## 7.2.2 Pro and cons of the different approaches

### 7.2.2.1 Computation requirements

Machine learning approaches, once calibrated, are obviously very little demanding in terms of computation. The inverse model is generally relatively simple and could be run very quickly. However, the calibration (or learning or training) process could require large computer resources, particularly for complex parametric model with a significant number of coefficients to be tuned and when the training dataset is large. The implementation of a LUT technique in algorithmic operational chains is very efficient because the radiative transfer model is run off-line. Conversely, iterative minimization methods require large computer resources because of its iterative nature. Improvements are possible using a meta-model. Further, automatic segmentation or discretization of the reflectance space (Pinty et al., 2011) will also reduce the number of inversions to be completed over a whole set of images.

### 7.2.2.2 Flexibility of the observational configuration

Iterative optimization methods allow retrieving canopy characteristics from several observational configurations. It is even possible to invert radiative transfer models concurrently over several pixels. This opens great potentials for exploiting additional temporal or spatial constraints Atzberger (2004); Lauvernet et al. (2008). LUT could theoretically cope with variable configurations at the expense of the dimensionality and thus the size of the tables, making them more difficult to manipulate. Conversely, machine learning methods require a fixed number of inputs. The characteristics of the configuration need thus to be used as inputs of the inverse parametric model, where the illumination and view directions are explicitly used. However, this increases the dimensionality of the system, making the calibration step more demanding and more difficult. One alternative is to calibrate several parametric models for each individual configuration and then select the proper calibrated inverse model.

### 7.2.2.3 Integration of prior information

The radiometric data driven approaches integrate the prior information directly in the cost function using a regularization term Combal et al. (2002). However, in the case of LUTs, it is also possible to restrict the simulations to the range of situations to be encountered as is done within the MODIS LAI and FAPAR algorithm that depends on the biome type considered Shabanov et al. (2005). For the machine learning approaches, the prior information is introduced through the distributions and co-distributions of the inputs of the radiative transfer model: when vegetation traits have to be estimated under situations where the type of canopies and their stage of development are known, it is more efficient to calibrate a specific inverse model for each individual situation. Note that Qu et al. (2008) proposed to use Bayesian networks where model simulations could be exploited along with a description of the distribution of the variables that may depend on growth stages or canopy types.

### 7.2.2.4 Associated uncertainties

The radiometric data driven approaches allow getting some estimates of the uncertainties associated to the solution by propagating the uncertainties associated to the measurements and to the model using the partial derivatives of the cost function with regards to the measurements (Lauvernet 2005). When using LUTs, uncertainties could be estimated by Monte-Carlo methods or approximated by the standard deviation of the ensemble of solutions defined by the uncertainties in the measurements (Knyazikhin et al. 1998). For machine learning methods, the error on the measurements may be assessed in different ways as proposed by Aires et al. (2004), which may also include the errors associated to the retrieval process itself. A more simple alternative solution is also proposed by Baret et al. (2013) based on the training dataset. Although the estimation of uncertainties of the retrievals is possible, it is generally limited by the poor knowledge on the input uncertainties associated to the reflectance measurements and radiative transfer models used. Knyazikhin et al. (1998) used a 20% relative uncertainty applied to MODIS top of canopy reflectance for LAI and FAPAR retrieval. Baret
et al. (2007) propose to use an additive uncertainty around 0.05 and a multiplicative uncertainty around 3%. This example shows that the uncertainties attached to each band is poorly known. Further, the structure of the uncertainties may also play an important role and is unfortunately very difficult to describe.

### 7.2.2.5 Robustness of the retrieval and quality assessment

A quality index needs to be associated to the retrieved values to inform about the status of the inversion process. For iterative optimization techniques, it could be criterions relative to the stop of the iterations (Gilbert 2002). As a matter of facts, the algorithm may sometimes encounter numerical problems occurring generally with very small values of J. No numerical problems are expected for LUT and machine learning approaches, and the quality index should mainly indicate whether the input reflectances were inside the definition domain and if the output solution is in the expected range of variation Baret et al. (2013). The performances of the approach will both depend on the minimization algorithm itself and on the level of illposedness of the inverse problem as a function of measurement configuration and model and measurement uncertainties.

### 7.2.3 Using VIs or reflectance for biophysical parameter retrieval?

The literature on the subject started about 30 years ago focusing on a restricted number of VIs available. At that time, VIs were the main approach used to characterize the vegetation. With the multiplication of satellite systems and the progress of our understanding of the physical principles governing the radiative transfer in vegetation canopies, other approaches were proposed as described in the previous section. VIs are combinations of reflectance in several bands, generally limited to 2 or three bands. They therefore exploit only the spectral dimension of the observation. They are using the differential absorption properties of the vegetation and background materials that govern the spectral variation of canopy reflectance. However, if canopy reflectance depends on the vegetation traits of interest, several other factors should be considered like soil background, atmosphere with main effects for the shorter wavelengths, sun/sensor geometry with the bidirectional effects, topography effects, scaling effects when the pixel is composed of patches containing contrasted objects, canopy architecture when the targeted trait is the biochemical composition or the biochemical composition when the targeted trait is related to canopy architecture.

VIs are therefore designed empirically with the objective to maximize the sensitivity to the given vegetation trait targeted, while minimizing the impact of confounding factors such as the soil background, atmospheric effects, topography effects, or the effect of other vegetation characteristics.

VIs have been used since the 1970's. At that time, even if the physical principles of the radiative transfer in vegetation canopies was well understood, there was no mature models and the required computer resources to use them to calibrate retrieval algorithms. VIs correspond to a parsimonious way to combine few spectral bands into a single quantity to be related to a vegetation variable of interest. Simple statistical models can be therefore calibrated over a small training dataset. For coarse resolution satellites for which the collection of ground data is difficult, VIs offered a convenient way to build relatively consistent long time series of vegetation observations.

The advancement of the radiative transfer models, the availability of high-quality and high spatial and temporal resolution satellites, the development of machine learning techniques, including deep-learning allows nowadays to develop more sophisticated approaches that will fully and optimally exploit the spectral information available to get more robust and accurate estimate of the few main vegetation traits. The reader is referred to WP6 report, section 4, Baret and Weiss (2024) for a full description of the development of VIs, how they are related and their sensitivity to vegetation traits, and how to mitigate with confounding effects to better assess biophysical variables from optical measurements.
# 7.3 Designing the algorithm for CWC retrieval from optical measurements

This study has the objective of proposing global and consistent CWC products at both decametric (Landsat-8, Sentinel-2) and 500m (MODIS) scale. For this purpose, we relied on former developments for operational algorithms Bacour et al. (2006) ; Baret et al. 2007b), including the one developed for the S2 SNAP toolbox Weiss et al. (2023) and developed a unified approach based on machine learning trained with RTM simulations, with particular attention to the following aspects: (1) Optimizing the representativeness of the training database by combining the knowledge extracted from two open databases (e.g. the newest version of TRY for plant traits and OSSL for soil background); (2) Proposing an appropriate validation scheme by firstly assessing CWC at decametric resolution that are much more compatible with the ground measurement footprint using ground data covering diverse biome types over the world, and then comprehensively evaluating the MODIS CWC using Landsat8 and Sentinel-2 CWC at the global scale; (3) Testing the effectiveness of optical indices and different band combinations for optimizing CWC retrievals more widely, and comparing the performances of CWC retrievals with those of VOD retrieved from the passive microwave sensors. The reader may also refer to the paper by Ma et al, submitted to Remote Sensing of Environment, August 2024.

## 7.3.1 Material

Here, we present the satellite optical data that are used as inputs to the CWC retrieval algorithm, then the databases employed to build a realistic training dataset using the PROSAIL radiative transfer model, and finally the site network that we used to evaluate our algorithm, including ground validation sites and those used for the inter-comparison exercise between different satellite products. In this way, we followed the recommendations of the CEOS (Committee on Earth Observation Satellite) LPV (Land Product Validation) group for vegetation biophysical variables.

### 7.3.1.1 Optical satellite data

