Global Change in Terrestrial Ecosystem Detected by Fusion of Microwave and Optical Satellite Observations

The detection of global land change via satellite observation is a major challenge in improving the understanding of global environmental change. In this study, we develop a new vegetation index which can be used as a proxy for the fractions of tree canopy and short vegetation, based on the simple linear regression between microwave vegetation optical depth (VOD) and optical leaf area index (LAI). Although we use no high-resolution reference data, the newly developed vegetation index successfully detects global land change which has been reported by previous estimations based on high-resolution reference data. We find that the relationship between VOD and LAI is non-stationary and the temporal change in the VOD-LAI relationship is an important signal for detecting global change in the terrestrial ecosystem.


Introduction
The detection of global change in terrestrial ecosystems is crucially important for understanding global environmental change. Vegetation indices generated from optical satellite observations (e.g., [1,2]) have been intensively used to monitor changes in terrestrial ecosystems. For example, Song et al. (2018) [3] developed a yearly global vegetation product in which every land pixel is characterized by its percent cover of tree canopy (>5 m in height), short vegetation, and bare ground, based on the optical-infrared Normalized Difference Vegetation Index (NDVI). They successfully detected long-term global land change. Wang et al. (2020) [4] analyzed the long-term trend of global Leaf Area Index (LAI) and developed a prognostic model to simulate it.
The limitation of optical satellite observation is that it is sensitive only to the photosynthetically active part of the terrestrial biomass (i.e., leaf biomass). It is not straightforward to retrieve the amount of the photosynthetically inactive part of the aboveground biomass (i.e., woody biomass) from the signal of optical satellite observation. To distinguish tree canopy (in which woody biomass is dominant in the total aboveground biomass) and short vegetation canopy (in which photosynthetically active biomass is dominant) in optical satellite observations, complex supervised machine learning with high spatial resolution reference data such as LANDSAT observation is usually required (e.g., [3]). To overcome this limitation, microwave remote sensing has been recognized as an alternative approach to monitor global change in the terrestrial ecosystem. Microwave Vegetation Optical Depth (VOD) is sensitive to vegetation water content, which includes information on both leaf and woody biomass (see [5,6] for the detailed description of VOD and the retrieval algorithms). Liu et al. (2015) [7] successfully monitored the long-term trend in global aboveground biomass carbon by microwave VOD data. Zhou et al. (2014) [8] indicated that microwave VOD was more sensitive to the degradation of the Congo rainforest by drought than the optical vegetation index.
Combining the microwave and optical vegetation indices has the potential to improve the understanding of vegetation dynamics. Jones et al. (2013) [9] found that the post-fire recovery of NDVI was much faster than that of microwave VOD, which indicated the faster recovery of short vegetation as compared to woody biomass after extreme wildfires in Alaska and Canada. Van Dijk et al. (2013) [10] found that in the Australian millennium drought, no decreasing trend could be detected in the optical vegetation indices, while microwave VOD significantly declined. By numerical simulation,   [11] attributed this difference to the higher resilience of short vegetation to drought than woody biomass. Tian et al. (2017) [12] detected change in woody biomass across global tropical drylands by combining microwave VOD and NDVI. These studies showed the high potential of fusion of microwave and optical satellite observation to distinguish the dynamics of tree canopy and short vegetation.
We aim to develop a new method to detect global land change by combining microwave VOD and optical LAI. We reveal that the simple linear regression between microwave VOD and optical LAI provides a useful proxy of the fractions of tree canopy and short vegetation. The proposed vegetation index successfully detects global change in the terrestrial ecosystem from 2003 to 2019 without any additional reference data.

