Revealing the widespread potential of forests to increase low level cloud cover

Assumptions for the space-for-time substitution

The main methodological concept in this study is the notion of a space-for-time substitution. Such approach has previously been used in various studies to estimate the effect of land cover change on temperature^20,22 or on the surface energy balance^22,72. The overarching assumption behind the method is that the difference in properties of neighbouring patches of land can serve as a surrogate for changes in time. While this main assumption largely holds for land surface properties, such as skin temperature, it requires a more detailed articulation into several underlying assumptions in order to apply the approach to atmospheric properties such as cloud cover. This is because atmospheric properties are prone to lateral movements, partially decoupling them from the land cover directly below them, and thus adding considerable complexity to the analysis.

The first underlying assumption is that the method will be mostly sensitive to low-level convective clouds generated in the boundary layer (i.e. cumulus clouds). These are typically formed under stable conditions of high pressure and low wind, and are thus expected to have a higher spatial correlation with the underlying landscape elements. Other types of low-level clouds, such as stratus clouds, are typically much more uniformly spread across the landscape, which would result in no difference in CFrC when comparing two distinct and neighbouring vegetation classes. For medium- or high-level clouds, their position will be determined mostly by the state of the atmosphere rather than by the land surface, resulting in a very low correlation with vegetation spatial patterns. The space-for-time substitution approach would thus similarly result in white noise.

The second assumption is that the boundary layer cumulus clouds will see very limited lateral advection between the moment of their formation and the satellite observation. Cumulus clouds over land show a very stable climatology where the peak formation is largely confined to the early afternoon (around 14:00), timing which remains very stable across space and season⁴³. Therefore this assumption should largely hold if the observations are made at this time.

The third assumption is that if we consider topographically flat terrain that is away from a coastline, general weather conditions are essentially the same at a local scale (i.e. a region of radius circa 25 km around a given point). Within such an area, we then assume that variations in low cloud cover are mostly determined by local differences in surface properties, themselves determined by the type and condition of the present land cover.

Preparation of input datasets

This study requires gridded geospatial datasets for two variables: cloud fractional cover and land fractional cover. Both datasets used here have been prepared in the frame of the European Space Agency’s (ESA) Climate Change Initiative (CCI)⁷³.

The Cloud CCI⁵⁵ provides a series of cloud properties derived from distinct satellite Earth observation platforms in a harmonized way. Here we use their cloud fractional cover variable (henceforth CFrC), which describes the fraction of a 0.05° × 0.05° pixel covered by clouds based on observations made at a finer spatial resolution at the given time of the satellite overpass. We chose to use Cloud CCI dataset based on the MODIS instrument on-board of the Aqua platform for two reasons. First, the timing of overpass of the Aqua platform (circa 13:30 local time at the Equator) coincides very well with the timing of peak of cumulus cloud formation⁴³, thus greatly limiting the extent of possible cloud advection between the moment of cloud formation and observation. Second, native spatial resolution of the MODIS instrument is superior to the alternative (AVHRR), and should result in a better sensitivity to the presence of small cumulus clouds. More specifically, out of the 5 spectral bands of the MODIS instrument used by the Cloud CCI to characterize cloud properties (bands 1, 2, 20, 31 and 32), two of them (bands 1 and 2) have a native spatial resolution of 250 m. While these are aggregated to 1 km (the spatial resolution of the other MODIS bands) prior to their ingestion in the cloud retrieval algorithm, their finer native granularity and quality should prove to be an asset for small cumulus cloud detection. The CCI MODIS-AQUA CFrC data is available for the period 2004–2014. The values are first averaged from daily to monthly scale, and then a single monthly value is calculated for every pixel over the period 2004–2014. The results are 12 layers each representing the multi-annual average CFrC for a given month.