In this study, two high resolution satellites including Sentinel-2 and Landsat-8 (20m-30m), as well as the moderate resolution satellite MODIS (1km) are used. For each sensor, we used surface reflectance L2A data at the time of acquisition date for Landsat-8 (https://www.usgs.gov/) and Sentinel-2A and 2B (https://scihub.copernicus.eu/). For MODIS, we considered the MCD43A4 Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR, V6.1) product at 500m spatial resolution. This product results from a kernel driven model normalization in the nadir direction over a 16-day period centered at the 9th day. For Sentinel-2A and 2B, the 10m green (B3) and red (B4) bands were resampled at 20m through spatial averaging, to be consistent with the spatial resolution of the other considered bands. We used a single NIR band for Sentinel-2, e.g., B8a which has a closer spectral range to Landsat-8 than B8 to keep better consistency between the sensors. It can be also noted that, in the SWIR domain (e.g. CWC sensitivity bands), MODIS has an additional band (i.e., SWIR3, 1230nm-1250nm) compared to Sentinel-2 and Landsat8. Additionally, the geometry information including sun zenith angle, view zenith angle and relative azimuth angle (i.e., solar azimuth minus view azimuth) were extracted for the three satellites. The three sensor surface reflectance products are derived with different algorithms due to both the inherent characteristics of the sensors (e.g., presence of a thermal band for the detection of clouds for Landsat-8 and MODIS) and choice made for the processing (e.g., atmospheric correction scheme, cloud screening, 16 days synthesis for MODIS). The largest differences are related to the cloud screening since different cloud definitions are used (e.g. thin cirrus, Skakun et al. (2022). Additionally, the compositing period and the coarse spatial resolution of MODIS may affect the surface reflectance product when small clouds occur within the 500m pixel Roy et al. (2008). To avoid such situations and inconsistencies between products, we filtered the data by considering that NDVI values should be consistent between the three sensors at the same acquisition date for a same location. For the same reasons as for the cloud influence, we also used the MODIS snow mask.
### 7.3.1.2 TRY and OSSL databases

As the biggest global effort to comprise a massive repository of global plant traits, TRY is capable of representing the actual leaf traits variations at an unprecedented spatial, temporal and climatic coverage KATTGE et al. (2011), Kattge et al. (2020). Since its establishment in 2007, TRY has grown continuously, thus improving representativeness for the terrestrial leaf traits, and we used the data available up to April, 2023. We extracted three key related leaf traits (leaf chlorophyll content, Cab, leaf dry matter content, Cdm, and LWC) that were used to determine the distribution of the PROSPECT model inputs after removing potential outliers by keeping only the data between the first and 99th percentile to 99% for each trait. With regard to the soil spectral data, the CWC retrieval algorithm relies on former developments, which used a soil spectral database collected by INRAE Avignon over different sites distributed all over France Weidong et al. (2002). However, some weaknesses were observed by validation studies of the SNAP algorithm for low vegetation areas Vinué et al. (2018), Estévez et al. (2020), which indicated that this database is not sufficiently representative of the diversity of the soil at the Earth surface. We therefore combined it with OSSL, another open access and massive database for soil spectral properties at the global scale Hengl et al. (2021). We resampled the OSSL soil spectra using a cubic convolution interpolation approach into 1nm spectral resolution, before combining it with the INRAE database. The resulting dataset of 35835 soil spectral samples was then exploited to improve the representativeness of the soil spectral background in our training database.

### 7.3.1.3 Ground CWC data

We evaluated the decametric CWC products using ground destructive measurements from five ground sites representing 3 biome types (2 grasslands, 2 forests and 1 cropland), distributed in diverse continents (Table 7.1, 7.1) by matching the Sentinel-2 and Landsat-8 acquisition dates with the date of measurements (+/- 2 days) during which we assume CWC exhibits very slight variations at the satellite overpass time. For the sites of Dahra, SMAPVEX16 and EucFACE, both Sentinel-2 and Landsat-8 could be used, while only Landsat-8 was present for VALSE2 and Qinghai with earlier observation periods. Regarding the spatial matching, the Dahra site can be considered as a homogeneous grassland, and the 28 CWC measurements of 1m² plot were averaged Tagesson et al. (2015); for EucFACE, the ground CWC was measured in rings of 25m diameter Quan et al. (2017) Griebel et al. (2023); while for Qinghai sites, the ground CWC was acquired at each 30m*30m plot. The size of sample range for these two sites are close to one pixel of Sentinel-2 (20m) or Landsat-8 (30m). In terms of temporal matching, only the decametric optical satellite data were extracted at the dates of the ground measurements (+/- 2 days). Here it should be noted that for SMAPVEX16, ground CWC of diverse crop types including soybean, wheat, canola, oat, and corn etc. were collected. As also pointed out in existing studies by Djamai et al. (2019) and Garcia-Haro et al. (2020), the ground data in this area were calculated by the difference of wet and dry weights of leaves (for soybean, canola and corn) or of combinations of leaves and stems (for wheat and oat). Moreover, the data for canola during the crop flowering stage (DOY from 180 to 203) were not considered for the CWC assessment, as the flowers may be absent in PROSAIL simulations adopted in this study. Unfortunately, we could not access additional sites or data via destructive approaches to complement our validation exercise at this spatial scale, as we tried to make our best efforts to guarantee the representativeness and broadness of the CWC assessment in this study.

<table><caption>Table 7.1: Ground CWC sites</caption><thead><tr><td>Ground Sites</td><td>Continent/Country</td><td>Time range used</td><td>Matching satellite</td><td>Landcover</td><td>Reference</td></tr></thead><tbody><tr><td>Dahra</td><td>Africa/Senegal</td><td>2013 July-2019 October</td><td>Landsat-8/Sentinel2</td><td>Grassland</td><td>Tagesson et al. 2015</td></tr><tr><td>SMAPVEX16</td><td>America/Canada</td><td>2016 June-2016 July</td><td>Landsat-8/Sentinel2</td><td>Cropland</td><td>McNairn et al. 2017</td></tr><tr><td>VALSE2</td><td>Europe/France</td><td>2013 June-2013 September</td><td>Landsat-8</td><td>Forest</td><td>Weiss and Baret 2013</td></tr><tr><td>EucFACE</td><td>Oceania/Australia</td><td>2019 September-2021 June</td><td>Landsat-8/Sentinel2</td><td>Forest</td><td>Griebel et al. 2023</td></tr><tr><td>Qinghai</td><td>Asia/China</td><td>2014 July-2015 August</td><td>Landsat-8</td><td>Grassland</td><td>Quan et al. 2017</td></tr></tbody></table>
<figure><img src="image_53.png" /><figcaption>Figure 7.2: Distribution of the LANDVALSUB and ground validation sites used in this study. The 500m resolution IGBP map is shown as biome type description (Friedl et al. 2002). For the specific presentations of each ground sites, the ESA WorldCover 10m 2022 was adopted in this study (Zanaga et al. 2022). DBF: Deciduous Broadleaf Forest; WS: Woody Savannas; ENF: Evergreen Needleleaf Forest; MF: Mixed Forests; EBF: Evergreen Broadleaf Forest; CS: Closed Shrublands; DNF: Deciduous Needleleaf Forest; OS: Open Shrublands. Selection of LANDVAL sites</figcaption></figure>

### 7.3.1.4 Selection of LANDVAL sites

The satellite product inter-comparison exercise uses the LANDVAL network of 725 sites which was specifically developed for benchmarking vegetation products such as LAI, fAPAR, or albedo Fuster et al. (2020). LANDVAL is a compilation of different site networks that well represent the combination of climate, vegetation types and environment at the Earth Surface. Additionally, these 3kmx3km sites were selected to be homogeneous in terms of biome type, not too close to a water body, flat, and do not correspond to urban areas. In this study, to reduce the extraction and computation requirements as well as guarantee the representativeness of diverse biome types in the world, we extracted a selection of 308 sites that we named LANDVALSUB sites Figure 7.2. Firstly, we divided the Earth into rectangles of 20° in latitude (from -60° to 80°) and 40° in longitude (from -180° to 180°). We identified the land covers of all the LANDVAL sites present in each rectangle, and then randomly selected half of them for each biome class, except for the desert and bare soil sites. Moreover, to further check the variations of satellite CWC in low water conditions, we additionally selected 7 sites for bare soil or deserts (i.e., CWC close to zero). For each site, we then spatially averaged the CWC products at their original resolution over the 3km by 3km area. Additionally, as the MODIS product results from a synthesis over a 16-day temporal window, and considering the low revisit frequency of Landsat (16 days), we considered a 5-days tolerance for the benchmarking of the optical CWC products over the LANDVALSUB sites.
## 7.3.2 Method

We developed a generic and consistent algorithm to derive an operational CWC from any optical sensor. This algorithm relies on former findings related to the use of machine learning trained on radiative transfer simulations to operationally derive different vegetation biophysical variables products from medium resolution (e.g., hectometric, and kilometric) or decametric sensors We developed a generic and consistent algorithm to derive an operational CWC from any optical sensor. This algorithm relies on former findings related to the use of machine learning trained on radiative transfer simulations to operationally derive different vegetation biophysical variables products from medium resolution (e.g., hectometric, and kilometric) or decametric sensors Weiss et al. (2002), Baret et al. (2007), Verger et al. (2009), Brown et al. (2021). We set up the training database of the machine learning algorithm through the PROSAIL radiative transfer model simulations run with a priori distributions input variables partly provided by the TRY and OSSL open databases (Figure 7.3, left). Then, by considering each sensor spectral responses, and adding a realistic noise to simulate the actual satellite TOC reflectance, we tested different band and optical indices combinations as inputs to the machine learning algorithm and selected the best one evaluated on a synthetic database made of independent PROSAIL simulations (Figure 7.3, middle). Finally, we evaluated the CWC retrieval algorithm by assessing the accuracy of the Landsat-8 and Sentinel-2 products against ground measurements, inter-comparison the different optical CWC product on the LANDVALSUB site network, and comparing the spatio-temporal trends of MODIS CWC and microwave VOD products (Figure 7.3, right).

<figure><img src="image_54.png" /><figcaption>Figure 7.3: lowchart of CWC retrieval. The whole processing is separated into three main parts: (i) in green dashed frame, the parameterizing for PROSAIL model, concerning soil background, leaf traits as well as other essential variables; (ii) in blue dashed frame, the relationship building between the simulated TOC by PROSAIL and CWC; (iii) in yellow dashed frame, the assessment and discussions for the retrieved CWC from optical satellites.</figcaption></figure>

### 7.3.2.1 Setting up the PROSAIL training database

By coupling the leaf optical properties model PROSPECT and the canopy reflectance model SAIL, PROSAIL has been adopted for several vegetation parameter retrievals from optical satellite observations Jacquemoud et al. (2009) and was successfully used for operational kilometric vegetation products at the global scale, despite some inconsistencies between some biome type structure and the model turbid medium assumption
Weiss et al. (2007). As compared to more complex models (e.g. 3DRTM) and besides the fast computing aspect, it presents the main advantage of having a small number of inputs which greatly eases the setting up of a training database representing the distribution of the vegetation surface reflectance at the Earth surface. Moreover, most of PROSAIL input parameters have been extensively studied and correspond to measurable variables. Additionally, CWC can be easily computed from PROSAIL inputs as the product of the green leaf area index (LAI or Green Area Index, GAI, cm2/cm2) and the green leaf water content per green area (Cw, g/cm2). We therefore took a particular attention to the distribution of PROSAIL inputs to improve the representativeness and the realism of the training database. As suggested by Garcia-Haro et al. (2020), we used the latest version of the TRY database to simulate 3 PROSPECT input variables: leaf chlorophyll content, leaf dry matter content, and leaf water content (LWC) instead of the relative water content (e.g. leaf water content divided by fresh biomass). Using LWC allows to keep consistency between the leaf dry matter and water content and directly simulate realistic co-distributions. Additionally, we did not consider weighting the extracted TRY distributions by PFT occurrence at the Earth surface as Garcia-Haro et al. (2020) showed little differences between both distributions. Therefore, we used directly the kernel density estimation (KDE) method to approximate the probability density function (PDF) of each leaf trait in TRY (Figure 7.3, Middle). The distribution of the other variables required in PROSAIL were kept the same as in the SNAP algorithm. The readers can refer to details in the ATBD document Weiss et al. (2023).

<figure><img src="image_55.png" /><figcaption>Figure 7.4: Distributions of PROSAIL inputs extracted from the TRY database: Cdm (g/cm2, top), LWC (unitless, middle) and Cab (µg/cm2, bottom). Left: former distribution in the SNAP used in the SNAP algorithm (green) compared to the distribution in TRY, middle: pdf distribution of each trait approximated by KDE, right: final distribution used for the training database (Green). Note V0 means the algorithm version used in SNAP.</figcaption></figure>

Similarly, we used the OSSL soil database to simulate the soil background. However, to avoid some over-representation of similar soil spectra, we first reduced its dimension by using the soil brightness concept,
considering that the reflectance of each soil of the database in a given band can be estimated as a brightness multiplied by a reference soil spectrum, selected among an ensemble of n reference soil spectra. The n soils then constituted the database of soil reference value for building the PROSAIL training database. We used an iterative approach to determine the n soils by applying a k-means classification (with k=1,2,...,n) and selecting the k soil whose spectrum was the closest from the centroid of each class. The optimized number of classes was determined by computing the RMSE between each soil spectrum and the one reconstructed using the brightness parameter and the reference value. We found that n=47 among the 35835 spectra was sufficient to well represent the range of spectral shapes within a good accuracy (RMSE around 0.01), and the related figures could be seen in Weiss et al (2023) (Weiss et al. 2023).

### 7.3.2.2 Training the machine learning algorithm

Many inter-comparison studies have been performed between different machine learning regression algorithms. While, the literature concurs on the importance of building a comprehensive training and realistic dataset, there is no clear consensus about which algorithm performs best. We therefore decided to keep the Neural Networks used in the SNAP algorithm as this is the only method implemented in operational contexts and comprehensively evaluated. However, as the neural networks convergence depends on the initialization of weights and biases that is done randomly, we considered an ensemble approach by training 12 neural networks and use the corresponding median value as the final estimate. The number of 12 was optimized over the testing datasets. This allows to regularize the inversion process and get more stable estimates over time and space, we set up the training database of the machine learning algorithm through the PROSAIL radiative transfer model simulations run with a priori distributions input variables partly provided by the TRY and OSSL open databases (Figure 4). Then, by considering each sensor spectral responses, and adding a realistic noise to simulate the actual satellite TOC reflectance, we tested different band and optical indices combinations as inputs to the machine learning algorithm and selected the best one evaluated on a synthetic database made of independent PROSAIL simulations (Figure 3, middle). For this part, the reader is referred to the paper as the theoretical performances of diverse band combinations were found very similar, whether using the TOC only or inclusion of optical indices. To keep spectral consistency between sensors, we used B4+B8a+B11+NDII+Angles for S2, B4+B5+ B6+NDII+Angles for LANDSAT8 and B1+B2+B6+NDII+Angles for MODIS. Finally, we evaluated the CWC retrieval algorithm by assessing the accuracy of the Landsat-8 and Sentinel-2 products against ground measurements, inter-comparison the different optical CWC product on the LANDVALSUB site network (Figure 3, right). For the evaluation of CWC, four error metrics including root mean square error (RMSE), unbiased RMSE (ubRMSE) or precision, Bias or accuracy, and correlation coefficient (R) were considered in this study.

## 7.3.3 Evaluation of the CWC algorithm performances

### 7.3.3.1 Direct validation of satellite CWC estimates against ground measurements

Based on the evaluation results against ground sites as shown in Figure 7.1, regarding the correlation coefficient with ground CWC, the satellite retrievals perform the best in grasslands (Dahra and Qinghai) both in terms of R and RMSE, then followed in croplands (SMAPVEX16), and forests (VALSE2 and EucFACE). As for the temporal variations of ground CWC, in spite of the same best performance in grasslands, we obtained smaller RMSE values in forests than those in croplands. In general, the satellite CWC retrievals in this study achieved acceptable performances in various vegetated conditions, which confirms the effectiveness of PROSAIL model for CWC estimations using optical data found by Djamai et al. (2019) and Garcia-Haro et al. (2020).
<figure><img src="image_56.png" /><figcaption>Figure 7.5: Satellite CWC estimates against ground measurements. The green triangle and blue snowflake represent Landsat8 and Sentinel-2A/B, respectively. (a) Dahra and (e) Qinghai are grassland sites; (b) SMAPVEX16 is cropland site; and (c) VALSE2 and (d) EucFACE are forest sites. (f) All data represents the scatter plot when putting data of all the sites together.</figcaption></figure>

### 7.3.3.2 Benchmarking Sentinel-2, Landsat-8 and MODIS CWC products

The MODIS CWC product at 500m cannot be compared to ground measurements because of its large spatial resolution that hampers ground measurements at this footprint. Therefore, we performed a comprehensive comparison over several biome types by using the LANDVALSUB sites for year 2019, by averaging all the pixels within the 3km*3km area for three satellites. Figure 7.1 and Figure 1 show good consistencies between MODIS CWC and Landsat-8 or Sentinel-2 CWC for most of the biome types. The correlation is higher than 0.8, except for bare soil area (BA, e.g., desert area) and Evergreen Broadleaf Forest (EBF). The comparison between the two decametric sensors shows a better result than MODIS (Figure 8), achieving a correlation coefficient around 0.9, even for all biomes, including the BA and EBF. This indicates the strong robustness of our algorithm between the different satellites. The differences and uncertainties that are observed may be attributed to three features: (i) the scaling effect due to the nonlinearity of the relationship between CWC and reflectance Landsat-8 (ii) the differences between the pre-processing chains of the satellites, especially regarding the cloud screening and the compositing period. This is quite well illustrated for EBF that are characterized by small clouds often mis-detected at coarse spatial resolution since the signal results from the mixing of clouds and surface reflectance. Additionally, EBF are located near the Equator and are affected by very high cloud occurrence that may impact the retrieval of Nadir reflectance from the kernel BRDF model inversion over the 16 day compositing period, while we considered the closest acquisition date when comparing Landsat-8 and Sentinel-2.
<figure><img src="image_57.png" /><figcaption>Figure 7.6: Scatter plots for CWC between high resolution Sentinel-2 and moderate resolution MODIS in diverse biome types. BA is Sparse and Bare Areas; CUL is Cultivated lands (agricultural lands); DBF is Deciduous Broadleaf Forest; EBF is Evergreen Broadleaf Forest; HER is Herbaceous (e.g., grasslands); NLF is Needle-leaf Forest; OF is Other forests, including Mixed Forest; and SHR is Shrublands.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>Figure 7.7: Scatter plots for CWC between high resolution Landsat8 and moderate resolution MODIS in diverse biome types. BA is Sparse and Bare Areas; CUL is Cultivated lands (agricultural lands); DBF is Deciduous Broadleaf Forest; EBF is Evergreen Broadleaf Forest; HER is Herbaceous (e.g., grasslands); NLF is Needle-leaf Forest; OF is Other forests, including Mixed Forest; and SHR is Shrublands.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 7.8: Scatter plots for CWC between high resolution Sentinel-2 and Landsat8 in diverse biome types. BA is Sparse and Bare Areas; CUL is Cultivated lands (agricultural lands); DBF is Deciduous Broadleaf Forest; EBF is Evergreen Broadleaf Forest; HER is Herbaceous (e.g., grasslands); NLF is Needle-leaf Forest; OF is Other forests, including Mixed Forest; and SHR is Shrublands.</figcaption></figure>

## 7.4 Conclusions

In contrast to other vegetation products that are more mature and identified as essential climate variables like LAI and fAPAR, there were only few attempts to develop algorithm for CWC estimation, and no product is available yet, with very few available ground validation measurements. However, we were able to conduct both indirect comparisons between the products we developed at two spatial resolutions, and to conduct some direct validation as recommended by the CEOS LPV group. These studies also consider assessing temporal profiles of LAI and fAPAR as the spatio-temporal variations of these variables are quite well known, smooth and continuous (apart from natural hazards or strong human intervention), and available through different vegetation functioning or land surface models as a function of climate, soil and environmental conditions. These temporal variations are much less documented for CWC, which is susceptible to vary less smoothly, depending on the density and frequency of precipitation events or water stress, related to soil conditions. However, we compared the seasonal behaviour of our product with other similar ones (e.g. VOD). More detailed is provided in the paper submitted to remote sensing of environment and these behaviours are analysed in more details within WP10. While the performances of the proposed CWC product appear satisfying for most biome
types in terms of expected spatio-temporal behavior the direct validation exercise, although limited to the
several available sites, show uncertainties that could be attributed to three main issues related to (i) the
representativeness of the training database (ii) the accuracy and footprint of the ground measurements, and
(iii) the saturation effect for dense canopies. For the representativeness, even the global efforts concerning leaf
traits of diverse biome types were delivered into PROSAIL model in this study, the possible limitations and
bias of TRY concerning diverse biome types, locations and temporal samplings, would also partially contribute
to the uncertainties. For example, with not very high representation for crop species in TRY as those for
other vegetation species Kattge et al. (2020), the uncertainties of CWC evaluations in croplands could also be
partially explained. Moreover, since this study considered all the kinds of biome together by using PROSAIL
model, in which the canopy is represented only by green components and a turbid medium. This may induce
a high simplification of biome representation. In terms of ground sites, the diverse measurement strategies
could clearly explain the different results shown in 7.5. The in situ data in Dahra, EucFACE and Qinghai, were
obtained within the area ideally matching with one Landsat8/Sentinel-2 pixel (30m/20m). In this case, less
uncertainties were taken into the CWC evaluations due to the relieved spatial mismatching between satellite
pixel and ground site, comparing to SMAPVEX16 designed for SMAP mission with coarse grids (McNairn et
al. 2017). Moreover, the environmental heterogeneity of ground data would also slightly contribute to the
uncertainties. For the SMAPVEX16 site, although fully covered with croplands, diverse crop types including
Oat, Wheat, Canola, Black bean, Soybean and Corn with different vegetation features (e.g., canopy structures)
would increase the spatial heterogeneity in the same study area (even within the same satellite pixel). This
may rise the difficulties of optical modeling for vegetation compared to that in homogenous area (e.g., the
full grasslands for Dahra). Based on the discussions concerning the spatial mismatching and heterogeneity
issue, it is fairly reasonable to understand the satisfying performances in Dahra, EucFACE and Qinghai, but
more uncertainties for SMAPVEX16. Finally, the definition of CWC with regards to what the optical satellite
is sensitive to should be also taken into account. In general, the spectral saturation issue in high-density
vegetation traits would inherently limit the optical signals to completely model the whole canopy (Allen and
Richardson 1968; Tucker 1977; Mutanga et al. 2023) and contribute to some uncertainties against the
ground destructive measurements of the whole vegetation canopy. The denser the vegetation, the stronger the
saturation issue as well as the more uncertainties would be introduced into the CWC assessments. This could
partially explain why we get the better results in grasslands, then followed by croplands, and forests. In other
words, the CWC retrievals based on PROSAIL model could satisfyingly represent the ground CWC of whole
canopy in low and most medium densely vegetated areas (more refer to all the water of the canopy green
components), but for high-dense vegetation traits, only exhibiting the leaf water of canopy considering green
components mainly come from canopy leaf, also with the saturation issues due to the penetration capacity of
optical signals. In this case, the integration with microwave observations should be promising.
# Chapter 8

## Optical Datacube (WP 8-9)

### 8.1 Introduction and purpose

The purpose of this section is to detail the setup of the Optical datacubes processing chains and their configurations. Biophysical products such as LAI, Fapar, Fcover have been downloaded from different sources and aggregated at the same spatial resolution to compare them to the VOD data sets derived from passive microwave measurements and to assess the complementary contribution of their content in the understanding of the VOD variability. The following sections represents the activities developed for the work package 8 and 9 of the project.

### 8.2 Input/Output data

Table 8.1 below summarizes the information about input optical datasets that were downloaded and used for the creation of the Optical Datacubes 1 and 2.

The output data are Optical Datacubes 1 and 2 with following specifications :

* **Datacube #1** : MODIS Climate Monitoring Grid, latitude/longitude projection, Global coverage [-180°W; +180°E] and [-90°S; +90°N], 0.05° spatial resolution (image size 3600H x 7200W)

* **Datacube #2** : EASE2 25km grid, latitude/longitude projection, coverage [-180°W; +180°E] and [-83.5°S; +83.5°N], 0.25° spatial resolution (image size 1388W x 584H)
<table><caption>Table 8.1: Input data for the Optical Datacubes #1 and #2</caption><thead><tr><th></th><th>Satellite</th><th>Product title</th><th>Provided by</th><th>Temporal range</th><th>Temporal resolution</th><th>Grid / projection, spatial resolution and coverage</th></tr></thead><tbody><tr><td>1</td><td>AVHRR</td><td>LAI, FAPAR, FCOVER</td><td>POSTEI, Theia, CNES</td><td>1981-2020</td><td>10 days</td><td>Geographic lat/lon projection, 0.05°, [-180°W; +180°E], [-90°S; +90°N]</td></tr><tr><td>2</td><td>Proba V, Sentinel-3 OLCI</td><td>LAI, FAPAR, FCOVER (v1), NDVI (v2)</td><td>Copernicus Land Service</td><td>2014 - present</td><td>10 days</td><td>Plate carrée, 1/112° (1km), [-180°W; +180°E], [-60°S; +75°N]</td></tr><tr><td>3</td><td>VGT1, VGT2, Proba V</td><td>LAI, FAPAR, FCOVER (v2), NDVI (v3), Top Of Canopy (TOC) Reflectance (v1)</td><td>Copernicus Land Service</td><td>1999 - 2020</td><td>10 days</td><td>Plate carrée, 1/112° (1km), [-180°W; +180°E], [-60°S; +75°N]</td></tr><tr><td>4</td><td>MODIS Terra</td><td>MOD13C1 Vegetation Index (VI) (v6.1)</td><td>United States Geological Survey (USGS)</td><td>2000 - present</td><td>Multi-day</td><td>Geographic lat/lon projection, 0.05°, [-180°W; +180°E], [-90°S; +90°N]</td></tr><tr><td>5</td><td>MODIS Terra and Aqua (combined)</td><td>MCD43C4 BRDF (v6.1)</td><td>United States Geological Survey (USGS)</td><td>2000 - present</td><td>Daily product</td><td>Geographic lat/lon projection, 0.05°, [-180°W; +180°E], [-90°S; +90°N]</td></tr><tr><td>6</td><td>MODIS Terra and Aqua (combined)</td><td>MCD15A2H LAI, FAPAR (v6.1)</td><td>United States Geological Survey (USGS)</td><td>2002 - present</td><td>Multi-Day</td><td>Sinusoidal grid (tiles of 10°x10° at the equator), 500m, [-180°W; +180°E], [-90°S; +90°N]</td></tr><tr><td>7</td><td>MODIS Terra and Aqua (combined)</td><td>MCD18C2 PAR (v6.1)</td><td>United States Geological Survey (USGS)</td><td>2002-present</td><td>Daily product</td><td>Geographic lat/lon projection, 0.05°, [-180°W; +180°E], [-90°S; +90°N]</td></tr></tbody></table>

Based on the Datacubes specifications and the source data projections, the collections that had to undergo processings were as follows:

* projection for the Datacube 1 : Copernicus Land Service vegetation variables and the MCD15A2H vegetation variables

* projection for the Datacube 2 : all of the collections

## 8.3 Datacube #1 reprojection algorithm

Based on the Table 8.1, certain optical datasets were identified for the reprojection in order to match the projection and spatial resolution of the Datacube #1. This first geometrical reprojection was done mainly
with the gdalwarp tool. gdalwarp is an image mosaicing, reprojection and warping utility (Fig. 8.1). This program can reproject data to any supported projection. Multiple resampling methods are offered by the tool, and the resampling kernel size also depends on the chosen method.

<figure><img src="image_60.png" /><figcaption>Figure 8.1: Example of the projection with gdalwarp from a more fine spatial resolution (left) to a more coarse spatial resolution (right)</figcaption></figure>

Amongst multiple resampling methods offered by gdalwarp, some of these methods can induce information loss dû to the kernel size and the resampling method itself (Fig. 8.2). Lost information is the information coming from pixels that do not fall into the resampling kernel.

<figure><img src="image_61.png" /><figcaption>Figure 8.2: gdalwarp 2x2 bilinear kernel with the hatched zones showing the pixels where the information has been lost during the resampling</figcaption></figure>

This problem can be solved either by:

1. changing the resampling method and taking a bigger kernel which would take into account a larger number of pixels at once (cubic resampling with 3x3 kernel, for example)
or

2. adding an averaging convolution step (as a low-pass filter) before the projection.

The second method is visualised on the Fig. 8.3. It shows that adding the convolution step with a 3x3 mean kernel allows to use the bilinear resampling algorithm without losing any pixel that falls in-between the resampling kernels.
<figure><img src="image_62.png" /><figcaption>Figure 8.3: Bilinear gdalwarp convolution with the preliminary mean-convolution step</figcaption></figure>

The example above is a demonstration of a particular example but the explanation of the convolution step can be generalized. When the difference between input and output spatial resolutions is too big, the change of the resampling method with a bigger kernel won't necessarily help, and aliasing artefacts may occur int final result. The robust solution in this case is to add a convolution step before projection. This averaging convolution step has to summarize data of the neighbourhood of each pixel in a mean value, and prevent an aliasing in the final image, as it will play the role of the low-pass filter. The convolution has to bring the input spatial resolution closer to the final one, but not over-filter the data. This anti-aliasing concept is the same that the one used in the processing of 1D signals. On top of the possible information loss, the aliasing can also be the reason of the reshaped or modified forms or curves which can falsify the original data (Fig. 8.4).

<figure><img src="image_63.png" /><figcaption>Figure 8.4: Example of aliased and anti-aliased results on a test image</figcaption></figure>

Considering the resolutions of the input data for the Datacube #1, the reprojection of the data contained two main steps:

1. convolution of the data to the "intermediate" spatial resolution and,

2. the reprojection itself

These steps were the same for every reprojected collection with some slight differences in the configuration of the convolution kernel size which depends on the gap between source and output spatial resolutions.

### 8.3.1 Convolution to the intermediate resolution

The averaging kernel is defined as follows (Eq. 8.1), where dimensions *N* and *M* are defined based on the input/output spatial resolutions:
$$N \left\{ \underbrace{\begin{bmatrix} 1 & \dots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \dots & 1 \end{bmatrix}}_{M} \cdot \frac{1}{NM} \right. \qquad (8.1)$$

After that, the convolution method is a usual kernel convolution algorithm, but with coefficients that are allows to compute the mean value of the neighbourhood (Fig. 8.5):

<figure><img src="image_64.png" /><figcaption>Figure 8.5: 2D image convolution concept</figcaption></figure>

Two kernel sizes were defined for the processing:

1. 5x5 kernel for the Copernicus 1km products

2. 7x7 kernel for the MODIS MCD15A2H 500m tiles

For the majority of the data, main variables are “physical” quantitative data (reflectances, vegetation indices), and the convolution can be applied as it is. However, for some other additional variables, like quality flags, the data is not something physical or continuous; it won't make much sense to compute a mean value or do a projection on such qualitative information. Thus, for this type of data, a special approach was developed in this project and is discussed below.

### 8.3.1.1 Special case : convolution of the non-physical qualitative variables

For the quality/classification flags, the convolution is not meant to compute the mean value of the neighbourhood, but the percentage of a chosen flag in the region. Fig. 8.6 shows an example of this approach. For a given flag description or mask value, the kernel compute the contribution percentage of this flag in the final projected pixel. Then, this percentage can be projected with gdalwarp the same way as other variables. The size of the kernel is the same as the one that was used for the physical variables of the same product collection.
<figure><img src="image_65.png" /><figcaption>(a) Original image with flag value LAND/WATER per pixel</figcaption></figure>

<caption>(b) Convolved image with LAND pixels percentage as the result</caption>

<caption>Figure 8.6: Convolution to the intermediate resolution in a case of a qualitative variable</caption>

## 8.3.2 Reprojection with gdalwarp

As it was already mentioned above, the projection step was done via the gdalwarp utility. This tool is also helpful for merging tiles and image parts together or for reprojection between different spatial references. In the scope of this project the functionality of merging tiles was used for the MODIS MCD15A2H collection, where original products are computed on the sinusoidal grid (Fig. 8.7).

<figure><img src="image_66.png" /><figcaption>Figure 8.7: MODIS sinusoidal tile grid</figcaption></figure>

gdalwarp is a command line tool, and it accepts multiple input options that allow to tune the projection to one's needs. Below is the verbatim of the gdalwarp help function.

```
Usage: gdalwarp [--help-general] [--formats]
       [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid]
       [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
       [-refineGCPS tolerance [minimumGCPS]]
       [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
       [-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...]
       [-wt Byte/Int16] [-srcnodata "value [value...]"] [-dstnodata "value [value...]"]
       [-dstalpha [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
       [-cutline datasource] [-cl layer] [-cwhere expression]
       [-csql statement] [-cblend dist_in pixels] [-crop_to cutline]
       [-of format] [-co "NAME=VALUE"]* [-overwrite]
       [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*
       [-doo NAME=VALUE]*
   srcfile* dstfile
```
Available resampling methods:

near (default), bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, Q1, Q3.

In our case the input parameters was defined as follows:

* `-t_srs` target spatial reference, was extracted from MODIS products which were already in a Datacube #1 spatial reference and defined by

```
GEOGCS["Unknown datum based upon the Clarke 1866 ellipsoid",
DATUM["Not specified (based on Clarke 1866 spheroid)",
SPHEROID["Clarke 1866",
6378206.4,
294.9786982138982,
AUTHORITY["EPSG","7008"]]]],
RIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]]
```

* `-te`, georeferenced extents of output file, was defined by "-180 -90 180 90" as the products are of global coverage

* `-tr`, output spatial resolution, was defined by "0.05 -0.05"

* `-r`, resampling method, was defined by "bilinear"

* `-dstnodata` and `-srcnodata`, no-data masking values respectively in the output and input data, was defined by "-999" for both because the convolution step provides the no-data pixels with -999 value and for the harmonization of the data this value was also kept for the Datacubes as the no-data value

* option `-co` defined with `WRITE BottomUp=NO`, this is a NetCDF driver option which allows to set the y-axis order for export. In order to export images with the "north-up" convention, this option has to be set to NO.

The input data for the `gdalwarp` utility is the convolved data, and the output data is the final projected data that is concatenated afterwards along time dimension into the Datacube #1. Source data that was already in the Datacube #1 projection and on the same grid was directly concatenated into the Datacube #1 format without `gdalwarp` or low-pass filter processing.

## 8.4 Datacube #1 validation

The first Optical Datacube contains the main data of the product, projected onto a grid of 0.05° spatial resolution (image size 3600H x 7200W), as well as some metadata and additional variables which give additional information on the products used for the Datacube aggregation. Global and variables' attributes are validated by manual checking. Variables' values are validated using global map visualisation and comparisons of time series "first Optical Data Cube VS source data" for given coordinates on the globe. In this section, we will describe the validation methodologies for two different scenarios:

1. same spatial resolution between source data and Datacube #1

2. source data is on the different grid and with different spatial resolution than Datacube #1
## 8.4.1 Products with grid shared by source data and Datacube #1

### 8.4.1.1 Methodology

For the case where no reprojection is needed, variables are validated using:

* comparison of global map from first optical datacube and source data, generally for one moment in January and one in August. Absence of artefacts and spatial coherence between source data and first optical data cube are checked.

* comparison of time series data (over a year) from the first optical datacube and from the source data for a specific position. In this case, when the grid is shared by source data and first optical datacube, a pixel in the source data corresponds to the same pixel in the first optical datacube. Therefore, for a given period and position, the variability of variables in this datacube must be the same as the one in the corresponding optical database.

For this second step, time series over three hundred LANDVAL sites were systematically created for all of the collections, in order to assess the Datacube #1 results. The distribution of these sites is showed on Fig. 8.8. These reference sites are well known for their high temporal stability and are used by CNES for the absolute calibration of remote sensing sensors, thus allowing us to have a strong basis of locations where we can test the results of our Datacube processing. For this validation five sites were selected amongst the three hundred in order to represent a diverse range of cases in terms of time series patterns and variable range values allowing a more quick and efficient study :

* Aguascalientes (21.7°N/-102.32°E)

* Aire Adour (43.7°N/0.25°E)

* Arm_Cart_SGP (36.64°N/-97.5°E)

* BAC_LIEU (9.28°N/105.73°E)

* Belmanip_00075 (41.58°N/-77.85°E)

<figure><img src="image_67.png" /><figcaption>Figure 8.8: Map of the 308 LANDVAL ground validation sites (left) and the five selected station for the validation (right)</figcaption></figure>

### 8.4.1.2 Results

* AVHRR LAI product

On Figure 8.9 no artefacts are observed, and there is a good match between the two global maps.

On the Figure 8.10, the variability of the LAI over the year 1983 is exactly the same for the first optical data cube and the source data as expected. Same validation methodology was applied to the FCOVER and FAPAR products of the AVHRR collection. The global maps comparison and the analysis of the time series provided the same conclusion as the one demonstrated for LAI product : good match is observed between maps and
the variability of the data is exactly the same between Datacube #1 and the source data for the mentioned variables of the AVHRR collection. The comparison was also declined for another date in the summer season to observe seasonal changes, and this comparison also demonstrated a good match between the source and the Datacube # 1 data. The projection of the data was validated.

<figure><img src="image_68.png" /><figcaption>Figure 8.9: Global visual comparison between AVHRR LAI original and harmonized data on January 5, 1983</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Figure 8.10: LAI time series (year 1983) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations.</figcaption></figure>

• MOD13C1 NDVI product

For the MOD13C1 NDVI product, good matches were observed on the global maps for the winter/summer seasons, and good match were also observed for the time series of the five selected stations. The projection of the data was validated.

<figure><img src="image_70.png" /><figcaption>Figure 8.11: Global visual comparison between MOD13C1 NDVI original and harmonized data on a February 18, 2000</figcaption></figure>
<figure><img src="image_71.png" /><figcaption>Figure 8.12: Global visual comparison between MOD13C1 NDVI original and harmonized data on a September 5, 2000</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 8.13: MOD13C1 NDVI time series (year 2000) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

• MCD18C2 Photosynthetically Active Radiation (PAR) product

Good matches were observed on the global maps for the winter/summer seasons, and good match were also observed for the time series of the five selected stations. For the illustration purposes, only plots of one of the variables are presented below. The projection of the data was validated.

<figure><img src="image_73.png" /><figcaption>Figure 8.14: Global visual comparison between MCD18C2 GMT 09:00 PAR original and harmonized data on an August 25, 2002</figcaption></figure>
<figure><img src="image_74.png" /><figcaption>Figure 8.15: MCD18C2 PAR time series (year 2000) comparison between first optical data cube (in red) and source data (in blue) for the selected stations</figcaption></figure>

• MCD43C4 Nadir Reflectance product

Good matches were observed on the global maps for the winter/summer seasons, and good match were also observed for the time series of the five selected stations. For the illustration purposes, only plots of one of the variables are presented below. The projection of the data was validated.

<figure><img src="image_75.png" /><figcaption>Figure 8.16: Global visual comparison between MCD43C4 Global Nadir Reflectance (Band 4) original and harmonized data on a February 25, 2000</figcaption></figure>
<figure><img src="image_76.png" /><figcaption>Figure 8.17: Global visual comparison between MCD43C4 Global Nadir Reflectance (Band 4) original and harmonized data on an July 29, 2000</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Figure 8.18: MCD43C4 Global Nadir Reflectance (Band 4) time series (year 2000) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

## 8.4.2 Products with different grids and projections

### 8.4.2.1 Methodology

In this specific case, quantitative variables (such as LAI, FCOVER...) and qualitative variables (such as quality flag) require validation using different methodologies:

* quantitative variables are validated using the same methodology as described in the previous section (global map and time series comparisons)

* qualitative variables are manually validated through pixel-to-pixel comparisons

As described in the section 8.4.1.1, a comparison is performed between the global map generated from the first optical datacube and the source data (when global map is available). Nothing changes for this method.

However, in order to perform a time series analysis and comparison with data with different spatial resolutions, following methods were applied:

1. find the nearest neighbour pixel to the ground validation site in the original source data (via minimal geodesic distance), this pixel is represented by the red square on Fig. 8.19

2. identify the Datacube #1 pixel that includes the source station pixel and retrieve its center; this pixel is indicated by the yellow circle on Fig. 8.19
3. retrieve from the original source data the 5x5 km² neighborhood around the DC1 found pixel. The size of this neighbourhood in pixels will depend on the original data spatial resolution. In the example of Fig. 8.19 the 5x5 neighborhood is 25 Copernicus pixels because the original Copernicus spatial resolution is 1km.

4. calculate the mean value and standard deviation of pixel values, as well as the number of pixels with available data *within the selected area around the DC1 central pixel*

5. compare time series from the DC1 center pixel, mean +/- std distribution and the original source data pixel

<figure><img src="image_78.png" /><figcaption>Figure 8.19: Diagram of time series comparison methodology (example for a Copernicus product with 1x1 km² spatial resolution)</figcaption></figure>

An example of such method visualization (applied to the Copernicus LAI product) is presented on Figure 8.20. On the left plot, the legend is as follows:

**red dashed line** : LAI value from the Datacube #1 pixel center

**green line** : LAI value from the source Copernicus data (nearest neighbour pixel to the station)

**blue line** : the average (mean) value of the 25 source pixels surrounding the DC1 center pixel

**shaded blue area** : mean +/- standard deviation of the 25 source pixels surrounding the DC1 center pixel

To validate the first Copernicus optical data cube we expect to have the red and green line within the blue area.

If it is not the case, the distribution of a variable is examined more closely to determine if it can be explained by the source data distribution in the 5x5 km² neighbourhood (for example, at the decimal year 1999.7 on Fig. 8.20). This case is illustrated on the right plot of the Figure 8.20. The blue curve represents LAI pixel values distribution from the 25 1km pixels that fall into the neighbourhood around the DC1 center pixel (green zone on Fig. 8.19). The central grey dashed line corresponds to the distribution mean and the grey shaded area represents the mean +/- std of the LAI distribution. The red point corresponds to the source data pixel containing the LANDVAL site (red square on Fig. 8.19). This pique explains the less strong match between the Datacube #1 pixel values and the station 1km pixel values around the dates of decimal year 1999.7 (left plot on Fig. 8.20).
<figure><img src="image_79.png" /><figcaption>Figure 8.20: Example of time series analysis for the Copernicus 1km LAI dataset</figcaption></figure>

In the following subsection, we present detailed comparisons between the first optical datacube and source data for each of the collections that have undergone the reprojection. Firstly, we compared global maps, followed by a comparison of time series. The results showed a strong agreement between the first optical data cube and source data.

### 8.4.2.2 Results

* LAI Copernicus Global Land Service product

<figure><img src="image_80.png" /><figcaption>Figure 8.21: Global visual comparison between Copernicus LAI original and harmonized data on January 10, 1999</figcaption></figure>
<figure><img src="image_81.png" /><figcaption>Figure 8.22: Copernicus LAI time series (year 1999) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

We can observe that on all five plots of Figure 8.22 the red curve is well within the shaded blue area. Upon examining the values of the 1 km resolution pixel, it is evident that for the "ARM CART SGP" and "BELMANIP_00075", there are data points that fall outside the blue zone. As discussed in section 8.4.2.1, this can be verified and analysed by plotting the distribution of the 25 Copernicus pixels included within the first optical data cube's pixel. As it was seen on Figure 8.20, the value of the 1 km Copernicus pixel (red point) can fall outside the mean of the distribution, causing the occasional mismatch between the 1km pixel value and the Datacube #1 pixel value for a given station.

Same validation was performed on the FCOVER and FAPAR products of the Copernicus collection. Similar results were observed whether studying the global maps distribution or analyzing time series plots. The projection was validated.

* Top-of-Canopy Reflectance (TOCR) Copernicus Global Land Service product

Good match and no artefacts were observed on the global maps. During the time series analysis, the Datacube #1 values are well within the mean ± std value range. The Copernicus 1 km pixel value sometimes falls outside of that range due to the data reprojection and the reasons explained in the section 8.4.2.1. The projection was validated.

<figure><img src="image_82.png" /><figcaption>Figure 8.23: Global visual comparison between Copernicus Global TOCR original and harmonized data on a January 3, 1999</figcaption></figure>
<figure><img src="image_83.png" /><figcaption>Figure 8.24: Global visual comparison between Copernicus Global TOCR original and harmonized data on a May 3, 1999</figcaption></figure>

<figure><img src="image_84.png" /><figcaption>Figure 8.25: Copernicus TOCR time series (year 1999) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

• NDVI Copernicus Global Land Service product

Good match and no artefacts were observed on the global maps. During the time series analysis, the Datacube #1 values are well within the mean ± std value range. The Copernicus 1 km pixel value falls sometimes outside of that range due to the data reprojection and the reasons explained in the section 8.4.2.1. The projection was validated.

<figure><img src="image_85.png" /><figcaption>Figure 8.26: Global visual comparison between Copernicus Global NDVI original and harmonized data on a January 1, 1999</figcaption></figure>
<figure><img src="image_86.png" /><figcaption>Figure 8.27: Global visual comparison between Copernicus Global NDVI original and harmonized data on a May 1, 1999</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>Figure 8.28: Copernicus NDVI time series (year 1999) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

• MCD15A2H LAI product

The MCD15A2H FPAR (Combined Fraction of Photosynthetically Active Radiation) and LAI (Leaf Area Index) product is an 8-day composite dataset with 500 meter pixel size. Due to the fine spatial resolution, the MCD15A2H products use the Sinusoidal grid tilling system (Fig. 8.7). In consequence, after merging the tiles together into a global map, we used the methodology presented in section 8.4.2 for this product. We will present global map and time series comparison. For each time series comparison the legend is the same as explained in the section 8.4.2.1.

The original MCD15A2H data is provided on the very fine Sinusoidal grid tilling system, hence only the analysis of the global map from the Datacube #1 was possible (Fig. 8.29). The use of the others collections' LAI products can help to assess the coherence of the observed tendencies. Based on the comparison on Fig. 8.29, we can see no artefacts and incoherent data for the MCD15A2H reprojected product. Zones with higher values are less homogeneous on the MODIS map compared to the same zones on the AVHRR map but it is due to the more fine resolution of the source MODIS data; the more fine resolution brings us more details while the AVHRR data gives a more averaging result with its more coarse resolution. Nevertheless, the tendencies between the two maps are coherent and follow the same pattern.

We can observe on Figure 8.30 that in all five figures, except for Aguascalientes station, the red curve is well within the shaded blue area. Upon examining the values of the 500 m resolution pixel, it is evident that there are data points that fall outside the blue zone. It can be explained by the 500m pixel distribution.

Based on the global map observations and time series analysis, no artefacts and incoherent data were detected for this product. The projection was validated.
<figure><img src="image_88.png" /><figcaption>Figure 8.29: Comparison of the global LAI maps from the MCD15A2H and AVHRR Datacubes #1</figcaption></figure>

<figure><img src="image_89.png" /><figcaption>Figure 8.30: MCD15A2H LAI time series (year 1999) comparison between first optical data cube (in red) and source data (in blue) for the five selected stations</figcaption></figure>

### 8.4.2.3 Conclusion

In this section, we have presented detailed comparisons between the first optical data cube and source data for each data cube. Firstly, we compared global maps, followed by a comparison of time series. As expected to validate the first optical data cube, the results showed a strong agreement between the first optical data cube and source data.

## 8.5 Datacube #2 reprojection algorithm

In order to use the project optical data presented above to improve our understanding of the vegetation opacity data sets retrieved from SMOS measurements and other low frequency passive microwave measurements, it is necessary to produce this optical data at 0.25° spatial resolution (SMOS spatial resolution). To do this, a Second Optical Datacube at 0.25° spatial resolution were created from the Datacube #1 presented above.

For this second datacube, the input is the first datacube itself, where all of the data has already been harmonized and projected onto the same intermediate grid. In order to obtain the Datacube #2, the first datacube data has to undergone a new re-projection where the objective is to smooth the data to the 0.25° spatial resolution via a kernel of a weighted average. This kernel is meant to represent the behaviour of the SMOS antenna. The weights of the kernel are calculated following a predefined formula elaborated by the CESBIO colleagues of this project. The section aims to detail different steps of this reprojection algorithm.
### 8.5.1 Kernel computation

The mathematical approximation of the SMOS antenna coefficients that multiply received signal follows the formula 8.2:

$$w = \frac{(\mathrm{sinc}(\frac{0.0733 \cdot d}{\pi}))^{1.4936}}{1 + d^{2.1030} \cdot 2.5748118104204e^{-04}} \quad (8.2)$$

Where d is the distance between the central point right underneath the antenna and any other point in the neighbourhood.

For a point with coordinates $(x_0, y_0)$ in 2D dimension, we define the distance to another point $(x, y)$ by the formula 8.3:

$$d = \sqrt{(x_0 - x)^2 + (y_0 - y)^2} \quad (8.3)$$

For a latitude/longitude point, we can define the distance respecting the curvature of the globe by the formula

$$d = R \cdot 2 \cdot \arctan \left( \frac{\sqrt{\sin^2(d_{\text{lat}}/2) + \cos(lat_0) \cdot \cos(lat_n) \cdot \sin^2(d_{\text{lon}}/2)}}{\sqrt{1 - (\sin^2(d_{\text{lat}}/2) + \cos(lat_0) \cdot \cos(lat_n) \cdot \sin^2(d_{\text{lon}}/2))}} \right) \quad (8.4)$$

where:

* $d_{\text{lat}}$ and $d_{\text{lon}}$ are respectively the latitude and longitude distance in radians between two points

* $lat_0$ and $lat_n$ are respectively latitude values of the reference point and the second point

The kernel is applied to the Datacube #1 data on the Datacube #1 grid. For any point of the grid we can compute a matrix with distances to this point's neighbours (Fig. 8.31).

<figure><img src="image_90.png" /><figcaption>Figure 8.31: Notion of distance matrix in the kernel computation methodology</figcaption></figure>

Then, we can apply the formula 8.2 from above to compute kernel weights for each of the neighbour points based on computed distances. An example of such distance matrix and associated kernel is presented on Figure 8.32. In this example the neighbourhood is 140 km² ($x$ and $y$ varying from -70 to 70km from the central point).
<figure><img src="image_91.png" /><figcaption>Figure 8.32: Example of the distance matrix and the associated kernel matrix with normalized weights</figcaption></figure>

Giving the formula 8.3, we can already notice that the kernel weights depends only on the latitude/longitude distance between points and their latitudes, and does not depend on points' longitude values. Hence, for two points of the same latitude on a regular grid, the distances between these points and their respective identical neighbourhood will be the same. Thus, the distance matrix and the kernel could only be computed once for each of the different latitude values of the grid in order to speed up the data processing.

Regarding the kernel coefficients formula 8.2, there are three important additional steps.

First of all, as it can be seen on the right plot of Figure 8.32, there are blank pixels where the weights are masked, even if the distance values are present. This is due to the cut-off distance that is defined by the user, and beyond which the weights are considered null. In our case this distance is 61.5 km (thus maximum allowed neighbourhood of 123x123 km²).

Second, the mathematical aspect of the formula allows for the imaginary numbers. It was decided to not take into account the eventual imaginary part of the coefficients, and resolve this problem by replacing all of the imaginary numbers by the minimum real-only number present in the kernel.

Lastly, the weights can be normalized to make the convolution easier and to adapt it for a real physical variable use-case. To do so, we divide each coefficient by the sum of all of the coefficients, so the sum of the final weights makes 1 in total.

In summary, the kernel computation should go through following steps in this order:

1. compute the distance matrix

2. compute the weights following the formula

3. replace imaginary numbers by the minimum of the real-only coefficients

4. replace coefficients beyond the cut-off distance by 0

5. normalize weights

## 8.5.2 Convolution algorithm

The convolution algorithm contains the kernel calculation and the convolution itself which follows the classical convolution method (Figure 8.33), also used in the Datacube #1 production chain. In our case, the input is the Datacube #1 data, the filter is the SMOSE kernel with normalized weights and the resulting values will constitute the Datacube #2.
<figure><img src="image_92.png" /><figcaption>Figure 8.33: 2d image convolution concept</figcaption></figure>

For each of the SMOS EASE grid points, the algorithm is as follows:

1. for a point of the SMOSE EASE grid (which is the grid of the Datacube #2) we find the nearest neighbour in the Datacube #1 grid (which is constant for every Datacube 1 product). This found nearest neighbour will then be the central starting point on the Datacube #1 grid for further computations (Fig. 8.34).

<figure><img src="image_93.png" /><figcaption>Figure 8.34: Datacube 1 (blue crosses) vs SMOSE (red points) grid, nearest neighbours are circled in green</figcaption></figure>

2. from the found central point, we compute the distance matrix via the formula 8.3 and cut-off any pixels with distance > cut-off distance as they won't be useful in the processing

3. based on the distance matrix, we compute the weights as described in the previous section (coefficient 0 if distance is still > cut-off distance, no imaginary numbers, normalized coefficients)

4. we perform the classical convolution between the computed kernel and the data sub-extract of the Datacube #1 which is overlayed by the kernel and centered on the found nearest neighbour

5. the resulting value corresponds to the product re-projected data value on the SMOSE point coordinates from which we started the step n°1

### 8.5.3 Invalid/no data pixels contribution

In case when the data overlaid by the kernel contains some no-data pixels or fill values, these pixels are not taken into account during the convolution, the resulting value contains only valid values multiplied by the corresponding coefficients. On Figure 8.35 no-data pixels from the left are propagated to the plot on the right as they must not be taken into account. Nevertheless, the kernel is always computed without taking into account valid/invalid pixels as its purpose is to simulate the instrument behaviour that does not adapt its functioning depending on the seen/received information.
<figure><img src="image_94.png" /><figcaption>Figure 8.35: Example of the data convolution with no data pixels</figcaption></figure>

### 8.5.4 Special case n°1 : uncertainties/error variables

In the case of the variables containing uncertainties or error estimation values, the convolution concept remains the same except that the kernel coefficients have been changed so as to calculate the classic average of the kernel covered data. The invalid data is treated in the same manner as explained in the section 8.5.3. The valid data that remains after applying kernel and filtering the invalid data (Fig. 8.36) is averaged to the final value. Kernel's averaging coefficients are calculated on the basis of the valid data only.

<figure><img src="image_95.png" /><figcaption>Figure 8.36: Data masking for the average computation of the uncertainties/error estimation variables</figcaption></figure>

### 8.5.5 Special case n°2 : quality flags

The case of quality flags and other qualitative variables were treated separately as it required a different approach from the physical/quantitative variables presented above. Two scenarios were studied for this part :

1. no projection were performed between source data and Datacube #1. In this case, qualitative variables have not yet undergone any processing.

2. projection between source data and Datacube 1 were performed

#### 8.5.5.1 No projection between source data and Datacube 1

In the case of source data being already in the Datacube #1 projection, no data processing were performed. This data did not undergo any changes. Thus the question of how to project the qualitative data onto the SMOSE grid was addressed.

The retained solution was the same as for the case where we had to reproject the source data into the Datacube 1 grid (section 8.3.1.1): for each flag or qualitative data that we are interested in, we want to compute the percentage of contribution of such pixels to the final SMOS pixel, using the same neighbourhood that in the case of a physical variable (as explained in previous paragraphs). An example is presented on Figure 8.37.

Here we have a flag that represents whether a pixel in the Datacube #1 AVHRR product is recognized as Evergreen Broadleaf Forest or not. When we couple the data extract with a kernel, we end up with the same
neighbourhood as in the physical variable reprojection case. Then, in the remaining data we calculate the contribution as :

$$c = \frac{N_{flag pixels}}{N_{valid pixels}} \qquad (8.5)$$

The resulting value is then affected to the SMOS current pixel, in the same manner as for the physical variable convolution.

<figure><img src="image_96.png" /><figcaption>Figure 8.37: Example of a quality flags vs kernel overlay</figcaption></figure>

### 8.5.5.2 Projection between source data and Datacube 1

In the case where the data needed reprojection between the source grid and the Datacube 1 grid, the quality flags processing has already been assessed. The variables in the Datacube 1 product already represent the contributions of the flags to the Datacube 1 pixels.

Here, we applied the same technique as for the uncertainties and errors data: we compute the classic average of the contributions in the same kernel neighbourhood as for the physical variables.

## 8.6 Datacube #2 validation

The second Optical Data Cube contains the main data of the product, projected onto a grid of 0.25° spatial resolution (image size 1388W x 584H), as well as some metadata and additional variables which give information on the products used for the aggregation of the data cube.

Global and variable attributes are validated by manual checking. Other variables are validated using global map visualisation and with comparisons of time series from the second Optical Data Cube and the first optical data cube for a given position.

### 8.6.1 Methodologies

#### 8.6.1.1 Global map comparison

A comparison is performed between the global map generated from the second optical data cube at 0.25° resolution and the first optical data cube at 0.05° resolution. The aim is to check the absence of artefacts.

#### 8.6.1.2 Time series comparison

For a more detailed and accurate comparison between the datacubes, following processing were performed :
1. find the DC2 nearest neighbour pixel to the ground validation site (via the minimal geodesic distance)

2. find the DC1 nearest neighbour pixel to the DC2 pixel center and retrieve its value with the 25 adjacent DC1 pixels

3. calculate the mean and standard deviation of the 25 DC1 pixels, as well as the number of pixels with available data within the selected area

4. find the DC1 nearest neighbour pixel to the ground station

5. compare the time series from DC2, DC1 and from the DC1 25 pixel neighbourhood

It's important to note that in this method we average the 25 nearest pixels and don't use the convolution kernel, which can lead to significant differences in heterogeneous areas. This point will be taken into account in the analysis of the results.

For each time series comparison, the legend is as follows:

**red dashed line**: variable value from the DC2 center pixel (validation site nearest neighbour)

**green dashed line**: variable value from the DC1 center pixel (validation site nearest neighbour)

**blue line**: average value of the 25 DC1 pixels surrounding the center of the DC2 pixel

**shaded blue area**: mean +/- std value of the 25 DC1 pixels surrounding the center of the DC2 pixel

For this validation activity the time series for the three hundred LANDVAL sites were also systematically extracted from the Datacube # 2 (see section 8.4.1.1). The analysis were performed especially on the same five selected ground sites as for the Datacube # 1 validation, but can be easily derived to other LANDVAL stations for which the time series were already generated.

## 8.6.2 Results

* LAI AVHRR products from POSTEL, Theia (CNES)

No artefacts or incoherence were observed between the Datacube # 1 and 2 plotted data.

<figure><img src="image_97.png" /><figcaption>Figure 8.38: Global visual comparison between AVHRR LAI Datacube #1 and 2 data on a January 5, 2000</figcaption></figure>
<figure><img src="image_98.png" /><figcaption>Figure 8.39: Global visual comparison between AVHRR LAI Datacube #1 and 2 data on a June 30, 2000</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>Figure 8.40: AVHRR LAI time series (year 2000) comparison for the five selected stations</figcaption></figure>

We can see that for each station the patterns are the same for each curve. The green curve, representing the DC1 pixel containing the validation site, is more variable and sometimes falls outside the blue shaded zone. This is explained by the fact that it represents a pixel and may therefore represent an extremum of the value distribution.

In the case of the "Bac Lieu" station, the blue curve is well above the red curve in the middle of the year. This is due to the use of the convolution kernel for the red curve and the use of a simple average over 25 pixels for the blue curve. Convolution is used on many more than 25 pixels, and if the area is very heterogeneous, this can imply significant differences between the simple average and the use of the convolution kernel.

To verify this hypothesis, let's take the following example:
<figure><img src="image_100.png" /><figcaption>Figure 8.41: LAI time series for the "Bac Lieu" validation site and the corresponding pixels values</figcaption></figure>

Figure 8.41 presents the LAI time series for the "Bac Lieu" validation site and the corresponding pixels values. The upper plots represent the time series and the red vertical bar represents the time for which we will plot the geographical area with LAI values. For the bottom plots the legend is:

**white cross** : DC1 grid pixel center

**red diamond** : DC2 grid pixel center

**red dot** : LANDVAL station

**red rectangle** : 25 pixel neighbourhood around the DC2 pixel containing the ground station

We can see that at t = 2000.3 the LAI values are homogeneous in the area of interest, so there will be little difference between using the simple average and using the convolution kernel. This is not the case for t = 2000.5. In this case, the simple average zone includes higher values, but low values are nearby. The latter are taken into account when using the convolution kernel and will therefore lower the average. This is why, in this case, the red curve is below the blue area. This type of analysis was used to study similar cases for other variables. However, for the sake of clarity and concision, this example will be the only one detailed in this document.

For the AVHRR collection, the similar analysis was also performed on the FAPAR and FCOVER products, resulting in a good match between the two Datacubes data. The projection was validated.

* MOD13C1 NDVI products from USGS

The MOD13C1 product introduces two types of variables to validate: physical variable (NDVI) and quality variables (status or quality flags). First of all, a good match was observed between the Datacubes plots for the physical NDVI variable (Figures 8.42 and 8.43 below).
<figure><img src="image_101.png" /><figcaption>Figure 8.42: Global visual comparison between MOD13C1 NDVI Datacube #1 and 2 data on a January 1, 2005</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>Figure 8.43: Global visual comparison between MOD13C1 NDVI Datacube #1 and 2 data on a June 26, 2005</figcaption></figure>

The time series for the NDVI were also analyzed in the scope of this product (Fig. 8.44). The curves are following similar tendencies, except some cases where the disparities can be explained as before by the different averaging kernels and the heterogeneity of the station surroundings.

<figure><img src="image_103.png" /><figcaption>Figure 8.44: MOD13C1 NDVI time series (year 2005) comparison for the five selected stations</figcaption></figure>

Second, the assessment of the quality flag reprojection was performed. In the Datacube #1 this MODIS data had not yet undergone any projection, so the quality flags in this Datacube are the native ones. Pixel values match certain flag or status in the original data, which is :

* 0 = good data
* 1 = marginal data

* 2 = snow/ice

* 3 = cloudy

* 4 = estimated data

In the Datacube # 2 on the other hand, these flags were projected into their respective contribution percentages in the final pixels. Thus, the Datacube # 2 has the "16 days pixel reliability" variable with 5 layers, corresponding to the 5 possible states. Each layer gives the percentage of the number of pixels in the corresponding state that contributed to the final Datacube # 2 pixel :

* 0 = percentage of good pixels

* 1 = percentage of marginal pixels

* 2 = percentage of snow/ice pixels

* 3 = percentage of cloudy pixels

* 4 = percentage of estimated pixels

On Figures 8.45 and 8.46 we can observe a good coherence between the original flags and the projected data. Indeed, for the January 1st, 2005 we can see that all of the yellow/green pixels are tagged as snow/ice pixels with a value of 2, and when projected (middle plot on the Fig. 8.45) this area has the biggest percentage of the snow/ice pixels. Same works for the pixels tagged as "good" pixels (blue pixels on the left plot); on the projected data (right plot) we see that the majority of this surface have significant values of the percentage of "good" pixels. Same observations and analysis have been applied for the data on the June 26, 2005 presented on Figure 8.46.

<figure><img src="image_104.png" /><figcaption>Figure 8.45: Comparison of the MOD13C1 native and projected quality flags in Datacubes # 1 and 2 on January 1, 2005</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>Figure 8.46: Comparison of the MOD13C1 native and projected quality flags in Datacubes # 1 and 2 on June 26, 2005</figcaption></figure>
Based on the observations and analysis from above, the projection of this collection was validated.

* MCD18C2 PAR products from USGS

The comparison of the global maps of the PAR variable were performed for every variable subclass (GMT0900, 1200 etc) but only one of the example is presented here for the sake of the concision. Here, a good match is observed between the global maps of the Datacube # 1 and # 2.

<figure><img src="image_106.png" /><figcaption>Figure 8.47: Global visual comparison between MCD18C2 Datacube #1 and 2 data on a January 1, 2005</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>Figure 8.48: Global visual comparison between MCD18C2 Datacube #1 and 2 data on a July 1, 2005</figcaption></figure>

For the time series analysis we can also see that for each station the patterns are the same for each curve. The green curve, representing the DC1 pixel containing the validation site, is more variable and sometimes falls outside the blue shaded zone. This is explained by the fact that it represents a pixel and may therefore represent an extremum of the value distribution.
<figure><img src="image_108.png" /><figcaption>(a) AIRE ADOUR site (43.7°, 0.25°)</figcaption></figure>

<figure><img src="image_109.png" /><figcaption>(b) BAC LIEU site (9.28°, 105.73°)</figcaption></figure>

Figure 8.49: MCD18C2 GMT0900 PAR time series (year 2005) comparison for two of five stations

• MCD43C4 Nadir Reflectance products from USGS

For the MCD43C4 collection no artefacts were observed while comparing global maps of the data. The very high or low latitude zones (Greenland and South pole) were not projected onto the Datacube #2 grid because the Datacube #2 grid does not have points in these areas.

<figure><img src="image_110.png" /><figcaption>(a) MCD43C4 Nadir Reflectance Datacube #1 data</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>(b) MCD43C4 Nadir Reflectance Datacube #2 data</figcaption></figure>

Figure 8.50: Global visual comparison between MCD43C4 Datacube #1 and 2 data on a January 1, 2005

<figure><img src="image_112.png" /><figcaption>(a) MCD43C4 Nadir Reflectance Datacube #1 data</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>(b) MCD43C4 Nadir Reflectance Datacube #2 data</figcaption></figure>

Figure 8.51: Global visual comparison between MCD43C4 Datacube #1 and 2 data on a July 1, 2005
Similar tendencies between the data were also observed on the time series plots (figures below). The projection of this collection was validated.

<figure><img src="image_114.png" /><figcaption>Figure 8.52: MCD43C4 Nadir Reflectance time series (year 2005) comparison for the five selected stations</figcaption></figure>

• MCD15A2H LAI products from USGS

No artefacts were observed on the global maps for the MCD15A2H collection; good matches between the time series of the data for the five selected stations. The projection was validated.

<figure><img src="image_115.png" /><figcaption>Figure 8.53: Global visual comparison between MCD15A2H LAI Datacube #1 and 2 data on a January 1, 2005</figcaption></figure>

<figure><img src="image_116.png" /><figcaption>Figure 8.54: Global visual comparison MCD15A2H LAI Datacube #1 and 2 data on a June 26, 2005</figcaption></figure>
<figure><img src="image_117.png" /><figcaption>Figure 8.55: MCD15A2H LAI time series (year 2005) comparison for the five selected stations</figcaption></figure>

* Copernicus Global Land Service LAI products

Global map and time series analysis were performed for all of the Copernicus products (LAI, FAPAR, FCOVER, NDVI, TOCR). For the sake of clarity and concision, only the LAI example is presented in the document. Good match and no artefacts were observed during the global map comparison; correct tendencies were observed on the data time series. Same results were demonstrated during other products validation. The projection of this collection was validated.

<figure><img src="image_118.png" /><figcaption>Figure 8.56: Global visual comparison between Copernicus LAI Datacube #1 and 2 data on a January 10, 2000</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>Figure 8.57: Global visual comparison between Copernicus LAI Datacube #1 and 2 data on a June 30, 2000</figcaption></figure>
<figure><img src="image_120.png" /><figcaption>Figure 8.58: Copernicus LAI time series (year 2000) comparison for the five selected stations</figcaption></figure>

### 8.6.3 Conclusion

In this section, we have presented detailed methodologies of production and validation of Optical Datacubes #1 and #2 for each of the studied optical collections. In the validation, we first compared global maps, followed by a comparison of time series. The results showed a strong agreement between the two datacubes, thus validating the projection of the data onto the two grids (intermediate grid and EASE2 grid).
# Chapter 9

# Comparative Data Analysis (WP 10)

## 9.1 Introduction

The WP10 consists in comparing several datasets used or produced in the project. Ideally, CWC should have been compared to VOD and optical indices but due to delays in the design and production this was not possible. Therefore, the comparison was done with different frequencies VOD, optical indices and soil moisture. The goal was to get more insight into how the remote sensing observations can provide information on the hydraulic behaviour of vegetation.

Vegetation Optical Depth (VOD), derived from microwave observations, has gained significant attention in remote sensing studies due to its sensitivity to water content, in contrast to optical indices. Previous crop-related studies using VOD have predominantly concentrated on phenological analysis at specific microwave frequencies and regional scales. This study examines the seasonal dynamics of multi-frequency VODs and NDVI in comparison to that of soil moisture (SM) across global croplands. The time lags between the peaking periods of different VODs, the NDVI, and the SM signal were used to describe the vegetation's response to water availability. An unsupervised machine learning algorithm (K-means clustering) was applied to classify these temporal relationships across global croplands. The analysis identified eleven distinct clusters, six of which could be related to specific regions and crop types, enabling the analysis of their response to microwave radiation as a function of their potential hydric behaviors. Regions such as North America's corn and wheat belts are distinctly identified. A systematic lag between higher frequency VODs and NDVI was observed across all clusters, ranging from 15 to 49 days. On the contrary the relationship between L-Band VOD (LVOD) and NDVI exhibited significant regional variations: In northern North America, LVOD peaked before NDVI, while in areas characterized by dry/warm weather, they were almost synchronous. Time series similarities and relative time shifts were also computed using the dynamic time warping method for a limited number of points. The comparison with the previous time lag analysis shows agreement with respect to the average time lags for the corresponding clusters but not for all of them. This WP provides new insights into cropland water-use strategies linking isohydric (constant water levels) and anisohydric (varying water levels) behaviors to specific crop types and local environmental conditions. This enhances our understanding of passive microwave signals over vegetated areas in agricultural systems worldwide.

## 9.2 Data

### 9.2.1 SMOS LVOD and SM

The SMOS satellite mission (Kerr et al., 2010) is the second Earth Explorer mission, and was launched in 2009. Its main objective is to measure soil moisture (SM) and sea surface salinity (SSS) worldwide. The satellite carries a 2D interferometric radiometer and measures Earth's microwave emission at 1.4 GHz (L-Band) at
multiple angles of incidence, ranging from 0 to around 60 degrees. The ascending (descending) orbit crosses the equator at 6:00 (18:00).

The SMOS LVOD and SM are derived from multiangular brightness temperature (TB) observations via the implementation of the L-MEB (L-Band Microwave Emission of the Biosphere) model (Wigneron et al., 2007b; Kerr et al., 2012b). The retrieved values for both variables are derived simultaneously through an iterative process aimed at minimizing the difference between the observed TB and modeled TB Kerr et al. (2012b). The level 3 products Al Bitar et al. (2017b) (version 330) for LVOD and SM are utilized in this study.

LVOD is a unitless measure that parameterizes the extinction effects (absorption and scattering) of the vegetation layer on surface microwave emissions; lower values represent sparser/lower vegetation and higher values the opposite. SM represents the water content of the superficial layer of soils (around 5cm depth) and is expressed as water volume per unit soil volume (m³/m³). Both products are provided by the CATDS (Centre Aval de Traitement de Données SMOS); the dataset covers daily observations since 2010 up to the present. They are projected onto the EASE (Equal-Area Scalable Earth) grid version 2 (Brodzik et al., 2012b), with pixels of 25x25 km at approximately 30 degrees latitude. In this study, datasets covering the period between 2011 and 2020 are utilized. Even though MW at the L-Band are protected for scientific use, they are still sensitive to Radio Frequency Interference (RFI). In that sense, both the LVOD and SM were filtered. Any pixel with a RFI probability appearance larger than 0.2 was excluded. Additionally, pixels where the χ² (a variable that measures the goodness of fit in modeled TB with respect to observed TB) surpasses 2 K are not taken into account.

## 9.2.2 VODCA VOD

The global long-term Microwave Vegetation Optical Depth Climate Archive (VODCA) (Moesinger et al., 2020b) is a harmonized multisensor time series of VOD products. The dataset comprises three time series for the C- (6.9 GHz), X- (10.7 GHz), and Ku-band (19 GHz), derived from level 2 multisensor products (SSM/I, TMI, AMSR-E, AMSR-2, and WindSat). All VOD values are retrieved using the Land Parameter Retrieval Model (LPRM) (Owe et al., 2008). The C-band VOD (CVOD) time series covers the period 2002-2018, the X-band (XVOD) 1997-2018, and the Ku-band (KuVOD) 1987-2017. The merging strategy is based on Cumulative Density Function (CDF) matching, with the AMSR-E product serving as the scaling reference. Temporal aggregation is achieved by overlapping the available observations. The resulting products have a spatial sampling of 0.25°. As part of the harmonization procedure, inaccurate data in the original time series were removed, considering: RFI (using the level 2 RFI flags from the original datasets), land surface temperature (all observations below 0°C were removed), and unrealistic VOD values (negative VOD values were removed). All VODCA time series were re-gridded to the reference SMOS L3 grid (EASE grid v2) by averaging the original pixels into 25x25 km EASE grid cells. In this study, observations from 2011 on-wards are used.

## 9.2.3 NDVI MODIS

The NDVI from MODIS (Moderate Resolution Imaging Spectro-radiomete, product MYD13C1) (Tucker, 1979; Huete et al., 2002) is used as proxy of vegetation health and photosynthetic activity. The original dataset is a 16 day composite with a spatial resolution of 1 km. Similar to the VODCA time series, the dataset is projected onto the SMOS L3 EASE grid v2. Data spanning from 2011 to 2020 are used in this study.

## 9.2.4 Crop mask

The global crop extent product (GFSAD1KCD) from the Global Food Security-support Analysis Data (GFSAD) (version v001) (Teluguntla et al., 2016) is used in this study. The dataset is produced by the USA's National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS) at a spatial resolution of 1 km. It provides cropland extent, water management practices (irrigated or rainfed), and crop dominance (circa 2010). Crop dominance is determined based on datasets from three separate studies on global crop type distribution. Water management practices are derived from the map of irrigated/rainfed
cropland area produced by the International Water Management Institute. The dataset is taken into the EASE grid v2; the dominant value from the finer resolution map (per pixel) is used as the representative value within the coarser EASE grid.

### 9.2.5 Climate classification mask

The Köppen-Geiger climate classification map (Version 2), provided by GloH2O and developed by (Beck et al., 2023) at 1 km spatial resolution, is used. The map is representative of the period between 1991 to 2020. It is derived primarily from three distinct temperature datasets and four precipitation datasets. The classification comprises five primary climatic groups: A) tropical, B) arid, C) temperate, D) cold (continental), and E) polar climates, each defined by specific temperature and precipitation thresholds.

