Impact of Large-Scale A ff orestation on Surface Temperature : A Case Study in the Kubuqi Desert , Inner Mongolia Based on the WRF Model

Afforestation activities in the Kubuqi Desert, Inner Mongolia, China, have substantially increased tree and shrub coverage in this region. In this study, the response of the surface temperature to afforestation is simulated with the Weather Research and Forecasting model. The surface temperature changes are decomposed into contributions from the intrinsic surface biophysical effect and atmospheric feedback, using the theory of intrinsic biophysical mechanism. The effect of afforestation on the surface temperature is 1.34 K, −0.48 K, 2.09 K and 0.22 K for the summer daytime, the summer nighttime, the winter daytime and the winter nighttime, respectively, for the grid cells that have experienced conversion from bare soil to shrub. The corresponding domain mean values are 0.15 K, −0.2 K, 0.67 K, and 0.06 K. The seasonal variation of surface temperature change is mainly caused by changes in roughness and Bowen ratio. In the daytime, the surface temperature changes are dominated by the biophysical effect, with albedo change being the main biophysical factor. In the nighttime, the biophysical effect (mainly associated with roughness change) and the atmospheric feedback (mainly associated with change in the background air temperature) contribute similar amounts to the surface temperature changes. We conclude that the atmospheric feedback can amplify the influence of the surface biophysical effect, especially in the nighttime.


Introduction
Afforestation in semi-arid areas is an important method to protect soil, combat densification and improve the local environment.According to the Food and Agriculture Organization (FAO) of the United Nations, the area of planted forests was 264 Mha worldwide in 2010, and the rate of afforestation was 4.3 Mha/year during the decade beginning in 1990 and 4.9 Mha/year during the decade starting in 2000 [1].Of the global total afforested area, about 22% is in semi-arid climates [1].Several countries in arid regions reported a net afforestation between 1990 and 2010, including Egypt, Tunisia, Cyprus, Iraq, Israel, Kuwait, Kyrgyzstan, Lebanon, Syria, Turkey, and the United Arab Emirates.The afforestation activities in Turkey can be traced to 1870, and 2.06 Mha of afforestation was completed by 2015 [2].In the semiarid Mediterranean Basin, afforestation activities have covered more than 2.5 Mha during the second part of the 20th century [3].The area of afforestation in North Africa and West Asia is 4 Mha from 1950 to 2000, representing 6% of the actual forest and woodland in these regions [4].In China, multiple afforestation programs have taken place since 1978, resulting in the largest afforested area in the world (about 69 Mha in 2013, [5]).
The Kubuqi Desert is the seventh largest desert in China and is a major dust source that affects the air quality in Inner Mongolia and the Beijing municipality.To reduce desertification and to control dust storms, local governments and private sectors have been implementing large-scale environmental improvement projects (including planting trees and shrubs).A total of 0.6 Mha of bare land was replaced with vegetation by 2010 [6].The environmental improvement in the Kubuqi Desert is an example of successful ecological restoration in China and was featured by the United Nations Convention to Combat Desertification and the United Nations Environment Programme during the Fifth Kubuqi International Desert Forum [7].
The change of land-use changes the biophysical properties of the land surface (surface roughness, evaporation, albedo).It has been long known that the regional climate is sensitive to even small changes in land surface properties [8].The biophysical effects of land-use change to climate constitute their regulation of the exchanges of energy, water, and momentum between the earth's surface and the lower atmosphere [9][10][11].The biogeochemical effects refer to greenhouse gas forcing due to the changing CO 2 flux to the atmosphere and the possible carbon storage in plants [10,12].More recent studies on forest management have shown that the biophysical changes can result in similar or stronger local climate changes than the effects arising from biogeochemical processes [10,[13][14][15][16][17][18].
How the large-scale afforestation influences the local and regional climate through biophysical processes is a subject of several previous studies.Using a space-for-time approach that compares the surface temperatures of forests with those of adjacent grasslands and croplands, Peng, et al. [19] found that afforestation in China has a cooling effect in the daytime and a warming effect in the nighttime.Cao, et al. [20] assessed the impacts of land-use change in the agro-pastoral transitional zone of North China between 2001 and 2010 with the Weather Research and Forecasting (WRF, [21]) model.Their results show that the conversion of cropland to grassland, a significant land-use change in this region, resulted in local cooling in the summer and warming in the winter.In their study, the air temperature change of the Kubuqi desert area between 2000 and 2010 is not significant, ranging from −0.2 to 0.2 K in the summer and −0.4 to 0.4 K in the winter.Ge,et al. [22] investigated the effects of vegetation fraction change from 1982 to 2000 on the temperature in Eastern China with the WRF model.They concluded that climate warming was slowed in the regions with increased vegetation and was enhanced in the regions with decreased vegetation.Using the National Center for Atmospheric Research (NCAR) second generation regional climate model (RegCM2) [23,24], Liu, et al. [25] found that the afforestation projects in northern China increase precipitation by 15% and reduce the air temperature by nearly 0.5 K. Some other studies also investigate the effects of afforestation on soil quality and precipitation (e.g., Zhang, et al. [26] and Kuriqi [27]).In the present study, we focus on the surface temperature change and the related energy flux changes.
The land surface temperature T s is an important climatic variable for understanding the interactions between land surface changes and atmospheric changes.This is because the T s is a sensitive indicator of perturbations brought by land-use changes to the surface energy balance.It is a quantity routinely monitored by remote sensing and calculated by land surface models.To generate mechanistic insights into the land-atmosphere interactions, the total T s changes brought on by afforestation and other Forests 2019, 10, 368 3 of 21 land-use activities are often partitioned into component contributions associated with individual biophysical factors.One theory used for the partitioning calculation is the intrinsic biophysical mechanism (IBPM) theory proposed by Lee, et al. [28].According to the IBPM theory, changes in the surface temperature are a consequence of a local longwave radiative feedback on perturbations of the surface albedo and energy redistribution caused by changes to the aerodynamic resistance (or surface roughness) and surface evaporation.This theory has been used in observational studies involving pairs of sites [13,28] and in climate modeling studies comparing subgrid T s variations [29][30][31].
An assumption underlying these applications is that the background atmosphere, as characterized by atmospheric forcing variables at the blending height, mainly the incoming shortwave radiation, the incoming longwave radiation and the air temperature, is unaffected by changes in the surface biophysical properties (Figure 1a,b).intrinsic biophysical mechanism (IBPM) theory proposed by Lee, et al. [28].According to the IBPM 99 theory, changes in the surface temperature are a consequence of a local longwave radiative feedback on perturbations of the surface albedo and energy redistribution caused by changes to the aerodynamic resistance (or surface roughness) and surface evaporation.This theory has been used in observational studies involving pairs of sites [13,28] and in climate modeling studies comparing subgrid Ts variations [29][30][31].An assumption underlying these applications is that the background atmosphere, as characterized by atmospheric forcing variables at the blending height, mainly the incoming shortwave radiation, the incoming longwave radiation and the air temperature, is unaffected by changes in the surface biophysical properties (Figure 1 a & b).In the real atmosphere however, feedback processes can amplify or diminish the biophysical effect.Possible atmosphere feedback mechanisms include changes of lapse rate, water vapor content, cloud amount and cloud altitude [32,33].In the context of the surface energy balance, these mechanisms are manifested in changes in the incoming longwave radiation (ΔL↓) and shortwave radiation (ΔK ↓ ) and the background air temperature (ΔTa).Green, et al. [34] found that the atmospheric feedback explains 30% of surface radiation variations in the Mediterranean region, where the enhancement of latent and sensible heat transfers due to denser vegetation increases the boundary layer height and convection, affecting cloudiness, and consequently the incident surface radiation.He,et al. [35] found that the atmospheric feedback caused by shrub encroachment creates a warmer microclimate in the northern Chihuahuan desert, reducing juvenile shrubs' vulnerability to freeze-induced mortality, which in turn favors shrub growth.The total Ts change should consist of a component arising from the intrinsic surface biophysical change and a component that accounts for the atmospheric feedback (Figure 1c).Splitting the forcing results into contributions of the intrinsic biophysical effect and the atmospheric feedback can determine the nature of the both terms.A deeper understanding of the relative importance of the two terms is required in the design of afforestation projects with limited resources.Explicit measures of each term can inform decisions for this afforestation and the associated land-use management.Policymakers can benefit from such studies to choose appropriate climate-effective mitigation actions.In the real atmosphere however, feedback processes can amplify or diminish the biophysical effect.Possible atmosphere feedback mechanisms include changes of lapse rate, water vapor content, cloud amount and cloud altitude [32,33].In the context of the surface energy balance, these mechanisms are manifested in changes in the incoming longwave radiation (∆L ↓ ) and shortwave radiation (∆K ↓ ) and the background air temperature (∆T a ).Green, et al. [34] found that the atmospheric feedback explains 30% of surface radiation variations in the Mediterranean region, where the enhancement of latent and sensible heat transfers due to denser vegetation increases the boundary layer height and convection, affecting cloudiness, and consequently the incident surface radiation.He, et al. [35] found that the atmospheric feedback caused by shrub encroachment creates a warmer microclimate in the northern Chihuahuan desert, reducing juvenile shrubs' vulnerability to freeze-induced mortality, which in turn favors shrub growth.The total T s change should consist of a component arising from the intrinsic surface biophysical change and a component that accounts for the atmospheric feedback (Figure 1c).Splitting the forcing results into contributions of the intrinsic biophysical effect and the atmospheric feedback can determine the nature of the both terms.A deeper understanding of the relative importance of the two terms is required in the design of afforestation projects with limited resources.Explicit measures of each term can inform decisions for this afforestation and the associated land-use management.Policymakers can benefit from such studies to choose appropriate climate-effective mitigation actions.
The published studies cited above have given an overall conclusion about how land-use change may affect the local temperature.To our best knowledge, no studies have isolated the roles of Forests 2019, 10, 368 biophysical forcing and atmosphere feedback, nor have they quantified the dominant contributors to the regional temperature change.In this study, the biophysical effect and atmospheric feedback of afforestation in the Kubuqi Desert are evaluated with the WRF model and the IBPM theory.Here, afforestation is defined as planting trees and shrubs in bare land.Our focus is on the surface temperature response to afforestation in this region.Specifically, we aim: (1) to evaluate the seasonal and diurnal surface temperature change caused by land-use and land-cover change in the Kubuqi Desert, and (2) to compare the relative contributions of the intrinsic biophysical effect and the atmospheric feedback to surface temperature changes.
This study extends the work we have recently published [36].In that study, we deployed a space-for-time strategy to investigate the T s response to afforestation.Specifically, we used the IBPM theory to decompose the T s difference observed between an afforested site and a shrub site in the Kubuqi Desert into component contributions.In the present study, we use the WRF model to quantify the T s change at the regional scale.The flux observations described in Wang, et al. [36] are used to verify the WRF model performance.An advantage of the WRF modeling over the pair-sites analysis is that we can quantify the effect of not only shrub-to-forest conversion but also other types of land-use change in the Kubuqi Desert.By coupling the atmosphere dynamically with the land surface in the model domain, the model allows us to investigate the role of atmospheric feedback via changes in the incoming shortwave radiation, the incoming longwave radiation and the background air temperature.In the space-for-time analysis, these quantities are assumed the same between the paired sites, but the assumption does not necessarily hold when the background atmosphere is allowed to respond to land-use changes.In the present study, we attempted to address whether the atmospheric feedback processes can amplify (or diminish) the biophysical effect using WRF simulation.