The second type of data needed for the analysis is the fraction of the 0.05° × 0.05° pixels that are covered by distinct vegetation types (essentially trees and grasses) and by other land cover classes (urban areas, bare soil, etc.). These are derived from the Land Cover CCI⁵⁴, a set of consistent annual maps describing, with a spatial resolution of 300 m, how the terrestrial surface is covered based on The United Nations Land Cover Classification Scheme⁷⁴. This information is aggregated both spatially and thematically using a specifically designed framework⁷⁵ to produce maps of general land fractional cover with a spatial resolution of 0.05° to match that of the cloud fractional cover data. The procedure is very similar to that done in a previous study²². For the context of this study, which has a focus on afforestation, the interest lies on transitions among three main vegetated classes, namely: deciduous forest, evergreen forests and herbaceous vegetation. Herbaceous vegetation is composed of both grasses and crops, irrespective of management practice such as irrigation. While irrigation has a clear biophysical effect of its own⁶⁰, we deemed the land cover product was not consistent enough for this specific class. For reasons that are explained in the respective methodological section below, the full compositional description of the landscape is necessary (i.e. beyond the classes of interest), and therefore land cover fractions of the following classes are also generated: shrublands, savannas, wetlands, water, bare or sparsely vegetated, snow or ice, and urban.

Retrieving potential cloud fractional cover change

Under the above-mentioned assumptions, we apply a space-for-time substitution algorithm developed in a previous study²² to the cloud fractional cover and land fractional cover datasets. We summarize the main aspects of the methodology, along with the few necessary adaptations, but the reader requiring more detail is redirected to the original papers^22,76. The approach consists in applying an un-mixing operation over a spatially moving window containing n pixels. Over each window we apply a linear regression based on a matrix X containing the explanatory variables, in which each column of X represents the fractional cover of a given land cover type for each of the n pixels. The response variable is a vector y containing the n values of CFrC for the n pixels, while the vector β represents the regression coefficients:

This is equivalent to solving the following system of equations:

in which the digits of the subscript of x, e.g. x_ij, represent the land cover fraction j in pixel i, for the n pixels in the moving window and the m classes that are considered. Once identified, we can use the β coefficients to predict the local y value corresponding to a given composition x, including that composed of a single land cover j by setting x_j = 1 and all other x values to zero. However, applying a regression directly on X carries a risk due to the compositional nature of the data (i.e. the sum of each row adds up to one), as the analysis of any given subset of compositional components can lead to very different patterns, results and conclusions⁷⁷. To avoid this, we reduce the dimensionality of X through singular value decomposition (SVD) after removing the mean of each column:

where M is the appropriate matrix of column means, U and V are the matrices containing, respectively, the left-hand and right-hand singular vectors, and D is a diagonal matrix containing the singular values representing the standard deviations of the ensuing dimensions. The squared values of D represent the variance explained by each dimension, and can thus serve to define z, a reduced subset of dimensions that conserves 100% of the original matrix’s variation. The corresponding z right-hand singular vectors, V_z, can then be used to find the appropriately transformed predictor matrix of reduced dimension Z as follows:

which can now be regressed onto the CFrC y:

where Z has been augmented with a leading column of 1s to accommodate an intercept term in the regression. We then use the standard method to obtain an estimate of β_z:

However, because of the matrix transformation from X to Z, the regression coefficients β_z do not provide direct information on the relationship between land fractional cover and cloud fractional cover (as in a normal regression). To identify the z values associated with a particular vegetation or land cover type (within the local analysis defined by the moving window), we define a ‘dummy pixel’ whose composition contains only a single class, with all other classes in its composition set to zero. This pixel’s composition is then transformed, and its y value predicted. This is the y associated with that vegetation type. To generalize this for all compositional components of interest, we define a matrix P with as many rows as these compositional components that we wish to predict. P is centred on the same column means as above (M, specific to each local analysis), and then multiplied by the correct number of transposed right-hand singular vectors (V_z, again, specific to each local analysis).

Predicted y_p values for each vegetation or land cover type (identified by predicting the appropriately transformed ‘dummy pixels’) are then calculated as:

The expected change in variable y associated with a transition from vegetation type A (e.g. herbaceous vegetation) to vegetation type B (e.g. deciduous forest) at the centre of the local window is then the difference between the y_p predicted for each ‘pure’ vegetation type:

The uncertainty in the estimation of Δy_A→B can be expressed as a standard deviation using the following expression:

where ({sigma }_{rm{A}}^{2}) and ({sigma }_{rm{B}}^{2}) are the variances in the estimates of y_A and y_B, and σ_AB is their covariance. These variances and covariances are in turn obtained from the covariance matrix, defined from the regression as:

