Assessing Global and Regional Effects of Reconstructed Land-Use and Land-Cover Change on Climate since 1950 Using a Coupled Land–Atmosphere–Ocean Model

: Land-useand land-coverchange(LULCC)is one of themost important forcingsaffectingclimatein the past century. This study evaluates the global and regional LULCC impacts in 1950–2015 by employing an annually updated LULCC map in a coupled land–atmosphere–ocean model. The difference between LULCC and control experiments shows an overall land surface temperature (LST) increase by 0.48K in the LULCC regions and a widespread LST decrease by 0.18K outside the LULCC regions. A decomposed temperature metric (DTM) is applied to quantify the relative contri- bution of surface processes to temperature changes. Furthermore, while precipitation in the LULCC areas is reduced in agreement with declined evaporation, LULCC causes a southward displacement of the intertropical convergence zone (ITCZ) with a narrowing by 0.5 8 , leading to a tripole anomalous precipitation pattern over the warm pool. The DTM shows that the temperature response in LULCC regions results from the competing effect between increased albedo (cooling) and reduced evaporation (warming). The reduced evaporation indicates less atmospheric latent heat release in convective processes and thus a drier and cooler troposphere, resulting in a reduction in surface cooling outside the LULCC regions. The southward shift of the ITCZ implies a northward cross-equatorial energy transport anomaly in response to reduced latent/sensibleheat of the atmosphere in the Northern Hemisphere, where LULCC is more intensive. Tropospheric cooling results in the equatorward shift of the upper-tropospheric westerly jet in both hemispheres, which, in turn, leads to an equatorward narrowing of the Hadley circulation and ITCZ.