The original 1 km map was reprojected onto the EASE Grid v2. This transformation involved selecting the predominant value from the higher resolution map to represent the area within each pixel of the EASE grid.

### 9.2.6 IGBP land surface classification

The International Geosphere-Biosphere Program (IGBP) land cover classification is used (Loveland et al., 2000). Specifically, the dataset produced by Broxton et al. (2014b) at 0.5km spatial resolution, which uses the MODIS MCD12Q1 products to create a single map representative of the period 2001-2010. This map is projected into the EASE grid v2 by assigning the predominant class of the original classification into each resulting SMOS L3 pixel. Here, we focus on IGBP cover class 12, which represents croplands. According to Loveland and Belward (1997) and Sulla-Menashe and Friedl, class 12 corresponds to regions 'with temporary crops followed by harvest and a bare soil period' where the cultivated area corresponds to ≥60 % of the pixels surface. Perennial woody crops are not included in this classification. To exclusively focus on regions dominated by croplands, only pixels classified as IGBP land cover class 12 (across all previously mentioned datasets) were used. This masking, in conjunction with the removal of spurious values in the SMOS dataset (as detailed in Subsection 9.2.1), results in approximately 8000 pixels containing available information classified under IGBP class 12.

## 9.3 Methods

### 9.3.1 Scaling and noise removal from time-series