The diagonal terms in Σ are the variances of individual predictions of (individual) classes. The off-diagonal parts of Σ hold the covariances between these predictions. As a reminder, the uncertainty σ_A→B calculated in this way is related to the methodological uncertainty and does not include the uncertainty in the input variables of land cover or cloud fractional cover.

In the default set-up for this study, we concentrate on two transitions: herbaceous vegetation to deciduous forest and herbaceous vegetation to evergreen forest. These are calculated using a spatial window of 7 × 7 pixels, each pixel being of 0.05°, resulting in a squared spatial window of circa 35 km in size. To ensure there are enough values to do the un-mixing over each window, we established that there must be a minimum of 60% of valid values in each window, and that at least 40% must have distinct compositions. The operation is applied to all 12 monthly layers of CFrC, resulting in 12 maps of Δy with a 0.05° spatial resolution for each of the two vegetation cover transitions.

Post-processing

A series of post-processing steps are required to ensure the results of the Δy maps can be used to evaluate the effect of land on cloud cover. The first step is to mask all pixels in which there is insufficient co-occurrence of the two vegetation classes involved in the transition. This co-occurrence is quantified by an index of vegetation co-occurrence⁷⁶, I_c, calculated from the land fractional cover layers using the same spatial moving window of 7 × 7 pixels as used before. This index is calculated pairwise, i.e. for 2 vegetation classes of interest A and B, using two vectors p_A and p_B, describing the presence of these two vegetation classes in each of the i pixels in the moving window. It also requires the definition of another i point evenly distributed along a hypothetical line B = 1 − A in the two-dimensional space describing the presences of vegetation class A and vegetation class B. These points, whose position in the 2-D space are labelled q_A and q_B, represent an ideal situation of maximum co-occurrence that serves as a reference to establish the index. The formal definition of the index is thus:

The minimum operator in the numerator selects the smallest distance that a given point p can have to any of the q points. The sum relates to the sum of this distance for all i points in the spatial moving window. The denominator characterizes the maximum distance that the point p can encounter. I_c will range from 0 to 1 corresponding to a gradient of ‘no presence of either class’ to ‘full and evenly balanced presence of both classes’. As in⁷⁶, we retain only pixels with I_c ≥ 0.5 where we consider that there is sufficient information at local scale concerning both vegetation types to derive meaningful information about the target land cover transition.

The second step is to remove the potential orographical effects, which can be especially problematic given that forests are more likely to be located over mountainous areas due to human action⁵⁶. Here we mask the areas where considerable topographical variation occurs within the 7 × 7 pixel moving window of interest using the same implementation described in⁷⁶. This involves using 3 different indicators, v₁, v₂ and v₃, calculated over the moving window based on μ_h and σ_h, which are, respectively, the mean and the standard deviation of elevation over each grid cell of the input cloud dataset. These are defined as follows:

For an interpretation of these metrics, readers are invited to read⁷⁶. These three indicators are combined together in a single layer depicting all pixels satisfying all of the following conditions: v₁ < 50 m, v₂ < 100 m and v₃ < 100 m. Pixels that fail any of these conditions are masked out from all the layers of results.

The third step is to aggregate the information to a coarser spatial support corresponding to the moving window. This step combines a decorrelation operation with a weighted averaging to result in a single value for a new pixel with a spatial resolution of 0.35°. The decorrelation is necessary as every spatial moving window overlaps, resulting in highly auto-correlated results. The degree of auto-correlation is characterized by an n × n matrix R_a containing the fraction of overlap between every pair of windows. By combining R_a with D_a, a diagonal matrix containing the estimation uncertainties for the central pixel of each window in its diagonal, we can build a new covariance matrix Σ_a (the subscript ‘a’ is used to differentiate these matrices involved in this aggregation step from those used before):

The vector of weights for the spatial aggregation from 0.05° to 0.35° that includes the decorrelation is then obtained as:

which can then be used to calculate the aggregated (overline{{{Delta }}y}) as:

while the aggregated uncertainty ({sigma }_{overline{{{Delta }}y}}^{2}) is given by:

