Final Report 10 YEARS OF SMOS – PASSIVE
MICROWAVE VEGETATION OPACITY STUDY
(PM-VO-S) OSMOSE

Arnaud Mialon, Simon Boitard, Nemesio Rodríguez-Fernández, Julio César Salazar-Neira, Beatrice Berthelot, Daria Malik, Mike Schwank, Yiwen Zhou, Frédéric Baret, Hongliang Ma, et al.

► To cite this version:

Arnaud Mialon, Simon Boitard, Nemesio Rodríguez-Fernández, Julio César Salazar-Neira, Beatrice Berthelot, et al. Final Report 10 YEARS OF SMOS – PASSIVE MICROWAVE VEGETATION OPACITY STUDY (PM-VO-S) OSMOSE. SO-TN-CB-GS-0119, CESBIO. 2024, pp.SO-TN-CB-GS-0119. hal-04840129

HAL Id: hal-04840129

https://hal.science/hal-04840129v1

Submitted on 16 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 10 YEARS OF SMOS – PASSIVE MICROWAVE VEGETATION OPACITY STUDY (PM-VO-S)
## OSMOSE

Contract No. 4000137990/22/NL/IA

Ref. number: SO-TN-CB-GS-0119

## Final Report

September 26, 2024

Prepared by:

<figure><img src="image_2.png" /></figure>
## Version History

<table><thead><tr><td>Issue</td><td>Date</td><td>Details</td><td>Editor</td></tr></thead><tbody><tr><td>1.0</td><td>26 Sept. 2024</td><td>version sent to ESA</td><td>Arnaud Mialon</td></tr><tr><td>0.1</td><td>21 Feb. 2024</td><td>Initial document template</td><td>Arnaud Mialon</td></tr><tr><td></td><td></td><td></td><td></td></tr></tbody></table>

Contacts: Arnaud Mialon (arnaud.mialon@univ-tlse3.fr) or Béatrice Berthelot (beatrice.berthelot@magellium.fr).

## Authors

<table><thead><tr><td>Partener</td><td>Author Name</td></tr></thead><tbody><tr><td>CESBIO</td><td>Arnaud Mialon<br>S. Boitard<br>N. Rodriguez-Fernandez<br>J.C. Salazar</td></tr><tr><td>Magellium</td><td>B. Berthelot<br>D. Malik</td></tr><tr><td>WSL</td><td>M. Schwank<br>Y. Zhou</td></tr><tr><td>EMMAH</td><td>F. Baret<br>H. Ma<br>M. Weiss</td></tr><tr><td>GlobEO</td><td>S. Mermoz</td></tr></tbody></table>

## Acknowledgements
# Contents

<table><tr><td>Part I: Project description</td><td>8</td></tr><tr><td>1 Introduction</td><td>8</td></tr><tr><td>2 Scientific background</td><td>9</td></tr><tr><td>3 OSMOSE Project overview</td><td>11</td></tr><tr><td>3.1 OSMOSE project objectives</td><td>11</td></tr><tr><td>3.2 OSMOSE project tasks</td><td>12</td></tr><tr><td>3.3 Project organization, workflow</td><td>12</td></tr><tr><td>3.4 Project schedule</td><td>14</td></tr><tr><td>4 Microwave Theoretical Analysis (WP 1)</td><td>16</td></tr><tr><td>4.1 Introduction</td><td>16</td></tr><tr><td>4.2 Methodology</td><td>16</td></tr><tr><td>4.2.1 Single-layer Discrete Scatter Model</td><td>16</td></tr><tr><td>4.2.2 Tree Structure Model (TSM)</td><td>19</td></tr><tr><td>4.2.3 Multi-layer Two Stream Model</td><td>24</td></tr><tr><td>4.2.4 Optimization</td><td>25</td></tr><tr><td>4.3 Results and Sensitivity Analysis</td><td>25</td></tr><tr><td>4.3.1 Effective scattering albedo and VOD based on the JRC structure</td><td>25</td></tr><tr><td>4.3.2 Sensitivity analysis</td><td>26</td></tr><tr><td>4.3.3 Impact of Tree Constituents and Their Dry-matter Fraction</td><td>27</td></tr><tr><td>4.3.4 Summary</td><td>29</td></tr><tr><td>4.4 Application in VOD retrieval</td><td>30</td></tr><tr><td>4.5 Conclusion and Future Work</td><td>31</td></tr><tr><td>5 Deriving VOD and AGB from passive microwaves (WP 2-3-4)</td><td>32</td></tr><tr><td>5.1 Introduction and purpose</td><td>32</td></tr><tr><td>5.2 Input data</td><td>32</td></tr><tr><td>5.2.1 General Considerations</td><td>32</td></tr><tr><td>5.2.2 SMOS</td><td>32</td></tr><tr><td>5.2.3 AMSRE</td><td>32</td></tr><tr><td>5.2.4 AMSR2</td><td>33</td></tr><tr><td>5.2.5 ERA 5</td><td>33</td></tr><tr><td>5.3 VOD Output Data file</td><td>34</td></tr><tr><td>5.4 Principle of the retrieval</td><td>36</td></tr><tr><td>5.4.1 Cost function minimization</td><td>37</td></tr><tr><td>5.4.2 Radiative Transfer Model</td><td>38</td></tr><tr><td>5.4.2.1 L-MEB</td><td>38</td></tr><tr><td>5.4.2.2 2S</td><td>38</td></tr><tr><td>5.4.2.3 Soil dielectric constant and emissivity</td><td>38</td></tr><tr><td>5.5 Operational configuration of the retrieval processing chains</td><td>38</td></tr><tr><td>5.6 VOD sensitivity analysis</td><td>40</td></tr><tr><td>5.6.1 Influence of the Radiative Transfer Model</td><td>40</td></tr><tr><td>5.6.2 Influence of the retrieval parameters</td><td>44</td></tr><tr><td>5.6.2.1 L-Band</td><td>44</td></tr></table>
<table><tbody><tr><td>5.6.2.2</td><td>C-Band</td><td>45</td></tr><tr><td>5.6.3</td><td>Influence of constraining the retrieval with upper and lower boundaries</td><td>45</td></tr><tr><td>5.6.4</td><td>Influence of the ω and Hr parametrization</td><td>46</td></tr><tr><td>5.6.5</td><td>Influence of retrieving one or two parameters</td><td>47</td></tr><tr><td>5.7</td><td>Validation of VOD</td><td>47</td></tr><tr><td>5.7.1</td><td>versus datacube2</td><td>47</td></tr><tr><td>5.7.1.1</td><td>L-Band</td><td>48</td></tr><tr><td>5.7.1.2</td><td>C-Band</td><td>50</td></tr><tr><td>5.7.1.3</td><td>X-Band</td><td>52</td></tr><tr><td>5.7.2</td><td>versus VODCA</td><td>53</td></tr><tr><td>5.7.2.1</td><td>L-band</td><td>54</td></tr><tr><td>5.7.2.2</td><td>C-band</td><td>54</td></tr><tr><td>5.7.2.3</td><td>X-band</td><td>55</td></tr><tr><td>5.7.3</td><td>Versus Wang X-VOD</td><td>57</td></tr><tr><td>5.8</td><td>AMSRE/AMSR2</td><td>57</td></tr><tr><td>5.8.1</td><td>Assessment of the AMSRE/AMSR2 bias for the C-band</td><td>57</td></tr><tr><td>5.8.2</td><td>Correction of the AMSRE/AMSR2 bias</td><td>58</td></tr><tr><td>5.8.3</td><td>Validation of the CDF matching</td><td>59</td></tr><tr><td>6</td><td>Analysis of AGB (WP5)</td><td>61</td></tr><tr><td>6.1</td><td>Introduction</td><td>61</td></tr><tr><td>6.2</td><td>Data</td><td>61</td></tr><tr><td>6.2.1</td><td>L-Band Level 3 Brightness Temperature and VOD/SM data-sets</td><td>61</td></tr><tr><td>6.2.2</td><td>C- and X-band VOD</td><td>62</td></tr><tr><td>6.2.3</td><td>Above Ground Biomass reference maps</td><td>62</td></tr><tr><td>6.2.4</td><td>Data for AGB evaluation</td><td>62</td></tr><tr><td>6.2.4.1</td><td>Lidar-based AGB data</td><td>62</td></tr><tr><td>6.2.4.2</td><td>In situ AGB data</td><td>63</td></tr><tr><td>6.2.4.3</td><td>FAO AGB data</td><td>63</td></tr><tr><td>6.2.4.4</td><td>AGB maps from Xu et al. (2021)</td><td>64</td></tr><tr><td>6.3</td><td>Methods</td><td>64</td></tr><tr><td>6.3.1</td><td>AGB estimation from multi-angular L-band TBs</td><td>64</td></tr><tr><td>6.3.2</td><td>AGB estimation from parametric function for L-band</td><td>65</td></tr><tr><td>6.3.3</td><td>AGB estimation from multi-frequency VOD</td><td>67</td></tr><tr><td>6.4</td><td>Results</td><td>67</td></tr><tr><td>6.4.1</td><td>Estimation of AGB from L-band brightness temperatures</td><td>67</td></tr><tr><td>6.4.1.1</td><td>Quality evaluation and predictors selection</td><td>67</td></tr><tr><td>6.4.1.2</td><td>Evaluation of AGB estimates from ANNs</td><td>69</td></tr><tr><td>6.4.1.3</td><td>Comparison of estimates from ANN and VOD fitted functions</td><td>72</td></tr><tr><td>6.4.2</td><td>L-AGB OSMOSE</td><td>72</td></tr><tr><td>6.4.3</td><td>AGB evaluation</td><td>74</td></tr><tr><td>6.4.3.1</td><td>AGB evaluation method</td><td>74</td></tr><tr><td>6.4.3.2</td><td>AGB evaluation results</td><td>75</td></tr><tr><td>6.4.4</td><td>Estimation of AGB from multi-frequency VODs</td><td>80</td></tr><tr><td>6.5</td><td>Conclusion</td><td>82</td></tr><tr><td>7</td><td>Optical indices (WP 6-7)</td><td>84</td></tr><tr><td>7.1</td><td>Introduction and purpose</td><td>84</td></tr><tr><td>7.2</td><td>Review of methods to estimate Vegetation traits from reflectance observations (WP6)</td><td>84</td></tr><tr><td>7.2.1</td><td>Biophysical parameter retrieval approaches from optical remote sensing</td><td>84</td></tr><tr><td>7.2.2</td><td>Pro and cons of the different approaches</td><td>86</td></tr><tr><td>7.2.2.1</td><td>Computation requirements</td><td>86</td></tr><tr><td>7.2.2.2</td><td>Flexibility of the observational configuration</td><td>86</td></tr><tr><td>7.2.2.3</td><td>Integration of prior information</td><td>86</td></tr><tr><td>7.2.2.4</td><td>Associated uncertainties</td><td>86</td></tr><tr><td>7.2.2.5</td><td>Robustness of the retrieval and quality assessment</td><td>87</td></tr></tbody></table>
<table><tbody><tr><td>7.2.3</td><td>Using VIs or reflectance for biophysical parameter retrieval?</td><td>87</td></tr><tr><td>7.3</td><td>Designing the algorithm for CWC retrieval from optical measurements</td><td>88</td></tr><tr><td>7.3.1</td><td>Material</td><td>88</td></tr><tr><td>7.3.1.1</td><td>Optical satellite data</td><td>88</td></tr><tr><td>7.3.1.2</td><td>TRY and OSSL databases</td><td>89</td></tr><tr><td>7.3.1.3</td><td>Ground CWC data</td><td>89</td></tr><tr><td>7.3.1.4</td><td>Selection of LANDVAL sites</td><td>90</td></tr><tr><td>7.3.2</td><td>Method</td><td>91</td></tr><tr><td>7.3.2.1</td><td>Setting up the PROSAIL training database</td><td>91</td></tr><tr><td>7.3.2.2</td><td>Training the machine learning algorithm</td><td>93</td></tr><tr><td>7.3.3</td><td>Evaluation of the CWC algorithm performances</td><td>93</td></tr><tr><td>7.3.3.1</td><td>Direct validation of satellite CWC estimates against ground measurements</td><td>93</td></tr><tr><td>7.3.3.2</td><td>Benchmarking Sentinel-2, Landsat-8 and MODIS CWC products</td><td>94</td></tr><tr><td>7.4</td><td>Conclusions</td><td>97</td></tr><tr><td>8</td><td>Optical Datacube (WP 8-9)</td><td>99</td></tr><tr><td>8.1</td><td>Introduction and purpose</td><td>99</td></tr><tr><td>8.2</td><td>Input/Output data</td><td>99</td></tr><tr><td>8.3</td><td>Datacube #1 reprojection algorithm</td><td>100</td></tr><tr><td>8.3.1</td><td>Convolution to the intermediate resolution</td><td>102</td></tr><tr><td>8.3.1.1</td><td>Special case: convolution of the non-physical qualitative variables</td><td>103</td></tr><tr><td>8.3.2</td><td>Reprojection with gdalwarp</td><td>104</td></tr><tr><td>8.4</td><td>Datacube #1 validation</td><td>105</td></tr><tr><td>8.4.1</td><td>Products with grid shared by source data and Datacube #1</td><td>106</td></tr><tr><td>8.4.1.1</td><td>Methodology</td><td>106</td></tr><tr><td>8.4.1.2</td><td>Results</td><td>106</td></tr><tr><td>8.4.2</td><td>Products with different grids and projections</td><td>110</td></tr><tr><td>8.4.2.1</td><td>Methodology</td><td>110</td></tr><tr><td>8.4.2.2</td><td>Results</td><td>112</td></tr><tr><td>8.4.2.3</td><td>Conclusion</td><td>116</td></tr><tr><td>8.5</td><td>Datacube #2 reprojection algorithm</td><td>116</td></tr><tr><td>8.5.1</td><td>Kernel computation</td><td>117</td></tr><tr><td>8.5.2</td><td>Convolution algorithm</td><td>118</td></tr><tr><td>8.5.3</td><td>Invalid/no data pixels contribution</td><td>119</td></tr><tr><td>8.5.4</td><td>Special case n°1: uncertainties/error variables</td><td>120</td></tr><tr><td>8.5.5</td><td>Special case n°2: quality flags</td><td>120</td></tr><tr><td>8.5.5.1</td><td>No projection between source data and Datacube 1</td><td>120</td></tr><tr><td>8.5.5.2</td><td>Projection between source data and Datacube 1</td><td>121</td></tr><tr><td>8.6</td><td>Datacube #2 validation</td><td>121</td></tr><tr><td>8.6.1</td><td>Methodologies</td><td>121</td></tr><tr><td>8.6.1.1</td><td>Global map comparison</td><td>121</td></tr><tr><td>8.6.1.2</td><td>Time series comparison</td><td>121</td></tr><tr><td>8.6.2</td><td>Results</td><td>122</td></tr><tr><td>8.6.3</td><td>Conclusion</td><td>131</td></tr><tr><td>9</td><td>Comparative Data Analysis (WP 10)</td><td>132</td></tr><tr><td>9.1</td><td>Introduction</td><td>132</td></tr><tr><td>9.2</td><td>Data</td><td>132</td></tr><tr><td>9.2.1</td><td>SMOS LVOD and SM</td><td>132</td></tr><tr><td>9.2.2</td><td>VODCA VOD</td><td>133</td></tr><tr><td>9.2.3</td><td>NDVI MODIS</td><td>133</td></tr><tr><td>9.2.4</td><td>Crop mask</td><td>133</td></tr><tr><td>9.2.5</td><td>Climate classification mask</td><td>134</td></tr><tr><td>9.2.6</td><td>IGBP land surface classification</td><td>134</td></tr><tr><td>9.3</td><td>Methods</td><td>134</td></tr><tr><td>9.3.1</td><td>Scaling and noise removal from time-series</td><td>134</td></tr></tbody></table>
<table><tbody><tr><td>9.3.2 Descriptors of crops hydraulics</td><td>135</td></tr><tr><td>9.3.3 Group Identification through K-Means Clustering</td><td>137</td></tr><tr><td>9.3.4 Attributes and cluster selection</td><td>137</td></tr><tr><td>9.3.5 Dynamic time warping</td><td>137</td></tr><tr><td>9.4 Results</td><td>138</td></tr><tr><td>9.4.1 Pixel Seasonality from pulse diagrams</td><td>138</td></tr><tr><td>9.4.2 Attribute selection</td><td>140</td></tr><tr><td>9.4.3 Unsupervised classification</td><td>140</td></tr><tr><td>9.4.4 Assessment of Unsupervised classification using external data</td><td>141</td></tr><tr><td>9.4.5 Dynamic time warping</td><td>142</td></tr><tr><td>9.5 Discussion</td><td>142</td></tr><tr><td>9.5.1 Cluster Interpretation</td><td>142</td></tr><tr><td>9.5.2 Temporal lags between PMWs and spectral indices</td><td>143</td></tr><tr><td>9.5.3 Physical interpretation of the vegetation parameters</td><td>144</td></tr><tr><td>9.5.4 Comparing to dynamic time warping</td><td>146</td></tr><tr><td>9.6 Conclusion</td><td>146</td></tr><tr><td><strong>10 Deliverables, progress reports, meetings</strong></td><td><strong>147</strong></td></tr><tr><td>10.1 Deliverables and reports</td><td>147</td></tr><tr><td>10.2 Project meetings</td><td>148</td></tr><tr><td><strong>11 Communication, outreach, and promotion</strong></td><td><strong>149</strong></td></tr><tr><td>11.1 Website</td><td>149</td></tr><tr><td>11.2 Outreach to other projects</td><td>150</td></tr><tr><td>11.3 Publications</td><td>150</td></tr><tr><td>11.4 Conference and workshop contributions</td><td>151</td></tr><tr><td><strong>12 Conclusions</strong></td><td><strong>152</strong></td></tr><tr><td><strong>References</strong></td><td><strong>153</strong></td></tr></tbody></table>
ACRONYMS

<table><tr><td>2S</td><td>Two Stream Radiative Transfer Model</td></tr><tr><td>AMSRE</td><td>Advanced Microwave Scanning Radiometer-EOS</td></tr><tr><td>AMSR2</td><td>Advanced Microwave Scanning Radiometer 2</td></tr><tr><td>ADF</td><td>Auxiliary Data File</td></tr><tr><td>ATBD</td><td>Algorithm Theoretical Baseline Document</td></tr><tr><td>CATDS</td><td>Centre Aval de Traitement des Données SMOS</td></tr><tr><td>CESBIO</td><td>Centre d'Etudes Spatiales de la Biosphère</td></tr><tr><td>DQX</td><td>Data Quality Index</td></tr><tr><td>EASE</td><td>Equal-Area Scalable Earth Grid</td></tr><tr><td>ECMWF</td><td>European Centre for Medium-range Weather Forecasting</td></tr><tr><td>FAPAR</td><td>Fraction of Absorbed Photosynthetically Active Radiation</td></tr><tr><td>IGBP</td><td>International Geosphere-Biosphere Programme</td></tr><tr><td>LAI</td><td>Leaf Area Index</td></tr><tr><td>L-MEB</td><td>L - band Microwave Emission of the Biosphere</td></tr><tr><td>MD</td><td>Cardioid Model</td></tr><tr><td>MEMLV</td><td>Microwave Emission Model for Layered Vegetation</td></tr><tr><td>NDVI</td><td>Normalized Difference Vegetation Index</td></tr><tr><td>RFI</td><td>Radiometric Frequency Interference</td></tr><tr><td>RTM</td><td>Radiative Transfer Model</td></tr><tr><td>SM</td><td>Soil Moisture</td></tr><tr><td>SMOS</td><td>Soil Moisture and Ocean Salinity</td></tr><tr><td>Tau/τ</td><td>Vegetation Optical thickness</td></tr><tr><td>TB</td><td>Brightness Temperature</td></tr><tr><td>VOD</td><td>Vegetation Optical Depth</td></tr><tr><td>WP</td><td>Work Package</td></tr></table>
# Part I: Project description
# Chapter 1

# Introduction

This technical note is the final report of the 10 YEARS OF SMOS – PASSIVE MICROWAVE VEGETATION OPACITY STUDY (PM-VO-S), which the consortium (CESBIO, Magellium, INRAE, WSL, GlobEO) carried out from June 2022 to September 2024, in the context of the ESA contract No. 4000137990/22/NL/IA. This final report describes the context and the activities done during these two years.

The first sections give the context of the project (Chapter 2 and Chapter 3). Chapter 4 then introduces an innovative model to describe better two essential microwave parameters representing the impact of vegetation that are τ and ω. Chapter 5 describes the method to derive the multi-frequency VOD, and an evaluation of the VOD. Chapter 6 details the Above Ground Biomass estimates obtained from the VOD. It also provides a validation of the derived AGB with in situ data. Chapter 7 starts the optical part, by giving the context and state of the art of the vegetation indices derived from optical data. Chapter 7 then describes a method of a new dataset which estimates the water content of the leaves components of the vegetation. Chapter 8 details the datacube 2 and the methodology to produce homogenized datasets. The analysis and the complementarity of all the indices to study the vegetation are analysed in Chapter 9. Chapter 12 concludes the project, and we lastly provide details of the project in terms of deliverables and meeting Chapter 10, outreach and communication Chapter 11.
# Chapter 2

# Scientific background

The thermal emission arising from the Earth surface at microwave frequencies depends on the soil characteristics such as its temperature, roughness and soil moisture content, which control the soil emissivity (Ulaby, 1976). In the presence of vegetation, part of the soil emission is absorbed and scattered. These effects can be parameterized using radiative transfer models such as the τ-ω model (Mo et al. (1982), Ulaby and Wilson (1985), Ferrazzoli and Guerriero (1996), Wigneron et al. (2007a), Liu et al. (2011)), where τ is the optical depth and ω is the single-scattering albedo. τ was shown to be linked to the vegetation water content (VWC, kg m⁻²) (Kirdiashev et al. (1979), Mo et al. (1982), Jackson and Schmugge (1991)) and to other vegetation properties such as the leaf area index (Jackson and Schmugge (1991), Van de Griend and Wigneron (2004), Wigneron et al. (2007a)). Therefore, τ is known as vegetation optical depth (VOD). VOD is also a function of the vegetation structure, which determines its dependence on the incidence angle and on the polarization of the radiation (Ulaby and Wilson (1985), Wigneron et al. (1995), Hornbuckle et al. (2003), Schwank et al. (2005)). Passive microwave radiometry is therefore a promising tool for monitoring vegetation on a global scale.

VOD samples the vegetation canopy, including woody vegetation, which uses root zone soil moisture (Andela et al. (2013)). VOD was used to study deforestation in South America (Van Marle et al. (2016)) and Africa (Brandt et al. (2018)). It has also been used to reveal teleconnections linking the state of the vegetation in Australia and El Niño Southern Oscillation (Liu et al. (2007)). In addition, Liu et al. (2015) showed the high potential of microwave VOD (X and C bands) to monitor the AGB dynamics on a large scale. Combining passive microwave VOD and optical data such as the NDVI, contributes to a more robust assessment of the vegetation characteristics (Liu et al. (2011)). The mentioned studies used VOD derived from different radiometers operating at different frequencies (Liu et al. (2011)) such as SSM/I at 19GHz (K-band), TRMM-TMI at 10.7 GHz (X-band) and the Advanced Microwave Scanning Radiometer - Earth Observing System (AMSR-E) at 10.7 and 6.9 GHz (C-band). It is worth noting that the VOD is a radiative parameter, and so intrinsically depends on the frequency of the electromagnetic radiation. VOD retrieved at different frequencies provide complementary information and should not be rescaled as if they were a single quantity as done by some studies (Liu et al. (2015); Moesinger et al. (2020a)). It is recommended that a specific VOD dataset is noted as B-VOD, where B stands for the microwave band (X-VOD, C-VOD, etc.). In principle, the lower the frequency, the lower the VOD for a given level of VWC (Wigneron et al. (1995, 2004); Ferrazzoli and Guerriero (1996)). Consequently, L-band (1.4 GHz, 21 cm) observations, which are less attenuated through the vegetation canopy, are capable of sampling the vegetation layer up to higher biomass values compared to higher-frequency observations.

