Introduction

Forests have the capacity to buffer global warming by storing large amounts of carbon from the atmosphere via photosynthesis1,2,3. Alongside the biochemical effects, forests can influence the local and regional climate through biophysical processes, including alterations in land surface water and energy balance4,5,6,7. On the local scale, the higher albedo and lower evapotranspiration (ET) following deforestation cause either surface cooling or warming, depending on which process holds dominance8,9,10. These cooling or warming impacts have the potential to offset or intensify, respectively, the warming effects connected to the released carbon caused by deforestation11,12,13,14,15,16. Extensive studies on the direct biophysical effects of deforestation on surface temperature have unveiled a latitudinal shift from tropical warming to boreal cooling8,9,17,18,19. Nevertheless, globally, alterations in surface albedo are more prevalent in the direct biophysical temperature response than ET because of its wider-scale impact17. This suggests that the global warming attributed to the biochemical effects of deforestation could potentially be mitigated by the cooling effects resulting from increased surface albedo and consequently altered radiative balance12,14,16. Yet, the impact of forest indirect biophysical processes on clouds and their associated radiative balance has not been well addressed, and the assessment of how changes in cloud radiative effects interact with the surface albedo effects remains unquantified. Understanding the response of clouds and their radiative effects to deforestation, however, is crucial due to the overwhelming effect clouds play on the Earth's energy budget. It stands as a major challenge in evaluating land-use-change-driven climate change20,21,22,23,24.

Observational studies allow for the conclusion that deforestation may predominantly reduce global cloud cover22,23,25, but with contrasting impacts across various regions21. These studies mostly compare clouds above forests and open land in adjacent geographical units (i.e., space-for-time substitution) and find larger cloudiness over forests. This commonly adopted method assumes that forests and neighboring land units share the same climate background, thereby deducing local effects through distinctions in land surface conditions. Apart from observations-based studies, general circulation models (GCMs) have been widely employed to quantify the impacts of deforestation26,27,28. GCMs show a global average enhancement in cloud cover with deforestation29. Unlike the observational studies that concentrate solely on local effects, GCMs probably possess the ability to encompass both local and non-local effects of deforestation. Hence, separating local and non-local effects could facilitate comparisons between these two distinct methods and enhance comprehension of the biophysical mechanisms of deforestation on clouds24.

Given the essential roles of cloud vertical structures in influencing radiative processes30,31, a sole concentration on overall cloud cover may be insufficient for a comprehensive analysis of the changes in cloud radiative effects from deforestation. Typically, low, highly reflective clouds have a cooling effect as they reflect solar radiation. In contrast, high, semi-transparent clouds contribute to warming by allowing shortwave radiation to pass through while impeding longwave radiation32,33. The alterations in cloud vertical profiles following deforestation have not received adequate attention, and addressing this gap is essential for gaining a deeper understanding of the consequent changes in cloud radiative effects.

In this study, we approach the evaluations of cloud profiles and associated radiative response to deforestation from two distinct viewpoints: the space-for-time substitution method from observations-driven reanalysis and the idealized deforestation experiments available from GCM simulations. Given that the outcomes from GCMs contain both local and non-local signals, we then isolate the local signals using a chessboard-like method24,34, enabling a comparative analysis between the two distinct ways. Using both methods, this work consistently indicates a global reduction in low-level clouds and a decline in high-level clouds over tropical regions in response to deforestation. In addition, we explore the potential physical mechanisms through which deforestation induces alterations in cloud profiles, suggesting that changes in turbulent heat flux could be a crucial factor. Finally, we quantify the impact of deforestation on cloud radiative forcing within the Earth-atmosphere system, with findings indicating that the warming effects of clouds, to a substantial extent, counterbalance the cooling effects of surface albedo at a global scale.

Results

Cloud profile changes