WRF Modeling
The WRF model is a mesoscale atmospheric model widely used in applications ranging from regional climate prediction to the investigation of land-atmosphere interactions (https://www.mmm.ucar.edu/weather-research-and-forecasting-model).In this study, we used the WRF dynamics solver described by Skamarock, et al. [21].The NOAH land surface model [37,38] was used to simulate the surface fluxes to the atmosphere.The other main physical parameterizations used for the simulations are presented in Table 1.Monin-Obukhov [40,41] a CAM3, the Community Atmosphere Model Version 3 for longwave and shortwave radiation [42,43]; b WSM3, the WRF Single-Moment 3 class microphysics scheme [44]; c YSU, the Yonsei University planetary boundary layer scheme [45,46].
The WRF model is configured with two nested domains (Figure 2).The coarse outer domain (D01) covers a total area of 1106 km × 966 km with a 14-km grid spacing in both horizontal directions.The inner domain (D02) is centered at 40.38 • N and 108.25 • E and covers the entire Kubuqi Desert (238 km × 196 km), using a 2 km grid spacing.Considering the computing time consuming and to keep the model away from collapse due to high resolution of the outer domain, a parent grid ratio of 7 was used to achieve a fine resolution of the inner domain.This ratio is higher than the recommended ratios of 3 and 5 [21].It may cause the nonequilibrium of weather and turbulence with grid spacing in the border grid cells [47].To avoid bias, seven layer grid cells on the boundaries were excluded in the statistical analysis.The projection is Lambert with two standard parallels of 30 • N and 60 • N and the central meridian of 109 • E. Two one-year continuous simulations, labeled as EXP_2000 and EXP_2008, were performed using identical settings except for the underlying land-use and land-cover data.Years 2000 and 2008 were chosen because the data from meteorological stations and eddy flux tower were available during these periods.In EXP_2000 the land-cover is that in the year 2000 and EXP_2008 is in the year 2008.Both EXP_2000 and EXP_2008 use the same initial and lateral boundary conditions in 2008.The experimental design of the two simulations is shown in Figure S1 as a flow chart.The integration time spans 0000 UTC on 1 January 2008 to 1800 UTC on 31 December 2008 for both EXP_2000 and EXP_2008, and the output time interval is 1 h.a CAM3, the Community Atmosphere Model Version 3 for longwave and shortwave radiation [42,43]; 167 b WSM3, the WRF Single-Moment 3 class microphysics scheme [44]; c YSU, the Yonsei University 168 planetary boundary layer scheme [45,46].

169
The WRF model is configured with two nested domains (Figure 2).The coarse outer domain 170 (D01) covers a total area of 1106 km × 966 km with a 14-km grid spacing in both horizontal directions.171 The inner domain (D02) is centered at 40.38°N and 108.25°E and covers the entire Kubuqi Desert (238 172 km × 196 km), using a 2 km grid spacing.Considering the computing time consuming and to keep 173 the model away from collapse due to high resolution of the outer domain, a parent grid ratio of 7 174 was used to achieve a fine resolution of the inner domain.This ratio is higher than the recommended 175 ratios of 3 and 5 [21].It may cause the nonequilibrium of weather and turbulence with grid spacing 176 in the border grid cells [47].To avoid bias, seven layer grid cells on the boundaries were excluded in  2).The area of bare land decreased from 15,711 km 2 in 2000 to 6176 km 2 in 2008, or about 20% of the study region.The reduction in the bare land area was mainly caused by conversion to shrubland.The spatial distribution of the transition zone matches the aerial seeding operations prior to 2008 [6].The dominant original shrub in Kubuqi desert is Artemisia Ordosica Krasch., which is also chosen to plant through aerial seeding.The grasses in this area mainly include licorice (Glycyrrhiza uralensis Fisch.) and alfalfa (Medicago Sativa L.).The principal crops are maize and sunflower (Helianthus annuus L.).