SMOS retrieval algorithm (Kerr et al. (2012a, 2020)) is based on the L-MEB (L-band Model Emission of the Biosphere) transfer radiative model to compute the emission of the surface and its interactions with a vegetation layer (through the Tau Omega model). The full-polarization and multi-angular capabilities of SMOS allow the simultaneous retrieval of the soil moisture content and L-VOD. Lawrence et al. (2014a); Grant et al. (2016) compared SMOS L-VOD to X-VOD and C-VOD measured by AMSR-E and to visible/infrared vegetation indices. In crop zones, such as the MODIS vegetation indices, L-VOD increases during the growing season and decreases during senescence (Lawrence et al. (2014a)). On a global scale, L-VOD is less correlated
to optical/visible vegetation indices than X/C-VOD, suggesting that L-VOD can add more complementary information with respect to optical/infrared indices than X/C-VOD (Grant et al. (2016)). For instance, Rahmoune et al. (2014) found a significant correlation between L-VOD and tree height estimates. Vittucci et al. (2016) also discussed this relationship and compared it to the one estimated with X/C-VOD, which shows higher values for low tree-height than SMOS L-VOD, as expected. Vittucci et al. (2016) also showed a close to linear relationship between L-VOD and AGB at 20 selected points over Peru, Columbia and Panama. During the SMOS+Veg ESA project, Rodríguez-Fernández et al. (2018) showed that SMOS L-VOD could be used as a proxy to Above Ground Biomass (AGB estimates), by showing its higher sensitivity to AGB than optical indices and VOD derived at higher passive microwave frequencies. This work was extended globally by Mialon et al. (2020). Rodríguez-Fernández et al. (2018) showed the possibility of computing simple L-VOD to AGB relationships. Using this approach, subsequent studies such as those of Brandt et al. (2018) and Fan et al. (2019) were able to study the evolution of carbon stocks for the last ten years. The contribution of passive radiometry to above-ground carbon stocks monitoring arises from the fact that despite a low spatial resolution, they give access to the temporal evolution thanks to a global coverage with high revisit.

Considering that L-VOD is sensitive to both biomass and its water content, L-VOD was averaged on a monthly and yearly basis to smooth out the daily variability. However, considering the duration of the SMOS mission, it is possible to study the trends of AGB over the fifteen years of existing data. This is of significant interest, knowing that current AGB maps derived from radar and optical data are computed as yearly estimates and are not available regularly.

Estimating AGB from the VOD is however not an easy task, as the assumptions to derive VOD and the uncertainties associated with the AGB(VOD) relationship imply large uncertainties in terms of derived AGB. The present project takes advantage of the neural network approach to directly estimate AGB from passive microwave observations. The VODs and AGB from VODs are used in this project to evaluate the strengths of the new proposed approach. Considering that optical indices are linked to the upper part of the green vegetation, which differs from what is observed at microwave frequencies, the complementarity of these indices is evaluated to fully characterize the vegetation layer.

Even if SMOS provides the longest consistent time series of L-band VOD, many applications would require a simultaneous use of VOD at other frequencies. However, currently existing multi-frequency VOD data sets (VODCA, Moesinger et al. (2020a), SMOS-IC, AMSR-E/2, Wigneron et al. (2021); Wang et al. (2021a)) have been derived using different algorithms, or extending algorithms specifically designed for one frequency such as the L-Band Microwave Emission of the Biosphere, to other frequencies. Some of them include a rescaling of VOD obtained from different frequencies using CDF-matching (Moesinger et al. (2020a)), ignoring that VOD is intrinsically frequency dependent. Current analysis of multifrequency VOD data is a complex task because it is not clear that the observed behavior arises from artifacts due to non-harmonized data sets or from a geophysical origin (see for instance some examples in Sect. 2.1.3.4 from the ESA Land Carbon Constellation study). Therefore, generating a vegetation data set based on low frequency microwave brightness temperatures is needed, and should cover the 1.4 to 37 GHz frequency range using a generic algorithm well adapted to all the frequencies.

The acquisition of vegetation parameter datasets generated from measurements in the optical domain and processed to a datacube consistent with the data cube containing the microwave data is also necessary. First, because there is a large amount of scientific literature studying vegetation with optical indices, and second, because the information contained is complementary to that in VOD, in particular for L-VOD. Using these microwave and optical data cubes, it is possible to investigate the mechanisms driving the VOD at different frequencies and their relationships with the hydric state of the vegetation. The synergies of multi-frequency VODs and optical indices, also allow us to to perform better estimations of the Above Ground Biomass.
# Chapter 3

## OSMOSE Project overview

### 3.1 OSMOSE project objectives

Considering the context described in 2, the OSMOSE has the following objectives, presented in Table 3.1 which also shows the work packages covering each objective.

<table><caption>Figure 3.1: Objectives of the project and the corresponding work packages (WP)</caption><thead><tr><td>Objective</td><td>Description</td><td>WP</td></tr></thead><tbody><tr><td>010</td><td>Generate a vegetation data set based on low frequency microwave brightness temperatures covering the 1.4 to 37 GHz frequency range</td><td>1-4</td></tr><tr><td>015</td><td>Acquire and analyse a set of data records based on measurements obtained in the optical domain, preferably including Sentinel-3 OLCI and / or Proba-V measurements</td><td>6-9</td></tr><tr><td>020</td><td>Exploit the complementary information of the passive microwave and optical data sets. Spatial and temporal variabilities on the local scale, i.e. individual pixels, shall be included in the analyses. Mismatches in spatial support shall be discussed and minimized. Generate Above Ground Biomass from multi-frequency VOD data complemented by optical indices.</td><td>10<br>5</td></tr><tr><td>025</td><td>Provide a roadmap for a continuous generation of vegetation data products based on optical measurements and low frequency passive microwave measurements. Outline potential application areas for the synergistic use of the observations.</td><td>11</td></tr></tbody></table>
## 3.2 OSMOSE project tasks

The OSMOSE project was organized into twelve work packages to address these objectives, detailed in the next chapters. Overviews of the organization and activities are displayed in Figures 3.2, 3.3 and 3.4. These WP are entitled:

* **WP 0:** is the management and the website

* **WP 1:** Theoretical Analysis of multi-frequency brightness temperatures and VOD

* **WP 2:** Multi-mission Multi-frequency Passive Microwave VOD Data Sets- Algorithm Design Definition and Development

* **WP 3:** Multi-mission Passive Microwave VOD Data Sets - Input Data Preparation

* **WP 4:** Multi-mission Passive Microwave VOD Data Sets - Data Set Generation and Analysis

* **WP 5:** Derivation of Above Ground Biomass

* **WP 6:** Optical Vegetation Phenology / Indices - Literature Review

* **WP 7:** Multi-mission Phenology / Indices Data Set - Pre-processing Algorithm Design Definition and Development

* **WP 8:** Multi-mission Vegetation Phenology / Indices Data Set - InputData Set Preparation

* **WP 9:** Multi-mission Optical Vegetation Data Set - Data Cube Generation and Analyses

* **WP 10:** Comparative Data Analysis

* **WP 11:** Roadmap for Operational / continuous Data Generation and Exploitation

## 3.3 Project organization, workflow

Matthias Drusch -ESA- is the technical officer of this project, led by Arnaud Mialon (CESBIO) with the support of Beatrice Berthelot (Magellium). The different tasks are linked to each other so that the outputs of the activities benefit the other tasks.
<figure><img src="image_3.png" /><figcaption>Figure 3.2: Overview of the projet.</figcaption></figure>

Each task is divided into activities as shown in Figure 3.3.
<figure><img src="image_4.png" /><figcaption>Figure 3.3: Organization of tasks.</figcaption></figure>

## 3.4 Project schedule

The PM-VOS (OSMOSE) project was planned for two years starting in June 2022. Figure 3.4 shows the overall project schedule.
<figure><img src="image_5.png" /><figcaption>Figure 3.4: Gantt chart of the OSMOSE project, for the period June 2022-September 2024.</figcaption></figure>
# Chapter 4

# Microwave Theoretical Analysis (WP 1)

## 4.1 Introduction

The primary aim of WP1 is to develop a physics-based microwave emission model (MEM) that elucidates the dependence of effective radiative transfer properties on frequency and vegetation state parameters. The developed model is named as the Microwave Emission Model for Layered Vegetation (MEMLV). It is applied to coniferous forest as an exemplary study, covering a wide frequency range from P- to Ka-band (0.4 GHz – 37 GHz). The simulation can be used to support and guide spaceborne multi-frequency VOD retrieval in WPs 2-4. This is achieved by enhancing the parameterizations of the retrieval algorithm and deepening the understanding of the physics principles behind VOD retrieval.

## 4.2 Methodology

Physics-based forward models describe the fundamental mechanisms of the interaction between electromagnetic waves and vegetation canopies Fung and Chen (2010); Chauhan et al. (1994); Schwank et al. (2021), providing valuable insights for VOD retrieval. The developed MEMLV is a compound model that consists of three major modelling components: (i) the single-layer discrete model (SL-DSM) for calculating the single scattering albedo and VOD of a single-layer of vegetation canopies; (ii) a new Tree Structure Model (TSM) for representing the forest architecture in multi-layer structure and defining the characteristics of tree constituents in each layer, based on the concept of radio realism; (iii) the Two-Stream microwave emission model (2S-MEM) for combining SL-DSM and TSM to calculate the brightness temperature (BT) of forested ground. Finally, based on the simulated BT, an optimization process is employed to obtain the effective scattering albedo, $\omega_{\text{eff}}$, and the effective VOD, $\tau_{\text{eff}}$ of an equivalent single-layer canopy above ground. This process is a key to bridge the physics-based forward modelling with retrieval algorithms.

### 4.2.1 Single-layer Discrete Scatter Model

The discrete scatter model (DSM) Lang and Sidhu (1983); Chauhan et al. (1994); Zhou et al. (2020, 2023) employs a coherent approach to simulate the radiative transfer properties of canopies. It offers great flexibility in adopting various arrangements of vegetation constituents by representing vegetation as an effective medium with embedded discrete scatterers. The DSM has been employed in modelling the emission of a single layer canopy with below ground Kurum et al. (2012); Kurum (2013). However, this structure may not be sufficient to represent realistic coniferous forests.

In this WP, forest architectures are modeled using a multilayer structure to achieve a more realistic representation. Each layer contains tree constituents that have distinct number densities and size distributions. The single-layer DSM (SL-DSM) Chauhan et al. (1994); Zhou et al. (2020) (simplified DSM without considering the ground effect) is used to calculate the single scattering albedo and VOD of each layer.
<figure><img src="image_6.png" /><figcaption>Figure 4.1: Sketch of the geometry and coordinate system used in the SL-DSM based on an example of a vegetation layer that includes only one single class of scatterers.</figcaption></figure>

The essential part of the SL-DSM is the scatterer's scattering amplitude, $F$, which is defined as the ratio between the amplitude of the far-field scattered electric field and the amplitude of the incident electric field. Assuming a cylindrical scatterer with one end located at the origin of a Cartesian coordinate, a wave with an amplitude $E_0^p$ and direction along the unit propagation vector $k_i$ is incident on the scatterer, as shown in Fig. 4.1. The electric field of the incident wave is:

$$E_i^p(r) = E_0^p \exp(-jk_0 k_i \cdot r), \quad (4.1)$$

where $k_0 = 2\pi/\lambda$ represents the free-space wavenumber; $\lambda$ is wavelength; $\mathbf{r}$ is the vector from the origin to the field point. The variables $p$ and $q$ denote the orthogonal polarizations of the incident and scattered waves, respectively, which can be either horizontal (H) or vertical (V).

The incident wave is then scattered in all directions with different amplitudes. Along a specific direction $\mathbf{k}_s$, the scattered electric field can be written as Seker and Schneider (1988):

$$E_s^q(\mathbf{r}) = F^{(p,q)}(\mathbf{k}_s, \mathbf{k}_i) \cdot E_0^p \frac{\exp(-jk_0 \mathbf{k}_s \cdot \mathbf{r})}{|\mathbf{r}|}, \quad (4.2)$$

where $F^{(p,q)}(\mathbf{k}_s, \mathbf{k}_i)$ is the scattering amplitude specified by $p, q, \mathbf{k}_s$, and $\mathbf{k}_i$. If $\mathbf{k}_i$ and $\mathbf{k}_s$ are respectively defined by the associated zenith and azimuthal angles $(\theta_i, \phi_i)$ and $(\theta_s, \phi_s)$, and the orientation of an individual cylindrical scatterer is defined by $(\theta_c, \phi_c)$, then its scattering amplitude $F^{(p,q)}(\mathbf{k}_s, \mathbf{k}_i, f)$ at frequency $f$ can be expressed as $F^{(p,q)}(\theta_s, \phi_s; \theta_i, \phi_i; \theta_c, \phi_c; f)$. Fig. 4.1 depicts the geometry of a cloud of a specific class of scatterers in a layer. Hereinafter, a “class” of scatterers is defined for scatterers having the same dimensions (radius and length) and the same orientational probability density functions (PDFs) in a medium.

The directional patterns of the scattering amplitude are further investigated for a cylinder having radius $r = 1$ cm, length $l = 26$ cm, $m_d = 0.5$ kg·kg⁻¹, $S = 5$ ppt, $T = 20$°C, as an example. The cylinder is vertically oriented, and the incident wave is along the direction $\{\theta_i, \phi_i\} = \{40°, 0°\}$. The polarizations of the incident and scattered waves are assumed as $\{p, q\} = \{H, H\}$, respectively. Fig. 4.2 shows the 3D patterns of the scattering amplitudes at L-, C-, X- and Ka-band together with their projections. The brown vertical lines represent the cylinders, the red and blue arrows represent the directions of the incidence and forward scattering, respectively. The magnitude of the forward scattering amplitude, $|FF|$, is also given in the figure as a quantitative reference for the scale of plots. It is seen that at L-band, the specular- and forward scattering have comparable amplitudes. Small lobes are also present in the directions of the back- and bistatic scattering. At C-band, the forward scattering amplitude $FF$ becomes greater, and it dominates the patterns at higher frequencies (X- and Ka-band). This trend is due to the transition in scattering mechanism from Rayleigh scattering to Mie scattering as frequency increases Robitaille (2007).
<figure><img src="image_7.png" /><figcaption>Figure 4.2: Three-dimensional spatial patterns of scattering amplitudes $F^{(H,H)}$ simulated for vertically oriented cylinders ($\theta_c = 0^\circ$) illuminated by a wave incident at $\theta_i = 40^\circ$, $\phi_i = 0^\circ$. From left to right: L-band ($f = 1.4$ GHz), C-band ($f = 6.8$ GHz), X-band ($f = 11$ GHz), and Ka-band ($f = 37$ GHz). $|FF|$ indicates the absolute values of forward scattering amplitudes.</figcaption></figure>

Assuming the zenith and azimuthal orientational PDFs are $P_{\theta_c}$ and $P_{\phi_c}$ for a class of scatterers in the medium, respectively, the averaged scattering amplitude for this class of scatterers is:

$$ \langle F^{(p,q)}(\theta_s, \phi_s; \theta_i, \phi_i; f) \rangle = \int_{\theta_c=0}^{\pi} \int_{\phi_c=0}^{2\pi} F^{(p,q)}(\theta_s, \phi_s; \theta_i, \phi_i; \theta_c, \phi_c; f) P_{\theta_c} P_{\phi_c} d\theta_c d\phi_c, \quad (4.3) $$

where angular brackets $\langle\rangle$ denote ensemble average over the angular statistics for a certain class of scatterers. Consequently, the averaged scattering cross-section $\hat{\kappa}_s^p$ and extinction cross-section $\hat{\kappa}_e^p$ for this class of scatterers are Zhou et al. (2023); Kurum (2013):

$$ \hat{\kappa}_s^p(f) = \sum_{q=\{H,V\}} \int_{\theta_s=0}^{\pi} \int_{\phi_s=0}^{2\pi} \langle |F^{(p,q)}(\theta_s, \phi_s; \theta_i, \phi_i; f)|^2 \rangle \sin \theta_s d\theta_s d\phi_s \quad (4.4) $$

and

$$ \hat{\kappa}_e^p(f) = \frac{4\pi}{k_0} \cdot \mathrm{Im}\{\langle FF^{(p,p)} \rangle\} = \frac{4\pi}{k_0} \cdot \mathrm{Im}\{\langle F^{(p,p)}(\theta_{s,f}, \phi_{s,f}; \theta_i, \phi_i; f) \rangle\}, \quad (4.5) $$

where $FF^{(p,p)}$ is the forward scattering amplitude along the forward direction defined by $\theta_{s,f} = \pi - \theta_i$ and $\phi_{s,f} = \pi + \phi_i$. Hereinafter, the symbol $\hat{\cdot}$ is used for the radiative transfer properties of a class of scatterers with the same dimensions and orientational PDFs.

In general, the fundamental relationships shown in (4.3) – (4.5) can be applied to scatterers having any arbitrary shapes. Furthermore, these equations assume the direction of the incident wave is fixed while the scattered wave is in all directions. In passive radiometry, however, the waves incident on a scatterer are from all directions, while a radiometer attached to an ideal pencil beam antenna only senses the scattered wave in one fixed zenith direction, denoted by $\theta_{RM}$. These two conditions are essentially the same, as can be explained by the reciprocity of Maxwell's equations. Therefore, it is inherent that $\theta_i = \theta_{RM}$, as shown in Fig. 4.1.

To simulate a realistic vegetation layer, a composite of scatterer classes representing different vegetation constituents must be considered. Assuming the canopy consists of $i = 1, \dots, M$ classes of cylindrical scatterers with different radius $r_i$ and length $l_i$ and orientational PDFs, the scattering- and extinction coefficients $\kappa_s^p, \kappa_e^p$ of the resultant canopy layer can be obtained by Zhou et al. (2023):

$$ \kappa_s^p = \sum_{i=1}^{M} \rho_i \cdot \hat{\kappa}_{s,i}^p \quad \text{and} \quad \kappa_e^p = \sum_{i=1}^{M} \rho_i \cdot \hat{\kappa}_{e,i}^p, \quad (4.6) $$

where $\rho_i$ is the number density of the $i$-th class of vegetation constituents presented in the layer, and $\hat{\kappa}_{s,i}^p, \hat{\kappa}_{e,i}^p$ are its associated averaged scattering- and extinction cross-sections.
<figure><img src="image_8.png" /><figcaption>Figure 4.3: (a) Photo of the JRC tree in the anechoic chamber Fortuny and Sieber (1999). (b) Vectorization of the tree architecture Zhao (2013)</figcaption></figure>

Finally, the single scattering albedo ($\omega^p$) and the VOD ($\tau^p$) of a single layer with a thickness of $d$, composed of composite scatterers, are given by Kurum (2013):

$$ \omega^p = \frac{\kappa_s^p}{\kappa_e^p} \quad \text{and} \quad \tau^p = d \cdot \kappa_e^p, \qquad (4.7) $$

The scattering and absorption behavior of discrete scatterers over frequencies forms one of the fundamental elements of the MEMLV. As natural canopies consist of scatterers with different sizes, number densities and orientations PDFs, it is important to have a realistic tree structure to support the modelling of VOD and scattering albedo for coniferous forests. This will be discussed in the next section.

### 4.2.2 Tree Structure Model (TSM)

The development and parameterization of the TSM are inspired by an experiment conducted for a 5m-high Caucasian fir tree at the Joint Research Centre (JRC), Ispra, Italy Lang et al. (1998); Fortuny and Sieber (1999). The sample tree, referred to as the JRC tree, was placed in an anechoic chamber for radar measurements, and the photograph of the tree in the chamber Fortuny and Sieber (1999) is presented in Fig. 4.3(a). The structure of the JRC tree was measured, and a specific tree vectorization technique Landry et al. (1997) was employed to reconstruct the tree architecture, including the dimensions and numbers of the tree constituents, as shown in Fig. 4.3(b) Zhao (2013). Furthermore, the dielectric constants of the JRC tree constituents were measured and documented in Zhao (2013). These experimental data for the JRC tree are used as guidelines for parameterizing the TSM presented in the next section.

The TSM allows to discretely describe a variety of tree architectures based on simple equations. It assumes an individual tree is rotationally symmetric, and the discretization involves three aspects:

(i) A tree is stratified by $j = 0, 1, \dots, N$ horizontal layers from the bottom to the top. If the height of the lower bound of tree crowns at which the first branch emerges from the trunk is $H_{low}$, and each layer in tree crowns has a uniform thickness of $d$, then the mean height above the ground of the j-th layer, $h_j$, can be represented as:
<figure><img src="image_9.png" /><figcaption>Figure 4.4: Illustration of the TSM. (a) Panoramic forest structure with marked layer heights $h_j$, thickness $d$ and surface area $S_j$ (b) Side view of a tree with marked macro-parameters $H_{low}$, $H_{top}$, $\alpha$, $\beta$, $L_1$, together with $T_{B,j}$ emitted from stacks of $j$ layers and the total upwelling $T_{B,j}^{ML}$.</figcaption></figure>

$$h_j = \begin{cases} \frac{H_{low}}{2}, & \text{for } j = 0 \\ H_{low} + \frac{(j-1)}{2} \cdot d, & \text{for } j = 1, 2, \dots, N \end{cases} \qquad (4.8)$$

where $h_{j=0}$ is the mean height of the trunk section that is below the tree crown.

(ii) The dimensions of tree constituents are discretized by categorizing them into a limited number of classes, each represented by cylinders with specified radii $r_i$ and lengths $l_i$ ($i = 1, 2, \dots, M$). The index $i$ ascends as the size of the scatterer decreases.

(iii) The tree constituents are categorized into three types: trunk, branches and needles. Each type can comprise multiple classes of cylinders. The numbers of cylinder classes for trunk, branches, and needles are denoted as $M_t$, $M_b$, and $M_n$, respectively.

According to these discretization rules, four major components of the TSM are introduced below:

(i) **Tree macro structure:** The cross-sectional area of tree crowns is assumed to reduce as $h_j$ increases based on the structure of the JRC tree and observations of most natural trees in coniferous forests. Fig. 4.4 illustrates an exemplary three-layer tree crown structure, where each crown layer is represented by a square cuboid. The edge lengths of these cuboids are assumed to decrease linearly with a crown-shape parameter, $\alpha$, as $h_j$ increases. If the edge length of the bottom cuboid is denoted as $L_1$, then the edge length of the $j$-th layer is given by:

$$L_j = L_1 - \alpha \cdot (h_j - h_1) \qquad (4.9)$$

The TSM also includes the trunk layer, $j = 0$, located below the bottom layer of the crown, as illustrated in Fig. 4.4. In addition, the trunk is assumed to be also present within the tree crown layers, and its effect is simulated alongside branches and needles.

To realistically simulate the radiative transfer properties of coniferous forests, the canopy gaps (i.e. the area between trees without tree crown coverage) have been considered. The coefficient $\beta$ is introduced as the canopy gap fraction, representing the average distance between the centers of neighboring trees in forests as $(\beta + 1) \cdot L_1$, as shown in Fig. 4.4(b).
The parameters $\alpha$, $\beta$, and $L_1$ characterize the macro vertical structure and spatial arrangement of multi-layer trees in forests to model the upwelling BT, $T_B^{ML}$. As illustrated in Fig. 4.4, $T_B^{ML}$ is simulated as the weighted sum of BT emitted from the stack of j canopy layers above ground, $T_{B,j}$. Therefore, $T_B^{ML}$ can be written as:

$$T_B^{ML} = \frac{\sum_{j=0}^{N} T_{B,j} S_j}{L_1^2(1+\beta)^2}, \qquad (4.10)$$

with

$$S_j = \begin{cases} L_1^2(1+\beta)^2 - L_1^2, & j=0 \\ L_j^2 - L_{j+1}^2, & 0 < j < N \\ L_N^2, & j=N \end{cases} \qquad (4.11)$$

where $S_j/[L_1^2(1+\beta)^2]$ represents the surface area fraction of the stack of j canopy layers. Specifically, $S_0$ is the surface area of canopy gaps, as shown in Fig. 4.4(a).

(ii) **Trunks**: The bottom section of the trunk (cylinder class $i=1$) has the greatest radius and is considered as the sole constituent in layer $j=0$. Starting from the first layer of the tree crown ($j=1$), the trunk section has decreasing radii as $h_j$ increases, mixing with branches and needles. On the uppermost crown layer ($j=N$), it is assumed that the trunk is not present. Therefore, there are a total of $M_t=N$ classes of trunk sections in the forest canopy, and their radii can be written as:

$$r_i = r_1 - A \cdot (i-1) \cdot d, \quad \text{with } i=1, \dots, N, \qquad (4.12)$$

where A is a parameter adjusting the changing rate of the trunk radius as the aboveground height $h_j$ increases. The trunk sections are assumed to be mostly vertical.

(iii) **Branches**: The cylinder class index for branches is $i=N+1, \dots, N+M_b$ since $i=1, \dots, N$ accounts for the trunk. It is assumed that both the mean radius, $\bar{r}(h_j)$, and the radius-distribution PDF, $P_r(r, \bar{r}(h_j))$, of branches are dependent on $h_j$. Since $P_r(r, \bar{r}(h_j))$ is a PDF, the normalization condition $\int_0^{+\infty} P_r(r, \bar{r}(h_j)) dr = 1$ must hold to ensure the probability of finding a branch with any possible radius is 100% at a given $h_j$. As often used in the statistical representation of randomly distributed positive quantities, a log-normal function is used to represent $P_r(r, \bar{r}(h_j))$:

$$P_r(r, \bar{r}(h_j)) = \frac{1}{r \cdot \sigma(h_j) \cdot \sqrt{2\pi}} \cdot \exp \left\{ -\frac{1}{2} \left( \frac{\ln(r/\bar{r}(h_j))}{\sigma(h_j)} \right)^2 \right\}, \qquad (4.13)$$