Data
The data we used are summarized in Table 1. We used the microwave VOD product based on the Land Parameter Retrieval Model (LPRM) provided by the National Aeronautics and Space Administration (NASA). The description of the LPRM algorithm can be found in [5,13]. Microwave VOD was retrieved from C-band (6.9 GHz) brightness temperature observed by the Advanced Microwave Scanning Radiometer for Earth observation system (AMSR-E) and AMSR2 (the reason why we selected C-band can be found in the next section). The VOD products from both AMSR-E (from June 2002 to September 2011) and AMSR2 (from July 2012 to December 2019) were used. We used only night scene data, in order to reduce the effect of a surface temperature bias in the retrieval algorithm [14]. The spatial resolution of this dataset is 0.25 • . The temporal resolution of this dataset is approximately 2-day, and we resampled it to 8-day by averaging data every 8 days to match the temporal resolution of the satellite LAI data. Note that the VOD data are completely independent from the LAI data. The retrieval algorithm does not use any information from optical remote sensing.
We used the optical LAI product processed by Ichii et al. (2017) [15]. Ichii et al. (2017) [15] applied detailed quality control to MODerate resolution Imaging Spectroradiometer (MODIS) onboard Terra and Aqua satellites LAI L4 data (MCD15A2H; [16]). The spatial and temporal resolutions are 0.25 • and 8-day, respectively. Although a finer temporal resolution is ideal, it is difficult to observe LAI in a finer temporal resolution since the optical satellite observation is strongly affected by atmospheric conditions [15,16].  [3,18] As land cover data, we used the MODIS land cover climate modeling grid (MCD12C1) version 6 data product [17]. This land cover product has 17 classes. The temporal resolution of this dataset is yearly. The original spatial resolution of this dataset is 0.05 • and we resampled it to 0.25 • by the nearest neighbor approach.
We used the Vegetation Continuous Fields (VCF) data product v001 [18]. In this data product, every land pixel is characterized by its per cent cover of tree canopy (>5 m in height), short vegetation, and bare ground. The VCF data were developed by Song et al. [3]. The spatial Remote Sens. 2021, 13, 3756 3 of 13 and temporal resolutions are 0.05 • and yearly, respectively. When we directly compared VCF data and our proposed index pixel by pixel, we resampled the data to 0.25 • by box averaging.

Methods
The flowchart of this study can be found in Figure 1. The Microwave and Optical Fusion Approach (MiOFA) was originally proposed by Sawada et al. [19] to improve the ability to retrieve surface soil moisture and vegetation water content from C-band brightness temperature. The MiOFA has been applied and evaluated in field observation sites [19][20][21]. In this paper, we used the equations proposed in [19] to generate a useful index to quantify changes in terrestrial ecosystems, while the original paper [19] used the same equations to improve the estimation of surface soil moisture. Although the equations used in the MiOFA are empirical, they have been used in many previous works and their parameters are well-defined. In the MiOFA, the relationship between VOD and total aboveground Vegetation Water Content (VWC) can be formulated as: where b is a parameter and σ quantifies the seasonally unchanged contribution to VOD, which is mainly attributed to surface soil roughness. Because it is difficult to distinguish the effects of vegetation and surface soil roughness on brightness temperature, the existing VOD data include the effect of surface soil roughness (e.g., [19][20][21][22][23]). Note that our "VWC" includes only the aboveground vegetation water content, which is affected by the seasonal cycle of vegetation dynamics. The seasonally unchanged VWC unaffected by the vegetation phenology as well as surface soil roughness is recognized as σ in our model. The important assumption of the MiOFA is that b-parameter for C-band VOD is time-invariant and speciesindependent.   [20] and   [21] supported this assumption by a field experiment. Although Jackson and Schmugge (1991) [24] originally formulated this parameter as a species-dependent variable, they also revealed that the species-dependency is minimal when the longer wavelength's microwave such as C-band is used. Although bparameter was assumed to be the function of land cover in some satellite soil moisture retrieval algorithms, the variability of their specified b-parameter is reasonably small (e.g., [25]).
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 15 temporal resolution of this dataset is yearly. The original spatial resolution of this dataset is 0.05° and we resampled it to 0.25° by the nearest neighbor approach. We used the Vegetation Continuous Fields (VCF) data product v001 [18]. In this data product, every land pixel is characterized by its per cent cover of tree canopy (>5 m in height), short vegetation, and bare ground. The VCF data were developed by Song et al. [3]. The spatial and temporal resolutions are 0.05° and yearly, respectively. When we directly compared VCF data and our proposed index pixel by pixel, we resampled the data to 0.25° by box averaging.

