Global terrestrial stilling: does Earth’s greening play a role?

Previous studies have documented that surface wind speed (u) has been increasing over the ocean but decreasing over land for the past several decades. The decreasing u at the surface over land has been referred to as terrestrial stilling. A plausible hypothesis for terrestrial stilling is an increase in surface roughness associated with changes in land surface (e.g. enhanced vegetation growth, landscape fragmentation or urbanization). One of the most widespread land surface changes is enhanced vegetation leaf area index (LAI) known as greening, particularly over the middle to high latitudes of the Northern Hemisphere where strong stilling is observed from weather station data. In this study, we examine the hypothesis that enhanced vegetation LAI is a key driver of global terrestrial stilling. We first characterized the trend in u over the ocean using long-term satellite altimeter measurements, and the trend in u over land using continuous wind records from 4305 in situ meteorological stations. We then performed initial condition ensemble Atmospheric Model Intercomparison Project-type simulations using two state-of-the-art Earth system models (IPSL-CM and CESM) to isolate the response of u to the historical increase in LAI (representing the greening) for the period 1982–2011. Both models, forced with observed sea surface temperature and sea ice and with LAI from satellite observation, captured the observed strengthening of Pacific trade winds and Southern Ocean westerly winds. However, these simulations did not reproduce the weakening of surface winds over land as significantly as it appears in the observations (−0.006 m s−1 versus −0.198 m s−1 during 1982–2011), indicating that enhanced LAI (greening) is not a dominant driver for terrestrial stilling.