The seasonal behavior of six variables across croplands was studied. The first four served as indicators of water content across vegetation's structural components: LVOD, CVOD, XVOD, and KuVOD. Given the greater sensitivity of lower-frequency microwaves to larger vegetation components, LVOD was considered as a proxy for water content in primary plant structures such as stems, while higher frequencies served as descriptors of smaller components such as small branches and leaves.

The fifth variable, the SM, obtained from the SMOS mission, provided information on water availability for plant uptake. The sixth variable, the NDVI, served as an indicator of croplands' phenology through vegetation greenness.

Before analysis, time-series were pre-treated. To minimize the influence of outliers in the seasonal signal, all values beyond the 0.05 and 0.95 quantiles were removed before smoothing. Given that the original time-series are affected by noise, a 30-day moving average per year was then computed to increase the signal-to-noise ratio.

Once the time-series of each year were smoothed individually, an interannual daily mean was computed. This involved aggregating the daily values across all available years to obtain a single time-series for each pixel; days with less than 3 observations in the period 2011-2020 were removed, just as yearly time-series with less than 180 days.

Finally, and given that the selected variables changed in scale and magnitude, a normalized version of each
variable was used to facilitate the comparison of their seasonal behaviors. The normalization was based on the maximum and minimum values of the interannual mean (per pixel), as follows:

$$ \text{Normalized } X_i = \frac{X_i - \max(X)}{\max(X) - \min(X)} $$

where X is the mean time-series and X_i the value for a given day. Even after normalization, time series can still exhibit noise. To preserve the primary features of the seasonal behavior for each variable and eliminate the influence of local minima or maxima, a smoothing spline is employed to represent the seasonality per pixel.

### 9.3.2 Descriptors of crops hydraulics

After normalizing the time series, the subsequent step involves identifying periods that reflect plant responsiveness to water availability. Periods characterized by optimal conditions for plant growth are determined using SM, a direct measure of water availability for plant uptake. Simultaneously, periods when vegetation is thriving and nearing peak development are assessed using the VODs and the NDVI. The simultaneous analysis of these periods aims to explore the potential synergies between water availability and plant development.

Considering that these critical moments of optimal conditions are dependent on local factors, such as soil type and meteorological conditions, a local threshold was established for all normalized time-series to standardize the global analysis. Liu et al. (2023) compared different approaches for estimating phenological metrics. Their study demonstrated that, when compared to ground data, local thresholds are more accurate for estimating the start of the season (and the end of the season) using NDVI.

This threshold method has been previously used to estimate land surface phenology metrics using vegetation indices - VIs (such as NDVI or LAI) (White et al., 2009; Garonna et al., 2014; Karkauskaite et al., 2017; Liu et al., 2023; Li et al., 2021). Karkauskaite et al. (2017) established a local threshold of 20% to determine the start of the growing season in northern latitudes, which is equivalent to 0.2 of the normalized NDVI time series in our study. Garonna et al. (2014) used a 50% threshold, while Liu et al. (2023) used both 20% and 50%. Additionally, Li et al. (2021) used multiple thresholds, defining 75% as the moment of peak growing season.

This document employs then a local threshold of 0.75 to ensure the capture of the period of peak vegetation growth/development conditions.

Hence, a normalized SM value exceeding 0.75 represents the periods where superficial moisture levels have surpassed at least three-fourths of their annual range. Similarly, for the VOD, it indicates that water retention for a given plant component is nearing its relative maximum, and for NDVI, represents the greening-up of vegetation, representing full development and health. Periods when the normalized time series of each variable exceed this threshold will be referred to as the peaking period (PP) for that variable.

Next, considering our objective to understand the seasonal responses of both passive microwaves and optical indices to water availability, the onset of each PP will be used as the primary descriptor of VODs and NDVI reactivity to water availability, hereafter called the 'Start of' (SO) the peaking period. Thus, the day of the year (DOY) when the SM surpasses the 0.75 threshold is defined as SMSO, for LVOD it will be LVODSO, and so forth.

Finally, given the spatial scale of the study area, which encompasses different climatic regions, the DOY for the SMSO, and the other variables, is expected to vary significantly from one location to another. To mitigate the inherent variability in the SO dates, our analysis focuses on the relative timing of VOD's and NDVI's responses to water availability, rather than their absolute calendar dates.

Consequently, the relative lag in days between SM and the different plant descriptors is computed. These lags will be identified as 'Δ + Descriptor'. For instance, the lag between the LVODSO and the SMSO is computed as LVODSO - SMSO and is referred to as ΔLVOD, and similarly for the other variables. The PP for the SM is defined as the longest period, of consecutive days, during which the threshold of 0.75 is exceeded. For the remaining variables, the PP is determined as the longest period, with a consecutive number of days above the threshold, closest to that of the SM.
<figure><img src="image_121.png" /><figcaption>Figure 9.2: Global mean behavior of Δ values per cluster (N=11), lags are given in days. The SM<sub>SO</sub> is represented as a green vertical line in each box-plot. The order of box-plots has been deliberately modified so cluster with similar behaviors are adjacent to each other.</figcaption></figure>
### 9.3.3 Group Identification through K-Means Clustering

Given the extensive scale of the study area and the possible representation of each pixel by multiple attributes (e.g., ΔLVOD, ΔCVOD, ΔNDVI), a straightforward pixel-wise classification approach is not feasible. To extract the relationships and patterns within our dataset, an unsupervised machine learning classification algorithm was employed. Specifically, the K-means clustering algorithm was selected for this purpose. K-means clustering operates on the principle of partitioning an input dataset into 'K' distinct clusters based on shared similarities within the data. Initially, the algorithm assigns centroids randomly. Subsequently, each data point is linked to the nearest centroid, thereby forming the initial clusters. The centroids are then recalculated based on the mean values of the data within each cluster, and the distances of each point to its new centroid are recomputed. This iterative process continues until the positions of centroids exhibit insignificant changes.

Finally, due to the sensitivity of the k-means method to the initial set of centroids, 1000 draws were performed, and the one resulting in the smallest intra-cluster sum of squared distances (CSSD) was selected.

### 9.3.4 Attributes and cluster selection

The K-means clustering algorithm needs two primary inputs: first, a dataset comprising relevant attributes that describe the phenomenon under investigation, and second, a predetermined number of groups (N) in which the dataset will be divided.