where the height-dependent parameter, $\sigma(h_j)$, determines the shape of $P_r(r, \bar{r}(h_j))$, allowing for different radius distributions of branches at different $h_j$. As a log-normal function, $P_r(r, \bar{r}(h_j))$ peaks near $r = \bar{r}(h_j)$, takes the value zero at $r = 0$ cm, and converges to zero as $r \to \infty$. The smaller $\sigma(h_j)$ is, the steeper the peak of $P_r(r, \bar{r}(h_j))$ is, indicating an increased likelihood of the presence of a branch with $r \approx \bar{r}(h_j)$. For the $i^{th}$ class of cylindrical scatterers with a radius of $r_i$, its number density $\rho_i(h_j)$ at height $h_j$ can be obtained by integrating $P_r(r, \bar{r}(h_j))$ over the radius range $r_i - \Delta r_i \le r \le r_i + \Delta r_i$:

$$\rho_i(h_j) = B \cdot \int_{r_i - \Delta r_i}^{r_i + \Delta r_i} P_r(r, \bar{r}(h_j)) dr, \qquad (4.14)$$

where $i = N+1, \dots, N+M_b$, B is a parameter used to adjust the number density of all branches, and $\Delta r_i$ is the half-bin width of the integration range.

For natural coniferous trees, as $h_j$ increases, the crown is more likely to be dominated by small branches. Therefore, $\bar{r}(h_j)$ and $\sigma(h_j)$ are represented by functions that decrease with increasing $h_j$:

$$\bar{r}(h_j) = C - D \cdot h_j^{n_0}, \qquad (4.15)$$
with

$$C = \frac{H_{\text{top}}^{n_0} \cdot R_{\text{low}} - H_{\text{low}}^{n_0} \cdot R_{\text{top}}}{H_{\text{top}}^{n_0} - H_{\text{low}}^{n_0}},$$

$$D = \frac{R_{\text{low}} - R_{\text{top}}}{H_{\text{top}}^{n_0} - H_{\text{low}}^{n_0}},$$

and

$$\sigma(h_j) = E - F \cdot \frac{h_j - H_{\text{low}}}{H_{\text{top}} - H_{\text{low}}}, \qquad (4.16)$$

where $R_{\text{low}}$ and $R_{\text{top}}$ denote the mean radius of branches present in the lower-bound, $H_{\text{low}}$, and upper-bound, $H_{\text{top}}$, of the tree crown, respectively. The parameter $n_0$ allows the TSM to describe various patterns of how the mean radius $\bar{r}(h_j)$ of branches decreases with increasing $h_j$. In (4.16), $E$ and $F$ are parameters that adjust how $\sigma(h_j)$ changes with $h_j$.

The TSM also includes the statistical description of the branches' orientations expressed by the zenith and azimuthal angles $\theta_{c,b}$ and $\phi_{c,b}$ (Fig. 4.1). Similar to other constituents, branches are assumed to have an isotropic distribution $P_{\phi_{c,b}}$ in azimuth:

$$P_{\phi_{c,b}} = \frac{1}{2\pi} \quad \text{for} \quad 0 \le \phi_{c,b} \le 2\pi \qquad (4.17)$$

The zenith orientation PDF, $P_{\theta_{c,b}}$, is assumed to be uniform within a limited range of zenith angles $\theta_{c}$. Furthermore, it is assumed that the smaller the branch is, the wider the range of its zenith orientation. This reflects that larger branches are more likely to orient along a particular zenith direction, whereas smaller branches tend to have more random zenith orientations. Therefore, $P_{\theta_{c,b}}$ is represented by:

$$P_{\theta_{c,b}} = \begin{cases} (\theta_{\text{b-max}} - \theta_{\text{b-min}})^{-1}, & \text{if } \theta_{\text{b-min}} \le \theta_{c,b} \le \theta_{\text{b-max}} \\ 0, & \text{otherwise} \end{cases} \qquad (4.18)$$

with $\theta_{\text{b-min}} = \pi/2 - \pi/2 \cdot (G \cdot r + 1)^{-1}$ and $\theta_{\text{b-max}} = \pi/2 + \pi/2 \cdot (G \cdot r + 1)^{-1}$, where $G$ is the parameter to adjust the shape of the radius-dependence in $P_{\theta_{c,b}}$. For cylinders with $r \to 0$ cm, the limits of their zenith angle approach $\theta_{\text{b-min}} = 0$ and $\theta_{\text{b-max}} = \pi$, corresponding to isotropic zenith orientation.

(iv) **Needles**: Since the needles of mature coniferous trees have relatively constant dimensions, they are represented by a single class of cylinders. This results in their cylinder class index $i = M$. The orientational distribution of the needles is considered as isotropic. The number density of needles, $\rho_{i=M}$, is strongly associated with the number density of branches. For pine trees, often the needles only grow on small branches. Therefore, the branches are further classified into small branches and large branches. If the total number of classes for small branches is $M_{b,s}$, the number density of needles can be represented as:

$$\rho_M(h_j) = \sum_{i=M-1-M_{b,s}}^{M-1} (K \cdot \rho_i(h_j) \cdot l_i), \qquad (4.19)$$

where $K$ represents the average number of needles growing per unit length on small branches.

Equations (4.8) to (4.19) define the TSM following the concept of radio-realism. It employs a minimal number of observable parameters to describe a coniferous forest from the perspective of its macro-geometric structure to its fine distribution of tree constituents. Although here the TSM is developed and used for coniferous forests, its formulation makes it flexible to adapt to other types of forests and vegetation.

The TSM is parameterized such that resulting numbers of tree constituents closely match the in-situ measurements of the JRC tree (Section 4.2.2). The JRC tree has a height of ~5 m and a diameter of ~3.4 m Zhao (2013). To achieve a similar crown shape, the TSM is parameterized by $H_{\text{top}} = 5.5$ m, $H_{\text{low}} = 0.5$ m, $L_1 = 3.4$ m, and $\alpha = 0.6$. The number of layers and their thicknesses are assumed as $N = 5$ and $d = 1$ m, respectively. With these macro parameters, the volume of each layer can be calculated.
<table><caption>Table 4.1: Comparison of the cylinder classes $i = 1, 2, \dots, M = 15$ categorized in TSM with corresponding JRC-Tree experiments</caption><thead><tr><th>Cylinder classes</th><th></th><th>r (cm)</th><th>l (cm)</th><th>number W<sub>i</sub></th><th>r - JRC (cm)</th><th>number JRC</th></tr></thead><tbody><tr><td rowspan="5">Trunk</td><td>i=1</td><td>9</td><td></td><td></td><td rowspan="5">9 (DBH)<br/>3 (avg)</td><td rowspan="5">1</td></tr><tr><td>i=2</td><td>7</td><td></td><td></td></tr><tr><td>i=3</td><td>5</td><td rowspan="3">100</td><td rowspan="3">1</td></tr><tr><td>i=4</td><td>3</td></tr><tr><td>i=5</td><td>1</td></tr><tr><td rowspan="5">Large branch</td><td>i=6</td><td>1.0</td><td>30</td><td>19</td><td>0.9~1.0</td><td>23</td></tr><tr><td>i=7</td><td>0.8</td><td>26</td><td>52</td><td>0.7~0.9</td><td>63</td></tr><tr><td>i=8</td><td>0.6</td><td>22</td><td>170</td><td>0.5~0.7</td><td>173</td></tr><tr><td>i=9</td><td>0.45</td><td>19</td><td>231</td><td>0.4~0.5</td><td>165</td></tr><tr><td>i=10</td><td>0.35</td><td>16</td><td>540</td><td>0.3~0.4</td><td>621</td></tr><tr><td rowspan="5">Small Branch</td><td>i=11</td><td>0.25</td><td>13</td><td>1526</td><td>0.2~0.3</td><td>1199</td></tr><tr><td>i=12</td><td>0.17</td><td>10</td><td>2484</td><td rowspan="3">0.1~0.2</td><td rowspan="3">4854</td></tr><tr><td>i=13</td><td>0.12</td><td>8</td><td>3620</td></tr><tr><td>i=14</td><td>0.07</td><td>6</td><td>9786</td></tr><tr><td>Needle</td><td>i=15</td><td>0.03</td><td>2.5</td><td>2.8E6</td><td>0~0.1</td><td>11760</td></tr></tbody></table>

For trunk sections, it is assumed their radii $r_i = \{9, 7, 5, 3, 1\}$ cm for $i = \{1, 2, 3, 4, 5\}$, associated with layers $j = \{0, 1, 2, 3, 4\}$ counted from bottom to top, as shown in Table 4.1. The maximal radius, $r_1 = 9$ cm, is the same as the diameter-at-breast-height (DBH) of the JRC tree, while the average of $r_i$ for trunk sections is slightly larger than the average trunk radius of the JRC tree Zhao (2013). The trunk sections are assumed to be mostly vertical with a uniform orientation PDF within the range $0^\circ \le \theta_c \le 5^\circ$ and $0^\circ \le \phi_c \le 360^\circ$.

The branches of the JRC tree generally have radii from 0 to 1 cm. Accordingly, the TSM represents branches by nine different classes ($i = 6, \dots, 14$) of cylinders with radii $r_i$ from 0.07 cm to 1 cm, as shown in Table 4.1. Furthermore, branches with $r_i \ge 0.45$ cm are considered as large branches ($i = 6, \dots, 9$), while the rest are small branches ($i = 10, \dots, 14$).

According to the parameterization given in Appendix A, the number densities of these classes of branches are computed for each layer using (4.13) and (4.14). Based on the layer volumes, the total number of the $i^{th}$ class of scatterers, $W_i$, present in the entire tree can be obtained. The resulting $W_i$ are in good agreement with the in-situ measurements of the JRC tree, as shown in Table 4.1.

The in-situ measurements for the JRC tree also included branch lengths Zhao (2013). The results showed that the branch lengths in general range from 0 to 30 cm and are inversely proportional to the number of branches, following a similar distribution of the radius, see Table 4.1. However, the reported lengths and radii of the branches were decoupled, meaning for a branch with a specific radius $r_i$, its corresponding average length $l_i$ was not documented. Thus, the TSM assumes $l_i$ in accordance with the elastic similarity hypothesis McMahon and Kronauer (1976):

$$l_i = Q \cdot r_i^{2/3} \quad \text{with} \quad i = 6, 7, \dots, 14. \qquad (4.20)$$

The parameter $Q$ is found to be approximately 33 to confine $l_i$ within the documented range of branch lengths (0 – 30 cm). The corresponding $l_i$ and $r_i$ for branches are given in Table 4.1.

For needles, the TSM is parameterized such that the total number of needles $W_{i=15}$ matches the corresponding count $2.8 \times 10^6$ available from the in-situ measurements of the JRC tree (Table 4.1). The radius and length of needles are chosen as $r_{i=15} = 0.03$ cm and $l_{i=15} = 2.5$ cm based on Fortuny and Sieber (1999), which is confirmed by our own measurements made for a pine tree.

The dielectric constant of tree constituents, $\varepsilon_c = \varepsilon'_c - j \cdot \varepsilon''_c$, is modelled based on Matzler (1994):

$$\varepsilon_c(S, T, m_d, f) = 0.522 \cdot (1 - 1.32 \cdot m_d) \cdot \varepsilon_{sw}(S, T, f) + 0.51 + 3.84 \cdot m_d, \qquad (4.21)$$
<figure><img src="image_10.png" /><figcaption>Figure 4.5: Comparison between simulated frequency dependence $\varepsilon_c(f) = \varepsilon'_c(f) - j \cdot \varepsilon''_c(f)$ of different types of tree constituents and corresponding permittivities measured at $f = \{1.9, 5.5, 9.1\}$ GHz for the JRC tree.</figcaption></figure>

where $m_d$ (kg · kg⁻¹) is the dry-matter fraction and $\varepsilon_{sw}$ is the dielectric constant of sap water. It is assumed that $\varepsilon_{sw}$ can be represented by the dielectric constant of saline water, which is obtained by the Meissner-Wentz model Meissner et al. (2014) as a function of salinity $S$, temperature $T$ and frequency $f$.

Equation (4.21) is validated with the dielectric measurements made for the JRC tree at L- (1.9 GHz), C- (5.5 GHz), and X- (9.1 GHz) bands (Table 4.2 of *Zhao* (2013)). The $m_d$ and $S$ are adjusted (with $T = 20°C$) to match the modeled $\varepsilon_c$ with the measurements. This yields $m_d = 0.15 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 10$ ppt for needles; $m_d = 0.30 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 10$ ppt for branches and trunks. Figure 4.5 demonstrates good agreement between the modelled $\varepsilon_c$ and measured dielectric constants for the needles, branches, and trunks at the frequencies of the measurements.

### 4.2.3 Multi-layer Two Stream Model

According to the SL-DSM and TSM introduced in Section 4.2, the cross-sections $\hat{\kappa}_{s,i}^p, \hat{\kappa}_{e,i}^p$ and number densities $\rho_i$ of the $i^{th}$ ($i = 1, 2, ..., M$) scatterer in the $j^{th}$ ($j = 1, 2, ..., N$) layer can be obtained. Then the VOD and single scattering albedo for the $j^{th}$ layer, $\tau_j^p(\theta_{RM}, f)$ and $\omega_j^p(\theta_{RM}, f)$, are simulated for two polarizations $p = \{H, V\}$, four observation angles $\theta_{RM} = \{0°, 20°, 40°, 60°\}$, and frequencies over $0.4 \text{ GHz} \le f \le 37 \text{ GHz}$. Then they are fed into the ML-2S-MEM to simulate upwelling BT of multi-layer coniferous forests, $T_B^{p,ML}(\theta_{RM}, f)$, with the following primary input parameters:

(i) Transmissivities $t_j^p(\theta_{RM}, f)$ and volume reflectivities $r_j^p(\theta_{RM}, f)$ of vegetation layers $j = 0, ..., N$. They are computed from the aforementioned $\omega_j^p$ and $\tau_j^p$ using equation (18) in Schwank et al. (2018).

(ii) Interface reflectivities $s_j^p$ between vegetation layers. They are assumed as $s_j^p = 0$, meaning that the reflection between vegetation layers are neglected. In addition, the change of the wave propagation direction due to refraction at these interfaces are also neglected. These assumptions can be made without restrictions because tree crowns inherently form a continuous medium.

(iii) Interface reflectivity $s_G^p(\theta_{RM}, f)$ of the ground. It is computed from the effective ground dielectric constant $\varepsilon_G$ using Fresnel's equations with surface roughness correction. The dielectric mixing model Mironov et al. (2009) is used to express $\varepsilon_G$ by means of clay fraction (set to 0.16), ground temperature (set to $T_G = 293 \text{ K}$), and volumetric soil moisture (set to $SM = 0.3 \text{ m}^3 \cdot \text{m}^{-3}$). The HQN model Wang and Choudhury (1981) is used for roughness correction with the same empirical parameters $\{h_G, q_G, n_G^V, n_G^H\} = \{0.3, 0, 0, 2\}$ as used in SMOS retrievals over forested areas Wigneron et al. (2017).

(iv) Down-welling sky radiance $T_{sky}(\theta_{RM})$. It is estimated with the empirical approach Pellarin et al. (2003) that requires the observation angle $\theta_{RM}$, air temperature 2-m above ground (set to $T_{air} = 293 \text{ K}$), and elevation above sea level set to 0.4 km). Although the sky radiance model was developed for L-band ($f = 1.4 \text{ GHz}$), it is applied in this study to estimate $T_{sky}$ for microwave frequencies up to Ka-band ($f = 37 \text{ GHz}$). The error introduced by this assumption is minimal as the reflection at the upper bound of the canopy is generally small.
(v) Temperatures of the vegetation layers j. They are assumed as constant $T_j = 293$ K.

With the aforementioned input parameters, the BT emitted from the stack of j canopy layers above ground, $T_{B,j}^p(\theta_{RM}, f)$, can be computed. Using (4.10) and (4.11), the $T_B^{p,ML}(\theta_{RM}, f)$ can be obtained by summing $T_{B,j}^p(\theta_{RM}, f)$ weighted by the corresponding surface area fraction.

### 4.2.4 Optimization

Due to the large number of input parameters in the ML-2S-MEM, it is difficult to implement it directly in retrieval algorithms. The previous study Li et al. (2020) has demonstrated the feasibility of VOD retrieval using the single-layer (SL) 2S-MEM, which is the simplest form of ML-2S-MEM that requires much fewer input parameters by considering only one vegetation layer and the below ground (see Appendix B of Li et al. (2020)). Thus, the goal is to link the ML-2S-MEM with SL-2S-MEM so that the results derived from this research can facilitate retrieval practice. This is achieved by finding the effective VOD, $\tau_{eff}$, and the effective scattering albedo, $\omega_{eff}$, to minimize the cost function (CF) between $T_B^{p,ML}$ and the BT emitted from the equivalent single vegetation layer with beneath ground, $T_B^{p,SL}$:

$$CF = \sum_{\{p, \theta_{RM}\}} \left[ T_{B}^{p,ML}(\theta_{RM}, f) - T_{B}^{p,SL}(\tau_{eff}, \omega_{eff}, \theta_{RM}, f) \right]^2, \quad (4.22)$$

where $\theta_{RM} = \{0^\circ, 20^\circ, 40^\circ, 60^\circ\}$ and $p = \{H,V\}$.

Formally (4.22) corresponds to a two-parameter retrieval of $\{\tau_{eff}(f), \omega_{eff}(f)\}$ derived from the modeled $T_B^{p,ML}(\theta_{RM}, f)$ considered as "synthetic measurements" under ideal conditions (i.e., without noise and interference). The optimization process assumes $\tau_{eff}$ and $\omega_{eff}$ are at nadir and are only dependent on $f$, while their dependence on $p$ and $\theta_{RM}$ is not considered. This assumption not only simplifies the representation of $\tau_{eff}$ and $\omega_{eff}$, but also aligns with the current widely used {SM, VOD} retrieval algorithms Li et al. (2022); Fernandez-Moran et al. (2017) that only retrieves one VOD value at nadir for a pixel at a certain frequency according to a pre-defined scattering albedo. It should be noted that, however, the MEMLV can certainly include the polarization- and angular-dependence of $\tau_{eff}$ and $\omega_{eff}$ in the optimization process for future studies.

## 4.3 Results and Sensitivity Analysis

This section presents the results simulated using the MEMLV. First the results based on the original TSM parameterized by the JRC tree is presented. Then sensitivity analyses are shown for further evaluation of the impact of different parameters.

### 4.3.1 Effective scattering albedo and VOD based on the JRC structure

Figure 4.6(a) shows that as frequency goes higher, $\tau_{eff}$ increases most rapidly from L- to S-band (up to $f \approx 3$ GHz), then it increases almost linearly up to $f \approx 15$ GHz, finally it stays consistent from $f \approx 15$ GHz to $f \approx 37$ GHz. The $\tau_{eff}$ is in general greater than or similar to $\tau_j$ up to $f \approx 10$ GHz, then $\tau_{eff}$ becomes smaller than $\tau_j$ (especially for $f > 20$ GHz) due to the saturation effect. To the authors' knowledge, this is the first time that the saturation at high frequency has been quantitatively simulated using radiative transfer modeling.

To better understand this saturation effect, it is important to distinguish the physical meaning of $\tau_j$ from $\tau_{eff}$. The $\tau_j$ corresponds to the definition of the extinction of a medium with a certain thickness $d$, as expressed in (4.7). If the medium has a fixed extinction coefficient $\kappa_e$ (i.e. constant number density of scatterers), $\tau_j$ is proportional to $d$ and can become infinitely large when $d \to \infty$. In contrast, $\tau_{eff}$ is obtained via an optimization process that minimizes the difference between $T_B^{p,ML}$ and $T_B^{p,SL}$, as shown in (4.22). Above a certain frequency, when an increase in $\tau_j$ can no longer effectively change $T_B^{p,ML}$, the saturation of $\tau_{eff}$ occurs. Another reason for $\tau_{eff}$ being smaller than $\tau_j$ is that the simulation of $\tau_{eff}$ considers canopy gaps and tree crown shapes, while $\tau_j$ does not take these aspects into account.
<figure><img src="image_11.png" /><figcaption>Figure 4.6: Frequency spectra of the effective- and layer- radiative transfer properties simulated for a forest stand consisting of JRC trees. (a) VOD $\tau_{eff}$ and $\tau_{j}$; (b) scattering albedo $\omega_{eff}$ and $\omega_{j}$.</figcaption></figure>

Figure 4.6(b) shows that $\omega_{eff} \approx 0$ at P-band. It then increases rapidly, reaching a peak at S-band ($f \approx 3$ GHz). Following this peak, $\omega_{eff}$ decreases until the X-band ($f \approx 10$ GHz) and then steadily increases again at higher frequency. This peak is partially due to the frequency spectrum of the tree constituents' dielectric constant $\epsilon_{c}$ simulated with (8). As shown in Fig. 4.5, the imaginary part $\epsilon_{c}''$ has a minimum at approximately the same frequency as the peak of $\omega_{eff}$. Accordingly, absorption contributes less to the total extinction, resulting in an increased $\omega_{eff}$ at S-band. For $f > 10$ GHz, the increasing trend of $\omega_{eff}$ is due to the dominating contribution from needles, whose single scattering albedo increases monotonically with frequency.

The single scattering albedos of individual layers, $\omega_{j}$, shown in Fig. 4.6(b) have a relatively large dispersion when $f < 10$ GHz. The $\omega_{j}$ of the bottom crown layer ($j = 1$) has the greatest value, particularly when $f < 10$ GHz, because it has the largest fraction of large scatterers (trunks and large branches). As the layer gets higher in the canopy (increasing $j$), the fraction of the small scatterers (small branches and needles) increases, resulting in a reduced value of $\omega_{j}$. For $f < 10$ GHz, $\omega_{eff}$ demonstrates a combined effect of $\omega_{j}$ from all layers due to their higher transmissivity at this frequency range. For $10$ GHz < $f < 20$ GHz, $\omega_{eff}$ is closer to the $\omega_{j}$ of higher layers (e.g. $j = 3, 4, 5$) due to the reduced sensing depth. For $f > 20$ GHz, all $\omega_{j}$ become closer to each other, and $\omega_{eff}$ is almost the same as $\omega_{j}$ of the uppermost layer ($j = 5$) due to the weak transmissivity.

### 4.3.2 Sensitivity analysis

As explained in Section 4.2.2, the TSM was developed based on a discretization of the measured architecture of the 5m-high JRC tree (Fig. 4.3). The discretization considers $N = 5$ layers with thickness $d = 1$ m (Fig. 4.4). To investigate the impact of the tree height $H$, $d$ was scaled by a parameter $\gamma$. Here, it is defined as $\gamma = \{1, 2, 3, 4, 5\}$, resulting in the respective $d = \{1, 2, 3, 4, 5\}$m and $H = \{5, 10, 15, 20, 25\}$m.

According to (4.7), height scaling is reflected in scaling of $\tau_{j}^{p}$ of each layer $j$. To better simulate the realistic structure of trees in coniferous forests, the self-shading effect Kramer (2012) due to light competition has been considered. As $H$ increases, the lower layers of trees become more shaded by the upper layers of trees, resulting in a lower number density of tree constituents, $\rho_{i}$. Since $\tau_{j}^{p}$ is also dependent on $\rho_{i}$ according to (4.6) and (4.7), the scaling factor of $\tau_{j}^{p}$ for lower layers would be smaller than $\gamma$. Therefore, the layer-VOD scaling factor, $\gamma_{\tau}$, is defined as,

$$ \gamma_{\tau} = \begin{cases} \gamma, & j = 0 \\ \gamma - (\gamma - 1) \cdot (1 - j/N), & j = 1, 2, \dots, N \end{cases} \qquad (4.23) $$
For the trunk layer ($j = 0$), the $\gamma_{\tau}$ is the same as $\gamma$, as the self-shading effect has almost no effect on the trunk. When $\gamma = 1$, $\gamma_{\tau}$ is always equal to $\gamma$ for self-consistency, as observed in the original 5m JRC tree. As $\gamma$ increases, the VOD of the top layer ($j = N$) always has an equivalent increase, since the self-shading effect is minimal for this layer. When the layer becomes lower, $\gamma_{\tau}$ decreases but still remains greater than 1. For instance, when the tree height is 25m (scaled by $\gamma = 5$), the VOD of the bottom layer ($j = 1$) increases to only 1.8 times that of the original 5m tree due to the self-shading effect.

It is assumed that the height scaling has no impact on the scatterers' size distribution for each layer, resulting in that single scattering albedo $\omega_j^p$ remain unchanged. This assumption significantly simplifies the whole process.

The crown-shape parameter, $\alpha$, also affects $\omega_{eff}$ and $\tau_{eff}$. As shown in (4.9) and Fig. 4.4, a tree with a lower $\alpha$ would have a more cube-like crown, whereas with a higher $\alpha$, it would have a more pyramid-like crown. Because changing $\alpha$ in the TSM also changes the volume of each layer and the total number of all constituents, further adjustment has been made to $B$ and $n_0$ in (4.14) and (4.15) to obtain comparable results. These adjustments are only effective when $\alpha \neq 0.6$ to remain the self-consistency of the original TSM. When $\alpha$ decreases, the adjustments made in $B$ and $n_0$ lower the overall number density and slightly change the size-distribution of constituents, making the number of constituents contained within an individual tree similar to the JRC tree.