Methods
The flowchart of this study can be found in Figure 1. The Microwave and Optical Fusion Approach (MiOFA) was originally proposed by Sawada et al. [19] to improve the ability to retrieve surface soil moisture and vegetation water content from C-band brightness temperature. The MiOFA has been applied and evaluated in field observation sites [19][20][21]. In this paper, we used the equations proposed in [19] to generate a useful index to quantify changes in terrestrial ecosystems, while the original paper [19] used the same equations to improve the estimation of surface soil moisture. Although the equations used in the MiOFA are empirical, they have been used in many previous works and their parameters are well-defined. In the MiOFA, the relationship between VOD and total aboveground Vegetation Water Content (VWC) can be formulated as: (1) where b is a parameter and quantifies the seasonally unchanged contribution to VOD, which is mainly attributed to surface soil roughness. Because it is difficult to distinguish the effects of vegetation and surface soil roughness on brightness temperature, the existing VOD data include the effect of surface soil roughness (e.g., [19][20][21][22][23]). Note that our "VWC" includes only the aboveground vegetation water content, which is affected by the seasonal cycle of vegetation dynamics. The seasonally unchanged VWC unaffected by the vegetation phenology as well as surface soil roughness is recognized as in our model. The important assumption of the MiOFA is that b-parameter for C-band VOD is timeinvariant and species-independent.   [20] and Sawada et al., (2017b) [21] supported this assumption by a field experiment. Although Jackson and Schmugge (1991) [24] originally formulated this parameter as a species-dependent variable, they also revealed that the species-dependency is minimal when the longer wavelength's microwave such as C-band is used. Although b-parameter was assumed to be the function of land cover in some satellite soil moisture retrieval algorithms, the variability of their specified b-parameter is reasonably small (e.g., [25]).  Optical LAI products essentially observe green and moist leaves, so that optical LAI can be the proxy of moisture in leaves [26,27]. However, since optical LAI is insensitive to water in woody biomass, species-dependent parameters are needed to directly relate optical LAI to total aboveground VWC. Paloscia and Pampaloni (1988) [26] proposed the empirical relationship between VWC and optical LAI, which was supported by the in-situ observations (e.g., [21]): where y is a species-dependent parameter which cannot be directly observed by satellite sensors. The approximation of Equation (2) is mathematically valid only when y is significantly greater than LAI. The quasi-linear relationship between satellite observed optical vegetation indices and in-situ VWC has been found in many previous works, which were synthesized in [27]. Although Paloscia and Pampaloni (1988) [26] and Gao et al. (2015) [27] analyzed only crops, Sawada et al. (2017) [21] revealed that the linear relationship shown in Equation (2) can be extended to tree canopy with moderate LAI. It should also be noted that "VWC" in the Equation (2) only includes the aboveground vegetation water content affected by leaf phenology, which justifies that VWC approaches to zero as LAI is close to zero and is consistent to the Equation (1). The y-parameter can be recognized as the contribution of leaf biomass to the total aboveground VWC. The lower y indicates the larger fraction of short vegetation since the small increase of LAI induces the large increase of total aboveground VWC. The higher y indicates that there is a larger amount of tree canopy. The y-parameter can be a good index of canopy biomass structure. By substituting Equation (2) to Equation (1), we obtain Equation (3).
The MiOFA proposes to estimate the parameters of canopy biomass structure (b/y) and surface soil roughness (σ) by the linear regression between microwave VOD and optical LAI. In this study, we propose to simply use b/y, the slope of the linear regression, as the newly developed vegetation index. In the original paper of the MiOFA [19], the σ-parameter was focused on to improve the retrieval of surface soil moisture. In this paper, we focused on the slope of the linear regression, b/y, as a vegetation index. Because we assumed that b is a constant parameter, we hypothesized that tree canopy (short vegetation) is dominant in pixels with lower (higher) b/y.
By performing the linear regression between microwave VOD and optical LAI, we developed the yearly global b/y dataset from 2003 to 2019. Since we had no C-band brightness temperature observations in the transition period from AMSR-E to AMSR2, we did not estimate b/y in 2011 and 2012 (see also Section 2). There are long-term and temporally continuous VOD datasets derived by merging the multi-sensor VOD products (e.g., [28]). However, the assumption of the constant b-parameter is correct only in the lower frequencies such as C-band, so we avoided the use of VOD retrieved and corrected by brightness temperature with higher frequencies in the main results. The characteristics of b/y retrieved by C-band VOD from [28] can be found in the Supplementary Materials. Although no significant inconsistency between AMSR-E and AMSR2 has been found in the previous calibration studies (e.g., [29][30][31]), we found that the VOD from AMSR2 was systematically higher than that from AMSR-E. We corrected this bias and the bias correction method can be found at Text S1 in the Supplementary Materials. It should be noted that the retrieved b/y and its trend were not significantly affected by this bias, so we could arrive at the same conclusion as shown below even if this bias is ignored. Because the unseasonal bias is absorbed into σ (see Equation (3)) in our algorithm, it does not affect the trend of b/y. We showed b/y in the pixels where the slope of the linear regression is positive and statistically significant with a 95% confidence level. In the excluded pixels, the assumptions to derive Equations (1)-(3) may not be applicable, so we did not show b/y in those pixels. However, those excluded pixels are mainly located in densely vegetated rainforest and desert, and we could obtain b/y in the other land cover types.
In summary, we performed the linear regression between MODIS LAI and biascorrected AMSR-E/2 VOD in every land pixel. The slopes and intercepts of this linear regression are recognized as b/y and σ (see Equation (3)), respectively. We will show and discuss the global distribution of the b/y index. Our b/y index quantifies the contribution of leaf biomass to total aboveground VWC. It is currently impossible to directly evaluate this leaf biomass contribution to total aboveground VWC by the independent data. Instead, we will evaluate the usefulness of our b/y index by comparing it with the global land Remote Sens. 2021, 13, 3756 5 of 13 cover data. We performed the Kolmogorov-Smirnov test to evaluate if the distribution of b/y in one land cover type can be distinguished from that in another land cover type using MODIS land cover data. We also compared the spatial distribution of our b/y index with VCF data to evaluate the ability of b/y to represent the land cover conditions. As we discussed above, tall trees may have smaller b/y (smaller leaf biomass contributions to total aboveground VWC) compared to short grasses, although there might be great diversity of the VWC distribution in the body of plants which are classified even in the same land cover types. We assumed that the b/y and σ is non-stationary. We performed the linear regression year by year to obtain the trend in the b/y index. We assumed there are no intra-annual changes in the b/y index. After evaluating the performance of our b/y index as the proxy of land cover, we will discuss the long-term global change in terrestrial ecosystems using the b/y index as well as the VCF data. All analyses were performed using Python. See also Figure 1 for the procedure of this study.