The proposed attribute matrix includes six variables. The first four, which are the ΔVODs at different frequencies, describe the reactivity of water storage in various plant components to SM availability. The fifth variable, the ΔNDVI, provides insights into the reactivity of the leafy components of plants to water availability. Finally, to adequately account for external environmental factors that significantly influence both water availability and plant development, the mean temperature (T) is added as the sixth and final attribute.

On the other hand, selecting the optimal N usually depends on prior knowledge of the dataset. However, to our knowledge, no comprehensive classification addressing the hydraulic behaviors of crops at a global scale exists. Consequently, defining an initial number of possible classes is difficult, making any number of clusters theoretically acceptable. In this study, we chose to assess clustering quality by varying the number of N within a specific range.

The climate classification by Köppen-Geiger (Beck et al., 2023) divides the study area into four main climatic zones: A - Tropical, B - Arid, C - Temperate, and D - Continental. Here, we opted to establish the minimum N at 4 (to cover at least the possible climatic regions of the study area) and observe the quality of clustering with increasing values. With too few clusters (2 or 3), the model might lose information and fail to capture essential patterns in the data. Conversely, with too many clusters, the resulting classification may lose interpretability. Clustering quality was evaluated based on the mean Silhouette value (a measure of cluster cohesion) and the geographical distribution of clusters.

### 9.3.5 Dynamic time warping

In addition to the methods discussed in Sect. 9.3.2, a more automatic method was used to compare different time series for a selected number of points. The goal is also to infer the complementarity of VOD (derived from L, X and C bands) and vegetation phenology indices (derived from optical data).

We choose to perform time-series analysis using a statistical approach for analyzing data that has been structured through time and detect patterns, trends, and anomalies. In time-series analysis, there are numerous approaches for searching for similarities, including the Euclidean distance, intercorrelation computation, dynamic time warping (DTW), and shape-based methods like the Fourier transform and Symbolic Aggregate Approximation (SAX).

In this activity, the degree of similarity of different VOD and optical indices time series will be analyzed through three metrics:

* Dynamic Time Warping (DTW) distance
* Pearson correlation

* time shift between time series

The DTW distance is a non-linear, non-time dependent measure of the alignment – similarity – between two time series. Contrary to the classical methods for time series comparisons (i.e. Euclidean distance), which perform a point-to-point comparison, the DTW is able to perform a suitable matching of the time series even if they are out of phase in the time axis or have different lengths (Jeong et al., 2011) and thus could be used to identify anomalies in time series. Pearson correlation and time shift estimation will bring additional information to the analysis about the correlation between time series and their similarities. The time shift estimation was performed via the cross-correlation function which is a measure of similarity of two series as a function of the displacement of one relative to the other. Thus, finding the maximum of the cross-correlation function and its corresponding displacement value means finding the best overlay between the time series and estimating the shift between them.

We describe the comparative analysis of the time series produced, for several locations on the earth's surface, from microwave VOD data (SMOS, AMSRE, AMSR2) and optical AVHRR LAI data variable available in the Datacube #2 archive. In our analysis we were targeting the LANDVAL in-situ sites that correspond to the CROPLAND class from the IGBP classification. The list of such sites was determined beforehand by cross-referencing the list of the LANDVAL sites and the global classification map provided by the IGBP.

To complete this study we performed following steps :

1. Create VOD time series from SMOS database for the filtered CROPLAND sites

2. Analyse obtained time series visually in order to select sites with consistent and continuous data only, sites with few data points or with no apparent patterns were set aside for a latter study

3. Process VOD time series with moving average method to smooth the data and to highlight the underlying pattern

4. Create VOD and AVHRR LAI time series for the final filtered CROPLAND sites

5. Overlay time series for each station, analyze it visually and perform statistical parameters computation

First, in order to filter down the list of CROPLAND sites and to select stations with only viable and consistent data, visual analysis of the time series was performed.

After the visual analysis, the list was filtered down from 120 to 27 CROPLANDs stations. On those stations a pre-processing step was done to prepare the data for the latter comparisons and statistics computation. To smooth the data and to highlight the underlying patterns in the time series, the moving average method was used. The smoothing step allows to make time series cleaner which then allows a more precise and coherent statistics computation when comparing time series between them. This step was then performed on all of the selected stations and for all of the VOD time series (SMOS L, AMSRE C/X, AMSR2 C/X). After performing the smoothing step, we can also plot time series from different satellites together to analyse their similarities and disparities if exist.

After the selection of stations and pre-processing of the VOD data, time series can be plotted with the AVHRR LAI variable for the comparison and statistics computation. To analyse similarities between time series in a more statistical way with quantitative metrics, we used metrics such as Dynamic Time Warping

## 9.4 Results

### 9.4.1 Pixel Seasonality from pulse diagrams

The seasonality of the selected variables was systematically analyzed across global croplands. A discernible seasonal pattern emerges across all variables, characterized by a peak occurring approximately between June and September for vegetation descriptors (VODs and NDVI); the temperature time-series on the contrary
<figure><img src="image_122.png" /><figcaption>Figure 9.3: Pulse diagrams for single cropland-pixels in USA, Argentina, Spain, Sudan, Morocco and Australia. Mean pulses computed from observations between 2011 and 2020. Red lines correspond to the PP for each variable.</figcaption></figure>

reaches its peak in the first months of the year. The first variable to ascend and reach its maximum after the SM is the NDVI, followed by the VODs. Overall, all frequencies of VOD exhibit synchrony throughout the year.

While figures of this nature offer a comprehensive depiction of the overall behavior of the selected variables, the inherent disparities in scale, ranges, and units prevent a more intuitive interpretation. As described in subsection 9.3.1, in this study we have chosen to represent these seasonal patterns as relative pulses (normalized time-series). Values closer to 1 indicate proximity to the annual maximum, while values close to 0 represent proximity to the minimum. In figure 9.3 (Australia) the size and color of each pulse represent the magnitude of the normalized value: brighter and larger dots indicate values closer to 1, while darker and smaller dots indicate the opposite. The PP of each variable are highlighted in red.

The new time series representation offers enhanced interpretability and facilitates graphical comparison across various geographic locations. Pulse diagrams for different locations are depicted in figure 9.3. Notably, consistent patterns emerge across different regions, as clearly observed in Sudan, Morocco and Australia figure 9.3.

In these locations, the dynamics of LVOD<sub>SO</sub> and NDVI<sub>SO</sub> closely track variations in SM<sub>SO</sub>. However, in Morocco, LVOD starts its peak development nearly synchronously with NDVI, whereas in Sudan, the NDVI<sub>SO</sub> occurs first. In both regions, the occurrence of higher frequency VODs lags behind SM by approximately 45 days. A comparable trend to that observed in Morocco is discernible in Australia, albeit with a delay of approximately six months. This temporal shift can be attributed to seasonal variations linked to the differing latitudinal positions of the two regions.

At northern latitudes (USA, Spain), the pattern changes drastically. Vegetation descriptors tend to lag several months behind SM<sub>SO</sub>. While clear similarities/differences in behavior can be observed in different regions, performing a pixel-by-pixel classification is impractical due to the data volume and complex, unique patterns
of individual pixels. Therefore, the utility of the classification algorithm lies in its capacity to group similar patterns efficiently without manual intervention.

## 9.4.2 Attribute selection

As mentioned in Section 9.3.4, the attributes selected to describe VODs and NDVI reactivity to SM include the delta values (the lag of each variable SO with respect to SM<sub>SO</sub>) and T. However, when examining the seasonal signal of higher frequency VODs across multiple locations, including but not limited to those shown in figure 9.3, a strong similarity was observed. The time series of these VODs demonstrated coordinated increases and decreases, indicating synchronized seasonal patterns.

Additionally, significant correlations (as R coefficients) between ΔCVOD, ΔXVOD, and ΔKVOD were observed at a global scale; ranging from 0.84 to 0.91. Consequently, only one of the higher frequency VODs was retained. Given that the C-Band provides information on medium-sized components of vegetation, and considering that larger components are already represented by LVOD and leafy vegetation by the NDVI, ΔCVOD was retained as representative of the higher frequency VODs.

Finally, to simplify the selected features and address the remaining correlations between ΔNDVI and ΔCVOD, Principal Component Analysis (PCA) was employed to reduce the dimensionality of the attributes, namely ΔLVOD, ΔCVOD, ΔNDVI, and T. The first principal component (PCA1) accounts for around 54% of the total variance in the dataset, followed by 25% for PCA2 and 14% for PCA3. Therefore, the first three principal components were selected as the final attributes for the clustering algorithm, collectively explaining 93% of the variance in the original dataset.

## 9.4.3 Unsupervised classification

The geographical distribution of the produced clusters (N=11) is shown in figure 9.2. The lower part of the figure presents the mean Δ values for the VODs and the NDVI per cluster, with the green vertical line representing SM<sub>SO</sub>. The order of clusters was deliberately adjusted so that clusters with similar Δ behaviors are adjacent.

The first group of clusters comprises K4 (yellow), K7 (golden), and K8 (blue). As observed in subsection ??, increasing the number of N during the classification process has little effect in their geographical distribution suggesting distinctive patterns and clear differences in their behaviors. Elements in this group are characterized by a distinct pattern: ΔLVOD tends to behave almost in synchrony with ΔNDVI, while higher frequency ΔVODs are consistently larger than the latter two (see corresponding box-plots in figure 9.2).

Both K4 and K7 are spatially close to each other, located in Australia, Africa, the southern Iberian Peninsula, and tropical America. As presented in Subsection ??, K4 is a branch of K7 that appears when N reaches 5. Among all eleven clusters, K4 and K7 present the lowest Δ values across all variables. For K4, NDVI<sub>SO</sub> is almost in synchrony with SM<sub>SO</sub>, with a mean ΔNDVI close to zero (-2±42 days). LVOD<sub>SO</sub> starts later in time with a ΔLVOD of approximately 19±42 days. In K7, the same pattern is observed, but both ΔLVOD and ΔNDVI increase by about 23 days (42±41 and 25±28 days, respectively). ΔCVOD is approximately 42±50 and 66±33 for K4 and K7, respectively.

K8 on the other hand is located at northern latitudes, mostly in North America, with lesser extents in Europe and Asia. Overall, K8 increases its Δ values to 54±30, 89±28, and 63±32 for LVOD, CVOD, and NDVI, respectively. It presents three main differences compared to K4 and K7: its geographical position (further north), an increase in Δ values of roughly 20 days (relative to K7), and the tendency of LVOD<sub>SO</sub> to be closer to SM<sub>SO</sub> than NDVI<sub>SO</sub>, which, as observed in figure 9.2, is not often the case for other clusters.

The second group corresponds to clusters K3 (green), K5 (magenta), and K9 (brown). All the clusters in this group can be considered subdivisions of the main four clusters at N=4. K3 is the most stable of the three and is the first to appear (at N=5) as a subdivision of the original blue and orange clusters (see N=4 in figure ??). K5 and K9, on the other hand, are the last clusters to be created, at N=11 and N=10, respectively. K5 is primarily a branch of both the blue and orange clusters at N=4, while K9 is a subdivision of the golden cluster at N=4.
All Δ values in the second group follow a stair-like pattern, with ΔNDVI at the top (56±24, 62±46, and 26±65 days for K3, K5, and K9, respectively), followed by ΔCVOD (89±27, 86±45, and 90±53 days), and finally ΔLVOD (121±26, 234±38, and 156±54 days). Though their Δ values for NDVI and CVOD are similar in magnitude, the main source of their differences lies in their geographical locations and the ΔLVOD behavior (which ranges from 121 to 234 days). K3 and K5 are both located at northern latitudes, while K9 is in the pan-tropical region.

Elements in the third group, comprising K6 (crimson) and K11 (pink), are derivations of the orange cluster at N=4 (see figure ??). These clusters are primarily located in western Europe and are close to each other, though K11 is situated further north. Two main attributes define the elements in G3: first, the Δ values for NDVI and CVOD surpass 120 days (which is not observed in the previously mentioned clusters), and second, the ΔLVOD can either start at the same time as SM<sub>SO</sub> (ΔLVOD = 0±57 days) or much later (245±34 days). Specifically, ΔNDVI values are 150±51 days for K6 and 140±31 days for K11, while ΔCVOD values are 170±56 days and 171±32 days, respectively.

The final group corresponds to clusters K1, K2 and K10. These clusters demonstrated stability with increasing N and are part of the original four clusters observed when N=4; except for K1 that appears when N reaches 7. Geographically, they are located in America and Northern China, although K10 can have some elements in Europe. In North America, these clusters are adjacent to each other, whereas in South America, K2 dominates the region, with K1 and K10 being less numerous. The main feature of these clusters is the short lag between the Δ values across all attributes, resulting in overall synchrony of all variables. The main differentiating factor between K2 and K1 with respect to K10 is the displacement of all Δ values with respect to SM<sub>SO</sub>: around 100 days in K10 and over 200 days in K2. K1 and K2 exhibit a similar behavior, differing only in lag direction. The negative lag in K1 could be considered equivalent to the positive lag in K2 due to the time series' cyclical nature.

### 9.4.4 Assessment of Unsupervised classification using external data

We used the crop mask by Teluguntla et al. (2016) and the Köppen-Geiger climate classification (Beck et al., 2023), both introduced in Section 9.2, to cross-validate our results. The geographical distribution of classes from these products aligns closely, likely due to the dependence of agricultural practices on regional climate conditions.

In North America, particularly between Canada and the USA, K8 (blue) and K10 (orange) are the dominant clusters. K8 predominantly represents wheat-barley crops [RA: Wh.Ba.(D)], while K10 represents corn-soy crops [RA: Co.So.(D)]. Both clusters share a similar continental (cold) climate without a dry season, but K8 experiences slightly colder summers than K10.

K8 exhibits an 80% agreement with both the GFSAD1KCD and Köppen-Geiger maps, while K10 shows around 70% agreement. These percentages are computed from the fraction of pixels of each cluster that fall into the crop/climate classifications.

In K8, LVOD is in synchrony with changes NDVI, whereas in K10, it is aligned with CVOD. The water column (represented by the VODs) and the greening of plants (represented by the NDVI) presented a quicker reaction to water availability in K8, with delta values ranging between 50 and 90 days, compared to over 100 days in K10.

Another major cluster in America is K2 (cyan), primarily located in the Mississippi River Basin and the Argentinian Pampas. It shows approximately 83% similarity with the climatic class [Cfa], characterized as a humid subtropical climate. In northern regions, the dominant crop is irrigated wheat-rice [IR: Wh.Ri.(D)], with about 41% agreement, while in the south, it consists mostly of a mix of rain-fed crops [RA: Mix3], accounting for 53%.

In these regions both the VODs and the NDVI seasonal behavior was in synchrony. The CVOD starts its season about 15 days later than the NDVI, this difference increases to 27 days with the LVOD. K2 exhibits the largest Δ values among all clusters (above 180 days).

The fourth main cluster is K3 (green), which appears in North America and Europe, with a small fraction in northern China. It represents areas with mixed crops [RA: Mix3] in continental climates, with warm summers
and no dry season [Dfb]. This cluster covers the transitional regions where the LVOD has a high reactivity to water availability to areas where its response tends to be slower.

K7 (golden) is the predominant cluster in southern latitudes, primarily found in Australia and Africa, with a small presence in the southern Iberian Peninsula (Figure 9.2). In Australia and Africa, it corresponds mainly to wheat-barley crops [RA: Wh.Ba.(D)], with approximately 69% and 55% agreement respectively, while in Europe, it is mainly related to mixed crops [RA: Mix3]. K7's climates feature hot, dry summers in Europe and Africa, and cold semi-arid weather in Australia. Its behavior is similar to K8, with LVOD nearly synchronous to NDVI (followed by CVOD). Overall, the Δ values are shorter than in K8, falling below 50 days (for LVOD and NDVI), indicating higher reactivity in the water column to SM variations.

The sixth of the main clusters is K4 (yellow), which is characteristic of warm weather and primarily located in the pan-tropical region. Approximately half of its pixels are associated with regions dominated by Wheat-Rice-Soy [RA: Wh.Ri.So.(D)], and the other half with areas having a low fraction of croplands [MinorFrac]. K4's behavior is similar to K7 and K8, but it exhibits the lowest Δ values among all clusters, indicating higher sensitivity in the water column to water availability.

Europe and Northern Asia, on the other hand, show the most significant changes when the number of clusters is increased, encompassing seven out of the eleven possible clusters (as shown in figure 9.2). This suggests a high degree of variability and complexity in the attributes used for the classification in these regions.

Some clusters in Europe and Asia exhibit behaviors close to those of the other regions, such as K5, K9, and K11 being similar to K3, possibly indicating areas of transition characterized by highly heterogeneous crops. However, other clusters present clear distinctive patterns. For example, in K6, LVOD is strongly synchronous with SM, but NDVI is not.

### 9.4.5 Dynamic time warping

The computed metrics of the CROPLANDs stations are presented on a latitude/longitude map on the figure 9.4. The distance and shift values are the median values computed for each station across the five VOD datasets (SMOS L band, AMSRE C and X bands, AMSR2 C and X bands), the AVHRR LAI dataset being the one to which we compare the VOD time series.

On this map we can see that the majority of the stations present same order of magnitude in their respective metrics. The distance values are around 20-40 in average, and the shifts are up to a maximum of 50 days in most of the cases. No apparent spatial variability was found analyzing given points.

We can also overlay these results to the CROPLANDs clusters that were found in Sect. 9.4.3 in order to see of our spatial variability matches cluster distributions. We can see that for some classes we observe similar metrics for the stations inside, such as BELMANIP_00089 and DIRECT_00001_AGRO stations from the class n°10. However we also have sites in the Australian region that have different metrics while being in the same n°7 cluster.

## 9.5 Discussion

In this study, we aimed to understand the seasonal behavior of PMW optical depth and NDVI in response to water availability over croplands. This analysis resulted in a well-structured geographical classification of vegetation based on the sensitivity of these variables to SM variations, revealing the potential hydric behavior of plants in different cropland areas.

### 9.5.1 Cluster Interpretation

As shown in Section 9.4, some of the clusters exhibit a clear agreement with the geographical distribution of the different classes from crop type and the Köppen-Geiger maps. Thus, it seems possible to relate the joint response of PMW and NDVI to SM as a function of agricultural practices and climate conditions. The classification found in Section 9.4.3, explained in terms of time lags, clearly depicts well-known regions.
<figure><img src="image_123.png" /><figcaption>Figure 9.4: CROPLANDs stations with their metrics overlaid on clusters in North America.</figcaption></figure>

For example, the cluster K10 presented strong spatial agreement with the USA's corn belt, as observed in figure 9.2. This region is characterized by strong synchronization among all VODs, which consistently lag behind the NDVI.

Furthermore, cluster K8 predominantly spans North America's wheat belt, while K7 encompasses the southern coast of Australia, the Mediterranean region, South Africa, and the majority of croplands in Ethiopia. Both clusters exhibit a distinctive synchrony between LVOD and NDVI (as shown in figure 9.2), and are predominantly found in regions where wheat is a major crop. A similar behavior is observed in K4, situated within the semi-arid Sahel belt; although it contains wheat crops, it also presents a significant fraction of non-agricultural low vegetation. We consider, therefore, that K4 and K7 represent warm arid ecosystems in Africa and Australia, predominantly featuring wheat crops. In contrast, K8 corresponds to cold continental regions in North America, which are also associated with wheat cultivation (USA's wheat belt).

Interesting behaviors were observed in other clusters, particularly K2 and K1, where the VODs and NDVI reached their peak development with a significant temporal offset with respect to SM's. The K2 cluster is found in three distinct regions: the Mississippi River basin (part of the USA's rice belt), the Argentinian pampas and China's Northern rice region. Regions like Europe, known for having a mosaic of different crop types, indeed showed the widest variety of clusters, underlying the high degree of variability in the attributes used for the clustering over this region.

While the clustering identified several classes, only the six mentioned before corresponded to well-defined geographical regions linked to specific cropland types. Therefore, in the next section, a potential interpretation of the response of PMW VODs and NDVI to water availability in those regions is given. Other regions exhibiting interesting patterns, like those represented by the remaining clusters require further research.

## 9.5.2 Temporal lags between PMWs and spectral indices

Previous studies on PMW applications for crops have mainly focused on the sensitivity of microwaves to plant phenology within the United States (Jones et al., 2012), particularly in the region known as the USA's corn-belt (Patton and Hornbuckle, 2013; Lawrence et al., 2014b; Hornbuckle et al., 2016; Guan et al., 2017). A systematic lag between the peaks of LAI and the LVOD, with LAI appearing around 19 days earlier, was observed by (Lawrence et al., 2014b) in some locations of this region. The mean seasonal behavior of the
XVOD over the corn-belt was studied by (Guan et al., 2017), they observed that the NDVI started its growing season before the XVOD by around 1 month. Our results align with (Guan et al., 2017) regarding the lags between XVOD and NDVI. The observed lags were around 25 days with a comparable response in LVOD of approximately 28 days. The observed lags are slightly longer than those reported by Lawrence et al. (2014b) by around 6 days.

Additionally, Vreugdenhil et al. (2017), using active CVOD over Australia, observed that LAI reached its growing season 1 to 2 months before CVOD over croplands; here we observed a similar behavior and estimated it at around 37 days. Field experiments with L-Band radiometers in Europe estimated that the NDVI appeared around two months before the LVOD over a mixed parcel in eastern France, where approximately 55% of the surface was identified as croplands (Barrée et al., 2021). The absence of data specifically corresponding to that region in our study prevents direct validation of these observations.

Other research has observed this lag between optical and PMW indices across different cover types, besides croplands (Jones et al., 2011, 2012). Jones et al. (2011) observed that in regions with low vegetation, the KVOD lags behind various VIs, although no estimation of actual time lags was provided. A similar approach was adopted by Jones et al. in a subsequent study (Jones et al., 2012), but focused exclusively on North America. In their study, across various cover types, they observed that in southern latitudes, the XVOD reached its growing season, after the NDVI by approximately 5.7 weeks (around 40 days). We estimated this temporal lag to be approximately 26 days. They also noted that in northern latitudes of North America, the XVOD appeared before the NDVI by about 7.5 weeks (around 53 days). Our observations, however, revealed contrasting patterns over croplands in this region: the NDVI peaked about 26 days before the high-frequency VODs. Possible explanations for the differences between our findings and those of Jones et al. include the methods used to estimate the start of the season and the exclusive focus on croplands in our study, as opposed to various cover types in theirs.

