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

Abstract

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.

1. 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₂ 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 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).

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 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.

2. Methods

2.1. 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. Model configuration and main physical parameterizations used for the two simulations.

Items	Description
Model version	WRF 3.7.1
Dynamics solver	Advanced Research WRF
Time step	30 s
Output interval	1 h
Vertical level	27
Radiation scheme	CAM3 a
Surface model	Noah land surface model [37,38]
Cumulus scheme	Grell-Freitas ensemble scheme [39]
Microphysic scheme	WSM3 b
PBL scheme	YSU c
Surface layer	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.

The major land-use changes in the Kubuqi Desert during 2000–2008 were conversions from bare land to shrubland (Figure 3c,

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⁻²), L_↓ is incoming longwave radiation (W m⁻²), T_a is air temperature (K) and the change Δ denotes the difference between EXP_2008 and EXP_2000. ΔT₁ denotes the surface temperature change caused by intrinsic biophysical effect and ΔT₂ denotes the surface temperature change caused by atmospheric feedback.

Land Change Type	Area (km2)	Time	ΔK↓	ΔL↓	ΔTa	ΔT1	ΔT2	Albedo Term1	Roughness Term2	Bowen Ratio Term3	Soil Heat Flux Term4	K↓ Term5	L↓ Term6	Ta Term7	Sum (Terms 1 through 7)
No change	30,552	winter daytime	−0.11	0.34	0.17	0.02	0.18	0.02	0.02	−0.02	0.00	0.00	0.01	0.17	0.20
		winter nighttime	0.00	0.11	0.00	−0.01	0	0.00	−0.02	0	0.01	0.00	0.00	0.00	−0.01
		summer daytime	−0.84	0.50	0.07	−0.03	0.07	0.01	0.04	−0.05	−0.03	−0.01	0.01	0.07	0.04
		summer nighttime	0.00	−0.28	−0.13	−0.06	−0.06	0.00	−0.06	−0.02	0.02	0.00	0.00	−0.06	−0.12
Bare land to shrub	10,268	winter daytime	−3.49	2.56	0.63	1.51	0.58	1.11	−0.05	0.37	0.08	−0.05	0.06	0.57	2.09
		winter nighttime	0.00	0.92	0.04	0.12	0.1	0.00	0.07	0.01	0.04	0.00	0.01	0.09	0.22
		summer daytime	−9.31	6.83	0.33	0.97	0.37	1.74	−0.44	−0.16	−0.17	−0.14	0.18	0.33	1.34
		summer nighttime	0.00	2.02	−0.39	−0.26	−0.22	0.00	−0.37	0.01	0.10	0.00	0.02	−0.24	−0.48
Shrub to grass	2332	winter daytime	−1.97	1.62	0.35	0.68	0.31	0.77	−0.28	0.19	0.00	−0.04	0.05	0.30	0.99
		winter nighttime	0.00	0.86	0.05	−0.02	−0.01	0.00	−0.06	0.01	0.03	0.00	0.04	−0.05	−0.03
		summer daytime	−4.69	1.95	0.00	−0.44	0	0.87	−0.13	−0.97	−0.21	−0.09	0.05	0.04	−0.44
		summer nighttime	0.00	1.34	−0.48	−0.16	−0.19	0.00	−0.25	−0.14	0.23	0.00	0.04	−0.23	−0.35
Cropland to grass	992	winter daytime	0.29	−0.58	0.13	0.03	0.12	−0.14	0.20	−0.02	−0.01	0.01	−0.02	0.13	0.15
		winter nighttime	0.00	−0.64	−0.01	−0.13	−0.07	0.00	−0.17	0.01	0.03	0.00	−0.03	−0.04	−0.20
		summer daytime	0.33	−2.67	0.01	0.13	−0.02	−0.24	0.24	0.11	0.02	0.01	−0.06	0.03	0.11
		summer nighttime	0.00	−2.81	−0.16	−0.06	−0.19	0.00	−0.06	−0.01	0.01	0.00	−0.09	−0.10	−0.25
