Contrasting impacts of forests on cloud cover based on satellite observations
Cloud cover and environmental datasetsThe 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 validated71,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 information30.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 changeTo 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 resolution53. 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.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 coverThe potential effect of forest on cloud (ΔCloud) was quantified as the mean cloud difference between unchanged forests and nearby non-forest as:$$Delta {{{{{rm{Cloud}}}}}}={{{{{{rm{Cloud}}}}}}}_{{{{{{rm{forest}}}}}}}-{{{{{{rm{Cloud}}}}}}}_{{{{{{rm{nonforest}}}}}}}$$
(1)
where Cloudforest and Cloudnonforest 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 globe73. 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 bodies74, 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. 58. 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 (ΔCloudloss) was quantified as the cloud cover difference between forest loss (Cloudloss) and nearby unchanged forest pixels (Cloudforest) using the same window searching strategy as the potential effect (Eq. 2).$$Delta {{{{{{rm{Cloud}}}}}}}_{{{{{{rm{loss}}}}}}}={{{{{{rm{Cloud}}}}}}}_{{{{{{rm{loss}}}}}}}-{{{{{{rm{Cloud}}}}}}}_{{{{{{rm{forest}}}}}}}$$
(2)
where ΔCloudloss is the actual impact of forest loss on cloud cover, Cloudloss and Cloudforest 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 200 W/m2).Scale-dependency of potential cloud effect of forestTo 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. More