The resulting $\tau_{eff}$ and $\omega_{eff}$ for trees with $H = \{5, 10, 15, 20, 25\}$m and $\alpha = \{0.1, 0.6\}$ are plotted in Fig. 4.7. The insets show the $\tau_{eff}$ at P- and L-band (P-$\tau_{eff}$ and L-$\tau_{eff}$), and $\omega_{eff}$ at L-band (L-$\omega_{eff}$). In general, $\tau_{eff}$ increases and $\omega_{eff}$ decreases as $H$ increases. There is a clear saturation for $\tau_{eff}$ at $f > 15$ GHz and/or $H > 10$ m due to the limited sensing depth at high frequency. Likewise, the $H$-dependence of $\omega_{eff}$ gets weaker as frequency goes higher. Both the L-$\tau_{eff}$ and P-$\tau_{eff}$ have strong dynamics as $H$ changes. When $H$ changes from 5m to 25m, P-$\tau_{eff}$ almost doubled, while L-$\tau_{eff}$ increases by 50% – 70%. These quantities are dependent on how the VOD is scaled. For most cases, the P-$\omega_{eff}$ is negligible.

As $\alpha$ decreases, the air volume between tree crowns, and the number density of constituents within the volume of tree crowns decreases. In other words, the constituents are more evenly and sparsely distributed among the canopies. This causes a general increase for $\tau_{eff}$, as shown by the comparison between Fig. 4.7(a) and (b). Regarding the scattering albedo, it is interesting to see from Fig. 4.7(d) that L-$\omega_{eff}$ reduces significantly as $H$ increases. For this case ($\alpha = 0.1$), the tree crown has a cube-like shape and most of the large scatterers in lower layers are obscured by above canopy layers that contain predominantly small scatterers. When $H$ increases, the effect from these large scatterers in lower layers become weaker due to the increased thickness of upper layers. Since large scatterers have relatively high single scattering albedo at low frequencies, this obscuring effect causes an obvious decrease in L-$\omega_{eff}$ as $H$ increases. This effect is not very obvious for the case with a more pyramid-like crown ($\alpha = 0.6$) shown in Fig. 4.7(c). For this case, regardless of the canopy height, there is always a relatively large area for lower layers exposing directly to the radiometer due to the pyramid-like shape of crown. Therefore, the scattering effect from those large scatterers in lower layers persistently affect $\omega_{eff}$ when $H$ increases.

Besides $H$ and $\alpha$, the canopy-gap fraction $\beta$ poses another important macroscopic parameter. Fig. 4.8 shows $\omega_{eff}$ and $\tau_{eff}$ for a forest consisting of scaled JRC trees with $H = 15$ m when $\beta$ varies from 0.05 to 0.25. It is seen that the $\tau_{eff}$ in general decreases significantly as $\beta$ increases because of the reduced area covered by trees. Compared to $\tau_{eff}$, $\omega_{eff}$ demonstrates greater robustness to variations in $\beta$, with only a slight decrease in the magnitude of $\omega_{eff}$ as $\beta$ increases.

### 4.3.3 Impact of Tree Constituents and Their Dry-matter Fraction

It is known that the VOD correlated with VWC. The impact of dry-matter fraction and sap-water salinity on $\tau_{eff}$ and $\omega_{eff}$ is evaluated and plotted in Fig. 4.9(a) and (c). The evaluation is made for a 15m-high forest for the “normal” and “dry” cases. The “normal” case follows the parameterizations used in Section 4.2.2: $m_d = 0.15 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 10 \text{ ppt}$ for needles, and $m_d = 0.3 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 10 \text{ ppt}$ for branches and trunks. The “dry” case assumes $m_d = 0.3 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 10 \text{ ppt}$ for needles, $m_d = 0.4 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 5 \text{ ppt}$ for branches, and $m_d = 0.7 \text{ kg} \cdot \text{kg}^{-1}$ and $S = 5 \text{ ppt}$ for trunks.

The effects of $m_d$ on $\tau_{eff}$ and $\omega_{eff}$ are shown in Fig. 4.9(a) and (c), respectively. The dry case has a lower $\tau_{eff}$ than the normal case up to $f \approx 13$ GHz. Specifically, the P-$\tau_{eff}$ and L-$\tau_{eff}$ for the dry case are about 35%
<figure><img src="image_12.png" /><figcaption>Figure 4.7: Frequency spectra of the effective radiative transfer properties of coniferous forests with scaled heights $H = \{5, 10, 15, 20, 25\}$ m and crown-shape parameter $\alpha = \{0.1, 0.6\}$. (a) $\tau_{\text{eff}}$- and (c) $\omega_{\text{eff}}$ for $\alpha = 0.6$; (b) $\tau_{\text{eff}}$- and (d) $\omega_{\text{eff}}$ for $\alpha = 0.1$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 4.8: Frequency spectra of the effective radiative transfer properties for coniferous forests with height $H = 15$ m and canopy-gap fractions $\beta = \{0.05, 0.1, 0.15, 0.2, 0.25\}$ (a) for $\tau_{\text{eff}}$; (b) for $\omega_{\text{eff}}$, where the representation of line colors is indicated in the legend of (a).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 4.9: Frequency spectra of the effective radiative transfer properties for coniferous forests with height $H = 15$ m under different conditions. (a) $\tau_{\text{eff}}$- and (c) $\omega_{\text{eff}}$ for “normal” and “dry” canopy; (b) $\tau_{\text{eff}}$- and (d) $\omega_{\text{eff}}$ for “normal” and “defoliated” canopy.</figcaption></figure>

and 25% lower than for the normal case, respectively. Compared with the normal case, the $\tau_{\text{eff}}$ of the dry case reaches saturation at higher frequency. After reaching saturation, $\tau_{\text{eff}}$ for the dry case is slightly larger (< 10%) than for the normal case. The comparison between $\omega_{\text{eff}}$ for the normal and dry cases follows the opposite trend. For the dry case, $\omega_{\text{eff}}$ is slightly higher at low frequency but lower at high frequency. At low frequency, when $m_d$ increases, the single scattering albedo $\hat{w}^p$ of branches and trunks increases, leading to an increase in $\omega_{\text{eff}}$ for the dry case. As frequency increases, the radiative transfer contribution from needles becomes dominated. The $\hat{w}^p$ of dry needles has lower increasing rate along $f$. This explains why $\omega_{\text{eff}}$ could become lower for the dry case at very high frequencies.

The effects of defoliation on $\tau_{\text{eff}}$ and $\omega_{\text{eff}}$ are shown in Fig. 4.9(b) and (d), respectively. For this analysis, the normal case is the same as before, while the “defoliated” case assumes the number density of needles as 0. The results show that both the P-$\tau_{\text{eff}}$ and L-$\tau_{\text{eff}}$ for the defoliated case are smaller than for the normal case. The impact of leaves on L-$\tau_{\text{eff}}$ is minor, agreeing with other studies Guglielmetti et al. (2007); Santi et al. (2009). From X- to Ku-band, $\tau_{\text{eff}}$ for the normal case is significantly higher than for the defoliated case, highlighting the suitability of these bands for detecting the defoliation states of forests. It is interesting to see that $\tau_{\text{eff}}$ of the defoliated case reaches saturation at much higher frequency (≈ 25 GHz) compared with the normal case (≈ 13 GHz). Since needles in general have small single scattering albedo, $\omega_{\text{eff}}$ for the defoliated case is higher than for the normal case, particularly at low frequencies.

### 4.3.4 Summary

Simulations demonstrate the dependence of $\tau_{\text{eff}}$ and $\omega_{\text{eff}}$ on the macro-parameters of coniferous forests (e.g., canopy height and -gap fraction), the structure of individual trees (e.g., crown shape), and the tree constituents (e.g., foliation stages and dry-matter fraction). The impacts of these factors are summarized in Table 4.2, where the sign “↑” indicates increasing, “↓” means decreasing, and “~” suggests staying relatively consistent. Table 4.2 offers a qualitative estimation of the dependence of $\tau_{\text{eff}}$ and $\omega_{\text{eff}}$ on different ancillary parameters to
<table><caption>Table 4.2: Summary of the impact from different parameters on the effective VOD and albedo</caption><thead><tr><th rowspan="2">f (GHz)</th><th colspan="3">ω<sub>eff</sub>(f)</th><th colspan="3">τ<sub>eff</sub>(f)</th></tr><tr><th>1-5</th><th>5-15</th><th>15-37</th><th>1-5</th><th>5-15</th><th>15-37</th></tr></thead><tbody><tr><td>Tree crown slope α ↑</td><td>↑</td><td>↑</td><td>~</td><td>↓</td><td>↑</td><td>↑</td></tr><tr><td>Tree height H ↑</td><td>↓</td><td>~</td><td>~</td><td>↑</td><td>↑</td><td>~</td></tr><tr><td>Canopy-gap fraction β ↑</td><td>↓</td><td>↓</td><td>↓</td><td>↓</td><td>↓</td><td>↓</td></tr><tr><td>Leaf number/volume ↑</td><td>↓</td><td>↓</td><td>↓</td><td>↑</td><td>↑</td><td>~</td></tr><tr><td>Dry-mat. m<sub>d</sub>↑ & Sal. S ↓</td><td>↑</td><td>~</td><td>↓</td><td>↓</td><td>~</td><td>~</td></tr></tbody></table>

<table><caption>Table 4.3: An example of effective VOD and scattering albedo for a forest stand at P- (0.4 GHz), L-(1.4 GHz), C-(6.8 GHz), X-(11 GHz), Ku-(19 GHz), and Ka-(37 GHz) bands.</caption><thead><tr><th></th><th>P-τ<sub>eff</sub> (P-ω<sub>eff</sub>)</th><th>L-τ<sub>eff</sub> (L-ω<sub>eff</sub>)</th><th>C-τ<sub>eff</sub> (C-ω<sub>eff</sub>)</th><th>X-τ<sub>eff</sub> (X-ω<sub>eff</sub>)</th><th>Ku-τ<sub>eff</sub> (Ku-ω<sub>eff</sub>)</th><th>Ka-τ<sub>eff</sub> (Ka-ω<sub>eff</sub>)</th></tr></thead><tbody><tr><td>5m</td><td>0.210 (0.000)</td><td>0.334 (0.194)</td><td>0.680 (0.196)</td><td>0.969 (0.181)</td><td>1.120 (0.212)</td><td>1.100 (0.294)</td></tr><tr><td>10m</td><td>0.262 (0.000)</td><td>0.381 (0.166)</td><td>0.744 (0.197)</td><td>1.030 (0.184)</td><td>1.136 (0.212)</td><td></td></tr><tr><td>15m</td><td>0.308 (0.000)</td><td>0.420 (0.150)</td><td>0.784 (0.199)</td><td>1.060 (0.185)</td><td>1.141 (0.212)</td><td></td></tr><tr><td>20m</td><td>0.351 (0.000)</td><td>0.454 (0.139)</td><td>0.810 (0.200)</td><td>1.075 (0.186)</td><td></td><td>1.103 (0.294)</td></tr><tr><td>25m</td><td>0.390 (0.000)</td><td>0.484 (0.133)</td><td>0.826 (0.201)</td><td>1.083 (0.186)</td><td></td><td></td></tr></tbody></table>

support improving the retrieval algorithms.

Quantitatively, Table 4.3 provides exemplary values of $τ_{eff}$ and $ω_{eff}$ for different frequencies $f$ and tree heights $H$ with a fixed canopy-gap fraction $β = 0.05$. The results further corroborate the conclusions drawn previously. It is seen that at C-band and higher frequencies, the $H$-dependence of $ω_{eff}$ and $τ_{eff}$ is very weak due to the reduced sensing depth. The frequency and tree height at which this effect becomes obvious depend on forest structures and dry-matter fraction, as demonstrated in Figs. 4.7-4.9.

## 4.4 Application in VOD retrieval

There are numerous possibilities of applying the MEMLV to enhance the multi-frequency VOD retrieval. The major aspects are outlined as follows:

### (i) Parameterization of retrieval algorithms:

The sensitivity study summarized in Table 4.2 provides guidelines for adjusting the parameterization of scattering albedo in retrieval based on canopy characteristics. The values of $ω_{eff}$ in Table 4.3 can facilitate the retrieval of variables such as VOD through the inversion of 2S-MEM. Although their direct applicability may be affected by natural forest variability, their frequency- and tree height-dependent trends are expected to hold for natural forest stands. Consequently, VOD retrieval could be optimized by using scaled $ω_{eff}$ values while maintaining the trends of frequency- and height-dependence observed in Table 4.3. This approach offers a more dynamic and representative parameterization strategy for retrieval algorithm, potentially in combination with data from other sensors. For instance: (1) The $ω_{eff}$ of forests can be adjusted based on tree height $H$ measured by LiDAR or InSAR, along with canopy-gap fraction $β$ and tree crown-shape parameter $α$ estimated using image pattern recognition techniques. This will enable data fusions from the modelling level to improve estimations of $ω_{eff}$ and enhance VOD retrieval for forests with various characteristics. (2) Seasonal variations in tree constituents could be factored into parameterizing $ω_{eff}$ for VOD retrieval. For example, defoliation or reduced liquid sap water during cold periods may result in higher $ω_{eff}$, particularly at low microwave frequencies. Conversely, during warming periods, one might expect a decrease in $ω_{eff}$ due to growth and associated increasing gravimetric water content of canopy constituents.
## (ii) Interpretation of retrieved VOD:

Properly interpreting the retrieved VOD is often challenging due to limited knowledge about the physical interactions between canopies and electromagnetic waves at various frequencies. The simulated $τ_{eff}$ using the MEMLV provides physics-based explanations, justifies assumptions and elucidates specific behaviors of the retrieved VOD. Furthermore, the simulated $τ_{eff}$ helps establish the relationships between VOD and the architecture and physiology of canopies, as presented in Figures 4.7-4/9. Consequently, this facilitates the quantification of parameters such as AGB and VWC based on the retrieved VOD.

## (iii) Definition of future campaigns and missions:

The crucial preparatory phase before launching microwave campaigns and missions requires the definitions of their scientific needs, advantages and goals. The MEMLV enables sensitivity analysis of microwave measurements across different frequencies and environmental conditions to support in such definitions. For instance, our simulations indicated a negligible P-$ω_{eff}$ for coniferous forests (Table 4.3). This finding is an advantageous for future P-band missions aimed at quantifying forest state parameters, as it simplifies the retrieval process and reduces uncertainties involved in the retrieval scheme.

# 4.5 Conclusion and Future Work

A new compound modelling tool, the MEMLV, is developed and presented along with sensitivity analyses for coniferous forests. It is designed as a toolbox to simulate the effective radiative transfer parameters across frequencies of land covered with different types of vegetation under a variety of conditions. The primary aim of the MEMLV is to bridge the gap between highly detailed physics-based forward emission models and retrieval processes over a wide frequency range. To achieve this goal, the MEMLV incorporates three major components:

(i) the SL-DSM for calculating the single scattering albedo and VOD of a single-layer of vegetation canopies;

(ii) the TSM for representing the forest architecture in multi-layer structure and defining the characteristics of tree constituents in each layer;

(iii) the 2S-MEM for calculating the BT of forested ground and obtaining the effective radiative transfer properties that can be used in retrieval practices.

The principal outcomes of this research are the frequency spectra of effective VOD, $τ_{eff}$, and effective scattering albedo, $ω_{eff}$. The $τ_{eff}$ generally increases as frequency increases but can saturate at high frequency. The $ω_{eff}$ increases rapidly at low frequency and peaks at around S-band, then it decreases until X-band (10 GHz), and finally, it increases almost linearly as frequency further increases.

Another major finding concerns the sensitivity analyses of $τ_{eff}$ and $ω_{eff}$ with respect to crucial model parameters. They are summarized below:

(i) Canopy height $H$: at low frequencies, $ω_{eff}$ slightly decreases and $τ_{eff}$ increases when $H$ increases. At high frequencies, $τ_{eff}$ and $ω_{eff}$ show minimal variation with $H$ due to saturation.

(ii) Canopy-gap fraction $β$: when $β$ increases, $ω_{eff}$ slightly decreases while $τ_{eff}$ significantly decreases.

(iii) Dry-matter fraction $m_d$: as $m_d$ increases, $ω_{eff}$ increases while $τ_{eff}$ decreases at low frequency.

(iv) Defoliation: it can cause an increase in $ω_{eff}$ across all frequencies, a slightly lower $τ_{eff}$ at P- and L-band, and a significant lower $τ_{eff}$ from X- to Ku-band.

Although the simulated radiative transfer properties have shown encouraging results when compared with satellite retrieval and literature values, rigorous ground-based microwave experiments are essential to fully validate MEMLV. Moving forward, it is crucial to perform close-range multi-frequency radiometer experiments on forest sites with well-documented ground-truth data. Additionally, experiments conducted under controlled environmental conditions (e.g., irrigated vs. non-irrigated) and across different growth seasons would be also very valuable for thoroughly assessing the performance of the MEMLV. In the long term, the MEMLV is anticipated to serve as a universal tool for the remote sensing community, constructively supporting and improving retrieval practices for satellites operating across frequencies from P- to Ka-band.
# Chapter 5

# Deriving VOD and AGB from passive microwaves (WP 2-3-4)

## 5.1 Introduction and purpose

The purpose of this part, is to detail the setup to derive VOD, the retrieval chains and their configurations. Data from the SMOS, AMSRE and AMSR2 sensors are used to derive the VOD at L, C and X bands. The following represents the activities developed for the tasks (work packages) 2, 3 and 4 of the project (see 3).

## 5.2 Input data

### 5.2.1 General Considerations

All the products discussed in the following are provided for the OSMOSE project as:

* **Format**: netcdf files, characterized by the .nc extension name

* **Grid system**: EASE Grid version 2 (Brodzik et al. (2012a), Brodzik et al. (2014)), that is 25km x 25km at 30deg of latitude

* **Grid dimensions**: global earth is covered by a matrix of 584x1388 pixels

### 5.2.2 SMOS

The Level 3 Brightness Temperature product (Al Bitar et al. (2017a)) is a daily global polarised brightness temperature product, arranged by incidence angle values, in full polarisations (H, V, stokes3, and stockes4). It includes all brightness temperatures acquired that day, transformed to ground polarisation reference frame (H and V polarisation), binned (5°) and averaged into fixed angle classes. Bin 9 corresponds to SMAP incidence angle (40° +/- 2°), and bin 15 is a placeholder for the time being. Ascending and descending orbits are processed separately.

### 5.2.3 AMSRE

The Advanced Microwave Scanning Radiometer for EOS (AMSR-E) aboard NASA's Aqua satellite, provides Brightness Temperatures at Vertical and Horizontal polarisations at an incidence angle of 55 deg, for the 6
frequencies 89.0 GHz, 36.5 GHz, 23.8 GHz, 18.7 GHz, 10.65 GHz and 6.925 GHz. The OSMOSE project focuses on the C and X bands, so only the 6.925 GHz and 10.65 GHz are kept.

* Radiometric accuracy: <0.5 K

* Spatial Resolution: 51 x 29 km (10.65 GHz, X-band), 75 x 43 km (6.925 GHz, C-band)

* Incidence Angle: 55 deg.

The purpose of the OSMOSE project is to set up an homogenized VOD data set, on the EASE Grid version 2 (25 km). The NSIDC provides the AMSRE data on such grid, called the "MEaSUREs dataset" (Brodzik et al. (2016), Brodzik and Long (2018)), but the data are split into two products, one on the global projection, but only covering latitudes lower than 70 deg. and one on the polar grid. The two projections are then combined daily (ascending/descending orbits processed separately) to ensure a global coverage, i.e. including the high latitudes. The MEaSUREs database encompasses different passive microwave sensors (e.g. SSM/I, SSMR...), which are inter-calibrated and for which the projection on the EASE Grid 2 is properly done, in case the methods developed for the OSMOSE project can be applied to the other passive microwave satellites.

### 5.2.4 AMSR2

The AMSR2 sensor is onboard the GCOM-W1 satellite, launched in July 2012 by the JAXA. It provides brightness temperatures at Vertical and Horizontal polarization, at an incidence angle of ~55°deg.. The data are split between Ascending/Descending orbits, corresponding to an overpassing time of 13:30/01:30, Maeda et al. (2016).

* Name of dataset: LEVEL3

* original spatial resolution: 0.25 deg

* Website: https://suzaku.eorc.jaxa.jp/

### 5.2.5 ERA 5

ERA 5 is the re-analysis of climate data from ECMWF, which provides various atmospheric and surface geophysical information every hour, on a 0.25deg spatial resolution grid Hersbach et al. (2023). For the OSMOSE project, we used a subset of the catalog available from the website https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=form, described in the next section. The ERA5 are pre-processed for the applications to 1/ use only the closest time slot to the satellite overpassing time. It means for each satellite used in this project, and for each orbit, the hourly ERA5 data may change.

Other auxiliary data are considered in the retrieval:

* **IGBP Landcover**: The IGBP LandCover map Broxton et al. (2014a) is used to describe the surface land classification, and was spatially resampled from its native resolution to the EASE v2.0 25km Grid. The 17 landcovers are gathered into main classes presented in Table 5.1. The VOD retrievals are performed over the two classes NO and FO, for which sets of parameters are defined (see Table below).

* **Soil Grid texture**: The clay content in the soil comes from the Food and Agriculture Organization map FAO (1988). This map was reprojected onto the EASE 2.0 Grid. This information is used in the dielectric model of the soil, to compute the emission of the surface soil.
<table><thead><tr><td>LANDCOVER CODE</td><td>Name</td><td>DT aggreg class</td></tr></thead><tbody><tr><td>0</td><td>Special Discard LC Marker</td><td>None:0</td></tr><tr><td>1</td><td>Evergreen Needle leaf Forest</td><td>CFO:2</td></tr><tr><td>2</td><td>Evergreen Broadleaf Forest</td><td>CFO:2</td></tr><tr><td>3</td><td>Deciduous Needle leaf Forest</td><td>CFO:2</td></tr><tr><td>4</td><td>Deciduous Broadleaf Forest</td><td>CFO:2</td></tr><tr><td>5</td><td>Mixed Forests</td><td>CFO:2</td></tr><tr><td>6</td><td>Closed Shrublands</td><td>CNO:1</td></tr><tr><td>7</td><td>Open Shrublands</td><td>CNO:1</td></tr><tr><td>8</td><td>Woody Savanas</td><td>CFO:2</td></tr><tr><td>9</td><td>Savanas</td><td>CNO:1</td></tr><tr><td>10</td><td>Grasslands</td><td>CNO:1</td></tr><tr><td>11</td><td>Permanent Wetland</td><td>CNO:1</td></tr><tr><td>12</td><td>Croplands</td><td>CNO:1</td></tr><tr><td>13</td><td>urban and Built-Up</td><td>CNO:1</td></tr><tr><td>14</td><td>Cropland/Natural Vegetation Mosaic</td><td>CNO:1</td></tr><tr><td>15</td><td>Snow and Ice</td><td>CEI:7</td></tr><tr><td>16</td><td>Barren or Sparsely Vegetated</td><td>CNO:1</td></tr><tr><td>17</td><td>Ocean water (Saline)</td><td>CWS:5</td></tr><tr><td>18</td><td>Open Water (Pure)</td><td>CWP:4</td></tr></tbody></table>

Figure 5.1: Correspondence IBGP land cover classes and SMOS Fractions as NO: Nominal, FO: Forest, EI: Ice/permanent SNOW, WS: Saline Water, WP: pure water. "C" stands for Classes and "F" for Fraction. From CESBIO technical note SO-TN-CB-GS-057

## 5.3 VOD Output Data file

The last products to be detailed are the output products containing the VOD. Every band is reported as separate files, as well as the orbits which are distinguished between ascending and descending. Note that the ascending/descending orbits do not correspond for SMOS and AMSR-E/2.