Two distinct approaches (see Methods) are employed in this study to assess the potential impact of deforestation on cloud fraction profiles. The first method draws upon five available GCMs (Table S1) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6)35. It entails analyzing the idealized global deforestation simulations (deforest-glob) conducted in the Land Use Model Intercomparison Project (LUMIP)36, and comparing them against the pre-industrial control simulations (piControl). The second method uses the space-for-time substitution to contrast the multi-year average cloud fraction profiles between the neighboring unaltered forested and unaltered open land grids. In this approach, the potential effects of deforestation on cloud profiles are measured by land cover data from the Moderate Resolution Imaging Spectroradiometer (MODIS) and cloud profiles from the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth reanalysis (ERA5). It should be noted that MODIS only provides data for specific times within the diurnal cycle (morning and noon), which may introduce a low bias on the estimate of forest-cloud impacts in the data23, compared to the model analysis. One significant drawback of the cloud profile data from active satellites is that they have relatively small footprints and sample sizes. As a result, data from numerous satellite passes must be averaged or combined to create a product with sufficient coverage. Given the finer spatial resolution of ERA5 cloud profiles, and their much larger coverage, in comparison to the available gridded data derived from active satellite observations, along with the strong correlation exhibited between ERA5 and the observations (Supplementary Fig. 1), we employ long-term ERA5 data instead. The moderate correlations between ERA5 and satellite-retrieved cloud profiles in the boundary layer are caused by two aspects: one is the data quality of ERA5 itself, and another is the limitation of active satellite sensors on the detection of low-level clouds, especially under conditions of thick upper clouds or strong surface returns37,38,39,40. However, the integration of denser and higher-quality observations over land enhances the accuracy of ERA5 boundary layer cloud data compared to over oceans (Supplementary Fig. 1), thereby better suiting this land-focused study. As GCMs contain both local and no-local effects, we extract the local effects from the total signals (see Methods). Isolating local effects can aid in understanding the biophysical mechanisms of deforestation on clouds. Despite the differing principles behind the two methods, it is noted that the space-for-time substitution also solely considers local effects, allowing for a comparison between these two approaches.

While ref. 28 outlined diverse spatial patterns in how cloud cover responds to deforestation across GCMs in LUMIP, once the local effects are isolated, they reveal consistent spatial patterns (Fig. 1a). This implies that the inconsistencies across models documented by ref. 28 primarily arise from discrepancies in non-local effects. For a quantitative comparison between the GCMs and ERA5, we quantify the sensitivity of the cloud fraction profile to deforestation by calculating the changes in cloud fraction per deforestation fraction. Even with distinct principles, both methods show consistency in this specific change across the spatial distribution regarding cloud vertical profile responses to deforestation (Fig. 1). Globally, cloud cover below 700 hPa decreases in response to deforestation, showing consistency with satellite observations21,22,23. The decrease in tropical cloud cover is restricted to relatively low altitudes according to the ERA5 space-for-time substitution method. The response to deforestation is most pronounced in low-level clouds, with additional reductions found for tropical high-level clouds (higher than 500 hPa). Low-level clouds, which are closely coupled with land surface, form and evolve in response to surface heating, moisture fluxes, and other boundary layer processes41,42,43. Therefore, on a global scale, the impacts of deforestation are expected to be more noticeable on low-level clouds than on high-level clouds. While deep-convection clouds are generally not well-coupled with the surface, surface conditions can influence their initiation44,45. Therefore, deforestation is not limited to affecting shallow clouds that are fully coupled to the surface but can also impact deep convective clouds to a certain extent. As a result, the height of deep convective cloud tops can roughly indicate the maximum altitude at which deforestation affects clouds. Since the cloud top height of deep convective clouds varies across regions, with those in tropical regions reaching higher altitudes than those in boreal zones (Supplementary Fig. 2), deforestation is more likely to affect high-level clouds in the tropics.

Fig. 1: Local sensitivity of cloud fraction profile to deforestation.
figure 1

a Zonal mean of the cloud fraction profile difference between the deforest-glob and piControl simulations (deforest-glob minus piControl) per deforestation fraction. The data were the ensemble mean of the local effect extracted from CMIP6 model simulations (see Methods). The stippling represents four or more of the five models showing the same sign. b Zonal mean ERA5 cloud fraction profile variations per deforestation fraction using the space-for-time substitution (open land minus forest; see Methods). Only latitudes possessing more than ten available samples are considered to ensure representativeness.

Full size image

Discussion of physical mechanisms of forest-cloud impacts

Various biophysical processes are engaged in the interactions between forests and clouds, yet identifying the factors that dictate where cloud enhancement or reduction occurs across global deforested areas has remained unclear21,22,29. In terms of the thermodynamics and moisture factors involved in cloud formation, cloud cover in certain areas might be restricted by the heating needed for uplift46,47. In others, it might be restricted by the availability of moisture48. In the following, we explore these two fundamental factors.