6
Krasch., which is also chosen to plant through aerial seeding.The grasses in this area mainly include licorice (Glycyrrhiza uralensis Fisch.) and alfalfa (Medicago Sativa L.).The principal crops are maize and sunflower (Helianthus annuus L.).Table 2. Comparison of biophysical effect and dynamic feedback for different land-use/land-cover transition and the contributions to changes in surface temperature from individual factors (Terms 1-7 in K).K ↓ is incoming shortwave radiation (W m −2 ), L ↓ is incoming longwave radiation (W m −2 ), T a is air temperature (K) and the change ∆ denotes the difference between EXP_2008 and EXP_2000.∆T 1 denotes the surface temperature change caused by intrinsic biophysical effect and ∆T 2 denotes the surface temperature change caused by atmospheric feedback.The land-use and land-cover data are primarily based on the Landsat Thematic Mapper digital images acquired in 2000 and 2008.Random Forest algorithm is used in the supervised classification [48,49].We upscale the classification results from the 30-m satellite pixel resolution to the 2-km resolution by assigning the land-use to the dominant type in each WRF grid cell.

Land
The initial and lateral boundary conditions for both simulations were provided by the weather forecast model Global Forecast System (ftp://nomads.ncdc.noaa.gov/GFS/analysis_only)product for 2008 with a 0.5 • by 0.5 • spatial resolution at 6-hourly intervals.The observation data used for evaluating the model performance were obtained from the National Meteorological Information Center, Chinese Meteorological Administration (CMA [50], http://data.cma.gov.cn/).These datasets contain the daily climatological variables (the pressure, the temperature at 1.5 m, the wind at 10 m and the precipitation) for 2008 for nine meteorological stations across the Kubuqi Desert (Figure 3c, black dots).In addition, the surface flux variables (sensible heat, latent heat, ground heat) were evaluated against the data obtained in 2008 at an eddy-covariance site located in a shrubland ecosystem in the center of the Kubuqi Desert (40 • 32 N, 108 • 41 E, Figure 3c, red dot; [36]).

Offline Calculations
According to the IBPM theory [28], the surface temperature is expressed as: where T s is surface temperature (K), T a is air temperature (K) at the blending height, ) is given by: where σ is the Stephan-Boltzmann constant, f is an energy redistribution factor (dimensionless), and R * n is apparent net radiation: where K ↓ is incoming shortwave radiation (W m −2 ), α is surface albedo, and L ↓ is incoming longwave radiation (W m −2 ), the proof of Equation ( 1) is given in the supporting information.The energy redistribution factor f is given by: where ρ is air density, C p is specific heat of air at constant pressure, β is Bowen ratio, and r t is total heat transfer resistance.In a site pair analysis, the surface temperature perturbation ∆T s is found by differentiating Equation (1) by holding K ↓ , L ↓ and T a constant between the paired sites.If we allow K ↓ , L ↓ and T a to change in response to land-use changes, we obtain a more general expression for ∆T s , as where ∆ denotes the difference in a variable between EXP_2008 and EXP_2000 (value in EXP_2008 minus value in EXP_2000).For example, ∆α is the difference in albedo between EXP_2008 and .The change in the energy redistribution factor consists of changes in surface roughness (∆f 1 ) and Bowen ratio (∆f 2 ): In Equation ( 5 In this study, the first four terms represent the intrinsic biophysical effect, and the last three terms represent the dynamic feedback.One specific objective of this study is to compare the total biophysical effect (sum of Terms 1 to 4) and the total effect arising from the dynamic feedback (sum of Terms 5 to 7).
Each term in Equation ( 5) is calculated offline using the variables saved from the online WRF model calculation.The term "offline" means that we calculated each term in Equation ( 5) after the WRF simulation using the model outputs.The "offline" ∆T s represents the summation of all terms in Equation ( 5) and the "online" ∆T s represents the change of surface temperature simulated in the WRF model.Two intermediate variables, the Bowen ratio and the total resistance to heat diffusion are computed as where H and LE are sensible and latent heat flux, respectively.The IBPM calculation is applied to hourly data in each grid.The hourly results are averaged according to daytime and nighttime periods.When the incoming shortwave radiation is larger than 0 W/m 2 , the results are considered as daytime; otherwise, they are considered as nighttime.The analysis is done separately for daytime and nighttime periods to account for the diurnal asymmetry of the biophysical effect [51].

Model Evaluation
Model evaluation was carried out as shown in Figure 4 (monthly mean) and Figure 5 (hourly scale).A general agreement is found between the monthly mean air temperatures (2 m level result) in the EXP_2008 simulation and the observations at the nine meteorological stations within the model domain (Figure 4), with the R 2 value of around 0.99 and a mean bias of −0.74 • C.This result is consistent with the WRF simulation by Zhang, et al. [52] which showed a cold bias between 1.0 • C and 1.5 • C for the whole of China.The root mean square error (RMSE) is 1.76 • C, slightly better than the result (1.92 • C) obtained by Cao, et al. [20].In the summer, the RMSEs, ranging from 0.94 (June) to 1.33 • C (July), are smaller than the results of Zhang, et al. [52] (ranging from 1 to 2.5 • C).In the winter, the RMSEs range from 1.64 (January) to 2.39 • C (November).According to Zhang, et al. [52], the RMSEs range from 2-3 • C in the winter.The simulated results from 14 February to 24 February and 4 August to 14 August 2008 were extracted from the nearest grid to the eddy covariance tower and compared with the observations.The results show that the simulated radiation and energy components (incoming shortwave radiation, incoming longwave radiation, sensible heat flux, latent heat flux and soil heat flux) are in agreement with those from observations (Figure 5).In the summer, the mean bias is The simulated sensible heat flux is close to the observation in the clear sky condition, but it is higher than the observation on 8 and 14 August.These two days are cloudy days.This overestimation may result from a random time lag in the radiation and sensible heat flux.However, this time lag is not systematic but highly depends on atmosphere condition.For example, the model successfully simulates the radiation and energy fluxes on another cloudy day (9 August).
Forests 2019, 10, x FOR PEER REVIEW 11 of 22 24th and August 4th to August 14th,2008 were extracted from the nearest grid to the eddy covariance tower and compared with the observations.The results show that the simulated radiation and energy components (incoming shortwave radiation, incoming longwave radiation, sensible heat flux, latent heat flux and soil heat flux) are in agreement with those from observations (Figure 5).In the summer, the mean bias is 2.4 W m −2 , −12.8 W m −2 , 2.5 W m −2 , 63.9 W m −2 and −14.5 W m −2 for sensible heat flux, latent heat flux, ground heat flux, incoming short wave radiation and incoming longwave radiation, respectively.The simulated sensible heat flux is close to the observation in the clear sky condition, but it is higher than the observation on August 8th and 14th.These two days are cloudy days.This overestimation may result from a random time lag in the radiation and sensible heat flux.However, this time lag is not systematic but highly depends on atmosphere condition.For example, the model successfully simulates the radiation and energy fluxes on another cloudy day (August 9th).

Online Versus Offline Surface Temperature Change
The spatial distributions of ΔTs calculated online and offline for the summer daytime (Figure 6a   and b) and nighttime (Figure 6e and f) were compared.The results for the winter season are given in Figure 7.The online results express ΔTs as the difference in the surface temperatures calculated by the EXP_2008 and the EXP_2000 model simulations.The offline ΔTs is the sum of the seven individual terms in Equation ( 5) calculated with the offline diagnostic variables.The spatial distributions of the online ΔTs and the offline ΔTs match each other, for both the daytime and the nighttime period and in both seasons.In the daytime, for the summer and the winter, both the online result and offline result show a warming effect, especially for the afforested area along the Yellow River and in the southwestern portion of the domain (Figure 6a&b and Figure 7a&b).The surface temperature does not change in the core area of the desert (in the middle of the domain).North of the Yellow River, both online and offline results show a fragmented pattern of warming and cooling effects.For the summer nighttime period, most of the southwestern area exhibiting land-cover change shows a cooling trend (Figure 6e), while in the winter night, the cooling trend is not clearly (Figure 7e).The differences between online and offline ΔTs can be found in Figure S2.The RMSE is 0.35 K, 0.07 K, 0.17 K and 0.07 K for the summer daytime, the summer nighttime, the winter daytime, and the winter nighttime, respectively.

Online versus Offline Surface Temperature Change
The spatial distributions of ∆T s calculated online and offline for the summer daytime (Figure 6a,b) and nighttime (Figure 6e,f) were compared.The results for the winter season are given in Figure 7.The online results express ∆T s as the difference in the surface temperatures calculated by the EXP_2008 and the EXP_2000 model simulations.The offline ∆T s is the sum of the seven individual terms in Equation ( 5) calculated with the offline diagnostic variables.The spatial distributions of the online ∆T s and the offline ∆T s match each other, for both the daytime and the nighttime period and in both seasons.In the daytime, for the summer and the winter, both the online result and offline result show a warming effect, especially for the afforested area along the Yellow River and in the southwestern portion of the domain (Figures 6a,b and 7a,b).The surface temperature does not change in the core area of the desert (in the middle of the domain).North of the Yellow River, both online and offline results show a fragmented pattern of warming and cooling effects.For the summer nighttime period, most of the southwestern area exhibiting land-cover change shows a cooling trend (Figure 6e), while in the winter night, the cooling trend is not clearly (Figure 7e).The differences between online and offline ∆T s can be found in Figure S2.The RMSE is 0.35 K, 0.07 K, 0.17 K and 0.07 K for the summer daytime, the summer nighttime, the winter daytime, and the winter nighttime, respectively.Figure 8 and Figure 9 present the comparison between online and offline results as scatter plots (Figure 8a for the summer daytime, Figure 8d for the summer nighttime, Figure 9a for the winter daytime and Figure 9d for the winter nighttime), with each data point denoting one grid point.The results follow the 1:1 line with high determination coefficients (R 2 > 0.9) except for the summer daytime (R 2 = 0.78).Figure 8 and Figure 9 present the comparison between online and offline results as scatter plots (Figure 8a for the summer daytime, Figure 8d for the summer nighttime, Figure 9a for the winter daytime and Figure 9d for the winter nighttime), with each data point denoting one grid point.The results follow the 1:1 line with high determination coefficients (R 2 > 0.9) except for the summer daytime (R 2 = 0.78).Figures 8 and 9 present the comparison between online and offline results as scatter plots (Figure 8a for the summer daytime, Figure 8d for the summer nighttime, Figure 9a for the winter daytime and Figure 9d for the winter nighttime), with each data point denoting one grid point.The results follow the 1:1 line with high determination coefficients (R 2 > 0.9) except for the summer daytime (R 2 = 0.78).

Seasonal Differences
Generally, the replacement of bare land with shrubland increases the temperature in the daytime both in the winter and in the summer.In the summer, the land-use change from 2000 to 2008 has a slight warming effect of about 0.06 K for the whole domain.If we focus on the afforestation grids where the replacement of bare land with shrubland had occurred, the warming effect is 0.43 K.In the winter, the average warming effect is 0.24 K for the whole domain and 1.12 K for the afforestation area.The calculated results in this part are mean values of the daytime results and the nighttime results, so they are not the same as those in Table 2.
The attribution analysis reveals that the seasonal variation is mainly caused by roughness change (Term 2) and Bowen ratio change (Term 3, Table 1, Figures S3-S6, partition results for the summer daytime, summer nighttime, winter daytime and winter nighttime, respectively).Take the afforestation area as an example: In the summer, Term 2 is negative (about −0.40 K), and the effect is negligible in the winter.As for Term 3, the effect is negative (about −0.08 K) in the summer, which offsets part of the total warming effect.In the winter, Term 3 is positive (about 0.20 K), which enhances the overall warming effect.The seasonal difference caused by these two terms is 0.68 K, which is almost the same as the total difference (0.75 K) for the afforestation grid cells.Besides, the soil heat flux changes (Term 4) have an effect about -0.02 K in the summer an 0.02 K in the winter across the model domain.The absolute values of Term 4 are higher in the summer (both daytime and nighttime) than those in the winter, mainly due to the seasonal variations of soil moisture.

Daytime versus Nighttime
When evaluating the influence of land-use change, it is necessary to consider daytime and nighttime periods separately because diurnal asymmetry exists in both the magnitude and the sign of ∆T s [28,51].In the summer daytime, the surface temperature of the whole domain increases by an average of 0.15 K from 2000 to 2008 (Figure 6a, Table 2).Conversion of bare land to shrubland causes a summer daytime warming effect of 1.34 K, with the main contributor of albedo change (1.74 K, Table 2, Figure S3).In the winter daytime, the surface temperature of the whole domain increases by 0.67 K, and the grid cells of bare land to shrubland conversion (the afforestation area) warms by an average of 2.09 K.Both the albedo change (1.11 K, Term 1), the Bowen ratio change (0.37 K, Term 3) and the air temperature change (0.57K, Term 7) contribute to the winter daytime warming effect (Table 2, Figure S5).
In the summer night, the whole domain cools by an average of 0.20 K (average ∆T s = −0.20 K).The afforestation grid cells have a substantial nighttime cooling effect (average ∆T s = −0.48K), with the main contributor of roughness change in these grid cells (−0.37 K, Table 2, Figure S4).In these grid cells, the background air temperature change also makes a significant contribution (−0.24K) to the nighttime surface cooling.
In the winter nighttime, the surface temperature of the whole domain shows a negligible change (∆T s = 0.06 K).The grid cells of bare land to shrubland conversion experience a slight warming effect (∆T s = 0.22 K) mainly due to roughness change (0.07 K, Figure S6) and air temperature change (0.09 K).We found that the diurnal asymmetry, the phenomenon in which the daytime ∆T s has the opposite sign with the nighttime ∆T s , is evident in the summer (Figure 6a,e), but not in the winter (Figure 7a,e).
For those grid cells where the land-use did not change between 2000 to 2008, the surface temperature change is minimal (Table 2), ranging from −0.12 K in the summer night to 0.2 K in the winter daytime.As the licorice planting is promoted in local plantation, three other kinds of land-use change, shrubland (Artemisia Ordosica) to grassland (licorice and alfalfa), cropland (maize and sunflower) to grassland and shrubland to cropland, account for 5%, 2% and 1% of the total area.The surface temperature change is minor for grid cells that experienced conversion of cropland to grassland, with the absolute values of ∆T s less than 0.25 K.As for those cells experienced conversion of shrubland to grassland and shrubland to cropland, the daytime ∆T s ranges from 0.86 K to 0.99 K in the winter, and −0.44 K to −0.7 K in the summer, and the nighttime ∆T s are smaller, ranging from −0.35 K to 0.47 K.A detailed information about these partitioning results and whole domain averaged results can be found in Table 2.

Warming Effect in the Daytime
As shown in Figures 6 and 7, we found the replacement of bare land with shrubland increases the temperature in the daytime, and the warming effect is stronger in the winter than in the summer.Peng, et al. [19] used the MODIS satellite product to compare the surface temperature of forests with those of adjacent croplands and grasslands.He found that afforestation in China has a cooling effect (−1.1 ± 0.6 K) in the daytime for those areas with mean annual precipitation between 400 and 600 mm y −1 .Although the forest absorbs more incoming radiation because of lower albedo, it dissipates even more energy through evaporative cooling.The results in Peng, et al. [19] are consistent with our previous study [36], which also shows that the planted forest is cooler than the adjacent shrubland (−0.5 ± 0.2 K) in the daytime.However, in the present study, the major land-use change on the regional scale is the replacement of bare land with shrubland.Thus, the warming signal reflects the difference between shrubland and bare land instead of the difference between forest and other vegetation types (cropland, grassland or shrubland).In Peng, et al. [19], the albedo of the forest is 2.8% lower than that of cropland and 1.7% lower than that of grassland.Forest has deeper root and it can easily acquire the groundwater, resulting to higher evapotranspiration.The energy loss through evapotranspiration cools the surface because it easily exceeds the extra absorbed solar radiation.Besides, the forest can increase the soil moisture and then increases the soil heat capacity, which can also cool the surface [19].While in this study, the albedo of shrubland is 12% lower than that of bare land, and the contribution of extra absorbed energy obviously exceeds the contribution of other terms which can cause a cooling effect (Figure 6a, Figure 7a).The sign and magnitude of the temperature change here are different with those in Wang, et al. [36] mainly due to the different albedo changes in the two processes.In the present study, the albedo of shrubland is lower than that of bare land (12%).While in Wang, et al. [36], the albedo of the planted forest is 30% higher than that of shrubland.The difference in albedo change, which can explain the contradiction of similar studies, is sensitive to the type of land-use change.
The Kubuqi Desert is a typical Northern desert belonging to the temperate continental climate.But in the tropical climate, the warming effect is often reported in the deforestation processes.For example, Abera, et al. [53] found the net warming effect is 1.3 K when the land cover converted from forest to cropland.Although the albedo change decrease the land surface temperature by 0.12 K, the dominant influence of non-radiative mechanisms (roughness change and evapotranspiration change) is about 10 times higher than that of radiative forcing.The study site of Abera, et al. [53] is the Horn of Africa, which is close to the Indian Ocean.The precipitation in the forest area can reach 1500 mm, about five times higher than that in the Kubuqi Desert.The higher soil moisture and evapotranspiration keep the forest cooler than the cropland, while in our study, the contribution of evapotranspiration is small due to the water deficit and the dominant influence is radiative forcing.

Intrinsic Biophysical Effect versus Atmospheric Feedback
In this section, we focus on the grid cells of bare land to shrubland conversion as these cells have experienced the largest changes in surface temperature in comparison to other grid cells.In the daytime, the contribution of the surface biophysical changes is larger than that of the atmospheric feedback, accounting for nearly 86% and 98% of the total ∆T s in the summer (Figure 8b) and in the winter, respectively (Figure 9b).
Specifically, the biophysical effect is dominated by the albedo change (Term 1, Equation ( 5)), which contributes 1.74 K to the summer daytime temperature increase and 1.11 K to the winter daytime temperature increase.The afforested area has a lower albedo (0.22) than bare land (0.25) in the model (default values in WRF version 3.7.1),and therefore absorbs more solar energy and warms up the surface in the daytime.The contribution of roughness change (Term 2) is −0.44 K in the summer, which offsets part of the albedo warming effect.The mechanism underlying this cooling is well documented: trees and shrubs enhance the convection between the land surface and the atmosphere, allowing more efficient heat dissipation via turbulent diffusion than bare land [54].In the winter, the influence of Term 2 is relatively small (−0.05K), because the shrub is dormant and the structure above ground is less than that in the summer.The shrubland has high canopy conductance and transpiration rate, which result in increased partitioning of net radiation into the latent heat flux.Thus, the Bowen ratio change (Term 3) has a cooling effect (−0.16 K) in the summer.
The conversion of bare land to shrubs and forests causes changes in the background atmosphere.The enhanced evapotranspiration changed the composition and properties of the clouds, and the result is that the incoming shortwave radiation (K ↓ ) decreases by −9.31 W m −2 (daytime mean) in the summer and −3.49W m −2 in the winter.Resulting from the cloud changes, the incoming longwave radiation (L ↓ ) increases by 6.83 W m −2 (daytime mean) in the summer and 2.56 W m −2 in the winter.But the changes in K ↓ and L ↓ do not contribute much to the surface temperature change since they counteract each other.Afforestation causes the domain-averaged surface incoming shortwave radiation to decrease by 1.9 W m −2 (daytime mean of summer and winter seasons).The land conversion causes the air temperature at the blending height to increase in both summer (0.33 K) and winter (0.63 K).The signs of background air temperature changes are consistent with the signs of surface temperature changes and it is the dominant atmospheric feedback process (Term 7; Table 2, Figures S3 and S5).We infer that the perturbation of boundary layer property which can cause air temperature change is sensitive to the land use types.
In the nighttime, the influence of atmospheric feedback is comparable with that of the biophysical effect (Figures 6e-h and 7e-h).In the summer, the total contribution of atmospheric feedback is −0.22 K, whereas the total biophysical effect is −0.26K.The corresponding values for the winter nighttime are 0.1 K and 0.12 K. Specifically, the atmospheric feedback is manifested mainly in the background air temperature change (Term 7), contributing −0.24 K and 0.09 K to the summer and winter nighttime surface temperature change, respectively.The main biophysical effect is roughness change (Term 2), contributing −0.37 K and 0.07 K to the summer and winter nighttime surface temperature changes, respectively.Afforestation causes the incoming longwave radiation (L ↓ ) to increase by 2.02 W m −2 (nighttime mean value) in the summer and 0.92 W m −2 in the winter, which causes a slight warming effect (0.02 K in the summer and 0.01 in the winter).As for the whole domain, the incoming longwave radiation increases by 0.3 W m −2 (nighttime mean of summer and winter seasons).
Figures 8a-c and 9a-c also show that the offline ∆T s in the daytime is dominant by biophysical effect (∆T 1 ).Most of the feedback effect (∆T 2 ) scatters in the range from −0.5 K to 0.5 K, or about 1/4 of the range of ∆T 1 .When the regression slope is closer to 1 (Figures 8b,c,e,f and 9b,c,e,f), the corresponding factor contributes more to the ∆T s .In the summer, the regression slope between ∆T 1 and online daytime ∆T s is 0.86 (R 2 = 0.70), but the slope value between ∆T 2 and online daytime ∆T s is 0.24 (R 2 = 0.58).The slope can be considered as an estimated fraction of the ∆T s caused by the biophysical effect/atmospheric feedback.It means the biophysical effect contribute much more to the ∆T s .In the winter, most of the variation in online daytime ∆T s can be explained by ∆T 1 ; the regression slope is 0.98, which is very close to 1, with a higher R 2 (0.90), indicating the biophysical effect is still the dominant factor.However, the regression slope between ∆T 2 and ∆T s is 0.28 (R 2 = 0.65) in the winter daytime.
In the nighttime, both ∆T 1 and ∆T 2 show a significant linear relationship with online ∆T s .The regression slopes in the summer are 0.57 (∆T 1 ) and 0.51 (∆T 2 ) with high R 2 values (0.95-0.98).The values in the winter are 0.54 (∆T 1 , R 2 = 0.99) and 0.47 (∆T 2 , R 2 = 0.97).All the slope values are close to 0.5, which means that the contribution of biophysical effect is about equal to that of atmospheric feedback.
We found the atmosphere feedback effects within the grid cells of afforestation are also not the same.The frequency histograms of atmosphere feedback effects were shown in Figure 10, and the frequency represents the number of grid cells.For example, the effects of atmosphere feedback range from 0 to 0.7 K with an average value of 0.37 K (∆T 2 , Table 2) in the summer daytime.The difference within the same land-use change type is related to the upwind land cover changed and the scale of afforestation.The prevailing wind in the Kubuqi Desert is from the northwest.The large scale afforestation along the yellow river and afforestation in the southwestern part of the domain might have an impact on the downwind region.The atmosphere feedback effect in the core desert area is about 0.2 K higher than the result in the unchanged desert area in the eastern part of the domain (Figure 6d) caused by the impact of afforestation along the yellow river in the upwind region.Moreover, the atmosphere feedback effect of the multiple adjacent conversion grid cells is higher than that of the isolated conversion grid cells.For example, the atmosphere feedback effect in the major afforestation area (southwestern part of the domain) is 0.6 K, while the result in the sparse afforestation area (40.13N, 109.17E, Figure 3c) is 0.3 K in the summer daytime.In the summer nighttime, the atmosphere feedback in the southwestern part of the domain contributes a cooling signal of −0.5 K, while the result in the sparse afforestation area is −0.1 K (Figure 6h).These results indicate that the atmospheric feedback in the downwind region can be influenced by the upwind area through horizontal transmission, and large scale conversion area has greater potential to change the atmosphere than isolated conversion area.
Forests 2019, 10, x FOR PEER REVIEW 18 of 22 18 within the same land-use change type is related to the upwind land cover changed and the scale of afforestation.The prevailing wind in the Kubuqi Desert is from the northwest.The large scale afforestation along the yellow river and afforestation in the southwestern part of the domain might have an impact on the downwind region.The atmosphere feedback effect in the core desert area is about 0.2 K higher than the result in the unchanged desert area in the eastern part of the domain (Figure 6d) caused by the impact of afforestation along the yellow river in the upwind region.Moreover, the atmosphere feedback effect of the multiple adjacent conversion grid cells is higher than that of the isolated conversion grid cells.For example, the atmosphere feedback effect in the major afforestation area (southwestern part of the domain) is 0.6 K, while the result in the sparse afforestation area (40.13 N, 109.17E, Figure 3c) is 0.3 K in the summer daytime.In the summer nighttime, the atmosphere feedback in the southwestern part of the domain contributes a cooling signal of −0.5 K, while the result in the sparse afforestation area is −0.1 K (Figure 6h).These results indicate that the atmospheric feedback in the downwind region can be influenced by the upwind area through horizontal transmission, and large scale conversion area has greater potential to change the atmosphere than isolated conversion area.

Limitations and Future Work
We found that for the 1:1 comparison of the online and offline ΔTs, the daytime results (Figures 8a & 9a) are more scattered than the nighttime results (Figures 8d & 9d) for both seasons.The derivation of Eq (1) assumes that nonlinear interactions in the change terms can be ignored.For example, concerning the changes in the net shortwave radiation, we assumed: where subscripts 1 and 2 represent EXP_2008 and EXP_2000, respectively, K↓ is the mean value of K ↓1 and K↓2, and α is the mean value of α1 and α2.When the difference between K↓1 and K↓2 or between α1 and α2 is large, as is the case in the summer daytime, this approximation may cause biases.In the winter daytime, the scatter is much reduced (R 2 = 0.96, Figure 9a, RMSE = 0.17 K) than in the summer

Limitations and Future Work
We found that for the 1:1 comparison of the online and offline ∆T s , the daytime results (Figures 8a  and 9a) are more scattered than the nighttime results (Figures 8d and 9d) for both seasons.The derivation of Equation (1) assumes that nonlinear interactions in the change terms can be ignored.For example, concerning the changes in the net shortwave radiation, we assumed: where subscripts 1 and 2 represent EXP_2008 and EXP_2000, respectively, K ↓ is the mean value of K ↓1 and K ↓2 , and α is the mean value of α 1 and α 2 .When the difference between K ↓1 and K ↓2 or between α 1 and α 2 is large, as is the case in the summer daytime, this approximation may cause biases.In the winter daytime, the scatter is much reduced (R 2 = 0.96, Figure 9a, RMSE = 0.17 K) than in the summer daytime (R 2 = 0.78, Figure 8a, RMSE = 0.35 K).In the future study, a more ingenious method should be proposed to furtherly improve the accuracy and reduce the RMSE.Another limitation of this study is that we only consider the real land-use change situation within two separate years (2000 and 2008).In future research, longer time horizon study or virtual afforestation scenarios should be investigated to achieve a deeper comprehension of the surface temperature change mechanism.

Conclusions
In this study, we evaluated the regional impact of afforestation on surface temperature in the Kubuqi desert.We found the effect of afforestation on the surface temperature is 1.34 K, −0.48 K, 2.09 K and 0.22 K for the summer daytime, the summer nighttime, the winter daytime, and the winter nighttime, respectively, for the grid cells that have experienced conversion from bare soil to shrubland (the major land use change in this region).Although the afforestation provides effective dust control and health benefits, the potential warming effect in the daytime should be considered.
The good agreement between the online and the offline calculations of ∆T s found in this study is good evidence that IBPM method performs well when it is connected with the WRF model.According to partitioning results of the IBPM theory, the seasonal variation of surface temperature change is mainly caused by the changes in roughness and Bowen ratio.Also, the daytime surface temperature change is dominated by the surface biophysical effect, in particular, the effect of albedo change.The nighttime temperature change can be attributed equally to surface biophysical effect and atmospheric feedback.
We conclude the atmospheric feedback can amplify the influence of the surface biophysical effect in the afforestation area.Our results showed that ∆T 1 (biophysical effect) and ∆T 2 (atmospheric feedback) have the same sign in the area of bare land to shrubland conversion in all situation.The ratio ∆T 2 /∆T 1 in the daytime is lower than the ratio in the nighttime, and the amplification is even bigger in the nighttime.This effect of atmospheric feedback should be considered in policy making for further afforestation projects, which can help better allocate the afforestation resources.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4907/10/5/368/s1, Figure S1: The experimental design of the two simulations.Figure S2.The difference between online and offline ∆Ts.(a), (c), (e) and (g) represent the results of summer daytime, summer nighttime, winter daytime and winter nighttime, respectively.The x-axis of the hist plots is the difference between offline ∆Ts and online ∆Ts.

Forests 2019 ,of 22 3 a
10, x FOR PEER REVIEW 3 quantity routinely monitored by remote sensing and calculated by land surface models.To generate 95 mechanistic insights into the land-atmosphere interactions, the total Ts changes brought on by 96 afforestation and other land-use activities are often partitioned into component contributions 97 associated with individual biophysical factors.One theory used for the partitioning calculation is the 98

Figure 1 .
Figure 1.Conceptual diagram of biophysical effect and atmospheric feedback.Ta is air temperature, K↓ is incoming shortwave radiation, L↓ is incoming longwave radiation, and the change Δ denotes the difference between different land cover types.ΔTs indicates the surface temperature change.(a) represents the original status before the land-use change; (b) represents the status after the land-use change only considering the biophysical effect; (c) represents the status after the land-use change considering the biophysical effect and atmospheric feedback.

Figure 1 .
Figure 1.Conceptual diagram of biophysical effect and atmospheric feedback.T a is air temperature, K ↓ is incoming shortwave radiation, L ↓ is incoming longwave radiation, and the change ∆ denotes the difference between different land cover types.∆T s indicates the surface temperature change.(a) represents the original status before the land-use change; (b) represents the status after the land-use change only considering the biophysical effect; (c) represents the status after the land-use change considering the biophysical effect and atmospheric feedback.

Figure 2 .Figure 2 .
Figure 2. Illustration of the two nested model domains used in the WRF simulations.

Figure 3 .
Figure 3. (a) Land-use and land-cover map for 2000 in the Kubuqi Desert; (b) land-use and land-cover map for 2008.(c)land-use and land-cover change from 2000 to 2008 (black dots represent the locations of meteorological stations, and the red dot represents the location of an eddy-covariance tower).

Figure 3 .
Figure 3. (a) Land-use and land-cover map for 2000 in the Kubuqi Desert; (b) land-use and land-cover map for 2008.(c) land-use and land-cover change from 2000 to 2008 (black dots represent the locations of meteorological stations, and the red dot represents the location of an eddy-covariance tower).

Forests 2019 ,
10, 368 9 of 21 EXP_2000 (albedo in 2008 minus albedo in 2000) ), Term 1 on the right-hand side of the equation represents the effect of albedo change, Terms 2 and 3 represents the role of energy redistribution associated with roughness change and with Bowen ratio change, respectively, Term 4 represents ground heat flux change, Term 5 is the effect of incoming shortwave radiation change, Term 6 is the effect of incoming longwave radiation change and Term 7 is related to air temperature change.