In this study, these observations were extended to croplands around the world, not just in North America or Australia, and to multiple frequencies. A systematic lag between NDVI and higher frequencies VOD in all clusters was observed, ranging from 15 to 49 days. Furthermore, the current analysis revealed that, over croplands, the LVOD does not consistently lag behind the NDVI. Unlike higher frequency VODs, LVOD exhibits variable timing in reaching its peak development relative to NDVI.

### 9.5.3 Physical interpretation of the vegetation parameters

The delayed response of high-frequency VODs compared to NDVI, can be attributed to the sensitivities of both variables to different aspects of vegetation. NDVI responds to leaf greening, while the VOD is sensitive to the water stored in the different components of vegetation. Across the entire study area, photosynthetic activity (NDVI) reaches its peak development initially, followed by water accumulation in the finer components of vegetation (high-frequency VODs) later in the growth cycle. This behavior is similar to that found in field measurements of C-Band radar over cornfields, which showed that LAI increases sooner and faster than canopy water content (Wigneron et al., 1999).

In contrast, the timing of LVOD relative to both NDVI and high-frequency VODs appears to vary with crop type. In the corn belt of the United States, LVOD reaches its peak growth in synchrony with high-frequency VODs, and both occur after the NDVI. This pattern suggests that water distribution within the plant's components—such as leaves, branches, and stem—remains relatively uniform throughout its development. Additionally, both NDVI and VODs reach their peak approximately 120 days after SM does. This delay between water availability and peak vegetation development indicates a slower rate of vegetative growth following the onset of water availability. In contrast, in regions where wheat is the dominant crop, such as Africa and Australia, VODs and NDVI respond more rapidly to water availability, peaking within 66 days after the SM signal. Unlike areas where corn is the primary crop, the LVOD in these regions closely follows the NDVI and appears before the high-frequency VODs. It appears that vegetation prioritizes photosynthetic activity first, followed by water storage in larger components, and only then stores water in smaller components.

It has been shown (Nguyen et al., 2022) that under water stress conditions, corn exhibits isohydric behavior (constant water potential), maintaining stable leaf water potential (LWP) values. In contrast, winter wheat displays anisohydric characteristics (variable water potential), showing more negative LWP values during water stress episodes. Moreover, a comparative analysis of these two crops (Wang et al., 2012) showed a disparity
in transpiration rates, with corn fields exhibiting consistently higher levels than winter wheat. Wang et al. (2012) observed that the transpiration in summer corn accounted for ≈71% of total evapotranspiration, while for winter wheat, it accounted for ≈61%. The stable LWP and higher transpiration rate suggest that corn is actively managing its water use to maintain more consistent water levels throughout its tissues during growth, which may explain the observed synchrony in PMW VODs. Wheat, on the other hand, has been shown to use water more efficiently (Wang et al., 2012) and exhibit distinct physiological strategies, as observed in the adjustment of LWP to water availability by Nguyen et al. (2022). The efficient use of available water could explain the faster reaction time to SM pulses. Changes of LWP and smaller transpiration rates may explain why LVOD peaks before high-frequency VODs. Vegetation in these regions seems to prioritize maintaining hydration in essential tissues, particularly larger structural components (indicated by LVOD). Only after, plants allocate water to smaller, more dynamic parts (captured by high-frequency VODs) under fluctuating water availability conditions.

Nonetheless, we also noted that LVOD precedes the NDVI in the USA's wheat-belt and in the Iberian peninsula. This behavior has been observed in other types of vegetation and for other VOD frequencies relative to NDVI. For instance, Jones et al. (2011) found that XVOD initiated its growing season before NDVI in vegetation over northern latitudes of North America. They proposed that cold temperatures inhibit vegetation's photosynthetic activity, explaining the temporal offset between VOD and NDVI signals. Their hypothesis suggests that plants may begin water uptake in non-photosynthetic tissue (represented by LVOD or XVOD) before temperatures are warm enough to initiate significant chlorophyll production, thus accounting for the delayed appearance of NDVI relative to VOD. This behavior could apply to perennial crops, where woody structures such as trunks and branches are maintained between growing seasons. However, it is less likely for annual crops, which perish and regrow within the same season and are the focus of our study. Moreover, Tong et al. (2019) observed a similar behavior (XVOD preceding NDVI) in Northern Australian shrublands, a region characterized by warm weather. This observation suggests that low temperatures may not be the sole driver behind this behavior. Further research is needed in regions where this particular pattern is observed.

Furthermore, we observed that in the USA's rice-belt, the Argentinian pampas and China's northern rice region, there is a considerable temporal offset in the seasonal signal of VODs and NDVI with respect to SM. Bousquet et al. (2021) found negative correlations between LVOD and SM in seasonally inundated areas, such as wetlands and floodplains. They suggested that this behavior is more pronounced in grasslands than in forests because low vegetation is often submerged and becomes invisible to sensors, while standing water becomes more prominent in the signal. This could explain the substantial Δ values observed in these regions which are prone to seasonal inundation, including potential areas of paddy cultivation. This is particularly evident in the Mississippi basin, where Bousquet et al. (2021) identified a highly dynamic water fraction. Although less prominent, their results also showed that the Argentinian Pampas and northern China exhibited relatively high dynamic water fractions compared to other regions. Similar phase opposition between VOD and VIs, such as the LAI and the Gross Primary Production (GPP), has been documented in tropical woodlands (Tian et al., 2018), evergreen broadleaf forests, regions with frequent inundations (Jones et al., 2011), and over wetlands (Teubner et al., 2018).

Moreover, a link between weather conditions and the response time of the different variables to SM was observed. Regions in Australia, Africa, and North America where the main crop is considered to be wheat exhibited shorter reaction times to SM pulses when located in drier and warmer regions. For instance, in Australia and Africa, VODs and NDVI are more responsive to water availability compared to North America. Specifically, LVOD reacts up to 35 days faster, NDVI up to 60 days faster, and CVOD up to 47 days faster than in North America. In crops near the Sahel region of Africa, which is known for its dry and warm conditions, the fastest reaction of VODs and NDVI to SM was observed ; the lag for ΔLVOD was approximately 19 days, while ΔNDVI was nearly synchronous with SM. This suggests that the water column stored in crops (captured by PWM VODs) in these regions is more sensitive to water availability, and that this sensitivity seems to be driven by harsher climatic and environmental conditions. This hypothesis is supported by studies like Huxman et al. (2004) which showed that photosynthetic activity in herbaceous vegetation in arid and semi-arid ecosystems responds rapidly and significantly to water availability.

More studies are needed to understand the PMW signals over vegetated regions with particular behaviors, like those where LVOD is strongly synchronized with SM but not with other VIs. As shown by Liu et al. (2021), classifying vegetation by hydric properties, such as those inferred from VOD, including vegetation water content or LWP, could enhance the parameterization of plant hydraulics in large-scale land surface
models.

### 9.5.4 Comparing to dynamic time warping

The dynamic time warping analysis was done for a small number of points and it is not possible to define a particular spatial pattern or apparent clusters. The comparison of DTW results with the pulses analysis cannot be done systematically for the same reason. The results are sometimes similar to those of the time pulses analysis. For instance for the points shown in figure 9.4 BELMANIP_00088 belong to cluster 3 (green, Table ??) for which the LVOD-NDVI lag is around 68 days on average, while the DTW analysis for this point suggest no shift. Regarding the two sites in the orange cluster (cluster 10), they also suggest no shift while on average the orange cluster has a LVOD-NDVI lag of 27 days. In contrast, the shift computed with DTW for the points in clusters 7 (gold) and 2 (cyan) are 10.35 and 26.7 days, respectively. This values are in better agreement with those of the pulse analysis which are 15 and 27 days, respectively.

## 9.6 Conclusion

By integrating data from passive microwave and optical sensors, this study presents a geographical classification of vegetation over croplands based on the seasonal responses of multi-frequency VODs and NDVI to SM pulses. Results from the clustering algorithm align with known agricultural practices and climatic classifications, highlighting the capability to capture significant ecological and agricultural patterns in vegetation. The analysis identified eleven distinct classes, of which six could be related to distinct (climate) regions and crop types.

Globally, we identified a systematic lag between higher frequency VODs and NDVI, ranging from 15 to 49 days. In contrast, the seasonality of the LVOD with respect to high-frequency VODs and NDVI varied with location. In wheat-dominated areas, the vegetation's water column and photosynthetic activity were highly reactive to SM pulses. We observed that the water content in vegetation's larger structural components, as indicated by LVOD, exhibited near synchrony with NDVI and a faster response to water availability compared to that in smaller, more dynamic plant parts, which are better captured by high-frequency VODs. Conversely, in corn-dominated regions, such as the USA's corn-belt, the PMW VODs and optical indices exhibited a more gradual response to water availability. In these areas, the LVOD was in synchrony with the high-frequency VODs and both lagged behind the NDVI. We hypothesize that this multi-frequency analysis reveals the isohydric nature of regions predominantly cultivating corn and the anisohydric characteristics of areas where wheat is the primary crop.

Additionally, the classification effectively isolated regions with particular seasonal behaviors, like the Mississippi River basin, the Argentinian pampas and China's northern rice region, where PMW VODs and NDVI lagged significantly with respect to SM. This behavior seems to be related to seasonal changes in the water fraction. Other locations, presenting particular patterns such as LVOD in synchrony with SM but not with NDVI/CVOD require additional investigation. This underscores the critical importance of studying PWM in vegetated regions to understand how surface conditions affect the signal.

These findings highlight the benefits of integrating passive microwave and optical indices to better understand the complex interactions between vegetation and water availability. However, the spatial resolution of the data, while suitable for large-scale analysis, may overlook finer-scale variations within croplands. The clear patterns observed in certain regions may be attributed to spatial resolution and extensive monocultural practices. In fragmented landscapes like Europe, the clustering exhibited reduced efficacy in identifying distinct patterns. Subsequent research could incorporate field measurements to validate and refine our classifications.
# Chapter 10

## Deliverables, progress reports, meetings

### 10.1 Deliverables and reports

Table 10.1 presents the OSMOSE project deliverables, and Table 10.2 gives the progress reports that were prepared regularly.

<table><caption>Table 10.1: OSMOSE deliverables table</caption><thead><tr><td>Task #</td><td>Deliverable</td><td>Long Name</td><td>Type</td><td>Reference</td></tr></thead><tbody><tr><td>T0</td><td>D1</td><td>Web Site</td><td>web</td><td>https://osmose.magellium.com/</td></tr><tr><td>T1</td><td>D2</td><td>Theoretical Analyses Microwave Radiative Transfer</td><td>TN 1</td><td>TN_OSMOSE_WP1_22Jun23.pdf</td></tr><tr><td>T2</td><td>D3</td><td>ATBD microwaves</td><td>TN 2</td><td>OSMOSE_D03TN02_ATBD_VODchain_20240109.pdf</td></tr><tr><td>T2</td><td>D4</td><td>Processing Algorithms for the generation of the VOD</td><td>Algo.</td><td></td></tr><tr><td>T3</td><td>D5</td><td>Data base for the generation of VOD</td><td>Data</td><td></td></tr><tr><td>T3</td><td>D6</td><td>Data Base Manual</td><td>TN 3</td><td>OSMOSE-NT-025-MAG-D6<br>OSMOSE_D06TN03_Manual DataBase_20240118</td></tr><tr><td>T4</td><td>D7</td><td>Output Database VOD</td><td>Data</td><td></td></tr><tr><td>T4</td><td>D8</td><td>VOD data set documentation and analysis report</td><td>TN 4</td><td>OSMOSE_D08TN04_VOD_Analysis_20240402.pdf<br>(ref. SO-TN-CB-GS-0118)</td></tr><tr><td>T5</td><td>D9</td><td>AGB data documentation and analyses report</td><td>TN 5</td><td>OSMOSE_WP5_report.pdf<br>(ref. SO-TN-CB-GS-0117)</td></tr><tr><td>T5</td><td>D10</td><td>AGB estimates from VOD and TB</td><td>Data</td><td></td></tr><tr><td>T6</td><td>D11</td><td>Optical data</td><td>TN 6</td><td>Review_VegTraits_Estimation.pdf</td></tr><tr><td>T7</td><td>D12</td><td>Pre processing ATBD Optical</td><td>TN 7</td><td>OSMOSE_ATBD_VWC_Optical_V3.0<br>OSMOSE-NT-007-MAG-v1.0_D12</td></tr><tr><td>T7</td><td>D13</td><td>ATBD Second Optical Data Cube:<br>- description<br>- validation manual</td><td>TN8</td><td>OSMOSE-NT-024-MAG-D13.pdf<br>OSMOSE-NT-006-MAG-D13.pdf<br>OSMOSE-NT-037-MAG.pdf</td></tr><tr><td>T7</td><td>D14</td><td>Pre processing algorithm</td><td>Algo.</td><td></td></tr><tr><td>T8</td><td>D15</td><td>Input data base for optical</td><td>Data</td><td>OSMOSE-NT-008-MAG-v1.0_D15</td></tr><tr><td>T8</td><td>D16</td><td>Update technical note 3</td><td>TN 3</td><td>OSMOSE-NT-025-MAG-D16_v1.0</td></tr><tr><td>T9</td><td>D17</td><td>Output database optical</td><td>Data</td><td></td></tr><tr><td>T9</td><td>D18</td><td>Optical Vegetation Data manual</td><td>TN 9</td><td>OSMOSE-NT-027-MAG-v1.0_D18</td></tr><tr><td>T10</td><td>D19</td><td>Data analysis report</td><td>TN 10</td><td>SO-TN-CB-GS-0125</td></tr></tbody></table>
<table><caption>Table 10.1: OSMOSE deliverables table</caption><thead><tr><td>Task #</td><td>Deliverable</td><td>Long Name</td><td>Type</td><td>Reference</td></tr></thead><tbody><tr><td>T11</td><td>D20</td><td>Road map</td><td>TN 11</td><td>SO-TN-CB-GS-0121</td></tr><tr><td></td><td>FR</td><td>Final report</td><td>TN</td><td>SO-TN-CB-GS-0119</td></tr><tr><td></td><td>SR</td><td>Summary report</td><td>TN</td><td>SO-TN-CB-GS-0120</td></tr><tr><td></td><td>FP</td><td>Final presentation</td><td>TN</td><td></td></tr></tbody></table>

<table><caption>Table 10.2: OSMOSE Progress Report table</caption><thead><tr><td>Number</td><td>Period</td><td>Long Name</td></tr></thead><tbody><tr><td>01</td><td>08/09/2022-15/11/2022</td><td>OSMOSE-PR-013-MAG-V1.0-ProgressReport_20221118_v1.docx.pdf</td></tr><tr><td>02</td><td>06/12/2022-24/01/2023</td><td>OSMOSE-PR-019-MAG-V1.0-ProgressReport_20230120_v1.docx</td></tr><tr><td>03</td><td>18/02/2023-20/04/2023</td><td>OSMOSE-PR-026-MAG-V1.0-ProgressReport_20230417_v1.pdf</td></tr><tr><td>04</td><td>06/10/2023-06/12/2023</td><td>OSMOSE-PR-038-MAG-V1.0-ProgressReport_20231214_v1.1.pdf</td></tr><tr><td>05</td><td>10/01/2024-05/03/2024</td><td>OSMOSE-PR-045-MAG-V1.0-ProgressReport_20240316_v1.pdf</td></tr><tr><td>06</td><td>09/04/2024-20/06/2024</td><td>OSMOSE-PR-049-MAG-V1.0-ProgressReport_20240618.pdf</td></tr></tbody></table>

## 10.2 Project meetings

The OSMOSE team met once per month to coordinate the activities and to share information, and progress meeting with ESA were organized regularly (see Table 10.3). These meetings took place by videoconference, expect the MTR which was hold in Toulouse at CESBIO, and the Final Review also in Toulouse at CESBIO.

<table><caption>Table 10.3: OSMOSE meetings and milestones table</caption><thead><tr><td>Meeting/Milestone</td><td>Long Meeting Name</td><td>Location</td><td>Date</td></tr></thead><tbody><tr><td>KOM</td><td>Kick-off Meeting</td><td>Teleconference</td><td>09-06-2022</td></tr><tr><td>PM1</td><td>Progress Meeting 1</td><td>Teleconference</td><td>08-09-2022</td></tr><tr><td>PM2</td><td>Progress Meeting 2</td><td>Teleconference</td><td>06-12-2022</td></tr><tr><td>PM3</td><td>Progress Meeting 3</td><td>Teleconference</td><td>13-02-2023</td></tr><tr><td>PM4</td><td>Progress Meeting 4</td><td>Teleconference</td><td>06-07-2023</td></tr><tr><td>MTR</td><td>Mid-Term Review</td><td>CESBIO Toulouse</td><td>05-10-2023</td></tr><tr><td>PM6</td><td>Progress Meeting 6</td><td>Teleconference</td><td>09-01-2024</td></tr><tr><td>PM7</td><td>Progress Meeting 7</td><td>Teleconference</td><td>09-06-2024</td></tr><tr><td>FR</td><td>Final Review</td><td>CESBIO, Toulouse</td><td>03-10-2024</td></tr></tbody></table>
# Chapter 11

## Communication, outreach, and promotion

### 11.1 Website

A dedicated website was set up to share information about the OSMOSE project as well as to make all the project documents available for the team and forESA. The link is the following one:
https://osmose.magellium.com/
<figure><img src="image_124.png" /></figure>

## 10 years of SMOS – Passive Microwave Vegetation Opacity Study

The overall objective of this project is to create two consistent data cubes containing long, consistent, times series of vegetation variables estimated with microwaves and visible observations. The microwave Vegetation Optical Depth (VOD) will be derived from at frequencies from 1.4GHz to 37 GHz. Vegetation variables such as the Fraction of Absorbed Photosynthetically Active Radiation or Leaf Chlorophyll Content will be derived from visible wavelengths. The VOD and visible data will be used to explore a number of applications, in particular to estimate the Above Ground Biomass (AGB). An alternative method using machine learning approaches to estimate AGB directly from the brightness temperatures measured at microwave wavelengths, without estimating VOD, will also be studied.

The two methods to estimate AGB will be evaluated to i) assess the most efficient and appropriate method to estimate AGB ii) improve current AGB estimates from remote sensing observations.

SMOS is the core of this project as it provides 12 years of L-band observations and most importantly, L-band observations are very sensitive to the biomass and its vegetation water content, as shown by recent studies."

<caption>Figure 11.1: Main page of the OSMOSE webpage.</caption>

## 11.2 Outreach to other projects

The outputs of the OSMOSE project are shared with the ESA SMOS ESL group, and in particular within the context of the SMOS 10 years SM project, led by CESBIO, as the change request n°4 ESL (Contract No. 4000130567).

## 11.3 Publications

The following publications were done in the context of the OSMOSE project:

* Zhou, Y. and Schwank, M. and Kurum, M. and Houtz, D. and Zhao, Q. and Lang, R. and Mialon, A. and Drusch, D., 2024. Microwave Emission Model for Layered Vegetation (MEMLV): An Exemplary Study for Coniferous Forests from P- to Ka- Band. In review TGRS, 2024. doi: https://doi.org/XXX.

* Salazar-Neira, J. C. and A. Mialon, and P. Richaume, and S. Mermoz, and Y.H. Kerr, and A. Bouvet, and T. Le Toan, and S. Boitard, and N. J. Rodríguez-Fernández. Above-ground biomass estimation based on multi-angular I-band passive microwaves brightness temperatures. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 16:5813-5827, 2023. doi: 10.1109/JSTARS.2023.3285288
* Salazar Neira, J. C. and Mialon, A. and Le Toan, T. and Bouvet, A. and Richaume, P. and Rodríguez-Fernández, N. J., 2024. The response of Microwave Optical Depth and Optical Indices to Water Availability in Croplands, in review Remote Sensing of Environment

* Ma, H. and Weiss, M. and Malik, D. and Berthelot, B. and Yebra, M. and Nolan, R. and Mialon, A. and Zeng, J. and Quan, X. and Tagesson, H. T. and Olioso, A. and Baret, 2024. Satellite vegetation canopy water from Sentinel-2, Landsat8 and MODIS: the principle, algorithm and assessment, submitted to RSE.

* Ma, H., Weiss M., et al., & Baret, F. Blending and comparison of satellite LAI and FAPAR products from multi-sensors using diverse fusion approaches, (to be submitted to RSE)

* Boitard, S. and Mialon, A. and Mermoz, S. and Rodríguez-Fernández, N.J., and Richaume, P. and Salazar Neira, J.C. and Tarot, S. and Yann H. Kerr, 2024. Above ground biomass dataset from SMOS L band vegetation optical depth and reference maps. in review ESSD.

## 11.4 Conference and workshop contributions

The OSMOSE team participated to various conferences to present the activities, as shown in Table 11.1. The OSMOSE project was regularly presented at the ESA SMOS ESL (Expert Support Laboratory), i.e. PM n°27 in Toulouse-France (05-October-2022), PM n°28 in Roma-Italy (18-January-2023), PM n°29 in Zurich-Switzerland (16-May-2023), PM n°30 in Toulouse-France (03-October-2023), PM n°31 in Sodankyla-Finland (24-January-2024), PM n°32 in Zurich-Switzerland (15-May-2024), PM n°33 in Toulouse-France (09-October-2024).

<table><caption>Table 11.1: Conferences and Workshops</caption><thead><tr><td>Conference/Workshop</td><td>Location</td><td>Date</td><td>Website</td></tr></thead><tbody><tr><td>IGARSS</td><td>Athens, Greece</td><td>7-12 July 2024</td><td>https://www.2024.ieeeigarss.org/index.php</td></tr><tr><td>EGU</td><td>Vienna, Austria</td><td>14-19 April 2024</td><td>https://www.egu24.eu/</td></tr><tr><td>IGARSS</td><td>Pasadena, USA</td><td>16-21 July 2023</td><td>https://2023.ieeeigarss.org/</td></tr><tr><td>ESA PollnSar</td><td>Toulouse, France</td><td>19-23 June 2023</td><td>https://polinsar-biomass2023.esa.int/</td></tr><tr><td>ESA LPVE</td><td>Frascati ESRIN, Italy</td><td>ESA LPVE</td><td>https://nikal.eventsair.com/lpve-2023/</td></tr><tr><td>EGU 2022</td><td>Vienne, Austria</td><td>01 May 2023</td><td>https://www.egu22.eu/</td></tr><tr><td>ESA 4th Carbon from Space</td><td>Frascati, ESRIN, Italy</td><td>25 Nov 2022</td><td>https://4thcarbonfromspace.esa.int/</td></tr><tr><td>IGARSS</td><td>Kuala Lumpur, Malaysia</td><td>17-22 July 2022</td><td>https://igarss2022.org/</td></tr><tr><td>ESA Living Planet Symposium</td><td>Bonn, Germany</td><td>23-27 May 2022</td><td>https://lps22.esa.int</td></tr></tbody></table>
# Chapter 12

