Revealing the widespread potential of forests to increase low level cloud cover
Assumptions for the space-for-time substitutionThe 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 temperature20,22 or on the surface energy balance22,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 season43. 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 datasetsThis 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)73.The Cloud CCI55 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 formation43, 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 CCI54, 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 Scheme74. This information is aggregated both spatially and thematically using a specifically designed framework75 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 study22. 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 own60, 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 changeUnder the above-mentioned assumptions, we apply a space-for-time substitution algorithm developed in a previous study22 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 papers22,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:$${bf{y}}={bf{X}}beta$$
(1)
This is equivalent to solving the following system of equations:$$left{begin{array}{ll}{y}_{1}=&{beta }_{1}{x}_{11}+{beta }_{2}{x}_{12}+…+{beta }_{m}{x}_{1m}\ {y}_{2}=&{beta }_{1}{x}_{21}+{beta }_{2}{x}_{22}+…+{beta }_{m}{x}_{2m}\ vdots &\ {y}_{n}=&{beta }_{1}{x}_{n1}+{beta }_{2}{x}_{n2}+…+{beta }_{m}{x}_{nm}end{array}right.$$
(2)
in which the digits of the subscript of x, e.g. xij, 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 xj = 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 conclusions77. To avoid this, we reduce the dimensionality of X through singular value decomposition (SVD) after removing the mean of each column:$$({bf{X}}-{bf{M}})={bf{U}}{bf{D}}{{bf{V}}}^{t}$$
(3)
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, Vz, can then be used to find the appropriately transformed predictor matrix of reduced dimension Z as follows:$${bf{Z}}=({bf{X}}-{bf{M}}){{bf{V}}}_{z}$$
(4)
which can now be regressed onto the CFrC y:$$y={bf{Z}}{beta }_{z}+varepsilon$$
(5)
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:$${beta }_{z}={left({{bf{Z}}}^{t}{bf{Z}}right)}^{-1}{{bf{Z}}}^{t}y$$
(6)
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 (Vz, again, specific to each local analysis).$${{bf{Z}}}_{{rm{p}}}=({bf{P}}-{bf{M}}){{bf{V}}}_{z}$$
(7)
Predicted yp values for each vegetation or land cover type (identified by predicting the appropriately transformed ‘dummy pixels’) are then calculated as:$${y}_{{rm{p}}}={{bf{Z}}}_{{rm{p}}}{beta }_{z}$$
(8)
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 yp predicted for each ‘pure’ vegetation type:$${{Delta }}{y}_{{rm{A}}to {rm{B}}}={y}_{rm{B}}-{y}_{rm{A}}$$
(9)
The uncertainty in the estimation of ΔyA→B can be expressed as a standard deviation using the following expression:$${sigma }_{{rm{A}}to {rm{B}}}=sqrt{{sigma }_{rm{A}}^{2}+{sigma }_{rm{B}}^{2}-2{sigma }_{rm{AB}}}$$
(10)
where ({sigma }_{rm{A}}^{2}) and ({sigma }_{rm{B}}^{2}) are the variances in the estimates of yA and yB, and σAB is their covariance. These variances and covariances are in turn obtained from the covariance matrix, defined from the regression as:$${mathbf{Sigma }}={{bf{Z}}}_{{rm{p}}}{rm{Var}}[beta ]{{bf{Z}}}_{{rm{p}}}^{t}$$
(11)
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-processingA 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-occurrence76, Ic, 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 pA and pB, 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 qA and qB, 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:$${I}_{{rm{c}}}=1-frac{{sum }_{i}min {sqrt{{left({q}_{A}-{p}_{A}right)}^{2}+{left({q}_{B}-{p}_{B}right)}^{2}}}}{{sum }_{i}sqrt{{q}_{A}^{2}+{q}_{B}^{2}}}$$
(12)
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. Ic 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 in76, we retain only pixels with Ic ≥ 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 action56. Here we mask the areas where considerable topographical variation occurs within the 7 × 7 pixel moving window of interest using the same implementation described in76. This involves using 3 different indicators, v1, v2 and v3, 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:$${v}_{1}=frac{1}{n}mathop{sum }limits_{i=1}^{n}{sigma }_{h,i}$$
(13)
$${v}_{2}=| {mu }_{h}-frac{1}{n}mathop{sum }limits_{i=1}^{n}{mu }_{h,i}|$$
(14)
$${v}_{3}=| {sigma }_{h}-{v}_{1}|$$
(15)
For an interpretation of these metrics, readers are invited to read76. These three indicators are combined together in a single layer depicting all pixels satisfying all of the following conditions: v1 More