In comparison to forests, open land typically exhibits higher surface albedo (Supplementary Fig. 3) and lower ET (Supplementary Fig. 4). Increased surface albedo from deforestation causes cooling by reflecting more shortwave radiation. This cooling effect is counterbalanced by lower ET8. Both the cooling caused by the surface albedo difference and the warming due to ET difference vary across latitudes, indicating that the magnitude and even the sign of local surface air temperature (SAT) changes resulting from alterations in forests differ across climate regions. When examining SAT changes in deforested areas, shifting from forests to open land induces surface warming in tropical regions (Supplementary Fig. 5). This is primarily due to the prevailing impact of ET on the temperature signal, although alterations in surface albedo partially counteract this surface warming. In contrast, the overall biophysical effect of deforestation leads to surface cooling in the boreal zone from GCMs (Supplementary Fig. 5), which agrees with previous studies from both observations8,49 and simulations24,50. Notably, the impact of surface albedo becomes more pronounced as latitude increases, while the influence of ET tends to diminish with higher latitudes. Hence, in boreal regions, increased surface albedo emerges as the predominant factor of surface cooling. However, surface cooling is not as evident from the space-for-time substitution approach, since SAT is nonlinearly influenced by both radiative and non-radiative processes. Previous observation-based studies also suggest small changes in SAT due to deforestation in boreal regions51,52.

Moreover, the reduction in incoming solar radiation and the drop in SAT caused by the higher surface albedo results in a substantial decrease in sensible heat flux (SH) within the boreal zone; however, in the tropics, the decline in the surface turbulent heat flux primarily stems from the reduction in latent heat flux (LH) due to the dominant role of ET (Fig. 2). Deforestation increases near-surface wind speed due to a decrease in surface roughness (Supplementary Fig. 6), enhancing the heat and water vapor exchange rate between the surface and the air, thereby increasing turbulent fluxes. Apart from being influenced by the near-surface wind speed, SH and LH, respectively, however, are also related to the temperature and humidity gradients between the surface and the air. The fluxes are proportional to the product of the wind speed and the gradient53,54. We find that the decreased temperature gradient between the surface and the air in the boreal zones resulting from deforestation outweighs the role of near-surface wind speed (Supplementary Fig. 7), leading to a reduction in SH. Since there is no proxy for the humidity gradient between the surface and the air, we infer changes in humidity gradient from changes in ET. An increase in wind speed is accompanied by a decrease in ET (Supplementary Fig. 4), suggesting that the decreased humidity gradient between the surface and the air caused by deforestation primarily drives the reduction in LH. Thus, when combining the alterations in LH and SH, the decrease in surface turbulent heat flux depicted in Fig. 2 is evident globally. In conclusion, the response of cloud cover to the reduction in turbulent heat flux is illustrated through the decrease in water vapor supply due to decreased LH in the tropics and the weakening in the uplifting process caused by decreased SH in the boreal regions.

Fig. 2: Changes in surface turbulent heat fluxes due to deforestation.
figure 2

a Global pattern of the surface turbulent heat flux (latent heat (LH) + sensible heat (SH)) difference between the deforest-glob and piControl simulations (deforest-glob minus piControl). Diagonal hashing indicates four or more of the five models showing the same sign of change. b Box plots of the CMIP6 surface turbulent heat flux (LH + SH, LH and SH) differences between the deforest-glob and piControl simulations over both tropical and boreal areas. c ERA5 surface turbulent heat flux (LH + SH) variations due to deforestation using the space-for-time substitution (see Methods). d Box plots of the ERA5 surface turbulent heat flux (LH + SH, LH and SH) variations due to deforestation. The data in (a, b) is the ensemble mean of the local effect extracted from CMIP6 model simulations (see Methods). Boxes in (b, d) show the 25th to 75th percentiles of the data, whiskers display the 5th to 95th percentiles, horizontal yellow lines in the boxes represent the median values, and red dots are the mean values. (e, f) Same as (a) but for LH and SH, respectively. (g, h) Same as (c) but for LH and SH, respectively.

Full size image

The local effects derived from ERA5 contain both the mean-state difference and the secondary mesoscale circulation (See Methods). While the mean-state difference mechanism has been discussed above, the secondary circulation processes can enhance or inhibit the cloud responses to deforestation21,23,55. This secondary circulation-induced cloud change is mainly driven by SH21, indicating that cloud enhancement occurs over open land with higher SH compared to adjacent forests in the tropics, whereas cloud inhibition occurs over open land with lower SH than surrounding forests in boreal areas. Additionally, the magnitude of secondary circulations and their impact on clouds may change diurnally, driven by differential heating contrast in the diurnally varying land surface heat fluxes between adjacent patches with different land covers55.

