Directional Climate Trend, Intensified Intraannual Variability, and Changes in Land Cover Drive the Dynamics of Vegetation Greenness in Peri‐Urban China During 2001–2015

Spatiotemporal dynamics of remote‐sensed vegetation indices and derived variables such as land surface phenology have often been studied as a function of climate with implied assumption of negligible land‐cover change. However, this assumption is not valid for areas with intensive human activities such as fast‐developing countries. To address the impact of land cover on vegetation greenness, we analyzed the trend of EVI59, enhanced vegetation index in May–September, within the 25‐km2 patches centered at 674 meteorological stations located mostly in peri‐urban areas of China, and the impacts of land cover and intraannual climate variability on the EVI, during 2001–2015 by means of linear mixed‐effect model. Impacts of land cover were assessed with the enhancement of sum of the squared difference with respect to the climate model. The climate models explained on average 60% sum of the squared difference of the full models (climate plus land cover), and this proportion reached 83% and 94% for the temperate grassland and the high‐cold Tibet, respectively. Including land cover enhanced on average 40% of the sum of the squared difference, and the enhancement is over 60% for the dense‐populated warm‐temperate‐deciduous forest and subtropical‐evergreen forest zones. The impact of land cover not only depends on the intensity but also on the type of land‐cover change. Forest regrowth in east China and vegetation growth from bare land in temperate desert significantly enhanced greenness but the shift from agriculture to shrubs in temperate grassland may not significantly alter the greenness. We showed that climate variability is important for EVI in all zones except Tibet, and the climate variability contributed on average 64% of climate impacts.


