Global Unsupervised Assessment of Multifrequency Vegetation Optical Depth Sensitivity to Vegetation Cover

Vegetation optical depth (VOD) has contributed to monitor vegetation dynamics and carbon stocks at different microwave frequencies. Nevertheless, there is a need to determine which are the appropriate frequencies to monitor different vegetation types. Also, as only a few VOD-related studies use multifrequency approaches, it is needed to evaluate their applicability. Here, we analyze the sensitivity of VOD at three frequencies (L-, C-, and X-bands) to different vegetation covers by applying a global-scale unsupervised classification of VOD. A combination of these frequencies (LCX-VOD) is also studied. Two land cover datasets are used as benchmarks and, conceptually, serve as proxies of vegetation density. Results confirm that L-VOD is appropriate for monitoring the densest canopies but, in contrast, there is a higher sensitivity of X-, C-, and LCX-VOD to the vegetation cover in savannahs, shrublands, and grasslands. In particular, the multifrequency combination is the most suited to sense vegetation in savannahs. Also, our study shows a vegetation–frequency relationship that is consistent with theory: the same canopies (e.g., savannahs and some boreal forests) are classified as lighter ones at L-band due to its higher penetration (e.g., as shrublands), but labeled as denser ones at C- and X-bands due their saturation (e.g., boreal forests are labeled as tropical forests). This study complements quantitative approaches investigating the link between VOD and vegetation, extends them to different frequencies, and provides hints on which frequencies are suitable for vegetation monitoring depending on the land cover. Conclusions are informative for upcoming multifrequency missions, such as the Copernicus Multifrequency Image Radiometer.


I. INTRODUCTION
R EMOTE sensing is a useful tool for the regular and global monitoring of the ecosystem's health, vegetation distribution and its dynamics, and changes in global carbon and water cycles. This is paramount to develop climate change mitigation strategies to reduce the global atmospheric CO 2 [1], [2]. The most widely used techniques for vegetation monitoring are based on visible -near infrared vegetation indices (VIS/NIR), such as the Normalized Difference Vegetation Index, which measures the photosynthetic activity and its spatial and temporal changes [3]. Still, VIS/NIR indices are limited by 1) the influence of clouds and aerosols, and 2) the fact that the relationship of these indices with biomass is limited by saturation at high biomass density, as it is only representative of the top layer of the vegetation [4].
Emerging as a complementary tool overcoming these issues, passive microwave remote sensing is nearly transparent to clouds and-although with a coarse resolution-is able to sense the vegetation at different layers and depths, depending on the frequency. In particular, microwave radiometers measure the radiation emitted by the Earth's surface, which is a function of several parameters, including its temperature, soil moisture, soil roughness, vegetation water content, and vegetation biomass and structure [5]. Vegetation effects are represented in radiative transfer models by the scattering albedo and by the attenuation of the vegetation over soil and plant microwave emissions. The latter is measured by the dimensionless parameter vegetation optical depth (VOD), being effective to monitor vegetation response to drought [6]. At the low frequencies (i.e., L-band: 1.4 GHz), the penetration depth of microwaves through the vegetation canopy is greater, sampling the vegetation for most of the canopy layer thickness [7], [8].
Several studies have used the VOD to analyze different vegetation properties, choosing the appropriate frequency depending on which characteristics were to be studied. X-band VOD (X-VOD) has been applied to study gross primary production and This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ evapotranspiration, as it is representative of the photosynthetically active biomass of plants (i.e., canopy leaves [9], [10]). Similarly, Konings and Gentine [11] provided estimates of the degree of isohydricity of plants by using the X-VOD in order to exclude the stem contribution to the retrievals and keep only a VOD signal that is sensitive to the leaf water potential. In contrast, when the full vegetation layer in dense canopies is the subject under study, the application of L-band VOD (L-VOD) is needed to ensure a larger penetration of the measured microwave emissions. L-VOD has been used to study deforestation in tropical forests of South America and Africa [12], [13], it has been related to vegetation height [14], [15], and it has been widely applied to map biomass and to analyze carbon trends (e.g., [16], [17], and [18]).
Regarding multifrequency VOD studies, research in [8] compared the sensitivity of L-, C-and X-VOD to above-ground carbon measured from airborne Lidar in South and Central American forests, showing that L-band is more sensitive to carbon density in the dense tropical forests. However, the authors also indicated that the synergy of multifrequency observations would be appropriate for measuring biomass in less dense canopies, such as grasslands, shrublands, or low forests, in the Andes range. Pringent and Jiménez [19] evaluated the synergy of satellite passive microwave observations between 1.4 and 36 GHz for vegetation characterization over the tropics also showing the potential of a multifrequency approach. Nevertheless, global analyses of vegetation characteristics from multifrequency VOD are still lacking. They are needed to understand which frequencies are appropriate to monitor the vegetation density and water content from the different vegetation types over the Earth surface. This would provide further knowledge on how to study vegetation properties with future multifrequency passive microwave missions, such as the Copernicus Imaging Microwave Radiometer (CIMR), which will operate at L-, C-, X-, Ku-, and Ka-bands [20].
In this study, we aim to qualitatively analyze, at global scale, the sensitivity of VOD at different frequencies (L-, C-, and Xbands) to the vegetation density. To achieve this, an unsupervised global-scale classification of VOD has been implemented by using these frequencies both separately and combined. Results have been compared to land cover classes, which serve here as a proxy of vegetation density. Our research questions are: Which is the qualitative relationship between VOD frequencies and land cover classes? Which VOD frequencies could be more appropriate to monitor vegetation for the different land cover classes? By answering these questions, we will clarify which VOD frequencies are more sensitive to the different land cover classes, in which regions, and to what extent the result is consistent with the fact that lower frequencies are more sensitive to denser canopies.

