# Deforestation-induced warming over tropical mountain regions regulated by elevation

Zhenzhong Zeng, Dashan Wang, Long Yang, Jie Wu, Alan D. Ziegler, Maofeng Liu, Philippe Ciais, Timothy D. Searchinger, Zong-Liang Yang, Deliang Chen, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Zhenzhong Zeng, Dashan Wang, Long Yang, Jie Wu, Alan D. Ziegler, et al. Deforestation-induced warming over tropical mountain regions regulated by elevation. Nature Geoscience, 2021, 14, pp.23-29.
10.1038/s41561-020-00666-0 . insu-03660177

HAL Id: insu-03660177

https://insu.hal.science/insu-03660177v1

Submitted on 7 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Deforestation-induced warming over tropical mountain regions regulated by elevation

Zhenzhong Zeng<sup>1,2✉</sup>, Dashan Wang<sup>1</sup>, Long Yang<sup>3</sup>, Jie Wu<sup>1,4</sup>, Alan D. Ziegler<sup>1,5</sup>, Maofeng Liu<sup>1,2</sup>, Philippe Ciais<sup>1,6</sup>, Timothy D. Searchinger<sup>1,7</sup>, Zong-Liang Yang<sup>1,8</sup>, Deliang Chen<sup>1,9</sup>, Anping Chen<sup>1,10</sup>, Laurent Z. X. Li<sup>1,11</sup>, Shilong Piao<sup>1,12</sup>, David Taylor<sup>1,13</sup>, Xitian Cai<sup>1,14</sup>, Ming Pan<sup>1,2</sup>, Liqing Peng<sup>2</sup>, Peirong Lin<sup>1,2</sup>, Drew Gower<sup>1,2</sup>, Yu Feng<sup>1</sup>, Chunmiao Zheng<sup>1</sup>, Kaiyu Guan<sup>1,15</sup>, Xu Lian<sup>1,12</sup>, Tao Wang<sup>1,16</sup>, Lang Wang<sup>1,7</sup>, Su-Jong Jeong<sup>1,18</sup>, Zhongwang Wei<sup>1,19,20,21</sup>, Justin Sheffield<sup>2,22</sup>, Kelly Caylor<sup>1,23</sup> and Eric F. Wood<sup>1,2</sup>

Agriculture is expanding in tropical mountainous areas, yet its climatic effect is poorly understood. Here, we investigate how elevation regulates the biophysical climate impacts of deforestation over tropical mountainous areas by integrating satellite-observed forest cover changes into a high-resolution land-atmosphere coupled model. We show that recent forest conversion between 2000 and 2014 increased the regional warming by 0.022 ± 0.002 °C in the Southeast Asian Massif, 0.010 ± 0.007 °C in the Barisan Mountains (Maritime Southeast Asia), 0.042 ± 0.010 °C in the Serra da Espinhaço (South America) and 0.047 ± 0.008 °C in the Albertine Rift mountains (Africa) during the local dry season. The deforestation-driven local temperature anomaly can reach up to 2 °C where forest conversion is extensive. The warming from mountain deforestation depends on elevation, through the intertwined and opposing effects of increased albedo causing cooling and decreased evaporation causing warming. As the elevation increases, the albedo effect increases in importance and the warming effect decreases, analogous to previously highlighted decreases of deforestation-induced warming with increasing latitude. As most new croplands are encroaching lands at low to moderate elevations, deforestation produces higher warming from suppressed evapotranspiration. Impacts of this additional warming on crop yields, land degradation and biodiversity of nearby intact ecosystems should be incorporated into future assessments.

Mountains and hills are generally less economically viable and environmentally suitable for intensive agriculture than low-lying and non-sloping lands. Based on this premise, the United Nations Food and Agriculture Organization (FAO) projected little to no additional net cropland expansion within global mountains<fnref n="1" />, a projection that has been incorporated into future scenarios for climate change impact assessments²⁻<fnref n="4" />. However, the reality of the past decade deviates from this expectation for many sub-tropical and tropical mountains<fnref n="1" />,⁵. High-resolution forest cover datasets<fnref n="4" />⁻⁶ reveal that farmers have been carving a new agricultural frontier, causing high rates of

mountain forest loss in several areas⁷. Forest conversion to crops in the tropics affects climate by causing a decrease in evapotranspiration (warming) and an increase in albedo (cooling), leading to an overall local warming effect⁸⁻¹¹. This local warming, superimposed on regional climate change², may cause a greater threat to the region than expected. Globally, the climate feedbacks of deforestation are regulated by latitude, ranging from warming in the tropics to cooling in the mid- to high latitudes¹²⁻<fnref n="14" />. By analogy, we expect deforestation-induced warming over tropical mountain regions to also be regulated by elevation, as the evapotranspiration effects possibly attenuate at higher altitudes.

<footnote n="1">¹School of Environmental Science and Engineering, Southern University of Science and Technology, Shenzhen, China. ²Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ, USA. ³School of Geography and Ocean Science, Nanjing University, Nanjing, China.</footnote>

<footnote n="4">⁴Department of Geoscience and Natural Resource Management, University of Copenhagen, Copenhagen, Denmark. ⁵Faculty of Fisheries Technology and Aquatic Resources, Maejo University, Chiang Mai, Thailand. ⁶Laboratoire des Sciences du Climat et de l'Environnement, CEA/CNRS/UVSQ, Gif-sur-Yvette, France. ⁷Woodrow Wilson School, Princeton University, Princeton, NJ, USA. ⁸Department of Geological Sciences, Jackson School of Geosciences, University of Texas at Austin, Austin, TX, USA. ⁹Regional Climate Group, Department of Earth Sciences, University of Gothenburg, Gothenburg, Sweden.</footnote>

<footnote n="10">¹⁰Department of Biology, Colorado State University, Fort Collins, CO, USA. ¹¹Laboratoire de Météorologie Dynamique, Centre National de la Recherche Scientifique, Sorbonne Université, Ecole Normale Supérieure, Ecole Polytechnique, Paris, France. ¹²Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing, China. ¹³Geography Department, National University of Singapore, Singapore, Singapore.</footnote>