Introduction
Surface wind speed (u) is the kinetic energy of the atmosphere that influences turbulent fluxes of water, energy, and momentum between the atmosphere and land or ocean surfaces. A change in u is an important manifestation of climate change, in addition to temperature and precipitation trends. While u over oceans has been increasing during the past 30 years (Wentz et al 2007, Young et al 2011, Zheng et al 2016, u measured at in situ weather stations at 10 m above ground has significantly decreased over land, particularly in the middle to high latitudes of the Northern Hemisphere (Vautard et al 2010, McVicar et al 2012b, Wu et al 2018. This phenomenon has been referred to as terrestrial stilling (Roderick et al 2007), and it has important scientific, socioeconomic, and environmental implications for fields including, but not limited to, agriculture and hydrology (McVicar et al 2012a), wind energy (Pryor et al 2006), the insurance sector (Sparks et al 1994), air pollution dispersion (McVicar and Roderick 2010), dust storms (Cowie et al 2013, Fan et al 2014, extreme temperatures (Horton et al 2015) and wind-related hazards and catastrophes (Kim and Paik 2015). For example, one of the most direct impacts of stilling is on long-term wind power industry that is rapidly developing and requires large infrastructural investments. So far, the drivers and physical mechanisms behind terrestrial stilling have not been determined, and thus it is unclear whether the observed stilling will continue (Vautard et al 2010, Peterson et al 2011, McVicar et al 2012b, Wu et al 2018. Decreasing u over land can be attributed to the driving force (i.e. a change in large-scale atmospheric circulation) and/or the drag force (i.e. an increase in surface roughness associated with land surface change). If the trend in u is mainly driven by changes in the force term, the decrease in u should be observable over both land and ocean and should also be present in climate reanalysis products (Pryor et al 2009, Vautard et al 2010. However, satellite altimeter measurements show a widespread increase of u over the ocean (Wentz et al 2007, Young et al 2011. Neither NCEP/NCAR nor ECMWF ERA-interim reanalysis products produce any trend in u over land , Pryor et al 2009, Vautard et al 2010. Therefore, the terrestrial stilling was hypothesized to be due to an increasing drag force on the atmosphere from the land surface, rather than a decreasing driving force (Vautard et al 2010).
One of the most widespread changes in the Earth's land cover is the enhanced vegetation growth during the past several decades, which is observed from satellites as the greening of the Earth (Zhu et al 2016). The global vegetation leaf area index (LAI, as a proxy for land greening) has enhanced by 8% during the past 30 years as indicated by satellite-observations (Zhu et al 2013(Zhu et al , 2016. The greening is particularly large over middle to high latitudes in the Northern Hemisphere (Myneni et al 1997, Xu et al 2013, Zhu et al 2016, where terrestrial stilling is observed (McVicar et al , 2012a(McVicar et al , 2012b. Enhanced vegetation LAI is associated with increased primary productivity and biomass, including tree height. Depending upon carbon allocation, an increment of net primary productivity is redistributed among root, stem and trunk. The latter will increase plant height and also surface roughness, which can decrease u. Vautard et al (2010) found a significant and negative spatial correlation between increased vegetation greenness (spring-summer normalized difference vegetation index (NDVI) trend) and stilling (u trend) across the stations in the Northern Hemisphere. They hypothesized that the cause of the terrestrial stilling could be an increase in surface roughness associated with enhanced vegetation activity, but NDVI is not directly related to LAI, and saturates with biomass. Further, NDVI can also change from vegetation structural shifts like the development of shrubs, as observed in dry regions and Arctic tundra (Frost andEpstein 2014, Tape et al 2006).
More importantly, correlation does not imply causation. The correlation could be derived from measuring a linear relationship due to factors that have nothing to do with NDVI and u. The role of enhanced vegetation growth (the greening of the Earth) in the terrestrial stilling remains unknown. The objective of this study is to test whether enhanced vegetation growth contributed to the terrestrial stilling observed over the past decades. This question is addressed by combining in situ and satellite observations of winds together with simulations of coupled land-atmosphere models.
In this study, we first examine the extent and patterns of changes in u during the period 1982-2011 using records from in situ meteorological stations over land and long-term satellite-derived observations over the ocean (Wentz et al 2007, Remote Sensing Systems 2015. We selected 4305 in situ stations that have more than 20 years of wind records from 26 247 meteorological stations (Global Surface Summary of the Day-GSOD). To isolate the wind response to changes in LAI, we ran two state-of-the-art coupled land-atmosphere models (the Community Earth System Model, CESM (Hurrell et al 2013) and the Institut Pierre Simon Laplace Coupled Model, IPSLCM (Marti et al 2010)) forced with satellite-observed LAI over a 30 years period . Note that LAI controls photosynthetic primary production, and is a key variable impacting vegetation growth in all models. We performed ensemble decadal Atmospheric Model Intercomparison Project (AMIP)-type simulations to make sure that the simulated climate is comparable with the historical climate (He and Soden 2016). The accuracy of the CESM and IPSLCM simulations in reproducing the historical climate change and the biophysical climatic feedbacks was evaluated by Zeng et al (2017Zeng et al ( , 2018a.

Observations
Satellite-derived observations of monthly average u over the ocean were obtained from the RSS V7 DISCOVER microwave radiometer data set on a 1°g rid for the period 1988-2011 (Wentz et al 2007, Remote Sensing Systems 2015. The satellite microwave radiometer, including SSM/I F08 through F15, SSMIS F16 and F17, WindSat, and AMSR2, are first inter-calibrated at the brightness temperature level, and are then used to construct the merged, 1°wind speed data set using a consistent processing methodology for all sensors (Remote Sensing Systems 2015). This data set has been carefully inter-calibrated and consistently processed. The in situ measurements of u (at 10 m height) over land were obtained from the GSOD (ftp://ftp.ncdc.noaa.gov/pub/data/gsod). This meteorological dataset is produced by National Centers for Environmental Information (NCEI) and includes 26 247 in situ stations based on the data exchanged under the World Meteorological Organization World Weather Watch Program. NCEI eliminated random errors in the original data, and performed extensive automated quality control. For this analysis, we only selected stations with more than 20 years of wind records during the period of 1982-2011. In determining whether a year of data is available, monthly mean u was calculated if more than 50% of the days in a month were available, and annual mean u was calculated only if all the months in the year were available. As a result, a total of 4305 meteorological stations over continents were selected. Current climate models simulations are very uncertain, particularly at regional scales, partly as a result of sea surface temperature and sea ice (SST) biases in coupled ocean-atmosphere models (He and Soden 2016). To improve the projections of anthropogenic climate change, we performed simulations with the land-atmosphere coupled models forced with observed SST during the past several decades. We obtained SST data for the period of 1982-2011 from the AMIP (http://pcmdi.llnl.gov/projects/amip/). This is a dataset grounded on observations and widely used in the AMIP-type simulations. There are significant increasing trends (significant level, p<0.05) of SST over most of the ocean except for the tropical Pacific Ocean and the Southern Ocean (figure S1(a) is available online at stacks.iop.org/ERL/13/124013/ mmedia). Overall, global SST significantly increased by 0.09°C per decade during the study period (p<0.01).
The greening of the Earth is represented by LAI observations from satellites (Zhu et al 2016). Monthly LAI was obtained from the advanced very high resolution radiometer 8 km global LAI product, which is derived from the third generation GIMMS NDVI3g from 1982 to 2011 (Zhu et al 2013). During the past 30 years, global land LAI has increased significantly by 8%, dominated by growth in the Amazon rainforest, the Sahel, Eurasian boreal forests, Europe, India, the eastern United States and southern China (figure S1(b)).

Numerical experiments
We performed several simulations with an ensemble of variable initial conditions using two state-of-the-art ESMs: the IPSLCM (version 4) from the Institute Pierre Simon Laplace at France (Marti et al 2010), and the CESM (version 1.3) from the National Center for Atmospheric Research at United Sates (Hurrell et al 2013). The horizontal spatial resolutions of IPSLCM and CESM are 3.75°×2.5°and 1.875°×2.5°, respectively, while the time step is 1.5 min. The methodology follows Zeng et al (2017Zeng et al ( , 2018a, where details of the simulation protocols can be found in the methods section and the supplementary information. We conducted two 30-year long AMIP-style experiments (AMIP_STD and AMIP_LAI) with the two models. Both experiments are constrained by observed SST from 1982 to 2011. In addition, while the LAI in AMIP_STD is forced with climatological monthly values, AMIP_LAI uses observed monthly LAI maps from 1982 to 2011. To reduce the uncertainty due to initial conditions (Lorenz 1963, Daron andStainforth 2013), all experiments contain a large ensemble with different initial conditions (n=30 for both models). The initial conditions of all realizations are identical, derived from a model output of an unperturbed 60 years run. The ensemble average was then used to characterize the response of u to prescribed changes in SST and LAI.
Surface u at 10 m was extracted from the model simulations. Similar to the reality, u is dependent upon drive and drag forces in the models. The drive term is pressure gradient force primarily owing to the simulated heterogeneity of temperature. It drives u over both the ocean and land surfaces. The drag term is caused by surface roughness associated with land surface characteristics. The drag term is negligible for u over the ocean, but not over land. SST change drives changes in climate (Kosaka and Xie 2013, He and Soden 2016), causing warming gradient that changes pressure gradient and u over land. LAI change has effects on both force and drag terms. On one hand, it triggers changes in land surface air temperature (Zeng et al 2017). The LAI-induced cooling has heterogeneity and this will superimpose on the force term.
On the other hand, surface roughness, similar to surface albedo, is computed from the simulated/prescribed state of the vegetation. Generally, models use LAI to represent vegetation growth and calculate photosynthetic primary production. Increased biomass is redistributed to root, stem and trunk, and is further used to calculate tree height. The simulated height of trees is then used to compute surface roughness (Su et al 2001), which influences the drag term. Below we take the land surface model in IPSLCM (details see Krinner et al 2005) as an example. As for photosynthesis at the canopy level, the assimilation rate (A c ) is integrated with LAI: where l is the cumulative LAI, A n is assimilation at the leaf level. Thus enhanced LAI increases photosynthetic primary production and carbon assimilation, which is further allocated to the biomass in different compartments (e.g. leaves, stems, roots, fruits) based on the allocation scheme in the model. The simulated biomass in the wood (B m ) is further used to calculate the wood's diameter (D): where π is the ratio of a circle's circumference to its diameter, p 1 is a parameter relevant to pipe density (40, dimensionless), p 2 is a tune parameter relevant to pipe (0.5, dimensionless). The equation for D and the behind parameters are developed based on the pipe model theory (Shinozaki et al 1964). D is further used to calculate vegetation height (H): where p 3 is a tune parameter relevant to pipe (0.3, dimensionless). Finally, the surface roughness length (z 0 ) is calculated based on the simulated vegetation height and vegetation fractions: where H PFT is vegetation height for each PFT, z _ 0 oh is a parameter to get z 0 from height of canopy (0.0625, dimensionless), z _ 0 bare is bare soil roughness length (0.01 m), z _ 0 nobio is roughness length for other nonvegetation surface types (e.g. ice, lake, ocean; nobio is the number of the types) (0.001 m), F PFT is vegetation fraction for each PFT, F bare is fraction for bare ground, F nobio is fraction for each non-vegetation surface types.  figure 1(g)). In North America and Europe, decreasing u primarily occurred during the last two decades (figures 1(c) and (d)). In Central Asia and India, the rapid decrease of u occurred during the first two decades and the decreasing rate slowed down during the last decade (figures 1(e) and (g)). Meanwhile, in Eastern Asia, mean annual u increased during the first two decades, and then decreased during 2002-2011 ( figure 1(f)). The widespread decreasing u at the surface over the northern mid-latitudes during the past 30 years is documented in several studies (e.g. Vautard et al 2010, McVicar et al 2012b. There are fewer stations in other regions, in particular the tropics and the Southern Hemisphere ( figure 1(a)). It is thus difficult to conclude the wind speed change over these regions.

Results and discussion
Overall, the trends in u for the four seasons follow similar patterns with those in mean annual u (figures 2 and S2). During the past three decades, mean u averaged over all weather stations significantly declined by  f2)). The reason why the declines of u in North America, Europe and Eastern Asia are much weaker than that in Central Asia and India is that there are more stations with positive trends in u in these regions ( figure S2). For example, many stations surrounding the Mediterranean show positive trends in u during the period 1982-2011 (figure S2), making a relative weak trend in u averaged over Europe (figure 2). During JJA, the trend in u averaged over Europe is even non-significant (p=0.56, figure 2(c2)).

Modelled response of u over ocean and land to SST changes
To isolate the climate response of u to the greening of the Earth, we ran the CESM and IPSLCM models forced with observed SST and variable/climatological LAI. We first tested the ability of the models to represent the response of u to SST-driven climate change. Note that simulated climate change driven by SST change compares well with the historical spatialtemporal variations in precipitation, temperature and evapotranspiration, as shown in Zeng et al (2018a, figures S3-5 in the paper). The DISCOVER satelliteretrieved u field over oceans during the past decades (Wentz et widespread negative trend in u over the Western Pacific and the Northern Atlantic ( figure S3(a)). This is in line with previous studies that documented a decreasing u along with weakening circulation under global warming (Held andSoden 2006, Vecchi andSoden 2007). However, strongly positive trends in u are observed over the Central Pacific and for the Southern Ocean westerly winds ( figure S3(a)). Averaged over Hereafter, all linear trends were calculated by robust regression (Rousseeuw and Leroy 2005) to reduce the impact of extreme mean annual u. Inset robust fit lines and numbers show trend of regional u during the study period. Dotting indicates a significant trend (p<0.05).
all the oceans, there is no significant trend in annual u during the period 1988-2011 ( figure S3(b)). However, Central Pacific trade winds significantly increased at a rate of 0.41 m s −1 per decade, while the Southern Ocean westerly winds increased by 0.22 m s −1 per decade. Both trends mainly occur during 1992-2011 (p<0.05; figures S3(c), (d)). It has been reported that the increasing wind over this region has a profound effect on global climate, as it is suspected to be the cause of the recent warming hiatus (Kosaka andXie 2013, England et al 2014).
The simulations showed that SST-forced u is generally decreasing in regions where SST are increasing, and increasing in the regions where SST are decreasing (figures 3(a), (d)). Spatially, SST increased over the entire ocean with the exception of the Eastern Pacific and Southern Ocean (figure S1(a)). Significant decreasing u is found over the western Pacific Ocean, the Indian Ocean, the Equatorial Atlantic and the Western North Atlantic, where SST have warmed during the past decades (figures 3(a), (d)). In contrast, SST over the Eastern Pacific and Southern Ocean has been decreasing. Both models simulate increasing u over these regions. Interestingly, although the warming rate is much greater than the cooling rate (figure S1(a)), the strengthening of u where SST decreased is much stronger than the weakening of u where SST increased (figures 3(a), (d)). Comparing modelled u with satellite-observed u, the pattern of the trends of u over the ocean are well reproduced by the models, particularly the strengthening of the Pacific trade winds  Combined with previous findings that the tropical Pacific cooling is a result of the unprecedented strengthening in Pacific trade winds (England et al 2014) and the results from this study that SST trends explain the strengthening in Pacific trade winds, we hypothesize the existence of a strong positive coupling between SST and Pacific trade winds. The positive feedback indicates that the climate system is sensitive to initial conditions (Lorenz 1963) over this region (approximately Nino 3.4 region): the strengthening in Pacific trade winds cools the ocean SST in the eastern Pacific and pushes the warm pool to the west, which reinforces Pacific trade winds and further cools the ocean SST in the eastern Pacific; globally, the cooling over this region feeds back to the global climate system and has led to the recent ongoing warming hiatus (Kosaka and Xie 2013).
If, like that of the change in u over the oceans, the change in u over land is driven by the global circulation change, the historical observations of u change over land should be reproduced by the SST-driven model simulations as well. However, in both model simulations (IPSLCM and CESM, static LAI simulations), the magnitude of the SST-forced u trend for the last 30 years is an order of magnitude smaller than the observed (negative) u trend over land (figures 4(a), (b) versus figure 1(a)). While the sign of SST-forced u trend is consistent with the observed u trend over Eurasia and India (negative trend), it is opposite to the observed u trend over North America (widespread positive trend; figures 4(a), (b) versus figure 1(a)). Similar results are found in the NCEP/NCAR and the ECMWF ERA-interim reanalysis products: the observed u trends over ocean have been well captured, while that over the northern mid-latitudes was not (figures S4-5) (Kalma et al 2008, Pryor et al 2009, Vautard et al 2010. According to these model results, the SST-driven circulation change seems not a determinant for the global terrestrial stilling. We note, however, more details about the model uncertainty should be investigated in future. 3.3. Modelled response of u over land to the greening of the Earth Next, we tested the hypothesis that increase of LAI contributed to the stilling. If this hypothesis is valid, forcing the models with the observed LAI change in addition to SST should reproduce the in situ observed surface u change over land. However, even though the models are constrained by both observed LAI and SST, there are substantial discrepancies between the modelled trend in u (figures 4(c), (d)) and the observed trend in u over land ( figure 1(a)). We isolated the contribution of LAI trends to u change over land from the difference (Δ) between ensemble simulations with observed LAI change (AMIP_LAI) and climatological  figure 1(b)).
As described in the Method (2.2), the greening of the Earth impacts on u over land through two mechanisms: (1) it triggers temperature heterogeneity and induces changes in pressure gradients; (2) it increases surface roughness length (z 0 ) and the drag force, which slows down u over land. To estimate the LAI-induced change in z 0 , we extracted the modelled z 0 from two equilibrium simulations (Zeng et al 2018a) which are two 60-year long simulations that share the same setup but are forced by different seasonal LAI maps (the 1980s: the average observed from 1982 to 1986 versus the 2010s: the average observed from 2007 to 2011). While the relative change of LAI is 8% between the 1980s and 2010s, the modelled relative change in z 0 is only 1.3%, changing from 0.2569 in the 1980s to 0.2601 in the 2010s. Assuming an unchanged pressure gradient, the change in z 0 leads to 0.35% decrease of surface wind speed at 10 m height (Kelly and Jørgensen 2017), which is −0.013 m s −1 during 1982-2011. It is still an order of magnitude smaller than the observed trend from in situ observations (−0.198 m s −1 , figure 1(b)).
Spatially, as shown in figures 4(e)-(f), modelled significant decreasing u occurs in the regions where LAI is significantly increasing ( figure S1(b)), including Europe, Central Asia, India, Southern China, Eastern United States and Southern Australia. There is no robust signal of the modelled response of u to greening in North America (−0.003 m s −1 per decade in IPSLCM and 0.012 m s −1 per decade in CESM; figures 5(b1), (b2)). In Europe, CESM produces a significant decrease of regional u as a response to LAI D at a rate of −0.018 m s −1 per decade (p<0.01,  ) and CESM (f). The response is calculated as a difference between the trends in AMIP_LAI and AMIP_STD. Dotting indicates a significant trend (p<0.05). Note that the colour axis scaling is only 20% of that for the observed u trend as shown in figure 1(a). figure 5(c2)), which is 28% of the observed stilling in this region ( figure 1(d)). In Central Asia, both IPSLCM and CESM produced a decrease in regional u as a response to LAI D (−0.007 m s −1 per decade; only 4% of the observed stilling in the region). In East Asia, there is no significant LAI driven trend in u in the simulations of both models (p>0.05; figures 5(e1), (e2)). In India, IPSLCM produced a significant decrease in regional u as a response of LAI, D but the rate is only 3% of the observed decrease in the region.
Therefore, the simulations indicate that the modelled responses of u to the greening remain orders of magnitude smaller than the observed u trends over land. It rejected the hypothesis that the greening of the Earth is a driver of the global terrestrial stilling. Model simulations depends strongly on how roughness is parameterized as a function of LAI in land surface models and how good is their aerodynamic coupling with the atmosphere. Given the large uncertainties behind the parameterization and the processes in the model, we examined this model results with more analyses directly based on the observations (see next section).

More evidence based on the coherence between the observed changes in u and LAI
The spatial coherence between satellite greenness (NDVI) and stilling over the middle to high latitudes in the Northern Hemisphere is a reason why we hypothesized that the greening is a determinant of the terrestrial stilling (also see Vautard et al 2010). Temporally, we did find a significant correlation between mean annual u and LAI averaged over the boreal Northern Hemisphere (north 30°) for the period 1982-2011 (R=−0.53, p<0.01). But this is an apparent correlation as it becomes non-significant when both time series are detrended (R=0.15, p=0.44). Furthermore, the correlation coefficients between detrended regional averaged u and LAI over all the individual regions with sufficient in situ stations are not significant (all with p>0.05; e.g. North America,p=0.86;Europe,p=0.54;Central Asia,p=0.20;East Asia,p=0.62;and India,p=0.95). Even for the long-term trend, the northern mid-high latitudes have been greening from 1981 to 1999 (Myneni et al 1997), yet during this period, mean annual u did not decrease over North America (figure 1(c)) and Europe ( figure 1(d)), both with p-value higher than 0.05. In East Asia, while the regionaverage LAI significantly increased for 1982-1999 (p<0.05; figure S1(b)), the region-average u is increased rather than decreased during the same period (0.035 m s −1 per decade, p=0.06; figure 1(f)).
Furthermore, the coherence between LAI and u trends does not exist everywhere if we examine it at a finer scale. For example, vegetation surrounding the Mediterranean is greening but u observed by many stations in this region is increasing during the past decades (figures S1(b) and 1(a)). In southern China and Japan, vegetation is greening but u is also increasing. Vegetation is browning (LAI decrease) in southern Alaska, central Canada and northern Russia, yet u in these regions is decreasing rather than increasing (figures S1(b) and 1(a)). In addition, if the surface stilling was mainly caused by the greening, the stilling should be larger in the growing season over deciduous forests. However, the observed decline of u in DJF (boreal winter) is 64% more than the decline in JJA (boreal summer). The rapidest decline in u occurred in Central Asia during boreal winter where and when widespread browning rather than greening is dominated (figures 2(d4), S6(d)). This decreasing rate (−0.227 m s −1 per decade, p<0.01) is much more than that in any other region and during any other seasons. Although widespread greening has occurred over Europe during the boreal summer (JJA; figure  S6(b)), there is no significant trend in u averaged over Europe during this season (p=0.56, figure 2(c2)).
Lastly, if the terrestrial stilling is mainly caused by an increase in surface roughness associated with the local vegetation change, there should be a significant correlation between the u observed by in situ station and LAI surrounding the station for the past decades . However, we found no significant correlation between the detrended mean annual u and the detrended average LAI surrounding each station (0.5°×0.5°) at the majority of the stations (3842/ 4305, 89%; figure 6). These results support the simulation results that the greening is not a dominant driver of the terrestrial stilling during the past decades.
3.5. Caveat and future studies to attribute the terrestrial stilling Models are imperfect, resulting in uncertainties in model simulations. As for the LAI-induced change in driving force, models are required to accurately simulate the responses of evapotranspiration and albedo, so that changes in air temperature and air pressure gradient can be reproduced by the models. However, because of systematic error in evapotranspiration partitioning, the CESM model failed to represent the response of evapotranspiration to the greening (Zeng et al 2017), and thus the simulated u change is uncertain. As for the LAI-induced change in drag force, models are supposed to correctly simulate the response of plant height to LAI change via vegetation grow models, and also the response of surface roughness to plant height change based on boundary theories. Yet, we cannot evaluate the responses of plant height and surface roughness to LAI change owing to the lack of global observations of plant height change and surface roughness. Regional long term measurements of plant height and surface roughness will be useful to evaluate the roughness parameterization in the model, and to reduce the uncertainty in attributing the terrestrial stilling to changes in surface land cover. Last, while the uncertainties exist in the model simulations, the results are supported by several evidence based on the coherence analyses between the observed changes in u and LAI (see the section 3.3).
Another caveat is that the large scale deforestation (e.g. Southeast Asia during 2000-2014 (Zeng et al 2018b)) has not been considered in this study. Largescale deforestation may intensify rather than weaken u, as it leads to a decrease of surface roughness. However, surface roughness could increase due to the deforestation induced forest fragmentation (Taubert et al 2018), and thus deforestation may also weaken u depending on its spatial pattern. Current high-resolution satellite observations of land cover change make it possible to detect historical change in forest fragmentation (e.g. Taubert et al 2018;Zeng et al 2018c). Yet, the current model horizontal resolution may be not fine enough to characterize the heterogeneous nature of the landscape (e.g. regional deforestation and forest fragmentation). In addition, fluffy cumulus clouds may act like anchors to create drag and slow the winds. But they are too small to register on the resolution for current climate models (Nuijens et al 2018). Future study using higher-resolution models could resolve the issue whether finer characterizations of land surface and/or clouds is an important element of the wind stilling.

Conclusion
We evaluated the global trend in u with 4305 in situ meteorological observations over land for the period 1982-2011, and the long-term satellite altimeter measurements over ocean for the period 1988-2011. The results presented here are consistent with previous results of a widespread negative trend in u over land, particularly the northern mid-latitudes, and a profound positive trend in u over Central Pacific Ocean and Southern Ocean. To isolate the response of u to the greening of the Earth for 1982-2011, we performed ensemble decadal AMIP-type simulations by prescribing the observed LAI and SST in two state-of-the-art ESMs (IPSL-CM and CESM). The SST-driven global circulation change can fully explain the observed change in u over the ocean, specifically the strengthening of Pacific trade winds and Southern Ocean westerly winds, proving the ability of the models to represent the response of u to climate change. The simulations showed that SST-driven global circulation change can partly explain the stilling over the Eurasia, but not North America. As for the u response to greening, regionally, the greening in Europe can explain 28% of the observed decrease in u over the region. However, globally, increasing LAI during the last 30 years (+8%) can only explain in the two models formulations to a decrease of terrestrial u by −0.006 m s −1 , an order of magnitude smaller than the observed trend from in situ stations (−0.198 m s −1 ). Therefore, the simulations by both models reject the hypothesis that the greening of the Earth is a key driver of the global terrestrial stilling. This conclusion is also supported by the incoherency between the observed changes in u and LAI (regional, temporal and seasonal variations). As one of the drastic land use changes surrounding all in situ stations, urbanization plays a vital role in shaping the landscape on the Earth. The decrease of near-surface wind speed is likely to be closely associated with the changing surface aerodynamic characteristics due to urbanization near the stations at regional scale (e.g. Zha et al 2017, Liu et al 2018. Future studies should be dedicated to exploring the connection of urbanization with the global wind stilling. A prerequisite is a high-resolution, temporalconsistent data product that can accurately capture the spatial and temporal patterns of urbanization over the globe.