The variables in the output files are summarized in the following table:
<table><caption>Table 5.1: VOD file content</caption><tr><td>lat</td><td>latitude in degrees<br>written as double, vector dim=1x584</td></tr><tr><td>lon</td><td>longitude in degrees<br>written as double, vector dim=1x1388</td></tr><tr><td>Soil_Moisture</td><td>Surface Soil Moisture, in m<sup>3</sup>.m<sup>-3</sup><br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Soil_Moisture_StdError</td><td>Soil moisture range due to modeling uncertainties, in m<sup>3</sup>.m<sup>-3</sup><br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Optical_Thickness_Nad</td><td>or VOD, transmittivity of vegetation layer at Nadir, unitless<br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Optical_Thickness_Nad_StdError</td><td>VOD range due to modeling uncertainties, unitless<br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Quality_Flag</td><td>0: data OK, 1: data not recommended, 2: missing data, unitless<br>written as byte, 2D matrix dim=584x1388</td></tr><tr><td>Scene_Flags</td><td>Scene Flags<br>written as byte, 2D matrix dim=584x1388</td></tr><tr><td>Processing_Flags</td><td>Processing Flags, unitless<br>written as byte, 2D matrix dim=584x1388</td></tr><tr><td>RMSE</td><td>RMSE between SMOS-TB and modeled TB, in Kelvin<br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Soil_Temperature-Level1</td><td>ECMWF Soil Temperature at surface level 1, in Kelvin K.<br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>Days</td><td>Number of Days since 1/1/2000<br>written as integer, 2D matrix dim=584x1388</td></tr><tr><td>UTC_Seconds</td><td>Number of Seconds of the day<br>written as integer, 2D matrix dim=584x1388</td></tr><tr><td>UTC_Microseconds</td><td>Number of Microseconds of the day<br>written as integer, 2D matrix dim=584x1388</td></tr><tr><td>RFI_ratio</td><td>Index to evaluate the impact of RFI, unitless<br>written as double, 2D matrix dim=584x1388</td></tr><tr><td>UTC_Minutes</td><td>Number of minutes since the start of the day<br>written as float, 2D matrix dim=584x1388</td></tr></table>

### RFI_ratio

RFI_ratio (Table 5.1) only concerns SMOS-derived VOD. It is computed by dividing the number of RFI flags per bin (from L1c data) by the number of TB per bin. Users can filter the derived SM/VOD based on this ratio of RFI by applying a user-defined threshold. For example, applying a threshold of 0.2 means that any SM/VOD retrievals associated with a RFI_ratio greater than 0.2 should be discarded from any analysis (i.e. more than 20% of the input brightness temperatures are tagged as contaminated by RFI).

### Scene_Flags

The Scene_Flags is a binary flag encoded on a byte that identifies moderate or strong topography, frozen soil (if soil temperature below 273K) or if a pixel is polluted by more than 10% of Water/Urban/Permanent Ice land cover types. The table 5.2 summarizes the flag signification per bit.

<table><caption>Table 5.2: Bit signification of the Scene_Flags</caption><tr><td>Bit</td><td>8</td><td>7</td><td>6</td><td>5</td><td>4</td><td>3</td><td>2</td><td>1</td></tr><tr><td>Flag</td><td>No Flag</td><td>No Flag</td><td>No Flag</td><td>No Flag</td><td>Frozen Soil</td><td>Polluted Scene</td><td>Strong Topography</td><td>Moderate Topography</td></tr></table>

### Processing_Flag

The Processing_Flag is set to one when the RMSE between the modeled TB and the measured TB is above 12K.
Quality_Flag

The Quality flag is set to 1 when a Scene_Flag is raised or the Processing_Flag is raised.

**5.4 Principle of the retrieval**

In this pixel-based approach, the scene is considered homogeneous in terms of land classification and the properties of the pixels (roughness, vegetation, scattering albedo, etc) are area-weighted averages of the properties of the IGBP land cover classes present in the footprint.

In most cases, on appropriate pixels, a couple (SM, VOD) is fed to a forward model computing the corresponding TB at all incidence angles and for both horizontal and vertical polarizations. The values of the (SM, VOD) couple are iteratively modified until the cost function computed from the sum of the squared weighted differences between the modeled and measured TB reaches a minimum. The values of (SM, VOD) which minimize the cost function are written out to the output products. The general principle of the retrieval is described on Figure 5.2 and the following sections describe some of the processing steps more thoroughly.
<figure><img src="image_15.png" /><figcaption>Figure 5.2: Retrieval principle of the processing platform</figcaption></figure>

### 5.4.1 Cost function minimization

The VOD retrieval is based on the minimization of a cost function (pink box on Figure 5.2) of the following form:

$$ FC(Par) = \sum_{\theta,p=\{H,V\}} \frac{(TB_{sat}(\theta, p) - TB_{model}(\theta, p, Par))^2}{\sigma_{TB_{sat}(\theta, p)}^2} + \sum_{Par=sm,tau_{nad}} \frac{(Par - Par_{prior})^2}{\sigma_{Par_{prior}}^2} \quad (5.1) $$
with $TB_{sat}(\theta, p)$ being the satellite brightness temperatures from passive microwave sensor (SMOS, AMSRE, AMSR2), at polarisation $p$ (V Vertical or H Horizontal), $\theta$ the incidence angles, and $\sigma_{TB_{sat}(\theta,p)}$ the radiometric accuracy. $TB_{model}(\theta, p, Par)$ are modeled TB obtained from a radiative transfer model and derived free parameters $Par$ constrained by their initial value $Par_{prior}$ and uncertainties $\sigma_{Par_{prior}}$ (Kerr et al. (2012a, 2020)).

The derived parameters $Par$ are the ones that minimize FC. In most cases, $Par$ are the surface soil moisture and the VOD - Vegetation Optical Depth. For the OSMOSE project, $TB_{model}(\theta, p, Par)$ are computed using two reference radiative transfer models: the Two-Stream model and the L-MEB (L-band Microwave Emission of the Biosphere).

## 5.4.2 Radiative Transfer Model

### 5.4.2.1 L-MEB

The L-MEB is a zero-order radiative transfer model. The multi-layer approach considers the emission of soil, with a rough surface. The contribution of the vegetation is computed using the Tau-Omega model. All the elements are reported in Wigneron et al. (2017), and for SMOS in Kerr et al. (2012a). The Tau/omega model relies on a fixed omega (scattering albedo) and the free parameter Tau ($\tau$), also called the VOD Vegetation Optical Depth.

### 5.4.2.2 2S

The two-stream (2S) is a first-order radiative transfer model accounting for multiple reflections, and developed by Mätzler (1998) and Wiesmann and Mätzler (1999). It is mentioned in this section as it was considered to perform the retrieval. Tests were done but could not be implemented in the processing chains.

### 5.4.2.3 Soil dielectric constant and emissivity

Various dielectric models exist to compute the dielectric constant of a soil. For the OSMOSE project, all cases consider the Mironov dielectric constant model Mironov et al. (2009), as it was designed and validated for the frequencies of interest for the OSMOSE project Montpetit et al. (2015). Then, the soil emissivity is computed using the fresnel law at the soil-atmosphere interface.

## 5.5 Operational configuration of the retrieval processing chains

This section details the setup defined in the radiative transfer model. The retrieval principle described in the previous part is applied to all frequency bands with some adjustments. The processing chain automatically identifies the sensor and the frequency based on the input TB product. The configuration files vary with the frequency, as some parameters need to be adapted. In order to tune the retrieval platform following each sensor's characteristics (e.g frequencies, multi-observation angles for SMOS or one observation angle for AMSRE and AMSR2), the values of several key parameters were tested out and optimized (see section 5.6).

The processing configuration can be summed up as follows:
<table><caption>Table 5.3: Configuration set up of the VOD retrieval algorithm. RTM: Radiative Transfer Model, *: LV low vegetation, **: FO forest</caption><thead><tr><th>Band</th><th>L</th><th>C</th><th>X</th></tr></thead><tbody><tr><td>RTM</td><td>L-MEB</td><td>L-MEB</td><td>L-MEB</td></tr><tr><td>Q</td><td>0</td><td>0</td><td>0.13</td></tr><tr><td>Hr</td><td>0.1 (LV*); 0.3 (FO**)</td><td>0.15</td><td>0.15</td></tr><tr><td>ω</td><td>0.1 (LV*); 0.06 (FO**)</td><td>0.05</td><td>0.05</td></tr><tr><td>derived parameters</td><td>SM-VOD</td><td>SM-VOD</td><td>SM-VOD</td></tr><tr><td>boundary retrieved SM</td><td>None</td><td>[0-1]</td><td>[0-1]</td></tr><tr><td>boundary retrieved VOD</td><td>None</td><td>[0-2]</td><td>[0-2]</td></tr><tr><td>Effective soil temperature</td><td>ERA5 soil level 1</td><td>ERA5 soil level 1</td><td>ERA5 soil level 1</td></tr><tr><td>Vegetation Temperature</td><td>ERA5 skin t.</td><td>ERA5 skin t.</td><td>ERA5 skin t.</td></tr><tr><td>Dielectric model</td><td>Mironov</td><td>Mironov</td><td>Mironov</td></tr><tr><td>SM initial</td><td>0.2</td><td>0.2</td><td>0.2</td></tr><tr><td>VOD init(*)</td><td>global map</td><td>global map</td><td>global map</td></tr></tbody></table>

Initial values of free parameters are needed to start the first iteration of the minimization. SM is fixed at 0.2 m³ m⁻³. The initial VOD (* in Table 5.3) are computed as an average of VOD from a three-year run.
## 5.6 VOD sensitivity analysis

To produce VOD from Brightness Temperatures (TB), the OSMOSE project capitalized on an existing algorithm for the L-band and adapted it to the C and X bands. The Table below (Table 5.4) summarizes the tested parameter values.

<table><caption>Table 5.4: Tested parameters and associated values for tuning the retrieval platform per frequency band</caption><thead><tr><th>Parameter</th><th colspan="3">Tested value</th></tr><tr><th>Band</th><th>L</th><th>C</th><th>X</th></tr></thead><tbody><tr><td>RTM</td><td>LMEB or 2S</td><td>LMEB or 2S</td><td>LMEB or 2S</td></tr><tr><td>Hr</td><td>[0 0.15 0.3 0.4]</td><td>[0 0.15 0.3 0.4]</td><td>[0 0.15 0.3 0.4]</td></tr><tr><td>RFI filtering</td><td>[0.2 0.8]</td><td>N/A</td><td>N/A</td></tr><tr><td>ω</td><td>ω = f(TH)</td><td>[0.01:0.01:0.15]<br>ω = f(TH)</td><td>[0.04 0.05 0.06 0.08]</td></tr><tr><td>Free parameters</td><td>VOD/SM<br>VOD/ACard</td><td>VOD/SM<br>VOD/ACard<br>VOD only (SM from ERA5)</td><td>VOD/SM<br>VOD/ACard<br>VOD only (SM from ERA5)</td></tr><tr><td>Retrieval limits</td><td>Free<br>or 0≤SM≤1 & 0≤VOD≤2</td><td>Free<br>or 0≤SM≤1 & 0≤VOD≤2</td><td>Free<br>or 0≤SM≤1 & 0≤VOD≤2</td></tr><tr><td>Veg. Temperature</td><td>Skin Temperature<br>2m Air Temperature</td><td>Skin Temperature<br>2m Air Temperature</td><td>Skin Temperature<br>2m Air Temperature</td></tr></tbody></table>

Here are some explanations of some tested parameters:

* **RFI filtering**: to assess the retrieval quality when applying different RFI filtering thresholds, for SMOS only. The L3TB impacted by RFI are filtered using the ratio Nb\_RFI\_Flags / Nviews, i.e. numbers of SMOS TB impacted by RFI over the number of TB. Two thresholds are tested, 0.2 being too restrictive removing too many TBs. A threshold of 0.8 allows us to deal with more TBs. It is justified as two more filters are applied on the VOD: removing VOD with an RMSE~8K and outliers when computing the yearly average.

* **ω**: was already fine-tuned for the L-band MEB. Tests were carried out with a Tree Height (TH) dependent omega based on the results from WSL analyses (Zhou et al. (2024)). For C and X bands, additional tests with fixed ω values were conducted.

* **Free parameters**: As it was unsure whether SM would be a relevant parameter for C and X-bands, some tests were conducted to retrieve VOD and ACard, instead of VOD and SM. ACARD is the cardioid model Waldteufel et al. (2004), which is a simple model to estimate the module of the dielectric constant. It can be used to characterize particular cases (frozen/bare soils) where the soil model (i.e. Mironov dielectric model) is not adapted. In the cases of X and C bands, we tested 2-parameter retrievals with ACard instead of the SM, as it was not evident that performing a SM/VOD retrieval would lead to satisfying results in terms of VOD and also in terms of RMSE (modeled TB model versus SMOS TB). Indeed, C and X is less sensitive to the soil in particular with a vegetation layer. Additionally, the retrieval of a single parameter (with SM from ERA5) versus the retrieval of two parameters was also tested.

* **Retrieval limits**: Assess the impact on the VOD, of constraining the free parameters in the cost function.

In the following section, we present the sensitivity of the VOD at X and C bands with respect to the various parameters, as described above.

### 5.6.1 Influence of the Radiative Transfer Model

In the frame of this project, two Radiative Transfer Models (RTM) were tested:
* The L-band Microwave Emission of the Biosphere Model (LMEB) which uses the classical τ – ω model
  *Wigneron et al. (2017)*

* The Two-Stream Model (2S) which is able to simulate higher degrees of interaction between different layers (soil, vegetation, atmosphere) *Schwank and Naderpour* (2018)

A test for the L-band was conducted over 1 month for both ascending and descending orbits and using both RTM. The ω value is converted from LMEB values to 2S values using Eq. 23 from *Schwank and Naderpour* (2018). It means that the 2S model is not used per se, but is rather fitted into the LMEB. The ω<sub>LMEB</sub> values are 0.06 over forest pixels and 0.1 over low vegetation pixels. The output VOD were then averaged for each orbit type. One of the final purposes of the microwave part of the OSMOSE project is the estimation of the AGB from multi-frequency VOD. When conducting the sensitivity analysis of different parameters on the VOD, the spatial correlation of this averaged VOD with respect to the CCI AGB is a relevant indicator. In this study, the Pearson correlation coefficient (R) and Spearman correlation coefficient (Rho) are used. Figure 5.3 shows the spatial density plot between the AGB from the CCI *Santoro and Cartus* (2023) and the 1-month average VOD for ascending and descending orbits using the 2S or the LMEB for the descending (top panel) or ascending orbits (bottom panel).
Comparison of the L-VOD distribution generated with 2S
(blue) and TO (red) RTM wrt AGB from CCI 2018.

<figure><img src="image_16.png" /></figure>

Figure 5.3: Density plot of the OSMOSE L-VOD for June 2011 against the AGB from CCI 2018. The red color refers to the use of the TO model to retrieve the VOD and the blue color refers to the use of the 2S model. In both cases, the couple (VOD/ACard) was retrieved, the retrieval was not constrained and the default values of ω were used (eq 23 of *Schwank and Naderpour* (2018) was used to convert from TO to 2S value). The shaded areas correspond to the interval in-between quantiles 16 and 84 and contains 68% of the distribution. (a) for descending orbits and (b) for ascending orbits.

Figure 5.3 shows that the RTM does not significantly impact the VOD correlation with the AGB for the L-band for both orbit types. In both cases, the 2S L-VOD saturates to a lower value than the TO L-VOD. The choice of the RTM does not greatly impact the distributions.

A similar test was conducted for C-band over July 2011. The ω value are derived from the MEMLV (outputs of OSMOSE WP1, WSL) and depends on the the TH. A map of the C-band omega for the TO RTM is displayed in Figure 5.4. The ω values from TO were converted to 2S using Eq. (23) of *Schwank and Naderpour* (2018). The results are displayed on Figure 5.5.
<figure><img src="image_17.png" /><figcaption>Figure 5.4: Map of ω for the C-band TO RTM based on the WSL Look-Up Tables</figcaption></figure>
<figure><img src="image_18.png" /></figure>

Figure 5.5: Density plot of the OSMOSE C-VOD for July 2011 against the AGB from CCI 2018. The red color refers to the use of the TO model to retrieve the VOD and the blue color refers to the use of the 2S model. In both cases, the couple (VOD/ACard) was retrieved, the retrieval was not constrained and the default values of ω were used (eq 23 of *Schwank and Naderpour* (2018) was used to convert from TO to 2S value). The shaded areas correspond to the interval in-between quantiles 16 and 84 and contain 68% of the distribution.

With the 2S configuration, a significantly large area over Eastern Russia exhibits unexpected high VOD. This area is highlighted in purple on panel (a) of Figure 5.5 and the corresponding points are also highlighted in purple on panel (b). These 2S VOD are higher than the VOD over the Tropics and do not seem realistic as they correspond to AGB values around 50-100 Mg ha⁻¹. The correlation also drops with the 2S RTM. Hence, the TO model seems to be better suited for the C-band. The same results apply to the X-band.

## 5.6.2 Influence of the retrieval parameters

The retrieval algorithm allows us to derive various parameters that drive the RTM, and the emission observed by the sensors. The most obvious ones are the SM or the dielectric constant for the soil layer, and the VOD for the vegetation layer, as described at the beginning of the section.

### 5.6.2.1 L-Band

For SMOS L-band, two couples of free parameters were tested:
* Retrieving (VOD, SM)

* Retrieving (VOD, ACard) and setting SM to the default initial value of 0.2 m⁻³/m⁻³

Comparison of the L-VOD distribution generated with TO RTM retrieving (VOD/ACard) (red) or (VOD/SM) (blue)
Asc & Desc orbits are mixed

<figure><img src="image_19.png" /></figure>

Figure 5.6: Density plot of the OSMOSE L-VOD for June 2011 against the AGB from CCI 2018. The red color refers to the retrieval of the (VOD/ACard) couple and the blue color to the (VOD/SM) couple. The TO retrieval was not constrained and the default values of ω were used. The shaded areas correspond to the interval in-between quantiles 16 and 84 and contains 68% of the distribution.

Figure 5.6 shows that retrieving (VOD/ACard) or (VOD/SM) does not significantly impact the VOD correlation with the AGB from CCI 2018. It is not surprising as the VOD are monthly averaged, which probably smooths the small differences that may occur. The derived Acard is the dielectric constant, which the SM mainly drives at L-band. SM and Acard are in that case linked and lead to similar VOD. All correlation coefficients are around 0.83. The dispersion, highlighted by the colored shaded area, is lower when using the classical (VOD/SM) couple. The (VOD/SM) couple is then considered in the OSMOSE processing chain for L-band.

### 5.6.2.2 C-Band

The same couples (VOD/SM or VOD/ACard) were tested for the C-band. As the observing frequency is higher than L-band, there was a possibility that the C-band observations were not sensitive to the ground underneath a vegetation layer, and that retrieving ACard would be more relevant than retrieving SM. Like for L-band, retrieving Acard rather than SM does not impact the correlation of the VOD with the AGB from CCI 2018. (VOD/SM) will hence be retrieved for the C-band. This also ensures homogeneity with L-band retrieval.

### 5.6.3 Influence of constraining the retrieval with upper and lower boundaries

Constraining the retrieval does not affect the L-VOD correlation to the AGB from CCI 2018. The processing chain is unconstrained for the L-band.

For the C-band, constraining the retrieval slightly improves the VOD Spearman correlation with the AGB (Pearson remains about the same). The dispersion of the C-VOD wrt AGB marginally decreases for VOD
values lower than 0.7 when setting lower and upper boundaries to the free parameters. Also, constraining the
retrieval greatly improves the convergence between the demo chain (coded in matlab and hosted at CESBIO)
and the operational chain (coded in Python and hosted at Magellium). In order to ensure results as close as
possible to the demo trials, the retrievals are constrained between [0-2] for the VOD and [0-1] for SM for the
C and X bands.

5.6.4 Influence of the ω and *Hr* parametrization

Another important aspect of the processing chain is the parametrization of the scattering albedo (ω and the roughness of the soil, *Hr*). This test was only conducted for the C and X bands as the ω and *Hr* parametrization for L-band had already been extensively studied and optimized, as described in the literature (Fernandez-Moran et al. (2017) and Wigneron et al. (2017)).