<footnote n="14">¹⁴Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA. ¹⁵Natural Resources and Environmental Sciences, University of Illinois at Urbana-Champaign, Urbana, IL, USA. ¹⁶Key Laboratory of Alpine Ecology and Biodiversity, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China. ¹⁷Department of Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong, China. ¹⁸Department of Environmental Planning, Graduate School of Environmental Studies, Seoul National University, Seoul, Republic of Korea. ¹⁹River and Environmental Engineering Laboratory, Department of Civil Engineering, University of Tokyo, Tokyo, Japan. ²⁰Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai, China. ²¹Guangdong Province Key Laboratory for Climate Change and Natural Disaster Studies, School of Atmospheric Sciences, Sun Yat-sen University, Guangzhou, China. ²²School of Geography and Environmental Sciences, University of Southampton, Southampton, UK. ²³Department of Geography, Bren School of Environmental Science and Management, University of California, Santa Barbara, Santa Barbara, CA, USA. ✉e-mail: zzeng@princeton.edu</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Simulated temperature change resulting from deforestation in tropical mountains. **a**, Spatial pattern of satellite-observed forest loss during 2000–2014 within 0.5°×0.5° grid cells and the inner domains for numerical simulations. **b–e**, Modelled temperature change at 5 km × 5 km grid cells triggered by forest cover change over SE1 (**b**), SE2 (**c**), SA (**d**) and AF (**e**). The effect was calculated as the difference between the simulations using the lower boundary with and without forest cover change. Dotting represents significance at *P* < 0.05.</figcaption></figure>

We addressed this research question over four tropical mountain regions (Fig. 1); namely, the Southeast Asian Massif (SE1), the Barisan Mountains in Maritime Southeast Asia (SE2), the Serra da Espinhaço in South America (SA) and the Albertine Rift mountains in Africa (AF). Supplementary Figs. 1 and 2 show the net forest loss in hectares and the percentage for each 5 km × 5 km grid cell during 2000–2014. Taking SE1 as an example, the total forest loss amounted to 3.2 million hectares, an area 20 times greater than that of Bangkok, the largest city within the region. Half of the forest loss occurred in the areas with low elevation concentrated around Cambodia, while the other half was distributed in the mountainous areas, mainly in Laos and northern Thailand around Nan Province (Supplementary Fig. 1). Earth system models and global climate models not only ignore mountain cropland expansion in their scenarios, but also have too coarse a resolution to represent topographic effects on climate²,⁴. For this reason, we coupled the satellite-observed, agriculturally driven forest conversion⁴ into a high-resolution regional climate model over the four regions experiencing cropland expansion. The model is a recent version of the Weather Research and Forecasting model coupled with the Noah Land Surface Model (WRF–Noah, version 3.8.1)¹⁵,¹⁶, which takes a mosaic approach to incorporating homogeneous, non-interacting sub-grid tiles¹⁷.

The grid spacing for our study areas was 5 km × 5 km, ensuring that the orographic variations were adequately represented (Supplementary Figs. 1 and 2). During the dry season, agricultural lands in many mountainous areas are bare ground (for example, Supplementary Fig. 3). Accordingly, the areas deforested during 2000–2014 were treated as cropland in the lowlands (<300 m) and as bare ground in the highlands (≥300 m) in the model based on the reported threshold to distinguish lowland and highland over

tropical mountains⁴,¹⁸. These simulations were performed during the local dry season for the year 2014 using lateral boundary conditions from ERA5 (the fifth-generation reanalysis of the European Centre for Medium-Range Weather Forecasts¹⁹) and forest cover in the year 2014 compared with a control run with forest cover in the year 2000 (refer to Supplementary Figs. 4–8 and the Methods for more details). We assessed the uncertainty associated with the climate background using simulations in other dry seasons during 2012–2016. We also assessed the uncertainty associated with the lateral boundary conditions using the simulations driven by FNL (the National Centers for Environmental Prediction Final Operational Global Analysis data)²⁰.

## Biophysical climate impacts of mountain deforestation

Consistent with reports of a net warming effect caused by tropical deforestation¹²,¹⁴, we found a general temperature increase triggered by agriculturally driven deforestation in the four mountain regions during 2000–2014 (Fig. 1b–e). In the SE1 region, the mean temperature increased significantly by 0.022 °C (95% confidence interval (CI)=0.020–0.024 °C; t-test statistics: t-value=22.16; degrees of freedom=88; *P* < 0.01). This increase in surface air temperature, if accounted for in the global climate models for the Intergovernmental Panel on Climate Change, would lead to an additional increase of 22.7% in regional warming in the twenty-first century under Representative Concentration Pathway 2.6 (RCP2.6), or 8.1% and 6.6% under RCP4.5 and RCP6.0, respectively. The pattern of temperature change is heterogeneous, primarily originating from highland deforestation in northern Thailand and northern Laos (Fig. 1b and Supplementary Fig. 1). For example, mountain forest conversion in a mountainous area of Nan Province of northern Thailand (–14.8%; for location and landscape information, see
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 2** | Temperature change versus percentage of forest loss for 2000–2014 at the grid level. **a–d** The temperature change at each grid cell was calculated using the simulated daily temperature during the local dry season from 1 December 2014 to 28 February 2015 over SE1 (**a**) and from 1 June to 31 August 2014 over SE2 (**b**), SA (**c**) and AF (**d**). The lowland grid cells (black crosses) and highland grid cells (red dots) are separated. The lines are the best fit of deforested grid cells by ordinary least squares regression.</figcaption></figure>

Supplementary Figs. 1a and 3) significantly increased the temperature by 0.17 °C (*P*<0.01; Supplementary Fig. 8). Similarly, the spatial patterns of temperature change match the distribution of forest loss in the SE2, SA and AF regions, in which the simulated temperature also increased significantly, up to 0.010 °C (95% CI=0.004–0.017 °C; *P*<0.01), 0.042 °C (95% CI=0.033–0.052 °C; *P*<0.01) and 0.047 °C (95% CI=0.039–0.055 °C; *P*<0.01), respectively (Fig. 1c–e and Supplementary Fig. 2). However, in some areas where little forest was converted, the effects of increased albedo overwhelmed those from the decreased evapotranspiration, leading to an overall local temperature decrease (for example, in parts of northern Laos and eastern Thailand; Fig. 1b).