Results
Climatologic b/y from 2003 to 2019 is shown in Figure 2. We use all 8-day data of microwave VOD and optical LAI in the study period to perform the linear regression. The estimated slope is shown in Figure 2. Although the linear relationship between microwave VOD and optical LAI is statistically insignificant in the rainforest regions due to there being no significant seasonal cycle of VOD and/or LAI in these regions, we can retrieve b/y in most of the vegetated pixels. The existence of the linear relationship between microwave VOD and optical LAI in many land pixels is consistent with the previous findings (e.g., [14,28]). Generally, b/y tends to be low in the temperate and boreal forest areas. In semiarid grassland b/y is relatively high, indicating the large contribution of LAI to VWC. The spatial distribution of the retrieved σ-parameter shown in Figure S1 in the Supplementary Materials is similar to the previous global estimation of a roughness parameter in a radiative transfer model [22]. To confirm that b/y can be used as a proxy of the fractions of tree canopy and short vegetation, we estimated yearly b/y and compared it with the yearly land cover type data. We obtained the b/y index each year and compared it with the land cover map of the respective years. Figure 3 and Figure S2 indicate that we obtain low b/y in broadleaf and needleleaf forest areas; b/y increases in savannas, cropland, and grassland areas. In pixels To confirm that b/y can be used as a proxy of the fractions of tree canopy and short vegetation, we estimated yearly b/y and compared it with the yearly land cover type data. We obtained the b/y index each year and compared it with the land cover map of the respective years. Figure 3 and Figure S2 indicate that we obtain low b/y in broadleaf and needleleaf forest areas; b/y increases in savannas, cropland, and grassland areas. In pixels with lower (higher) b/y tree canopy (short vegetation) is dominant, so that our proposed b/y can be a proxy of fractions of tree canopy and short vegetation. We performed the Kolmogorov-Smirnov test to evaluate if the distribution of b/y in one land cover type can be distinguished from that in another land cover type. Table S1 in the Supplementary Materials shows that our b/y index is useful for classifying the satellite pixels into tree canopy, short vegetation, and bare ground. It is not straightforward to distinguish tree canopy dominant pixels from short vegetation dominant pixels by simply analyzing yearly optical LAI (Figure S3), although the clumping factors are estimated based on land cover types in the LAI retrieval algorithm. Boxplots in Figure 3 show the large variance of b/y in the specific land cover types, probably due to the coarse resolution of microwave VOD, which prevents the direct matchup of b/y with the land cover data (see also Section 2). In addition, our b/y index only quantifies the contribution of leaf biomass to the total aboveground VWC, which cannot be completely correlated to the land cover types. The large variance of b/y implies the diversity of the VWC distribution in the body of plants classified in the same land cover types. The characteristics of yearly b/y shown in Figure 3 could also be found when we compared it with the VCF data. Figure 4a reveals that the b/y index decreases as the percent cover of tree canopy increases. Yearly b/y and the VCF's tree canopy per cent cover are negatively correlated (R = −0.52). Figure 4b reveals that the b/y index increases as the percent cover of bare soil increases. Yearly b/y and the VCF's bare soil percent cover are positively correlated (R = 0.67). As we found in Figure 3, there is the large variance of b/y in the specific percent cover. First, a one-to-one relationship between b/y and VCF cover can-  Figure 4b reveals that the b/y index increases as the percent cover of bare soil increases. Yearly b/y and the VCF's bare soil percent cover are positively correlated (R = 0.67). As we found in Figure 3, there is the large variance of b/y in the specific percent cover. First, a one-to-one relationship between b/y and VCF cover cannot be expected since VCF has three-dimensional variables in each pixel. For example, the large variance of b/y with low percent cover of tree canopy can be found in Figure 4. This is because there are many combinations of short vegetation and bare soil percent covers with low percent cover of tree canopy, which induces the wide variety of the b/y index. Second, as we discussed above, our b/y index quantifies the contribution of leaf biomass to the total aboveground VWC, which cannot be completely correlated to the percent cover of tree canopy, short vegetation, and bare soil.  We can detect global land change using our yearly b/y. The linear trend of yearly b/y from 2003 to 2019 is shown in Figure 5. The decrease of b/y (shown in blue) indicates the reduction of the contribution of leaf biomass to the total aboveground VWC, which corresponds to the transition from short vegetation to tree canopy or the transition from bare ground to short vegetation (see also Figures 3 and 4). Most of the statistically significant trends shown in Figure 5 are consistent with previous findings. We can detect the increase of tree canopy in the Sahel region, which was also apparent in the data of [3]. Previous studies attributed this change to the increase of rainfall [32,33]. The increase of tree canopy in eastern Europe detected in our dataset is consistent with previous studies (e.g., [34,35]). We can detect the widespread increase of tree canopy in China, which resulted from reforestation and afforestation programs (e.g., [36]). The decreasing trend of b/y is also detected in western India, which can be attributed to the transition from bare ground to short vegetation due to agricultural activities [37]. The decrease of tree canopy in the southeastern Amazon (northeastern Brazil) detected in our dataset is also consistent with the previous findings [35], although higher spatial resolution data [35] showed much more widespread forest losses. We can detect global land change using our yearly b/y. The linear trend of yearly b/y from 2003 to 2019 is shown in Figure 5. The decrease of b/y (shown in blue) indicates the reduction of the contribution of leaf biomass to the total aboveground VWC, which corresponds to the transition from short vegetation to tree canopy or the transition from bare ground to short vegetation (see also Figures 3 and 4). Most of the statistically significant trends shown in Figure 5 are consistent with previous findings. We can detect the increase of tree canopy in the Sahel region, which was also apparent in the data of [3]. Previous studies attributed this change to the increase of rainfall [32,33]. The increase of tree canopy in eastern Europe detected in our dataset is consistent with previous studies (e.g., [34,35]). We can detect the widespread increase of tree canopy in China, which resulted from reforestation and afforestation programs (e.g., [36]). The decreasing trend of b/y is also detected in western India, which can be attributed to the transition from bare ground to short vegetation due to agricultural activities [37]. The decrease of tree canopy in the southeastern Amazon (northeastern Brazil) detected in our dataset is also consistent with the previous findings [35], although higher spatial resolution data [35] showed much more widespread forest losses.
The global land change detected by our newly developed b/y index is consistent with the state-of-the-art VCF dataset provided by [3] (Figure S4). Song et al. (2018) [3] estimated the long-term trend of the fractions of tree canopy, short vegetation, and bare ground from NDVI by supervised machine learning with high-resolution reference data. Our proposed method can detect a similar trend in the terrestrial ecosystem to [3] by the simple linear regression between microwave VOD and optical LAI, without any reference data. However, there are some differences between the trends detected by this study and the VCF data. The VCF detected short vegetation gain in the Congo rainforest, which we did not detect.
The VCF detected tree canopy gain in the eastern U.S., which cannot be detected by our b/y data. The significant tree canopy gain in Siberia, Scandinavia, and Canada detected by the VCF cannot be found in our b/y dataset. There are several reasons for these differences. First, the linear regression between microwave VOD and optical LAI is insignificant in densely vegetated areas, as shown in Figure 2, so that our b/y index cannot accurately quantify vegetation dynamics in the dense rainforest. Second, the spatial resolution of our b/y data (0.25 • ) is much lower than that of the VCF data (0.05 • ), so many impacts of human activities on vegetation dynamics in small scales may be invisible in our data. Third, it is difficult to accurately retrieve b/y in cold regions whose snow season is long. We could not retrieve VOD with snow cover, while the optical NDVI based VCF data include vegetation observation in snowy seasons. Note that the change in our b/y index is caused by the changes in both microwave VOD and optical LAI. Figure S4 also indicates that the linear trend of b/y is identical to neither LAI nor VOD, so that the fusion of microwave VOD and optical LAI provides unique information on global land change. We can detect global land change using our yearly b/y. The linear trend of yearly b/y from 2003 to 2019 is shown in Figure 5. The decrease of b/y (shown in blue) indicates the reduction of the contribution of leaf biomass to the total aboveground VWC, which corresponds to the transition from short vegetation to tree canopy or the transition from bare ground to short vegetation (see also Figures 3 and 4). Most of the statistically significant trends shown in Figure 5 are consistent with previous findings. We can detect the increase of tree canopy in the Sahel region, which was also apparent in the data of [3]. Previous studies attributed this change to the increase of rainfall [32,33]. The increase of tree canopy in eastern Europe detected in our dataset is consistent with previous studies (e.g., [34,35]). We can detect the widespread increase of tree canopy in China, which resulted from reforestation and afforestation programs (e.g., [36]). The decreasing trend of b/y is also detected in western India, which can be attributed to the transition from bare ground to short vegetation due to agricultural activities [37]. The decrease of tree canopy in the southeastern Amazon (northeastern Brazil) detected in our dataset is also consistent with the previous findings [35], although higher spatial resolution data [35] showed much more widespread forest losses.  To deepen the understanding of the long-term trend of b/y, we compare b/y with the other optical satellite products (i.e., LAI and VCF) in some specific regions ( Figure 6). It should be noted that all linear trends without performing statistical tests are shown in Figure 6 to clearly show how each index works; see Figure S4 for statistically significant trends. In the Sahel region, the VCF indicates the transition from short vegetation to tree canopy and from bare ground to short vegetation (Figure 6b,d), which is consistent with the decrease of b/y (Figure 6a). It is not straightforward to detect this land change by simply analyzing the linear trend of LAI (Figure 6c). The VCF reveals the significant transition from bare ground to short vegetation in the western part of India (Figure 6f,h). The linear trend of b/y in the western part of India is consistent with this land change (Figure 6e), while LAI increases in whole India (Figure 6g). In Europe, the VCF reveals the widespread increase of tree canopy and decrease of short vegetation (Figure 6j,l), which is consistent with the decrease of b/y (Figure 6i). The trends of b/y and LAI are not identical (Figure 6i,k), as we found in the other regions. The linear trend of VCF indicates widespread tree canopy gain with decrease in the fraction of short vegetation in central and southern China (Figure 6n,p). This land change can be detected by our method as the decrease The long-term trends of the VCF data and our b/y index are consistent (Figure 7). In Figure 7, the vertical and horizontal axis shows the trends in vegetation fractions quantified by the VCF dataset, and the color of the dots shows the trend of our b/y index. Each dot has the information on both VCF and b/y long-term trends from 2003 to 2016. For example, dots in the upper left corner of Figure 7a show changes from short vegetation to tree canopy in the VCF datasets. Most dots in this region shows the decrease of b/y. Dots in the lower right corner of Figure 7a show changes from tree canopy to short vegetation in the VCF datasets, which is consistent with the increase of b/y shown in the color of the dots. Many dots in the upper right of Figure 7a show the increase in the fraction of short vegetation with little change in the fraction of tree canopy, which implies the decrease of bare ground. In this case, our b/y index decreases in most pixels, which is also consistent with our theoretical hypothesis that the decrease of b/y corresponds to the transition from short vegetation to tree canopy or the transition from bare ground to short vegetation (see also Figures 3 and 4) The similar consistency between the VCF data and our b/y index could also be found when we sorted the pixels by the other two of the three fractions (Figure 7b,c). The long-term trends of the VCF data and our b/y index are consistent (Figure 7). In Figure 7, the vertical and horizontal axis shows the trends in vegetation fractions quantified by the VCF dataset, and the color of the dots shows the trend of our b/y index. Each dot has the information on both VCF and b/y long-term trends from 2003 to 2016. For example, dots in the upper left corner of Figure 7a show changes from short vegetation to tree canopy in the VCF datasets. Most dots in this region shows the decrease of b/y. Dots in the lower right corner of Figure 7a show changes from tree canopy to short vegetation in the VCF datasets, which is consistent with the increase of b/y shown in the color of the dots. Many dots in the upper right of Figure 7a show the increase in the fraction of short vegetation with little change in the fraction of tree canopy, which implies the decrease of bare ground. In this case, our b/y index decreases in most pixels, which is also consistent with our theoretical hypothesis that the decrease of b/y corresponds to the transition from short vegetation to tree canopy or the transition from bare ground to short vegetation (see also Figures 3 and 4) The similar consistency between the VCF data and our b/y index could also be found when we sorted the pixels by the other two of the three fractions (Figure 7b,c).
To evaluate the robustness of our b/y retrieval, we used the VOD Climate Archive (VODCA) [28] dataset and performed the linear regression shown in Equation (3). Here we discuss the effect of the uncertainty in the VOD data on our newly developed b/y index. Figure S5 reveals that the VODCA-based b/y index has the same characteristics as the original b/y index: in pixels with lower (higher) VODCA-based b/y, tree canopy (short vegetation) is dominant. Even if we replace the VOD dataset, the important characteristics of our b/y index hold. However, the absolute value of the VODCA-based b/y index tends to be smaller than the original b/y index (Figures 3 and S5). The long-term trend of the VODCA-based b/y index also tends to be smaller than that of the original b/y index (Figures 5 and S6). Although the retrieval algorithms of these two datasets are similar, many different satellite data were blended by the cumulative distribution function matching to generate the long-term VODCA data. It probably reduces the dynamic range of VOD, which makes the magnitude and trend of b/y smaller than the original b/y index. This result implies that the selection of VOD data affects the retrieval of the b/y index, although the characteristics of our b/y index are qualitatively unchanged with the different VOD datasets. To evaluate the robustness of our b/y retrieval, we used the VOD Climate Archive (VODCA) [28] dataset and performed the linear regression shown in Equation (3). Here we discuss the effect of the uncertainty in the VOD data on our newly developed b/y index. Figure S5 reveals that the VODCA-based b/y index has the same characteristics as the original b/y index: in pixels with lower (higher) VODCA-based b/y, tree canopy (short vegetation) is dominant. Even if we replace the VOD dataset, the important characteristics of our b/y index hold. However, the absolute value of the VODCA-based b/y index tends to be smaller than the original b/y index (Figure 3 and Figure S5). The long-term trend of the VODCA-based b/y index also tends to be smaller than that of the original b/y index ( Figure  5 and Figure S6). Although the retrieval algorithms of these two datasets are similar, many different satellite data were blended by the cumulative distribution function matching to generate the long-term VODCA data. It probably reduces the dynamic range of VOD, which makes the magnitude and trend of b/y smaller than the original b/y index. This result implies that the selection of VOD data affects the retrieval of the b/y index, although the characteristics of our b/y index are qualitatively unchanged with the different VOD datasets.