Figure 4 .
Figure 4. Comparison of monthly mean air temperatures in the EXP_2008 simulation and the observations of 9 meteorological stations.The number on the top of each panel is the station number.

Figure 4 .
Figure 4. Comparison of monthly mean air temperatures in the EXP_2008 simulation and the observations of 9 meteorological stations.The number on the top of each panel is the station number.

Figure 5 .
Figure 5. Comparisons between simulated air temperature, radiation and energy components with observations of eddy covariance tower in a shrub ecosystem from February 14th to February 24th (a) and August 4th to August 14th (b).Ta represents air temperature, H represents sensible heat, LE represents latent heat, G represents ground heat flux.

Figure 5 .
Figure 5. Comparisons between simulated air temperature, radiation and energy components with observations of eddy covariance tower in a shrub ecosystem from 14 February to 24 February (a) and 4 August to 14 August (b).Ta represents air temperature, H represents sensible heat, LE represents latent heat, G represents ground heat flux.

Figure 6 .
Figure 6.Partition of the daytime (a-d) and nighttime (e-h) ∆T s in the summer.(a) daytime ∆T s calculated online by WRF, (b) daytime ∆T s calculated offline with the IBPM theory, (c) biophysical effect (∆T 1 ), (d) atmospheric feedback (∆T 2 ).(e-h) represent the nighttime results corresponding to (a-d).