The type of land cover notably influences cloud formation processes by affecting surface heating, moisture availability, and atmospheric stability56,57,58. Forests generally promote cloud formation due to high moisture levels and low albedo22, while over deserts, typically fewer clouds form, due to low moisture availability and high albedo59. Grasslands have a moderate effect on cloud formation that is intermediate between forests and deserts23. Urban areas, with their unique heat island effect and pollution, potentially influence cloud properties and formation60. The varying impacts of various land covers on cloud formation also explain why the results of the two methods differ. The CMIP6 models only consider deforestation into grassland, whereas the diverse land covers between adjacent ERA5 grids disrupt the signals.

Implications for radiation and climate

Previous studies have concentrated on alterations in surface albedo following deforestation, yet there is a lack of quantitative analysis on changes in cloud albedo subsequent to deforestation22. Clouds, on average, exert a cooling effect on climate61. The decrease in cloud cover with deforestation, therefore, implies a warming effect on climate. The increase in surface albedo resulting from deforestation, in turn, contributes to a cooler climate17. Hence, clarifying the competitive relationship between these two elements is essential to the area of forest biophysical effects.

For a complete analysis, we also examine the disturbance of the outgoing radiation at the top of the atmosphere (TOA). The perturbation of outgoing radiation under all-sky conditions reflects the combined impacts of alterations in both surface and cloud properties from deforestation. Under clear-sky conditions, the radiation perturbations solely arise from alterations in surface properties. Thus, the alterations in TOA outgoing radiation due to cloud cover changes can be obtained through the difference between all-sky and clear-sky conditions (all-sky minus clear-sky, also known as cloud radiative effect). As denoted in Fig. 3, a universal pattern prevails worldwide: alterations in surface properties largely govern the overall outgoing radiation changes, with changes in cloud cover acting as a buffer. On average, from the CMIP6 idealized deforestation experiments, reduced cloud cover offsets approximately 44% of the surface albedo cooling effect; while from the space-for-time substitution method based on ERA5, the relative offset is about 26% (Fig. 3g). The disparities in numerical outcomes primarily result from methodological differences. Nonetheless, both methods lead to consistent conclusions. Given the saturation of CMIP6 latitudinal data, we proceed to examine the zonal disparities (Fig. 3h). The discernible result reveals that the compensatory impact of cloud cover compared to the surface albedo change is stable across latitudes, at roughly 50%.

Fig. 3: Changes in outgoing radiation at the top of atmosphere (TOA) due to deforestation.
figure 3

a, c, e Global pattern of the TOA outgoing radiation (shortwave + longwave) difference between the deforest-glob and piControl simulations (deforest-glob minus piControl), respectively, under all-sky, clear-sky, and all-sky minus clear-sky circumstances. Diagonal hashing indicates four or more of the five models showing the same sign of change. b, d, f ERA5 TOA outgoing radiation (shortwave + longwave) variations due to deforestation using the space-for-time substitution (see Methods). Global mean values and standard errors for (a–f) are shown in (g). The offset ratio is the proportion of all-sky minus clear-sky to the all-sky value. h CMIP6 zonal mean of the TOA outgoing radiation (shortwave + longwave) difference between the deforest-glob and piControl simulations under both clear-sky and all-sky minus clear-sky circumstances. The black line indicates the zonal mean offset ratio and the dashed yellow line is the ratio equal to −0.5. The CMIP6 data were the ensemble mean of the local effect extracted from multi-model simulations (see Methods).

Full size image

When comparing the shortwave and longwave components (Supplementary Figs. 8, 9), however, it becomes evident that the perturbations to the climate come mainly from shortwave, further indicating that changes in surface and cloud albedo are the most main causes. From a global average standpoint, the quantitative competition between clouds and surface albedo becomes apparent, showing that deforestation-induced reduction in cloudy-sky albedo partially counteracts the increased surface albedo by nearly half (Supplementary Fig. 10). Considering that alterations in cloud cover following deforestation approximately counterbalance half of the cooling effect caused by changes in surface albedo, neglecting the shifts in cloud-climate interactions introduces a large bias when investigating the biophysical effects of forests in the future.

Methods

CMIP6 simulations

