Contrasting impacts of forests on cloud cover based on satellite observations

Cloud cover and environmental datasets

The monthly mean MODIS cloud fraction at 0.05° used in this study was computed from the daily cloud mask data (“cloudy” label for the bits 0–1 of “state_1 km” band) included in the MODIS Surface Reflectance product (MYD09GA.006, overpass at local time of 13:30) of Aqua from 2002 to 2018, using the reduceResolution function with “mean” aggregation method on Google Earth Engine (https://earthengine.google.com/). The 1-km cloud mask was produced based on the MOD35_L2 cloud mask product, which had been extensively validated^71,72. Before computing cloud fractions, a snow/ice flag (the bit 12 of “state_1km” band) was used to remove snow or ice pixels in the cloud record because the high reflectivity of snow/ice degrades the accuracy of cloud detection, especially during winter in the northern hemisphere. Therefore, the estimated cloud effect would have larger uncertainty in boreal winter than in summer.

To complement MODIS-based cloud analyses, we used the Meteosat Second Generation (MSG) hourly cloud fraction data of 2004–2013 at a spatial resolution of 0.05°. The Coordinated Universal Time (UTC) of the raw MSG hourly cloud cover data was converted to local time before being used for analysis.

The cloud fraction from Sentinel-5P Near Real-Time (NRTI) data product was used in this analysis. This dataset is available from 2018-07-05 at a spatial resolution of 0.01° and it has an overpass time of 13:30 similar to MODIS. The Sentinel-5P cloud data, although having a short period of 2 years, allows for the separation of cloud effects into different cloud types, with the help of a cloud classification scheme based on cloud top pressure and cloud optical depth information³⁰.

Environmental variables include evapotranspiration (ET, MOD16A2 V6), land surface temperature (LST, MYD11A1 V6) from MODIS, and soil moisture (SM) from the TerraClimate dataset. All these environmental variables were averaged into monthly means at 0.05° resolution.

Elevation data are from SRTM Digital Elevation Data at 0.05° resolution. Land cover data include MODIS (MOD12C1) and European Space Agency (ESA) global land cover products, which were aggregated to 0.05°.

Defining forest cover change

To define forest/non-forest and forest cover change, we used the Global forest cover (GFC) product which provides global tree cover for the year 2000 (baseline), yearly forest loss from 2001 to 2018, and forest gain from 2000–2012 at 30 m resolution⁵³. The GFC data were aggregated to fractions at 0.05°. Net forest cover change was calculated as the sum of the loss and gain accumulated throughout the study period. Pixels with net forest cover change fractions smaller than 0.05 are considered to be “unchanged” and greater than 0.05 are considered to be “changed”. Unchanged forests and unchanged non-forest were defined as pixels with baseline tree cover fraction greater or less than 0.5 and with net forest change <0.05. For unchanged non-forest, pixels classified as water, snow/ice, or wetland were excluded using the major composite of MODIS land cover from 2002 to 2005 with the International Geosphere-Biosphere Program (IGBP) classification scheme. For “changed” forest pixels, forest loss was identified as those with a net forest loss >0.15. Forest loss defined this way is expected to pose a stronger signal on clouds than that with a lower threshold, and thus improves the detectability of cloud impact against natural variability of cloud cover.

Estimating potential and actual impacts of forest loss on cloud cover

The potential effect of forest on cloud (ΔCloud) was quantified as the mean cloud difference between unchanged forests and nearby non-forest as:

where Cloud_forest and Cloud_nonforest are multiyear or yearly mean cloud fractions averaged over unchanged forest and unchanged non-forest pixels, respectively. ΔCloud defined this way, with the reversed sign, represents the potential impact of forest loss on cloud cover at a given location. The methodology is designed to isolate the cloud effects of land surface conditions from those caused by meteorological conditions. It refers to local cloud impact (caused by land surface conditions) because effects from synoptic conditions and large-scale circulation changes/climate changes (meteorological conditions) are shared by both forest and non-forest and are therefore minimized through subtraction. If there is no effect of forests on cloud cover, the resulting ΔCloud would show random patterns with mixed positive and negative values instead of a systematic pattern, which indicates a cloud preference over forests or non-forest.

To implement Eq. 1, we used a moving window approach to search for comparison samples between forest and nearby non-forest pixels at locations that underwent “forest change” (i.e., net forest change >0.05) across the globe⁷³. Each moving window was sized at 9 × 9 pixels (0.45° × 0.45°) and two adjacent windows were half-overlapped with a distance of 5 pixels (i.e., the centers of two windows were 5 pixels apart along latitudinal and longitudinal direction). To avoid cloud inhibition effects from water bodies⁷⁴, water pixels and their one-pixel buffer zone were masked out in the window searching strategy for ΔCloud. Therefore, ΔCloud can be calculated using unchanged forest and non-forest pixels within each moving window. This window searching strategy ensures the proximity of the forest and non-forest pixels to pixels that underwent forest change, making the estimated potential effect more representative of the actual forest change impact. To test the sensitivity of ΔCloud to window size and time period, ΔCloud was also estimated using alternative window sizes: 11 × 11 (0.55° × 0.55°), 21 × 21 (1.05° × 1.05°), 51 × 51 (2.55° × 2.55°) pixels and different periods (2002–2007, 2008–2013, and 2014–2018). The resulting ΔCloud was similar to results with the window size of 9 × 9 (0.45° × 0.45°) and among split time periods (Supplementary Figs. 2, 3). Unlike using direct comparison in cloud cover (and other biophysical variables) between forest and non-forest, an alternative method is to utilize the regression coefficients of cloud cover (dependent variable) to land cover fraction (independent variable) and estimate cloud effects assuming 100% land conversion, as adopted by ref. ⁵⁸. The alternative regression-based approach is mathematically more complicated, and its implementation involves non-trivial post-processing compared with our method while producing qualitatively similar results.

A similar window searching strategy was applied to estimate the differences between forests and non-forest in LST (ΔLST), ET (ΔET), and soil moisture (ΔSM) (Supplementary Fig. 10).

The cloud impact estimated as the cloud differences between forest and non-forest could be confounded by their differences in topography, which is known to be an important factor for cloud formation. To minimize the topographic influence, we calculated the standard deviation (s.d.) of elevation within each moving window and removed samples with s.d. >100 m from the analysis. This filtering effectively excluded comparison samples over complex terrain such as mountainous regions so that the retained samples came from relatively flat areas.

The actual effect of forest loss on cloud (ΔCloud_loss) was quantified as the cloud cover difference between forest loss (Cloud_loss) and nearby unchanged forest pixels (Cloud_forest) using the same window searching strategy as the potential effect (Eq. 2).

where ΔCloud_loss is the actual impact of forest loss on cloud cover, Cloud_loss and Cloud_forest are the multiyear or yearly mean cloud cover averaged over forest loss and unchanged forest pixels, respectively. The actual impact (deforested vs. forests) shows good spatial resemblance to the potential effect (non-forest vs. forests, ΔCloud with the reversed sign), suggesting that the potential effect can provide a priori prediction of possible cloud change induced by forest loss (the correlation of the spatial pattern is 0.44, p < 0.05).

To quantify the progressive tree cover changes caused by forest loss, we calculated tree cover differences between forest loss and unchanged forest pixels following Eq. 3,

where ΔTree_year is the tree cover difference between forest loss and unchanged forest pixels at a given year. It is the sum of the tree cover difference in the baseline year 2000 (Tree2000_loss − Tree2000_forest) and the accumulated yearly forest loss differences from 2001 until a given year (the sigma term of Eq. 3).

The comparison samples obtained from the window searching strategy for potential and actual impacts were aggregated to 0.5° for display and further analysis.

Cloud effects of forests separated into different cloud types

By using cloud top pressure and cloud optical depth from the daily Sentinel-5P NRTI data, nine cloud types were classified according to the ISCCP (International Satellite Cloud Climatology Project) cloud classification scheme³⁰. The classified cloud types were 1-cirrus, 2-cirrostratus, 3-deep convection, 4-altocumulus, 5-altostratus, 6-nimbostratus, 7-cumulus, 8-stratocumulus, and 9-stratus. Cloud types 1–3, 4–6, and 7–9 corresponded to low, mid-, and high-clouds, respectively. Cloud types 3, 7, and 8 were convective clouds and the latter two were shallow convective clouds. The multiyear mean JJA total cloud fraction and fraction of each cloud type were calculated during the available time period and were aggregated to 0.05° from the original 0.01° resolution. We then applied the same moving window method to estimate the cloud effects of forests for total cloud cover as well as for different cloud types. The summed cloud effects of each cloud type equaled the total cloud cover effects. We expected convective cloud types (types 3, 7, and 8) to be influenced by forests, while other non-convective cloud types would not, so that their ΔCloud would show a more random pattern. The dominant cloud type for cloud effects of forests was determined by the cloud type whose ΔCloud had the same sign as the total cloud effect and had the largest magnitude (Supplementary Fig. 7).

We noted that there were regional differences in the cloud effects estimated from Sentinel-5P and the magnitude of the effect was also smaller than the other two datasets. For example, the Southeast US in MODIS was dominated by negative ΔCloud (64.67%) whereas in Sentinel-5P it showed more positive ΔCloud (57.09%) (Supplementary Fig. 13). The large spatial coverage of positive ΔCloud in Europe in MODIS and MSG was slightly reduced with Sentinel-5P. These regional differences might be linked to potential bias in cloud fractions of Sentinel-5P, because we found that cloud fractions of Sentinel-5P were systematically lower than that of both MODIS and MSG. However, the cloud effects of Sentinel-5P in the Amazon were consistent with MODIS (80.95%) in terms of coverage, showing a prevailing cloud inhibition (81.45%) (Supplementary Fig. 13). Cloud inhibition in Central Africa with a spatial coverage of 51.84% was slightly more in line with the widespread negative ΔCloud in MSG (67.56%) than in MODIS (36.40%).

Given these differences in the cloud effects among datasets, the results from Sentinel-5P still provided strong support that convective clouds dominated the cloud effects of forests at both global and regional scales (Supplementary Fig. 7).

Attribution of cloud effect of forests

Since cloud effects of forests may result from contributions of both vegetation properties and orography, we used tree cover and elevation as indicators to represent each of their effects. Elevation was selected as an indicator of the orographic lifting mechanism. We acknowledge that the reality is much more complicated than this highly simplified representation of the orographic cloud effect. However, for a global-scale analysis, elevation could still provide a first-order approximation of the orographic effect.

To decompose the potential cloud effect of forests into contributions from tree cover and elevation, we first estimated sensitivities of cloud cover to tree cover and elevation respectively, following a linear regression model defined in Eq. 4.

where S_tree and S_ele were the sensitivities of cloud cover to tree cover and elevation respectively, and the intercept c was unused in this study. The sensitivity parameters were estimated for each moving window separately if it had a nonzero tree cover. The estimated slope of cloud cover to elevation (S_ele) was positive in the majority of the world (Supplementary Fig. 6d), suggesting that a higher elevation indeed promotes cloud formation. Next, we calculated tree cover differences (ΔTree) and elevation differences (ΔEle) between unchanged forests and non-forest pixels similarly to ΔCloud. Then the cloud differences induced by tree cover (ΔCloud_tree) and by elevation (ΔCloud_ele) can be obtained by multiplying their sensitivities by the corresponding differences as Eqs. 5 and 6. The sensitivity and differences parameters were averaged to 0.5° resolution before using in Eqs. 5 and 6 (Supplementary Fig. 6).

The reconstructed ΔCloud given by the sum of ΔCloud_tree and ΔCloud_ele explained about 70% of the original ΔCloud.

To attribute ΔCloud to tree cover and elevation-induced cloud changes, we compared the sign and magnitude of original ΔCloud, ΔCloud_tree, and ΔCloud_ele. If ΔCloud_tree and ΔCloud_ele both had the same sign as ΔCloud, the one with greater magnitude was classified as the dominant factor. If only one of ΔCloud_tree and ΔCloud_ele had the same sign as ΔCloud, the factor with the same sign was classified as the dominant factor. If neither ΔCloud_tree nor ΔCloud_ele had the same sign as ΔCloud, the dominant factor was classified as others. As a result, the potential cloud effects could be attributed to five classes: tree cover induced cloud increase (Tree+) and decrease (Tree−), orography induced cloud increase (Orography+) and decrease (Orography−), and others.

Linking cloud effect with sensible heat flux

Sensible heat data were obtained from three independent sources: satellite estimate⁴, a Community Land Model version 5 simulation⁷⁵, and 30 paired forest and non-forest flux sites³⁵.

Satellite estimates provide changes in the combined sensible heat and ground heat fluxes (H+G) under different land cover conversions at 1° spatial resolution based on MODIS data (a total of 45 pairs of land conversions for “HG_IGBPdet”). The combined fluxes of H+G were estimated as the residual of surface energy components as described in ref. ⁴. Due to the small contribution of G to H+G, we referred to “H+G” as “H” for simplicity in the following text and the main text. To obtain sensible heat differences between forest and non-forest (ΔH) that are compatible with ΔCloud, we extracted the dominant land cover type for unchanged forest (e.g., evergreen broadleaf) and non-forest pixels (e.g., crop) within each moving window from the ESA land cover product. The dominant land cover types for forests and non-forest were upscaled to 1° resolution with the “major” method (figure not shown for 1°, but a similar one for 0.5° is shown in Supplementary Fig. 9). For each one-degree grid box with a dominant forest type (e.g., evergreen broadleaf) and non-forest type (e.g., crop), ΔH can be extracted from the corresponding sensible heat change value that matches the specific land conversion (e.g., evergreen broadleaf to crop) at the same grid box from the 45 pairs of land cover conversions defined within the “HG_IGBPdet” dataset.

CLM5 is the land component of a state-of-the-art earth system model Community Earth System Model 2 (Ref. ³⁴). The CLM5 simulation was conducted at the spatial resolution of 0.5° from 1997 to 2010, driven by a revised climatology GSWP3 as the atmospheric forcing (http://hydro.iis.utokyo.ac.jp/GSWP3/), with the plant phenology prescribed from satellite products, the land cover of 2000, and the separated soil columns configuration^76,77. The years 1997 to 2001 were the spinup period and excluded from the analysis (please see detailed description in ref. ⁷⁵). In CLM, different types of vegetation within a grid cell are represented as separated tiles of different plant functional types (PFTs). We used subgrid PFT-level model outputs to calculate sensible heat differences between different land cover types within the same model grid. The subgrid tiles within a model grid cell share the same atmospheric forcing, therefore replicating the assumption of similar meteorological conditions of the space-for-time approach¹². To match the CLM5 model resolution, the dominant land cover types for forests and non-forest of each moving window were upscaled to 0.5° using the ESA land cover data (Supplementary Fig. 9). Because CLM adopted a different land classification scheme, we created a look-up table to convert CLM land cover to the IGBP classification scheme (Supplementary Table 3). The differences in the sensible heat flux (ΔH) between a specific forest and a non-forest type can be extracted from the sensible heat values of the corresponding PFTs.

A total of 30 paired flux sites were used in this study to calculate sensible heat differences between forest and non-forest (ΔH). Twenty-eight site pairs were processed by ref. ³⁵ using FLUXNET data and two additional Amazon site pairs were from the ORNL archive⁷⁸ (Supplementary Table 1). ΔH was calculated as the mean sensible heat flux difference between the paired forest and non-forest site during the daytime (8:00 to 16:00). ΔCloud for each site pair was extracted from the central location of the line linking two sites. Unlike ΔCloud used in the main analysis which was aggregated to 0.5°, we here used ΔCloud aggregated to 1° without the elevation s.d. criteria and the one-pixel water buffer removal to increase available ΔCloud value for each site pair. When analyzing the relationship between ΔH and ΔCloud, two flux pairs were excluded because the matched ΔCloud was missing (pair 29) and an outlier in ΔH (pair 22 with ΔH > 200 W/m²).

Scale-dependency of potential cloud effect of forest

To investigate how the potential cloud effect varies with spatial scale, we reprocessed the MODIS cloud cover and GFC data into different spatial resolutions to emulate the scale change (using “mean” for cloud cover and “major” method for forest cover). Specifically, the 0.05° cloud and GFC data used in the main analysis were aggregated to coarser resolutions (0.1°, 0.25°, 0.5°, and 1°) and ΔCloud was re-estimated with the window searching strategy of slightly different configurations to accommodate the resolution change (Supplementary Fig. 12). The specific parameters of the window searching strategy under different resolutions are provided in Supplementary Table 2, including raw data resolution, window size, window distance, and display resolution. For a given resolution, ΔCloud was estimated with two-parameter combinations to ensure the robustness of the results.

Source: Ecology - nature.com