Shrub to cropland	644	winter daytime	−2.53	1.78	0.32	0.59	0.27	0.88	−0.58	0.23	0.06	−0.04	0.05	0.26	0.86
		winter nighttime	0.00	1.39	0.24	0.22	0.25	0.00	0.25	0.01	−0.04	0.00	0.04	0.21	0.47
		summer daytime	−5.83	4.34	−0.01	−0.69	−0.01	0.95	−0.52	−1.02	−0.10	−0.09	0.09	−0.01	−0.7
		summer nighttime	0.00	4.39	−0.36	−0.24	−0.10	0.00	−0.39	−0.06	0.21	0.00	0.09	−0.19	−0.34
Whole domain	46,648	winter daytime	−0.95	0.86	0.29	0.38	0.29	0.30	−0.02	0.08	0.02	−0.01	0.02	0.28	0.67
		winter nighttime	0.00	0.30	0.03	0.03	0.03	0.00	0.01	0.00	0.02	0.00	0.01	0.02	0.06
		summer daytime	−2.88	1.91	0.12	0.01	0.14	0.28	−0.08	−0.12	−0.07	−0.04	0.05	0.13	0.15
		summer nighttime	0.00	0.28	−0.11	−0.1	−0.1	0.00	−0.12	−0.02	0.04	0.00	0.00	−0.10	−0.2

). The area of bare land decreased from 15,711 km² in 2000 to 6176 km² 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.).

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.

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]).

2.2. Offline Calculations

According to the IBPM theory [28], the surface temperature is expressed as:

$$T_s=T_a+\frac{{\lambda }_0}{1+f}\left(R_n^{*}-G\right)$$

(1)

where T_s is surface temperature (K), T_a is air temperature (K) at the blending height, G is ground heat flux (W m⁻²), λ₀ (K W⁻¹ m⁻²) is given by:

$${\lambda }_0=\frac{1}{4\sigma T_s^3}$$

(2)

where σ is the Stephan-Boltzmann constant, f is an energy redistribution factor (dimensionless), and $R_n^{*}$ is apparent net radiation:

$$R_n^{*}=K_{\downarrow }\left(1-\alpha \right)+L_{\downarrow }-\sigma T_a^4$$

(3)

where K_↓ is incoming shortwave radiation (W m⁻²), α is surface albedo, and L_↓ is incoming longwave radiation (W m⁻²), the proof of Equation (1) is given in the supporting information.

The energy redistribution factor f is given by:

$$f=\frac{\rho C_p}{4r_t\sigma T_s^3}\left(1+\frac{1}{\beta }\right)$$

(4)

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

$$\begin{array}{cccc}\Delta T_s\approx & \frac{{\lambda }_0}{1+f}\left(-K_{\downarrow }\Delta \alpha \right)+& \frac{-{\lambda }_0}{{\left(1+f\right)}^2}\left(R_n^{*}-G\right)\Delta f_1+& \frac{-{\lambda }_0}{{\left(1+f\right)}^2}\left(R_n^{*}-G\right)\Delta f_2+\\ & \mathrm{Term} 1& \mathrm{Term} 2& \mathrm{Term} 3\\ \frac{-{\lambda }_0}{1+f}\Delta G& +\frac{{\lambda }_0}{1+f}\left(1-\alpha \right)\Delta K_{\downarrow }& +\frac{{\lambda }_0}{1+f}\Delta L_{\downarrow }& +\left(\Delta T_a-\frac{{\lambda }_0}{1+f}\sigma \Delta T_a^4\right)\\ \mathrm{Term} 4& \mathrm{Term} 5& \mathrm{Term} 6& \mathrm{Term} 7\end{array}$$

(5)

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 EXP_2000 (albedo in 2008 minus albedo in 2000). The change in the energy redistribution factor consists of changes in surface roughness (Δf₁) and Bowen ratio (Δf₂):