A. Vegetation Optical Depth
L-band VOD (L-VOD: 1.4 GHz) is derived from the SMOS-IC version 2 product (produced by INRA-CESBIO from the SMOS mission [21]). L-VOD is shown in Fig. 1(a). In this product, both soil moisture and VOD are retrieved simultaneously by using the radiative transfer model L-band microwave emission of the biosphere (L-MEB), where the vegetation layer contributes to the radiative emission at L-band by attenuating and scattering the soil emission and by adding its own contribution to the total radiation measured above the canopy. The SMOS-IC product [22] has the advantage of being as independent as possible of auxiliary data, as it considers the footprints to be homogeneous in order to avoid uncertainties and errors linked to inconsistent auxiliary datasets [21], making it more suitable to perform vegetation studies, such as vegetation seasonality [23], crop modeling [24], and biomass estimation [18], [25]. The SMOS-IC dataset is provided on the Equal-Area Scalable Earth Grid version 2 (EASE2) [26] with a spatial resolution of 25 x 25 km 2 at 30°of latitude. C1-band VOD (C-VOD: 6.9 GHz) and X-band VOD (X-VOD: 10.7 GHz) products, shown in Fig. 1(b) and (c), are derived from the Advanced Microwave Scanning Radiometer 2 carried on the Global Change Observation Mission 1st Water (GCOM-W1) satellite. Soil moisture and VOD are retrieved by using the land parameter retrieval model through a nonlinear iterative procedure by applying the microwave polarization index [27]. The ground resolutions of C-and X-VOD are 35 x 62 km 2 and 24 x 42 km 2 , respectively. The dataset is provided on a 25-km grid [28]. The period covered in this study spans from 2016 to 2018.

B. Land Cover Maps
Two different land cover products are used to study the VODfrequency-land cover relationship as well as to understand how different land cover products can impact the results and their interpretation. On the one hand, the International Geosphere-Biosphere Program (IGBP) has been applied [see Fig. 2(a)]. This is a 17-class land cover dataset obtained with unsupervised classification using data from the Moderate Resolution Imaging Spectroradiometer, and with postclassification refinement. Its spatial resolution is 500 m [29]. On the other hand, the Copernicus Climate Change Service (C3S) provides global 22-class land cover maps at 300 m spatial resolution for 2016 to 2019. Here, the map for 2018 has been applied [see Fig. 2

A. VOD Processing
C-and X-VOD datasets have been linearly interpolated to match the EASE2 25-km grid of L-VOD. Also, four screening steps have been applied: 1) Pixels with a fraction of open water bodies, ice/snow, and/or urban areas larger than 5% have been screened out. 2) L-VOD values have been filtered by removing strong topography as it may impact the angular signature of radiometers brightness temperatures [30]. 3) Since the presence of radio frequency interferences can affect the quality of the retrievals, the RMSE between the modeled brightness temperature (obtained with the L-MEB model) and the SMOS measured brightness temperature was used as an indicator of the retrieval quality. Values with RMSE > 6 K have been screened out [22]. Note that the second and third steps were only applied to SMOS data, as the AMSR2 product is already filtered by those parameters. 4) Outliers for the three VOD products have been removed by 1) computing the differences between raw VOD data and 30-day smoothed VOD data (moving average) and 2) removing values lower/higher than the 10th/90th percentiles of this result. The VOD values have been yearly averaged using both the ascending and descending orbits to remove the VOD diurnal variations due to their sensitivity to the vegetation water content and canopy rain interception. The coefficient of variation of the year time-series has also been calculated, where all the pixels with a coefficient of variation higher than 1 were excluded due to their high dispersion. Moreover, only pixels with a number of samples higher than 50 days per year have been considered in this analysis.