Implications to the Existing VOD Retrieval Algorithms
The results shown above imply that the empirical equations and assumptions of the MiOFA described in Section 2.2 are accurate on a global scale. Since these empirical equations and assumptions are closely related to the algorithms for retrieving surface soil moisture and VOD from brightness temperature (e.g., [5,13,19,38]) and for assimilating brightness temperature into land surface models (e.g., [39][40][41][42]), our retrieved b/y and parameters may greatly contribute to improving these algorithms and the understanding of microwave radiative transfer in the canopy. The most important implication of this study is that the relationship between microwave VOD (or VWC) and the optical vegetation index is non-stationary, although it is assumed to be stationary in previous studies on the development and validation of the VWC retrieval algorithms (e.g., [27]). In many cases, the

Implications to the Existing VOD Retrieval Algorithms
The results shown above imply that the empirical equations and assumptions of the MiOFA described in Section 2.2 are accurate on a global scale. Since these empirical equations and assumptions are closely related to the algorithms for retrieving surface soil moisture and VOD from brightness temperature (e.g., [5,13,19,38]) and for assimilating brightness temperature into land surface models (e.g., [39][40][41][42]), our retrieved b/y and σ parameters may greatly contribute to improving these algorithms and the understanding of microwave radiative transfer in the canopy. The most important implication of this study is that the relationship between microwave VOD (or VWC) and the optical vegetation index is non-stationary, although it is assumed to be stationary in previous studies on the development and validation of the VWC retrieval algorithms (e.g., [27]). In many cases, the newly developed VOD products were evaluated by the correlation coefficient between VOD and optical vegetation indices (e.g., [6,14,28]), assuming the stationary relationship between them. The evaluation and interpretation of VOD datasets can be improved by explicitly considering the non-stationary relationship between VOD and optical vegetation indices. The trend in the VOD-LAI relationship is an important signal for detection of global environmental change, so it should not be neglected in the development, validation, and analysis of satellite land observation products.