Introduction
Remote-sensed vegetation greenness can indicate many properties of terrestrial ecosystems (Ollinger et al., 2008;Soininen et al., 2015). The high absorbance of red light and the high reflectance of near-infrared light by green vegetation provide the spectral contrast used to construct vegetation indices, such as normalized vegetation index (NDVI) and enhanced vegetation index (EVI; Huete et al., 1997). The remotely sensed vegetation indices have long been used as a spatiotemporally consistent measurement of canopy structure (Gamon et al., 1995), leaf area index (LAI; Viña et al., 2011), leaf nitrogen content (Knyazikhin et al., 2013), and leaf biomass (Anderson et al., 1993;Tucker, 1979). These vegetation indices helped to derive important parameters of ecosystem processes such as carbon assimilation (Veroustraete et al., 1996) and evapotranspiration (Zhang et al., 2016). Vegetation greenness has also been used to quantify the ecosystem feedback to the climate system (Findell et al., 2017), to monitor the land-cover change (Joshi et al., 2016), and to assess the disasters such as fire and flood (Bello & Aina, 2014). The dynamics of vegetation greenness, which follows essentially the pace of land surface phenology, responds to the variability/fluctuation and changes in climate (Cleland et al., 2007;Menzel et al., 2006).
Modeling of vegetation greenness and its derived variables with respect to the changes of climate and other environmental features is an important aspect of remote-sensing applications (Pettorelli et al., 2005) and involved many studies. Vegetation indices were directly regressed on precipitation and temperature to quantify the impact of climate change on vegetation greenness (Hao et al., 2012;Xu et al., 2016). To reconcile the effects of drought-induced green-up and drought-induced tree death on the carbon flux in the Amazon area, a set of relationships between EVI and climate variables, such as vapor pressure deficit and plant available water, were examined (Brando et al., 2010). The remote-sensing recorded die-off of woodland vegetation in the southwest of North America was found to be caused by drought (Breshears et al., 2005). The changes in NDVI of advanced very high resolution radiometer were found as results of shifted rainfall and temperature (Piao et al., 2004). The relationship between advanced very high resolution radiometer NDVI and climate variables was also used to parameterize the terrestrial ecosystem simulation models (Gao et al., 2000(Gao et al., , 2007. Phenology has also been analyzed to show diverse responses to climate (Cleland et al., 2007;Menzel et al., 2006;Zhang et al., 2004). The factor of land-cover change has been assumed to be a minor signal and neglected in the above analyses, so that the greenness or variable derived from vegetation indices was solely regarded as a function of climate. The assumption is valid for areas where population density is low or where land-cover change is a much weaker signal compared to climatic change.
However, this assumption of negligible land-cover change may not be valid, especially for regions such as China, India, and other developing countries where rapid economic development has induced substantial land-cover changes in and around urban areas (Bren d'Amour et al., 2017;Pandey & Seto, 2015;Zhang et al., 2014). Ignoring the impact of land-cover change may substantially bias the relationship between vegetation greenness and climate change. To exclude the effects of land-cover change on vegetation greenness, Yu et al. (2017) derived the relationship between the nighttime warming and the EVI based on the Moderate Resolution Imaging Spectroradiometer (MODIS) imagery for patches without land-cover change in a tropic region. Similarly, the land surface phenology derived from MODIS EVI on patches without landcover change was regressed against climate variables . Phenoregion was proposed to detect areas with minimal nonclimate forcing for global phenological monitoring (White et al., 2005), and impacts of compositional changes in agricultural ecosystems on satellite-derived phenology trend were explored in the Midwest of United States . Efforts have also been made to disentangle signals of climate fluctuation from those of land-cover changes (Dieguez & Paruelo, 2017).
The land-cover changes in China over this century feature intensive urban expansion into the previous periurban agricultural land. Reforestation or agroforestry also takes place in the abandoned agricultural land after the rural residents flow to urban areas . It is expected that the most intensive land-cover change would be in the peri-urban areas in recent decades. In such circumstances, land-cover change needs to be included explicitly as explanatory variable in models of the spatiotemporal dynamics of vegetation greenness.
Climate variability is another important driver to vegetation greenness. It is predicted that the greenhouse warming increases the frequency of extreme climate events such as El Niño or La Niña (Cai et al., 2014(Cai et al., , 2015. The dynamics of vegetation indices also reflects temporal ecosystem stability and resilience (de Keersmaecker et al., 2014). Dynamics of climate can feature both trend and fluctuation. Comparing to the impact of the trend component, the warming-intensified fluctuation, including extreme climate events such as drought and storms, may have abrupt devastating impact on the vegetation and ecosystem function (Thornton et al., 2014). The question how the intensified climate variability translates into the dynamics of vegetation greenness and ecosystem functions still lacks a general solution.
This study is aimed to address the spatiotemporal dynamics of vegetation greenness for areas with most intensive land-cover change and for time period with great climate variability. Assuming that most meteorological stations in China are located in peri-urban areas with intensive land-cover changes, we analyzed the climate variability and land-cover changes and their impacts on vegetation greenness for areas surrounding the meteorological stations during 2001-2015 period. Including both land cover and climate as explanatory variables for the spatiotemporal dynamics of MODIS EVI in peri-urban areas will allow us to place a cap on the effects of land-cover change on vegetation greenness. The period of 2001-2015 has two strong La Niña events in 2007-2008 and 2010-2011 and one very strong El Niño event in 2015-2016 (https://ggweather. com/enso/oni.htm, https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php), thus features great climate variability. This feature allowed us to explore how climate variability, together with land-cover changes, impacts the vegetation greenness.
We hypothesized that the relative importance of climate versus land cover on vegetation greenness depends not only on the intensity but also on the types of land-cover changes in the particular climate regime. The impact of land-cover change is stronger than that of climate in areas where intensive land-cover changes cause great difference in greenness, such as reforestation in east China and vegetation regrowth from bare land in the northwestern desert. On the other hand, the impact of climate is stronger than that of land cover in areas where vegetations are strongly sensitive to climate change, such as the vegetation in the semiarid and arid zones, or where land-cover change did not cause significant changes in greenness. We also hypothesized that changes in climate variability account for most, that is, more than half, of the impact of climate change in most regions in China.
We acquired daily meteorological records from the China Meteorological Administration for 674 stations, and extracted MODIS EVI and annual land cover in patches of 25 km 2 centered on the stations. Intraannual variability of the climate variables was computed as the standard deviation of the monthly data. Linear mixed-effect model (LME) was used to derive the linear temporal trends of climate variables and land cover, and to quantify the response of MODIS EVI to the climate changes, climate variability, and land-cover changes.

Data Acquisition and Preparation
Daily meteorological data consisting of 674 stations ( Figure 1) with total precipitation, P (mm), maximum (T max ) and minimum (T min ) temperatures (°C), mean wind speed (m/s), total sunshine (hr), and mean relative humidity (%) in the period of 2000-2015 were obtained from the China Meteorological Administration. Daily reference evapotranspiration (ET 0 ) was computed according to FAO documents (Allen et al., 1998) based on the available meteorological data. The daily data were then averaged for temperatures or summed for precipitation and evapotranspiration for each month and each wet year (from 1 October to 30 September next year; Peters et al., 2014). Number of dry days during May-September (D 59 ) with daily precipitation less than 0.1 mm was counted for each year. The intraannual variabilities of the climate variables are characterized by their standard deviations computed from monthly data.
Annual land-cover data for China were classified by the authors group  based on MODIS imagery with the virtual ground-truth interpreted from high-resolution Google Earth images. For this analysis, we combined the land cover into six classes: forest (F c ), shrub (S c ), grass (G c ), sparse vegetation or bare land (B c ), agriculture (A c ), and urban land (U c ).
We extracted MODIS enhanced vegetation index (product MOD13Q1 V006) of a square area of 25 km 2 centered at the geographical coordinates of each meteorological station by means of Google Earth Engine (Gorelik et al., 2017). The mean EVI from May to September (EVI 59 ) was computed for each year. For most areas of China, May to September is the growing season of the deciduous plants. Annual land covers were also extracted from the land-cover maps of China for the patches, and the cover fractions were computed for the six land-cover types.

Linear Trend of EVI, Climate, Land Cover, and Climate Variability
We intended to analyze the peri-urban vegetation greenness in various zones of China. According to the 1:6,000,000 vegetation zonation map of China, the continental China is divided into eight vegetation zones: cold-temperate coniferous forest (CCON) and temperate coniferous-broadleaved mixed forest (MIXD) in the northeast, warm-temperate deciduous broadleaved forest (DBRD) in the north, subtropical evergreen broadleaved (GBRD) and seasonal rainforest (RAIN) in the southeast, temperate grassland (GRAS) in the north and northwest, desert in the vast northwest (DSRT), and high-cold vegetation in the Tibet Plateau (TIBT). Considering that the cold-temperate coniferous forest and the rainforest zones are relatively small in proportions of area, and have far fewer meteorological stations available, we grouped the cold-temperate coniferous forest zone into the mixed forest zone, and the seasonal rainforest zone into the subtropical evergreen broadleaved forest zone. Hereafter, we refer these zones as mixed forest (MIXD), deciduous (DBRD) and evergreen forest (GBRD), grassland (GRAS), desert (DSRT), and Tibet (TIBT). The distribution of the meteorological stations within the six zones ( Figure 1 and Table S1) shows the decreasing density of the stations from the heavily populated east to the sparsely populated western desert. Mean annual precipitation, temperature, and land covers for stations in each zone were illustrated in Table 1 as the background information.
The layered data structure, that is, each zone has multiple meteorological stations and patches, and each station/patch has multiple years of measurement, requires mixed/hierarchical model approaches. We applied LME to derive the linear slopes of EVI 59 and the major climate and land-cover variables by regressing them on time (year) for each vegetation zone with meteorological station as the grouping variable. The slope in the fixed effects signifies the mean trend in each vegetation zone, whereas the slopes of the random effects represent the deviation of trends of meteorological stations. Changes in EVI and climate/land cover were obtained as the slopes times the time interval from 2001 to 2015.

Analysis of the Impacts of Climate Variability and Land-Cover Change on Vegetation Greenness
LME was used to investigate the response of EVI 59 to land cover, climate changes, and intraannual climate variability, for each vegetation zone at annual scale with meteorological station as the grouping variable. The function lmer() in R packages "Lme4" and "LmerTest" (Bates & Machler, 2015;Kuznetsova et al., 2017) were used to take the advantage of the capability of stepwise selection of the independent variables to minimize the restricted maximum likelihood convergence criterion (similar to Akaike Information Criterion). The function lmer() also handles reasonably the problem of degrees of freedom associated with restricted maximum likelihood (Bates & Machler, 2015). LME also requires the variables to be in approximately the same order of magnitude to make it less error-prone in the extraction of eigenvectors and eigenvalues. To accommodate, we divided P and ET 0 by 1,000, S P (standard deviation of precipitation) and S ET0 (standard deviation of ET 0 ) by 100, D 59 by 365.25, and T max and T min by 10 to make them in the order of 1 or 0.1. All the land- cover variables are in fractions. After the variables are selected, we used function "lme()" in "nlme" package to rerun the models with first order of serial correlation.
For each vegetation zone, we fitted the following models for mean EVI 59 : a null model with only fixed and random intercepts (Model 0), a climate model with climate and its variability (Model 1), a model with only climate variability (Model 1.v), and a full model with both climate and land cover (Model 2), as explanatory variables. The coefficients of the fixed effects reflect the mean response of vegetation greenness of the patches to the climate and land cover at the level of vegetation zones, whereas the coefficients of the random effects discern the deviation of the response among the patches. LME does not offer a measure to quantify the variation in the data explained by the model. Several studies tried to derive the statistics similar to R 2 for ordinary multilinear regression. We used the "rsquared()" function in "piecewiseSEM" package to find R 2 for models in each zone. We also used the approach of Xu (2003) to compute the sum of the squared difference (SSD) in predicted EVI 59 between the models with climate, with and without land cover (Models 1 and 2), and the null model (Model 0). The relative importance of climate can be assessed with the fraction of the full model SSD explained by the climate model (ratio of SSD of Model 1 over that of Model 2). The relative importance of land cover is regarded as the fraction of SSD enhanced by incorporating land cover. To assess the impact of climate variability, we divided the SSD of the model of climate variability (Model 1.v) by that of the climate model (Model 1) and by that of the full model (Model 2).

Change of EVI 59 , Climate and Its Intraannual Variability, and Land Cover in 2001-2015
Our analysis is based on the vegetation patches in peri-urban areas. In the following description, when we say anything for a particular vegetation zone, we mean the mean trend/change or mean relationship for all the vegetation patches in this zone, not the continuous area. At most, our results and conclusions may reflect the discontinuous peri-urban areas of cities and towns.
The average changes of EVI 59 , climate, and land cover for each of the six vegetation zones were computed as the slopes of the fixed effects of the linear mixed-effect model times the time interval ( Figure 2). EVI 59 in all the six zones were significantly increased (p < 0.05). The increases of 0.029 and 0.033 in the mixed forest and the grassland zones, respectively, are the strongest, and the increase of 0.006 in Tibet is the weakest (Figure 2a). This direction of EVI 59 change does not match up with the direction of changes in precipitation P of the stations in the zones (Figure 2b). The patches in the grassland zone had the greatest increase in mean EVI 59 in correspondence with the average 59.6-mm increase in precipitation (p < 0.05), whereas in the deciduous forest zone the increase in EVI 59 was associated with the average −56.4-mm decrease in precipitation (p < 0.05). The mixed forest zone had the greatest increase in precipitation of 80.5 mm (p < 0.05), but the system is not as responsive as the grassland. The change in the number of dry days from May to September, D 59 , is strongly inversely correlated with that of the P (Figure 2d). It was significantly increased in Tibet (eight days), but significantly decreased in the mixed forest and the grassland (p < 0.05). Annual mean of daily maximum temperature, T max , significantly increased in Tibet, desert, and deciduous forest (Figure 2a; p < 0.05), which corresponded to the decreased P (Figures 2a and 2b). The two zones with increased P (mixed forest and grassland) showed decreased T max . The minimum temperature, T min , increased across six zones and the increases were significant in Tibet, evergreen forest, and desert ( Figure 2a; p < 0.05). Among the six vegetation zones Tibet had the greatest increase in both the maximum (0.37°C) and minimum (0.52°C) temperatures, exhibiting stronger warming in the high altitude. The reference evapotranspiration (ET 0 ) significantly increased in Tibet, desert, and evergreen forest (p < 0.05), and the steepest increase of 55.2 mm was found in Tibet, which corresponded to the greatest increases in both T min and T max (Figures 2a and 2b). On the other hand, ET 0 significantly decreased in mixed forest and grassland (p < 0.05), which may be caused by the decreased T max . Note. P, annual precipitation; T max , annual average of daily maximum temperature; T min , annual average of daily minimum temperature; F c , forest cover; S c , shrub cover; G c , grass cover; A c , agriculture cover; and U c , urban cover. MIXD, temperate mixed forest and cold-temperate coniferous forest; DBRD, warm-temperate deciduous-broadleaved forest; GBRD, subtropical evergreen-broadleaved forest; GRAS, temperate grassland; DSRT, temperate desert; and TIBT, Tibet high-cold vegetation.
The trend analysis also indicated increased intraannual climate variabilities, which reflect both seasonality and extreme events, as standard deviations computed from the monthly data (Figures 2c and 2e). The grassland had a significantly increased variability of precipitation (ΔS P = 4.9 mm; p < 0.05). The intraannual variability of maximum temperature, S Tmax , significantly increased for all the vegetation zones except grassland (p < 0.05), for which S Tmax was decreased by 0.16°C. The evergreen forest zone in southeast had the greatest increase in S Tmax (0.37°C), and the zones of Tibet and mixed forest also had large increases in S Tmax (0.32°C and 0.3°C, respectively). All the zones had increases in the variability of the minimum temperature, S Tmin , and the increase was significant in all zones (p < 0.05) except that in grassland and Tibet. The mixed forest in northeast had the greatest increase in the S Tmin (0.44°C). The changes in variability of the ET 0 seemed to correlate much with the change of ET 0 itself (Figures 2b and 2c). The mixed forest and grassland had the greatest decrease in S ET0 (−3.1 mm), whereas Tibet had the significant increase of 1.3 mm in S ET0 (Figure 2c; p < 0.05).
Agricultural land cover, A c , significantly retreated in all zones (p < 0.05), and the reduction was more than 10% in deciduous forest (−15.9%), evergreen forest (−12.5%), and grassland (−11.0%) ( Figure 2f and Table  S1). Analysis of land-cover change  indicated that there are large areas of newly reclaimed agriculture in the northeast China (the mixed forest zone) and in the desert in the west to compensate the overall agriculture reduction. For the peri-urban patches, we found that in the desert zone there are 23 out of the 94 patches with an average rate of 5.2% increase in agriculture cover. In the mixed forest zone, we found 9 out of 51 patches with an average rate of 4% increase in agriculture cover. Patches with increased agriculture are negligible in the other zones.
In contrast to the reduction of agricultural land, urban cover, U c , increased significantly across all zones (p < 0.05), and the increases were the strongest in deciduous forest (11.9%), evergreen forest (8.2%), and grassland (5.3%). Forest cover, F c , was increased in all the zones, and the change was significant (p < 0.05) except that in mixed forest and desert. Evergreen and deciduous forest zones had the steepest increases of 4.8 and 2.9% in F c , respectively (Figure 2f). Shrubs cover, S c , also increased in most zones, and the change of 7.0% in grassland was the greatest. Bare land cover, B c , showed a maximum reduction of 3.8% in the desert zone ( Figure 2f and Table S1). Table 2 shows the coefficients of fixed effects and the variables selected for the random effects in the climate model (Model 1). Each column of the table represents an equation for a vegetation zone. For example, the first column represents the equation for mixed forest zone as the following:

The Climate Model
Variables and coefficients out of the parenthesis are the fixed effects describing the mean relationship between EVI 59 and climate variables across all vegetation patches in the zone, whereas variables within the parentheses are the random effects and were evaluated for each station (i) to obtain the i th of a and b. Table 2 indicates that precipitation P appeared in the fixed effects and positively impacts the EVI 59 for all zones except evergreen forest where moisture is more abundant than other zones (Table 1), and the response of EVI 59 to precipitation is the strongest in desert and grassland zones with the greatest coefficient of 0.115 and 0.109, respectively. D 59 negatively impacts EVI 59 in evergreen forest. The minimum temperature, T min , positively impacts EVI 59 in deciduous forest, evergreen forest, and Tibet zones, and Tibet showed the strongest response with coefficient of 0.093 in correspondence to the strongest increase in T min (Figure 2a). The maximum temperature, T max , only appeared in the equation for the evergreen forest and negatively impacts EVI 59 with coefficient of −0.046. The variability of precipitation, S p , negatively impacts EVI 59 in all zones except Tibet, and the temperate grassland showed the strongest impact of S p with coefficient of −0.041. Similarly, S ET0 appeared in all zones except Tibet and negatively impacts EVI 59 , and the impact is strongest for the temperate grassland with coefficient of −0.146. The effect of the variability of minimum temperature, S Tmin , on greenness is not consistent: positive in evergreen forest and desert but negative in grassland. The variability of maximum temperature, S Tmax , appeared in the equation for deciduous forest and positively impacts EVI 59 with coefficient of 0.051.
Variables in random effects were preselected based on their ability to reduce the restricted maximum likelihood convergence criterion of LME model, and then reselected by the LmerTest. We found that in most of the cases the effective reduction of the criterion is caused by S ET0 , P, and D 59 in the descending order, as they appeared in four or five zones (Table 2). ET 0 appeared in the equations for evergreen forest and desert zones, and T min appeared in the equations for evergreen forest and Tibet zones. ET 0 and S ET0 are synthetic variables including radiation, temperature, relative humidity, and wind effects. Precipitation and evapotranspiration control the water budget of terrestrial ecosystems. D 59 describes the variation of daily rainfall. The random effects thus showed dominance by the water budget and its intraannual variations. The serial correlation coefficients varied from 0.31 for mixed forest to 0.77 for desert with the mean of 0.53 (Table 2).

The Full Model With Both Climate and Land Cover
The variables of the climate and land cover were combined in the regression of EVI 59 on climate and land cover (Model 2), and reselected by the LmerTest. Table 3 shows the coefficients of the fixed effects and the variables selected for the random effects for mean EVI 59 . Comparing to the climate model (Model 1), we found that some variables of the fixed effects in the previous climate model were eliminated when land-cover variables were added. For example, S p was eliminated from the equation of deciduous forest zone. T min only Note. ρ stands for first-order autocorrelation. P, annual precipitation; D 59 , number of dry days during May-September; ET 0 , reference evapotranspiration; T max , annual average of daily maximum temperature; and T min , annual average of daily minimum temperature. S p , S Tmax , S Tmin , and S ET0 are standard deviations of precipitation, maximum temperature, minimum temperature, and reference evapotranspiration, respectively. MIXD, temperate mixed forest and cold-temperate coniferous forest; DBRD, warm-temperate deciduous-broadleaved forest; GBRD, subtropical evergreen-broadleaved forest; GRAS, temperate grassland; DSRT, temperate desert; and TIBT, Tibet high-cold vegetation. Significant indicators: *for p value less than 0.05 and † for p value less than 0.1. appeared in the equation for Tibet, and S Tmin only remained for grassland zone. The sign of T max in the equation of evergreen forest zone became positive. However, the main theme of climate model on ecosystem water budget still persists in the full model (Model 2).
Forest cover, F c , appeared in the equations for all three forest zones in east China with positive effects on greenness (Table 3). Agricultural cover, A c , appeared in mixed and deciduous forest zones, as well as grassland, with positive effects on greenness. Bare ground cover, B c , was selected with negative effects for desert and Tibet zones where bare land decreased. Shrub cover, S c , was selected for deciduous forest zone. Urban cover, U c , was selected for evergreen forest zone with negative effect.
For the random effect, climate variable S ET0 remained for mixed forest, deciduous forest, and grassland zones, P presented in grassland, desert, and Tibet zones. D 59 remained for mixed forest, evergreen forest, and Tibet zones (Table 3). T min appeared only in the equation for Tibet, A c presented in all six zones except Tibet, and U c presented in the two mostly developed zones, that is, deciduous and evergreen forest zones. S c appeared in deciduous forest, grassland, and desert zones, in which shrub got proliferated (Figure 2f). B c presented in desert where bare ground largely decreased and F c appeared in evergreen forest where forest expansion is the highest (Figure 2f). Figure 3 shows how the full models fit to the mean EVI 59 . We found the R 2 as 0.93 for mixed forest, deciduous forest, and evergreen forest zones, 0.94 for grassland, 0.99 for desert, and 0.98 for Tibet.
The plot of SSD of the climate model (Model 1), climate variabilityonly model (Model 1.v), and full model with both climate and land cover (Model 2) showed that for annual mean EVI 59 (Figure 4), the SSD of the climate model accounted for 94% of SSD of the full model for Tibet, 83% for grassland, 53% for desert, 52% for mixed forest, 39% for deciduous forest, and 38% for evergreen forest. Therefore, land cover accounted for more than 60% of SSD of the full models for evergreen (62%) and deciduous forest (61%) zones, the two mostly developed zones with the most extensive land cover changes (Figure 2f). Land cover accounted for around half of SSD of the full models for mixed forest (48%) and desert (47%) zones, where the impacts of land cover and climate changes were comparable. The average percentage of the climate model is 60% over all zones. The impact of climate is much stronger than that of land cover, thus dominates in the high-cold Tibet as well as the grassland. On the other hand, the impact of land cover is much higher than that of climate for the evergreen forest and the deciduous forest zones where dense population is associated with intense economic activities.

Impact of Climate Variability on Mean EVI 59
The standard deviation of the climate variables computed from monthly data reflected mostly the climate seasonality, but it may also include the effects of extreme climate events. Intraannual variability appeared in the equations for all zones except Tibet (Table 3). The effects of intraannual variabilities are all negative except S Tmax in the equation for evergreen forest (Table 3); hence, climate variability in general negatively impacts vegetation greenness.
To look into the overall impacts of climate variability, we obtained the SSD of the model with climatevariability only (Model 1.v) which accounted for 91, 52, 48, 64, and 68% of SSD of the climate model (Model 1) for mixed forest, deciduous forest, evergreen forest, grassland, and desert zones, respectively Note. ρ stands for first-order autocorrelation. P, annual precipitation; D 59 , number of dry days during May-September; ET 0 , reference evapotranspiration; T max , annual average of daily maximum temperature; and T min , annual average of daily minimum temperature. F c , forest cover; S c , shrub cover; G c , grass cover; B c , bare ground cover; A c , agriculture cover; and U c , urban cover. S p , S Tmax , S Tmin , and S ET0 are standard deviations of precipitation, maximum temperature, minimum temperature, and reference evapotranspiration, respectively. MIXD, temperate mixed forest and cold-temperate coniferous forest; DBRD, warm-temperate deciduous-broadleaved forest; GBRD, subtropical evergreen-broadleaved forest; GRAS, temperate grassland; DSRT, temperate desert; and TIBT, Tibet high-cold vegetation. Significant indicators: *for p value less than 0.05 and † for p value less than 0.1.
( Figure 4). Referenced to the full model (Model 2), intraannual climate variability accounted for 47, 20, 18, 53, and 36% of SSD for the above zones. On average, the climate variability accounted for 64% of SSD referenced to the climate model, and 35% of the SSD referenced to the full model ( Figure 4b). Therefore, climate variability plays important roles in vegetation greenness, and climate variability accounted for more than 60% SSD of the climate model, thus dominated climate impacts on vegetation greenness in the zones of cold-temperate mixed forest, temperate grassland, and temperate desert.

Changes of Climate and Its Variability in the Patches Centered at the Meteorological Stations
For the first decade of this century, increased precipitation was reported for northern China (including the northeast), but decreased precipitation for southern China based on the observed data (Wang, Gao, Liu, et al., 2016). This shifted pattern might be a result of warming-induced northward rainfall shift (Putnam & Broecker, 2017). Our analysis showed that among the six vegetation zones, station-based precipitation in the temperate grassland in northern China and the temperate-cold-temperate mixed forest in northeastern China significantly increased, whereas precipitation in other regions decreased (Figure 2). The regions with increased precipitation have decreased maximum temperature partly due to the increased cloudy sky associated with rainfall ( Figure 2); on the other hand, the regions with decreased precipitation all have increases in maximum temperature. The increase in minimum temperature is found in all zones and the increases are greater than in maximum temperature except in warm-temperate deciduous-broadleaved forest, implying more nighttime warming than daytime warming (Reza, 2015). We found that direction of changes in intraannual standard deviation of precipitation is mostly associated with the changes in precipitation in all the zones except temperate desert. Thus, the increase in annual precipitation also means more concentrated precipitation profile due to increased seasonality and/or extreme climate events. Increasing standard deviation of precipitation in the northern grassland zone implies more concentrated rainfall in the growing season as the average precipitation profile of the area is that most rainfall concentrates in summer (Na et al., 2018). On the other hand, decreased precipitation in warm-temperate deciduous-broadleaved

Journal of Geophysical Research: Biogeosciences
forest, subtropical evergreen-broadleaved forest, and Tibet also implies less rainfall in summer. Increase in minimum temperature (Richard et al., 2017) on one hand contributed to elevated evapotranspiration and maintenance respiration; on the other hand, it may also prolong the vegetation growing season by increasing thermal time and delaying the end of the growing season. Number of dry days has been reported to increase over the period from 1961 to 2012 for most areas of China (Duan et al., 2017). For the period of 2001-2015, we found that the changes of the number of dry days in May to September are inversely correlated to that of the precipitation: that is, reduced number of dry days is associated with increased precipitation in the temperate grassland and mixed forest zones, whereas increased number of dry days is associated mostly with the decreased rainfall in other zones.
We also found that the increases in the intraannual standard deviation of the minimum temperature are greater than those in the mean minimum temperature in all zones except Tibet where the increase in mean minimum temperature is the greatest. Hence, the changes in seasonality and/or extremes of minimum temperature are greater than changes in mean. Increase in S Tmin means greater warming summer than warming winter or warming summer or cooling winter. We also found decreased minimum temperature in winter for northeastern China , and the winter chilling is assumed to reduce thermal time requirement and to advance the growing season. Climate models suggest that steep winter warming is temporary, and summer will eventually warm more than winter (Egan & Mullin, 2016).

Importance of the Moisture and the Temperature on Vegetation Greenness
Vegetation greenness reflects mostly the green leaf area of ecosystem which is mostly supported by the availability of water (Chapin et al., 2011;Cowling & Field, 2003). Therefore, precipitation and reference evapotranspiration are major controlling climate variables. This has been demonstrated in our study. Except in the moisture-abundant zone in subtropical south China, precipitation appeared as an important fixed effect on vegetation greenness (Table 3). As a synthetic climate variable, standard deviation of reference evapotranspiration, together with the precipitation and the number of dry days, dominate the random effects of the linear mixed-effect model for the six zones to differentiate the response of vegetation greenness among the vegetation patches. As to the minimum temperature, the effect of prolonging the growing season dominates over promotion of respiration and evapotranspiration; therefore, T min showed positive effect on mean EVI 59 in Tibet where low temperature is the major limiting factor. Although the daily maximum temperature is important for canopy photosynthesis which helps the leaf growth, its effect did not appear for most zones except the evergreen forest zone partly because of its correlation with the minimum temperature and precipitation.

Importance of Land-Cover Change on Vegetation Greenness
The main features of land-cover/land-use changes in the patches centered at the meteorological stations are agriculture reduction, urbanization, and reforestation/afforestation (Figure 2f). Agriculture cover is an important random effect to discern the responses of greenness among the patches in all the vegetation zones except Tibet in China (Table 3). Agriculture cover also appeared in three zones of mixed forest, deciduous forest, and grassland as fixed effect where agriculture retreat is extensive. However, urban cover only appeared as fixed and/or random effects in the heavily populated deciduous forest and evergreen forest zones where urban expansions are most intensive (Table 3 and Figure 2f). Reforestation contributed to the greening in China (Chen et al., 2019) as we found forest cover in the three forest zones in east China helped the increases in the peri-urban greenness (Table 3 and Figure 2f).
Land cover was included in the trend analysis of greenness for the city of Guangzhou in south China to illustrate the effects of abrupt land-cover changes (Zhu et al., 2016). Using the annual land-cover data for China, we found that for peri-urban areas of densely populated southeast China, land-cover changes explained much more variations in greenness than climate (Figure 4). Our result thus implied that land cover cannot be neglected in the analysis of vegetation greenness and derived quantities. This is especially important for densely populated areas with heavy economic activities.

Relative Importance of Climate Versus Land Cover
It is worth pointing out that the trend of EVI 59 , an important measure of canopy leaf area or biomass, does not follow pattern of precipitation ( Figure 2) as we found the greenness increased in all zones, whereas precipitation decreased in four of the six zones. Aside from the moisture balance, there are a few factors that tend to affect the canopy greenness: CO 2 fertilization, increased growing season due to raised temperature, and reforestation as one of the important part of land-cover changes. We did not include the CO 2 fertilization in this analysis due to the lack of common understanding of long-term impact of raised CO 2 concentration on canopy greenness. Some studies indicated that under the elevated CO 2 scenarios, the total leaf area would decrease and the root biomass increases (Körner et al., 2005), but others found that CO 2 enrichment did not affect the partition of carbon among plant biomass (McCarthy et al., 2010). All CO 2 enrichment experiments use step increase of CO 2 . It is not clear how the reality of gradual rise of CO 2 will affect the leaf biomass. The minimum temperature appeared in the equation for Tibet, and the positive coefficient may imply a lengthened growing season with minimum temperature, a result similar to what found in a previous study (Shen et al., 2016). Also, recent studies found that the area of surface water has increased substantially in Tibet due to warming which may alleviate the moisture stress on plants Wang et al., 2019).
The relative importance of climate versus land cover depends not only on the intensity of economy-driven land-cover/land-use changes indicated by the stacked bar height in Figure 2f but also on the types of land-cover changes and the particular climate regime. The deciduous forest and the evergreen forest zones are the major fast-developing regions in China, and the economy growth drove most extensive land-use/ land-cover changes as reflected with the highest stacked bars of two zones in Figure 2f. The impact of climate change dominates over that of the land-cover change in Tibet because the land-cover change is light (Figure 2f) and the limiting factor is mostly the minimum temperature.
There is also high rate of agriculture reduction and shrubland and urban expansion in the grassland (Figure 2f). The increased shrub cover in grassland is results of the Grain for Green policy (https://forestsnews.cifor.org/52964/grain-for-green-how-china-is-swapping-farmland-for-forest?fnl=en) and shrub encroachment into grassland, a process happening in many place of the world (D'Odorico et al., 2012).
Research showed that the shrubs in grassland have very similar greenness to the grass or forbs (Shimada et al., 2012). Hence, the replacement of grass and agriculture by shrubs may not cause significant change in greenness. The weak land-cover change effect on greenness in grassland is reflected in the ANOVA of the full model, in which agriculture cover, as the only land-cover type appeared in the fixed effects, contributed only 6% of total F-statistic. In addition, the temperate grassland in China lies in the semiarid and arid climate zone, and was found much more sensitive to rainfall increase than forests (Peng et al., 2012). All these factors might explain why the impact of climate dominates over that of land cover in the grassland zone.
The west China desert is also sensitive to rainfall. However, the insignificant change in rainfall does not explain the significant increased greenness (Figure 2). Bare ground in the peri-urban patches in desert zone significantly reduced and the shrubland and grassland expansion into the bare land might promote the greenness (Figure 2f). This is reflected in the ANOVA of the full model for the desert zone, which showed that bare land contributed to 75% of the total F-statistic for the fixed effects. Moreover, urban areas in the desert zone are mostly found within or near the oases. Increases in temperature accelerate the melting of glaciers in the mountains and provide more water for the oases which help vegetation growth and promote the greenness in the peri-urban areas.
The mixed forest zone has relative low land-cover change rate (Figure 2f). However, the land cover accounted for 48% of the full model SSD. The agriculture cover appeared in both fixed and random effects, and reclamation of wetland in some patches compensated the agriculture reduction in other 10.1029/2019JG005336

Journal of Geophysical Research: Biogeosciences
patches . The 27% average forest cover (Table 1) explains much of the greenness, which might lead to overestimation of the importance of land-cover change. Moreover, previous studies based on monthly climate data captured the changes in phenology in northeast China . The annual climate variables in this study may not be able to capture the phenological change that happened in this zone, thus might underestimate the climate impact on vegetation greenness.

Importance of Climate Variability
The impact of climate variability can be much stronger than mean climate change, and an extreme event may damage or remove the aboveground organs of vegetation. Forest ecosystem function was found to respond to climate variation (Goulden et al., 1996). Drought-induced defoliation and die-off of trees have been observed worldwide (Camarero et al., 2015). The mechanism of the die-off has been found primarily due to the hydraulic failure of xylem rather than the carbon starving (Saiki et al., 2017). Even without defoliation, tropical dry forest is much less green in dry season than wet season (Chapin et al., 2011).
The period of our study includes one very strong ENSO event (2015) and two strong La Niña events of 2007-2008 and 2010-2011, thus featuring great intraannual and interannual climate variability. We found that intraannual climate variability of the precipitation and the reference evapotranspiration negatively impacted annual mean greenness of the patches in the three forest zones in east China. Greater variability in precipitation and evapotranspiration means more concentrated seasonal profiles of the moisture input and output of ecosystems. More concentrated precipitation implies more runoff-loss in the period with ample or excessive precipitation and prolonged droughts in the period with less rainfall. Concentrated evapotranspiration tends to create great deficit in moisture budget and may cause defoliation and decreases in greenness. Aside from the water-budget related variability, the seasonality of minimum temperature also negatively impacted the greenness of the patches in grassland zone. The increased S Tmin makes the winter colder or the summer hotter. Hotter summer implies higher maintenance respiration to bring down the growth of leaves. In contrast to the negative impact of the abovementioned variabilities, the variability of the maximum temperature positively impacted the greenness of patches in evergreen forest zone (Table 3). An increase in intraannual variability of the maximum temperature for the evergreen forest zone may help summer photosynthesis to boost the leaf growth. The strength of the impact of the climate variability, measured as the percent of SSD referenced to the climate impact (Figure 4), peaked for patches in the temperate-cold-temperate mixed forest zone to explain 91% of the climate impact. The climate impact in Tibet is the strongest among the six zones (Figure 4), explaining 94% of SSD of the full model. However, the climate variability variables did not appear in the equation for patches in Tibet. One possible explanation is that increases in the mean minimum temperature is strong enough to explain the variation of greenness (Table 3 and Figure 4). The strength of the impact of the climate variability was also high for patches in the desert (68%) and the grassland (64%) zones in north and northwest China with cold weather (Table 1), but moderate for patches in the deciduous (52%) and the evergreen forest (48%) zones in east China. Hence, the contributions of climate variability in the regions with low temperature (mixed forest, grassland, and desert) are much greater than those in deciduous forest and evergreen forest zones. The reason why climate variability affects more the low-temperature mixed forest, grassland, and northwest desert zones is probably the plants in these regions live in the boundary of their niches and are more sensitive to intraannual variability than plants living in the middle of the niches.

Concluding Remarks
In this study, we applied linear mixed-effect model to investigate the impacts of climate and land-cover changes and climate variability on vegetation greenness in peri-urban patches around the meteorological stations in China. We found that changes in temperature variability of the meteorological stations are much stronger than the changes in the mean climate, and the changes in land cover of these patches feature great agriculture retreat, fast urban expansion, and moderate forest regrowth. We found that the impact of land cover is much stronger in patches located in the densely populated zones with fast economic growth. On the other hand, the impact of climate change dominates for patches within the high-cold Tibet and less densely populated temperate grassland. The variabilities of precipitation and reference evapotranspiration 10.1029/2019JG005336

Journal of Geophysical Research: Biogeosciences
significantly, negatively impacted the greenness in the patches of the mixed forest, deciduous forest, evergreen forest, grassland, and desert zones.
China will continue the urbanization in the next decades. Our results of the land-cover change in combination with the intensified climate variability will continue to increase their impacts on the services of the urban and peri-urban ecosystems.