B. Unsupervised Classification Analysis of VOD
A K-means clustering for the three VOD frequency bands individually, and for the combination of the three frequencies, has been applied. For the latter (LCX-VOD), the three frequencies have the same weight, meaning that are equally important for  Table I, to distinguish both datasets (the original IGBP land cover has 17 categories, whereas C3S has 22). Categories not mentioned in Table I were modified based on the closest correspondence between IGBP and C3S. the clustering process. Each VOD cluster is compared with each land cover class. The land cover is used as a qualitative proxy of the density of vegetation. The unsupervised machine learning K-means algorithm allocates each data point to the nearest cluster by finding the smallest Euclidean distance between the input vector and the centroid vector [31]. The number of clusters was selected from a silhouette analysis [32] on K-mean clustering computed from 5 to 9 labels. Silhouette analysis is a goodnessof-clustering index that studies the separation distance between the resulting clusters. It ranges between -1 and 1. Coefficients near to +1 indicate that the sample is far away from the neighbor clusters (i.e., it can be only assigned to one cluster), whereas coefficients close to 0 indicate that the cluster is very close to the decision boundary between two clusters (i.e., its classification is not clear). Negative values indicate that those samples might have been assigned to the wrong cluster. The silhouette analysis displayed in Fig. 3 shows that using five clusters provides an appropriate clustering, with all their silhouette coefficients over 0.6, being greater than those found for 6 to 9 clusters divisions (see Supplementary Material). For the five-cluster configuration, the combination of the three frequencies only shows few negative values. For these reasons, finally, five different clusters were applied to study the relationship between single frequency and multifrequency VODs and the different land cover classes.

C. Reclassification of Land Cover
IGBP and C3S land cover maps have been resampled to the EASE2 25-km grid by assigning to each pixel the dominant class. Only the pixels with a dominant fraction of a single class higher than 60% have been included in the analysis, to guarantee a more representative and homogeneous sample. To compare the same number of clusters and land cover classes, both LC datasets were reclassified into five categories, which encompass all major vegetation types on Earth. Table I tabulates the aggregation of land cover classes according to their vegetation canopy density. Fig. 4(a) and (b) shows the maps of the resulting reclassifications. Note that the homogeneity filter mentioned above removed more pixels in the C3S dataset than in the IGBP land cover dataset. Therefore, the fact that C3S raw dataset is more heterogeneous than IGBP dataset, representing the land with five more classes, causes a larger loss of data when filtered [e.g., no data are available in large regions of North America and Australia; see Fig. 4 Concerning to differences between land cover classifications and the accuracy of the products, some studies [33], [34] have shown that the accuracy of the different land cover maps is below 60%. Part of the differences between IGBP and C3S are also due to their different spatial resolutions. The higher resolution Fig. 3. Silhouette analysis for K-means clustering on sample data with five labels. Each cluster is represented by a horizontal "fin shark" shape. It indicates a decreasing number of pixels from the upper, widest part, to the lower, thinnest part (e.g., in cluster 2, a lot of pixels are closer to cluster 1 than to cluster 3). Refer to the Supplementary Material for a detailed description of the silhouette analysis from 6 to 9 labels. of C3S can partially explain its higher heterogeneity. Therefore, results will be also interpreted and discussed according to differences between land cover products.

D. Performance of the Classification
The resulting clusters (see Section III-B) have been matched to each land cover class (see Section III-C), and interpreted according to vegetation density. The performance of the classification algorithm has been analyzed by comparing the VOD clustering with the land cover types in two steps.
First, the performance has been assessed globally by doing an overall cluster-class fitting analysis (i.e., without evaluating the specific cluster-class pairs performances) for each frequency and for the combination of frequencies. To that goal, the following three metrics have been used.

1) Homogeneity (of VOD clusters):
This measures the normalized conditional entropy of the class distribution given the proposed clustering (h k ). Thus, it serves to evaluate how homogeneous the proposed clustering is. It is computed as 1-h k to fulfill the convention of 1 being desirable (full homogeneity) and 0 undesirable (full heterogeneity).
Here, we express it as a percentage to ease the interpretation. 2) Completeness (of land cover classes): This measures the normalized conditional entropy of the clusters distribution given the proposed land cover classes (h c ). Thus, it serves at evaluating how complete the proposed land cover classes are. It is computed as 1-h c to fulfill the convention of 1 being desirable (full completeness) and 0 undesirable (full incompleteness). Here, we express it in percentage to ease the interpretation. 3) V-measure: It corresponds to an entropy-based measure that evaluates the accuracy by using a combination of homogeneity and completeness. It is computed as the harmonic mean of distinct scores of these two metrics. Here, it is expressed as a percentage to ease its interpretation. Further description of the three overall performance metrics is found in [35].
Second, the specific performance of each cluster-class pair has been calculated for each frequency and for the combination of frequencies. This has been done 1) for both land cover classifications separately, and 2) only considering the pixels with the same land cover label in both classifications in order to test the consistency of the results given some inaccuracies and mismatching in land cover datasets (see Section III-C). To do so, matching matrices between clusters and classes have been computed, with diagonal cells indicating the expected land cover class-VOD cluster matchings. Then, the following two metrics have been used to test the cluster-class performances.
1) The cluster consistency (i.e., for VOD) measures the percentage of a VOD cluster that is formed by the same land cover class. It ranges from 0%, when the cluster is totally inconsistent, to 100%, when the result is totally consistent. 2) The class consistency (i.e., for land cover classes) measures the percentage of a land cover class that is assigned to the same VOD cluster. It ranges from 0%, when the class is totally inconsistent, to 100%, when the class is totally consistent. Finally, note that the goal of the performance analysis is not to determine how accurate is the land cover classification. This is out of context given the low spatial resolution of passive microwave measurements: VIS/NIR sensors are the appropriate tools to this task. Instead, we aim at determining a qualitative correspondence between land cover classes and VOD clusters and frequencies (see Section I).