Limitations
There are several limitations in this study. First, the C-band VOD is affected by radio frequency interference (RFI) in many parts of the world, such as the United States and the Middle East (see Figure S7 of [28]). Since we cannot use microwave observations with higher frequency to perform MiOFA (see Section 2.2), our future work will focus on the long-term L-band VOD data to avoid the effect of RFI.
Second, optical remote sensing is strongly affected by clouds, while microwave remote sensing has all weather capability. The influence of clouds negatively impacts the quality of the optical LAI product as well as our proposed b/y index, especially in rainforest areas.
Third, the C-band VOD may have limited penetration depth of the canopy to observe the whole structure of densely vegetated rainforests, although many previous works successfully monitored Amazon and Congo rainforests (e.g., [7,8]). This potential limitation of our method will also be addressed by the future development of long-term L-band VOD data.
Fourth, the spatial resolution of our b/y data is too coarse to be directly evaluated by the existing ground truth datasets and to be applied to assess the local-scale environments. In addition, combining the two datasets (i.e., VOD and LAI datasets) with significantly different spatial resolutions may bring additional uncertainty into our b/y data.
Fifth, further field-supported verification is necessary. Our newly developed b/y index, which is the proxy of the contribution of leaf biomass to total aboveground VWC, is so unique that it is difficult to find the de facto satellite observation products useful to validate our index. Note that the information in the VCF data do not fully correspond to what we observed with the b/y index, although we found our b/y index behaves similarly to the VCF data. It is currently impossible to get direct observation of the contribution of leaf biomass to the total aboveground VWC in the satellite footprint scales. Since the destructive sampling of vegetation is needed to directly observe it, an intensive field observation campaign is needed to deepen our understanding of the newly developed b/y index.