The temperature change was positively correlated with the percentage of forest loss in the grid cells in both the highlands and the lowlands, but with different magnitudes (Fig. 2). The correlation coefficients ranged from 0.78–0.86 (*P*<0.01) in the highlands over all four regions, whereas the correlation coefficient was no more than 0.36 in the lowlands over the three regions except for AF. Moreover, the sensitivities of dry season temperature change to forest loss in the highlands and the lowlands were different. Taking SE1 as an example, there was an increase of 0.16 °C per 10% of forest loss in the highlands. Full forest loss in a grid cell (100% forest loss; hereafter, deforestation) therefore resulted in a temperature increase of 1.63 °C (95% CI=1.61–1.64 °C; *P*<0.01). However, in the lowlands, where the sensitivity was weaker, the temperature change from deforestation was only 0.04 °C (Fig. 2a). The sensitivities of temperature change to highland deforestation over SA and AF were 1.53 and 1.50 °C, respectively (Fig. 2c,d)—values that are close to that for SE1. In the coastal mountain region SE2, the sensitivity was smaller in comparison (0.88 °C; Fig. 2b).

Irrigation in lowland areas allows for intensive year-round agriculture, but crops cannot be grown intensively throughout an entire year at high elevations. Rain-fed croplands therefore have minimal surface vegetation during the dry season (for example, Supplementary Fig. 3). The contrast in biophysical properties between forest and bare ground in the highlands was greater than that between forest and quasi-permanent cropland in the lowlands, resulting in different sensitivities of temperature change to deforestation (Fig. 2). The temperature change caused by afforestation in the highlands depends on the type of land use change. In this study, forest gain was mainly the result of conversion from cropland/grassland. In this scenario, the opposing effects of decreased albedo causing warming can neutralize the increased evapotranspiration causing cooling, explaining the lack of sensitivity of temperature to forest gain (Fig. 2).

## Deforestation-induced warming regulated by elevation

Our analyses further showed an elevational effect on temperature change caused by agriculturally driven deforestation. Highland deforestation caused a net increase in local temperature that significantly weakened with increasing elevation (*P*<0.01; Fig. 3) in all four regions. In SE1, deforestation caused a highland heat island effect by increasing the local air temperature by 1.84±0.03 °C at an elevation of 350±50 m (*P*<0.01). As the elevation increased, the effect decreased to 1.56±0.02 °C at 750±50 m (*P*<0.01; inset in Fig. 3a), falling to 0.67±0.05 °C at 1,400±50 m (*P*<0.01). The effect became insignificant at elevations greater than 1,600 m, with few samples above this threshold (Fig. 3a).

The climatic impacts of mountain deforestation over the other three regions were also regulated by elevation, although the
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 3** | Regulation by elevation of the deforestation-induced warming over mountains. **a–d**, Results for SE1 (**a**; the total highland forest loss amount: 1.20×10⁶ ha), SE2 (**b**; 0.54×10⁶ ha), SA (**c**; 3.58×10⁶ ha) and AF (**d**; 1.81×10⁶ ha). For each elevation bin (±50 m), ordinary least squares regression was applied to regress the temperature change versus the forest loss fraction, defining a sensitivity of temperature change to deforestation. The inserts show the sensitivity of temperature change to deforestation in 750 ± 50 m elevation bins as examples. Error bars show the standard error of the sensitivity. The red dot size indicates the forest loss within each elevation bin.</figcaption></figure>

sensitivities varied (Fig. 3b–d). In SA, the effect was greatest at 450 ± 50 m (1.81 ± 0.01 °C; *P* < 0.01) and decreased at higher elevations. In AF, deforestation caused a temperature change of more than 1.48 °C between 450 and 1,000 m; it then decreased gradually (Fig. 3d), eventually becoming insignificant for elevations >3,000 m. The pattern over SE2 was generally similar to the other regions, but the magnitude of sensitivity was lower (Fig. 3b). Nevertheless, deforestation was more widespread in low-elevation mountains than at high elevations (Fig. 3). Taking SE1 as an example, 83% of highland forest area loss was at elevations ranging from 300 to 1,000 m, with 17% between 1,000 and 1,600 m and very limited forest loss above 1,600 m. Thus, the mean regional temperature increases were dominated by deforestation from low-elevation mountainous areas.

To examine further deforestation feedbacks with increasing elevation, we compared the albedo and evapotranspiration changes between the simulations with and without satellite-observed forest loss (Fig. 4 and Supplementary Fig. 10). Taking SE1 as an example, an increased area of forest loss resulted in a systematic increase in albedo and a decrease in evapotranspiration (Fig. 4a,b) for both lowland and highland areas. However, the sensitivities of albedo and evapotranspiration changes to percentage forest loss were greater in highlands compared with lowlands (Fig. 4a,b), explaining the contrasting temperature change induced by lowland versus highland deforestation (Fig. 2a). In the highlands at an elevation of 350 ± 50 m, we found that deforestation significantly increased

albedo by 0.249 ± 0.003 (P < 0.01; 56.43 ± 0.68 W m⁻²) and remarkably decreased evapotranspiration by 2.44 ± 0.04 mm d⁻¹ (P < 0.01; 70.61 ± 1.16 W m⁻²; Fig. 4c,d). The net change was an increase of 14.18 W m⁻² in the energy budget, implying that the warming induced by the reduction in evapotranspiration dominated the cooling by higher albedo values, resulting in a net warming effect.