Cloud fraction profile, tree cover fraction, surface LH, surface SH, surface temperature, surface radiation fluxes, surface ET, near-surface wind speed, near-surface air temperature, as well as radiation fluxes at the TOA from five GCMs (Table S1) participating in the CMIP6 are adopted in this study35. The idealized global deforestation simulations (deforest-glob) from the LUMIP36 are analyzed in comparison to the pre-industrial control simulations (piControl). The deforest-glob setup assumes that a total forest area of 20 million km2 is linearly removed from the top 30% forested area with a fixed rate of 400,000 km2 yr−1 over a period of 50 years across the globe. This is then followed by at least a 30-year simulation with a constant land cover to achieve stable conditions. The last 30 years of the deforest-glob and piControl simulations are compared (deforest-glob minus piControl) to derive the mean response to deforestation28. Due to differences in resolution among GCMs, the ensemble mean statistics are calculated by bilinear remapping of diagnostics from individual GCMs to a 2° × 2° grid, and vertically to 27 pressure levels from 1000 to 100 hPa.

Reanalysis datasets

From the ECMWF ERA562, we utilize the cloud fraction profiles data alongside elevation, surface LH, surface SH, surface temperature, surface radiation fluxes, surface ET, 10-m wind speed, 2-m air temperature, and TOA radiation fluxes to examine the impacts of deforestation. Datasets spanning from 2001 to 2021, featuring a spatial resolution of 0.25° × 0.25° and encompassing 28 vertical pressure levels from 1000 to 100 hPa, are employed for the analysis.

Observed land cover

For delineating forested and open land areas, we use land cover data from the MODIS dataset (MCD12C1, version 6.1)63, relying on the International Geosphere-Biosphere Program (IGBP) classification layer to define the land cover types. Annual data for the years 2001−2021 with a spatial resolution of 0.05°×0.05° are adopted. Here, five forest types (evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, and mixed forest) are merged into a single forest classification. The forest fraction is bilinearly gridded spatially into 0.25° × 0.25° to align with the ERA5 data.

Observed cloud profile and cloud top pressure

In assessing the accuracy of ERA5 cloud profiles, we analyse active satellite-observed cloud profiles. The cloud profile retrievals from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) and CloudSat between 2007 and 2010, are aggregated to a spatial resolution of 2° × 2° and a vertical resolution of 480 m64. The fusion of data from both sensors facilitates an extensive depiction of the vertical cloud structure. This comprehensive view is achieved by leveraging the distinct wavelengths each sensor employs (CloudSat: ~2 mm, CALIPSO: 532 and 1064 nm), catering to various cloud and precipitation particles in both liquid and solid phases.

The International Satellite Cloud Climatology Project (ISCCP)-HGM65 monthly average data with a spatial resolution of 1° × 1° from 2001 to 2016 is employed to obtain the observed cloud top pressure for iced convective clouds. The zonal average ISCCP cloud top pressure is interpolated to GCMs or ERA5 pressure levels.

Climate zones

In this study, climate zones are defined according to the global maps of the Köppen–Geiger climate classification (Version 1)66. The Köppen–Geiger historical map contains 30 climate zones at a resolution of 1 km. Tropical and boreal regions are each merged from corresponding subdivided climate zones.

Extracting local effect from GCMs

Deforestation exerts a local impact on the climate within deforested areas (local effect) by modifying land surface characteristics such as albedo, roughness, and ET. Additionally, it affects both deforested and open land grids by altering the advection of heat and moisture, as well as influencing atmospheric circulation (non-local effect)67. Distinguishing between local and non-local effects within GCMs is crucial as coupled models encompass the complete climate response to deforestation, incorporating both local and non-local impacts. Moreover, it allows to develop a more profound insight into the mechanisms influencing the local effects in comparison to those governing the non-local effects. Mesoscale processes typically have spatial scales between 10 and 1000 km. In the CMIP6 models, mesoscale circulations within the analysis resolution (10–200 km) are partly local, while the larger ones (200–1000 km) are considered non-local. Therefore, we assume that the term “non-local effect” in the GCMs refers to both the large-scale circulations (>1000 km) and the mesoscale circulations beyond the analysis resolution (200–1000 km).