Figure 7 .
Figure 7. Partition of the daytime (a-d) and nighttime (e-h) ∆T s in the winter.(a) daytime ∆T s calculated online by WRF, (b) daytime ∆T s calculated offline with the IBPM theory, (c) biophysical effect (∆T 1 ), (d) atmospheric feedback (∆T 2 ).(e-h) represent the nighttime results corresponding to (a-d).

Figure 8 .
Figure 8.(a) Relationship between daytime ΔTs in the summer calculated offline with the IBPM theory

Figure 8 .
Figure 8.(a) Relationship between daytime ∆T s in the summer calculated offline with the IBPM theory and the ∆T s calculated online by WRF; (b) the biophysical effect (∆T 1 ) versus the online ∆T s ; (c) the atmospheric feedback (∆T 2 ) versus the online ∆T s .(d-f) represent the nighttime results corresponding to (a-c).The red dash lines represent the 1:1 line and the solid red lines represent the regression results.

Forests 2019 , 14 Figure 8 .
Figure 8.(a) Relationship between daytime ΔTs in the summer calculated offline with the IBPM theory and the ΔTs calculated online by WRF; (b) the biophysical effect (ΔT1) versus the online ΔTs; (c) the atmospheric feedback (ΔT2) versus the online ΔTs.(d)-(f) represent the nighttime results corresponding to (a)-(c).The red dash lines represent the 1:1 line and the solid red lines represent the regression results.