As elevation increased, the deforestation-driven reduction in evapotranspiration significantly weakened (Fig. 4d) while the increase in albedo remained constant (Fig. 4c), explaining why deforestation-driven warming tends to decrease with elevation (Fig. 3a). The smaller evapotranspiration reduction from deforestation at higher elevations (Fig. 4d) could have been a result of the temperature lapse rate (0.72 °C decrease per 100 m elevation; Supplementary Fig. 11). Higher elevations lower the evaporative demand from temperature and water vapour pressure deficit, and thus, potential evapotranspiration, which dampens the sensitivity of evapotranspiration to forest cover losses. The deforestation-driven evapotranspiration and albedo change patterns over the other three regions were comparable with those from SE1 (Fig. 4 and Supplementary Fig. 10). In short, the overall temperature increase was mainly driven by deforestation in low-elevation mountains, where the evapotranspiration warming effect outweighed the albedo cooling effect.

To test the robustness of the elevation regulation concept, we focused on SE1 and repeated the analysis using simulations
<figure><img src="image_13.png" /><figcaption>**Fig. 4** | Mechanisms causing the elevation regulation of the deforestation-driven warming effect in SE1. **a,b,** Scatter plots of albedo change (a) and evapotranspiration (ET) change (b) versus percentage forest loss for 2000–2014 at the grid level, partitioned into lowland grids (black crosses) and highland grids (red dots). **c,d,** Sensitivity of albedo change (c) and evapotranspiration change (d) to deforestation against elevation. Error bars show the standard error of the sensitivity.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

during the four closest dry seasons to check whether variations in yearly monsoon and dry season characteristics influence the result (Supplementary Figs. 12 and 13). We also performed the analysis with FNL instead of the ERA5 lateral boundary conditions to test the influence of large-scale conditions outside the domain (Supplementary Figs. 14 and 15). For elevations lower than 1,500 m where >99.85% of deforestation occurred, our findings were robust regardless of the year of simulation (Supplementary Figs. 12 and 13) or lateral boundary conditions used (Supplementary Figs. 14 and 15). The results were more uncertain at elevations higher than 1,500 m where very little deforestation occurred (n < 30).

We further used the change in equivalent temperature ($T_E$) instead of surface air temperature to examine the elevation-dependent warming effects of mountain deforestation on surface air moist enthalpy. By including both dry and moist terms, $T_E$ provides a good indicator of changes in the near-surface total heat content. For example, the relationship between $T_E$ change and forest loss in the SE1 highlands was overall consistent with the deforestation-induced local surface air temperature change (Supplementary Fig. 16a). $T_E$ increased with the percentage of forest loss in the highlands, whereas in the lowlands $T_E$ decreased instead, an effect that may be related to a decrease in specific humidity. The $T_E$ warming effect of highland deforestation significantly attenuated with elevation (Supplementary Fig. 16b), confirming the regulating role of elevation on the biophysical climate impacts of mountain deforestation (Fig. 3).

## Comparison of the regulation by elevation versus latitude

Focusing on deforestation-prone tropical mountain regions$^{4,21}$, we demonstrated an attenuation of deforestation-induced warming with elevation. This is caused by deforestation-led evapotranspiration decline with rising elevation, along with the compensation by a constant albedo cooling effect. The intertwined and opposing biophysical climate feedbacks of deforestation result in an elevation-dependent overall local warming. This diagnosis of elevation-dependent biophysical climate impacts of deforestation is analogous to the latitudinal variation of the climate effects from forest cover change$^{12,14}$. Nonetheless, forest cover change in mid- to high latitudes is dominated by afforestation$^{21,22}$.

Previous research has revealed a critical latitude of 45°N that delineates the crossover of albedo-induced warming versus evapotranspiration-induced cooling for afforestation climatic feedbacks$^{12,22,23}$. For latitudes above this limit, the increased shortwave radiation and lower albedo increase from more leaf area result in a net warming effect. However, here, the critical elevation at which the sign of tropical mountain deforestation-related climate feedbacks reverses is unknown. This may be due to an absence of snow in tropical mountain regions. In high-latitude mountainous areas, where snow (and its albedo effect) is frequently present, trees generally remain snow free on snowy days, leading to lower albedo in forested landscapes compared with open land. Furthermore, even if snow may cover some of the forest canopy, shortwave radiation is scattered rather than reflected. Hence, deforestation in boreal mountains, but not in snow-free tropical mountains, may incur an
additional albedo cooling effect that would offset warming from reduced evapotranspiration.

## Implications

The deforestation-induced highland heat island effect impacts not only deforested areas but also the surrounding areas within a radius of 50 km (ref. ²⁴). This influence has several implications for regional environmental risks, including biodiversity loss, human and livestock health threats and lower crop yields. The SE1, SE2, SA and AF regions fall within the Indo-Burma biodiversity hotspot, the Sundaland biodiversity hotspot. Brazil's Cerrado biodiversity hotspot and the Eastern Afromontane biodiversity hotspot, respectively²⁵,²⁶, which harbour many endemic species that are highly sensitive to climate change²⁷,²⁸. If mountain deforestation continues, endemic species in these regions will probably suffer greater extinction risks with the combined effects of habitat loss/degradation and concomitant warming²⁹,³⁰. In addition, local warming may increase the risk of extreme heat that threatens human health³¹, and worsen the living conditions of rural communities located at lower mountain elevations³²,³³. Higher temperatures potentially increase the risks related to fire, as well as decreased crop quality and increased harvest failure³⁴. Moreover, regional warming may reduce thermal barriers to the dispersal of insects and the potential pathogens they carry to highland areas, increasing the exposure of people, crops and livestock to greater risks of infection³⁵. Meanwhile, as the temperature increases, insect pests consume more crops, further impacting agricultural yields³⁶.