E. Changes in Vegetation Density and Wetness Between Seasons
Further study of vegetation patterns has been conducted by including seasonal analyses. In that sense, the K-means clustering has been computed separately for VOD averages of the periods December-February, March-May, June-August, and September-November. The changes between March-May and June-August, and between June-August and September-November have been evaluated both globally and regionally. The following five regions, including different land covers and seasonal patterns, have been studied: 1) the US Corn Belt and southern Canada, 2) the Sahel, 3) the Iberian Peninsula, 4) the Miombo woodlands in southern Africa, and 5) boreal forests in Russia.
Results have been compared with the expected vegetation patterns and literature and have been interpreted not only according to vegetation density, but also to vegetation wetness/dryness, as seasonal patterns of VOD are indicative of both magnitudes.
Transitions including the December-February period have been excluded after checking that the screening of snow and frozen ground regions in the Northern Hemisphere involved losing too much information, thus making clustering results not comparable with other seasons. Table II tabulates the three metrics applied to evaluate the performance of the classification method for each VOD frequency and their combination. Results show that K-means clustering performs better at L-VOD for both IGBP and C3S classifications (V-measures = 45.5% and 50.4%, respectively), followed by LCX-VOD (34.8% and 37.4%), and decreasing as frequency increases (C-and X-VOD; V-measures ∼34% and ∼30%, respectively; see Table II). These results are consistent with the enhanced sensitivity (due to low saturation) of the L-VOD in the densest vegetation, in contrast to the VOD at higher frequencies (e.g., [7] and [14]). In addition, the fact that metrics' results for LCX-VOD are in the same range as those for single frequencies (see Table II) suggests that complementarity between VOD frequencies exists, as reported in previous research [8], although it needs further study. Still, it is important to mention that the performance analysis reported in Table II includes the five VOD clusters at global scale. Hence, these results are not able to indicate the best-suited band for monitoring each type of vegetation.  Fig. 4). These patterns are general and show relevant exceptions in some frequency-land cover correspondences which are discussed in next sections (e.g., northernmost boreal regions at L-band are linked to the shrubland-like cluster 4 in the IGBP, whereas the area is partially covered by forests in the C3S). Fig. 5 shows the matching matrices between VOD clusters and land cover classes together with their consistencies. It also shows that cluster 3, which has a fragmented distribution in Fig. 4, is mostly linked to savannas according to the IGBP classification. Fig. 5 confirms the expected relationship between clusters and land covers. In that sense, the diagonal cells in Fig. 5 match with the highest number of pixels of each cluster in 85% of cases (17 out of 20) for the IGBP classification. In the case of C3S, this percentage is lower (12 out of 20; 60%), mainly because there is a low number of samples of savannah pixels in the C3S land cover classification. Nevertheless, a qualitative trend is appreciated with two low vegetation density clusters (3 and 4) being assigned also to a low vegetation density land cover (grasslands/croplands). Therefore, a general cluster-land cover linkage can be defined as in Table III.