Figure 9 .
Figure 9. (a) Relationship between daytime ΔTs in the winter calculated offline with the IBPM theory and the ΔTs calculated online by WRF; (b) the biophysical effect (ΔT1) versus the online ΔTs; (c) the atmospheric feedback (ΔT2) versus the online ΔTs.(d)-(f) represent the nighttime results corresponding to (a)-(c).The red dash lines represent the 1:1 line and the solid red lines represent the regression results.

Figure 9 .
Figure 9. (a) Relationship between daytime ∆T s in the winter calculated offline with the IBPM theory and the ∆T s calculated online by WRF; (b) the biophysical effect (∆T 1 ) versus the online ∆T s ; (c) the atmospheric feedback (∆T 2 ) versus the online ∆T s .(d-f) represent the nighttime results corresponding to (a-c).The red dash lines represent the 1:1 line and the solid red lines represent the regression results.

Figure 10 .
Figure 10.The frequency histograms of atmosphere feedback effects (K) in the afforestation area.The frequency represents the number of grid cells.(a)-(d) represents the result in the summer daytime, winter daytime, summer nighttime and winter nighttime, respectively.

Figure 10 .
Figure 10.The frequency histograms of atmosphere feedback effects (K) in the afforestation area.The frequency represents the number of grid cells.(a-d) represents the result in the summer daytime, winter daytime, summer nighttime and winter nighttime, respectively.