Again, for more details on these operations, see^22,76.

The methodological procedure is still susceptible to the generation of some unrealistic data. A reason for this might be uncertainties and errors in the input data, which could not be explicitly taken into account in the methodology. To mitigate their effect, we decided to filter for outliers. We first removed the aggregated (overline{{{Delta }}y}) values that are calculated based on less than 20% of the underlying 0.05° pixel values. We also removed values in which the aggregated methodological uncertainty (({sigma }_{overline{{{Delta }}y}}^{2})) is above 0.1. We then removed data falling outside the 1st and 99th percentiles of the entire dataset for each transition.

The final post-processing step is to merge together the two separate datasets corresponding to the two different forest types (deciduous and evergreen). This is done using weighted averages based on their respective presence in each grid cell. A simplified version is also produced at a reduced resolution of 1°, useful for visualization purposes (e.g. in Fig. 1).

Considerations on the space-for-time method

The indirect biophysical effects of land cover change discussed in this paper refer exclusively to local cloud cover (so-called ‘local effect’). Our methodology relies on spatial gradients in a local moving window, and it is therefore intrinsically incapable of quantifying large-scale non-local effects. However, it has been shown with modelling studies that changes in forest cover can also have relevant non-local biophysical effects. For instance, albedo-induced cooling following deforestation is mainly a non-local effect through mechanisms of teleconnections that may ultimately compensate the opposite local warming effect⁷⁸. It is unclear how potential non-local effects would affect the local effects we observe. Modelling experiments with dedicated coupled runs would need to be established to investigate this further, but they would require a fine granularity to adequately represent the local changes in land cover.

A remark is warranted to discuss particular situations where there might be a strong influence of water on the space-for-time substitution signal. For example, trees may be expected to behave differently if planted over wetlands. We deemed that cases involving afforestation in wetland would deserve a specific attention that is beyond the scope of this more general global analysis. For this reason, we kept wetlands as a separate vegetation cover class, whose effect is therefore factored out when considering the transitions of afforestation from herbaceous vegetation to either evergreen forests or deciduous forests. Irrigation, on the other hand, is included in the herbaceous vegetation class. Therefore, if we consider potential afforestation in areas where irrigated crops are presently common, the resulting estimated CFrC change will relate to that of a transition from irrigated cropland to forest. Finally, another place where water could affect the method includes coastlines. Their effect are mitigated by the experimental set-up of the moving window, which requires that more than 60% of the pixels be valid non-water pixels. This generates a buffer around coastlines for which the values are ignored, precisely to avoid these areas where the cloud formation may be dominated by the sea. Local residual coastline effects might still remain is some places, and further methodological refinements could be foreseen to where the method be applied for a regional scope.

Despite having observations coincide with the peak formation of cumulus clouds, clouds might still move. Several methodological points should mitigate the consequences of this effect on our results. First, we work on monthly averaged CFrC, which is then averaged over the 11-year period to describe a multi-annual cloud regime for every grid cell irrespective of all wind conditions (in terms of both strength and direction) encountered during this period. In the unlikely event that at some places the wind climatology has indeed a systematic direction and a sufficient strength to displace the low-level boundary clouds far beyond the underlying forest in the short amount of time they have, the most likely result is a lack of correlation between the land and the atmosphere. This would lead to the same result as what we described in the assumptions for stratus or high clouds: the space-for-time method would lead to random noise. Cloud movements can actually be expected to effectively decrease the spatial gradients (among different neighbouring land cover classes) that we are studying, since these would be stronger in the absence of cloud movements. As a result, we expect that our methodology underestimates the actual magnitude of the changes in CFrC.