Our findings draw attention to the need for sustainable land management strategies to reduce degradation in mountains for climate change mitigation. Such strategies need to reconcile conservation with development. One feasible way, as recommended by the Intergovernmental Panel on Climate Change³⁷, is to consider sustainability criteria in the global trade of land and land-based commodities. Another complementary bottom-up strategy, which is currently being considered in Nan Province, Thailand, is: (1) intensify agriculture in lower-elevation hills where irrigation can be employed (around 22% of the world's mountainous areas, as estimated by the FAO³⁸); and (2) in tandem, reforest/rehabilitate previously converted areas at higher elevations. In support, the Nan Province government is planning to build 14 reservoirs over the next 10 years, partly to improve agriculture efficiency (Supplementary Fig. 17). These changes will not only conserve upland forests and the carbon they store, but will also reduce local warming in low-elevation mountains by enabling land to remain vegetated throughout the dry season.

In summary, this study shows that there is a profound heat island in deforestation-associated highlands and highlights the importance of using more realistic projections of land use change in future assessments. However, caveats and limitations remain. First, our work has not considered the biophysical climate effects of burned areas. Many farmers set fire to clear fields of agricultural waste. These fires directly increase local temperature by burning, and indirectly increase local temperature via reducing the albedo cooling effect as black charcoal from biomass burning absorbs more solar radiation³⁹. Second, hillslope redistribution and lateral flow of groundwater can feedback to the atmosphere; thus, forest-to-crop conversion in hillslopes could involve multiple groundwater-atmosphere interacting processes⁴⁰. Yet, the current models have not incorporated lateral groundwater flow. Another caveat is the indirect effects of clouds. Although clouds are parameterized in models, their effects are largely unclear. In addition, soil erosion can decrease the water holding capacity of soil, which in turn amplifies evapotranspiration reduction caused by forest conversion. To improve the projections of land use change impacts, the next generation of high-resolution regional climate models should include these processes.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-020-00666-0.

Received: 2 April 2020; Accepted: 11 November 2020; Published online: 14 December 2020

## References

1. Chini, L. P., Hurtt, G. C. & Frolking, S. Harmonized Global Land Use for Years 1500–2100, V1 (ORNL DAAC, 2014).

2. IPCC Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) (Cambridge Univ. Press, 2013).

3. Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. *Science* **342**, 850–853 (2013).

4. Zeng, Z. et al. Highland cropland expansion and forest loss in Southeast Asia in the twenty-first century. *Nat. Geosci.* **11**, 556–562 (2018).

5. Zeng, Z., Gower, D. & Wood, E. F. Accelerating forest loss in Southeast Asian Massif in the 21st century: a case study in Nan Province, Thailand. *Glob. Change Biol.* **24**, 4682–4695 (2018).

6. Alexandratos, N. & Bruinsma, J. World Agriculture Towards 2030/2050: The 2012 Revision Working Paper No. 12-03 (FAO, 2012).

7. Grogan, K., Pflugmacher, D., Hostert, P., Mertz, O. & Fensholt, R. Unravelling the link between global rubber price and tropical deforestation in Cambodia. *Nat. Plants* **5**, 47–53 (2019).

8. Mahmood, R. et al. Land cover changes and their biogeophysical effects on climate. *Int. J. Climatol.* **34**, 929–953 (2014).

9. Lawton, R. O., Nair, U. S., Pielke, R. A. Sr & Welch, R. M. Climatic impact of tropical lowland deforestation on nearby montane cloud forests. *Science* **294**, 584–587 (2001).

10. Ray, D. K., Nair, U. S., Lawton, R. O., Welch, R. M. & Pielke, R. A. Sr Impact of land use on Costa Rican tropical montane cloud forests: sensitivity of orographic cloud formation to deforestation in the plains. *J. Geophys. Res.* **111**, D02108 (2006).

11. Lawrence, D. & Vandecar, K. Effects of tropical deforestation on climate and agriculture. *Nat. Clim. Change* **5**, 27–36 (2015).

12. Bonan, G. B. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. *Science* **320**, 1444–1449 (2008).

13. Lee, X. et al. Observed increase in local cooling effect of deforestation at higher latitudes. *Nature* **479**, 384–387 (2011).

14. Alkama, R. & Cescatti, A. Biophysical climate impacts of recent changes in global forest cover. *Science* **351**, 600–604 (2016).

15. Chen, F. The Noah Land Surface Model in WRF: A Short Tutorial (RAL, TIIMES, NCAR, 2007); https://go.nature.com/3lcWcsS

16. Skamarock, W. C. et al. A Description of the Advanced Research WRF Version 3. NCAR Tech. Note NCAR/TN-475+STR (UCAR, 2008).

17. Li, D., Bouzeid, E., Barlage, M., Chen, F. & Smith, J. A. Development and evaluation of a mosaic approach in the WRF–Noah framework. *J. Geophys. Res. Atmos.* **118**, 11918–11935 (2013).

18. Scott, J. C. The Art of Not Being Governed: An Anarchist History of Upland Southeast Asia (Yale Univ. Press, 2014).

19. Hersbach, H. & Dee, D. ERA5 reanalysis is in production. *ECMWF Newsl.* **147**, 7 (2016).

20. NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids (National Centers for Environmental Prediction, National Weather Service, NOAA & US Department of Commerce, accessed 1 August 2018); https://doi.org/10.5065/D65Q4T4Z

21. Song, X. P. et al. Global land change from 1982 to 2016. *Nature* **560**, 639–643 (2018).

22. Betts, R. A. Offset of the potential carbon sink from boreal forestation by decreases in surface albedo. *Nature* **408**, 187–190 (2000).

23. Peng, S. S. et al. Afforestation in China cools local land surface temperature. *Proc. Natl Acad. Sci. USA* **111**, 2915–2919 (2014).

24. Cohn, A. S. et al. Forest loss in Brazil increases maximum temperatures within 50 km. *Environ. Res. Lett.* **14**, 084047 (2019).

25. Myers, N., Mittermeier, R. A., Mittermeier, C. G., Fonseca, G. A. & Kent, J. Biodiversity hotspots for conservation priorities. *Nature* **403**, 853–858 (2000).

26. Mittermeier, R. A. et al. in *Biodiversity Hotspots* (eds Zachos, F. E. & Habel, J. C.) 3–22 (Springer, 2011).