# Conclusions

The ESA PM-VOS project aimed at studying the vegetation layer by remote sensing data from passive microwaves and optical/Infra red observations. It was the continuity of the ESA "SMOS+veg" project (Contract No. 4000137990/22/NL/IA) with the main goals to: i) improve AGB estimates from passive microwave sensors, such as SMOS, AMSR-E, and AMSR-2, and ii) develop harmonized datasets of microwave band-specific VOD to assess evolutions of vegetation parameters (AGB, for instance) over the last 20 years from passive microwave sensors; iii) provide a 20-year data-cube of optical indices (NDVI, LAI,FAPAR); iv) study the complementarity of all these indices to characterize and monitor state parameters of vegetation.
This project lasted two years from June 2022 to October 2024, and provided the following important outputs:

* First, a microwave emission model for layered vegetation (MEMLV) was developed to estimate the effective VOD and scattering albedo of forest canopies from P- to Ka-band (0.4 GHz – 37 GHz). The MEMLV provided a physics-based approach to assess the impact from different parameters, such as frequency, forest structure and vegetation water content, on the retrieved VOD and albedo.

* Second, leaf water content was estimated from optical data

* Third, estimations of AGB were performed from passive microwave observations, including their first assessment based on in situ measurements. AGB were derived from VOD and from a neural network approach.

* Finally, the synergy of all these indices was investigated for croplands. Passive microwave VOD were also used in synergy with vegetation indices from optical/near-infrared data, showing the complementarity of these indices to study the seasonal dynamics of the vegetation.
References

Aires, F., P. Weston, P. de Rosnay, and D. Fairbairn (2021), Statistical approaches to assimilate ascat soil moisture information—i. methodologies and first assessment, *Quarterly Journal of the Royal Meteorological Society*, **147**(736), p. 1823–1852.

Al Bitar, A., A. Mialon, Y. Kerr, F. Cabot, P. Richaume, E. Jacquette, A. Quesney, A. Mahmoodi, S. Tarot, M. Parrens, A. Al-yaari, T. Pellarin, N. Rodriguez-Fernandez, and J.-P. Wigneron (2017a), The globl smos level 3 daily soil moisture and brightness temperature maps, *Earth System Science Data*, **9**, p. 293–315, doi:10.5194/essd-2017-1.

Al Bitar, A., A. Mialon, Y. H. Kerr, F. Cabot, P. Richaume, E. Jacquette, A. Quesney, A. Mahmoodi, S. Tarot, M. Parrens, A. Al-Yaari, T. Pellarin, N. Rodriguez-Fernandez, and J.-P. Wigneron (2017b), The global SMOS Level 3 daily soil moisture and brightness temperature maps, *Earth System Science Data*, **9**(1), p. 293–315, doi:10.5194/essd-9-293-2017.

Andela, N., Y. Liu, A. Van Dijk, R. De Jeu, and T. McVicar (2013), Global changes in dryland vegetation dynamics (1988–2008) assessed by satellite remote sensing: comparing a new passive microwave vegetation density record with reflective greenness data, *Biogeosciences*, **10**(10), p. 6657–6676.

Atzberger, C. (2004), Object-based retrieval of biophysical canopy variables using artificial neural nets and radiative transfer models, *Remote Sensing of Environment*, **93**, p. 53–67.

Baccini, A., S. J. Goetz, W. S. Walker, N. T. Laporte, M. Sun, D. Sulla-Menashe, J. Hackler, P. S. A. Beck, R. Dubayah, M. A. Friedl, S. Samanta, and R. A. Houghton (2012), Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps, *Nature Climate Change*, **2**(3), p. 182–185, doi:10.1038/nclimate1354.

Bacour, C., F. Baret, D. Béal, M. Weiss, and K. Pavageau (2006), Neural network estimation of lai, fapar, fcover and lai×cab, from top of canopy meris reflectance data: Principles and validation, *Remote Sensing of Environment*, **105**, p. 313–325, doi:10.1016/j.rse.2006.07.014.

Baret, F., and M. Weiss (2024), Review of methods to estimate vegetation traits from reflectance observations, *Report*, INRAE.

Baret, F., O. Hagolle, B. Geiger, P. Bicheron, B. Miras, M. Huc, B. Berthelot, F. Nino, M. Weiss, O. Samain, J. L. Roujean, and M. Leroy (2007), Lai, fapar and fcover cyclopes global products derived from vegetation: Part 1: Principles of the algorithm, *Remote Sensing of Environment*, **110**(3), p. 275–286.

Baret, F., M. Weiss, R. Lacaze, F. Camacho, H. Makhmara, P. Pacholczyk, and B. Smets (2013), Geov1: Lai and fapar essential climate variables and fcover global time series capitalizing over existing products. part1: Principles of development and production, *Remote Sensing of Environment*, **137**(0), p. 299–309, doi:http://dx.doi.org/10.1016/j.rse.2012.12.027.

Barrée, M., A. Mialon, T. Pellarin, M. Parrens, R. Biron, F. Lemaître, S. Gascoin, and Y. Kerr (2021), Soil moisture and vegetation optical depth retrievals over heterogeneous scenes using LEWIS L-band radiometer, *International Journal of Applied Earth Observation and Geoinformation*, **102**, p. 102,424, doi:10.1016/j.jag.2021.102424.
Beck, H. E., T. R. McVicar, N. Vergopolan, A. Berg, N. J. Lutsko, A. Dufour, Z. Zeng, X. Jiang, A. I. J. M. van Dijk, and D. G. Miralles (2023), High-resolution (1 km) Köppen-Geiger maps for 1901–2099 based on constrained CMIP6 projections, *Scientific Data*, **10**(1), p. 724, doi:10.1038/s41597-023-02549-6.

Bousquet, E., A. Mialon, N. Rodriguez-Fernandez, C. Prigent, F. H. Wagner, and Y. H. Kerr (2021), Influence of surface water variations on VOD and biomass estimates from passive microwave sensors, *Remote Sensing of Environment*, **257**, p. 112,345, doi:10.1016/j.rse.2021.112345.

Brandt, M., J.-P. Wigneron, J. Chave, T. Tagesson, J. Penuelas, P. Ciais, K. Rasmussen, F. Tian, C. Mbow, A. Al-Yaari, N. Rodriguez, G. Schurgers, W. Zhang, J. Chang, Y. Kerr, A. Verger, C. Tucker, A. Mialon, and R. Fensholt (2018), Satellite passive microwaves reveal deforestation and drought-induced carbon losses in Africa, *Nature Ecology Evolution*, **Vol. 2 Issue 5**, p. 827–835, doi:10.1038/s41559-018-0530-6.

Brocca, L., S. Hasenauer, T. Lacava, F. Melone, T. Moramarco, W. Wagner, W. Dorigo, P. Matgen, J. Martínez-Fernández, P. Llorens, J. Latron, C. Martin, and M. Bittelli (2011), Soil moisture estimation through ascat and amsr-e sensors: An intercomparison and validation study across europe, *Remote Sensing of Environment*, **115**(12), p. 3390–3408, doi:https://doi.org/10.1016/j.rse.2011.08.003.

Brodzik, M. J., and D. G. Long (2018), NASA measures calibrated passive microwave daily ease-grid 2.0 brightness temperature esdr (cetb) algorithm theoretical basis document, revision 1.0 26 April 2018.

Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. Savoie (2012a), EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets, *ISPRS International Journal of Geo-Information*, **1**(1), p. 32–45, doi:10.3390/ijgi1010032, ISSN 2220-9964.

Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. Savoie (2012b), EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets, *ISPRS International Journal of Geo-Information*, **1**(1), p. 32–45, doi:10.3390/ijgi1010032.

Brodzik, M. J., B. Billingsley, T. Haran, B. Raup, and M. H. Savoie (2014), Correction: Brodzik, M.J., et al. EASE-Grid 2.0: Incremental but Significant Improvements for Earth-Gridded Data Sets. *ISPRS International Journal of Geo-Information* 2012, 1, 32–45, *ISPRS International Journal of Geo-Information*, **3**(3), p. 1154–1156, doi:10.3390/ijgi3031154, ISSN 2220-9964.

Brodzik, M. J., D. G. Long, M. A. Hardman, A. Paget, and R. Armstrong (2016), Measures calibrated enhanced-resolution passive microwave daily ease-grid 2.0 brightness temperature esdr, version 1.

Brown, L. A., R. Fernandes, N. Djamai, C. Meier, N. Gobron, H. Morris, F. Canisius, G. Bai, C. Lerebourg, and C. Lanconelli (2021), Validation of baseline and modified sentinel-2 level 2 prototype processor leaf area index retrievals over the united states, *ISPRS Journal of Photogrammetry and Remote Sensing*, **175**, p. 71–87.

Broxton, P. D., X. Zeng, D. Sulla-Menashe, and P. A. Troch (2014a), A global land cover climatology using modis data, *Journal of Applied Meteorology and Climatology*, **53**(6), p. 1593 – 1605, doi:https://doi.org/10.1175/JAMC-D-13-0270.1.

Broxton, P. D., X. Zeng, D. Sulla-Menashe, and P. A. Troch (2014b), A Global Land Cover Climatology Using MODIS Data, *Journal of Applied Meteorology and Climatology*, **53**(6), p. 1593–1605, doi:10.1175/JAMC-D-13-0270.1.

Ceccato, P., S. Flasse, S. Tarantola, S. Jacquemoud, and J.-M. Grégoire (2001), Detecting vegetation leaf water content using reflectance in the optical domain, *Remote Sensing of Environment*, **77**(1), p. 22–33, doi:https://doi.org/10.1016/S0034-4257(01)00191-2.

Chauhan, N. S., D. M. Le Vine, and R. H. Lang (1994), Discrete scatter model for microwave radar and radiometer response to corn: comparison of theory and data, *IEEE Transactions on Geoscience and Remote Sensing*, **32**(2), p. 416–426.

Combal, B., F. Baret, and M. Weiss (2002), Improving canopy variable estimation from remote sensing data by exploiting ancillary information. case study on sugar beet canopies., *Agronomie*, **22**, p. 205–215, in press.
Djamai, N., R. Fernandes, M. Weiss, H. McNairn, and K. Goïta (2019), Validation of the sentinel simplified level 2 product prototype processor (sl2p) for mapping cropland biophysical variables using sentinel-2/msi and landsat-8/oli data, *Remote Sensing of Environment*, **225**, p. 416-430, doi:https://doi.org/10.1016/j.rse.2019.03.020.

Dong, M., S. Liu, R. Jiang, J. Qi, B. de Solan, A. Comar, L. Li, W. Li, Y. Ding, and F. Baret (2024), Comparing and combining data-driven and model-driven approaches to monitor wheat green area index with high spatio-temporal resolution satellites, *Remote Sensing of Environment*, **305**, p. 114-118, doi: https://doi.org/10.1016/j.rse.2024.114118.

Estévez, J., J. Vicent, J. P. Rivera-Caicedo, P. Morcillo-Pallarés, F. Vuolo, N. Sabater, G. Camps-Valls, J. Moreno, and J. Verrelst (2020), Gaussian processes retrieval of lai from sentinel-2 top-of-atmosphere radiance data, *ISPRS Journal of Photogrammetry and Remote Sensing*, **167**, p. 289-304, doi:10.1016/j.isprsjprs.2020.07.004, export Date: 25 January 2021.

Fan, L., J.-P. Wigneron, P. Ciais, J. Chave, M. Brandt, R. Fensholt, S. S. Saatchi, A. Bastos, K. Al-Yaari, A. ad Hufkens, Y. Qin, X. Xiao, C. Chen, R. B. Myneni, R. Fernandez-Moran, A. Mialon, N. Rodriguez-Fernandez, Y. Kerr, F. Tian, and J. Peñuelas (2019), Satellite-observed pantropical carbon dynamics, *Nature plants*, **5**(9), p. 944-951.

FAO, F. (1988), Soil map of the world: revised legend, *World Soil Resources Report*, **60**.

Fernandez-Moran, R., A. Al-Yaari, A. Mialon, A. Mahmoodi, A. Al Bitar, G. De Lannoy, N. Rodriguez-Fernandez, E. Lopez-Baeza, Y. Kerr, and J.-P. Wigneron (2017), SMOS-IC: An Alternative SMOS Soil Moisture and Vegetation Optical Depth Product, *Remote Sensing*, **9**(5), doi:10.3390/rs9050457, article Number 457.

Ferrazzoli, P., and L. Guerriero (1996), Passive microwave remote sensing of forests: A model investigation, *IEEE Transactions on Geoscience and Remote Sensing*, **34**(2), p. 433-443.

Fortuny, J., and A. J. Sieber (1999), Three-dimensional synthetic aperture radar imaging of a fir tree: First results, *IEEE Transactions on Geoscience and Remote Sensing*, **37**(2), p. 1006-1014.

Fung, A. K., and K.-S. Chen (2010), *Microwave scattering and emission models for users*, Artech house.

Fuster, B., J. Sánchez-Zapero, F. Camacho, V. García-Santos, A. Verger, R. Lacaze, M. Weiss, F. Baret, and B. Smets (2020), Quality assessment of proba-v lai, fapar and fcover collection 300 m products of copernicus global land service, *Remote Sensing*, **12**(6), doi:10.3390/rs12061017.

Garcia-Haro, F. J., M. Campos-Taberner, A. Moreno, H. T. Tagesson, F. Camacho, B. Martínez, S. Sánchez, M. Piles, G. Camps-Valls, and M. Yebra (2020), A global canopy water content product from avhrr/metop, *ISPRS Journal of Photogrammetry and Remote Sensing*, **162**, p. 77-93.

Garonna, I., R. De Jong, A. J. De Wit, C. A. Mücher, B. Schmid, and M. E. Schaepman (2014), Strong contribution of autumn phenology to changes in satellite-derived growing season length estimates across Europe (1982-2011), *Global Change Biology*, **20**(11), p. 3457-3470, doi:10.1111/gcb.12625.

Goel, N. (1989), *Inversion of canopy reflectance models for estimation of biophysical parameters from reflectance data*, p. 205-251, asrar, g. ed., Wiley Interscience.

Grant, J. P., J.-P. Wigneron, R. De Jeu, H. Lawrence, A. Mialon, P. Richaume, A. Al Bitar, M. Drusch, M. Van Marle, and Y. Kerr (2016), Comparison of smos and amsr-e vegetation optical depth to four modis-based vegetation indices, *Remote Sensing of Environment*, **172**, p. 87-100, doi:http://dx.doi.org/10.1016/j.rse.2015.10.021.

Griebel, A., M. M. Boer, C. Blackman, B. Choat, D. S. Ellsworth, P. Madden, B. Medlyn, V. Resco de Dios, A. Wujeska-Klause, and M. Yebra (2023), Specific leaf area and vapour pressure deficit control live fuel moisture content, *Functional Ecology*, **37**(3), p. 719-731.

Guan, K., J. Wu, J. S. Kimball, M. C. Anderson, S. Frolking, B. Li, C. R. Hain, and D. B. Lobell (2017), The shared and unique values of optical, fluorescence, thermal and microwave satellite data for estimating large-scale crop yields, *Remote Sensing of Environment*, **199**, p. 333-349, doi:10.1016/j.rse.2017.06.043.
Guglielmetti, M., M. Schwank, C. Mätzler, C. Oberdörster, J. Vanderborght, and H. Flühler (2007), Measured microwave radiative transfer properties of a deciduous forest canopy, *Remote Sensing of Environment*, **109**(4), p. 523–532, doi:https://doi.org/10.1016/j.rse.2007.02.003.

Hansen, M. C., P. V. Potapov, R. Moore, M. Hancher, S. A. Turubanova, A. Tyukavina, D. Thau, S. V. Stehman, S. J. Goetz, T. R. Loveland, et al. (2013), High-resolution global maps of 21st-century forest cover change, *science*, **342**(6160), p. 850–853.

Hengl, T., J. Sanderman, and L. Parente (2021), Open soil spectral library (training data and calibration models), doi:https://zenodo.org/record/5805138.

Hersbach, H., B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut (2023), Era5 hourly data on single levels from 1940 to present. copernicus climate change service (c3s) climate data store (cds), *Tech. rep.*, ECMWF, doi:10.24381/cds.adbb2d47, accessed on 2022.

Hornbuckle, B. (2023), Vegetation water content... or vegetation canopy water?, in *IGARSS 2023 - 2023 IEEE International Geoscience and Remote Sensing Symposium*, p. 3163–3166, doi:10.1109/IGARSS52108. 2023.10283364.

Hornbuckle, B. K., A. W. England, R. D. De Roo, M. A. Fischman, and D. L. Boprie (2003), Vegetation canopy anisotropy at 1.4 ghz, *IEEE Transactions on Geoscience and Remote Sensing*, **41**(10), p. 2211–2223.

Hornbuckle, B. K., J. C. Patton, A. VanLoocke, A. E. Suyker, M. C. Roby, V. A. Walker, E. R. Iyer, D. E. Herzmann, and E. A. Endacott (2016), SMOS optical thickness changes in response to the growth and development of crops, crop management, and weather, *Remote Sensing of Environment*, **180**, p. 320–333, doi:10.1016/j.rse.2016.02.043.

Huete, A., K. Didan, T. Miura, E. Rodriguez, X. Gao, and L. Ferreira (2002), Overview of the radiometric and biophysical performance of the MODIS vegetation indices, *Remote Sensing of Environment*, **83**(1-2), p. 195–213, doi:10.1016/S0034-4257(02)00096-2.

Huxman, T. E., K. A. Snyder, D. Tissue, A. J. Leffler, K. Ogle, W. T. Pockman, D. R. Sandquist, D. L. Potts, and S. Schwinning (2004), Precipitation pulses and carbon fluxes in semiarid and arid ecosystems, *Oecologia*, **141**(2), p. 254–268, doi:10.1007/s00442-004-1682-4.

Jackson, T., and T. Schmugge (1991), Vegetation effects on the microwave emission of soils, *Remote Sensing of Environment*, **36**(3), p. 203–212.

Jackson, T., T. Schmugge, and J. Wang (1982), Passive microwave sensing of soil moisture under vegetation canopies, *Water Resour. Res.*, **18**, p. 1137–1142.

Jacquemoud, S., W. Verhoef, F. Baret, C. Bacour, P. J. Zarco-Tejada, G. P. Asner, C. François, and S. L. Ustin (2009), Prospect+sail models: A review of use for vegetation characterization, *Remote Sensing of Environment*, **113**, p. S56–S66, doi:https://doi.org/10.1016/j.rse.2008.01.026.

Jarque, C. M., and A. K. Bera (1987), A test for normality of observations and regression residuals, *International Statistical Review/Revue Internationale de Statistique*, p. 163–172.

Jones, M. O., L. A. Jones, J. S. Kimball, and K. C. McDonald (2011), Satellite passive microwave remote sensing for monitoring global land surface phenology, *Remote Sensing of Environment*, **115**(4), p. 1102–1114, doi:10.1016/j.rse.2010.12.015.

Jones, M. O., J. S. Kimball, L. A. Jones, and K. C. McDonald (2012), Satellite passive microwave detection of North America start of season, *Remote Sensing of Environment*, **123**, p. 324–333, doi:10.1016/j.rse.2012.03.025.

Karkauskaite, P., T. Tagesson, and R. Fensholt (2017), Evaluation of the Plant Phenology Index (PPI), NDVI and EVI for Start-of-Season Trend Analysis of the Northern Hemisphere Boreal Zone, *Remote Sensing*, **9**(5), p. 485, doi:10.3390/rs9050485.
KATTGE, J., S. DÍAZ, S. LAVOREL, I. C. PRENTICE, P. LEADLEY, G. BÖNISCH, E. GARNIER, M. WESTOBY, P. B. REICH, I. J. WRIGHT, J. H. C. CORNELISSEN, C. VIOLLE, S. P. HARRISON, P. M. Van BODEGOM, M. REICHSTEIN, B. J. ENQUIST, N. A. SOUDZILOVSKAIA, D. D. ACKERLY, M. ANAND, O. ATKIN, M. BAHN, T. R. BAKER, D. BALDOCCHI, R. BEKKER, C. C. BLANCO, B. BLONDER, W. J. BOND, R. BRADSTOCK, D. E. BUNKER, F. CASANOVES, J. CAVENDER-BARES, J. Q. CHAMBERS, F. S. CHAPIN III, J. CHAVE, D. COOMES, W. K. CORNWELL, J. M. CRAINE, B. H. DOBRIN, L. DUARTE, W. DURKA, J. ELSER, G. ESSER, M. ESTIARTE, W. F. FAGAN, J. FANG, F. FERNÁNDEZ-MÉNDEZ, A. FIDELIS, B. FINEGAN, O. FLORES, H. FORD, D. FRANK, G. T. FRESCHET, N. M. FYLLAS, R. V. GALLAGHER, W. A. GREEN, A. G. GUTIERREZ, T. HICKLER, S. I. HIGGINS, J. G. HODGSON, A. JALILI, S. JANSEN, C. A. JOLY, A. J. KERKHOFF, D. KIRKUP, K. KITAJIMA, M. KLEYER, S. KLOTZ, J. M. H. KNOPS, K. KRAMER, I. KÜHN, H. KUROKAWA, D. LAUGHLIN, T. D. LEE, M. LEISHMAN, F. LENS, T. LENZ, S. L. LEWIS, J. LLOYD, J. LLUSIÀ, F. LOUAULT, S. MA, M. D. MAHECHA, P. MANNING, T. MASSAD, B. E. MEDLYN, J. MESSIER, A. T. MOLES, S. C. MÜLLER, K. NADROWSKI, S. NAEEM, U. NIINEMETS, S. NÖLLERT, A. NÜSKE, R. OGAYA, J. OLEKSYN, V. G. ONIPCHENKO, Y. ONODA, J. ORDOÑEZ, G. OVERBECK, W. A. OZINGA, S. PATIÑO, S. PAULA, J. G. PAUSAS, J. PEÑUELAS, O. L. PHILLIPS, V. PILLAR, H. POORTER, L. POORTER, P. POSCHLOD, A. PRINZING, R. PROULX, A. RAMMIG, S. REINSCH, B. REU, L. SACK, B. SALGADO-NEGRET, J. SARDANS, S. SHIODERA, B. SHIPLEY, A. SIEFERT, J.-F. SOSINSKI, E. SOUSSANA, E. SWAINE, N. SWENSON, K. THOMPSON, P. THORNTON, M. WALDRAM, E. WEIHER, M. WHITE, S. WHITE, S. J. WRIGHT, B. YGUEL, S. ZAEHLE, A. E. ZANNE, and C. WIRTH (2011), Try – a global database of plant traits, *Global Change Biology*, **17**(9), p. 2905–2935, doi:https://doi.org/10.1111/j.1365-2486.2011.02451.x.