Linked to the previous point, a special consideration is needed to the question of scale. The fact that convective clouds are susceptible to move laterally between the moment they are formed above the forest and the moment they are detected by the satellite also means that choosing the right spatial resolution to do the analysis is not directly evident. The finer the pixel, the more likely it is that the clouds have moved to neighbouring pixels. On the other hand, if the pixel is too coarse it may dilute the overall effect, leading further to underestimation of the magnitude. To partly explore the sensitivity of the methodology to scale, we realized some variant experiments. Over Europe, the Cloud CCI project additionally provides the same MODIS-based data product with a finer spatial resolution of 0.02° (instead of the standard 0.05° resolution). Using this specific dataset, two experiment variants have been made. The first consists of applying the space-for-time methodology with the same configuration: i.e. a moving window of 7 × 7 pixels, but here resulting in a finer scale with aggregated pixels of circa 14 km instead of circa 35 km. The second instead considers a 17 × 17 pixel moving window covering an area of circa 35 km, similar to the original experiment using 0.05° data. To do both of these specialized experiments, the land fractional covers had to be generated again from the CCI Land cover data to result in land fractional cover layers at the corresponding spatial resolution of 0.02°. The results of this analysis are summarized in Fig. 6, where it appears clear that the general patterns remain constant, but that the magnitude of the effect can depend on scale.

Confrontation with an alternative method

In order to strengthen confidence in our results based on the space-for-time substitution approach, we also explore an alternative method. This other approach investigates the relationship between land cover change and biophysical effects by quantifying their effects after real land cover changes occur²¹. This method requires factoring out the local effect from climate variability to isolate the effect from the land cover using a similar moving window approach. For this purpose, the difference in monthly CFrC (Δy) between two years at a given location is expressed as the effect of forest cover change (Δy_fcc) plus the residual signal (Δy_res) due to climate variability, so the former can be estimated as such:

To increase the number of samples that can be used with this approach, the calculation is done based on all possible pairwise comparisons of different years within the 2004–2014 period, not just between the first and last year. This is applied independently for each month. We here use the same implementation as a previous study²¹, but we use a forest cover percentage derived from the CCI Land Cover to remain consistent with the space-for-time approach.

While these results would arguably be a more direct estimation of cloud cover, there are two disadvantages that preclude it from using it as the main approach. First, the number of actual changes in land cover is much smaller and harder to identify than the potential changes used in the space-for-time approach, resulting in a strongly reduced number of observations from which to derive an estimation (and thus a necessity to have larger and less reliably stable moving windows). Second, land cover change is a gradual process and thus the effect of the change (e.g. change in cloud cover regime in this case) may not be fully expressed in the temporal interval considered. These reasons explain why the results are not expected to be identical to those obtained from using the space-for-time substitution. Their role here is to add confidence to the main analysis by showing that two different approaches do generally show the same direction in the response of cloud cover, as shown in Fig. 3.

Confrontation with ground observations

To consolidate our CFrC change estimations from satellite we reproduce comparable estimations from ground-based observations of cloud cover from synoptic weather stations. These observations are made by visual examination of the sky under a specific protocol at regular sub-daily timesteps⁷⁹. The global SYNOP database originated from the archive of the European Centre for Medium-Range Weather Forecasts (ECMWF) and assembles a harmonized set of such measurements across the world. However, the density and completeness of records is very variable geographically and temporally, with many more usable records in Europe than in the rest of the world.

To be usable for our purposes, we first need to associate a corresponding forest cover value to each station and filter out those for which cloud formation may be dominated by other factors. We assume here that the CFrC measurement at every SYNOP location relates to the fraction of forest cover within a region of influence defined by a radius of 15 km around each station. We thus associate to each point a value of forest cover fraction based on the Global Forest Change maps⁸⁰ using the Google Earth Engine platform⁸¹. We prefer to use these simpler binary maps since their finer spatial resolution (30 m) is more appropriate for this fine-scale station-level analysis than the thematically richer CCI Land cover maps used for the space-for-time analysis. We similarly attribute to each station a value of water fractional cover related to the same 15 km radius using the global surface water maps⁸². This information is then used to filter out stations with more than 5% of water fractional cover where the micro-climate conditions may be strongly impacted by the presences of water bodies. A second filtering is based on the topography, using the same criteria as the one for the space-for-time analysis above. Finally, we only retain stations that have a minimum number of records corresponding approximately to 7 years of CFrC records between 2004 and 2014 in order to ensure we are not focusing on the situation of a single year.