27. Freeman, B. G. & Freeman, A. M. Rapid upslope shifts in new Guinean birds illustrate strong distributional responses of tropical montane species to global warming. *Proc. Natl Acad. Sci. USA* **111**, 4490–4494 (2014).

28. Freeman, B. G., Scholer, M. N., Ruizgutierrez, V. & Fitzpatrick, J. W. Climate change causes upslope shifts and mountaintop extirpations in a tropical bird community. *Proc. Natl Acad. Sci. USA* **115**, 11982–11987 (2018).
29. Thomas, C. D. et al. Extinction risk from climate change. *Nature* **427**, 145–148 (2004).

30. Sala, O. E. et al. Global biodiversity scenarios for the year 2100. *Science* **287**, 1770–1774 (2000).

31. Mora, C. et al. Global risk of deadly heat. *Nat. Clim. Change* **7**, 501–506 (2017).

32. Kohler, T. & Maselli, D. *Mountains and Climate Change: From Understanding to Action* (Geographica Bernensia, 2009).

33. Huber, U. M. et al. *Global Change and Mountain Regions: An Overview of Current Knowledge* (Springer Science & Business Media, 2006).

34. Beniston, M. in *Climate Variability and Change in High Elevation Regions: Past, Present & Future* (ed. Diaz, H. F.) 5–31 (Springer, 2003).

35. *Practical Chemotherapy of Malaria: Report of a WHO Scientific Group* (WHO, 1990).

36. Deutsch, C. et al. Increase in crop losses to insect pests in a warming climate. *Science* **361**, 916–919 (2018).

37. *Special Report on Climate Change and Land* (IPCC, 2019).

38. Egan, P. A. & Price, M. F. *Mountain Ecosystem Services and Climate Change: A Global Overview of Potential Threats and Strategies for Adaptation* (UNESCO, 2017).

39. Liu, Z., Ballantyne, A. P. & Cooper, L. A. Biophysical feedback of global forest fires on surface temperature. *Nat. Commun.* **10**, 214 (2019).

40. Forrester, M. M. & Maxwell, R. Impact of lateral groundwater flow and subsurface lower boundary conditions on atmospheric boundary layer development over complex terrain. *J. Hydrometeorol.* **21**, 1133–1160 (2020).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2020
## Methods

**Coupled Land-Atmosphere Regional Weather Model.** The WRF model, developed by the National Center for Atmospheric Research, is a fully compressible, nonhydrostatic, mesoscale model that has been widely used for regional climate studies$^{41–45}$. The Advanced Research version of the WRF model (ARW, version 3.8.1) was used in this study. We coupled the Noah Land Surface Model into the WRF model to accurately represent interactions between the land surface and the lower atmosphere$^{5,16}$. The mosaic approach was adopted in the Noah Land Surface Model to account for heterogeneities over sub-grid scales (Noah Mosaic$^{17}$). Instead of using the dominant land use over the grid, the mosaic approach considers surface fluxes over N categories (N=3) of land use type, with weightings based on fractional coverage. Sub-grid heterogeneity was of prime importance in the analyses as deforestation was rather sporadic over the regions studied and may not have changed the dominant land use type over the entire grid$^{6}$.

All simulations in this study shared the same physical options (Supplementary Table 1). The Kain-Fritsch scheme$^{46}$ was applied to simulate the effects of convective clouds at a sub-grid scale on the grid variables (for example, precipitation and temperature). The WRF single-moment 6-class microphysics scheme$^{47}$ was used to simulate ice microphysical processes more realistically. Longwave radiation was simulated using the Rapid Radiative Transfer Model, in which longwave fluxes and heating are efficiently and accurately calculated via the correlated-k approach$^{48}$. Shortwave radiation was simulated by the Dudhia scheme$^{49}$, which accounted for the effects of terrain slope and shadowing on the surface solar flux. Land surface turbulent fluxes were simulated by a widely used framework based on the Monin-Obukhov similarity theory$^{50}$. The Yonsei University planetary boundary layer scheme was used to parameterize the turbulent vertical fluxes of heat, constituents and momentum in the planetary boundary layer and throughout the atmosphere$^{51}$.

**Experimental design.** We investigated the biophysical climate impacts of agriculturally driven forest loss over four tropical mountain regions that are undergoing rapid deforestation in the twenty-first century (that is, SE1, SE2, SA and AF). Our analyses were constrained to the tropical mountain regions because recent deforestation and agricultural expansion dominate in the tropics$^{1,21}$ (Fig. 1). We designed two one-way nested domains with horizontal grid spacings of 25 and 5 km, respectively (Supplementary Fig. 4). The simulations over the four study areas shared the same configurations except for the domain dimensions, as summarized in Supplementary Table 1. We used two reanalysis products (ERA5$^{19}$ and FNL$^{20}$) for initial and lateral conditions for the numerical simulations. To better isolate the signal of the deforestation-driven climate impacts, we simulated the climate during the local dry season when the climate noise caused by external forcing was small (that is, from 15 November to 28 February in the SE1, and from 16 May to 31 August over the other three areas). We performed simulations for the year 2014 over the four study areas and treated the first 16 d of each model run as spin up. Our results are based on the simulations from 1 December to 28 February for SE1, and from 1 June to 31 August for the other areas.

We implemented two sets of numerical experiments with contrasting land use and land cover. For the control experiments, we used the Moderate Resolution Imaging Spectroradiometer (MODIS) 30-s land cover product$^{52}$ in both model domains to represent the current landscape features. MODIS 30-s provides global maps at 1-km spatial resolution of the Noah-modified 20-category International Geosphere-Biosphere Programme land cover classifications, the dominant land use category, 16-category top-layer soil type, the dominant soil category, topography height and sub-grid scale orography$^{53}$. We used these static maps to represent the land cover for the year 2000. The 20 land use categories comprise evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, mixed forests, closed shrublands, open shrublands, woody savannas, savannas, grasslands, permanent wetlands, croplands, urban and built up, cropland/natural vegetation mosaic, snow and ice, barren or sparsely vegetated, water, wooded tundra, mixed tundra and barren tundra. Values of parameters describing land surface characteristics (for example, albedo, emissivity, fractional vegetative cover, leaf area index (LAI) and surface roughness) were determined for the categories of land use type in each grid point as a function of the land use category via the default lookup tables (for example, GENPARM.TBL, MPTABLE.TBL, VEGPARM.TBL and SOILPARM.TBL). LAI parameters were updated for the tropical ecosystems using the MODIS Terra and Aqua combined LAI product$^{53}$ (Supplementary Table 2).