B. Global Relationship Between VOD Clusters and Land Cover Classes
Importantly, this relationship is only used for a general interpretation of the results and requires further investigation as it is not homogeneous for all land cover-frequency pairs. Hence, results shown in Figs. 4 and 5 are hereafter detailed and Bold numbers indicate which land cover class is dominant in each VOD cluster. Light blue is used to highlight which diagonal cells contain the dominant land cover (i.e., in which cases the VOD cluster and the expected land cover class match). The row summary (far right column) displays the consistency of VOD clusters (Clust. cons., i.e., the percentage of a VOD cluster which is formed by the dominant land cover). The column summary (bottom row) displays the consistency of land cover classes (Class cons., i.e., the percentage of a land cover class which is assigned to the corresponding VOD cluster).
interpreted. The explanation is structured in two sections ordered by increasing vegetation density to ease the interpretation.

C. Multifrequency and High-Frequency Approaches are Suited to Sense Low-Vegetation Canopies and Savannahs
Qualitatively, the spatial distribution of clusters 1 and 4 is linked to low vegetation density land cover classes. They are mainly found in grassland and open shrubland regions in central Australia, the Sahel, southern Africa, central and southern Argentina, central Asia, and-partially-the Great Plains in North America (see Fig. 4). Also, at L-band, the shrubland-like cluster 4 is found in African subtropical savannahs, as well as in northernmost boreal regions dominated by either forests or shrublands depending on the land cover map studied. These two cases are discussed in Section IV-D.
According to Fig. 5, the grass-like cluster 1 is mostly built by pixels of the class cropland/grassland (∼50% for IGBP and ∼60% for C3S), and by pixels of savanna and shrubland in a lower proportion. The cluster consistency of cluster 1 is very similar among frequencies, being moderate to high (52% to 68%, Fig. 5). It is the highest for the LCX-VOD combination, according to the IGBP, and for the L-VOD, according to the C3S (see Fig. 5). The low density of grasslands and croplands suggests that an improved sensitivity of high and multifrequency approaches should be expected, and that L-band should be less prone to study this kind of vegetation. The seasonal changes shown in Figs. 6 and 7 and the Supplementary Material confirm this fact. On the one hand, grass-like (cluster 1) and shrub-like (cluster 4) vegetation increase their density in the Sahel both for the beginning of the rainy season (June) and until the last rain dates in the region (October). This is detected in a much larger extension by X-band than by L-band (see Figs. 6 and 7). Spatially, similar patterns to those of X-band are found for C-and LCX-bands (see Supplementary Material). On the other hand, in the Iberian Peninsula, the summer drying trends captured by L-band are specifically found in the southwestern region, dominated by savannahs. In contrast, X-, C-, and LCX-bands detect losses of vegetation density/wetness in the entire Iberian Peninsula, including the central regions of Spain where large extensions of wheat crops are harvested in June and July (see Fig. 6). At the X-band frequency, this is detected as a transition from a savannah-like cluster 3 to a grass-like cluster 1, likely due to some saturation of X-band when the crop fields are wet during spring. At L-band, the crop harvesting regions are not detected [see Fig. 6(a)].
Still, the behavior of croplands is complex and depends on their density. In the large crop extensions of the US Corn belt and southern Canada, all frequency bands detect the growth of vegetation between spring and summer (see Figs. 6 and Supplementary Material). Changes detected by X-band (35%) are larger in extension than those detected by L-band (22%), but the latter has shown a good capacity to capture crop phenology and to model crop yield in the region [36], [37].
Also, it should be noted that grasslands show a complex behavior at L-band, since the presence of a litter layer of dead grass below the green vegetation can have a disproportionate effect on the L-band emission if wet [38]. This could impact the interpretation of results in grasslands, which are also limited by uncertainties between IGBP and C3S products and by lack of homogeneous data for the latter [especially in grasslands and shrublands in Australia and Argentina for the C3S classification; see Fig. 4(a) and (b)].
Concerning to cluster 4 (shrub-like vegetation), for C-, X-, and LCX-VOD, it is dominated (>55%) by shrublands in the IGBP classification, and by croplands/grasslands in the C3S classification (where shrublands account for ∼30% of the cluster). For these frequencies, according to the IGBP, the spatial patterns of clusters 1 and 4 distinguish well the distribution of grasslands and shrublands, respectively, in Australia, southern Africa, and Argentina. At these frequencies, cluster 4 is also found in the driest areas of the Asian steppes and the Sahel (see Fig. 4). In that sense, note that the shrub-like cluster 4 is linked to the driest and sparsest vegetation (closer to the deserts), whereas the grass-like cluster 1 is its natural continuation in the biogeographic gradient toward more humid climates (e.g., from north to south in the Sahel, and from inner to outer regions in Central Asia and Australia). Cluster 4 is, thus, mostly representing open shrublands, with the lowest vegetation density.
At L-band, instead, shrubland pixels in Central Asia, Argentina, Australia, and the Sahel are aggregated into the widely extended, grassland-dominated cluster 1 (see Figs. 4

and 5).
This is relevant as it shows a poor capacity for discriminating shrublands from grasslands at L-band (shrubland class consistency of 22% for IGBP and 6% for C3S; see Fig. 5). In contrast, the consistency of cluster 4 is the highest for the X-VOD product in the IGBP analysis. In the C3S classification, a lower agreement (∼30%) is found for the high-frequency bands and the combination, but it is still much higher than that for L-band (6%).
To sum up, the better suitability of high frequencies (and especially X-VOD) in open shrublands is consistent with their lower penetration capacity. In contrast, L-VOD can penetrate deeply the vegetation and remains unsensitive to differences between grasslands and shrublands, which suggests that the higher penetration of this frequency may reduce its sensitivity to vegetation properties (e.g., density) in short canopies.
Finally, cluster 3 is partially linked to savannahs according to the IGBP (∼45% in the class consistency; see Fig. 5). Nevertheless, its distribution is fragmented and dominates savannahs and shrublands in northeastern Brazil, some regions in Africa and Europe that are linked to low forests, savannahs, and/or transition regions and, mainly at X-band, an important part of the US Corn Belt (see Fig. 4). No analysis is feasible with the C3S land cover classification due to a very low consistency of cluster 3 at all frequencies (≤3%). This is due to a much lower number of savannah samples if compared with other land covers for C3S (see Fig. 5 and Section IV-E). For the IGBP product, the consistency of cluster 3 increases when the clustering algorithm combines the three frequency bands (LCX-VOD: 48%; see Fig. 5). Slightly lower cluster consistency values are found at C-band (42%) and at L-and X-bands (39% and 37%, respectively).
The analysis of cluster changes in African savannahs (Miombo woodlands) between March and August reports drying  trends at X-band, and no changes at L-band, consistently with the coupling of water storage and leaf phenology in the region (see [23, Fig. 3.c.i]). At X-band, the changes from cluster 2 (very dense forest) to cluster 5 (other forest), and from cluster 5 (other forest) to cluster 3 (savannahs), show how drying patterns lead to the reduction of saturation at this frequency.
Concerning the June-August to September-November variations, X-band clusters keep showing vegetation density/wetness loss, which is unexpected according to the patterns reported in [23]. This might be caused by no recovery of leaf greenness in these periods due to extreme droughts in the region (values of the Standard Precipitation Evapotranspiration Index, SPEI, are close to -2 at the 3-month time scale during the study periods [39]).
Overall, results are coherent with the fact that high-frequency data are more sensitive to leaves than L-VOD, which is more sensitive to stems and other woody components [40]. These results suggest that multifrequency approaches and high frequencies could improve our capability to estimate vegetation properties from passive microwave sensors in low vegetation, especially in shrublands. The use of a combined LCX-VOD data synergy is promising, with its highest sensitivity being found in savannahs. Fig. 4 shows that cluster 2 is qualitatively linked to the densest forests in tropical regions, regardless of the frequency band and the land cover. Instead, cluster 5 is linked mainly to IGBP boreal forests at L-band, and to African subtropical savannahs at LCX-, C-, and X-bands. In addition, it must be noted that the denseforest cluster 2 is also observed in boreal forests at LCX-, C-, and X-bands, with larger extension with increasing frequency, and that this cluster also dominates the forest-to-tundra transition of the northernmost latitudes at X-band (see Fig. 4). Also, at X-band, boreal forests show increasing cluster 2 extension when the vegetation becomes denser (June-August; see Fig. 6).