Here, we use a chessboard method as outlined by ref. 34 to assess the local effect. This method assumes that the unaltered and adjacent deforested grids share the same non-local effect21,67. The unaltered grids indicate that the forest cover within these grids remains unchanged in the deforest-glob simulation. Since the cloud cover changes within unaltered forest cover grids are entirely due to non-local effects, to generate a global map of the non-local effect, we determine the non-local effects within the deforestation grid by interpolating the signals from the surrounding unaltered forest cover grids. As spatial interpolation might alter existing values, we maintain the non-local signals within the unchanged forest cover grids as they are and only derive the non-local signals for the deforested grids. The local effect over the deforested grids thus can be derived by subtracting the interpolated non-local effect from the total effect. Notably, employing a chessboard-like method introduces horizontal interpolation errors, given that the local effect relies solely on interpolation from neighboring, unaltered grids. However, our study is centered on idealized deforestation scenarios, and prior study24, has demonstrated the possibility of isolating local effects using similar methodologies and datasets. Winckler, et al.34 conducted comparisons between simulations involving both sparse and extensive idealized deforestation, finding small differences in derived local effects from spatial interpolation. The difference between the two local effects of sparse and extensive deforestation simulations is of secondary importance as compared to the local effect itself. Additionally, by comparing the local and non-local effects of SAT as an example (Supplementary Fig. 11), we find that both effects are non-negligible relative to the total effect. Therefore, the calculated local effect is unlikely to introduce large uncertainties due to discrepancies in magnitudes with the non-local effect.

To mitigate uncertainties, we use ensemble mean results from five available GCMs. Grid points where four or more of the five models exhibited the same sign are highlighted to demonstrate where there is a high consistency among the models. However, it should be noted that since the model results are derived from idealized deforestation experiments, they may appear overly simplistic compared to observations.

Space-for-time substitution

In addition to idealized deforestation simulations, this study employs a space-for-time substitution method to assess the impacts of deforestation, combining MODIS land cover and ERA5 reanalysis datasets. Such an approach has previously been applied in various studies to evaluate the effect of alterations in land cover on temperature8,26, the surface energy budget5,68, or cloud cover21,22. The fundamental premise of this method is that neighboring land patches share the same climatic background, and variations in their characteristics can act as a proxy for temporal changes. This method exclusively includes the local effects, comprising two components: one is the mean-state difference—variations in land cover conditions that result in spatial disparities; another is the secondary mesoscale circulation within the moving window pixels —differential heating of adjacent land cover patches can create sea-breeze-like secondary mesoscale circulations21,23,55. However, it should be acknowledged that the space-for-time substitution approach also carries uncertainties as it is an indirect method to calculate local biophysical effects.

Areas designated as unaltered forested (or unaltered open land) are identified as pixels where the initial (in 2001) tree cover fraction exceeds 60% (or is below 40%) and with a net change in forest cover <10% from 2001 to 2021. Pixels with water coverage >10% are excluded. We use a moving window approach to search for comparison samples between unaltered forested and unaltered open land pixels. We choose for each moving window a size of 7 × 7 pixels (1.75° × 1.75°). To reduce the influence of topography69,70, we calculate the standard deviation (s.d.) of elevation within specific moving windows and omit samples where this s.d. exceeds 100 m following ref. 21. Finally, the potential effect of deforestation on a specific variable (\(\Delta {{{\rm{Var}}}}\)) is quantified as:

$$\Delta {{{\rm{Var}}}}={{{{\rm{Var}}}}}_{{{{\rm{open}}}}\, {{{\rm{land}}}}}-\overline{{{{{\rm{Var}}}}}_{{{{\rm{surrounding}}}}\, {{{\rm{forests}}}}}}$$
(1)

or

$$\Delta {{{\rm{Var}}}}=\overline{{{{{\rm{Var}}}}}_{{{{\rm{surrounding}}}} \, {{{\rm{open}}}} \, {{{\rm{lands}}}}}}-{{{{\rm{Var}}}}}_{{{{\rm{forest}}}}}$$
(2)

where Eqs. (1) and (2) are applicable to the case where the central pixel of the moving window is unaltered open land and unaltered forest, respectively. \({{{{\rm{Var}}}}}_{{{{\rm{open}}}} \, {{{\rm{land}}}}}\) and \({{{{\rm{Var}}}}}_{{{{\rm{forest}}}}}\) are multi-year mean variables over unaltered open land and unaltered forest pixels, respectively. \(\overline{{{{{\rm{Var}}}}}_{{{{\rm{surrounding}}}} \, {{{\rm{forests}}}}}}\) and \(\overline{{{{{\rm{Var}}}}}_{{{{\rm{surrounding}}}} \, {{{\rm{open}}}} \, {{{\rm{lands}}}}}}\) are the average values of the surrounding \({{{{\rm{Var}}}}}_{{{{\rm{forest}}}}}\) and \({{{{\rm{Var}}}}}_{{{{\rm{open}}}} \, {{{\rm{land}}}}}\) within a moving window when the central pixel is unaltered open land and unaltered forest, respectively.