For other (simulated) experiments, land use over the inner domain was replaced with a modified land use dataset to allow the difference between the simulated and control to accurately represent the forest loss during 2000–2014 in the study regions. Most satellite-based products of land cover change, such as the European Space Agency Climate Change Initiative$^{54}$, GlobeLand30 (from the National Geomatics Center of China)$^{55}$ and MCD12Q1 (from the National Aeronautics and Space Administration of the United States)$^{56}$, failed to capture the large-scale forest loss in Southeast Asia during the early twenty-first century$^{4}$. Therefore, we did not directly use the land cover maps in 2014 from these products. Instead, we superimposed the forest cover change of the high-resolution global maps of twenty-first-century forest cover (hereafter referred to as

HANSEN$^{2)}$ between 2000 and 2014 onto the original MODIS 30-s land cover product; the modified maps represent the land cover for 2014. According to the validation work based on FAO statistics, LiDAR detection and other satellite measurements, the HANSEN product shows an overall accuracy greater than 99% at the global scale$^{1}$. In previous work$^{4}$, we validated the accuracy of HANSEN forest loss pixels for the lowlands and uplands in Southeast Asia using the latest high-resolution satellite imagery (for example, RapidEye (5 m), Doves (3 m), EarlyBird-1 (0.8–3.0 m), IKONOS (3.2 m), QuickBird (0.6–2.4 m), Pleiades-1A (0.5–2.0 m) and Pleiades-1B (0.5–2.0 m)), further strengthening the accuracy assessment of the HANSEN dataset.

To generate the modified land use map for 2014, we first aggregated the HANSEN-derived net forest change (that is, gain minus loss) at 30-m resolution into total forest cover change for each 5 km × 5 km grid cell. The total forest cover change was then converted to the percentage forest cover change based on the land area of each 5 km × 5 km grid cell. In the lowland areas in tropical regions, irrigation allowed for intensive year-round cropping, while the rain-fed, highland croplands could not be grown throughout an entire year and had minimal surface vegetation during the dry season (for example, Supplementary Fig. 3). The elevation threshold of lowland and highland cropland was around 300 m (refs. $^{4,18}$). Therefore, the HANSEN-derived forest cover changes during 2000–2014 were superimposed onto the original MODIS 30-s land cover product to represent the satellite-observed deforestation, in which we increased the percentage of cropland (bare ground) in lowland (highland) accordingly. If there was net forest loss in the grid cell, we decreased the percentages of forest-type categories of land use proportionally (that is, evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest and mixed forests), and we increased the percentage of croplands if the elevation was lower than 300 m, or that of bare ground (that is, barren or sparsely vegetated ground) if the elevation was higher than 300 m. Inversely, if there was net forest gain in the grid cell, we then proportionally increased the percentages of forest-type categories of land use and proportionally decreased the percentages of non-forest-type categories of land use (excluding water, wooded tundra, mixed tundra and barren tundra).

The two sets of experiments (control and simulated) used the same model configurations as introduced above but were driven by the lower boundary conditions with and without the satellite-observed forest loss. Thus, the differences between the simulated and control experiments were used to elucidate the biophysical climate impacts of mountain deforestation over the four tropical regions.