Figure S3 .
Partition of the daytime intrinsic biophysical effect and atmospheric feedback in summer according to the IBPM theory.(a) effect of albedo change, (b) energy distribution associated with changes in roughness, (c) energy distribution associated with changes in Bowen ratio and (d) effect of soil heat flux changes, (e) effect of incoming shortwave radiation change, (f) effect of incoming longwave radiation change, (g) air temperature related term, (h) ∆Ts calculated offline with the IBPM theory.Figure S4.Partition of the nighttime intrinsic biophysical effect and atmospheric feedback in summer according to the IBPM theory.(a) effect of albedo change, (b) energy distribution associated with changes in roughness, (c) energy distribution associated with changes in Bowen ratio and (d) effect of soil heat flux changes, (e) effect of incoming shortwave radiation change, (f) effect of incoming longwave radiation change, (g) air temperature related term, (h) ∆Ts calculated offline with the IBPM theory.Figure S5.Partition of the daytime intrinsic biophysical effect and atmospheric feedback in winter according to the IBPM theory.(a) effect of albedo change, (b) energy distribution associated with changes in roughness, (c) energy distribution associated with changes in Bowen ratio and (d) effect of soil heat flux changes, (e) effect of incoming shortwave radiation change, (f) effect of incoming longwave radiation change, (g) air temperature related term, (h) ∆Ts calculated offline with the IBPM theory.Figure S6.Partition of the nighttime intrinsic biophysical effect and atmospheric feedback in winter according to the IBPM theory.(a) effect of albedo change, (b) energy distribution associated with changes in roughness, (c) energy distribution associated with changes in Bowen ratio and (d) effect of soil heat flux changes, (e) effect of incoming shortwave radiation change, (f) effect of incoming longwave radiation change, (g) air temperature related term, (h) ∆Ts calculated offline with the IBPM theory.Author Contributions: L.W. and D.F. generated the data.Y.Y. (Yizhou Yin) and Y.L. provided the basic data.All authors analyzed and discussed the data.The manuscript was written by L.W. and X.L. C.F., Z.W., Y.Y. (Yanzheng Yang) and G.L. revised the manuscript.Funding: This work is supported jointly by the Department of Hydraulic Engineering, the Department of Earth System Science, Tsinghua University and the School of Forestry and Environmental Studies, Yale University.We are grateful for the financial support from Ministry of Science and Technology of P.R.China (2016YFC0402701), National Science Foundation of China (NSFC 51825902), National Basic Research Program of China (2013CB956601) and China Scholarship Council.

Table 1 .
Model configuration and main physical parameterizations used for the two simulations.