Conclusions
We propose a new vegetation index by combining microwave VOD and optical LAI. Our new index can provide a good proxy of the respective fractions of tree canopy and short vegetation, and is useful for detecting global change in the terrestrial ecosystem. The advantage of our new method is that the vegetation index can be obtained by the simple linear regression between microwave VOD and optical LAI, so that no high spatial resolution reference data are required. Therefore, our proposed method makes it possible to easily quantify the condition of terrestrial ecosystems in regions where no reliable ground truth data can be obtained. The limitation is the lower spatial resolution compared to other existing methods. We show that our b/y index is consistent with the existing land cover datasets which were thoroughly evaluated by previous studies. However, a crucial need is the further evaluation of our b/y index by field observation of the contribution of leaf biomass to the total aboveground VWC in the satellite footprint scale. Future work will focus on the extensive verification of the new vegetation index by reference to field data. We will also focus on the applications of this new vegetation index in deepening the understanding of vegetation dynamics in extreme climate conditions such as drought.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13183756/s1, Text S1: Bias correction of VOD from AMSR-E and AMSR2, Figure S1: Climatologic σ [-] (see Section 3). Pixels which show a statistically significant slope between microwave VOD and optical LAI (student's t-test, p < 0.05) are depicted, Figure S2: Histograms of yearly b/y [-] (see Section 3) stratified by land cover types, Figure S3: Same as Figure S2 but for LAI, Figure S4: Linear trend of (a) b/y, (b) fraction of tree canopy, (c) LAI, (d) fraction of short vegetation, (e) VOD, and (f) fraction of bare ground from 2003 to 2016. For (b-f), the data generated by [3] were used. Pixels which show a statistically significant trend (stu-dent's t-test, p < 0.05) are depicted, Figure S5: Same as Figure 2 in the main manuscript but for the VODCA dataset, Figure S6: Same as Figure 4 in the main manuscript but for the VODCA dataset, Table S1: p-values of the Kolmogorov-Smirnov test for b/y values (in 2013) of the two paired land use types. Funding: This study was supported by JAXA grant ER2GWF102 and JSPS KAKENHI grant 18H03800 and 21H01430.