**Model performance.** After using the first 16 d of each model run as spin up time, the remaining periods (that is, 1 December 2014 to 28 February 2015 in SE1 and 1 June to 31 August 2014 over the other three areas) were used for validation and analyses. The model performance was first evaluated by comparing the simulated experiments with the surface air temperature from in situ weather stations in the four tropical mountainous areas (download on 15 April 2020 from ftp://ftp.ncdc.noaa.gov/pub/data/gsod). There were 125 Global Surface Summary of the Day (GSOD) ground-observed weather stations in SE1, while the stations were relatively sparse over the other three mountain regions (30, five and nine stations in total for SE2, SA and AF, respectively). It is noteworthy that the density and quality of ground observations may have added uncertainties to the results of model validations. For SE1 and SE2 in Southeast Asia, we also evaluated the model performance using a gridded surface air temperature dataset reconstructed from a dense network of daily gauge data around Asia (APHRODITE)$^{57}$.

In SE1, 40 of the 125 GSOD stations (32%) were located in the highlands (Supplementary Fig. 5b); this proportion was in line with the overall ratio of highland areas (41%; Supplementary Fig. 1). The model accurately simulated the spatial-temporal variability of mean daily surface air temperature (Supplementary Fig. 5). Compared with the ground-observed 2-m air temperature records from 125 weather stations (dots in the inset plot of Supplementary Fig. 5b), the WRF-Noah model accurately simulated the spatial-temporal variability of near-surface air temperature in the study region, regardless of which reanalysis product was used for initial and lateral boundary conditions (temporal variability: $R > 0.96$; $P < 0.01$ (Supplementary Fig. 5a); spatial variability: $R > 0.93$; $P < 0.01$ (Supplementary Fig. 5b)). The dynamic downscaling of reanalysis by the WRF-Noah model was almost as accurate as the APHRODITE product (green line in Supplementary Fig. 5a). In the other three tropical mountain regions, the model also generally captured the temporal variability of mean daily surface air temperature ($R = 0.83$, 0.72 and 0.66 for SE2, SA and AF, respectively; $P < 0.01$; Supplementary Fig. 6). Besides, our simulations provided finer spatial information in mountain regions (Supplementary Figs. 7 and 8). Because of the similarity of the results driven by the two sets of initial and lateral boundary conditions, we report the results of the simulations driven by the ERA5 lateral boundary conditions, and use the simulations driven by FNL as a test for robustness.

**Calculation of $T_E$.** As surface air temperature alone does not capture the change of heat content in the atmosphere, we calculated $T_E$, which represents the surface air moist enthalpy and contains both dry and moist heat content$^{38}$. In each grid cell, the daily $T_E$ (in °C) is given by:

$$T_E = T + T_M \quad (1)$$
where *T* is the simulated daily mean air temperature (°C) and *T*<sub>M</sub> is the moisture heat content, which is given by:

$$T_M = \lambda q/C_p = 1000 \times (2.501 - 0.002361 T)q/C_p \quad (2)$$

where λ is the latent heat of vapourization (MJ kg⁻¹), q is the simulated daily mean specific humidity (kg kg⁻¹) and C<sub>p</sub> (=1.013 kJ kg⁻¹ °C⁻¹) is the specific heat at constant pressure for air.

## Data availability

Data on satellite-observed high-resolution forest cover change in the twenty-first century are available at http://earthenginepartners.appspot.com/science-2013-global-forest. GSOD surface air temperature data are available at ftp://ftp.ncdc.noaa.gov/pub/data/gsod. The ERA5 reanalysis product is available at https://cds.climate.copernicus.eu/. The FNL reanalysis product is available at https://rda.ucar.edu/datasets/ds083.2. All of the datasets are also available on request from Z.Z.

## Code availability

We used the programmes MATLAB (R2014a) and ArcGIS (10.4) to generate all of the results. Analysis scripts are available at https://doi.org/10.6084/m9.fishshare.13280150.

## References

41. Leung, L. R., Kuo, Y. H. & Tribbia, J. Research needs and directions of regional climate modeling using WRF and CCSM. *Bull. Am. Meteorol. Soc.* **87**, 1747–1751 (2006).

42. Bukovsky, M. S. & Karoly, D. J. Precipitation simulations using WRF as a nested regional climate model. *J. Appl. Meteorol. Clim.* **48**, 2152–2159 (2009).

43. Caldwell, P., Chin, H. N. S., Bader, D. C. & Bala, G. Evaluation of a WRF dynamical downscaling simulation over California. *Clim. Change* **95**, 499–521 (2009).

44. Wang, J. L. & Kotamarthi, V. R. Downscaling with a nested regional climate model in near-surface fields over the contiguous United States. *J. Geophys. Res. Atmos.* **119**, 8778–8797 (2014).

45. Wi, S. et al. Climate change projection of snowfall in the Colorado River basin using dynamical downscaling. *Water Resour. Res.* **48**, W05504 (2012).

46. Kain, J. S. The Kain-Fritsch convective parameterization: an update. *J. Appl. Meteorol. Clim.* **43**, 170–181 (2004).

47. Hong, S. Y. & Lim, J. O. J. The WRF single-moment 6-class microphysics scheme (WSM6). *Asia Pac. J. Atmos. Sci.* **42**, 129–151 (2006).

48. Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J. & Clough, S. A. Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *J. Geophys. Res.* **102**, 16663–16682 (1997).

49. Dudhia, J. Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model. *J. Atmos. Sci.* **46**, 3077–3107 (1989).

50. Monin, A. S. & Obukhov, A. M. Basic laws of turbulent mixing in the surface layer of the atmosphere. *Tr. Akad. Nauk SSSR Geophiz. Inst.* **24**, 163–187 (1954).

51. Hong, S.-Y., Noh, Y. & Dudhia, J. A new vertical diffusion package with explicit treatment of entrainment processes. *Mon. Weather Rev.* **134**, 2318–2341 (2006).

52. Friedl, M. A. et al. Global land cover mapping from MODIS: algorithms and early results. *Remote Sens. Environ.* **83**, 287–302 (2002).

53. Chen, Y. et al. Generation and evaluation of LAI and FPAR products from Himawari-8 Advanced Himawari Imager (AHI) data. *Remote Sens.* **11**, 1517 (2019).

54. Land Cover CCI: Product User Guide Version 2.0 (ESA, 2017); https://go.nature.com/3nTzp6Y

55. Chen, J. et al. Global land cover mapping at 30 m resolution: a POK-based operational approach. *ISPRS J. Photogramm. Remote Sens.* **103**, 7–27 (2015).

56. Friedl, A. et al. MODIS Collection 5 global land cover: algorithm refinements and characterization of new datasets. *Remote Sens. Environ.* **114**, 168–182 (2010).

57. Yasutomi, N., Hamada, A. & Yatagai, A. Development of a long-term daily gridded temperature dataset and its application to rain/snow discrimination of daily precipitation. *Glob. Environ. Res.* **15**, 165–172 (2011).

58. Pielke, R. A. Sr, Davey, C. & Morgan, J. Assessing “global warming” with surface heat content. *Eos* **85**, 210–211 (2004).

## Acknowledgements

This study was supported by Lamsam-Thailand Sustain Development (B0891), the National Natural Science Foundation of China (42071022, 42001321), the China Postdoctoral Science Foundation (2020M672693), the start-up fund provided by the Southern University of Science and Technology (29/Y01296122) and the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA20060401). P.C. acknowledges support from the European Research Council Synergy project SyG-2013-610028 IMBALANCE-P and the ANR CLAND Convergence Institute. D.C. was supported by Swedish BECC and MERGE. We thank Della Research Computing at Princeton University and the Taiyi Supercomputer at the Southern University of Science and Technology for providing computing resources. We sincerely appreciate D. S. Wilcove for constructive comments on this paper.

## Author contributions

Z.Z. designed the research and wrote the draft. Z.Z. and D.W. performed the analysis.

Z.Z., D.W., L.Y. and M.L. performed the numerical simulations. All of the authors contributed to interpretation of the results and writing of the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41561-020-00666-0.

Correspondence and requests for materials should be addressed to Z.Z.

Peer review information Primary Handling Editors: Tamara Goldin; Xujia Jiang.

Reprints and permissions information is available at www.nature.com/reprints.