$$\Delta f_1=\frac{-\rho C_p}{4r_t\sigma T_s^3}\left(1+\frac{1}{\beta }\right)\frac{\Delta r_t}{r_t}$$

(6)

$$\Delta f_2=\frac{-\rho C_p}{4r_t\sigma T_s^3}\left(\frac{\Delta \beta }{{\beta }^2}\right)$$

(7)

In Equation (5), 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.

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

$$\beta =\frac{H}{LE}$$

(8)

$$r_t=\rho C_p\frac{T_s-T_a}{H}$$

(9)

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², 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].

3. Results

3.1. 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² 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 2.4 W m⁻², −12.8 W m⁻², 2.5 W m⁻², 63.9 W m⁻² and −14.5 W m⁻² 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 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).

3.2. 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 (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 Δ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² > 0.9) except for the summer daytime (R² = 0.78).

3.3. 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.

3.4. 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.48 K), 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.24 K) 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.

4. Discussion

4.1. Warming Effect in the Daytime

As shown in Figure 6 and Figure 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⁻¹. 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.

4.2. 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.05 K), 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⁻² (daytime mean) in the summer and −3.49 W m⁻² in the winter. Resulting from the cloud changes, the incoming longwave radiation (L_↓) increases by 6.83 W m⁻² (daytime mean) in the summer and 2.56 W m⁻² 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⁻² (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 (Figure 6e–h and Figure 7e–h). In the summer, the total contribution of atmospheric feedback is −0.22 K, whereas the total biophysical effect is −0.26 K. 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⁻² (nighttime mean value) in the summer and 0.92 W m⁻² 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⁻² (nighttime mean of summer and winter seasons).

Figure 8a–c and Figure 9a–c also show that the offline ΔT_s in the daytime is dominant by biophysical effect (ΔT₁). Most of the feedback effect (ΔT₂) scatters in the range from −0.5 K to 0.5 K, or about 1/4 of the range of ΔT₁. When the regression slope is closer to 1 (Figure 8b,c,e,f, and Figure 9b,c,e,f), the corresponding factor contributes more to the ΔT_s. In the summer, the regression slope between ΔT₁ and online daytime ΔT_s is 0.86 (R² = 0.70), but the slope value between ΔT₂ and online daytime ΔT_s is 0.24 (R² = 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₁; the regression slope is 0.98, which is very close to 1, with a higher R² (0.90), indicating the biophysical effect is still the dominant factor. However, the regression slope between ΔT₂ and ΔT_s is 0.28 (R² = 0.65) in the winter daytime.

In the nighttime, both ΔT₁ and ΔT₂ show a significant linear relationship with online ΔT_s. The regression slopes in the summer are 0.57 (ΔT₁) and 0.51 (ΔT₂) with high R² values (0.95–0.98). The values in the winter are 0.54 (ΔT₁, R² = 0.99) and 0.47 (ΔT₂, R² = 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₂, 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.

4.3. Limitations and Future Work

We found that for the 1:1 comparison of the online and offline ΔT_s, the daytime results (Figure 8a and Figure 9a) are more scattered than the nighttime results (Figure 8d and Figure 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:

$$K_{\downarrow 1}\left(1-{\alpha }_1\right)-K_{\downarrow 2}\left(1-{\alpha }_2\right)\approx -K_{\downarrow }\Delta \alpha +\left(1-\alpha \right)\Delta K_{\downarrow }$$

(10)

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 α₁ and α₂. When the difference between K_↓1 and K_↓2 or between α₁ and α₂ 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² = 0.96, Figure 9a, RMSE = 0.17 K) than in the summer daytime (R² = 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.

5. 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₁ (biophysical effect) and ΔT₂ (atmospheric feedback) have the same sign in the area of bare land to shrubland conversion in all situation. The ratio ΔT₂/ΔT₁ 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.