D. Relationship Between Frequencies and Vegetation Density Through Forests and Savannahs
In contrast, it is worth noting that, according to the C3S classification, northernmost boreal forests at L-band are dominated by a low-density vegetation cluster (cluster 4; see Fig. 4). The mismatch between land cover maps in this transition region is likely due to the presence of sparser and lower trees if compared with the densest taiga regions. Also, the shrub-like cluster 4 (low-density vegetation) is clearly dominant in subtropical African savannahs at L-band.
These contrasting patterns provide evidence of the relationship between frequency and vegetation density: boreal forests (which are not as dense as tropical ones) are sensed as very dense canopies with increasing frequencies, due to saturation, but the forest-to-tundra transition is sensed homogeneously as low-density vegetation by the lowest frequency, due to its higher penetration capacity (see Fig. 4). Similarly, high frequencies saturate in savannah regions sensing forest-density vegetation, whereas the highest penetration of L-band ignores the presence of tree canopies in these areas (see Fig. 4).
A global-scale quantitative analysis of this relation is reported in Fig. 5. Cluster 2 is dominated by tropical forests. Decreasing cluster consistency is observed for increasing frequencies (99% for L-VOD, ∼77% for LCX-VOD, ∼72% for C-VOD, and ∼52% for X-VOD). This is regardless of the land cover map applied. In that sense, at C-and X-VOD, the construction of this cluster is mixed with other types of forests, and with shrubland to a lesser extent in X-VOD. As in the qualitative analysis of Fig. 4, the effect of saturation at high frequencies is also self-evident here, confirming the need of L-VOD for vegetation monitoring in tropical rainforests [7], [8].
In contrast, the consistency of the tropical forest land cover class is higher for X-VOD and LCX-VOD (∼73%) and decreases for L-VOD (69%) and C-VOD (66%). Hence, the spatial patterns of the highest frequency and the combination of frequencies are similar to those of optical-infrared-based land cover maps. Note that, when comparing optical-infrared vegetation indices to X-VOD, they show a high correlation, and a lack of saturation. In contrast, when comparing these indices to L-VOD, they saturate and show lower correlation (e.g., see [41,Fig. 4]).
Concerning to the cluster 5 (lighter forest density), Fig. 5 shows that it has the highest cluster consistency for LCX-VOD according to the C3S land cover, and for L-VOD according to the IGBP. Cluster consistency decreases with increasing frequencies for both land cover maps, in line with the penetration depth (see Fig. 5).
In addition, note that nontropical forests are formed by a mix between evergreen and deciduous forests, and by both broadleaf and needleleaf forests. Future studies should address time and spatial dynamics in the relationship between vegetation density and microwave frequencies in these forests, which could reflect time-variations in absorption and volume scattering effects on VOD [42], [43].