Introduction
Land-use and land-cover change (LULCC) consists of a wide range of land surface conversions including the conversion from forests to crops and pasturelands, reforestation of formerly agricultural areas, afforestation, and all kinds of urbanization (Mahmood et al. 2014). From the years 1700 to 2000, 42%-68% of the global land surface has been transformed from natural vegetation into agriculture, rangeland, and woodland (Hurtt et al. 2006). By the end of the twentieth century, most LULCC had taken place over the temperate regions of the Northern Hemisphere (NH), such as East Asia, South Asia, Europe, and North America (Ramankutty and Foley 1999). After a rapid increase in the rate of deforestation during the 1980s and 1990s, recent satellite observations indicate a slowdown of LULCC in the past decade (FAO 2007). Reforestation and afforestation are observed in western Europe, North America, and China as a result of land abandonment and afforestation efforts (Nagendra and Southworth 2009), while deforestations are still concentrated in the tropics.
Previous observation studies agree that LULCC modifies the surface properties through changing soil and vegetation characteristics, including albedo, leaf area index (LAI), and roughness length (Gash and Nobre 1997;Silva Dias et al. 2002;Lyons 2002;Chagnon et al. 2004; Davin and de Noblet-Ducoudre 2010;Nair et al. 2011;Guo et al. 2016;Liu et al. 2016). When LULCC takes place, conversions from forests to crops/grasslands normally cause a higher albedo because the multiple reflections within the canopy are diminished in low vegetation. A reconstructed land-cover and biophysical parameter map by Steyaert and Knox (2008) showed that the albedo of the eastern United States increased by ;5% from 1650 to 1992, caused by the extensive LULCC that replaced deciduous forests and native grasslands with agricultural crops and pastures. Albedo contrast due to LULCC is especially amplified when snow is present: open land is normally entirely snow-covered, meaning the albedo can be as high as 0.9, whereas tree crowns will remain exposed above the snow in the forest (Betts 2001). The increase in albedo reflects more incident radiation and leads to a change in radiative forcing by 20.15 6 0.10 W m 22 globally, according to the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5) (Myhre et al. 2013). The reduction in net radiation tends to reduce the local land surface temperature (LST) via the radiation budget, referred to as the albedo effect (Davin et al. 2007).
Meanwhile, changes in LAI, vegetation cover, and surface roughness decrease the evapotranspiration on the leaf surface (Forster et al. 2007;Strack et al. 2008;Peng et al. 2014;Bright et al. 2017;Liu et al. 2019). A conversion from forest with deeper root depth to grassland/crop with shallower root depth results in a depressed latent heat (LH) within the planetary boundary layer (PBL) and potentially modifies atmospheric humidity and canopy air temperature. The change in aerodynamic roughness alters the surface wind speed and thus changes the momentum and water vapor exchange between land and atmosphere (Lee et al. 2011). Therefore, a higher Bowen ratio (the ratio of sensible heat to latent heat fluxes) and a modification in water cycles is found due to LULCC. The reduced latent and sensible heat (SH) fluxes tend to increase the LST by modifying the surface energy balance, referred to as the evaporation effect.
The LST response to LULCC is highly heterogeneous and can be warming or cooling, depending on its location and magnitude associated with the relative magnitude of albedo effect and evaporation effect (Betts 2000;Davin and de Noblet-Ducoudre 2010;Arora and Montenegro 2011;Li et al. 2015). Various methods have been proposed to evaluate the relative contributions of radiation and biophysical properties such as LH, SH, and aerodynamic resistance to the temperature response (Lee et al. 2011;Luyssaert et al. 2014;Chen and Dirmeyer 2016;Rigden and Li 2017). One of the widely used methods is the intrinsic biophysical mechanism (IBM), which attributes changes in LST to changes in net radiation, aerodynamic resistance, and the Bowen ratio (Lee et al. 2011). This approach has been widely used to investigate changes of temperature and the contributing factors (Lee et al. 2011;Peng et al. 2014;Zhao et al. 2014;Bright et al. 2017) and has been further developed to include atmospheric feedbacks (Chen and Dirmeyer 2016). However, it is reported that the IBM overestimates the contribution of aerodynamic resistance by assuming independence between the aerodynamic resistance and the Bowen ratio (Rigden and Li 2017). As the Bowen ratio is directly influenced by albedo and aerodynamic resistance, the evaluation of its contributions to temperature changes also includes those from radiation and resistance.
Another decomposition metric decomposes changes in LST to changes in incoming radiation, surface albedo, ground heat, sensible heat, and latent heat fluxes, providing a detailed breakdown of all components in the surface energy budget (Juang et al. 2007;Luyssaert et al. 2014;Xu et al. 2015;Chen and Dirmeyer 2016;Winckler et al. 2017;Hirsch et al. 2018). However, it does not link changes in turbulent fluxes with changes in biophysical properties such as surface roughness (Rigden and Li 2017). Moreover, it does not include the effect of atmospheric temperature and moisture conditions on LST. This method is further developed in this study to link latent and sensible heat fluxes to surface conditions (temperature and vapor pressure), atmospheric feedbacks (temperature and vapor pressure), and surface roughness, thus providing a diagnostic framework for evaluating the relative contribution of each variable to LST changes in general circulation model (GCM) simulations.
The climate effect of LULCC on local and regional precipitation has been thoroughly observed and simulated in the past few years. Many studies agree that LULCC causes a decrease in local precipitation due to local land-atmosphere interactions, such as a decrease in evapotranspiration, soil moisture, and clouds (Xue 1996;Pitman et al. 2004;Xue et al. 2016;Quesada et al. 2017). Some studies point out that the effects of afforestation and deforestation on monsoon precipitation depend on the relative location of the perturbed area and largescale circulation (Xue and Shukla 1996;Boone et al. 2016). Meanwhile, there remains a lack of a comprehensive understanding of how LULCC influences remote and global climate because the LULCC effects on precipitation and temperature outside the LULCC areas are inconsistent among previous model simulations (Pielke et al. 2007;Pitman et al. 2009;Pielke et al. 2011;Swann et al. 2012;Devaraju et al. 2015;Lejeune et al. 2017). Some early GCM simulations report that drastic LULCC in the tropical regions can affect remote regions through atmospheric teleconnections, although their mechanisms need further investigation (Chase 1995;Zhao et al. 2001;Snyder 2010;Schneck and Mosbrugger 2011). The Land-Use and Climate, IDentification of robust impacts (LUCID; Pitman et al. 2009) models are atmosphere-only GCMs with fixed sea surface temperatures (SST), which show no remote impact of LULCC. Using model results from phase 5 of the Climate Model Intercomparison Project (CMIP5), recent studies show that LULCC induces warming in low latitudes and cooling in high latitudes (Winckler et al. 2019a,b). However, there is less agreement about the latitude at which deforestation starts to have a cooling effect.
A few early atmospheric general circulation model (AGCM) studies report that LULCC could change the largescale circulation by shifting the ITCZ in the North African continent (Xue and Shukla 1993), displacing the African easterly jet (Li et al. 2007) and changing the strength and position of the Hadley and Walker circulation (Zhao et al. 2001;Snyder 2010). Recent literature shows that the global impact of LULCC by previous modeling studies can be reviewed and explained from the perspective of the global energy budget and energy transport (Swann et al. 2012;Devaraju et al. 2015;Laguë and Swann 2016;Devaraju et al. 2018). The results consistently show that continental deforestation/afforestation can exert a significant impact on global-scale circulations, tropical precipitation, and cloud cover through atmospheric teleconnections. This energy transport approach has been widely used in different studies to explore the global impact of large-scale forcings, such as aerosol radiative effects (Ming and Ramaswamy 2009;Ming et al. 2011;Haywood et al. 2016) and Arctic sea ice loss (Chiang and Bitz 2005;Broccoli et al. 2006;Deser et al. 2015;Tomas et al. 2016;Cvijanovic et al. 2017). It seems reasonable that drastic LULCC may also exert a global impact as one of the most important forcings in the past century (Mahmood et al. 2014). However, current studies of the global impact of LULCC are based on idealized experiment designs: either they convert all forests to grasses (Devaraju et al. 2015) or they replace all C3 grasses/agricultures with deciduous trees (Swann et al. 2012;Laguë and Swann 2016). Moreover, these studies applied slab ocean models, which may overestimate the local thermodynamic impact of air-sea coupling and may even modify the spatial structure of the response compared to the full-depth ocean (Deser et al. 2016). Therefore, a fully coupled ocean model that includes both dynamical and thermodynamic processes is needed for a proper representation of oceanic feedbacks. In this study, the remote impact of LULCC will be identified, and the key processes controlling the remote impact will be analyzed by using a coupled land-atmosphere-ocean GCM. Our study mainly focuses on how LULCC would invoke a global impact on precipitation and large-scale circulation, especially on the ITCZ and Hadley circulation.
In this study, we aim to investigate the global biophysical impacts of LULCC and identify the processes that affect landatmosphere-ocean interactions. A LULCC map with interannual variations has been developed based on historical and future land-use data from Hurtt et al. (2006Hurtt et al. ( , 2011. By implementing the yearly-updated vegetation map into a GCM, we provide an assessment of large-scale LULCC on the global climate since 1950. In section 2, we provide a brief introduction to the annually updated LULCC map used in the simulation, a description of the CFSv2 GCM, experiment design, temperature decomposition metric, and the statistical method. The responses of global and regional temperature and precipitation, along with the key mechanisms that control those responses, are presented in section 3. Discussion and conclusions are given in section 4.
2. Model, methodology, and experiment design a. Model, LULCC map, and experiment design In this research, we adopt the second version of the National Centers for Environmental Prediction (NCEP) Climate Forecast System (CFSv2), a fully coupled atmosphere-oceanland model (Saha et al. 2014). The atmospheric model is NCEP's Global Forecast System (GFS) with the horizontal resolution of the model set at T126, approximately 100-km grid resolution. Sixty-four vertical levels are used, most of which are in the troposphere. The ocean model is from the Geophysical Fluid Dynamics Laboratory Modular Ocean Model version 4 (MOM4; Griffies et al. 2004). The NCEP CFS is coupled with the second generation of the Simplified Simple Biosphere Model (SSiB2) as the land surface model (Xue et al. 1991;Zhan et al. 2003; hereafter CFSv2/SSiB2). The CFSv2/SSiB2 has been extensively used to study different aspects of the climate, including the influence of atmosphere-ocean interactions (Lee et al. 2019), aerosol indirect effects (Huang et al. 2019), and the potential role of spring LST anomaly on seasonal prediction (Xue et al. 2018;Diallo et al. 2019). A detailed description of the parameterizations of this version of CFSv2/ SSiB2 can be found in Huang et al. (2019).
The SSiB2 is a state-of-the-art vegetation biophysical model that includes climate-vegetation biophysical processes by explicitly considering the exchanges of energy and water between the land and the atmosphere (Xue et al. 1991;Sellers et al. 1996;Zhan et al. 2003). The vegetation map applied in CFSv2/ SSiB2 includes 13 types (see Fig. S1 in the online supplemental material) as in Xue et al. (2004), and each grid cell has one vegetation type. Each vegetation type includes a set of parameters, including roughness length, LAI, displacement height, greenness, and vegetation fractional coverage with climatological seasonal variations. SSiB2 has been implemented for global and regional climate simulations and highlights the importance of land-atmosphere interaction studies (Kang et al. 2007;Li et al. 2007;Ma et al. 2013a,b;Boone et al. 2016).
The LULCC map is generated based on the Land-Use Harmonization 2 (LUH2) datasets from Hurtt et al. (2006Hurtt et al. ( , 2011. The global gridded land-use data in LUH2 include 600 years of annual land-use transition mapping for the period of 1500-2100, which have been widely used as forcing data for global LULCC simulations including the CMIP5/CMIP6 (Findell et al. 2009;Brovkin et al. 2013;Lawrence et al. 2016). In this study, the crop and pasture fractions in LUH2 have been added to obtain an estimate of the total land-cover change. The combined crop and pasture fraction changes are defined as the LULCC fraction. Figure 1 shows the snapshots of the LULCC fraction relative to 1948 in the years 1950, 1970, 1990, and 2015. Significant land conversions occur after 1970 and continue for several decades, except for over the southern United States and western Europe where there is a clear decrease in the amount of land converted to crops or pastures .
A methodology has been developed to convert the LULCC fraction map to the annually updated vegetation map. Once the LULCC fraction exceeds a threshold value (a) or the fraction change compared to 1948 exceeds a threshold value (b), the area is degraded annually based on the default SSiB2 land-cover classification map in Fig. S1. The forests are degraded to low vegetation (grassland and shrub), and low vegetation is degraded to bare soil regions. By degrading plant functional types (PFTs), the associated surface parameters (LAI, transmittance, greenness, fractional coverage, vegetation height, displacement height, roughness length, soil parameters) are changed to show the LULCC effects. In this study, we only consider the simplest land degradation from forests to low vegetation and from low vegetation to bare soil and neglect other types of LULCC such as reforestation, irrigation, and urbanization.
The methodology for developing the annually updated vegetation map is used in Chilukoti and Xue (2020) to assess global impact of LULCC. After comparing with observed temperature and precipitation, Chilukoti and Xue (2020) concluded that LULCC simulations with the annually updated vegetation map had a smaller bias and a better interannual variability compared to control simulations with the potential vegetation map. The method in Chilukoti and Xue (2020) is further developed in this study by applying different criteria a and b in different time periods and 10 degraded regions in Table 1. For each region, criteria a and b are chosen to make the fraction of land degradation grids in the SSiB2 vegetation map comparable to the LULCC fraction in LUH2 (Fig. S2). Intensive LULCC takes place before 1990 because of agricultural activity. There are a few regions (East China, Tibet, India, Mexico, and Australia) in which no further degradation is found after 1980 and some regions even show recovery. However, degradation continues in tropical regions, which may result in significant changes in the tropical climate. Compared with the potential vegetation map, about 20% of global land areas have land degradation in 2015 in the LULCC simulations.
GCM ensemble simulations are conducted to investigate the effect of LULCC on the global climate ( Table 2). The sensitivity cases (LULCC) apply the vegetation map that is degraded annually based on the LUH2 and are compared with the control cases (CTL) using the original SSiB2 vegetation map in Xue et al. (2004). Both CTL and LULCC include three ensemble members whose initial conditions are obtained from years 8, 9, and 10 from the spinup simulations. Previous LULCC studies have demonstrated that conducting ensemble simulations can effectively screen the model variability, which would otherwise incorrectly hint at a global teleconnection (Lorenz et al. 2016). The ensemble simulations start from 1949 and run with atmospheric CO 2 concentrations in 1949-2015 from the World Meteorological Organization (WMO) Global Atmospheric Watch.
Since the CFS reanalysis (CFSR; Saha et al. 2010) starts from the year 1979, spinup runs are required to ensure a thermodynamically balanced initial condition for the simulations

b. Decomposition metric
We applied the decomposed temperature metric (DTM) developed by Juang et al. (2007) and Luyssaert et al. (2014) to quantify the relative dynamic contribution of albedo, incoming shortwave radiation, incoming longwave radiation, latent heat, sensible heat, and ground heat flux to temperature changes. The energy balance in the DTM can be arranged as follows: where a s is surface albedo, SW in is incoming shortwave radiation, LW in is incoming longwave radiation, LST is the surface temperature, s is the Stephan-Boltzmann constant, LH is latent heat flux, SH is sensible heat flux, and G is ground heat flux. The changes in the energy fluxes and LST can be acquired by applying the first-order derivative of the left-hand side and right-hand side of Eq. (1) as in Luyssaert et al. (2014): (2) The delta (D) indicates the LULCC minus CTL in our experiment. The residual term DI represents the residual due to the high-order derivative, which includes the covariance between attributing variables, which could not be quantified in the decomposition metric. As such, the contribution of energy budget changes to LST changes can be evaluated using Eq.
(2). The DTM is developed to further relate LH and SH with surface biophysical property and atmospheric condition changes induced by LULCC. Based on the bulk transfer equation (Verhoef et al. 1997), LH and SH over the vegetated surface are related to surface conditions, atmospheric variables, and aerodynamic resistance (Xue et al. 1991) as follows: LH 5 rC p l e s * 2 e a r a 1 r c , where r is the air density, C p is the specific heat of air at constant pressure, T a is the air temperature at the reference height, r a is aerodynamic resistance, l is the psychrometric constant (66.5 Pa 8C 21 ), e s * is saturated vapor pressure at the surface, e a is vapor pressure at the reference height, and r c is bulk stomatal (canopy) resistance. We assume r a 1 r c is the surface resistance for water vapor to distinguish it with aerodynamic resistance for heat (r a ).
The change of SH is decomposed to changes in LST, T a , and r a by implicit first-order derivative as in Eq. (2): The terms on the right-hand side of the equation indicate that the changes of SH are positively related to changes in LST and is negatively related to changes in atmospheric temperature and aerodynamic resistance. The LH response is decomposed in a similar way: It is assumed that the change of LH is positively related to change in surface moisture and is negatively related to changes in atmospheric moisture and resistance. Equations (4a) and (4b) will be employed in the analysis to evaluate the contribution of surface temperature/vapor pressure, atmospheric temperature/vapor pressure, and surface resistance to turbulent flux changes.
c. The energy flux framework The energy flux framework developed by Kang et al. (2009) has been widely used to understand changes of Hadley circulation and ITCZ caused by atmospheric thermal forcing, such as aerosol radiative effects, Arctic sea ice-cover changes, asymmetric In the long-term average when the heat storage can be neglected, the energy transport across latitude u, F a (u), can be calculated from the energy budget within the atmosphere (Hartmann 1994): where a is Earth's radius. The atmospheric energy (R ATM ) is the thermal heating within the atmosphere calculated using top-of-atmosphere radiation (R TOA ), surface radiation (R SFC ), LH, and SH [see Eq. (S2) in the online supplemental material]. When atmospheric thermal forcing occurs in one hemisphere (i.e., a cooling anomaly in the Northern Hemisphere), the energy imbalance between hemispheres necessitates a cross-equatorial energy transport toward the cooled hemisphere (the NH). In the tropics, the anomalous energy transport across the equator is realized through anomalous mass transport toward the Northern Hemisphere by the upper branch of the Hadley circulation , which would be accomplished by the southward shift of the rising branch of Hadley circulation and ITCZ (Kang et al. 2009). A detailed description of how the Hadley circulation and ITCZ position are related to the equatorial energy transport can be found in supporting information. Please note that we do not analyze the possible mechanisms that cause seasonal changes in ITCZ and Hadley circulation in this paper. This is because the energy flux framework we use to diagnose the shift of ITCZ is valid on the annual average when the rate of energy change in the atmosphere can be considered as negligible [i.e., ›E a /›t 5 0. in Eq. (S1)]. For a specific season, the rate of energy change cannot be neglected, which is beyond the scope of this study.

d. Statistical methods
For the three-member ensembles, we calculate the ensemble mean first before the significance tests. Then the statistical significance of differences between LULCC and CTL is evaluated using the Student's t test (Zwiers and von Storch 1995). To account for the correlation in space, we perform a field significance test by estimating the false discovery rate (FDR; Wilks 2006). The FDR not only tests field significance but also identifies the local grid points with p values small enough to reject the global null hypothesis (Lorenz et al. 2016).

Results
a. Global and regional changes of temperature, precipitation, and Hadley circulation due to large-scale LULCC In this study, we focus on annual mean results averaged in 1950-2015 over all ensemble members in LULCC and CTL. Figure 2 shows the difference of annual mean surface temperature between LULCC and CTL averaged over 1950-2015. We found a warming by 0.48 K averaged over the LULCC grids and a weak cooling (20.18 K) averaged over the widespread non-LULCC land areas (Fig. 2b). Over the globe, a slight cooling is found by 20.08 K. The warming regions are concentrated in low-latitude degraded regions, such as West Africa, southern Africa, South America (Brazil and southern Argentina), and the Tibetan Plateau; in the midlatitude regions, such as central Asia, colder LST is found over the regions of land degradation (Fig. 2a). Meanwhile, LULCC causes a widespread cooling by an average of 0.1 K over the ocean. The three ensemble members have similar temperature responses to LULCC in all regions in Fig. 2b. As shown in previous studies, the response of LST is heterogeneous over different land regions and is influenced by the competing albedo (cooling) effect and evaporation (warming) effect (Chen and Dirmeyer 2016). In section 3b, we will apply the decomposition method described in section 2c to quantify the contribution of each surface variable to temperature changes and discuss the opposing LST responses in lower and midlatitudes and the cooling in the nondegraded regions. For seasonal response, LULCC has caused temperature increases by 1.5 K in local warm seasons and by less than 0.5 K in local cold seasons (Fig. S3) at the low-latitude degraded regions (West Africa, southern Africa, South America, and tropical Australia). In midlatitude central Asia, the surface temperature is cooled by about 1.5 K in winter and spring and warmed by 1 K in summer. Similar to the annual results, a widespread cooling is found in non-LULCC land regions and widespread ocean areas. Compared to the temperature change in LULCC regions, the temperature responses over non-LULCC regions are relatively small but consistent throughout the year.
LULCC also has substantial effects on precipitation. As shown in Fig. 3a, annual precipitation reduction is found over degraded areas, including northwest and southern Africa, South America, East and South Asia, and southern North America, as well as over the adjacent ocean regions. The precipitation response over the degraded region is strongest during the local monsoon seasons: the change peaks in DJF in Southern Hemisphere while it peaks in JJA in Northern Hemisphere (Fig. S4). Around the western tropical Pacific, a tripole anomaly pattern can be identified: precipitation increases along the equator along with significant decreases occurring to the south and the north of the equatorial region. The ITCZ boundary narrows by about 0.58 following the definition in Byrne and Schneider (2016) (see section 3 of the supplemental material). Moreover, the rising branch of the Hadley circulation is displaced southward compared to its climatology location (Fig. 3b), relocating the centroid of the convergence zone southward. In DJF and MAM (Figs. S4a,b), the contraction of the ITCZ is also shown in the tropical Pacific: the precipitations increase along the equator but decrease in the southern and northern boundaries of ITCZ. The signal is weak in JJA and SON. The changes of ITCZ and largescale circulation on annual scale will be analyzed in section 3c from the perspective of energy budget and meridional energy transport within the climate system. b. LST response and its contribution from radiation and turbulent fluxes Figure 4 shows the differences of LST between LULCC and CTL and decomposed LST change as a function of radiation, turbulent fluxes, and heat storage as introduced in Eq. (2). LULCC results in an increase in surface albedo because of the higher reflectivity of bare soil and low vegetation compared with tall vegetation. The changes in albedo contribute to a decrease of LST restricted to LULCC regions (Fig. 4b). The widespread increase in SW in tends to warm the surface (Fig. 4c) and the extensive decrease in LW in has a cooling effect on the surface (Fig. 4d). The changes in SW in and LW in are related to climate feedback and will be explained later. Overall, both LH and SH decrease after land-cover changes, indicating a general warming effect on the land surface (Figs. 4e,f). Ground heat flux change plays a minor role in LULCC-induced LST change (Fig. 4g). The minor differences between LULCC-induced LST change (Fig. 4a) and the summed LST changes from each contributor (Fig. 4h) show that the DTM can reasonably delineate the contribution of each process. The high-order terms are larger over the semiarid regions such as central Asia, inland China, and the Sahel (Fig. 4i).
The changes in SW in and LW in are related to atmospheric moisture/temperature conditions (Fig. 5) and cloud feedback (Fig. S5). LULCC causes less evaporation to the PBL (Xue et al. 1997). The deficit in low-level specific humidity can extend to about 500 hPa and spread outside the LULCC regions (Fig. 5a). The decrease in water vapor (Fig. S5a) is expected to reduce the latent heat release during convective processes and cool the tropospheric temperature, consistent with previous studies showing a reduction of convective heating due to LULCC (Ma et al. 2013a,b;Song 2013). The cooling signal is amplified with increasing elevation, with the maximum cooling anomalies occurring around 300 hPa (Fig. 5b), which may be due to the larger moist adiabatic lapse rate in a drier atmosphere. The decrease in tropospheric temperature, in turn, is expected to emit less downward longwave radiation and reduce the longwave heating received by the surface (Fig. 4d). Also, the reduction in latent heat release implies anomalous subsidence that results in decreased cloud cover in most LULCC regions (Fig. S5b). Therefore, the surface receives more solar radiation, which tends to increase the LST as shown in Fig. 4c.
The DTM metric is applied in both degraded regions (Fig. 6a) and nondegraded regions (Fig. 6b) to quantify the contribution of each component to LST changes. The degraded area is further divided into lower latitudes (408S-408N) and midlatitudes (408-608N/S), as LST responds differently in these two areas. At lower latitudes, LULCC causes an increase in LST by 0.67 K, while at the midlatitudes a decrease in LST by 0.23 K is found (Fig. 6a). The opposite LST responses between lower and middle latitudes are also reported in observational LULCC studies (Zhang et al. 2014;Li et al. 2015;Alkama and Cescatti 2016;Bright et al. 2017;Duveiller et al. 2018). As shown in Fig. 6a, the LST response is mainly determined by the relative magnitude between the albedo effect (ABD) and the evaporation effect (LH), among which LH contribution has the largest discrepancy at lower and midlatitudes. The difference in LH contribution to LST changes between low latitudes and midlatitudes can be as large as 0.58 K, the reason for which will be discussed later.
The application of DTM to LST decrease in nondegraded regions is shown in Fig. 6b. The surface cooling (20.18 K) is dominated by the reduction in LW in , which is partly counterbalanced by the increased SW in . As discussed before, land degradation reduces surface evapotranspiration and therefore atmospheric water vapor (Fig. 5a), leading to cooler tropospheric temperature and a reduction in cloud cover ( Fig. 5b; see also Fig. S5b). As a result, negative DLW in and positive DSW in are found over the non-LULCC regions.
We further investigate the response of SH and LH and link them to the changes in biophysical properties due to LULCC. The SH changes can be decomposed to three major components: surface temperature (DLST), atmospheric temperature (DT a ), and aerodynamic resistance (Dr a ). The contribution of surface and atmospheric temperature changes is generally spread over the globe (Figs. 7b,c). As the sum of these two contributions produces a positive contribution to SH changes, an increase in the temperature gradient between surface and atmosphere is found, which tends to increase SH. The aerodynamic resistance increase due to lower surface roughness length and displacement height plays a dominant role in SH change and results in an overall SH decrease in degraded areas. The calculated SH changes from the decomposition components (Fig. 7e) reproduce the simulated SH changes (Fig. 7a) and the high-order terms account for less than 5% in most LULCC regions (Fig. 7f).
The decrease in LH is influenced by changes in LAI, soil moisture, aerodynamic resistance, and low-level moisture conditions (Nair et al. 2011;Wang et al. 2017). The change of LH is attributed to changes in surface vapor pressure (e s ), atmospheric vapor pressure (e s ), and surface resistance (r a 1 r c ) as shown in Fig. 8. The changes in surface vapor pressure caused by LULCC slightly increase LH in the tropical degraded region and decrease LH in other areas (Fig. 8b). Meanwhile, atmospheric water vapor decreases globally (De a , 0), which contributes to an increase of LH ({2rC p / [l(r a 1 r c )]}De a . 0; Fig. 8c). The vapor pressure gradient between surface and atmosphere is thus enhanced and causes an overall increase in LH. In contrast, the increase in surface resistance contributes to significant LH decrease over degraded regions (Fig. 8d). Since this process is dominant, LH decrease is found as a result of land degradation. The surface resistance (r a 1 r c ) for LH is closely related to soil moisture (Fig. S6), which controls the stomatal resistance. The changes in soil moisture are consistent with those in precipitation  (Fig. 3a). As in the SH analyses, the high-order terms are very small for LH decomposition (Fig. 8f).
The decomposition metric is applied to identify the important surface property that leads to different LH responses at lower and midlatitudes (Fig. 6a). As shown in Fig. 9, the change of LH is much larger at lower latitudes (29.84 W m 22 ) than at midlatitudes (25.54 W m 22 ). Our decomposition metric shows that the contribution from surface resistance to LH change is significantly different at lower latitudes (214.78 W m 22 ) and midlatitudes (26.44 W m 22 ). Surface resistance plays a dominant role in changing LH in the LULCC experiment, and the contributions from the surface and atmospheric vapor pressure are relatively small. It should be noted that the ensemble members have a similar result in the contribution of surface vapor pressure, atmosphere vapor pressure, and surface resistance to LH changes.
c. The effect of LULCC on precipitation and large-scale circulation from the perspective of meridional energy transport LULCC causes changes in radiation, turbulent fluxes, moisture, and temperature in the PBL and also indirectly influences cloud cover, convection, and precipitation. Figures 10a and 10b show the zonal-mean difference of precipitation, evapotranspiration, and moisture flux convergence (MFC) over land and ocean between LULCC and CTL during 1950-2015. The zonal-mean land precipitation changes are mostly consistent with evapotranspiration changes (Fig. 10a), indicating that the local effect in the degraded areas is dominant. When LULCC occurs, land covers with more vegetation and larger LAI are degraded to those with less vegetation and smaller LAI, and therefore the moisture supply for precipitation from evapotranspiration is decreased. As a result, precipitation is suppressed (Fig. 3a) and cloud cover is decreased mainly in the degraded areas (Fig. S5b). Over the Pacific Ocean in the northern midlatitudes (Fig. 10b), the precipitation change is affected by both evaporation and MFC changes. Ocean precipitation change around 408N is associated with evaporation decrease, in keeping with a cooling SST (Fig. 2a). In contrast, precipitation changes over the tropical oceans are mostly consistent with the MFC changes (Fig. 10b). LULCC causes an increase in precipitation south of the equator but decreases at the southern and northern boundaries of the ITCZ (Figs. 3a and 10b), which is associated with anomalous moisture flux convergence around the equator and divergence on both sides (Fig. 10c). The precipitation changes over tropical ocean indicate that the continental LULCC has a remote effect on tropical climate, mainly through changes in atmospheric moisture transport.
We investigate the LULCC impact on tropical precipitation using the energy flux framework in section 2c. The atmosphere energy balance and energy transport in CTL is shown in Fig. 11. The atmosphere has an energy surplus between about 358N/S (R ATM . 0 in Fig. 11a), which is transported out of the tropics FIG. 5. Differences of annual (a) specific humidity (g kg 21 ) and (b) temperature (K) between LULCC and CTL during 1950-2015 (LULCC 2 CTL). The contours represent the climatological mean from CTL and the shaded colors represent LULCC 2 CTL. Stippling indicates that the response is statistically significant at the 95% confidence level. through energy flux divergence (positive slope in Fig. 11b). Meanwhile, the atmospheric thermal cooling poleward of 358N/S (R ATM , 0 in Fig. 11a) necessitates a meridional energy flux convergence (negative slope in Fig. 11b). The negative F a (southward energy transport) at the equator in CTL indicates asymmetric atmospheric heating with respect to the two hemispheres in the climatology, associated with the shape of continents and atmosphere-ocean interactions (Xie 2004).
The occurrence of LULCC results in intensive cooling anomalies, as indicated by the negative R ATM difference in four latitude bands (Fig. 11a). Three of the cooling bands are located in the Northern Hemisphere, with one in the tropical region and the others in the extratropics. One intensive LULCC region is located in the Southern Hemisphere tropical region. There is also an atmospheric energy surplus in regions without LULCC (Fig. 11a), which is induced by climate feedbacks. The LULCC-induced R ATM change is further attributed to changes in R TOA , R SFC , LH, and SH, as shown in Fig. 11c. The anomalous R TOA and R SFC offset each other, indicating that reflected shortwave radiation from the surface travels through the atmosphere with little absorption and produces a rather small radiative effect on the atmosphere. In contrast, the anomalous latent heat flux from the surface to the atmosphere is consistent with the change of atmospheric energy (DR ATM ) at almost all latitudes, and the sensible heat flux plays a similar, albeit secondary, role. Our results emphasize the importance of nonradiative cooling of the atmosphere by latent and sensible heat fluxes in the PBL for the atmospheric energy balance, which is confirmed by the evaporation effect of LULCC (Figs. 4e,f).
Due to the nonradiative cooling in LULCC bands, an anomalous energy flux convergence is found at the LULCC latitudes (Fig. 11b), indicated by the negative slope of atmospheric energy flux. In regions without LULCC, the anomalous energy surplus (Fig. 11a) induced by climate feedbacks corresponds to an energy flux divergence (Fig. 11b). As the stronger atmospheric cooling occurs in the Northern Hemisphere, an anomalous northward energy transport of 0.05 petawatts (PW) takes place across the equator (Fig. 11b). According to the energy flux framework, the anomalous northward energy transport at the equator has been accomplished by the anomalous northward mass transport by the upper branch of Hadley circulation, necessitating a southward displacement of the rising branch of Hadley circulation (Fig. 3b). The southward shift of the rising branch of Hadley circulation relocates the precipitation centroid of the ITCZ southward, and the intensified equatorial updraft causes an increase in moisture flux convergence around the equator, producing more precipitation there (Figs. 3a and 10b).
A stronger descent is also found at the southern and northern boundaries of the Hadley circulation (Fig. 3b), which gives rise to the narrowing of the ITCZ as shown in Fig. 3a. The equatorward contraction is associated with the shift of the subtropical jets toward the equator as compared to the CTL (Fig. 12). Previous studies have suggested that the increase in upper-tropospheric temperature gradient associated with tropospheric warming under greenhouse increases can lead to a stronger jet that moves poleward (e.g., Chen and Held 2007;Chen et al. 2008), whereas the opposite is found for tropospheric cooling in response to more anthropogenic aerosols at the Northern Hemisphere midlatitudes (Ming and Ramaswamy 2009;Ming et al. 2011). The change of subtropical jet can be viewed from the perspective of the thermal wind equation: where (›T/›y) p is the meridional temperature gradient in the layer between pressure levels p 0 and p 1 , R is the gas constant for dry air, and f is the Coriolis force. In the LULCC experiment, the tropospheric cooling strengthens the meridional temperature gradient on the equatorward side of the jet while weakening it on the poleward side (Fig. 5b). According to the thermal wind relationship [Eq. (6)], the westerly jets are enhanced on the equatorward side and weakened on the polar side in both hemispheres, which, in turn, causes an equatorial contraction of the Hadley circulation due to the interaction between the subtropical jet and the descending branch of the Hadley cell. It narrows the Hadley circulation, and the resulting anomalous moisture flux divergence (Fig. 10c) causes a decrease in precipitation at the southern/ northern boundary of the ITCZ (Fig. 3a).

a. LST response in observational studies
Our simulations show that LULCC can cause surface warming in low latitudes and cooling in high latitudes due to the competing Observational studies provide independent information on the effects of deforestation on LST using satellite-retrieved products (Li et al. 2015;Alkama and Cescatti 2016;Duveiller et al. 2018) and semi-empirically derived LST based on the FLUXNET observations (Bright et al. 2017). By comparing the LST differences between forests and grasses, observational studies confirm our conclusion that LULCC has caused a warming in lower latitudes and a cooling in boreal regions, although the latitudes with cooling vary among datasets (Li et al. 2015;Alkama and Cescatti 2016;Bright et al. 2017;Duveiller et al. 2018). However, since these observation-based LULCC studies exclude the nonlocal LULCC effects in their methodologies and do not consider the conversions from grass/ shrub to bare ground that are included in our study, we are unable to make a quantitative comparison between simulated and observation-derived LST changes.
b. The regional and global effects of LULCC In this study, we investigate both regional/local and global/ nonlocal effects of LULCC during 1950-2015 using a coupled land-atmosphere-ocean Earth system model. Over the degraded regions, the ensemble GCM simulations show that LULCC has caused a warming effect by 0.67 K in lower latitudes and a cooling effect by 0.23 K in middle to high latitudes. The DTM analysis shows that the temperature response is dominated by changes in the latent heat and sensible heat fluxes in the tropical regions and by albedo changes in midlatitudes. We further conclude that the response of latent heat and sensible heat fluxes are dominated by surface roughness changes in a further developed DTM analysis. The precipitation is reduced due to the decrease in moisture supply from evapotranspiration associated with smaller LAI and vegetation fraction.
In the non-LULCC regions and oceans, we find a widespread cooling in surface temperature with a small seasonal variability throughout the year. Recent GCM studies on LULCC effects (Devaraju et al. 2018;Winckler et al. 2017Winckler et al. , 2019a also reported that nonlocal effects dominated the global mean temperature change or temperature change in a few regions. Using the developed DTM analysis, we are able to attribute the temperature decrease in the non-LULCC regions to the reduction in surface longwave heating corresponding to the drier and cooler troposphere after LULCC. A significant shrinking of ITCZ is simulated over tropical ocean associated with LULCC-induced tropospheric cooling and the resultant weakening and equatorward shift of the westerly jets. Meanwhile, a stronger nonradiative cooling in the Northern Hemisphere is caused by LULCC, which has been balanced by anomalous northward energy transport across the equator through the southward shift of Hadley cell. The responses of ITCZ have also been reported in a number of LULCC studies, which investigate the climate effects of idealized afforestation/deforestation in the Northern Hemisphere with slab ocean models (Swann et al. 2012;Devaraju et al. 2015;Laguë and Swann 2016). However, slab ocean models were found to be insufficient to simulate the global/nonlocal impact with projected thermal forcing (Deser et al. 2016;Tomas et al. 2016), as they overestimate the local thermodynamic impact of air-sea coupling and modify the spatial structure of the response compared to the fully coupled ocean model. Our study complements the previous LULCC studies by implementing annually updated LULCC processes in both hemispheres and by using a fully coupled ocean model, in which the oceans are dynamically and thermodynamically coupled to the atmosphere.

c. Uncertainties in LULCC simulation
Previous studies Pitman et al. 2009) have shown that simulated temperature and precipitation responses to LULCC have large discrepancies among GCMs. The uncertainties lie in several modeling aspects, such as different land surface schemes (Koster et al. 2006), different LULCC implementations (Boone et al. 2016), and also whether prescribed SSTs with interannual and seasonal variations, slab ocean models, or coupled ocean models were employed Pielke et al. 2011). Pielke et al. (2011) pointed out that the coupling strength between land and atmosphere could be one of the reasons that cause large differences in LULCC simulations. Koster et al. (2006) compared the land-atmosphere coupling strength using 12 GCMs with various land schemes. Our model (NCEP GFS/SSiB2) has a medium coupling strength, close to the ensemble mean, indicating that the land surface influence on the atmosphere in SSiB2 is within the range of current land surface models.
The uncertainties in simulated LULCC effects in this study come from several aspects. First, we apply the simplest LULCC implementation by changing the dominant vegetation type in one grid box in the LULCC map. Further implementation of vegetation tiling to improve the characterization of landscape complexity is needed. In addition, our simulations mainly consider deforestation and desertification while neglecting other types of LULCC such as irrigation, harvest, grazing, and urbanization. Furthermore, many countries have adopted better landuse management policies since the start of the twenty-first century, such as reforestation and afforestation (Hua et al. 2016) to mitigate climate changes, which have not been considered in our LULCC map. Although we have applied local and field significance tests, our conclusions are made based on a single model result. More analysis on LULCC effects, especially the nonlocal effects, should be carried out in the climate modeling community using land-atmosphere-ocean coupled models.

Conclusions
In this study, we investigated the global and regional effects of the LULCC in 1950-2015 on climate by conducting an ensemble of numerical simulations using the CFSv2/SSiB2. We incorporate the estimated land conversions by applying the LUH2 global land-use reconstructions from Hurtt et al. (2006Hurtt et al. ( , 2011. The vegetation fraction in the designed LULCC areas is degraded annually to a lower fraction, and by 2015 about 20% of global land areas have land degradation compared with CTL simulation. The GCM ensemble simulations show that LULCC results in a warming effect at lower latitudes and a cooling effect at midlatitudes in degraded regions. Meanwhile, the LST decreases in non-LULCC regions due to reduced longwave heating of the surface by a drier and cooled troposphere. Extensive LULCC activities in 1950-2015 are found to reduce the regional rainfall over LULCC regions and adjacent oceans. LULCC also narrows the ITCZ by about 0.58 and displaces it southward, forming a tripole precipitation anomaly around the FIG. 11. The 1950-2015 annual zonal-mean (a) atmospheric heating in CTL (dashed) and its difference between LULCC and CTL (solid). (b) The 1950-2015 atmospheric energy transport in CTL (dashed blue) and LULCC (dashed red) and its difference between LULCC and CTL (solid). The inset in the bottom right of (b) shows the enlarged box of atmosphere energy transport in CTL and LULCC near the equator. (c) The 1950-2015 annual-mean difference of zonal atmospheric heating (R ATM ; black) and its components: R TOA (red), 2R SFC (orange), LH (blue), and SH (green) between LULCC and CTL.
FIG. 12. The 1950-2015 annual-mean difference of zonal wind (m s 21 ) between LULCC and CTL (color shading) and climatological zonal wind in CTL (black contours). Stippling indicates statistical significance at 95% confidence level. warm pool. All ensemble members show a very similar magnitude of change, and both local and field significant tests are carried out to examine the statistical significance of the response.
The mechanism that causes a heterogeneous LST response in different latitudes is examined. The DTM analyses indicate that the temperature response is determined by the relative importance of increased albedo (cooling) and reduced evaporation (warming), and the magnitude of evaporation contribution decreases with latitude. We further decompose LH changes to those in surface vapor pressure, atmospheric vapor pressure, and surface resistance, and the results identify the dominant role of surface resistance in LH changes. The contribution of surface resistance to LH changes is significantly larger in the lower latitudes than in higher latitudes.
The reduction in local precipitation after LULCC results from the decrease in moisture supply from evapotranspiration associated with smaller LAI and vegetation fraction. The reduced nonradiative heating (SH and LH) from the surface to the atmosphere after LULCC leads to an asymmetric atmospheric cooling and results in a northward heat transport across the equator, accomplished by a southward shift of the rising branch of Hadley cell. The narrowing of the ITCZ is consistent with the equatorward contraction of the Hadley circulation, which is expected from the tropospheric cooling and the resultant weakening and equatorward shift of the westerly jet.
This study has also demonstrated that to appropriately assess the local and remote responses to LULCC, it is critical to use multiple methodologies. By using the DTM metric, we conclude that evaporation effects and albedo effects are the leading factors that dominate the local LULCC effects on surface temperature in the low latitudes and midlatitudes, respectively. The energy flux framework, which shows that LULCC has caused changes in atmospheric nonradiative heating due to the decrease of LH and SH from the surface, has been introduced to analyze the LULCC's remote effects on the ITCZ and Hadley circulation.