The next step is to group stations in pairs that are comparable but which relate to variable percentages of forest cover. The main criteria we use to find such pairs is the distances between stations. A minimum distance is set to 30 km to ensure that there is no overlap between the regions of influences of the two stations, assuming that this equally prevents auto-correlation in terms of the sky viewed at each point. A maximum distance of only 100 km is tolerated to avoid the fact that the large-scale weather situation differs considerably in terms of the impact on cloudiness above the two locations. With these constraints, the resulting pair sites are only found in Europe where the network is dense enough (see Fig. 4a). The difference in CFrC Δy_synop between both stations within each pair is computed, along with the differences in forest cover fraction Δx_synop.

The final step is to make the Δy_synop compatible with the Δy from satellite. The challenge is that the satellite estimation Δy_A→B refers to a full transition from a land cover class to another (i.e. 100% cover of A to 100% of B), while there are no single pair of SYNOP points showing a transition of more than 40% of forest cover. To harmonize the two concepts, we extrapolate the theoretical value of a full cover transition from no forest to full forest by fitting a linear regression between Δy_synop and Δx_synop using the values from all pairs of points across Europe. As we force the relationship to pass through the origin, the slope represents the expected change in cloud cover for a full transition from no forest to forest. This regression is applied separately for every hour and for every month (resulting in Fig. 4b). To compare with the satellite estimates, the SYNOP values corresponding to 14:00 (close to the overpass time of the Aqua satellite) are compared to the total average of the satellite estimates falling within 50 km of the centroid of each station pair. This is done both with the original satellite data at 0.05° spatial resolution and the dataset with a refined 0.02° spatial resolution, both of which are summarized in Fig. 4c.

Changes in the surface energy balance

We propose a final analysis confronting the potential changes in CFrC with equivalent changes in different components of the surface energy balance following potential afforestation. This is possible by leveraging on a previous study²², where the same space-for-time framework was used to quantify the changes in net radiation (R_n), in latent heat flux (LE) and in the combination of sensible and ground heat fluxes (H + G). These datasets are fully described in⁷⁶, in which two different levels of vegetation transitions are provided. The single transition between the class ‘forest’ and the class ‘crops and grasses’ of the more generic dataset (IGBPgen) is used here. We retain only the values where there are common spatio-temporal records of both datasets, i.e. those for changes in components of the surface energy balance (LE, R_n and H + G) and for changes in CFrC. By looking at how changes in these variables associated with afforestation co-vary with our estimated changes in CFrC, we can catch some insights on the underlying processes behind cloud formation above forests.

This analysis brings all estimates of CFrC change in a common system of physical coordinates irrespective of their geographical location or their seasonal behaviour (Fig. 7). To facilitate the interpretation, we exclude snow-covered areas, by using a threshold on the reduction in surface albedo following afforestation (Δα ≥ −0.1; such large differences in surface albedo occurs predominantly between dark emerging trees and snow-covered open lands). Once the snow effect is masked out, the general picture shows a clear increase in CFrC, with higher values when the afforestation entails a larger amount of net radiation (see gradients in Fig. 7a, b). When mapping the changes in CFrC with respect to changes in the non-radiative fluxes (Fig. 7c), we see that the large majority of cases where afforestation would occur (i.e. the points shown in the plot), there is either an increase in LE, in H + G or in both, and it is systematically associated to an increase in CFrC. This comforts the notion that cloud formation above forests are correlated to both the injection of extra moisture in the boundary layer and the increased bulk displacement of air generated by the fact the forest surface is darker. It also suggests that depending on where the potential afforestation occurs, one process might dominate over the other.

There are some shortcomings to this analysis that prevent us from drawing strong conclusions from it. First, unlike the CFrC product, the energy balance products are not specific to the early afternoon but instead resume the situation integrated over the entire day. This mismatch in the timing complicates the interpretation of a causal link between changes in surface energy balance and cloud formation, as the former are affected by other processes throughout the day (and night). Second, the original study²² assumed that the local indirect biophysical effects of vegetation cover change, such as the change in cloud cover we observe here, could be neglected in order to close the surface energy balance. Therefore, there is an inconsistency in this analysis, since the consequences on the energy balance of the change in CFrC are not reflected in the coordinates of the plots in Fig. 7. Despite these issues, we still think it is valuable to provide this analysis for completeness as a complementary description of the CFrC signal we have derived.

Source: Ecology - nature.com