E. Limitations and Future Work
The main limitations found in this study correspond to the differences between both land cover datasets, and to the nonperfect accuracy of land cover products. Further, the generalization of IGBP and C3S from 13 and 19 classes (water bodies, bare soils, and snow were excluded), respectively, to five classes, naturally led to the loss of the ability to describe detailed land cover characteristics. For example, open shrublands are much less dense than closed shrublands. The same occurs with savannahs and woody savannahs, where the latter may be more sensitive to L-VOD as the volume of woody structure is higher and because woody savannas might be covered with forest canopy up to 60%. On the other hand, nonwoody savannas are more sensitive to higher frequencies, as these areas are regions of transition between forests and low-density vegetation zones, such as grasslands.
In addition, various studies [33], [34] have shown that the accuracy of different land cover maps is below 60%. Some mismatching regions can be found, for instance, in the north of Russia, where IGBP-shrubland pixels are classified as forests by the C3S dataset, or in northeastern Brazil between savannah (IGBP) and shrublands (C3S), to name a few. Finally, the aggregation needed due to the different spatial resolutions of the land cover and the VOD products mixed different types of vegetation, explaining a reduced accuracy of the matching matrices. It must be recall again that, for this study, the absolute accuracy of the land cover-VOD matches is not relevant, whereas the focus is instead on the relative comparison across frequencies and land cover types.
To assess the impact of the inconsistencies in the results, Fig. 5 has been replicated by only considering those pixels where IGBP and C3S land cover labels are the same (see Supplementary  Material). Results are consistent in general terms. In that sense, diagonal cells in the Supplementary Material match with the highest number of pixels of each cluster in 60% of cases (12 out of 20) compared with 85% for IGBP alone and to 60% for C3S alone (see Section IV-B). Importantly, the results show again the higher cluster consistency for tropical forests at L-band, the decreasing consistency with increasing frequency in this forest type, or the similar behavior of all bands in the grasslands/croplands category (see Supplementary Material), thus driving to the same main conclusions that have been derived from the separate IGBP and C3S analyses. Still, note that the main inconsistency between both land cover classifications relies on the low number of pixels simultaneously classified as savannahs (n ∼ 500) due to a low sample in this land cover for C3S (n ∼ 1000) in front of IGBP (n ∼ 20 000).
Also, a lower sample is found for shrublands (n ∼ 4000) if compared with C3S (n ∼ 12 000) and IGBP (n ∼ 23 000). Consequently, a low number of pixels are assigned to the VOD clusters matching these land cover types, thus suggesting that the separate analysis presented in Fig. 5 shows greater completeness than that using only the coincident IGBP-C3S pixels (see Supplementary Material), which has been used complementarily to confirm the consistency of the results.
This study has shed light on the qualitative relationship between vegetation types/density and VOD frequencies at a global scale, showing that distinct density-frequency patterns emerge in most biomes and suggesting the application of highfrequency and multifrequency approaches to sense low and mid-density vegetation covers. This is, therefore, a starting point for globally assessing which frequencies and in which regions are more appropriate to sense different vegetation canopies. To that goal, future work could take advantage of biomass [44], new global canopy height models [45], and/or future missions. In that regard, active microwave measurements from the NISAR (L-band, [46]) and BIOMASS (P-band, [47]) missions, which are planned to be launched in 2023, will provide further information to determine the capability of different low frequencies to study forests worldwide. Furthermore, passive multifrequency microwave information will be available from the Copernicus Imaging Multifrequency Radiometer (CIMR, [48], planned for 2029) to precisely quantify the suitability of different frequencies and of multifrequency approaches to sense different vegetation canopies.