For this test, the OSMOSE retrieval platform was run with fixed values for ω and *Hr*. For the first runs, the ω values spanned [0.01-0.15] and the *Hr* values spanned [0 0.15 0.3 0.4]. The AMSRE C and X band VOD were derived for the month of July 2011 for ascending and descending orbits. Then, the temporal and spatial mean RMSE between the modeled (*TB*<sub>mod</sub>) and the measured *TB* (*TB*<sub>meas</sub>) was computed. This single mean RMSE value allowed us to compare the performance of the different (ω/*Hr*) configurations. The results of these tests for C-band over the global region are displayed in Figure 5.7. This figure shows that *Hr* has little influence on the RMSE compared to the value of ω. The RMSE slightly decreases from ω=0.01 to ω=0.05. The RMSE then increases with increasing values of ω. The configuration ω=0.05 and *Hr*=0.15 is selected for the C-band archive. The test was also conducted for the month of September 2011 and led to the same conclusion.

<figure><img src="image_20.png" /><figcaption>Figure 5.7: Mean RMSE between C-band *TB*<sub>mod</sub>*TB*<sub>meas</sub> over the global region for different *ω* and *Hr* couples and a 2 parameters (VOD-SM) retrieval. The RMSE was averaged over the course of July 2011, mixing ascending and descending orbits</figcaption></figure>

The same test for the X-band leads to the identical ω and *Hr* configuration as for the C-band.
### 5.6.5 Influence of retrieving one or two parameters

As AMSRE and AMSR2 only have one incidence angle, only two TB are available per observation, one for horizontal and one for vertical polarizations. Moreover, these two TB are not completely independent. Retrieving two parameters with two correlated TB may lead to erroneous results. This is why tests with a one-parameter retrieval were conducted with the same configuration described in the previous section. The results of these tests for the global region and the C-band are displayed in Figure 5.8. When comparing this figure with Figure 5.7, it is clear that a one-parameter retrieval leads to higher RMSE between the C-band $TB_{mod}$ and $TB_{meas}$. The same conclusion applies to the X-band. Hence, a two-parameter retrieval is adopted for the C and X-band archives. As a side note, Figure 5.8 is concordant with Figure 5.7 and shows the little influence of the value of $Hr$ compared to the value of $\omega$. Again, the lowest global RMSE is achieved with $\omega=0.05$.

<figure><img src="image_21.png" /><figcaption>Figure 5.8: Mean RMSE between C-band $TB_{mod}$ and $TB_{meas}$ over the global region for different $\omega$ and $Hr$ couples and a 1 parameter (VOD) retrieval. The RMSE was averaged over the course of July 2011, mixing ascending and descending orbits</figcaption></figure>

## 5.7 Validation of VOD

### 5.7.1 versus datacube2

Datacube 2 is a collection of datasets from visible sensors, projected on the ESA Grid version 2 (Brodzik et al. (2012a, 2014)). Refer to chapter 8 for a complete description of the datacube2.
5.7.1.1 L-Band

OSMOSE L-VOD vs CCI AGB

<figure><img src="image_22.png" /><figcaption>Figure 5.9: Density plot of the yearly L-VOD from OSMOSE wrt the CCI AGB for the years 2017-2020</figcaption></figure>

Figure 5.9 and Table 5.5 show that the OSMOSE L-VOD is very well correlated to the CCI AGB for all years. The relation is stable from one year to the other.

<table><caption>Table 5.5: R between the L-VOD from OSMOSE and AGB values from CCI for the years 2017-2020.</caption><thead><tr><th></th><th>2017</th><th>2018*</th><th>2019</th><th>2020</th></tr></thead><tbody><tr><th scope="row">R</th><td>0.86</td><td>0.86</td><td>0.86</td><td>0.86</td></tr><tr><th scope="row">Rho</th><td>0.9</td><td>0.9</td><td>0.9</td><td>0.9</td></tr></tbody></table>
OSMOSE L-VOD vs LAI Copernicus

<figure><img src="image_23.png" /><figcaption>Figure 5.10: Density plot of the yearly L-VOD from OSMOSE wrt the LAI from Copernicus for the years 2010, 2013, 2016 and 2020</figcaption></figure>

Figure 5.10 shows that the OSMOSE L-VOD is very well correlated to the LAI for all years. Only 4 years are displayed in Figure 5.10, but the relationship is stable from one year to the other.

Even though they are not displayed in this technical note, regressions between the L-VOD and other optical indices have been conducted:

* FAPAR (2010-2016): $0.82 \le R \le 0.86$ and $0.82 \le rho \le 0.86$

* FCOVER (2010-2020): $0.82 \le R \le 0.86$ and $0.81 \le rho \le 0.85$

* NDVI (2010-2020): $0.49 \le R \le 0.66$ and $0.51 \le rho \le 0.7$
5.7.1.2 C-Band

OSMOSE C-VOD vs CCI AGB

<figure><img src="image_24.png" /></figure>

Figure 5.11: Density plot of the yearly C-VOD from OSMOSE wrt to the AGB from CCI for the years 2017-2020

Figure 5.11 shows that the OSMOSE C-VOD is well correlated to the AGB from the CCI for all years and the relation is stable from one year to the other. In agreement with the literature (Rodríguez-Fernández et al. (2018)), the general correlation with the AGB is not as good as with the L-VOD (see Figure 5.9).
OSMOSE C-VOD vs LAI Copernicus

<figure><img src="image_25.png" /><figcaption>Figure 5.12: Density plot of the yearly C-VOD from OSMOSE wrt to the LAI from Copernicus for the years 2010, 2013, 2016 and 2020</figcaption></figure>

Figure 5.12 shows that the OSMOSE C-VOD is very well correlated to the LAI for all years. Only 4 years are displayed in Figure 5.12 but the relationship is stable from one year to the other. The linear correlation R is not as good as the one with the LVOD but the spearman correlation coefficients are very close (see Figure 5.10).

Even though they are not displayed in this technical note, regressions between the C-VOD and other optical indices have been conducted:

* FAPAR (2002-2016): $0.84 \le R \le 0.88$ and $0.84 \le rho \le 0.87$

* FCOVER (2002-2020): $0.83 \le R \le 0.86$ and $0.83 \le rho \le 0.86$

* NDVI (2002-2020): $0.6 \le R \le 0.76$ and $0.62 \le rho \le 0.77$
5.7.1.3 X-Band

OSMOSE X-VOD vs CCI AGB

<figure><img src="image_26.png" /><figcaption>Figure 5.13: Density plot of the yearly X-VOD from OSMOSE wrt the AGB from CCI for the years 2017-2020</figcaption></figure>

Figure 5.13 shows that the OSMOSE X-VOD is well correlated to the AGB from the CCI for all years and the relation is stable from one year to the other. In agreement with the literature (Rodríguez-Fernández et al. (2018)), the general correlation with the AGB is not as good as with the L-VOD (see Figure 5.9) but is very close to the correlation of the C-VOD with AGB (see Figure 5.11).
OSMOSE X-VOD vs LAI Copernicus

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>Figure 5.14: Density plot of the yearly X-VOD from OSMOSE wrt the LAI from Copernicus for the years 2010, 2013, 2016 and 2020</figcaption></figure>

Figure 5.14 shows that the OSMOSE X-VOD is very well correlated to the LAI for all years. Only 4 years are displayed in Figure 5.14 but the relationship is stable from one year to the other. The linear correlation R is not as good as the one with the LVOD but the spearman correlation coefficients are very higher (see Figure 5.10). In general, the correlations (R and Rho) are better than the C-VOD ones (see Figure 5.12).

Even though they are not displayed in this technical note, regressions between the X-VOD and other optical indices have been conducted:

* FAPAR (2002-2016): $0.86 \le R \le 0.90$ and $0.86 \le rho \le 0.89$

* FCOVER (2002-2020): $0.87 \le R \le 0.89$ and $0.85 \le rho \le 0.89$

* NDVI (2002-2020): $0.66 \le R \le 0.79$ and $0.68 \le rho \le 0.8$

5.7.2 versus VODCA

The version 2 (v2) of the VODCA dataset Zotta et al. (2024) was just released. In this new version of the dataset, the C, X and Ku bands were merged into one daily data archive. It includes the use of AMSRE and AMSR2 descending orbits processed in the LPRM. The L-band created from SMOS and SMAP ascending orbits processed in the LPRM algorithm, is distributed as a 10-day composite data archive. The full data archive is provided on a 0.25° grid, close to the EASE 2 grid used in the OSMOSE project. The VODCA v2 data archive was reprojected to the OSMOSE grid.

The yearly averaged OSMOSE VOD were compared to the yearly averaged VODCA dataset.
### 5.7.2.1 L-band

Figure 5.15 shows that the correlation between OSMOSE and VODCA for the L-band is excellent and uniform over the years. The R (Pearson coefficient) and Rho (Spearman coefficient) correlation coefficients are over 0.89 for all years between 2010 and 2020.

<figure><img src="image_29.png" /><figcaption>Figure 5.15: Density plot of the yearly L-VOD from OSMOSE wrt to VODCA v2 for (a) 2010, (b) 2013, (c) 2016 and (d) 2020</figcaption></figure>

Generally, VODCA VOD is lower than OSMOSE VOD, except for a region where the two have more similar values (i.e. 0.2-0.5, above the linear fit black line in Figure 5.15). These points are mainly located over boreal forests, and further investigations are needed to understand better these differences. The VODCA VOD also saturates at a lower value (~0.8) compared to the OSMOSE VOD which saturates at around 1.2 being more able to capture the vegetation dynamic.

### 5.7.2.2 C-band

Figure 5.16 shows that the correlation between OSMOSE and VODCA for the C-band is also excellent and uniform over the years. The ratio between OSMOSE C-band and VODCA CXKu is close to 1. The R and Rho correlation coefficients are over 0.87 for all years between 2003 and 2020.
<figure><img src="image_30.png" /><figcaption>Figure 5.16: Density plot of the yearly C-VOD from OSMOSE wrt to VODCA CXKu v2 for (a) 2010, (b) 2013, (c) 2016 and (d) 2020</figcaption></figure>

5.7.2.3 X-band

Figure 5.17 shows that the correlation between OSMOSE and VODCA for the X-band is also excellent and uniform over the years. The ratio between OSMOSE X-band and VODCA CXKu is close to 1. The R and Rho correlation coefficients are over 0.88 for all years between 2003 and 2020.
<figure><img src="image_31.png" /><figcaption>Figure 5.17: Density plot of the yearly X-VOD from OSMOSE wrt to VODCA CXKu v2 for (a) 2010, (b) 2013, (c) 2016 and (d) 2020</figcaption></figure>

A similar evaluation is conducted per land cover classes. Table 5.6 gathers all the correlation coefficients per IGBP class for OSMOSE L and C VOD and VODCA L and CXKu VOD.

<table><caption>Table 5.6: Correlation Coefficient (R) between the CCI AGB (year 2018) and VODs from OSMOSE (L/C/X) and VODCA.</caption><thead><tr><th></th><th>OSMOSE<br>L-VOD</th><th>VODCA<br>L-VOD</th><th>OSMOSE<br>C-VOD</th><th>OSMOSE<br>X-VOD</th><th>VODCA<br>CXKU-VOD</th></tr></thead><tbody><tr><td>Evergreen needleleaf forests</td><td><strong>0.58</strong></td><td>0.42</td><td>0.45</td><td>0.56</td><td>0.50</td></tr><tr><td>Evergreen broadleaf forests</td><td><strong>0.71</strong></td><td>0.66</td><td>0.55</td><td>0.55</td><td>0.52</td></tr><tr><td>Deciduous needleleaf forests</td><td>0.43</td><td><strong>0.53</strong></td><td>0.3</td><td>0.44</td><td>0.34</td></tr><tr><td>Deciduous broadleaf forests</td><td><strong>0.44</strong></td><td>0.3</td><td>-0.18</td><td>-0.04</td><td>0.29</td></tr><tr><td>Mixed forests</td><td><strong>0.49</strong></td><td>0.15</td><td>0.18</td><td>0.36</td><td>0.16</td></tr><tr><td>Closed shrublands</td><td>0.52</td><td><strong>0.68</strong></td><td>0.08</td><td>0.42</td><td>0.66</td></tr><tr><td>Open shrublands</td><td>0.56</td><td><strong>0.59</strong></td><td>0.51</td><td>0.49</td><td>0.54</td></tr><tr><td>Woody savannas</td><td><strong>0.54</strong></td><td>0.53</td><td>0.38</td><td>0.42</td><td>0.44</td></tr><tr><td>Savannas</td><td><strong>0.74</strong></td><td>0.72</td><td>0.52</td><td>0.65</td><td>0.62</td></tr><tr><td>Grasslands</td><td><strong>0.54</strong></td><td>0.52</td><td>0.44</td><td>0.48</td><td>0.47</td></tr><tr><td>Permanent wetlands</td><td>0.43</td><td><strong>0.52</strong></td><td>0.32</td><td>0.39</td><td>0.51</td></tr><tr><td>Croplands</td><td>0.48</td><td><strong>0.50</strong></td><td>0.39</td><td>0.37</td><td>0.50</td></tr><tr><td>Urban</td><td>0.50</td><td><strong>0.54</strong></td><td>0.29</td><td>0.16</td><td>0.38</td></tr><tr><td>Cropland and vegetation mix</td><td>0.60</td><td><strong>0.61</strong></td><td>0.4</td><td>0.44</td><td>0.57</td></tr><tr><td>Snow Ice</td><td><strong>0.67</strong></td><td>0.20</td><td>0.34</td><td>0.41</td><td>0.04</td></tr><tr><td>Barren</td><td><strong>0.16</strong></td><td>0.14</td><td>0.07</td><td>0.09</td><td>0.12</td></tr></tbody></table>
### 5.7.3 Versus Wang X-VOD

The OSMOSE X-VOD produced with the TO model with $\omega=0.05$ and $Hr=0.15$ using a constrained 2 parameter retrieval has been compared with the estimation available from Wang et al. (2024). In this paper, the authors conducted a regression between the 2017 X-VOD and the AGC from CCI 2017 over Africa. Figure 5.18 shows that the OSMOSE X-VOD averaged in 2017 over Africa shows a similar relationship with AGB from CCI as what Wang et al. (2024) found.

<figure><img src="image_32.png" /><figcaption>Figure 5.18: (a) Density plot of the 2017 X-VOD from Fig1 of Wang et al. (2024) wrt AGC CCI over Africa. (b) Density plot of the July 2011 OSMOSE X-VOD wrt to AGB CCI over Africa.</figcaption></figure>

## 5.8 AMSRE/AMSR2

AMSR2 is the follow-up mission of AMSRE and there is a nine-month gap between both sensors. One of the objectives of the OSMOSE project being to build a 20-year time series of C and X VOD, the continuity of the AMSRE/AMSR2 data needs to be assessed and potentially corrected. Several studies have already dealt with AMSRE/2 continuity and put up different methodologies. Wang et al. (2021b) established a linear relationship between the AMSRE and AMSR2 TB (in H and V polarization respectively) over stable and homogeneous land covers (Evergreen Broadleaf forest and snow and ice). This relationship was then applied globally to calibrate the AMSRE TB to the AMSR2 ones. Finally, the calibrated TB were in turn processed by the X-MEB (Wang et al. (2021b)) to produce a consistent 20-year-long X-VOD time series. The authors of Moesinger et al. (2020a) and Zotta et al. (2024) worked directly on the VOD and employed CDF matching techniques to scale AMSR2 to AMSRE. The following sentences are directly extracted from Zotta et al. (2024): VODCA v1 scaled AMSR2 to TMI if enough overlap was available or directly to AMSR-E, without temporal overlap, above and below 35° latitude N and S, respectively Moesinger et al. (2020a). This led to spatial inconsistencies in VODCA v1 [...]. Therefore, we changed the approach in VODCA v2 and used SSMI F17 Ku-band observations (scaled to AMSR-E X-band) as reference to bridge the gap between AMSR-E and AMSR.

For the OSMOSE project, the bias and scaling are studied and performed directly on the C and X-VOD of AMSRE and AMSR2 and the datacube #2 is used to check the consistency of the correction.

### 5.8.1 Assessment of the AMSRE/AMSR2 bias for the C-band

Since there is no overlap between AMSRE and AMSR2, the continuity of the VOD between both sensors was studied by taking into account the last year and a half of AMSRE (1st Jan 2010 to 27th Sep 2011) and the first year and a half of AMSR2 (3rd Jul. 2012 to 31st Dec. 2014), following the method described in Moesinger et al. (2020a). The mean of the C-VOD was computed for both sensors over their respective period. The difference between the mean VOD is then computed. This difference is referred to as a bias even though it does not purely represent a bias as there is no overlap between the sensors.
To better understand the bias between AMSRE and AMSR2 C-VOD, it is compared to the mean difference of various optical indices from the datacube #2 (FAPAR, FCOVER, LAI, NDVI) and also to the mean difference of the OSMOSE L-VOD over the same period. The global spatial correlation between the C-VOD bias and the optical indexes/L-VOD is also computed. Results for FAPAR appear in Figure 5.19.

<figure><img src="image_33.png" /><figcaption>Figure 5.19: (a) Mean difference between the AMSRE C-VOD (01/01/2010-09/27/2011) and the AMSR2 C-VOD (07/03/2012-12/31/2013). (c) Mean difference of the FAPAR between 2012/2013 and 2010/2011. (d) Global spatial correlation between FAPAR and C-VOD bias.</figcaption></figure>

Globally, the C-VOD bias is not correlated to any of the optical indices or L-VOD bias. *R* is lower than 0.41 for all cases and *Rho* is lower than 0.34. This indicates that globally, the C-VOD bias is not natural. Interestingly, the C-VOD bias is much more dynamic than the L-VOD bias. However, looking up Figure 5.19 and other indices show the same trend, 3 areas have similar C-VOD biases and optical/L-VOD biases: the eastern tip of Brazil, northern South Africa, and eastern Australia. For these three regions, the observed C-VOD bias is also observed in the optical indices/L-VOD bias as well, meaning that the C-VOD bias there is (at least) partly natural.

## 5.8.2 Correction of the AMSRE/AMSR2 bias

These comparisons highlight that, globally, the bias between AMSRE and AMSR2 needs to be corrected pixel-wise. One well-established method to correct a bias is the Cumulative Distribution Function (CDF) matching as described and applied in section III.B. of Madelon et al. (2022). The following sentences are directly extracted from Madelon et al. (2022): The CDF is a specific way to describe the distribution of a discrete or continuous variable X. The purpose of this function is to give the probability that X will take a value less than or equal to a certain threshold. [...] CDF matching consists of transforming the CDF of one variable X (source data) to mimic that of another variable Y (reference data) by using a function *f*.

For this study, AMSRE is rescaled to AMSR2 as AMSR2 is the most recent sensor. As a consequence, the source data are the last two years of AMSRE (28th Sept. 2009 to 27th Sept. 2011) and the reference data are the first two years of AMSR2 (3rd Jul. 2012 to 2nd Jul 2014). Different CDF matching functions can be used. The most simple one is a linear function using the mean and the standard deviation of the source and reference data. To use the linear CDF matching, the data time series need to be Gaussian. Several statistical tests (Chi-square goodness-of-fit test, the Lilliefors test *Lilliefors* (1967) and the Jarque-Bera test
Jarque and Bera (1987)) were applied on the AMSRE and AMSR2 C-VOD time series. These stastistical tests showed that the 2-year long time series could not be considered as Gaussian. Hence, the Polynomial fitting CDF matching from Brocca et al. (2011) was used to rescale AMSRE to AMSR2. Madelon et al. (2022) summarizes the principle: It uses a function approximation technique to find the relationships between the source and the reference data CDF. First, the source and reference CDFs are computed with the number of percentiles equal to the number of SM observations in their respective time series. Second, the source CDF is linearly interpolated on the percentiles of the reference CDF, and the difference between the two is computed. Then, the difference is plotted against the interpolated source CDF, and a polynomial fit is performed. Finally, the polynomial function is used to compute the correction to apply to the source data. As mentioned in Aires et al. (2021), this method has the advantage to provide a continuous transformation, but it is more complex to control. Indeed, it can be subjected to instabilities and extreme events are not sure to be well preserved.

A polynomial of degree three was used to rescale the C-VOD time series of each pixel. If there are fewer than 50 points in either the source or reference time series, the pixel is dropped (4 points in total).

### 5.8.3 Validation of the CDF matching

Time series of VOD with/without CDF matching correction are checked against SMOS L-VOD and optical indices from OSMOSE data cube 2. An example is displayed below:

<figure><img src="image_34.png" /></figure>

Figure 5.20: For the point [-11.41° N, -39.77°E] (eastern tip of Brazil): 2003-2014 time series of the AMSRE/AMSR2 C-VOD before CDF matching (red), 2003-2011 time series of AMSRE C-VOD after CDF matching with AMSR2 (green), 2010-2014 time series of the SMOS L-VOD (yellow), 2003-2014 time series of FAPAR from Copernicus (cyan). The L/C-VOD time series are 10 day moving averages.
<figure><img src="image_35.png" /><figcaption>Figure 5.21: Yearly VOD Bias for AMSRE after minus before CDF matching for C (top) and X (bottom) band</figcaption></figure>

Correction of the VOD from the AMSR sensors is mandatory. Figure 5.21 shows the significant differences when applying the CDF matching correction. For our analysis, we chose to correct the AMSRE VOD, as there is no absolute method to clearly state which sensor (i.e. AMSRE or AMSR2) provides the reference TB. The correction could also have been done at the TB level as in Wang et al. (2021a). It constitutes the most appropriate method. However, we chose to keep the "raw" VOD, so that we could perform an analysis with optical data from the datacube2. It seems that on a few nodes, the correction is not justified. One has to keep in mind that it may be the case for other grid nodes. Both VODs (i.e. not corrected and corrected with CDF matching) are then produced and delivered for the OSMOSE project, as one can further investigate the need to apply a correction, whether on AMSRE or AMSR2 sensors.
# Chapter 6

# Analysis of AGB (WP5)

## 6.1 Introduction

The purpose of this section is to detail the setup of the AGB estimation. Data from the SMOS, AMSRE and AMSR2 sensors are used to derive the AGB. The following represents the activities developed for the work package 5 of the project.

## 6.2 Data

### 6.2.1 L-Band Level 3 Brightness Temperature and VOD/SM data-sets

Level 3 TBs (L3TB) are provided in fourteen angle bins of 5º width (excluding bin number nine). Hereafter bin No.1 denotes incidence angles from 0 to 5º, bin No.2 from 5º to 10º and so on until bin No.14; exceptionally, the complementary bin number 9 is centered at 40º with a 2º width. L3TB TBs are provided in horizontal and vertical polarizations, hereafter H and V polarizations, respectively.

The Level 3 L-VOD/SM product (L3SM version 339), hereafter VOD (and SM), and TBs (version 330) are produced through the implementation of the L-MEB (L-Band Microwave Emission of the Biosphere) radiative transfert model *Wigneron et al. (2007a)* on the SMOS Level 3 retrieval algorithm *Al Bitar et al. (2017a)*.

Both the L3SM and L3TB data sets are provided by the CATDS (Centre Aval de Traitement de Données SMOS) from 2010 until 2023 on a daily time step, and are projected on the EASE (Equal-Area Scalable Earth) grid version 2. For this study, only the observations corresponding to the years of the AGB reference maps (2010, 2017 and 2018) were used. TBs from the ascending orbits were used. To remove the dependence to seasonal variations of VWC on the signal, annual averages of each variable were calculated for the aforementioned years.

A daily filtering was performed on the TB and VOD/SM products before calculating the annual averages. The probability of appearance of Radio Frequency Interference (RFI) was used as a quality criterion. In the case of TBs, any observation whose probability of appearance of RFIs (number of RFIs detected/number of daily observations) by angle of incidence was greater than 0.7 was discarded. For VOD products, in addition to the RFI probability filter (>0.2), a $\chi^2$ coefficient (which measures the goodness of fit from TB estimates and TB observations during the minimization process of the VOD/SM) lower than 2 was selected.

Additionally, and taking into account that the SMOS retrievals over coastal pixels tend to be strongly affected by the fraction of water in the observed scene (especially at high incidence angles), all observations one pixel away from the coasts (defined by the USGS land sea mask used by the SMOS algorithm) were removed.
### 6.2.2 C- and X-band VOD

Producing AGB with a multi-frequency approach requires a harmonized multi-frequency VOD dataset. The frequencies considered in the OSMOSE project are the L-band from SMOS and the C and X-bands from the Advanced Microwave Scanning Radiometer (AMSRE) and its follow-up mission AMSR2. The multi-frequency VOD dataset results from the processing of L, C and X-bands TB through a RTM (the so called τ-ω model) which parametrization varies from one frequency to the other. For now, L and C-VOD have been produced in the frame of the OSMOSE project. The L-VOD dataset spans 2010-2020 and the C and X-VOD product span 2002-2020. The products and the retrieval chain are described in the Chapter 5.

### 6.2.3 Above Ground Biomass reference maps

Three AGB maps from ESA's CCI biomass project with global coverage (Santoro et al. (2021)) were used as vegetation biomass references for the years 2010, 2017 and 2018 (datasets version 3.0), from now on referred to as CCI2010, CCI2017 and CCI2018, respectively. Reference maps are provided with a spatial resolution of 100 m, and representing biomass density in Mg/ha. The CCI2010 is produced from ALOS-1 (PALSAR-1) and Envisat ASAR data-sets, the CCI2017-18 are produced from ALOS-2 (PALSAR-2) and Sentinel-1 observations of SAR backscatter. According to the validation protocol applied by the CCI biomass project (comparing the CCI maps with true AGB observations from forest field inventories), the CCI2017 and CCI2018 were in better agreement with AGB references than the one from 2010. Though the CCI2010 showed a lower agreement with field estimations, it still offers an approximation of the biomass for that reference year, and therefore it was used in this study. In spite of their possible uncertainties, the CCI AGB maps will be referred to as 'references', assuming that they are a good approximation to the actual state of vegetation biomass globally.

In addition to the maps from the CCI project, we also tested as reference the AGB map from Avitabile et al. (2016). The latter developed a 1-km resolution AGB map, specifically for the pantropical region, by merging two earlier AGB maps created by Saatchi et al. (2011) and Baccini et al. (2012), targeting the year 2010. This combination utilized a novel collection of field data and adjusted high-resolution AGB maps. These elements were standardized and combined, yielding nearly 15,000 reference AGB data points at a 1-km resolution. To ensure accuracy, a data merging method that included bias correction and weighted linear averaging was applied across specific zones.

For this study, all AGB reference maps have been re-gridded to the common grid of SMOS L3 products (EASE grid v2) by averaging the original 100 m pixels within the ~25 km EASE grid cells.

### 6.2.4 Data for AGB evaluation

We used for AGB validation lidar-based AGB data, in situ AGB data and FAO data at the national scale. We also compared our maps with AGB maps from Xu et al. (2021).

#### 6.2.4.1 Lidar-based AGB data

Four Lidar-based AGB dataset were used in this study for evaluating our AGB maps. The characteristics of these AGB maps are summarized in Table 6.1. They have been processed in the frame of the ESA CCI Biomass project by Nicolas Labrière from the Évolution et Diversité Biologique Laboratory, University Toulouse 3 Paul Sabatier, Toulouse, France. The four dataset were acquired in Australia, USA, Brasil and Indonesia, and include 499 lidar-based AGB maps at 100m pixel size, representing in total an area of 24 340 km², acquired from 2011 to 2019. This timeframe is particularly suited to our analysis, in which we provided maps from 2011 to 2022.

The work of Nicolas Labrière was similar for the Australian TERN (Terrestrial Ecosystem Research Network), the National Ecological Observatory Network (NEON), the Sustainable Landscape Brazil (SLB) and the Carbon Monitoring System (CMS) consisted in the processing of tree inventory and airborne LiDAR data. After data cleaning, tree inventory data included approximately:
* TERN: 28 800 living trees from 9 inventories at 5 different sites scattered across Australia. They were converted into plot- and subplot-level AGB estimates and associated uncertainty, for a total of 188 units (34, 1 and 153 of 1ha, 0.5ha and 0.25ha, respectively),

* NEON: 58 800 living trees from 126 inventories at 37 different sites scattered across the conterminous United States, Alaska, Hawaii and Puerto Rico. They were converted into plot-level AGB estimates, for a total of 3 386 plots of 0.04ha (20×20m square plots)

* SLB: 32,500 living trees from 31 inventories at 22 different sites scattered across the Brazilian Amazon. They were converted into plot-level AGB estimates and associated uncertainty, for a total of 625 plots ranging from 0.16ha (40×40m square plots) to 1ha (20×500m transects),

Regarding CMS data, tree-by-tree data was not accessible but plot-level information was, for a total of 104 plots ranging from 0.1ha (20×50m rectangular plots) to 1ha (40×250m transects) covering both drylands (n=82) and wetlands (n=22).

Then, for each of the four dataset, LiDAR-derived canopy height models (CHMs) were produced for all the lidar point cloud data. Mean top-of-canopy height (TCH) was calculated within the footprint of each tree inventory unit from the corresponding CHM, and a power law model was built between plot AGB and TCH using all the plots units. This model was applied to all LiDAR footprints, and AGB maps and associated uncertainty were produced by Nicolas Labrière.

<table><caption>Table 6.1: Lidar-based AGB data characteristics</caption><thead><tr><td>Source</td><td>Country</td><td>Area in km²</td><td>Nb. lidar footprints</td><td>Nb. sites</td><td>Year</td></tr></thead><tbody><tr><td>TERN</td><td>Australia</td><td>370</td><td>13</td><td>13</td><td>2012-2015</td></tr><tr><td>NEON</td><td>USA</td><td>21 650</td><td>125</td><td>-</td><td>2013-2019</td></tr><tr><td>SLB</td><td>Brasil</td><td>820</td><td>286</td><td>69</td><td>2011-2018</td></tr><tr><td>CMS</td><td>Indonesia</td><td>1 500</td><td>75</td><td>-</td><td>2014</td></tr><tr><td></td><td>Total</td><td>24 340</td><td>499</td><td>-</td><td>2011-2019</td></tr></tbody></table>

### 6.2.4.2 In situ AGB data

We used for evaluation the AGB plots dataset from Ploton et al. (2020a), accessible online to all for free. This dataset consists in the AGB estimations and associated uncertainty for nearly 60,000 1-km pixels aggregated from nearly 100,000 ha of in situ forest management inventories for the 2000-early 2010s period. The data specifically targets five central African countries in the African humid tropical forests, providing detailed insights into this ecologically important region. The uncertainty on AGB estimations derived from 0.5-ha inventory plots (8.6-15.0%) is only moderately higher than the error obtained from 1-ha scientific sampling plots (8.3%). This uncertainty is acceptable, considering that larger plots (here 1-ha versus 0.5-ha) exhibit smaller uncertainties, and that management inventories present higher uncertainties that scientific sampling plots.

### 6.2.4.3 FAO AGB data

The FAO Forest Resource Assessment (FRA) provides a worldwide analysis of forest resources, including data on forest extent, condition, and changes over time, with a specific section dedicated to forest biomass FAO (1988). FAO thus offers a comprehensive overview of the world's forest resources, with key insights into forest AGB. Estimates of AGB and carbon stocks are provided for forests defined as vegetation patches with a minimum tree cover of 10% and minimum area of 0.5 ha. It is important to note that the FAO estimates are based on reports provided by each country, and some countries have limited capacities to estimate their biomass stocks accurately. However, the comparison between our estimates and the FAO estimates provide an idea of the amount of biomass estimated nationally.
### 6.2.4.4 AGB maps from Xu et al. (2021)

Xu et al. (2021) calculated the above- and below-ground biomass of global terrestrial ecosystems on an annual basis, covering the years 2000 to 2019. Reference data, including approximately 100,000 plots coupled with airborne LIDAR data covering more than 1 Mha of tropical forests globally and satellite lidar survey detailing the height structure of global vegetation across over 8 million sample footprints, were used to feed a machine learning model. Reference data were converted into estimates of biomass using established allometric models and root-to-shoot ratio, then used as training data for the machine learning algorithm, together with microwave and optical satellite imagery collected from 2000 to 2019. Note that these maps were used for comparison purpose only.

## 6.3 Methods

### 6.3.1 AGB estimation from multi-angular L-band TBs

The proposed methodology for estimating AGB maps from ANN is based on a purely statistical relationship between inputs, in this case multiangular TBs, and targets, the CCI reference maps. Yet, each polarization of the TBs is grouped into 1 of 14 bins (according to their angle of incidence - See Section 6.2.1). As a result, the number of possible predictors can be as large as 28. In order to determine the optimal set of predictors to retrieve biomass from ANNs, an evaluation of the estimated AGBs was performed. Taking advantage of the multi-angle capabilities of the SMOS mission, the quality of the produced biomass maps was evaluated in two ways. First, the influence of the TBs' angle of incidence on the quality of the AGB inversions was assessed. Then, the influence of different predictors was analyzed, being TBs' H and V the main ones. Other variables such as the Polarization Ratio (PR, Eq. 6.1) and soil temperatures (TS) were used as complementary predictors. Soil temperature was also used to estimate emissivities (E) from the TBs, calculated as $E = TB/TS$; Henceforth, E in horizontal and vertical polarizations will be referred to as $E_h$ and $E_v$, respectively. The polarization ratio is defined as:

$$PR = \frac{TB_V - TB_H}{TB_V + TB_H} \qquad (6.1)$$

Scattering events within the vegetation strongly depolarize the ground emission. Therefore, low PR values are associated with dense vegetation and higher values with sparsely vegetated soils

In addition, a hybrid approach for the estimation of AGB was evaluated. In this approach, the outputs of a physically based radiative transfer model (VOD or SM) were used along with the TBs as complementary predictors in the ANN input vector. To summarize, AGB estimates will be presented from two ANN approaches: the first referred to as 1) the data-driven approach (ANN_D), which is based solely on TBs (and their derived variables), and the second referred to as 2) the hybrid approach, which uses both TBs and VOD (or SM) as inputs for the networks (ANN_H), see Fig 6.1. An additional case using parametric functions based solely on the VOD will be used and discussed later in the document (VOD function).

The learning sets used by the networks are divided in three subsets. The first one, the training subset, is used to train the network and adjust its parameters (weights and biases) in order to reduce the difference between estimates and targets. The validation subset is used to monitor the performances of the ANN at each step of the training process and to detect overtraining if the performances continue to increase on the training subset but start decreasing on the validation subset. The last subset, the test one, is used to evaluate the ability of the ANN to generalize the patterns learned during the training stage to a set of new data. The training, validation and test subsets were built from a random selection of the global observations contained in the learning set, and represent respectively 50 (~58000 pixels/year), 20 (~23000 pixels/year) and 30% (~35000 pixels/year) of the original set.
Figure 6.1: Sketch of pure data-driven approach (AGB(ANN<sub>D</sub>)) using TBs or derived variables (like E or PR) and the hybrid approach (AGB(ANN<sub>H</sub>)) using TBs and variables from radiative transfert models (VOD or SM).

<figure><img src="image_36.png" /></figure>

The optimization of the ANN parameters is done using the Levenberg-Marquardt algorithm, the performance of the network is calculated using the Mean Square Error (MSE). A simple architecture is used: a single hidden layer with 10 neurons with a logistic-sigmoidal activation function, together with an output layer with a linear activation function. Initially, the training is carried out on input-target pairs of the same year. However, thanks to the three available years, a temporal cross-validation is applied. In other words, once a model is trained on data from a specific year, it is then applied to data from another year in order to compare the resulting inversion to the reference map of that period, e.g. an ANN trained with the input-target pair: TBs 2017-CCI2017, is subsequently applied to TBs of 2018, to compare their inversion to the CCI2018 map, and so on for all other possible combinations. An alternative training scheme, based on a multi-year learning-set (a concatenation of inputs and references from the three available years) was explored to assess the flexibility of such an ANN for a temporal cross-validation; the same ANN architecture mentioned before is used, the construction of the learning subsets (training, validation and test) is carried out from a random selection of the three-year concatenated data set.

### 6.3.2 AGB estimation from parametric function for L-band

Currently, the common approach to estimate AGB from passive microwaves is based on determining a relationship between VOD and AGB using a functional fit. Previous studies have explored this methodology and have proposed different parametric functions to estimate AGB from VOD (Liu et al. (2015); Rodríguez-Fernández et al. (2018)). Here, the one proposed by Rodríguez-Fernández et al. (2018) was used.
<figure><img src="image_37.png" /><figcaption>Figure 6.2: Overview of the methodology for estimating yearly AGB maps from SMOS L-VOD and an AGB reference map.</figcaption></figure>

The workflow is divided into three steps numbered from 1 to 3 in Fig. 6.2. The step 1 (pre-processing) refers to the aggregation and interpolation of the AGB calibration map. This was done by Magellium and it is described in Chapter 8.

In step 2 (2- Processing) the SMOS L-VOD daily measurements are masked, filtered, and temporally averaged. Only pixels over continental surfaces are considered, as provided in the 1km USGS (US Geological Survey) land-sea mask aggregated into the EASE 2 grid. Low-quality VOD are removed based on the RMSE between the modeled and measured TB. The SMOS footprints with a RMSE larger than 8K are filtered out. The VOD temporal series of each footprint are also checked for potential outliers: the values outside an interval of two sigmas around the yearly average are discarded. The cleaned VOD time series, are then temporally averaged on a yearly basis to iron out the effects on L-VOD of the diurnal and seasonal variations of the vegetation water content. At this stage, footprints with an AGB reference value equal to 0 are discarded, as they are not useful to our AGB estimation.
In step 3 (3- Fitting in Fig. 6.2), the AGB reference map is then compared pixel-wise against the annual SMOS L-VOD map (red points in the 3-Fitting part of Fig. 6.2) for the same year (for example SMOS L-VOD in 2018 against ESA Biomass CCI AGB map for the year 2018). The annual L-VOD is binned into 0.05-width bin. In each bin, the mean AGB from the reference map is computed (black points in the 3 - Fitting part of Fig. 6.2) and the set of parameters of the logistic function that best fits the mean AGB/L-VOD distribution is estimated. This logistic function is defined in Eq. 6.2:

$$ AGB = \frac{a}{1 + e^{-b(VOD-c)}} + d \qquad (6.2) $$

where a, b, c and d are the free parameters to obtain the logistic function that best fits the mean AGB distribution in the L-VOD bins (black curve in the 4 - Fitting part of Fig. 6.2). In equation 6.2, AGB is in Mgha⁻¹ and the L - VOD is dimensionless. Hence a and d are in Mg ha⁻¹ and b and c have no dimension.

### 6.3.3 AGB estimation from multi-frequency VOD

Multi-frequency approaches to estimate AGB from passive microwaves VOD was studied using a multi-linear regressions, neural networks and random forest.

## 6.4 Results

### 6.4.1 Estimation of AGB from L-band brightness temperatures

#### 6.4.1.1 Quality evaluation and predictors selection

Variations in the R² scores between AGB references and estimates, with respect to the angle of incidence and different predictors, are shown in Fig. 6.3. Several combinations of ANN predictors were tested, using single or grouped variables. Higher R² values were found in retrievals from higher incidence angles (bin >8, incidence angles >35°). In single predictors, the H-polarization and PR presented the highest R² values (around 0.72 and 0.80, respectively). The V and emissivities alone presented non-significant values, around 0.1 for the most part and 0.45 as the highest value (for E_h between 60 to 65° [Bin No. 14]). The highest R² were obtained when multiple predictors were used at the same time. The most performant combinations were the cases j, k and l in Fig. 6.3, all with similar values close to 0.85; the use of E_v and V did not significantly increased the coefficient' scores. The case f, which did not use TBs but only variables derived from TB and TS, and case i, which used only TBs and PR, showed similar scores for angles on incidence >35° (R² ~0.82). In inversions with multiple predictors, the less performant approach was that of case g, with slightly lower R² scores, at around 0.81.
Figure 6.3: R² values for AGB estimates depending on angle of incidence (Y axis) and predictors (X axis). Bins from 1 to 14 covering angles of incidence between 0 to 65º, in intervals of 5º -see Section 6.2.1-. For E the subindex represents the polarization: E<sub>h</sub> corresponds to the horizontal polarization, and E<sub>v</sub> to the vertical one. The inversions presented were produced from 2017 data-sets using the CCI2017 map as a reference.

<figure><img src="image_38.png" /></figure>

Thus, from two of the most performant groups of predictors: case *j*, with the highest values of R², and case *i* with relatively lower results but using exclusively TBs, the quality of AGB estimates (in terms of R²) was evaluated for different ranges of angles of incidence on a global and continental scale (see Fig. 6.4). AGB estimates obtained from ANN<sub>D</sub>(H, PR, E<sub>h</sub>) were evaluated for: all angles (14 angle bins and 3 variables, therefore 42 input predictors), angles between 35 and 65º [bins 8 to 14], and between 40 and 45º [bin number 10] (inputs a, b and c respectively in Fig. 6.4). Inversion from ANN<sub>D</sub>(H, V, PR) were evaluated for: angles between 35 and 65º [bins 8 to 14], and between 40 and 45º [bin number 10] (inputs d and e in Fig. 6.4). Three additional sets of predictors were studied, two for the hybrid approach (using VOD and SM as complementary predictors to H and PR) and one using TS due to its strong influence in TB observations (inputs f, g and h in Fig. 6.4).

Estimates from a single angle of incidence had the lowest R² values among all combinations. However, using one angle with the predictors (H,PR,E<sub>h</sub>) presented similar scores to those produced from angles 35 and 65º [bin 8 to 14] with (H,V,PR), around 0.86 in both cases. On the global scale, the most performant values come from the combination a, b and f in Fig. 6.4, with an R² of ~0.89. Finally, the hybrid case using SM (and TS as complementary predictor) showed similar coefficients (close to 0.87).

Figure 6.4: Influence of different angles of incidence on the R² of the best predictors. The label '8:14' represents the TBs between 35 and 65º (see Section 6.2.1). The '1:14' tag is equivalent to 0 to 65º. The inversions presented were produced using data-sets from 2017 and the CCI2017 map as a reference.

<figure><img src="image_39.png" /></figure>

A similar analysis to that of Fig. 6.3 and Fig. 6.4 was done for 2010 and 2018 (at a global scale and for Africa and South America), similar results were obtained. Consequently, henceforth the following groups of
predictors will be used: i) the optimal combination of predictors (H,PR,E<sub>h</sub>), as it gives the best compromise between angles of incidence and higher R<sup>2</sup> scores of the estimates, ii) the combination (H,V,PR), that despite having relatively low values for the coefficient of determination, exhibits the results that can be expected from multi-angular TBs alone, and iii) both hybrid cases, as they presented considerable R<sup>2</sup> values without the use of emissivities; Despite having similar results to the latter, the combination that includes TS as complementary predictor will not be used, since it is considered that the effect of the surface temperature on the AGB estimates was better captured from the emissivities (derived from TB and TS).

Model inversions from the selected cases were evaluated in both their ability to explain the variability of the reference AGB maps (R²) and their dispersion with respect to them (RMSE). The first three matrices (from top to bottom) of Fig. 6.5, present the R² values obtained for the AGB inversions produced for the three reference years. In Fig. 6.5 the title of each matrix represents the year of the data-set that was used for the training of models, therefore, the second matrix titled '2017' contains the R² of all the models trained using input variables from 2017 as predictors and the CCI2017 as reference. The columns, on the other hand, represent the year of the data-set [reference map] to which the trained model was applied to [compared to], e.g. the value of the upper right corner of the '2017 matrix' (so 0.87), represents the R² between: 1) the CCI2018 reference map and 2) the AGB estimate from the ANN<sub>D</sub>(H,PR,E<sub>h</sub>) trained with 2017 data and applied to data-sets from 2018.

On average, when the AGB estimates produced by ANN<sub>D</sub>(H,V,PR) were evaluated with respect to their respective reference map, a R² of around 0.85 units was obtained; this approach presented slightly lower R² values than the other cases. When looking at the remaining ANNs and the temporal extrapolation of models, two patterns were observed. First, both ANN<sub>H</sub>(H,PR,VOD) and ANN<sub>D</sub>(H,PR,E<sub>h</sub>) presented similar performances (R² ~0.87) for each year. However, the use of E<sub>h</sub> or VOD as a complementary predictor yielded relatively (but not significantly) higher values than the use of SM; the inclusion of horizontal polarized emissivities as an additional predictor of ANNs slightly increased the R² value (in ~0.03 units) with respect to the ANN<sub>D</sub>(H,V,PR). Second, the temporal cross-validation of models showed a decrease in R² values. Models trained in 2010 and applied to 2018 reduced their R² by about 0.2 units, when applied to 2017 data-sets the difference was not significant. A similar decrease was observed when applying models trained in 2018 to data from 2010 (~0.2). The difference in R² between 2017 and 2018 models was not significant.

Similar patterns were observed in the RMSE values (not shown here), i.e. a) an improvement from the inclusion of E<sub>h</sub> (and similar results from the VOD hybrid case) and b) a deterioration when applying models to data-sets from years different than those used during training. AGB estimates from ANNs trained using (H,V,PR) from 2018, 2017 and 2010 showed RMSE in the order of 38.7 [54%], 38 [54%] and 38.6 Mg/ha [55%] respectively. The ANN<sub>D</sub>(H,PR,E<sub>h</sub>) case reduced the previous values by 3.7 [4%], 3.8 [4%] and 2.4 Mg/ha [5%] respectively (similar decreases were observed with the ANN<sub>H</sub>(H,PR,VOD)).

Finally, ANNs were trained with learning-sets composed of the three reference years (multi-year learning-set mentioned at the end of Section 6.3.1); the temporal cross-validation of models trained under these conditions is shown in the '3 years' matrix in Fig. 6.5 (bottom matrix). The evaluation is done using only test subset (Sect. 6.3.1). The ANNs trained under this configuration, as in the previously mentioned cases, showed slightly higher results from both the ANN<sub>D</sub>(H,PR,E<sub>h</sub>) and the ANN<sub>H</sub>(H,PR,VOD) approach (increase of ~0.03 units) than from the ANN<sub>D</sub>(H,V,PR). Additionally, a significant improvement in the cross-validation of models trained with multiple years was observed. The aforementioned interannual differences in R<sup>2</sup> (0.2 for models trained in 2010 [2018] and applied to 2018 [2010] data-sets), were not significant under the 3-years training configuration (<0.01 units). When comparing the R<sup>2</sup> (and RMSE) from multi-year trained models with the results from single-year models (trained and applied on data from the same year), a slight (but not significant) decrease in R<sup>2</sup> was observed in some cases. However, to carry out the evaluation of AGB estimates, the models trained with a multi-year learning set (concatenation of three years) were chosen, since they are the ones that showed greater stability for inter-annual analyses.

### 6.4.1.2 Evaluation of AGB estimates from ANNs

Once the quality metrics for the different cases were evaluated, and the optimal training configuration defined, the spatial distribution of residuals was explored (that is, the differences between references and estimates). This analysis was performed exclusively on three cases of interest, more precisely: ANN<sub>D</sub>(H,PR,E<sub>h</sub>),
<figure><img src="image_40.png" /><figcaption>Figure 6.5: Temporal cross-validation of models (in terms of R²). The title of each matrix represents the year of the data-set used for training the models. On the Y axis of each matrix: The models used to estimate AGB. On the X axis : The years of the datasets to which the trained models were applied to.</figcaption></figure>
$ANN_D(H,V,PR)$ and $ANN_H(H,PR,VOD)$; the results by the $ANN_H(H,PR,SM)$ presented lower quality metrics to those produced by using VOD as an additional predictor, only the results of the latter are shown.

<figure><img src="image_41.png" /><figcaption>Figure 6.6: Comparison of the spatial structure of AGB maps. Inversions are produced by models trained with a multi-year data-set and applied to 2018 observation, the reference is the CCI2018 map. AGB is regrouped in bins of 35 Mg/ha of width.</figcaption></figure>

The evaluation of the spatial patterns of AGB estimates produced from the selected $ANN_D$ and $ANN_H$ over tropical Africa and America are shown in Fig. 6.6. In order to represent the spatial distribution of vegetation biomass, the AGB was regrouped into bins of 35 Mg/ha of width. In the left column, and in descending order, are the AGB maps from: CCI2018, the data-driven cases $ANN_D(H,V,PR)$ and $ANN_D(H,PR,E_h)$, and the hybrid case $ANN_H(H,PR,VOD)$ (the final panel will be discussed in Section 6.4.1.3). In the right column the count per AGB bin. All inversions are from models trained with the multi-year data-set and applied to 2018 observations.

Just like shown by the other quality metrics (bottom box in Fig. 6.5), both the $ANN_D(H,PR, E_h)$ and the $ANN_H(H,PR,VOD)$ presented similar results. However, none of the three cases was capable of reproducing
the highest AGB values, e.g. the redish spots in the Amazon rain-forest (around [Lat: -10, Lon: -70], or [Lat: 0, Lon: -55]), where models reproduced values around 385 Mg/ha whilst the reference show AGB values around ~500 Mg/ha. Similarly, even if the main structural features of reference maps were captured by ANNs, e.g. lower values in northern and southern tropical Africa, and higher values in the more equatorial region of both America (Amazon rainforest) and Africa (Congo forest), subtle spatial variations like the transition from densely to sparsely vegetated areas in the Amazon rainforest (region between Lat: [-10, 0], Lon: [-70, 60]) were not reproduced accurately. Nevertheless, the cluster of high AGB values from the reference map, located near the Atlantic Ocean on the tropical America, were more closely reproduced by $ANN_D(H,PR,E_h)$ and $ANN_D(H,V,PR)$, than by the hybrid approach. Regardless, ANNs tended to correctly replicate the main features of the distribution of the reference map (histograms on the right side of Fig. 6.6).

### 6.4.1.3 Comparison of estimates from ANN and VOD fitted functions

In view of the demonstrated potential of the methodology proposed here to estimate AGB, a comparison between AGB estimates produced from ANNs with those produced from the fitted VOD function was finally explored. The same evaluation criteria applied to the ANNs were applied to the VOD-based estimates. To ensure a fair comparison, the training of models and the comparison of estimates were performed on the same data-sets for both methods.

The fifth row of each matrix in Fig. 6.5 ($R^2$ results), gives a comparison of the performance of the VOD function with respect to ANNs. The $R^2$ scores between references and AGB estimates from the VOD functions presented similar results independently of year used for the training. $R^2$ values from $ANN_D(H,PR, E_h)$ are generally ~0.03 units higher than those from the VOD parametric function; a similar difference was observed when compared to $ANN_H(H,PR,VOD)$. The RMSE values associated with the inversion of the VOD function were around 2.8 Mg/ha larger than those on the $ANN_H(H,PR,VOD)$, the difference with respect to the $ANN_D(H,PR,E_h)$ case was on the order of ~2 Mg/ha.

The lower panels of Fig. 6.6 show the spatial patterns and histogram of AGB estimates produced from VOD. Estimates of the parametric function tended to simplify the spatial variations of biomass. In regions with higher AGB densities, the estimates of the parametric function deviated from the references more strongly than those from ANNs, resulting in inversions larger than the reference AGB in the Amazon rain forest while lower than the AGB reference in tropical Africa. Furthermore, while it is true that both methods (the ANNs and the parametric functions) underestimated the highest AGB values (Fig. 6.6), it is also true that the ANNs were closer to the latter than the VOD function. In general, the AGB reference distribution was reproduced more accurately by the ANNs (right-side Figures in Fig. 6.6). For example, the Amazon rain-forest was mostly represented by values between 315 and 385 Mg/ha under the VOD-based models, while according to the estimates of the reference, it can encompass densities between 245 and >500 Mg/ha.

One might think that a more complex approach such as using an ANN with VOD as only input instead of a parametric function could give better performances, however an ANN with VOD as input did not change significantly the results obtained with the parametric function.

## 6.4.2 L-AGB OSMOSE

Following the methodology described in Sect. 6.3.2, a yearly AGB time series has been estimated from the project L-VOD dataset. An example of the L-AGB in 2018 is displayed in Fig. 6.7. This L-AGB has been generated using the year 2018 as a reference. This means that the CCI calibration map reprojected to EASE 2 for the year 2018 and the 2018 L-VOD average, mixing ascending and descending orbits, have been used to calibrate the parameters of the logistic function (Eq. 6.2). The maximum AGB estimates are in the equatorial belt and reach 300 Mg ha⁻¹. The boreal forests, in the northern high latitudes, present AGB estimates around 100-150 Mg ha⁻¹. Temperate and arid regions show an estimated AGB lower than 50 Mg ha⁻¹. The spatial distribution is coherent with CCI 2018 (see Fig. 6.8) and with the input L-VOD.
<figure><img src="image_42.png" /><figcaption>Figure 6.7: AGB estimated from the OSMOSE L-Band VOD, mixing ascending and descending orbits</figcaption></figure>

Panels (a) to (d) in Fig. 6.8 show the spatial distribution of the differences across the four years (2017-2020) of the CCI reference maps. The maps are very similar from one year to the other. The region with the highest differences (more than 75 pixels with differences greater than 200 Mg ha⁻¹) is the equatorial part of Africa. The statistics (see Table 6.2) are also steady over the years. The AGB estimated from OSMOSE L-VOD is very well correlated to the CCI maps from 2017 to 2020 (R above 0.88). On the global scale, the OSMOSE derived L-AGB shows similar distribution, presents more variability, and slightly lower estimations compared to the CCI. Table 6.2 shows a global bias of ≈ -2.5 Mg ha⁻¹ and Fig. 6.7 highlights that this negative bias is particularly high over the Congo basin and the boreal forests of eastern Europe and Russia.

<figure><img src="image_43.png" /><figcaption>Figure 6.8: (a)-(d) Difference in Mg ha⁻¹ between the yearly AGB estimated from the 2018 global logistic fit minus CCI AGB values for the years 2017-2020.</figcaption></figure>
<table><caption>Table 6.2: R, bias (Mg ha⁻¹) and RMSE (Mg ha⁻¹) between the estimated AGB from OSMOSE LVOD and AGB values from CCI for the years 2017-2020. AGB was estimated using the global logistic fit computed from the VOD and the CCI AGB for the year 2018. *: year of reference.</caption><thead><tr><th></th><th>2017</th><th>2018*</th><th>2019</th><th>2020</th></tr></thead><tbody><tr><th>R</th><td>0.88</td><td>0.89</td><td>0.89</td><td>0.88</td></tr><tr><th>bias</th><td>-2.46</td><td>-2.42</td><td>-2.64</td><td>-2.26</td></tr><tr><th>RMSE</th><td>36.30</td><td>35.89</td><td>35.41</td><td>35.93</td></tr></tbody></table>

Overall, on a yearly basis, the AGB generated from the OSMOSE L-VOD compares well with the optical indices from the data cube 2.

### 6.4.3 AGB evaluation

#### 6.4.3.1 AGB evaluation method

The major problem inherent in comparing SMOS AGB maps and in situ data lies in the drastic difference in spatial resolution between these dataset. In fact, we are comparing SMOS AGB maps at 25 m pixel size with in situ data of pixel size ranging from 100 m (lidar-based AGB maps) to 1km (data from Ploton et al. (2020a)). The surface area of in situ data located in a single SMOS pixel therefore represents a very small fraction of the SMOS pixel surface area itself. In this particular case, in situ data are not necessarily representative of the SMOS pixel. For example, a SMOS pixel embedded with 50% of dense forest (mean AGB of 300 Mg/ha) and open forest (mean AGB of 100 Mg/ha) has a mean biomass of 200 Mg/ha, to be compared with in situ data located either in dense forest or in open forest. To account for this bias, we need to account for the representativity of the in situ plot on the ground.

To do so, we corrected in situ data AGB values using the 2010 tree cover density map at 30m pixel size from Hansen et al. (2013), as detailed in Figure 6.9. We first updated the 2010 tree cover map to 2020 by setting to 0 all the pixels considered as forest loss from 2010 to 2020 in the forest loss maps from Hansen et al. (2013). From the 30 m pixel size 2020 tree cover map, we then derived 100 m and 1 km tree cover maps (average interpolation) to fit with reference in situ data, and 25 km pixel size tree cover map by averaging all the 30m pixels located in the 25 km SMOS AGB pixels. We finally computed the adjusted in situ AGB pixels as follows:

$$AGB_{ref,adj_i} = AGB_{ref_i} \cdot \frac{E(tc_{smos})}{tc_{ref_i}} \qquad (6.3)$$

with AGB<sub>ref,i</sub> the AGB reference (in situ) value at pixel i, E(tc<sub>smos</sub>) the 25 km tree cover pixel and tc<sub>ref,i</sub> the 100m or 1km tree cover value at pixel i.

We finally performed the validation process by comparing AGB pixel values at 25km with adjusted reference in situ data values AGB<sub>ref,adj_i</sub>.
<figure><img src="image_44.png" /><figcaption>Figure 6.9: Overview of the methodology for the evaluation of yearly AGB maps from SMOS L-VOD and AGB in situ data described in section 6.2.4.</figcaption></figure>

### 6.4.3.2 AGB evaluation results

#### Lidar-based AGB data

The validation results associated with all lidar-based AGB maps described in table 6.1, compared with 2020 SMOS AGB from VOD calibrated with CCI AGB map are presented in figure 6.10. Results associated to CMS (Indonesia), NEON (USA), SLB (Brazil) and TERN (Australia) data are shown in green, red, blue and yellow respectively. The accuracy assessment resulted in an R value of 0.59, a bias of -58.5 Mg/ha and a RMSE of 104.8 Mg/ha. Considering AGB values from lidar maps greater than 100 Mg/ha, 150 Mg/ha and 200 Mg/ha lowered the validation results. The R value decreased from 0.56 to nearly 0, the bias from -58.5 Mg/ha to -122 Mg/ha and the RMSE increased from 104.8 Mg/ha to 158.1 Mg/ha when considering AGB values from lidar maps greater than 200 Mg/ha.

The TERN network demonstrated the best performance in terms of correlation, achieving an R value of 0.86. However, the bias and RMSE were found to be large, that is to say -145.5 Mg/ha and 205.5 Mg/ha respectively, showing a consequent under-estimation of AGB estimated from SMOS. The CMS network exhibited a R value of 0.74, a small bias of -1.9 Mg/ha and a RMSE of 67.2 Mg/ha, the NEON network a R value of 0.58, a bias of -68.4 Mg/ha and a RMSE of 93.8 Mg/ha, and the SLB network a R value of 0.51, a bias of -65.8 Mg/ha and a RMSE of 111.6 Mg/ha. As an example, the detailed results related to the CMS network, for all years, is shown in figure 6.11.
<figure><img src="image_45.png" /><figcaption>Figure 6.10: Validation using all lidar-based AGB vs 2020 SMOS AGB calibrated with CCI maps.</figcaption></figure>

<table><caption>Table 6.3: R, bias, MAE and RMSE between the estimated AGB from OSMOSE LVOD calibrated with CCI maps, and lidar-based AGB maps >0 Mg/ha, >100 Mg/ha, >150 Mg/ha and >200 Mg/ha.</caption><thead><tr><th></th><th>All</th><th>AGB>100 Mg/ha</th><th>AGB>150 Mg/ha</th><th>AGB>200 Mg/ha</th></tr></thead><tbody><tr><td>nb. sample</td><td>441</td><td>317</td><td>232</td><td>147</td></tr><tr><td>R</td><td>0.56</td><td>0.35</td><td>0.22</td><td>-0.04</td></tr><tr><td>bias</td><td>-58.5</td><td>-79</td><td>-97.9</td><td>-122</td></tr><tr><td>MAE</td><td>76.4</td><td>91.2</td><td>105.4</td><td>127.4</td></tr><tr><td>RMSE</td><td>104.8</td><td>119.8</td><td>134.8</td><td>158.1</td></tr></tbody></table>
<figure><img src="image_46.png" /></figure>

<table><thead><tr><th>CMS - Indonesia</th><th>R</th><th>R<sup>2</sup></th><th>bias (%)</th><th>MAE (%)</th><th>RMSE (%)</th><th>bias (Mg/ha)</th><th>MAE (Mg/ha)</th><th>RMSE (Mg/ha)</th></tr></thead><tbody><tr><td>CCIv4_2011</td><td>0.70</td><td>0.49</td><td>38.10</td><td>59.33</td><td>125.21</td><td>-0.86</td><td>55.51</td><td>71.94</td></tr><tr><td>CCIv4_2012</td><td>0.72</td><td>0.52</td><td>43.03</td><td>62.75</td><td>133.55</td><td>6.24</td><td>56.35</td><td>70.11</td></tr><tr><td>CCIv4_2013</td><td>0.75</td><td>0.56</td><td>42.60</td><td>59.65</td><td>126.61</td><td>10.54</td><td>53.24</td><td>67.54</td></tr><tr><td>CCIv4_2014</td><td>0.74</td><td>0.55</td><td>48.72</td><td>64.23</td><td>142.72</td><td>12.79</td><td>54.06</td><td>68.83</td></tr><tr><td>CCIv4_2015</td><td>0.72</td><td>0.51</td><td>40.06</td><td>59.92</td><td>123.89</td><td>8.18</td><td>56.03</td><td>70.32</td></tr><tr><td>CCIv4_2016</td><td>0.75</td><td>0.57</td><td>26.90</td><td>49.28</td><td>92.19</td><td>5.11</td><td>54.10</td><td>67.37</td></tr><tr><td>CCIv4_2017</td><td>0.73</td><td>0.54</td><td>34.27</td><td>55.09</td><td>111.36</td><td>5.59</td><td>54.38</td><td>68.89</td></tr><tr><td>CCIv4_2018</td><td>0.72</td><td>0.51</td><td>42.73</td><td>61.26</td><td>131.14</td><td>9.39</td><td>55.23</td><td>70.49</td></tr><tr><td>CCIv4_2019</td><td>0.72</td><td>0.51</td><td>44.12</td><td>61.75</td><td>135.64</td><td>8.66</td><td>54.23</td><td>70.44</td></tr><tr><td>CCIv4_2020</td><td>0.74</td><td>0.55</td><td>26.00</td><td>50.31</td><td>98.18</td><td>-1.91</td><td>52.69</td><td>67.23</td></tr><tr><td>CCIv4_2021</td><td>0.75</td><td>0.56</td><td>35.15</td><td>53.97</td><td>111.61</td><td>5.55</td><td>51.90</td><td>66.51</td></tr><tr><td>CCIv4_2022</td><td>0.74</td><td>0.54</td><td>33.76</td><td>53.55</td><td>109.29</td><td>4.49</td><td>52.60</td><td>67.92</td></tr></tbody></table>

Figure 6.11: Validation using lidar-based AGB from the CMS network in Indonesia vs SMOS AGB calibrated with CCI maps.

## In situ AGB data from Ploton et al. (2020a)

Validation results using in situ AGB data and SMOS AGB calibrated with CCI are presented in Figure 6.12. Compared to lidar-based AGB data, results showed very low correlation with $r^2$ from 0.03 to 0.08, a strong under-estimation of AGB estimates (bias ranging from -66.8 to -90.1 Mg/ha) and a RMSE from 112.8 up to 133.6 Mg/ha. Note that validation results when SMOS VOD is calibrated with Avitabile's map does not show better correlation ($r^2$ < 0.8), although the bias is found a bit smaller (from -25.6 to -48.5 Mg/ha). Finally, SMOS maps directly based on Tbs showed slightly better results ($r^2$ = 0.9 and bias from -22.6 to -1.6 Mg/ha). We don't think that SMOS AGB estimates are poorer over Africa. Instead, AGB estimates from commercial inventories in Ploton et al. (2020a), at 1km pixel size but originally with a 0.5ha resolution, are probably not representative enough. This hypothesis highlights the need of lidar-based AGB data calibrated with in situ tree inventories for an accurate statistical assessment of SMOS AGB maps.
<figure><img src="image_47.png" /></figure>

<table><thead><tr><th></th><th>r</th><th>r²</th><th>bias (%)</th><th>MAE (%)</th><th>RMSE (%)</th><th>bias (Mg/ha)</th><th>MAE (Mg/ha)</th><th>RMSE (Mg/ha)</th></tr></thead><tbody><tr><td>CClv4_2011</td><td>0.26</td><td>0.07</td><td>-25.31</td><td>34.77</td><td>41.43</td><td>-90.11</td><td>104.15</td><td>130.54</td></tr><tr><td>CClv4_2012</td><td>0.27</td><td>0.07</td><td>-23.93</td><td>34.16</td><td>40.93</td><td>-86.73</td><td>101.90</td><td>127.62</td></tr><tr><td>CClv4_2013</td><td>0.28</td><td>0.08</td><td>-21.51</td><td>32.27</td><td>38.87</td><td>-81.00</td><td>96.89</td><td>122.21</td></tr><tr><td>CClv4_2014</td><td>0.29</td><td>0.08</td><td>-17.87</td><td>30.78</td><td>38.02</td><td>-71.53</td><td>90.97</td><td>115.59</td></tr><tr><td>CClv4_2015</td><td>0.24</td><td>0.06</td><td>-18.48</td><td>31.52</td><td>38.96</td><td>-73.58</td><td>93.36</td><td>119.28</td></tr><tr><td>CClv4_2016</td><td>0.27</td><td>0.08</td><td>-16.71</td><td>30.78</td><td>38.50</td><td>-68.66</td><td>90.06</td><td>114.91</td></tr><tr><td>CClv4_2017</td><td>0.18</td><td>0.03</td><td>-17.50</td><td>32.79</td><td>41.00</td><td>-72.42</td><td>95.96</td><td>122.10</td></tr><tr><td>CClv4_2018</td><td>0.17</td><td>0.03</td><td>-24.04</td><td>35.78</td><td>43.27</td><td>-88.71</td><td>106.21</td><td>133.60</td></tr><tr><td>CClv4_2019</td><td>0.26</td><td>0.07</td><td>-21.75</td><td>33.49</td><td>40.57</td><td>-82.12</td><td>99.22</td><td>123.97</td></tr><tr><td>CClv4_2020</td><td>0.25</td><td>0.06</td><td>-20.21</td><td>32.38</td><td>39.67</td><td>-78.33</td><td>96.08</td><td>121.33</td></tr><tr><td>CClv4_2021</td><td>0.26</td><td>0.07</td><td>-18.86</td><td>31.75</td><td>39.16</td><td>-74.22</td><td>93.60</td><td>118.83</td></tr><tr><td>CClv4_2022</td><td>0.28</td><td>0.08</td><td>-15.92</td><td>30.28</td><td>37.85</td><td>-66.80</td><td>88.58</td><td>112.84</td></tr></tbody></table>

Figure 6.12: Validation using in situ AGB data from Ploton et al. (2020a) vs SMOS AGB calibrated with CCI maps.

FAO AGB data

Comparison of SMOS AGB calibrated with CCI data with FAO national statistics is presented in Figure 6.13. Correlation coefficient *r*² is comprised between 0.85 to 0.90. Estimates are slightly over-estimated for some countries, in particular where biomass stocks are very low in North Africa (Lybia, Algeria, Tunisia, Egypt, Niger, Mali) and Central Asia (Kazakstan, Turkmenistan, Afghanistan). Our estimates suggest a somewhat higher overall carbon stock than what the FAO reports, which can be attributed to the fact that we provided AGB estimates for entire countries. In contrast, the FAO focuses solely on carbon stocks within forests, defined as areas with at least 10% tree cover and spanning no less than 0.5 ha. This discrepancy highlights the potential carbon stock overlooked when limiting the scope to forests alone. Additionally, it is worth mentioning that FAO's Figures are derived from national reporting, and the capability to precisely gauge carbon stocks varies significantly across countries. In general, countries with limited capabilities for AGB estimation have no means of estimating their carbon stocks, and declare inaccurate values. Countries showing the highest values of AGB estimates from SMOS VOD and FAO are Brazil (almost the same estimates), Russia, Democratic Republic of the Congo, Indonesia, Colombia and China. SMOS AGB estimates calibrated using Avitabile's map are more correlated to FAO data (*r*² is comprised between 0.83 to 0.94) but are a bit more biased. Finally, SMOS AGB estimates using Tbs are the most correlated to FAO data (*r*² = 0.95) but also the most biased.
<figure><img src="image_48.png" /><figcaption>Figure 6.13: Validation using national FAO data vs SMOS AGB calibrated with CCI maps.</figcaption></figure>

## AGB maps from Xu et al. (2021)

We compared SMOS AGB estimates from VOD calibrated with CCI maps with AGB maps from Xu et al. (2021) at the national scale, and showed comparison results in Figure 6.14. We found that both maps were very correlated with $r^2$ from 0.86 to 0.96 depending on the year. Estimates from SMOS are a bit over-estimated compared with the maps from Xu et al. (2021). We found similar results when SMOS AGB were calibrated using Avitabile's map ($r^2$ from 0.91 to 0.94 with a higher positive bias) or estimated using Tbs ($r^2$ from 0.95 to 0.96 with a higher positive bias as well). When comparing directly the map from Xu et al. (2021) with FAO, comparison results were the best with $r^2$ from 0.91 to 0.92 and a negative bias. According to these results, it is difficult to assess which one of AGB maps from SMOS or from Xu et al. (2021) is the best one, but both provide excellent comparison results with FAO data.
<figure><img src="image_49.png" /><figcaption>Figure 6.14: Comparison of SMOS AGB calibrated with CCI maps with AGB maps from Xu et al. (2021) at the national scale.</figcaption></figure>

### 6.4.4 Estimation of AGB from multi-frequency VODs

The three methods mentioned in Sect. 9.3 (multi-linear regression, a feed-forward neural network and a random forest algorithm) were applied to multi-frequency VOD data with and without NDVI.

The Pearson correlation between estimated AGB and the reference increases from multi-linear regression (0.87) to neural networks (0.9) and random forest (0.92). A similar trend is found for Spearman $\rho$. Numerically, very similar results in terms of correlation were obtained when also adding NDVI for the multi-linear regression and the neural networks.

AGB estimations for different years were done using the multi-frequency approach establishing the relationship in 2018 using the CCI2018 AGB map. Maps were also computed with the different methodologies and the different predictors sets with and without NDVI. A summary of the results is shown in Fig. 6.15-6.16. The best results are found for the random forest approach including NDVI as predictor. Table 6.4 shows different quality metrics for other years than 2018, confirming these results.
<figure><img src="image_50.png" /><figcaption>Figure 6.15: Maps of the AGB estimation for year 2018 and differences with respect to the CCI map using L-, C- and X-VOD a multi-linear regression (left), a neural network approach (right).</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 6.16: Maps of the AGB estimation for year 2018 using a random forest approach with L-, C-, X-VOD as predictors (left) and also NDVI (right).</figcaption></figure>
<table><caption>Table 6.4: Comparison of the performances of different machine learning approaches for the estimation of AGB using multi-frequency data. Multi-linear regression (MultiReg), Neural Network (NN) and Random Forest (RF).</caption><thead><tr><th>Method</th><th>2017</th><th>2018</th><th>2019</th><th>2020</th></tr></thead><tbody><tr><td>MultiReg</td><td>0.88</td><td>0.88</td><td>0.88</td><td>0.88</td></tr><tr><td>NN</td><td>0.91</td><td>0.91</td><td>0.91</td><td>0.90</td></tr><tr><td>RF</td><td>0.92</td><td>0.95</td><td>0.92</td><td>0.90</td></tr></tbody></table>

<table><caption>Pearson correlation coefficient of estimated AGB LCXNDVI with CCI</caption><thead><tr><th>Method</th><th>2017</th><th>2018</th><th>2019</th><th>2020</th></tr></thead><tbody><tr><td>MultiReg</td><td>36.35</td><td>37.17</td><td>36.81</td><td>36.69</td></tr><tr><td>NN</td><td>32.26</td><td>32.59</td><td>32.57</td><td>34.04</td></tr><tr><td>RF</td><td>30.22</td><td>24.5</td><td>30.55</td><td>33.19</td></tr></tbody></table>

<table><caption>RMSE (Mg ha⁻¹) of estimated AGB LCXNDVI with CCI</caption><thead><tr><th>Method</th><th>2017</th><th>2018</th><th>2019</th><th>2020</th></tr></thead><tbody><tr><td>MultiReg</td><td>2.85</td><td>3.43</td><td>3.30</td><td>3.92</td></tr><tr><td>NN</td><td>0.03</td><td>0.31</td><td>0.43</td><td>-1.08</td></tr><tr><td>RF</td><td>-0.37</td><td>-0.04</td><td>-0.36</td><td>-2.1</td></tr></tbody></table>

<caption>Bias (Mg ha⁻¹) of estimated AGB LCXNDVI with CCI</caption>

## 6.5 Conclusion

This chapter discusses different ways to estimate AGB from L-band data:

* Using SMOS L-VOD and fitting a logistic function

* Using SMOS brightness temperatures and training a neural network

* Using SMOS L-VOD complemented by X- and C-band VOD and NDVI by multi-linear regression, neural networks and random forest

The first two approaches were evaluated extensively against a large number of AGB data from in situ plots, LIDAR estimations and national inventories.

Most of the published works using AGB estimations from L-band data use L-VOD and a reference AGB map to fit a parametric function. The first alternative methodology for the estimation of AGB from L-band passive microwave observations, which directly exploits the statistical relationships between biomass and multi-angle TB observations using Artificial Neural Networks. The results show that the estimation of AGB from the proposed method presented certain advantages with respect to the method using parametric functions used by most of the works published to date on the subject. The estimated AGB maps, presented a slightly higher correlation than the current methods, with the reference AGB maps (R² 0.87 instead of 0.84), a stronger agreement with the structural patterns of the reference maps over tropical regions, residual maps with lower absolute values and smaller uncertainties when compared to the state of the art AGB maps. The results suggest that accurate biomass estimates can be produced directly from SMOS L-Band multi-angular TBs. Being the optimal set of predictors for the estimation of AGB directly from TBs: the TBs in H polarization, the PR and the E<sub>h</sub> (for angles of incidence >35°).

In addition, it has been shown that using a single year to establish the relationship of the input variables and AGB introduces significant uncertainties when this relationship is used to study the evolution of AGB/Carbon Stocks. An alternative training scheme was proposed based on multi-year learning sets, which accounted for the differences mentioned above and produced more stable estimates for different years. Despite all of the
above, a more detailed study should be carried out on the possible synergies that the observations of different instruments (frequencies) could have with other canopy components, and on the potential applications that this new methodology may have in the future, like the estimation of carbon stocks or AGB temporal evolution.

The evaluation against a large number of in situ estimations of AGB show a good agreement (R ~ 0.6) with the in situ estimations extended using LIDAR tree height data. The estimations using SMOS brightness temperatures give slightly better results than those starting from LVOD. In contrast, the correlation with respect to in situ plots estimations by Ploton et al. (2020b) is very low. This is likely due to the fact that those estimations provided at a 1 km spatial scale come actually from 0.5 ha plots and they are probably not representative of the spatial scale of SMOS observations. It is interesting to remark that the comparison of AGB estimations from SMOS at the national level show a very good agreement with FAO inventories (R ~ 0.9), probing that when compared at large scale, AGB estimations from SMOS are very accurate.

Finally, different methods were studied to estimate AGB from multi-frequency VOD with or without NDVI. Random forest showed the best performances of the three methods studied. They also showed hints of providing better AGB estimates than L-VOD alone. More research on this point would be needed to explore the full potential of this promising approach.
# Chapter 7

## Optical indices (WP 6-7)

### 7.1 Introduction and purpose

The purpose of this section is to describe the development of the algorithm used to estimate the vegetation water content from optical sensors, as well as the corresponding validation results. The amount of vegetation water content can be assessed through several variables, such as leaf water content (LWC, unitless), which is defined as the total mass of water related to the total fresh or dry biomass of the leaves, and the Equivalent Water Thickness (EWT, g.cm-2) that is defined as the mass of water contained in a unit of leaf area *Ceccato et al.* (2001). At the canopy level, the amount of water in the vegetation can be related to any vegetative living organ (stem, trunk, leaves, reproductive organs etc.) or solely green leaves. Vegetation water content can thus be assessed through several variables such as live fuel moisture content (LFMC, %) *Yebra et al.* (2013) or canopy water content (CWC, gcm-2) also named vegetation water content (VWC, gm-2) *Hornbuckle* (2023). CWC correspond to the product of the canopy Plant Area Index (PAI, e.g. one-sided vegetation element area per unit horizontal ground area, unitless) and EWT, and are related to all the living vegetation elements *Jackson et al.* (1982).

We used MODIS, SENTINEL-2 and LANDSAT-7 Data to derive the CWC products. The validation lies in two steps, first validate the CWC product at the SENTINEL-2 scale against ground measurements that are easier to obtain, and then inter-compare the CWC MODIS product with aggregated S2 CWC. The following represents the activities developed for the tasks (work packages) 6 and 7 of the project (see Chapter 3).

### 7.2 Review of methods to estimate Vegetation traits from reflectance observations (WP6)

WP6 consisted in reviewing methods for estimating different vegetation traits, including CWC. Since the 80's, considerable improvement in the quality of vegetation traits derived from satellite or airborne systems has been achieved due to the advances of measurement capability of satellite instruments and to our understanding of the radiation regime within vegetation canopies *Liang* (2005). However, remote sensing observations sample a limited fraction of the radiation field reflected or emitted by the surface, and thus do not provide directly estimates of the targeted traits. It is therefore necessary to transform the radiance values recorded by the sensor into vegetation trait estimates using dedicated methods. The methods should ideally be accurate, precise and computationally efficient. Most importantly, they should require minimal calibration since they are supposed to be applied over diverse locations, seasons and conditions *Walthall et al.* (2004).

#### 7.2.1 Biophysical parameter retrieval approaches from optical remote sensing

Many methods have been proposed to retrieve vegetation traits from remote sensing observations *Myneni et al.* (1988), *Goel* (1989), Pinty and Verstraete 1991, XX Kimes et al. 2000; *Baret et al.* (2007), Houborg
and Boegh 2008; Verger et al. (2011), Laurent et al. (2013). They include empirical methods with calibration over experimental data sets. These simple methods are limited by the size and diversity of the calibration data set as well as by the uncertainties attached to the ground measurements. More complex ones based on the use of radiative transfer models have been proposed where no in situ calibration data set is required. Radiative transfer models describe the physical processes involved in the photon transport within vegetation canopies. They simulate the radiation field reflected by the surface for a given observational configuration once the vegetation and the background are known. Retrieving canopy characteristics from the radiation field as sampled by the sensor aboard satellite needs to 'invert' the radiative transfer model, i.e. to estimate some input variables from the measurement of the outputs of the model. This can be done explicitly with a radiative transfer model, or based on datasets where the satellite measurements are matching the canopy traits targeted as measured or estimated from an independent method. Note that the use of ground-based measurement is generally not possible for coarse resolution sensors because of the very demanding spatial sampling required. The retrieval methods can be split into two main approaches depending whether the emphasis is put on the inputs (the canopy biophysical variables driven approach, e.g. vegetation indices or machine learning approach) or the outputs (radiometric data driven approach, e.g. iterative optimization methods that minimize a cost function, and Look-Up-Tables) of the radiative transfer model Figure 7.1.

<figure><img src="image_52.png" /><figcaption>Figure 7.1: The two main approaches used to estimate vegetation traits (exemplified by LAI) from remote sensing data. On the left side, the approach focuses on radiometric data showing the solution search process leading to the estimated LAI value, LAI*. On the right side, the approach focuses on the biophysical variables showing the calibration of the inverse model (top) and the application using the inverse model with its calibrated coefficients (bottom). Δ represents the cost function to be minimized over the biophysical variables (right) or over the radiometric data (left).</figcaption></figure>

In a recent study, Wang et al. (2022) propose to use a Bayesian scheme to estimate some vegetation traits including GAI from observations of a range of satellites, allowing to assess a better prior information regarding the distributions and co-distributions of the radiative transfer model input variables. Domain adaptation techniques based on using GAN (Generative Adversarial networks) could also be investigated to make radiative transfer model simulation better reflect the actual measurements. Finally, hybrid methods Dong et al. (2024) and ensemble products are also interesting solutions to achieve better accuracy on the products Xiao et al. (2014).