Kattge, J., G. Bönisch, S. Díaz, S. Lavorel, I. C. Prentice, and et al (2020), Try plant trait database – enhanced coverage and open access, *Global Change Biology*, **26**(1), p. 119–188, doi:https://doi.org/10.1111/gcb. 14904.

Kerr, Y., P. Richaume, P. Waldteufel, P. Ferrazzoli, J. Wigneron, M. Schwank, and K. Rautiainen (2020), Algorithm theoretical basis document (ATBD) for the smos level 2 soil moisture processor, Technical Report TN-ESL-SM-GS-0001-4b SM-ESL (CBSA), p. 145p.

Kerr, Y. H., P. Waldteufel, J.-P. Wigneron, S. Delwart, F. Cabot, J. Boutin, M.-J. Escorihuela, J. Font, N. Reul, C. Gruhier, et al. (2010), The SMOS mission: New tool for monitoring key elements of the global water cycle, *Proceedings of the IEEE*, **98**(5), p. 666–687.

Kerr, Y. H., P. Waldteufel, P. Richaume, J. P. Wigneron, P. Ferrazzoli, A. Mahmoodi, A. Al Bitar, F. Cabot, C. Gruhier, S. E. Juglea, D. Leroux, A. Mialon, and S. Delwart (2012a), The SMOS Soil Moisture Retrieval Algorithm, *Geoscience and Remote Sensing, IEEE Transactions on*, **50**(5), p. 1384–1403, doi:10.1109/ TGRS.2012.2184548.

Kerr, Y. H., P. Waldteufel, P. Richaume, J. P. Wigneron, P. Ferrazzoli, A. Mahmoodi, A. Al Bitar, F. Cabot, C. Gruhier, S. E. Juglea, D. Leroux, A. Mialon, and S. Delwart (2012b), The SMOS Soil Moisture Retrieval Algorithm, *IEEE Transactions on Geoscience and Remote Sensing*, **50**(5), p. 1384–1403, doi:10.1109/ TGRS.2012.2184548.

Kirdiashev, K. P., A. A. Chukhlantsev, and A. M. Shutko (1979), Microwave radiation of the earth's surface in the presence of vegetation cover, *Radiotekhnika i Elektronika*, **24**, p. 256–264.

Knyazikhin, Y., J. V. Martonchik, R. B. Myneni, D. J. Diner, and S. W. Running (1998), Synergistic algorithm for estimating vegetation canopy leaf area index and fraction of absorbed photosynthetically active radiation from modis and misr data, *J. Geophys. Res.*, **103**(D24), p. 32,257–32,275, doi:10.1029/98jd02462.

Kramer, P. (2012), *Physiology of woody plants*, Elsevier.

Kurum, M. (2013), Quantifying scattering albedo in microwave emission of vegetated terrain, *Remote Sensing of Environment*, **129**, p. 66–74.

Kurum, M., P. E. O'Neill, R. H. Lang, A. T. Joseph, M. H. Cosh, and T. J. Jackson (2012), Effective tree scattering and opacity at L-band, *Remote Sensing of Environment*, **118**, p. 1–9.
Landry, R., R. Fournier, F. Ahern, and R. Lang (1997), Tree vectorization: a methodology to characterize fine tree architecture in support of remote sensing models, *Canadian journal of remote sensing*, **23**(2), p. 91–107.

Lang, R., R. Landry, A. Franchois, G. Nesti, and A. Sieber (1998), Microwave tree scattering experiment: comparison of theory and experiment, in *1998 IEEE International Geoscience and Remote Sensing. Symposium Proceedings*, vol. 5, p. 2384–2386, IEEE.

Lang, R. H., and J. S. Sidhu (1983), Electromagnetic backscattering from a layer of vegetation: A discrete approach, *IEEE Transactions on Geoscience and Remote Sensing*, **GE-21**(1), p. 62–71.

Laurent, V. C., W. Verhoef, A. Damm, M. E. Schaepman, and J. G. Clevers (2013), A bayesian object-based approach for estimating vegetation biophysical and biochemical variables from apex at-sensor radiance data, *Remote Sensing of Environment*, **139**, p. 6–17.

Lauvernet, C., F. Baret, L. Hascoët, S. Buis, and F.-X. Le Dimet (2008), Multitemporal-patch ensemble inversion of coupled surface-atmosphere radiative transfer models for land surface characterization, *Remote Sensing of Environment*, **112**(3), p. 851–861, doi:http://dx.doi.org/10.1016/j.rse.2007.06.027.

Lawrence, H., J.-P. Wigneron, P. Richaume, N. Novello, J. Grant, A. Mialon, A. Al Bitar, O. Merlin, D. Guyon, D. Leroux, S. Bircher, and Y. Kerr (2014a), Comparison between smos vegetation optical depth products and modis vegetation indices over crop zones of the usa, *Remote Sensing of Environment*, **140**, p. 396–406, doi:https://doi.org/10.1016/j.rse.2013.07.021.

Lawrence, H., J.-P. Wigneron, P. Richaume, N. Novello, J. Grant, A. Mialon, A. Al Bitar, O. Merlin, D. Guyon, D. Leroux, S. Bircher, and Y. Kerr (2014b), Comparison between SMOS Vegetation Optical Depth products and MODIS vegetation indices over crop zones of the USA, *Remote Sensing of Environment*, **140**, p. 396–406, doi:10.1016/j.rse.2013.07.021.

Li, C., L. Zhou, and W. Xu (2021), Estimating Aboveground Biomass Using Sentinel-2 MSI Data and Ensemble Algorithms for Grassland in the Shengjin Lake Wetland, China, *Remote Sensing*, **13**(8), p. 1595, doi: 10.3390/rs13081595.

Li, X., A. Al-Yaari, M. Schwank, L. Fan, F. Frappart, J. Swenson, and J.-P. Wigneron (2020), Compared performances of SMOS-IC soil moisture and vegetation optical depth retrievals based on Tau-Omega and Two-Stream microwave emission models, *Remote Sensing of Environment*, **236**, p. 111,502.

Li, X., J.-P. Wigneron, L. Fan, F. Frappart, S. H. Yueh, A. Colliander, A. Ebtehaj, L. Gao, R. Fernandez-Moran, X. Liu, M. Wang, H. Ma, C. Moisy, and P. Ciais (2022), A new SMAP soil moisture and vegetation optical depth product (SMAP-IB): Algorithm, assessment and inter-comparison, *Remote Sensing of Environment*, **271**, p. 112,921, doi:https://doi.org/10.1016/j.rse.2022.112921.

Liang, S. (2005), *Quantitative remote sensing of land surfaces*, John Wiley & Sons.

Lilliefors, H. W. (1967), On the kolmogorov-smirnov test for normality with mean and variance unknown, *Journal of the American statistical Association*, **62**(318), p. 399–402.

Liu, X., Y. Chen, Z. Li, and Y. Li (2023), Evaluating the Consistency of Vegetation Phenological Parameters in the Northern Hemisphere from 1982 to 2015, *Remote Sensing*, **15**(10), p. 2559, doi:10.3390/rs15102559.

Liu, Y., R. A. de Jeu, A. I. van Dijk, and M. Owe (2007), Trmm-tmi satellite observed soil moisture and vegetation density (1998–2005) show strong connection with el niño in eastern australia, *Geophysical Research Letters*, **34**(15).

Liu, Y., N. M. Holtzman, and A. G. Konings (2021), Global ecosystem-scale plant hydraulic traits retrieved using model-data fusion, *Hydrology and Earth System Sciences*, **25**(5), p. 2399–2417, doi: 10.5194/hess-25-2399-2021.

Liu, Y. Y., R. A. De Jeu, M. F. McCabe, J. P. Evans, and A. I. Van Dijk (2011), Global long-term passive microwave satellite-based retrievals of vegetation optical depth, *Geophysical Research Letters*, **38**(18).

Liu, Y. Y., A. I. Van Dijk, R. A. De Jeu, J. G. Canadell, M. F. McCabe, J. P. Evans, and G. Wang (2015), Recent reversal in loss of global terrestrial biomass, *Nature Climate Change*, **5**(5), p. 470–474.
Loveland, T., and A. Belward (1997), The International Geosphere Biosphere Programme Data and Information System global land cover data set (DISCover), *Acta Astronautica*, **41**(4-10), p. 681-689, doi:10.1016/S0094-5765(98)00050-2.

Loveland, T. R., B. C. Reed, J. F. Brown, D. O. Ohlen, Z. Zhu, L. Yang, and J. W. Merchant (2000), Development of a global land cover characteristics database and IGBP DISCover from 1 km AVHRR data, *International Journal of Remote Sensing*, **21**(6-7), p. 1303-1330, doi:10.1080/014311600210191.

Madelon, R., N. J. Rodríguez-Fernández, R. van der Schalie, T. Scanlon, A. Al Bitar, Y. H. Kerr, R. de Jeu, and W. Dorigo (2022), Toward the removal of model dependency in soil moisture climate data records by using an I-band scaling reference, *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, **15**, p. 831-848, doi:10.1109/JSTARS.2021.3137008.

Maeda, T., Y. Taniguchi, and K. Imaoka (2016), Gcom-w1 amsr2 level 1r product: Dataset of brightness temperature modified using the antenna pattern matching technique, *IEEE Transactions on Geoscience and Remote Sensing*, **54**(2), p. 770-782, doi:10.1109/TGRS.2015.2465170.

Matzler, C. (1994), Microwave (1-100 GHz) dielectric model of leaves, *IEEE Transactions on Geoscience and Remote Sensing*, **32**(4), p. 947-949.

Mätzler, C. (1998), Improved born approximation for scattering of radiation in a granular medium, *Journal of Applied Physics*, **83**(11), p. 6111-6117.

McMahon, T. A., and R. E. Kronauer (1976), Tree structures: deducing the principle of mechanical design, *Journal of theoretical biology*, **59**(2), p. 443-466.

Meissner, T., F. J. Wentz, and L. Ricciardulli (2014), The emission and scattering of L-band microwave radiation from rough ocean surfaces and wind speed measurements from the uppercaseAquarius sensor, *Journal of Geophysical Research: Oceans*, **119**(9), p. 6499-6522.

Mialon, A., N. J. Rodríguez-Fernández, M. Santoro, S. Saatchi, S. Mermoz, E. Bousquet, and Y. H. Kerr (2020), Evaluation of the sensitivity of smos I-vod to forest above-ground biomass at global scale, *Remote Sensing*, **12**(9), p. 1450.

Mironov, V. L., L. G. Kosolapova, and S. V. Fomin (2009), Physically and mineralogically based spectroscopic dielectric model for moist soils, *IEEE Transactions on Geoscience and Remote Sensing*, **47**(7), p. 2059-2070.

Mo, T., B. Choudhury, T. Schmugge, J. R. Wang, and T. Jackson (1982), A model for microwave emission from vegetation-covered fields, *Journal of Geophysical Research: Oceans*, **87**(C13), p. 11,229-11,237.

Moesinger, L., W. Dorigo, R. De Jeu, R. Van Der Schalie, T. Scanlon, I. Teubner, and M. Forkel (2020a), The global long-term microwave vegetation optical depth climate archive (vodca), *Earth System Science Data*, **12**(1), p. 177-196.

Moesinger, L., W. Dorigo, R. de Jeu, R. van der Schalie, T. Scanlon, I. Teubner, and M. Forkel (2020b), The global long-term microwave Vegetation Optical Depth Climate Archive (VODCA), *Earth System Science Data*, **12**(1), p. 177-196, doi:10.5194/essd-12-177-2020.

Montpetit, B., A. Royer, J.-P. Wigneron, A. Chanzy, and A. Mialon (2015), Evaluation of multi-frequency bare soil microwave reflectivity models, *Remote Sensing of Environment*, **162**, p. 186 – 195, doi:https://doi.org/10.1016/j.rse.2015.02.015.

Myneni, R., V. Gutschick, G. Asrar, and E. Kanemasu (1988), Photon transport in vegetation canopies with anisotropic scattering part ii. discrete-ordinates/exact-kernel technique for one-angle photon transport in slab geometry., *Agricultural and Forest Meteorology*, **42**, p. 17-40.

Nguyen, T. H., M. Langensiepen, T. Gaiser, H. Webber, H. Ahrends, H. Hueging, and F. Ewert (2022), Responses of winter wheat and maize to varying soil moisture: From leaf to canopy, *Agricultural and Forest Meteorology*, **314**, p. 108,803, doi:10.1016/j.agrformet.2021.108803.
Owe, M., R. De Jeu, and T. Holmes (2008), Multisensor historical climatology of satellite-derived global land surface moisture, *Journal of Geophysical Research: Earth Surface*, **113**(F1), p. 2007JF000769, doi: 10.1029/2007JF000769.

Patton, J., and B. Hornbuckle (2013), Initial Validation of SMOS Vegetation Optical Thickness in Iowa, *IEEE Geoscience and Remote Sensing Letters*, **10**(4), p. 647-651, doi:10.1109/LGRS.2012.2216498.

Pellarin, T., J.-P. Wigneron, J.-C. Calvet, M. Berger, H. Douville, P. Ferrazzoli, Y. H. Kerr, E. Lopez-Baeza, J. Pulliainen, L. P. Simmonds, et al. (2003), Two-year global simulation of L-band brightness temperatures over land, *IEEE Transactions on Geoscience and Remote Sensing*, **41**(9), p. 2135-2139.

Ploton, P., F. Mortier, N. Barbier, G. Cornu, M. Réjou-Méchain, V. Rossi, A. Alonso, J.-F. Bastin, N. Bayol, F. Bénédet, et al. (2020a), A map of african humid tropical forest aboveground biomass derived from management inventories, *Scientific Data*, **7**(1), p. 221.

Ploton, P., F. Mortier, M. Réjou-Méchain, N. Barbier, N. Picard, V. Rossi, C. Dormann, G. Cornu, G. Viennois, N. Bayol, A. Lyapustin, S. Gourlet-Fleury, and R. Pélissier (2020b), Spatial validation reveals poor predictive performance of large-scale ecological mapping models, *Nature Communications*, **11**(1), p. 4540, doi:10.1038/s41467-020-18321-y.

Qu, Y., J. Wang, H. Wan, X. Li, and G. Zhou (2008), A bayesian network algorithm for retrieving the characterization of land surface vegetation, *Remote Sensing of Environment*, **112**(3), p. 613-622, doi: http://dx.doi.org/10.1016/j.rse.2007.03.031.

Quan, X., B. He, M. Yebra, C. Yin, Z. Liao, X. Zhang, and X. Li (2017), A radiative transfer model-based method for the estimation of grassland aboveground biomass, *International Journal of Applied Earth Observation and Geoinformation*, **54**, p. 159-168, doi:https://doi.org/10.1016/j.jag.2016.10.002.

Rahmoune, R., P. Ferrazzoli, Y. K. Singh, Y. H. Kerr, P. Richaume, and A. Al Bitar (2014), Smos retrieval results over forests: Comparisons with independent measurements, *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, **7**(9), p. 3858-3866.

Robitaille, P.-M. (2007), On the Earth microwave background: absorption and scattering by the atmosphere, *Progress in Physics*, **3**(3), p. 3-4.

Rodríguez-Fernández, N. J., A. Mialon, S. Mermoz, A. Bouvet, P. Richaume, A. Al Bitar, A. Al-Yaari, M. Brandt, T. Kaminski, T. Le Toan, Y. H. Kerr, and J.-P. Wigneron (2018), An evaluation of smos I-band vegetation optical depth (I-vod) data sets: high sensitivity of I-vod to above-ground biomass in Africa, *Biogeosciences*, **15**(14), p. 4627-4645, doi:10.5194/bg-15-4627-2018.

Roy, D. P., J. Ju, P. Lewis, C. Schaaf, F. Gao, M. Hansen, and E. Lindquist (2008), Multi-temporal modis-landsat data fusion for relative radiometric normalization, gap filling, and prediction of landsat data, *Remote Sensing of Environment*, **112**(6), p. 3112-3130, doi:https://doi.org/10.1016/j.rse.2008.03.009.

Saatchi, S. S., N. L. Harris, S. Brown, M. Lefsky, E. T. A. Mitchard, W. Salas, B. R. Zutta, W. Buermann, S. L. Lewis, S. Hagen, S. Petrova, L. White, M. Silman, and A. Morel (2011), Benchmark map of forest carbon stocks in tropical regions across three continents, *Proceedings of the National Academy of Sciences*, **108**(24), p. 9899-9904, doi:10.1073/pnas.1019576108.

Santi, E., S. Paloscia, P. Pampaloni, and S. Pettinato (2009), Ground-based microwave investigations of forest plots in Italy, *IEEE Transactions on Geoscience and Remote Sensing*, **47**(9), p. 3016-3025, doi: 10.1109/TGRS.2009.2021613.

Santoro, M., and O. Cartus (2023), Esa biomass climate change initiative (biomass cci): Global datasets of forest above-ground biomass for the years 2010, 2017, 2018, 2019 and 2020, v4, NERC EDS Centre for Environmental Data Analysis, doi:10.5285/af60720c1e404a9e9d2c145d2b2ead4e.

Santoro, M., R. Lucas, H. Kay, and S. Quegan (2021), CCI BIOMASS PRODUCT USER GUIDE v3.

Schwank, M., and R. Naderpour (2018), Snow density and ground permittivity retrieved from I-band radiometry: Melting effects, *Remote Sensing*, **10**(2), p. 354.
Schwank, M., C. Matzler, M. Guglielmetti, and H. Fluhler (2005), L-band radiometer measurements of soil water under growing clover grass, *IEEE Transactions on Geoscience and Remote Sensing*, **43**(10), p. 2225–2237.

Schwank, M., R. Naderpour, and C. Mätzler (2018), Tau-omega and Two-Stream emission models used for passive L-band retrievals: Application to close-range measurements over a forest, *Remote Sensing*, **10**(12), p. 1868.

Schwank, M., A. Kontu, A. Mialon, R. Naderpour, D. Houtz, J. Lemmetyinen, K. Rautiainen, Q. Li, P. Richaume, Y. Kerr, et al. (2021), Temperature effects on L-band vegetation optical depth of a boreal forest, *Remote Sensing of Environment*, **263**, p. 112,542.

Seker, S., and A. Schneider (1988), Electromagnetic scattering from a dielectric cylinder of finite length, *IEEE Transactions on antennas and propagation*, **36**(2), p. 303–307.

Shabanov, N., D. Huang, W. Yang, B. Tan, Y. Knyazikhin, R. Myneni, D. Ahl, S. Gower, and A. Huete (2005), Analysis and optimization of the modis leaf area index algorithm retrievals over broadleaf forests, *IEEE Transactions on Geoscience and Remote Sensing*, **43**(8), p. 1855–1865.

Skakun, S., J. Wevers, C. Brockmann, G. Doxani, M. Aleksandrov, M. Batič, D. Frantz, F. Gascon, L. Gómez-Chova, O. Hagolle, D. López-Puigdollers, J. Louis, M. Lubej, G. Mateo-García, J. Osman, D. Peressutti, B. Pflug, J. Puc, R. Richter, J.-C. Roger, P. Scaramuzza, E. Vermote, N. Vesel, A. Zupanc, and L. Žust (2022), Cloud mask intercomparison exercise (cmix): An evaluation of cloud masking algorithms for landsat 8 and sentinel-2, *Remote Sensing of Environment*, **274**, p. 112,990, doi:https://doi.org/10.1016/j.rse.2022.112990.

Sulla-Menashe, D., and M. A. Friedl (), User Guide to Collection 6 MODIS Land Cover (MCD12Q1 and MCD12C1) Product.

Tagesson, T., R. Fensholt, I. Guiro, M. O. Rasmussen, S. Huber, C. Mbow, M. Garcia, S. Horion, I. Sandholt, and B. Holm-Rasmussen (2015), Ecosystem properties of semiarid savanna grassland in west Africa and its relationship with environmental variability, *Global change biology*, **21**(1), p. 250–264.

Teluguntla, P., P. Thenkabail, J. Xiong, M. Gumma, C. Giri, C. Milesi, M. Ozdogan, R. Congalton, J. Tilton, T. Sankey, R. Massey, A. Phalke, and K. Yadav (2016), NASA Making Earth System Data Records for Use in Research Environments (MEaSURES) Global Food Security Support Analysis Data (GFSAD) Crop Mask 2010 Global 1 km V001, doi:10.5067/MEASURES/GFSAD/GFSAD1KCM.001.

Teubner, I. E., M. Forkel, M. Jung, Y. Y. Liu, D. G. Miralles, R. Parinussa, R. Van Der Schalie, M. Vreugdenhil, C. R. Schwalm, G. Tramontana, G. Camps-Valls, and W. A. Dorigo (2018), Assessing the relationship between microwave vegetation optical depth and gross primary production, *International Journal of Applied Earth Observation and Geoinformation*, **65**, p. 79–91, doi:10.1016/j.jag.2017.10.006.

Tian, F., J.-P. Wigneron, P. Ciais, J. Chave, J. Ogée, J. Peñuelas, A. Ræbild, J.-C. Domec, X. Tong, M. Brandt, A. Mialon, N. Rodriguez-Fernandez, T. Tagesson, A. Al-Yaari, Y. Kerr, C. Chen, R. B. Myneni, W. Zhang, J. Ardö, and R. Fensholt (2018), Coupling of ecosystem-scale plant water storage and leaf phenology observed by satellite, *Nature Ecology & Evolution*, **2**(9), p. 1428–1435, doi:10.1038/s41559-018-0630-3.

Tong, X., F. Tian, M. Brandt, Y. Liu, W. Zhang, and R. Fensholt (2019), Trends of land surface phenology derived from passive microwave and optical remote sensing systems and associated drivers across the dry tropics 1992–2012, *Remote Sensing of Environment*, **232**, p. 111–307, doi:10.1016/j.rse.2019.111307.

Tucker, C. J. (1979), Red and photographic infrared linear combinations for monitoring vegetation, *Remote Sensing of Environment*, **8**(2), p. 127–150, doi:10.1016/0034-4257(79)90013-0.

Ulaby, F. T., and E. A. Wilson (1985), Microwave attenuation properties of vegetation canopies, *IEEE Transactions on Geoscience and Remote Sensing*, **GE-23**(5), p. 746–753.

Van de Griend, A. A., and J.-P. Wigneron (2004), The b-factor as a function of frequency and canopy type at h-polarization, *IEEE Transactions on Geoscience and Remote Sensing*, **42**(4), p. 786–794.