V. CONCLUSION
This research provides an unsupervised classification of L-, C-, and X-VOD bands and of the multifrequency combination of them, and qualitatively compares its results with land cover maps. Because land cover maps are used here as proxies for vegetation density, a qualitative assessment of the suitability of different VOD frequencies to sense different vegetation densities is provided at a global scale. The results derived from the present work complement quantitative approaches investigating the link between VOD and biomass and extend them to different frequencies, to potential multifrequency synergies, and to all major vegetation types in Earth.
The clustering analysis confirms that using L-VOD is more reliable for vegetation density monitoring at global scale. Nevertheless, L-band is distinctly best-suited for vegetation monitoring only in dense canopies. In contrast, the spatial distribution of C-and X-VOD clusters shows that these frequencies are tightly linked to areas with very low vegetation density, such as grasslands and shrublands. Medium-low density vegetation areas, such as savannas, can be sensitive to either L-VOD or C-VOD, but the performance of the combination of three frequency bands is considerably better, suggesting that multifrequency approaches are the most indicated in this biome.
Importantly, a relationship between vegetation density, microwave frequencies, and their penetration capacity is found across most biomes. On the one hand, at the lowest frequency (L-band), some denser canopies are classified (i.e., sensed) as lighter ones (e.g., some boreal forests, savannahs and grasslands are classified as open shrublands). On the other hand, at the highest frequencies (C-and X-bands), similar canopies are classified as even denser ones (e.g., boreal forests as tropical ones). The physical explanation of these results is based on the fact that lower frequency bands capture the attenuation of soil emissivity due to vegetation as it passes through the whole canopy, whereas higher frequencies capture the emission of the upper layers of the vegetation canopy, such as leaves and stems, and saturate in dense vegetation conditions. The higher penetration at L-band entails a reduced discrimination of vegetation densities in short canopies. These results are also confirmed by the increased capacity of high frequencies to sense short vegetation greening in the Sahel, by crop harvesting in the Iberian Peninsula being detected by X-VOD but not by L-VOD, and by saturation of the X-VOD when the summer biomass peak in boreal forests occurs.
Overall, this study shows that the use of different microwave frequency bands improves or complements the estimation of vegetation properties, such as density. The results provide hints on which frequency is more suitable for such estimations depending on the land cover. This is especially relevant in semiarid regions, where the applicability of multifrequency approaches seems more appropriate, because these regions account for 40% of global vegetation and drive important global carbon cycle variations. The results presented are informative for future vegetation studies relying on upcoming multifrequency missions, such as the CIMR.