# Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype

Jean-Christophe Calvet, Bertrand Bonan, Oscar Rojas-Munoz, Anna

Agusti-Panareda, Patricia de Rosnay, Peter Weston, Philippe Peylin, Cédric

Bacour, Vladislav Bastrikov, Fabienne Maignan, et al.

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

## To cite this version:

Jean-Christophe Calvet, Bertrand Bonan, Oscar Rojas- Munoz, Anna Agusti-Panareda, Patricia de Rosnay, et al. Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype. ECMWF. 2023. meteo-04541420

HAL Id: meteo-04541420

https://meteofrance.hal.science/meteo-04541420

Submitted on 10 Apr 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

CoCO2
Prototype system for a
Copernicus CO2 service

# Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype

Jean-Christophe Calvet, Bertrand Bonan, Oscar Rojas-Munoz, Anna Agusti-Panareda, Patricia de Rosnay, Peter Weston, Philippe Peylin, Cédric Bacour, Vladislav Bastrikov, Fabienne Maignan, Thomas Kaminski, Wolfgang Knorr, Michael Vossbeck, Marko Scholze

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

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

Co-ordinated by
ECMWF

<figure><img src="image_4.png" /></figure>
<figure><img src="image_3.png" /></figure>

CoCO2
Prototype system for a
Copernicus CO2 service

# Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype

Jean-Christophe Calvet, Bertrand Bonan, Oscar Rojas-Munoz, Anna Agusti-Panareda, Patricia de Rosnay, Peter Weston, Philippe Peylin, Cédric Bacour, Vladislav Bastrikov, Fabienne Maignan, Thomas Kaminski, Wolfgang Knorr, Michael Vossbeck, Marko Scholze

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

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

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

Co-ordinated by
ECMWF
<figure><img src="image_5.png" /></figure>

CoCO2
Prototype system for a
Copernicus CO2 service

# D3.4 Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype

**Dissemination Level:** Public

**Author(s):** Jean-Christophe Calvet, Bertrand Bonan, and Oscar Rojas-Munoz (Meteo-France), Anna Agusti-Panareda, Patricia de Rosnay and Peter Weston (ECMWF), Philippe Peylin, Cedric Bacour, Vladislav Bastrikov, and Fabienne Maignan (CEA, LSCE), Thomas Kaminski, Wolfgang Knorr, Michael Vossbeck (iLab), Marko Scholze (ULUND)

**Date:** 18/10/2023

**Version:** 0.1

**Contractual Delivery Date:** 30/06/2023

**Work Package/ Task:** WP3/ T3.3b

**Document Owner:** MF

**Contributors:** Meteo-France, ECMWF, iLab, ULUND

**Status:** Draft/ for Review/ Final

ECMWF

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

Empa

Materials, Science and Technology

ICOS

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

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

CE2

Max Planck Institute
for Biogeochemistry

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

TNO innovation for life

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

LUND UNIVERSITY

VU

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

VEU
UNIVERSITEIT
AMSTERDAM

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

WAGENINGEN

UNIVERSITY & RESEARCH

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

*CICERO

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

cmcc

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

Deutscher Wetterdienst
Wetter und Klima auf einer Hand

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

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

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

FMI

FORTH

iLab

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

MERCATOR
OCEAN
INTERNATIONAL

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

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

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

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

THE CYPRUS
INSTITUTE

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

Deutsches Zentrum
für Luft- und Raumfahrt
<figure><img src="image_7.png" /></figure>

CoCO2
Prototype system for a
Copernicus CO2 service

# D3.4 Demonstrator systems for using remote sensing data (LAI, VOD, SIF) in online global prior fluxes for the CO2MVS prototype

**Dissemination Level:** Public

**Author(s):** Jean-Christophe Calvet, Bertrand Bonan, and Oscar Rojas-Munoz (Meteo-France), Anna Agusti-Panareda, Patricia de Rosnay and Peter Weston (ECMWF), Philippe Peylin, Cedric Bacour, Vladislav Bastrikov, and Fabienne Maignan (CEA, LSCE), Thomas Kaminski, Wolfgang Knorr, Michael Vossbeck (iLab), Marko Scholze (ULUND)

**Date:** 18/10/2023

**Version:** 0.1

**Contractual Delivery Date:** 30/06/2023

**Work Package/ Task:** WP3/ T3.3b

**Document Owner:** MF

**Contributors:** Meteo-France, ECMWF, iLab, ULUND

**Status:** Draft/ for Review/ Final

ECMWF

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

Empa

Materials, Science and Technology

ICOS

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

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

CE2

Max Planck Institute
for Biogeochemistry

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

TNO innovation for life

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

LUND UNIVERSITY

VU

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

VEU
UNIVERSITEIT
AMSTERDAM

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

WAGENINGEN

UNIVERSITY & RESEARCH

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

Barcelona
Computing
Center
Catalunya
Institut de Recerca i Tecnologia de Catalunya

*CICERO

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

cmcc

CITTS

Deutscher Wetterdienst
Wetter- und Klima aus einer Hand

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

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

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

FMI

FORTH

iLab

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

MERCATOR
OCEAN
INTERNATIONAL

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

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

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

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

THE CYPRUS
INSTITUTE

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

Deutsches Zentrum
für Luft- und Raumfahrt
<figure><img src="image_25.png" /></figure>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_43.png" /></figure>
<figure><img src="image_26.png" /></figure>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_44.png" /></figure>
# CoCO2: Prototype system for a Copernicus CO2 service

Coordination and Support Action (CSA)
H2020-IBA-SPACE-CHE2-2019 Copernicus evolution –
Research activities in support of a European operational
monitoring support capacity for fossil CO2 emissions

**Project Coordinator:** Dr Richard Engelen (ECMWF)

**Project Start Date:** 01/01/2021

**Project Duration:** 36 months

Published by the CoCO2 Consortium

**Contact:**

ECMWF, Shinfield Park, Reading, RG2 9AX,
richard.engelen@ecmwf.int

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

The CoCO2 project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 958927.

ECMWF

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

Empa

Material Science and Technology

ICOS

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

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

cea

Max Planck Institute
for Biogeochemistry

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

TNO innovation for life

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

LUND UNIVERSITY

VU

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

VELESTER
UNIVERSITEIT
AMSTERDAM

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

WAGENINGEN

UNIVERSITY & RESEARCH

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

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

Groninger Supercomputing Center
Universiteit Groningen

CICERO

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

cmcc

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

Deutscher Wetterdienst
Wetter und Klima aus einer Hand

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

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

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

FMI

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

FORTH

iLab

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

MERCATOR
OCEAN
INTERNATIONAL

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

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

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

FCNOS

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

THE CYPRUS
INSTITUTE

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

Deutsches Zentrum
für Luft- und Raumfahrt
# CoCO2: Prototype system for a Copernicus CO2 service

Coordination and Support Action (CSA)
H2020-IBA-SPACE-CHE2-2019 Copernicus evolution –
Research activities in support of a European operational
monitoring support capacity for fossil CO2 emissions

**Project Coordinator:** Dr Richard Engelen (ECMWF)

**Project Start Date:** 01/01/2021

**Project Duration:** 36 months

Published by the CoCO2 Consortium

**Contact:**

ECMWF, Shinfield Park, Reading, RG2 9AX,
richard.engelen@ecmwf.int

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

The CoCO2 project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 958927.

ECMWF

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

Empa

Material Science and Technology

ICOS

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

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

cea

Max Planck Institute
for Biogeochemistry

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

TNO

Innovation
for life

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

LUND UNIVERSITY

VU

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

VELE
UNIVERSITE
AMSTERDAM

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

WAGENINGEN

UNIVERSITY & RESEARCH

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

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

Sahel
Supercomputing
Center

CICERO

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

cmcc

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

Deutscher Wetterdienst
Wetter und Klima auf einer Hand

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

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

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

FMI

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

FORTH

iLab

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

MERCATOR
OCEAN
INTERNATIONAL

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

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

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

FCNOS

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

THE CYPRUS
INSTITUTE

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

Deutsches Zentrum
für Luft- und Raumfahrt
Table of Contents

<table><tr><td>1</td><td>Executive Summary</td><td>9</td></tr><tr><td>2</td><td>Introduction</td><td>9</td></tr><tr><td>2.1</td><td>Background</td><td>9</td></tr><tr><td>2.2</td><td>Scope of this deliverable</td><td>10</td></tr><tr><td>2.2.1</td><td>Objectives of this deliverables</td><td>10</td></tr><tr><td>2.2.2</td><td>Work performed in this deliverable</td><td>10</td></tr><tr><td>2.2.3</td><td>Deviations and counter measures</td><td>10</td></tr><tr><td>3</td><td>ISBA</td><td>11</td></tr><tr><td>3.1</td><td>Methods</td><td>11</td></tr><tr><td>3.1.1</td><td>Offline simplified extended Kalman filter (SEKF)</td><td>11</td></tr><tr><td>3.1.2</td><td>Rescaling</td><td>11</td></tr><tr><td>3.1.3</td><td>Experiment setup</td><td>11</td></tr><tr><td>3.2</td><td>Data</td><td>13</td></tr><tr><td>3.3</td><td>Results and evaluation</td><td>13</td></tr><tr><td>3.3.1</td><td>Global GPP (GLOBE experiments)</td><td>13</td></tr><tr><td>3.3.2</td><td>VOD assimilation (CONUS experiments)</td><td>15</td></tr><tr><td>3.3.3</td><td>SIF consistency analysis (EUMED experiments)</td><td>16</td></tr><tr><td>3.4</td><td>Summary</td><td>18</td></tr><tr><td>4</td><td>ECLand</td><td>19</td></tr><tr><td>4.1</td><td>Methods</td><td>19</td></tr><tr><td>4.1.1</td><td>Offline simplified extended Kalman filter (SEKF)</td><td>19</td></tr><tr><td>4.1.2</td><td>Modifications for dynamic LAI analysis</td><td>20</td></tr><tr><td>4.1.3</td><td>Experiment setup</td><td>21</td></tr><tr><td>4.2</td><td>Data</td><td>22</td></tr><tr><td>4.2.1</td><td>Observations</td><td>22</td></tr><tr><td>4.2.2</td><td>Quality control</td><td>24</td></tr><tr><td>4.2.3</td><td>CDF-matching</td><td>25</td></tr><tr><td>4.3</td><td>Results</td><td>27</td></tr><tr><td>4.4</td><td>Evaluation</td><td>29</td></tr><tr><td>4.4.1</td><td>Impact on NWP</td><td>29</td></tr><tr><td>4.4.2</td><td>Impact on carbon cycle</td><td>34</td></tr><tr><td>5</td><td>ORCHIDEE</td><td>36</td></tr><tr><td>5.1</td><td>Methods</td><td>36</td></tr><tr><td>5.1.1</td><td>ORCHIDEE model and optimized parameters</td><td>36</td></tr><tr><td>5.1.2</td><td>Data assimilation system</td><td>37</td></tr><tr><td>5.1.3</td><td>Assimilated data</td><td>38</td></tr><tr><td>5.1.4</td><td>Other data for evaluation and comparison</td><td>39</td></tr></table>
# Table of Contents

<table><tr><td>1</td><td>Executive Summary</td><td>9</td></tr><tr><td>2</td><td>Introduction</td><td>9</td></tr><tr><td>2.1</td><td>Background</td><td>9</td></tr><tr><td>2.2</td><td>Scope of this deliverable</td><td>10</td></tr><tr><td>2.2.1</td><td>Objectives of this deliverables</td><td>10</td></tr><tr><td>2.2.2</td><td>Work performed in this deliverable</td><td>10</td></tr><tr><td>2.2.3</td><td>Deviations and counter measures</td><td>10</td></tr><tr><td>3</td><td>ISBA</td><td>11</td></tr><tr><td>3.1</td><td>Methods</td><td>11</td></tr><tr><td>3.1.1</td><td>Offline simplified extended Kalman filter (SEKF)</td><td>11</td></tr><tr><td>3.1.2</td><td>Rescaling</td><td>11</td></tr><tr><td>3.1.3</td><td>Experiment setup</td><td>11</td></tr><tr><td>3.2</td><td>Data</td><td>13</td></tr><tr><td>3.3</td><td>Results and evaluation</td><td>13</td></tr><tr><td>3.3.1</td><td>Global GPP (GLOBE experiments)</td><td>13</td></tr><tr><td>3.3.2</td><td>VOD assimilation (CONUS experiments)</td><td>15</td></tr><tr><td>3.3.3</td><td>SIF consistency analysis (EUMED experiments)</td><td>16</td></tr><tr><td>3.4</td><td>Summary</td><td>18</td></tr><tr><td>4</td><td>ECLand</td><td>19</td></tr><tr><td>4.1</td><td>Methods</td><td>19</td></tr><tr><td>4.1.1</td><td>Offline simplified extended Kalman filter (SEKF)</td><td>19</td></tr><tr><td>4.1.2</td><td>Modifications for dynamic LAI analysis</td><td>20</td></tr><tr><td>4.1.3</td><td>Experiment setup</td><td>21</td></tr><tr><td>4.2</td><td>Data</td><td>22</td></tr><tr><td>4.2.1</td><td>Observations</td><td>22</td></tr><tr><td>4.2.2</td><td>Quality control</td><td>24</td></tr><tr><td>4.2.3</td><td>CDF-matching</td><td>25</td></tr><tr><td>4.3</td><td>Results</td><td>27</td></tr><tr><td>4.4</td><td>Evaluation</td><td>29</td></tr><tr><td>4.4.1</td><td>Impact on NWP</td><td>29</td></tr><tr><td>4.4.2</td><td>Impact on carbon cycle</td><td>34</td></tr><tr><td>5</td><td>ORCHIDEE</td><td>36</td></tr><tr><td>5.1</td><td>Methods</td><td>36</td></tr><tr><td>5.1.1</td><td>ORCHIDEE model and optimized parameters</td><td>36</td></tr><tr><td>5.1.2</td><td>Data assimilation system</td><td>37</td></tr><tr><td>5.1.3</td><td>Assimilated data</td><td>38</td></tr><tr><td>5.1.4</td><td>Other data for evaluation and comparison</td><td>39</td></tr></table>
<table><tr><td>5.2</td><td>Results</td><td>39</td></tr><tr><td>5.2.1</td><td>Model - data improvement</td><td>39</td></tr><tr><td>5.2.2</td><td>Spatial distribution of GPP</td><td>41</td></tr><tr><td>5.3</td><td>Evaluation using other GPP products</td><td>42</td></tr><tr><td>5.4</td><td>Summary</td><td>44</td></tr><tr><td>6</td><td>SDBM</td><td>45</td></tr><tr><td>6.1</td><td>Methods</td><td>45</td></tr><tr><td>6.2</td><td>Data</td><td>46</td></tr><tr><td>6.3</td><td>Results and evaluation</td><td>47</td></tr><tr><td>6.4</td><td>Summary</td><td>52</td></tr><tr><td>7</td><td>Synthesis of results</td><td>53</td></tr><tr><td>8</td><td>Conclusion</td><td>56</td></tr><tr><td>9</td><td>References</td><td>57</td></tr></table>

Figures

<table><tr><td>Figure 1: Impact of assimilation on ISBA global GPP estimates</td><td>14</td></tr><tr><td>Figure 2: Seasonal rescaling of X-band VOD over CONUS: (top) deciduous forests, (bottom) C3 crops, (left) VOD vs. CGLS LAI, (right) rescaled VOD vs. CGLS LAI. Adapted from Mucia et al. 2022</td><td>15</td></tr><tr><td>Figure 3: Effect of the assimilation of LAI and VOD on ISBA (left) GPP and (right) evapotranspiration (ET) monthly correlation with independent estimates (FLUXCOM-V1 and ALEXI, described in Jung et al. 2020 and Anderson et al. 2011, respectively) over CONUS. VODX10 is VOD in X-band with 10-day undersampling. Adapted from Mucia et al. 2022</td><td>16</td></tr><tr><td>Figure 4: Daily time series of S5-P SIF and OL_EUMED and SEKF_LAI_EUMED GPP (green, dark, and red lines, respectively) over the region of Toulouse in France from 1 May 2018 to 31 December 2019</td><td>16</td></tr><tr><td>Figure 5: Correlation maps between the TROPOMI SIF observations and the GPP analysis for the 2018–2020 period (top, left) and the difference map (top, right), where red means an improvement after assimilation, and comparison (bottom) with the same maps obtained with GOME-2 for the 2007-2015 period (adapted from Leroux et al. 2018)</td><td>17</td></tr><tr><td>Figure 6: Correlation maps between the TROPOMI SIF observations and the LAI analysis for the 2018–2020 period (left column) and the difference map (right column, where red means an improvement after assimilation)</td><td>18</td></tr><tr><td>Figure 7: Schematic of the ECMWF LDAS within the IFS (de Rosnay et al., 2022)</td><td>20</td></tr><tr><td>Figure 8: Correlation between collocated SMOS L-band VOD and AMSR2 X-band VOD observations, for 1st January 2016</td><td>23</td></tr><tr><td>Figure 9: Map of daily coverage of AMSR2 6.9GHz VOD observations on 1st July 2016</td><td>24</td></tr><tr><td>Figure 10: Maps of AMSR2 observations contaminated by RFI for the 6.9GHz channel (top) and 7.3GHz channel (bottom), data is from 1st July 2016</td><td>25</td></tr><tr><td>Figure 11: Correlations between VOD observations and CONFESS LAI for high vegetation types for July 2018. The different panels show: upper left SMOS L-band VOD; upper right rescaled SMOS-L band VOD; middle left AMSR2 C-band VOD; middle right rescaled AMSR2 C-band VOD; lower right AMSR2 X-band VOD; lower left rescaled AMSR2 X-band VOD</td><td>26</td></tr></table>
<table><tr><td>5.2</td><td>Results</td><td>39</td></tr><tr><td>5.2.1</td><td>Model – data improvement</td><td>39</td></tr><tr><td>5.2.2</td><td>Spatial distribution of GPP</td><td>41</td></tr><tr><td>5.3</td><td>Evaluation using other GPP products</td><td>42</td></tr><tr><td>5.4</td><td>Summary</td><td>44</td></tr><tr><td>6</td><td>SDBM</td><td>45</td></tr><tr><td>6.1</td><td>Methods</td><td>45</td></tr><tr><td>6.2</td><td>Data</td><td>46</td></tr><tr><td>6.3</td><td>Results and evaluation</td><td>47</td></tr><tr><td>6.4</td><td>Summary</td><td>52</td></tr><tr><td>7</td><td>Synthesis of results</td><td>53</td></tr><tr><td>8</td><td>Conclusion</td><td>56</td></tr><tr><td>9</td><td>References</td><td>57</td></tr></table>

Figures

<table><tr><td>Figure 1: Impact of assimilation on ISBA global GPP estimates</td><td>14</td></tr><tr><td>Figure 2: Seasonal rescaling of X-band VOD over CONUS: (top) deciduous forests, (bottom) C3 crops, (left) VOD vs. CGLS LAI, (right) rescaled VOD vs. CGLS LAI. Adapted from Mucia et al. 2022</td><td>15</td></tr><tr><td>Figure 3: Effect of the assimilation of LAI and VOD on ISBA (left) GPP and (right) evapotranspiration (ET) monthly correlation with independent estimates (FLUXCOM-V1 and ALEXI, described in Jung et al. 2020 and Anderson et al. 2011, respectively) over CONUS. VODX10 is VOD in X-band with 10-day undersampling. Adapted from Mucia et al. 2022</td><td>16</td></tr><tr><td>Figure 4: Daily time series of S5-P SIF and OL_EUMED and SEKF_LAI_EUMED GPP (green, dark, and red lines, respectively) over the region of Toulouse in France from 1 May 2018 to 31 December 2019</td><td>16</td></tr><tr><td>Figure 5: Correlation maps between the TROPOMI SIF observations and the GPP analysis for the 2018–2020 period (top, left) and the difference map (top, right), where red means an improvement after assimilation, and comparison (bottom) with the same maps obtained with GOME-2 for the 2007-2015 period (adapted from Leroux et al. 2018)</td><td>17</td></tr><tr><td>Figure 6: Correlation maps between the TROPOMI SIF observations and the LAI analysis for the 2018–2020 period (left column) and the difference map (right column, where red means an improvement after assimilation)</td><td>18</td></tr><tr><td>Figure 7: Schematic of the ECMWF LDAS within the IFS (de Rosnay et al., 2022)</td><td>20</td></tr><tr><td>Figure 8: Correlation between collocated SMOS L-band VOD and AMSR2 X-band VOD observations, for 1st January 2016</td><td>23</td></tr><tr><td>Figure 9: Map of daily coverage of AMSR2 6.9GHz VOD observations on 1st July 2016</td><td>24</td></tr><tr><td>Figure 10: Maps of AMSR2 observations contaminated by RFI for the 6.9GHz channel (top) and 7.3GHz channel (bottom), data is from 1st July 2016</td><td>25</td></tr><tr><td>Figure 11: Correlations between VOD observations and CONFESS LAI for high vegetation types for July 2018. The different panels show: upper left SMOS L-band VOD; upper right rescaled SMOS-L band VOD; middle left AMSR2 C-band VOD; middle right rescaled AMSR2 C-band VOD; lower right AMSR2 X-band VOD; lower left rescaled AMSR2 X-band VOD</td><td>26</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>Figure 12: Monthly mean analysis increments of leaf area index for high vegetation types for January 2018 (upper left), April 2018 (upper right), July 2018 (lower left), October 2018 (lower right)</td>
      <td>27</td>
    </tr>
    <tr>
      <td>Figure 13: As Figure 12 but for leaf area index for low vegetation types</td>
      <td>28</td>
    </tr>
    <tr>
      <td>Figure 14: Monthly mean analysis increments of leaf area index for low vegetation types for July 2021 from LDAS-X-VOD (left) and LDAS-L-VOD (right)</td>
      <td>28</td>
    </tr>
    <tr>
      <td>Figure 15: Global change in RMSE of 2 metre temperature forecasts from T+48 to T+120 hours between IFS-X-VOD and IFS-CTRL. Statistics cover January 2018 to December 2021</td>
      <td>29</td>
    </tr>
    <tr>
      <td>Figure 16: As figure 15 but for 850hPa temperature forecasts</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Figure 17: As figure 15 but for 850hPa relative humidity forecasts</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Figure 18: Mean temperature (left) and relative humidity (right) forecast errors at 850hPa (upper) and mean 2 metre temperature forecast errors (lower) for IFS-X-VOD (black) and IFS-CTRL (red). Statistics cover January 2018 to December 2021</td>
      <td>31</td>
    </tr>
    <tr>
      <td>Figure 19: Change in RMSE for 2 metre temperature forecasts between IFS-X-VOD (left), IFS-C-VOD (centre), IFS-L-VOD (right) and IFS-CTRL. Statistics cover January 2018 to December 2021</td>
      <td>32</td>
    </tr>
    <tr>
      <td>Figure 20: Change in RMSE of 2 metre temperature forecasts between IFS-L-VOD and IFS-CTRL for December-January-February (top), March-April-May (2<sup>nd</sup>), June-July-August (3<sup>rd</sup>), September-October-November (bottom). Statistics cover January 2018 to December 2021</td>
      <td>33</td>
    </tr>
    <tr>
      <td>Figure 21: Differences in accumulated CO<sub>2</sub> gross primary production absolute bias between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</td>
      <td>34</td>
    </tr>
    <tr>
      <td>Figure 22: Differences in accumulated CO<sub>2</sub> gross primary production correlation between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</td>
      <td>34</td>
    </tr>
    <tr>
      <td>Figure 23: Comparison of the mean seasonal cycles of GPP for each PFT, between in situ data (black) and ORCHIDEE simulations (prior in blue and posterior in orange). Plain lines correspond to the mean values over the sites considered, and the shaded areas represent the corresponding standard deviation. Units are in gC/m<sup>2</sup>/day. The RMSD of the fit before and after assimilation are provided for each vegetation type, as well as the reduction in RMSD (in %). The PFTs indicated in green correspond to those for which FLUXSAT GPP data are used instead of in situ flux observations.</td>
      <td>40</td>
    </tr>
    <tr>
      <td>Figure 24: Same legend as figure 23 for SIF data (over 2018-2020)</td>
      <td>40</td>
    </tr>
    <tr>
      <td>Figure 25: Reduction in RMSD (in %) between ORCHIDEE simulations and SIF (left) and GPP (right) data following the data assimilation.</td>
      <td>41</td>
    </tr>
    <tr>
      <td>Figure 26: Yearly mean map for the simulations performed with the ORCHIDEE land surface model prior (left) and posterior (right) to data assimilation over the period 2018 (from May) - 2020. The global minimum, maximum, and mean values are provided (gC/m<sup>2</sup>/day), as well as the global budget (in PgC/m<sup>2</sup>/day).</td>
      <td>41</td>
    </tr>
    <tr>
      <td>Figure 27: Comparison of the average GPP time series for the simulations performed with ORCHIDEE (prior and optimized), the mean and standard deviation of the TRENDY simulations (black, grey shaded area) and the FluxCOM-V2 and FluxSat data products, for the northern hemisphere temperate latitudinal band (60°N / 30°N).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 28: Similar caption to Figure 27 for the northern tropics latitudinal band (30°N / 0°).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 29: Similar caption to Figure 27 for the southern tropics latitudinal band (0° / -30°N).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 30: Maps of annual GPP from SDBM for the year 2021 using the ERA5-derived water stress factor (top) and the BETHY-derived water stress factor (middle) compared to the annual GPP from FLUXCOM v2.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>Figure 31: Time series of globally aggregated, monthly modelled GPP from SDBM (blue line) using the ERA5-derived water stress factor (left) and the BETHY-derived water stress factor (right) compared against the globally aggregated FLUXCOM v2 GPP product (red line).</td>
      <td>49</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>Figure 12: Monthly mean analysis increments of leaf area index for high vegetation types for January 2018 (upper left), April 2018 (upper right), July 2018 (lower left), October 2018 (lower right)</td>
      <td>27</td>
    </tr>
    <tr>
      <td>Figure 13: As Figure 12 but for leaf area index for low vegetation types</td>
      <td>28</td>
    </tr>
    <tr>
      <td>Figure 14: Monthly mean analysis increments of leaf area index for low vegetation types for July 2021 from LDAS-X-VOD (left) and LDAS-L-VOD (right)</td>
      <td>28</td>
    </tr>
    <tr>
      <td>Figure 15: Global change in RMSE of 2 metre temperature forecasts from T+48 to T+120 hours between IFS-X-VOD and IFS-CTRL. Statistics cover January 2018 to December 2021</td>
      <td>29</td>
    </tr>
    <tr>
      <td>Figure 16: As figure 15 but for 850hPa temperature forecasts</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Figure 17: As figure 15 but for 850hPa relative humidity forecasts</td>
      <td>30</td>
    </tr>
    <tr>
      <td>Figure 18: Mean temperature (left) and relative humidity (right) forecast errors at 850hPa (upper) and mean 2 metre temperature forecast errors (lower) for IFS-X-VOD (black) and IFS-CTRL (red). Statistics cover January 2018 to December 2021</td>
      <td>31</td>
    </tr>
    <tr>
      <td>Figure 19: Change in RMSE for 2 metre temperature forecasts between IFS-X-VOD (left), IFS-C-VOD (centre), IFS-L-VOD (right) and IFS-CTRL. Statistics cover January 2018 to December 2021</td>
      <td>32</td>
    </tr>
    <tr>
      <td>Figure 20: Change in RMSE of 2 metre temperature forecasts between IFS-L-VOD and IFS-CTRL for December-January-February (top), March-April-May (2<sup>nd</sup>), June-July-August (3<sup>rd</sup>), September-October-November (bottom). Statistics cover January 2018 to December 2021</td>
      <td>33</td>
    </tr>
    <tr>
      <td>Figure 21: Differences in accumulated CO<sub>2</sub> gross primary production absolute bias between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</td>
      <td>34</td>
    </tr>
    <tr>
      <td>Figure 22: Differences in accumulated CO<sub>2</sub> gross primary production correlation between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</td>
      <td>34</td>
    </tr>
    <tr>
      <td>Figure 23: Comparison of the mean seasonal cycles of GPP for each PFT, between in situ data (black) and ORCHIDEE simulations (prior in blue and posterior in orange). Plain lines correspond to the mean values over the sites considered, and the shaded areas represent the corresponding standard deviation. Units are in gC/m<sup>2</sup>/day. The RMSD of the fit before and after assimilation are provided for each vegetation type, as well as the reduction in RMSD (in %). The PFTs indicated in green correspond to those for which FLUXSAT GPP data are used instead of in situ flux observations.</td>
      <td>40</td>
    </tr>
    <tr>
      <td>Figure 24: Same legend as figure 23 for SIF data (over 2018-2020)</td>
      <td>40</td>
    </tr>
    <tr>
      <td>Figure 25: Reduction in RMSD (in %) between ORCHIDEE simulations and SIF (left) and GPP (right) data following the data assimilation.</td>
      <td>41</td>
    </tr>
    <tr>
      <td>Figure 26: Yearly mean map for the simulations performed with the ORCHIDEE land surface model prior (left) and posterior (right) to data assimilation over the period 2018 (from May) - 2020. The global minimum, maximum, and mean values are provided (gC/m<sup>2</sup>/day), as well as the global budget (in PgC/m<sup>2</sup>/day).</td>
      <td>41</td>
    </tr>
    <tr>
      <td>Figure 27: Comparison of the average GPP time series for the simulations performed with ORCHIDEE (prior and optimized), the mean and standard deviation of the TRENDY simulations (black, grey shaded area) and the FluxCOM-V2 and FluxSat data products, for the northern hemisphere temperate latitudinal band (60°N / 30°N).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 28: Similar caption to Figure 27 for the northern tropics latitudinal band (30°N / 0°).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 29: Similar caption to Figure 27 for the southern tropics latitudinal band (0° / -30°N).</td>
      <td>43</td>
    </tr>
    <tr>
      <td>Figure 30: Maps of annual GPP from SDBM for the year 2021 using the ERA5-derived water stress factor (top) and the BETHY-derived water stress factor (middle) compared to the annual GPP from FLUXCOM v2.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>Figure 31: Time series of globally aggregated, monthly modelled GPP from SDBM (blue line) using the ERA5-derived water stress factor (left) and the BETHY-derived water stress factor (right) compared against the globally aggregated FLUXCOM v2 GPP product (red line).</td>
      <td>49</td>
    </tr>
  </tbody>
</table>
# 1 Executive Summary

The purpose of this work is to ingest satellite observations of leaf area index (LAI), microwave vegetation optical depth (VOD), and solar induced fluorescence (SIF) to constrain the terrestrial biogenic carbon uptake. Numerical experiments performed by 4 modelling approaches (ISBA, ECLand, ORCHIDEE, SDBM) are compared. ISBA and ECLand use an LDAS approach to perform a sequential assimilation of satellite observations to analyse model state variables, ORCHIDEE uses satellite data to optimize the value of model parameters, SDBM is directly driven by satellite FAPAR. The VOD assimilation is evaluated using ISBA and ECLand. The consistency of SIF with LAI and gross primary production (GPP) is evaluated using ISBA. Global FLUXCOM-V2 GPP estimates are used to benchmark models.

# 2 Introduction

## 2.1 Background

Following the CO2 Human Emissions (CHE) project, and the Observation-based system for monitoring and verification of greenhouse gases (VERIFY) project, CoCO2 aims to create a prototype anthropogenic CO2 emission Monitoring and Verification Support (i.e., the CO2MVS) systems at global, regional, and local scales in the framework of low-carbon economy and implemented commitments under the Paris Agreement (Delbeke et al., 2019) and the European Green Deal (EU Commission, 2021). This effort will support countries with an independent assessment of their emissions and reduction targets. Moreover, the resulting information, being consistent and reliable, will provide support to policy- and decision-making processes both at national and European level.

The global modelling and data assimilation work package 3 (WP3) of CoCO2 is focused on the development of the CO2MVS at a global scale, building on existing capacity, and will deliver a system that can provide information of CO2 emissions and other surface exchanges in near-real time, as well as in past reanalysis mode. The European Centre for Medium-range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) is at the core of the CO2MVS with planned developments in several areas, such as data assimilation, tracer advection, representation of anthropogenic impact on land-surface exchanges, simulations of biogenic fluxes, among others.

Biogenic fluxes play a fundamental role in the carbon cycle with a key component centred on the capacity of plants to absorb carbon via photosynthesis. This process is also crucial for the land-surface water and energy cycles. These three cycles, water/energy/carbon, are coupled and interact on time-scales ranging from minutes to centuries. Among different aspects of the processes involved, Land Use and Land Cover (LULC) is extremely relevant in the estimation of biogenic carbon fluxes (Quaife et al., 2008) and their changes are extensively studied in the context of climate change. Moreover, the errors found in the models' representation of LULC effects on the lower troposphere have also been shown to limit the progress in weather and climate predictability (e.g. Guo et al., 2011; Orth et al., 2016).

In the current ECMWF land surface model (ECLand) the LAI is based on a 2000-2008 climatology derived from Moderate Resolution Imaging Spectroradiometer (MODIS) collection 5 (MOD15A2) data and rescaled using a static LAI dataset to guarantee neutral impact on the ECMWF model (Boussetta et al., 2013). The interannual variability of LAI is not represented. Therefore, the representation of vegetation in ECLand does not fully benefit from the developments of remote sensing vegetation data sets during the past 20 years and developing techniques to integrate such data into ECLand is needed.
# 1 Executive Summary

The purpose of this work is to ingest satellite observations of leaf area index (LAI), microwave vegetation optical depth (VOD), and solar induced fluorescence (SIF) to constrain the terrestrial biogenic carbon uptake. Numerical experiments performed by 4 modelling approaches (ISBA, ECLand, ORCHIDEE, SDBM) are compared. ISBA and ECLand use an LDAS approach to perform a sequential assimilation of satellite observations to analyse model state variables, ORCHIDEE uses satellite data to optimize the value of model parameters, SDBM is directly driven by satellite FAPAR. The VOD assimilation is evaluated using ISBA and ECLand. The consistency of SIF with LAI and gross primary production (GPP) is evaluated using ISBA. Global FLUXCOM-V2 GPP estimates are used to benchmark models.

# 2 Introduction

## 2.1 Background

Following the CO2 Human Emissions (CHE) project, and the Observation-based system for monitoring and verification of greenhouse gases (VERIFY) project, CoCO2 aims to create a prototype anthropogenic CO2 emission Monitoring and Verification Support (i.e., the CO2MVS) systems at global, regional, and local scales in the framework of low-carbon economy and implemented commitments under the Paris Agreement (Delbeke et al., 2019) and the European Green Deal (EU Commission, 2021). This effort will support countries with an independent assessment of their emissions and reduction targets. Moreover, the resulting information, being consistent and reliable, will provide support to policy- and decision-making processes both at national and European level.

The global modelling and data assimilation work package 3 (WP3) of CoCO2 is focused on the development of the CO2MVS at a global scale, building on existing capacity, and will deliver a system that can provide information of CO2 emissions and other surface exchanges in near-real time, as well as in past reanalysis mode. The European Centre for Medium-range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) is at the core of the CO2MVS with planned developments in several areas, such as data assimilation, tracer advection, representation of anthropogenic impact on land-surface exchanges, simulations of biogenic fluxes, among others.

Biogenic fluxes play a fundamental role in the carbon cycle with a key component centred on the capacity of plants to absorb carbon via photosynthesis. This process is also crucial for the land-surface water and energy cycles. These three cycles, water/energy/carbon, are coupled and interact on time-scales ranging from minutes to centuries. Among different aspects of the processes involved, Land Use and Land Cover (LULC) is extremely relevant in the estimation of biogenic carbon fluxes (Quaife et al., 2008) and their changes are extensively studied in the context of climate change. Moreover, the errors found in the models' representation of LULC effects on the lower troposphere have also been shown to limit the progress in weather and climate predictability (e.g. Guo et al., 2011; Orth et al., 2016).

In the current ECMWF land surface model (ECLand) the LAI is based on a 2000-2008 climatology derived from Moderate Resolution Imaging Spectroradiometer (MODIS) collection 5 (MOD15A2) data and rescaled using a static LAI dataset to guarantee neutral impact on the ECMWF model (Boussetta et al., 2013). The interannual variability of LAI is not represented. Therefore, the representation of vegetation in ECLand does not fully benefit from the developments of remote sensing vegetation data sets during the past 20 years and developing techniques to integrate such data into ECLand is needed.
## 2.2 Scope of this deliverable

### 2.2.1 Objectives of this deliverables

This deliverable is integrated in Task 3.3 aiming to improve global simulations of biogenic fluxes and to integrate land-surface remote sensing observations and ancillary data related to vegetation carbon exchanges in the IFS prototype.

Key objectives of this document are to

* assess to what extent vegetation satellite products such as LAI, microwave Vegetation Optical Depth (VOD), and Solar-Induced Fluorescence (SIF) could be integrated into ECLand

* demonstrate the use of Land Data Assimilation Systems (LDAS) at ECMWF and Meteo-France

* present a potential observation operator for VOD

* assimilate VOD to analyse LAI

* intercompare Gross Primary Production (GPP) estimates from the different modelling approaches, with and without the assimilation of vegetation products: ECLand, ISBA, ORCHIDEE, and SDBM, developed by ECMWF, Meteo-France, LSCE, and ULund/iLab, respectively.

### 2.2.2 Work performed in this deliverable

This deliverable was produced through the collaboration between all authors from the different institutions. This included ECLand offline and online model simulations, ISBA simulations within the SURFEX modelling platform, ORCHIDEE simulations, SDBM.

ECMWF, Meteo-France and LSCE have assimilated vegetation variables in their data assimilation system and started exploring the use of new satellite data. ULund and iLab forced SDBM with satellite FAPAR vegetation data.

Each partner did their own data processing, analysis, evaluation, in their own data processing environment, and wrote a chapter of this document describing the methods, the data used, and the results.

Each chapter describes the progress achieved by a given partner within its own modelling system, with a self-assessment section (pros and cons), and yearly global GPP estimates derived from their model.

Chapters 3, 4, 5, and 6, correspond to ISBA (Meteo-France), ECLand (ECMWF), ORCHIDEE (LSCE), and SDBM (ULund and iLab). The FLUXCOM-V2 dataset produced by WP5 was used in the evaluation. In Chapter 7, pros and cons of data and modelling approaches are presented, together with a comparison of simulated global yearly GPP with FLUXCOM-V2.

### 2.2.3 Deviations and counter measures

There were no major deviations from the work plan.
## 2.2 Scope of this deliverable

### 2.2.1 Objectives of this deliverables

This deliverable is integrated in Task 3.3 aiming to improve global simulations of biogenic fluxes and to integrate land-surface remote sensing observations and ancillary data related to vegetation carbon exchanges in the IFS prototype.

Key objectives of this document are to

* assess to what extent vegetation satellite products such as LAI, microwave Vegetation Optical Depth (VOD), and Solar-Induced Fluorescence (SIF) could be integrated into ECLand

* demonstrate the use of Land Data Assimilation Systems (LDAS) at ECMWF and Meteo-France

* present a potential observation operator for VOD

* assimilate VOD to analyse LAI

* intercompare Gross Primary Production (GPP) estimates from the different modelling approaches, with and without the assimilation of vegetation products: ECLand, ISBA, ORCHIDEE, and SDBM, developed by ECMWF, Meteo-France, LSCE, and ULund/iLab, respectively.

### 2.2.2 Work performed in this deliverable

This deliverable was produced through the collaboration between all authors from the different institutions. This included ECLand offline and online model simulations, ISBA simulations within the SURFEX modelling platform, ORCHIDEE simulations, SDBM.

ECMWF, Meteo-France and LSCE have assimilated vegetation variables in their data assimilation system and started exploring the use of new satellite data. ULund and iLab forced SDBM with satellite FAPAR vegetation data.

Each partner did their own data processing, analysis, evaluation, in their own data processing environment, and wrote a chapter of this document describing the methods, the data used, and the results.

Each chapter describes the progress achieved by a given partner within its own modelling system, with a self-assessment section (pros and cons), and yearly global GPP estimates derived from their model.

Chapters 3, 4, 5, and 6, correspond to ISBA (Meteo-France), ECLand (ECMWF), ORCHIDEE (LSCE), and SDBM (ULund and iLab). The FLUXCOM-V2 dataset produced by WP5 was used in the evaluation. In Chapter 7, pros and cons of data and modelling approaches are presented, together with a comparison of simulated global yearly GPP with FLUXCOM-V2.

### 2.2.3 Deviations and counter measures

There were no major deviations from the work plan.
# 3 ISBA

At the Meteorological Research Division of Meteo-France, CNRM (Centre National de Recherches Météorologiques), LDAS-Monde (Albergel et al., 2017) has been developed as an offline LDAS capable of sequentially and simultaneously assimilating LAI and surface soil moisture (SSM) into the ISBA (Interactions between Soil-Biosphere-Atmosphere) land surface model (Noilhan and Mahfouf, 1996; Calvet et al., 1998, 2004; Gibelin et al., 2006; Barbu et al., 2014). LDAS-Monde has the ability to integrate satellite data into ISBA at global and regional scales and can be used to monitor and predict the state of land surface variables (Albergel et al., 2017, 2019, 2020; Tall et al., 2019; Bonan et al., 2020; Mucia et al., 2022). Knowing the strong positive influence of LAI assimilation on the simulation, and also recognising the weakness of assimilating these observations only every 10 d, Mucia et al. (2022) investigated the assimilation of seasonally linearly rescaled VOD as a proxy for LAI over contiguous United States (CONUS). In addition, the consistency between simulated and analysed LAI and GPP was assessed over the Euro-Mediterranean (EUMED) region using S5-P TROPOMI SIF data.

## 3.1 Methods

The LDAS-Monde tool uses a version of ISBA capable of simulating plant growth together with the diffusion of water and heat into the soil. A parsimonious approach is used to represent plant growth and leaf phenology (Calvet et al. 1998, Gibelin et al. 2006). The simulated leaf phenology is entirely driven by photosynthesis and depends on only two static parameters: the maximum leaf span time, which ranges from 150 d for low vegetation to 365 d for evergreen trees, and the minimum LAI, which ranges from 0.3 m²m⁻² for low vegetation to 1.0 m²m⁻² for evergreen trees. The resulting simulated LAI is flexible and can be analysed at any time when an observation is available. This allows for a sequential assimilation of LAI (Barbu et al. 2014). In addition to LAI analysis, LAI assimilation provides an analysis of root-zone soil moisture. Note that LDAS-Monde does not guarantee conservation of mass.

### 3.1.1 Offline simplified extended Kalman filter (SEKF)

Within the SURFEX (Surface Externalisée Version 8.1) modelling platform (Masson et al., 2013) developed by CNRM (http://www.umr-cnrm.fr/surfex), a package allows the offline assimilation of satellite-derived products into the ISBA land surface model using the LDAS-Monde data assimilation system. In the LDAS-Monde system (Albergel et al., 2017), observed satellite LAI and SSM data are assimilated into the ISBA LSM using a simplified extended Kalman filter (SEKF) data assimilation technique (Mahfouf et al., 2009). The SEKF technique uses finite differences to compute the flow dependency between the assimilated observations, SSM and LAI, and the control variables (soil moisture from soil layers 2 to 8 (0.01 to 1.00 m) and LAI). The eight control variables are updated directly by the observed variables according to their sensitivity as given by the SEKF Jacobian matrices (Barbu et al., 2014). Other model variables are updated indirectly through feedbacks and other biophysical processes as related to the control variables.

### 3.1.2 Rescaling

A seasonal linear rescaling technique is used to (1) fit assimilated SSM observations to the soil characteristics used in the ISBA model, (2) use VOD as a proxy for LAI, (3) use THEIA AVHRR LAI as a proxy for CGLS LAI. Assimilated LAI is not rescaled to fit the model LAI.

### 3.1.3 Experiment setup

Meteo-France carried out 8 numerical experiments, listed in Table 1.
# 3 ISBA

At the Meteorological Research Division of Meteo-France, CNRM (Centre National de Recherches Météorologiques), LDAS-Monde (Albergel et al., 2017) has been developed as an offline LDAS capable of sequentially and simultaneously assimilating LAI and surface soil moisture (SSM) into the ISBA (Interactions between Soil-Biosphere-Atmosphere) land surface model (Noilhan and Mahfouf, 1996; Calvet et al., 1998, 2004; Gibelin et al., 2006; Barbu et al., 2014). LDAS-Monde has the ability to integrate satellite data into ISBA at global and regional scales and can be used to monitor and predict the state of land surface variables (Albergel et al., 2017, 2019, 2020; Tall et al., 2019; Bonan et al., 2020; Mucia et al., 2022). Knowing the strong positive influence of LAI assimilation on the simulation, and also recognising the weakness of assimilating these observations only every 10 d, Mucia et al. (2022) investigated the assimilation of seasonally linearly rescaled VOD as a proxy for LAI over contiguous United States (CONUS). In addition, the consistency between simulated and analysed LAI and GPP was assessed over the Euro-Mediterranean (EUMED) region using S5-P TROPOMI SIF data.

## 3.1 Methods

The LDAS-Monde tool uses a version of ISBA capable of simulating plant growth together with the diffusion of water and heat into the soil. A parsimonious approach is used to represent plant growth and leaf phenology (Calvet et al. 1998, Gibelin et al. 2006). The simulated leaf phenology is entirely driven by photosynthesis and depends on only two static parameters: the maximum leaf span time, which ranges from 150 d for low vegetation to 365 d for evergreen trees, and the minimum LAI, which ranges from 0.3 m²m⁻² for low vegetation to 1.0 m²m⁻² for evergreen trees. The resulting simulated LAI is flexible and can be analysed at any time when an observation is available. This allows for a sequential assimilation of LAI (Barbu et al. 2014). In addition to LAI analysis, LAI assimilation provides an analysis of root-zone soil moisture. Note that LDAS-Monde does not guarantee conservation of mass.

### 3.1.1 Offline simplified extended Kalman filter (SEKF)

Within the SURFEX (Surface Externalisée Version 8.1) modelling platform (Masson et al., 2013) developed by CNRM (http://www.umr-cnrm.fr/surfex), a package allows the offline assimilation of satellite-derived products into the ISBA land surface model using the LDAS-Monde data assimilation system. In the LDAS-Monde system (Albergel et al., 2017), observed satellite LAI and SSM data are assimilated into the ISBA LSM using a simplified extended Kalman filter (SEKF) data assimilation technique (Mahfouf et al., 2009). The SEKF technique uses finite differences to compute the flow dependency between the assimilated observations, SSM and LAI, and the control variables (soil moisture from soil layers 2 to 8 (0.01 to 1.00 m) and LAI). The eight control variables are updated directly by the observed variables according to their sensitivity as given by the SEKF Jacobian matrices (Barbu et al., 2014). Other model variables are updated indirectly through feedbacks and other biophysical processes as related to the control variables.

### 3.1.2 Rescaling

A seasonal linear rescaling technique is used to (1) fit assimilated SSM observations to the soil characteristics used in the ISBA model, (2) use VOD as a proxy for LAI, (3) use THEIA AVHRR LAI as a proxy for CGLS LAI. Assimilated LAI is not rescaled to fit the model LAI.

### 3.1.3 Experiment setup

Meteo-France carried out 8 numerical experiments, listed in Table 1.
<table><caption>Table 1: Numerical experiments using ISBA and the offline LDAS-Monde tool*</caption><thead><tr><th>Experiment (time period) (area)</th><th>Assimilated satellite observations</th><th>Model equivalent</th><th>Model control variables</th><th>Atmospheric forcing</th></tr></thead><tbody><tr><td>OL_GLOBE (2002-2019) (global)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_LAI_SSM_GLOBE (2002-2019) (global)</td><td>GEOV2 CGLS LAI and ESA-CCI v6.1 SSM (rescaled)</td><td>LAI, WG2 (0.01-0.04 m)</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>OL_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_LAI_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>GEOV2 CGLS LAI</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_VOD_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>VODCA X-band VOD (rescaled)</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>OL_EUMED (2018-2020) (30-72 N, 25W-45E)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>HRES re-interpolated at 0.10°</td></tr><tr><td>SEKF_LAI_EUMED (2018-2019) (30-72 N, 25W-45E)</td><td>GEOV1 CGLS LAI</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>HRES re-interpolated at 0.10°</td></tr><tr><td>SEKF_LAI_EUMED (2020) (30-72 N, 25W-45E)</td><td>THEIA AVHRR LAI (rescaled to match GEOV1 CGLS LAI)</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>HRES re-interpolated at 0.10°</td></tr></tbody></table>

Different periods and regions were considered depending on the availability of atmospheric forcing data and assimilated satellite products. Two global experiments (GLOBE) were performed for the period 2002-2019, with and without LAI and SSM assimilation in the ISBA model: SEKF_LAI_SSM_GLOBE and OL_GLOBE, respectively. Three experiments were performed over CONUS from 2003 to 2018, with and without the assimilation of LAI, and with the assimilation of VOD: SEKF_LAI_CONUS, OL_CONUS, and SEKF_VOD_CONUS, respectively. Two experiments were performed over EUMED from 2018 to 2020, with and without assimilation of LAI: SEKF_LAI_EUMED and OL_EUMED, respectively. In 2020, the rescaled THEIA AVHRR LAI was assimilated instead of the CGLS LAI to mitigate the end of the PROBA-V satellite time series. All experiments used ERA5 as the atmospheric forcing, except the EUMED experiments, which are based on the higher spatial resolution HRES forcing produced by the IFS. The CONUS experiments were designed to evaluate the

<footnote>* OL, SEKF, CONUS, EUMED, LAI, SSM, VOD, CGLS, VODCA, WG2, WG8, ERA5, HRES stand for Open-Loop, Simplified Extended Kalman Filter, Continental United States, Euro-Mediterranean, Leaf Area Index, Surface Soil Moisture, microwave Vegetation optical Depth, Copernicus Global Land Service, Vegetation Optical Depth Climate Archive (Moesinger et al. 2020), volumetric soil moisture for model soil layers 2 (0.01-0.04 m) and 8 (0.80-1.00 m), ECMWF reanalysis (Hersbach et al., 2020), ECMWF high resolution operational high resolution weather analysis and forecasts, respectively. n/a stands for "not applicable"</footnote>
<table><caption>Table 1: Numerical experiments using ISBA and the offline LDAS-Monde tool*</caption><thead><tr><th>Experiment (time period) (area)</th><th>Assimilated satellite observations</th><th>Model equivalent</th><th>Model control variables</th><th>Atmospheric forcing</th></tr></thead><tbody><tr><td>OL_GLOBE (2002-2019) (global)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_LAI_SSM_GLOBE (2002-2019) (global)</td><td>GEOV2 CGLS LAI and ESA-CCI v6.1 SSM (rescaled)</td><td>LAI, WG2 (0.01-0.04 m)</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>OL_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_LAI_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>GEOV2 CGLS LAI</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>SEKF_VOD_CONUS (2003-2018) (20-55 N, 130W-60W)</td><td>VODCA X-band VOD (rescaled)</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>ERA5 re-interpolated at 0.25°</td></tr><tr><td>OL_EUMED (2018-2020) (30-72 N, 25W-45E)</td><td>n/a</td><td>n/a</td><td>n/a</td><td>HRES re-interpolated at 0.10°</td></tr><tr><td>SEKF_LAI_EUMED (2018-2019) (30-72 N, 25W-45E)</td><td>GEOV1 CGLS LAI</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>HRES re-interpolated at 0.10°</td></tr><tr><td>SEKF_LAI_EUMED (2020) (30-72 N, 25W-45E)</td><td>THEIA AVHRR LAI (rescaled to match GEOV1 CGLS LAI)</td><td>LAI</td><td>LAI, WG2 to WG8 (0.01-1 m)</td><td>HRES re-interpolated at 0.10°</td></tr></tbody></table>

Different periods and regions were considered depending on the availability of atmospheric forcing data and assimilated satellite products. Two global experiments (GLOBE) were performed for the period 2002-2019, with and without LAI and SSM assimilation in the ISBA model: SEKF_LAI_SSM_GLOBE and OL_GLOBE, respectively. Three experiments were performed over CONUS from 2003 to 2018, with and without the assimilation of LAI, and with the assimilation of VOD: SEKF_LAI_CONUS, OL_CONUS, and SEKF_VOD_CONUS, respectively. Two experiments were performed over EUMED from 2018 to 2020, with and without assimilation of LAI: SEKF_LAI_EUMED and OL_EUMED, respectively. In 2020, the rescaled THEIA AVHRR LAI was assimilated instead of the CGLS LAI to mitigate the end of the PROBA-V satellite time series. All experiments used ERA5 as the atmospheric forcing, except the EUMED experiments, which are based on the higher spatial resolution HRES forcing produced by the IFS. The CONUS experiments were designed to evaluate the

<footnote>* OL, SEKF, CONUS, EUMED, LAI, SSM, VOD, CGLS, VODCA, WG2, WG8, ERA5, HRES stand for Open-Loop, Simplified Extended Kalman Filter, Continental United States, Euro-Mediterranean, Leaf Area Index, Surface Soil Moisture, microwave Vegetation optical Depth, Copernicus Global Land Service, Vegetation Optical Depth Climate Archive (Moesinger et al. 2020), volumetric soil moisture for model soil layers 2 (0.01-0.04 m) and 8 (0.80-1.00 m), ECMWF reanalysis (Hersbach et al., 2020), ECMWF high resolution operational high resolution weather analysis and forecasts, respectively. n/a stands for "not applicable"</footnote>
assimilation of VOD. The EUMED experiments were designed to investigate the relationship between SIF, GPP and LAI.

## 3.2 Data

Several data sources are used to perform the numerical experiments listed in Table 1:

* Land cover is provided by ECOCLIMAP within SURFEX (http://www.umr-cnrm.fr/surfex/spip.php?article136),

* Atmospheric forcing: the meteorological variables of air temperature, wind speed, air specific humidity, atmospheric pressure, shortwave and longwave downwelling radiation, and liquid and solid precipitation are ingested into the ISBA model; this study uses atmospheric reanalyses from ECMWF's ERA5 (Hersbach et al. 2020) for GLOBE and CONUS experiments and ECMWF's operational high resolution weather analysis and forecasts (HRES) for EUMED experiments; the data are re-interpolated over the model grid (0.25° x 0.25° for GLOBE and CONUS, 0.10° x 0.10° for EUMED),

* Assimilated satellite observations: LAI derived from the SPOT-VGT and PROBA-V satellites by the Copernicus Global Land Service (CGLS) (https://land.copernicus.eu/global/products/lai), ESA-CCI SSM from combined active and passive microwave sensors (https://climate.esa.int/en/odp/#/project/soil-moisture), VOD from Moesinger et al. 2019 (https://zenodo.org/record/2575599).

* SIF satellite observations: this product is derived from Sentinel 5-P TROPOMI observations (https://s5p-troposif.noveltis.fr/data-access/) in the 743-758 nm window. We re-gridded this product on a regular grid at a spatial resolution of 0.5 x 0.5°.

## 3.3 Results and evaluation

### 3.3.1 Global GPP (GLOBE experiments)

Global GPP estimates are shown in Table 2 for the OL_GLOBE and SEKF_LAI_SSM_GLOBE experiments and for FLUXCOM-V2 from 2002 to 2019.

<table><caption>Table 2: Global GPP estimations (in PgC/year) from the ISBA land surface model from 2002 to 2019, and from FLUXCOM-V2.</caption><thead><tr><td>Year</td><td>OL_GLOBE</td><td>SEKF_LAI_SSM_GLOBE</td><td>FLUXCOM-V2</td></tr></thead><tbody><tr><td>2002</td><td>116.1</td><td>111.0</td><td>130.6</td></tr><tr><td>2003</td><td>118.3</td><td>113.1</td><td>130.6</td></tr><tr><td>2004</td><td>119.8</td><td>114.8</td><td>132.5</td></tr><tr><td>2005</td><td>118.5</td><td>113.9</td><td>132.2</td></tr><tr><td>2006</td><td>121.2</td><td>116.9</td><td>133.5</td></tr><tr><td>2007</td><td>119.9</td><td>115.3</td><td>133.3</td></tr><tr><td>2008</td><td>120.9</td><td>116.0</td><td>133.3</td></tr><tr><td>2009</td><td>121.1</td><td>116.2</td><td>133.9</td></tr><tr><td>2010</td><td>124.6</td><td>119.2</td><td>134.3</td></tr><tr><td>2011</td><td>124.2</td><td>119.3</td><td>135.3</td></tr><tr><td>2012</td><td>122.7</td><td>117.9</td><td>133.9</td></tr><tr><td>2013</td><td>123.1</td><td>118.3</td><td>134.7</td></tr><tr><td>2014</td><td>125.2</td><td>119.8</td><td>136.2</td></tr></tbody></table>
assimilation of VOD. The EUMED experiments were designed to investigate the relationship between SIF, GPP and LAI.

## 3.2 Data

Several data sources are used to perform the numerical experiments listed in Table 1:

* Land cover is provided by ECOCLIMAP within SURFEX (http://www.umr-cnrm.fr/surfex/spip.php?article136),

* Atmospheric forcing: the meteorological variables of air temperature, wind speed, air specific humidity, atmospheric pressure, shortwave and longwave downwelling radiation, and liquid and solid precipitation are ingested into the ISBA model; this study uses atmospheric reanalyses from ECMWF's ERA5 (Hersbach et al. 2020) for GLOBE and CONUS experiments and ECMWF's operational high resolution weather analysis and forecasts (HRES) for EUMED experiments; the data are re-interpolated over the model grid (0.25° x 0.25° for GLOBE and CONUS, 0.10° x 0.10° for EUMED),

* Assimilated satellite observations: LAI derived from the SPOT-VGT and PROBA-V satellites by the Copernicus Global Land Service (CGLS) (https://land.copernicus.eu/global/products/lai), ESA-CCI SSM from combined active and passive microwave sensors (https://climate.esa.int/en/odp/#/project/soil-moisture), VOD from Moesinger et al. 2019 (https://zenodo.org/record/2575599).

* SIF satellite observations: this product is derived from Sentinel 5-P TROPOMI observations (https://s5p-troposif.noveltis.fr/data-access/) in the 743-758 nm window. We re-gridded this product on a regular grid at a spatial resolution of 0.5 x 0.5°.

## 3.3 Results and evaluation

### 3.3.1 Global GPP (GLOBE experiments)

Global GPP estimates are shown in Table 2 for the OL_GLOBE and SEKF_LAI_SSM_GLOBE experiments and for FLUXCOM-V2 from 2002 to 2019.

<table><caption>Table 2: Global GPP estimations (in PgC/year) from the ISBA land surface model from 2002 to 2019, and from FLUXCOM-V2.</caption><thead><tr><td>Year</td><td>OL_GLOBE</td><td>SEKF_LAI_SSM_GLOBE</td><td>FLUXCOM-V2</td></tr></thead><tbody><tr><td>2002</td><td>116.1</td><td>111.0</td><td>130.6</td></tr><tr><td>2003</td><td>118.3</td><td>113.1</td><td>130.6</td></tr><tr><td>2004</td><td>119.8</td><td>114.8</td><td>132.5</td></tr><tr><td>2005</td><td>118.5</td><td>113.9</td><td>132.2</td></tr><tr><td>2006</td><td>121.2</td><td>116.9</td><td>133.5</td></tr><tr><td>2007</td><td>119.9</td><td>115.3</td><td>133.3</td></tr><tr><td>2008</td><td>120.9</td><td>116.0</td><td>133.3</td></tr><tr><td>2009</td><td>121.1</td><td>116.2</td><td>133.9</td></tr><tr><td>2010</td><td>124.6</td><td>119.2</td><td>134.3</td></tr><tr><td>2011</td><td>124.2</td><td>119.3</td><td>135.3</td></tr><tr><td>2012</td><td>122.7</td><td>117.9</td><td>133.9</td></tr><tr><td>2013</td><td>123.1</td><td>118.3</td><td>134.7</td></tr><tr><td>2014</td><td>125.2</td><td>119.8</td><td>136.2</td></tr></tbody></table>
<table><tr><td>2015</td><td>123.2</td><td>118.1</td><td>134.6</td></tr><tr><td>2016</td><td>126.8</td><td>120.6</td><td>136.2</td></tr><tr><td>2017</td><td>128.2</td><td>123.0</td><td>136.8</td></tr><tr><td>2018</td><td>127.0</td><td>122.5</td><td>137.1</td></tr><tr><td>2019</td><td>125.5</td><td>120.9</td><td>137.3</td></tr><tr><td>Mean values</td><td>122.6</td><td>117.6</td><td>134.2</td></tr></table>

All GPP estimates show a rapid increase of GPP during this period, of about 0.6 PgC/yr in both OL and SEKF ISBA simulations, and 0.4 PgC/yr on average for FLUXCOM-V2.

The assimilation of LAI and SSM satellite products tends to reduce the GPP simulated by ISBA, on average by 5 PgC/year. The mean GPP value of the analysis, from 2002 to 2019, is 117.6 PgC/year. This is about 12% smaller than the FLUXCOM-V2 mean GPP for the same period. The mean GPP estimate from OL_GLOBE is about 9% lower than FLUXCOM-V2.

The spatial distribution of the OL_GLOBE mean GPP bias is shown in Figure 1.b, and it can be observed that the areas with the largest negative bias are found in Europe, Brazil (Pantanal, Nordeste), Pakistan, northern India (Pendjab, Uttar Pradesh), southwestern China (Sichuan), forested areas of the northwestern USA, and Russian boreal forests.

<figure><img src="image_66.png" /><figcaption>Figure 1: Impact of assimilation on ISBA global GPP estimates.</figcaption></figure>

Figure 1.c shows that assimilation reduces the bias mainly in areas where the bias is positive: rainforests, eastern USA, southern Brazil, Paraguay, Uruguay.

Figure 1.d shows that the ISBA and FLUXCOM-V2 GPP correlate well ($R > 0.8$) at high latitudes. Poor correlation values are observed over areas with sparse vegetation (e.g., Kazakhstan).

Figure 1.e shows that assimilation generally improves the temporal correlation between ISBA and FLUXCOM-V2 GPP, except in some areas with sparse vegetation. In such areas, LAI observations may be lower than the minimum LAI value used in the model for sparse vegetation ($0.3 m^2 m^{-2}$) and the assimilation reduces the seasonal and inter-annual variability of the simulated LAI.
<table><tr><td>2015</td><td>123.2</td><td>118.1</td><td>134.6</td></tr><tr><td>2016</td><td>126.8</td><td>120.6</td><td>136.2</td></tr><tr><td>2017</td><td>128.2</td><td>123.0</td><td>136.8</td></tr><tr><td>2018</td><td>127.0</td><td>122.5</td><td>137.1</td></tr><tr><td>2019</td><td>125.5</td><td>120.9</td><td>137.3</td></tr><tr><td>Mean values</td><td>122.6</td><td>117.6</td><td>134.2</td></tr></table>

All GPP estimates show a rapid increase of GPP during this period, of about 0.6 PgC/yr in both OL and SEKF ISBA simulations, and 0.4 PgC/yr on average for FLUXCOM-V2.

The assimilation of LAI and SSM satellite products tends to reduce the GPP simulated by ISBA, on average by 5 PgC/year. The mean GPP value of the analysis, from 2002 to 2019, is 117.6 PgC/year. This is about 12% smaller than the FLUXCOM-V2 mean GPP for the same period. The mean GPP estimate from OL_GLOBE is about 9% lower than FLUXCOM-V2.

The spatial distribution of the OL_GLOBE mean GPP bias is shown in Figure 1.b, and it can be observed that the areas with the largest negative bias are found in Europe, Brazil (Pantanal, Nordeste), Pakistan, northern India (Pendjab, Uttar Pradesh), southwestern China (Sichuan), forested areas of the northwestern USA, and Russian boreal forests.

<figure><img src="image_67.png" /><figcaption>Figure 1: Impact of assimilation on ISBA global GPP estimates.</figcaption></figure>

Figure 1.c shows that assimilation reduces the bias mainly in areas where the bias is positive: rainforests, eastern USA, southern Brazil, Paraguay, Uruguay.

Figure 1.d shows that the ISBA and FLUXCOM-V2 GPP correlate well ($R > 0.8$) at high latitudes. Poor correlation values are observed over areas with sparse vegetation (e.g., Kazakhstan).

Figure 1.e shows that assimilation generally improves the temporal correlation between ISBA and FLUXCOM-V2 GPP, except in some areas with sparse vegetation. In such areas, LAI observations may be lower than the minimum LAI value used in the model for sparse vegetation ($0.3 m^2 m^{-2}$) and the assimilation reduces the seasonal and inter-annual variability of the simulated LAI.
3.3.2 VOD assimilation (CONUS experiments)

Satellite LAI observations are available every 10 d at best and can suffer from missing data for months over regions and seasons with heavy cloud cover, such as winter or monsoon conditions. In this context, Mucia et al. (2022) investigated the assimilation of low-frequency passive microwave vegetation optical depth (VOD), which is available in almost all weather conditions, as a proxy for LAI.

A generally positive relationship was found between observations of LAI and VOD in X-band (VODX). This relationship still contains variability strongly related to the dominant vegetation. These results are consistent with previous work comparing vegetation cover and VOD. Coniferous forests consistently had the weakest correlation, while C3 and C4 vegetation typically had the best. Seasonality also dramatically changed the LAI-VODX relationship, with winter score values typically the lowest and summer and fall correlations typically the strongest. Rescaled VODX was much more strongly correlated with LAI observations (Figure 2). In Figure 2, linear regression and correlation scores have been plotted over the data. Logarithmic transformations to the VODX data were also applied but with no significant increase in regression correlation or regression shape.

<figure><img src="image_67.png" /><figcaption>Figure 2: Seasonal rescaling of X-band VOD over CONUS: (top) deciduous forests, (bottom) C3 crops, (left) VOD vs. CGLS LAI, (right) rescaled VOD vs. CGLS LAI. Adapted from Mucia et al. 2022.</figcaption></figure>

The assessment over CONUS showed improved representation of evapotranspiration and GPP compared to independent observations over some months of VOD assimilation. The improvements in evapotranspiration and GPP correlations seen by assimilating rescaled VODX instead of LAI are almost entirely due to the more frequent observations of VOD. This is demonstrated by the fact that the SEKF VODX10 experiment, which assimilates adjusted VODX at the same frequency as LAI observations, performs significantly worse than the natural VODX observation frequency (Figure 3).
3.3.2 VOD assimilation (CONUS experiments)

Satellite LAI observations are available every 10 d at best and can suffer from missing data for months over regions and seasons with heavy cloud cover, such as winter or monsoon conditions. In this context, Mucia et al. (2022) investigated the assimilation of low-frequency passive microwave vegetation optical depth (VOD), which is available in almost all weather conditions, as a proxy for LAI.

A generally positive relationship was found between observations of LAI and VOD in X-band (VODX). This relationship still contains variability strongly related to the dominant vegetation. These results are consistent with previous work comparing vegetation cover and VOD. Coniferous forests consistently had the weakest correlation, while C3 and C4 vegetation typically had the best. Seasonality also dramatically changed the LAI-VODX relationship, with winter score values typically the lowest and summer and fall correlations typically the strongest. Rescaled VODX was much more strongly correlated with LAI observations (Figure 2). In Figure 2, linear regression and correlation scores have been plotted over the data. Logarithmic transformations to the VODX data were also applied but with no significant increase in regression correlation or regression shape.

<figure><img src="image_68.png" /><figcaption>Figure 2: Seasonal rescaling of X-band VOD over CONUS: (top) deciduous forests, (bottom) C3 crops, (left) VOD vs. CGLS LAI, (right) rescaled VOD vs. CGLS LAI. Adapted from Mucia et al. 2022.</figcaption></figure>

The assessment over CONUS showed improved representation of evapotranspiration and GPP compared to independent observations over some months of VOD assimilation. The improvements in evapotranspiration and GPP correlations seen by assimilating rescaled VODX instead of LAI are almost entirely due to the more frequent observations of VOD. This is demonstrated by the fact that the SEKF VODX10 experiment, which assimilates adjusted VODX at the same frequency as LAI observations, performs significantly worse than the natural VODX observation frequency (Figure 3).
<figure><img src="image_68.png" /><figcaption>**Figure 3:** Effect of the assimilation of LAI and VOD on ISBA (left) GPP and (right) evapotranspiration (ET) monthly correlation with independent estimates (FLUXCOM-V1 and ALEXI, described in Jung et al. 2020 and Anderson et al. 2011, respectively) over CONUS. VODX10 is VOD in X-band with 10-day undersampling. Adapted from Mucia et al. 2022.</figcaption></figure>

### 3.3.3 SIF consistency analysis (EUMED experiments)

The OL_EUMED and SEKF_LAI_EUMED GPP and LAI simulations were compared with the TROPOSIF product derived from the S5-P TROPOMI observations. Figure 4 shows that SIF and GPP compare well during the growing season but a mismatch appears during the senescence.

<figure><img src="image_69.png" /><figcaption>**Figure 4:** Daily time series of S5-P SIF and OL_EUMED and SEKF_LAI_EUMED GPP (green, dark, and red lines, respectively) over the region of Toulouse in France from 1 May 2018 to 31 December 2019.</figcaption></figure>

Figure 5 shows a map of the temporal correlation between TROPOMI SIF and SEKF_LAI_EUMED GPP for the period 2018-2020, together with the correlation difference between SEKF_LAI_EUMED and OL_EUMED.
<figure><img src="image_69.png" /></figure>

**Figure 3:** Effect of the assimilation of LAI and VOD on ISBA (left) GPP and (right) evapotranspiration (ET) monthly correlation with independent estimates (FLUXCOM-V1 and ALEXI, described in Jung et al. 2020 and Anderson et al. 2011, respectively) over CONUS. VODX10 is VOD in X-band with 10-day undersampling. Adapted from Mucia et al. 2022.

### 3.3.3 SIF consistency analysis (EUMED experiments)

The OL_EUMED and SEKF_LAI_EUMED GPP and LAI simulations were compared with the TROPOSIF product derived from the S5-P TROPOMI observations. Figure 4 shows that SIF and GPP compare well during the growing season but a mismatch appears during the senescence.

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

**Figure 4:** Daily time series of S5-P SIF and OL_EUMED and SEKF_LAI_EUMED GPP (green, dark, and red lines, respectively) over the region of Toulouse in France from 1 May 2018 to 31 December 2019.

Figure 5 shows a map of the temporal correlation between TROPOMI SIF and SEKF_LAI_EUMED GPP for the period 2018-2020, together with the correlation difference between SEKF_LAI_EUMED and OL_EUMED.
<figure><img src="image_70.png" /><figcaption>Figure 5: Correlation maps between the TROPOMI SIF observations and the GPP analysis for the 2018–2020 period (top, left) and the difference map (top, right), where red means an improvement after assimilation, and comparison (bottom) with the same maps obtained with GOME-2 for the 2007-2015 period (adapted from Leroux et al. 2018).</figcaption></figure>

A comparison with the previous results of Leroux et al. (2018) is presented. In the latter study, a GOME-2 SIF product is used for the period 2007-2015. Overall, the correlation between SIF and GPP is lower for TROPOMI SIF than for GOME-2. The increase in correlation due to assimilation is also often smaller (e.g. France, Ukraine). On the other hand, Figure 6 shows that the temporal correlation between TROPOMI SIF and SEKF_LAI_EUMED LAI is very good at mid-latitudes (from France to Ukraine) and a large increase in correlation due to assimilation is observed. A different result is observed at high latitudes, where the correlation of TROPOMI SIF with LAI is not as good as that with GPP.
<figure><img src="image_71.png" /><figcaption>Figure 5: Correlation maps between the TROPOMI SIF observations and the GPP analysis for the 2018–2020 period (top, left) and the difference map (top, right), where red means an improvement after assimilation, and comparison (bottom) with the same maps obtained with GOME-2 for the 2007-2015 period (adapted from Leroux et al. 2018).</figcaption></figure>

A comparison with the previous results of Leroux et al. (2018) is presented. In the latter study, a GOME-2 SIF product is used for the period 2007-2015. Overall, the correlation between SIF and GPP is lower for TROPOMI SIF than for GOME-2. The increase in correlation due to assimilation is also often smaller (e.g. France, Ukraine). On the other hand, Figure 6 shows that the temporal correlation between TROPOMI SIF and SEKF_LAI_EUMED LAI is very good at mid-latitudes (from France to Ukraine) and a large increase in correlation due to assimilation is observed. A different result is observed at high latitudes, where the correlation of TROPOMI SIF with LAI is not as good as that with GPP.
<figure><img src="image_71.png" /><figcaption>Figure 6: Correlation maps between the TROPOMI SIF observations and the LAI analysis for the 2018–2020 period (left column) and the difference map (right column, where red means an improvement after assimilation).</figcaption></figure>

3.4 Summary

The ISBA land surface model takes a parsimonious approach to simulating plant growth and phenology, using only a few parameters (e.g. only 2 for phenology). At the same time, rather complex soil-vegetation processes are represented (e.g. heat and water diffusion into the soil, drought-avoiding and drought-tolerant plant response to soil moisture). The resulting simulated LAI is flexible and can be analysed at any time when an observation is available. This allows sequential assimilation of LAI. In addition to LAI analysis, LAI assimilation provides a direct analysis of root zone soil moisture. Note that LDAS-Monde does not guarantee mass conservation.

This work shows that ISBA and FLUXCOM-V2 GPP correlate well ($R > 0.8$) at high latitudes. Poor correlation values are observed over areas with sparse vegetation (e.g. Kazakhstan). Assimilation generally improves the temporal correlation between ISBA and FLUXCOM-V2 GPP, except for sparse vegetation areas.

The assimilation of VOD improves the representation of evapotranspiration and GPP, and the improvements are almost entirely due to the more frequent observations of VOD compared to LAI.

The temporal correlation between TROPOMI SIF and the analysed LAI is very good at mid-latitudes (e.g. from France to Ukraine) and a strong increase of the correlation due to the assimilation is observed. SIF and GPP are in good agreement during the growing season, but a discrepancy appears during the senescence. A different result is observed at high latitudes, where the correlation of TROPOMI SIF with LAI is not as good as that with GPP. This result could be due to the rather poor representation of LAI of conifers in the ISBA model version we used, since coniferous forests are more common at high latitudes.
<figure><img src="image_72.png" /><figcaption>**Figure 6:** Correlation maps between the TROPOMI SIF observations and the LAI analysis for the 2018–2020 period (left column) and the difference map (right column, where red means an improvement after assimilation).</figcaption></figure>

## 3.4 Summary

The ISBA land surface model takes a parsimonious approach to simulating plant growth and phenology, using only a few parameters (e.g. only 2 for phenology). At the same time, rather complex soil-vegetation processes are represented (e.g. heat and water diffusion into the soil, drought-avoiding and drought-tolerant plant response to soil moisture). The resulting simulated LAI is flexible and can be analysed at any time when an observation is available. This allows sequential assimilation of LAI. In addition to LAI analysis, LAI assimilation provides a direct analysis of root zone soil moisture. Note that LDAS-Monde does not guarantee mass conservation.

This work shows that ISBA and FLUXCOM-V2 GPP correlate well ($R > 0.8$) at high latitudes. Poor correlation values are observed over areas with sparse vegetation (e.g. Kazakhstan). Assimilation generally improves the temporal correlation between ISBA and FLUXCOM-V2 GPP, except for sparse vegetation areas.

The assimilation of VOD improves the representation of evapotranspiration and GPP, and the improvements are almost entirely due to the more frequent observations of VOD compared to LAI.

The temporal correlation between TROPOMI SIF and the analysed LAI is very good at mid-latitudes (e.g. from France to Ukraine) and a strong increase of the correlation due to the assimilation is observed. SIF and GPP are in good agreement during the growing season, but a discrepancy appears during the senescence. A different result is observed at high latitudes, where the correlation of TROPOMI SIF with LAI is not as good as that with GPP. This result could be due to the rather poor representation of LAI of conifers in the ISBA model version we used, since coniferous forests are more common at high latitudes.
# 4 ECLand

Currently, vegetation parameters such as LAI are specified as seasonally-varying climatological monthly mean maps in the ECMWF Numerical Weather Prediction (NWP) system. This climatology is in the process of being updated using the latest CGLS LAI dataset from the CONFESS project (Boussetta & Balsamo, 2021) and has been shown to have a significant impact on the quality of NWP forecasts (Boussetta, personal communication). This change is aimed for inclusion in Integrated Forecasting System (IFS) cycle 49r1 to go operational in 2024. Despite the updates to the source of vegetation data contributing to the monthly climatology, the fundamental methodology will not change in 49r1.

The significant NWP impacts of the planned upgrade suggest that changes to the methodology of how the vegetation information is passed to the NWP system may have a substantial additional impact. Sub-monthly changes are accounted for by doing a temporal interpolation of the monthly fields based on how many days through the month the forecast is started from. However, one of the main weaknesses of the current approach is the inability to respond to inter-annual variability in vegetation because, although the climatology varies seasonally, it has no inter-annual variability. In fact, the current climatology used in the IFS is built using a weighted mean of the previous 20 years of available LAI observations. Inter-annual differences in vegetation can be large because of meteorological events such as droughts, above average rainfall and variations in 2 metre temperature. The effects of climate change have exaggerated these variations in recent years. The differences can be particularly large in the transition seasons (spring and autumn) when the vegetation characteristics undergo the largest changes during the annual cycle, especially in the mid-latitudes where the seasonal conditions vary the most.

In this study a daily LAI analysis is produced using the ECMWF land data assimilation system (LDAS) and this analysis is passed to the IFS NWP system instead of the monthly climatology. Using this approach has the advantage of more dynamic updates to the vegetation parameters and should be better able to respond to inter-annual variations than using the fixed monthly climatology. This follows similar work of Mucia et al (2022), Kumar et al (2020), although in this study the LAI of high vegetation types (e.g. evergreen and deciduous trees, mixed and interrupted forest) and low vegetation types (e.g. crops and mixed farming, grass, shrubs etc.) are considered separately and the LAI analyses are produced over the entire global domain.

## 4.1 Methods

### 4.1.1 Offline simplified extended Kalman filter (SEKF)

The ECMWF LDAS within the IFS is made up of several different component parts, see Figure 7. This system is weakly coupled to the atmospheric assimilation system whereby the land and atmospheric assimilation systems run separately to produce their respective analyses. The coupling is achieved within the forecast model where the land and atmospheric components of the Earth system interact and influence one another (de Rosnay et al., 2022). The forecast model provides the background for the next assimilation cycle and so on.
# 4 ECLand

Currently, vegetation parameters such as LAI are specified as seasonally-varying climatological monthly mean maps in the ECMWF Numerical Weather Prediction (NWP) system. This climatology is in the process of being updated using the latest CGLS LAI dataset from the CONFESS project (Boussetta & Balsamo, 2021) and has been shown to have a significant impact on the quality of NWP forecasts (Boussetta, personal communication). This change is aimed for inclusion in Integrated Forecasting System (IFS) cycle 49r1 to go operational in 2024. Despite the updates to the source of vegetation data contributing to the monthly climatology, the fundamental methodology will not change in 49r1.

The significant NWP impacts of the planned upgrade suggest that changes to the methodology of how the vegetation information is passed to the NWP system may have a substantial additional impact. Sub-monthly changes are accounted for by doing a temporal interpolation of the monthly fields based on how many days through the month the forecast is started from. However, one of the main weaknesses of the current approach is the inability to respond to inter-annual variability in vegetation because, although the climatology varies seasonally, it has no inter-annual variability. In fact, the current climatology used in the IFS is built using a weighted mean of the previous 20 years of available LAI observations. Inter-annual differences in vegetation can be large because of meteorological events such as droughts, above average rainfall and variations in 2 metre temperature. The effects of climate change have exaggerated these variations in recent years. The differences can be particularly large in the transition seasons (spring and autumn) when the vegetation characteristics undergo the largest changes during the annual cycle, especially in the mid-latitudes where the seasonal conditions vary the most.

In this study a daily LAI analysis is produced using the ECMWF land data assimilation system (LDAS) and this analysis is passed to the IFS NWP system instead of the monthly climatology. Using this approach has the advantage of more dynamic updates to the vegetation parameters and should be better able to respond to inter-annual variations than using the fixed monthly climatology. This follows similar work of Mucia et al (2022), Kumar et al (2020), although in this study the LAI of high vegetation types (e.g. evergreen and deciduous trees, mixed and interrupted forest) and low vegetation types (e.g. crops and mixed farming, grass, shrubs etc.) are considered separately and the LAI analyses are produced over the entire global domain.

## 4.1 Methods

### 4.1.1 Offline simplified extended Kalman filter (SEKF)

The ECMWF LDAS within the IFS is made up of several different component parts, see Figure 7. This system is weakly coupled to the atmospheric assimilation system whereby the land and atmospheric assimilation systems run separately to produce their respective analyses. The coupling is achieved within the forecast model where the land and atmospheric components of the Earth system interact and influence one another (de Rosnay et al., 2022). The forecast model provides the background for the next assimilation cycle and so on.
<figure><img src="image_72.png" /><figcaption>Figure 7: Schematic of the ECMWF LDAS within the IFS (de Rosnay et al., 2022)</figcaption></figure>

In this study an offline implementation of the LDAS is used (Rodríguez-Fernández et al., 2019). This is a simplified, uncoupled system which reads in the atmospheric forcing as well as the pseudo 2 metre temperature and 2 metre relative humidity observations from an external source - here the ERA5 reanalysis (Hersbach et al., 2020). The SEKF is run and the land surface analysis produced is passed to an offline version of the land surface forecast model H-TESSEL (Balsamo et al., 2009). This has the disadvantage that there is only a one-way coupling between the atmosphere and the land but allows for developments to be tested at low computing cost for multi-year periods before being implemented in the weakly coupled IFS LDAS.

The offline SEKF analyses 3 layers of soil moisture (0-7cm, 7-28cm, 28-100cm) using the aforementioned pseudo-observations of 2 metre temperature and relative humidity in combination with surface soil moisture observations from the Advanced Scatterometer (ASCAT) instrument onboard the MetOp series of satellites. The analysis is calculated using the Kalman filter equation:

$$ \mathbf{x}_a = \mathbf{x}_b + \mathbf{K} (\mathbf{y} - \mathbf{H}(\mathbf{x}_b)) \quad (4.1) $$

$$ \mathbf{K} = \mathbf{B}\mathbf{H}^{\mathrm{T}} (\mathbf{R} + \mathbf{H}\mathbf{B}\mathbf{H}^{\mathrm{T}})^{-1} \quad (4.2) $$

where $\mathbf{x}_a$ is the analysis, $\mathbf{x}_b$ is the background, $\mathbf{y}$ are the observations, $\mathbf{H}$ is the observation operator, $\mathbf{R}$ contains the Jacobians from the ensemble of data assimilations (EDA) to link the model variables to the observed variables, $\mathbf{B}$ is the background error covariance matrix and $\mathbf{R}$ is the observation error covariance matrix.

### 4.1.2 Modifications for dynamic LAI analysis

A number of modifications were made to the existing offline SEKF to produce a dynamic daily LAI analysis. LAI is not a prognostic variable. A daily LAI analysis is produced using the offline LDAS by assimilating the VOD observations (after Cumulative Distribution Function (CDF)-matching to LAI observations) and using the climatological LAI as the background. This analysed LAI is used in place of the climatology in IFS experiments. The control vector was extended to include high and low vegetation LAI in addition to the existing 3 layers of soil moisture. The background vector was extended to include the high and low vegetation LAI climatology as a first-guess. The background errors for LAI are assumed to be 0.2m²/m² and uncorrelated with the background errors for soil moisture in all layers. This is a simplistic but pragmatic assumption given that in the current system the background errors of soil moisture
<figure><img src="image_73.png" /><figcaption>Figure 7: Schematic of the ECMWF LDAS within the IFS (de Rosnay et al., 2022)</figcaption></figure>

In this study an offline implementation of the LDAS is used (Rodríguez-Fernández et al., 2019). This is a simplified, uncoupled system which reads in the atmospheric forcing as well as the pseudo 2 metre temperature and 2 metre relative humidity observations from an external source - here the ERA5 reanalysis (Hersbach et al., 2020). The SEKF is run and the land surface analysis produced is passed to an offline version of the land surface forecast model H-TESSEL (Balsamo et al., 2009). This has the disadvantage that there is only a one-way coupling between the atmosphere and the land but allows for developments to be tested at low computing cost for multi-year periods before being implemented in the weakly coupled IFS LDAS.

The offline SEKF analyses 3 layers of soil moisture (0-7cm, 7-28cm, 28-100cm) using the aforementioned pseudo-observations of 2 metre temperature and relative humidity in combination with surface soil moisture observations from the Advanced Scatterometer (ASCAT) instrument onboard the MetOp series of satellites. The analysis is calculated using the Kalman filter equation:

$$ \mathbf{x}_a = \mathbf{x}_b + \mathbf{K} (\mathbf{y} - \mathbf{H}(\mathbf{x}_b)) \quad (4.1) $$

$$ \mathbf{K} = \mathbf{B}\mathbf{H}^{\mathrm{T}}(\mathbf{R} + \mathbf{H}\mathbf{B}\mathbf{H}^{\mathrm{T}})^{-1} \quad (4.2) $$

where $\mathbf{x}_a$ is the analysis, $\mathbf{x}_b$ is the background, $\mathbf{y}$ are the observations, $\mathbf{H}$ is the observation operator, $\mathbf{R}$ contains the Jacobians from the ensemble of data assimilations (EDA) to link the model variables to the observed variables, $\mathbf{B}$ is the background error covariance matrix and $\mathbf{R}$ is the observation error covariance matrix.

### 4.1.2 Modifications for dynamic LAI analysis

A number of modifications were made to the existing offline SEKF to produce a dynamic daily LAI analysis. LAI is not a prognostic variable. A daily LAI analysis is produced using the offline LDAS by assimilating the VOD observations (after Cumulative Distribution Function (CDF)-matching to LAI observations) and using the climatological LAI as the background. This analysed LAI is used in place of the climatology in IFS experiments. The control vector was extended to include high and low vegetation LAI in addition to the existing 3 layers of soil moisture. The background vector was extended to include the high and low vegetation LAI climatology as a first-guess. The background errors for LAI are assumed to be 0.2m²/m² and uncorrelated with the background errors for soil moisture in all layers. This is a simplistic but pragmatic assumption given that in the current system the background errors of soil moisture
in the different layers are assumed uncorrelated when significant correlations are expected to exist.

The VOD observations (see section 4.2 for details) are ingested into the system once per day. In the absence of a physically-based observation operator to transform the background LAI into simulated VOD a rescaling from VOD to LAI is required. To do this a CDF-matching approach is used which is similar to the approach used to transform ASCAT surface soil moisture to the model soil moisture as part of the existing ASCAT soil moisture assimilation (Scipal et al., 2008). The CDF-matching is done using 2 years of training data (2016-17) to match the VOD observations to the climatology of the time-varying CONFESS LAI data.

The VOD observations are collocated with the model grid points and the CDF-matching parameters are calculated separately for each model grid point to account for geographical variations, and separately for each month of the year using a three-month rolling average to account for seasonal variations whilst preventing large jumps in the conversions for consecutive months.

For ASCAT surface soil moisture the CDF is matched both using the mean and standard deviation but for VOD/LAI matching this led to poor results due to a lack of month-to-month variability in the CONFESS LAI dataset, so the VOD/LAI CDF-matching just matches the means as shown in the equation below:

$$VOD' = (E(LAI) / E(VOD)) * VOD \quad (4.3)$$

where VOD' is the rescaled VOD, E(LAI) is the mean of the LAI, E(VOD) is the mean of the VOD and VOD represents the original VOD observations.

Once the CDF-matching is performed the transformed VOD observations are in the same units and have the same dynamic range as the model LAI so the observation operator (H) is simply the identity matrix. For simplicity, no correlations are assumed between the LAI and the other soil moisture observations are assumed and, vice versa, no correlations are assumed between the soil moisture and the transformed VOD observations. The specification of non-zero correlations in these blocks of the Jacobian matrix could be considered in future work.

### 4.1.3 Experiment setup

The experiments are run for 4 years between 2018 and 2021. This period is motivated by the availability of the observations (not to overlap with the CDF-matching training period 2016-17) and to sample multiple annual cycles of vegetation. The experiments are run at a resolution of T<sub>co</sub>319 (~35km) which is significantly coarser than the operational resolution of T<sub>co</sub>1279 (~9km). This resolution is purely chosen for pragmatic reasons to enable the running of multiple experimental configurations for such a long period of time, and is very close to the standard ECMWF testing resolution of T<sub>co</sub>399. Three separate experiment configurations are run in addition to the control run as detailed in Table 3.

<table><caption>Table 3: Summary of offline ECLand LDAS experiments</caption><thead><tr><td>Experiment ID</td><td>VOD observations assimilated</td><td>Period</td></tr></thead><tbody><tr><td>LDAS-CTRL</td><td>None</td><td>January 2018 to December 2021</td></tr><tr><td>LDAS-L-VOD</td><td>SMOS L-VOD (1.4GHz)</td><td>January 2018 to November 2021</td></tr><tr><td>LDAS-C-VOD</td><td>AMSR2 C-VOD (6.9GHz)</td><td>January 2018 to December 2021</td></tr><tr><td>LDAS-X-VOD</td><td>AMSR2 X-VOD (10.7GHz)</td><td>January 2018 to December 2021</td></tr></tbody></table>
in the different layers are assumed uncorrelated when significant correlations are expected to exist.

The VOD observations (see section 4.2 for details) are ingested into the system once per day. In the absence of a physically-based observation operator to transform the background LAI into simulated VOD a rescaling from VOD to LAI is required. To do this a CDF-matching approach is used which is similar to the approach used to transform ASCAT surface soil moisture to the model soil moisture as part of the existing ASCAT soil moisture assimilation (Scipal et al., 2008). The CDF-matching is done using 2 years of training data (2016-17) to match the VOD observations to the climatology of the time-varying CONFESS LAI data.

The VOD observations are collocated with the model grid points and the CDF-matching parameters are calculated separately for each model grid point to account for geographical variations, and separately for each month of the year using a three-month rolling average to account for seasonal variations whilst preventing large jumps in the conversions for consecutive months.

For ASCAT surface soil moisture the CDF is matched both using the mean and standard deviation but for VOD/LAI matching this led to poor results due to a lack of month-to-month variability in the CONFESS LAI dataset, so the VOD/LAI CDF-matching just matches the means as shown in the equation below:

$$VOD' = (E(LAI) / E(VOD)) * VOD \quad (4.3)$$

where VOD' is the rescaled VOD, E(LAI) is the mean of the LAI, E(VOD) is the mean of the VOD and VOD represents the original VOD observations.

Once the CDF-matching is performed the transformed VOD observations are in the same units and have the same dynamic range as the model LAI so the observation operator (H) is simply the identity matrix. For simplicity, no correlations are assumed between the LAI and the other soil moisture observations are assumed and, vice versa, no correlations are assumed between the soil moisture and the transformed VOD observations. The specification of non-zero correlations in these blocks of the Jacobian matrix could be considered in future work.

### 4.1.3 Experiment setup

The experiments are run for 4 years between 2018 and 2021. This period is motivated by the availability of the observations (not to overlap with the CDF-matching training period 2016-17) and to sample multiple annual cycles of vegetation. The experiments are run at a resolution of T<sub>co</sub>319 (~35km) which is significantly coarser than the operational resolution of T<sub>co</sub>1279 (~9km). This resolution is purely chosen for pragmatic reasons to enable the running of multiple experimental configurations for such a long period of time, and is very close to the standard ECMWF testing resolution of T<sub>co</sub>399. Three separate experiment configurations are run in addition to the control run as detailed in Table 3.

<table><caption>Table 3: Summary of offline ECLand LDAS experiments</caption><thead><tr><td>Experiment ID</td><td>VOD observations assimilated</td><td>Period</td></tr></thead><tbody><tr><td>LDAS-CTRL</td><td>None</td><td>January 2018 to December 2021</td></tr><tr><td>LDAS-L-VOD</td><td>SMOS L-VOD (1.4GHz)</td><td>January 2018 to November 2021</td></tr><tr><td>LDAS-C-VOD</td><td>AMSR2 C-VOD (6.9GHz)</td><td>January 2018 to December 2021</td></tr><tr><td>LDAS-X-VOD</td><td>AMSR2 X-VOD (10.7GHz)</td><td>January 2018 to December 2021</td></tr></tbody></table>
Note that the SMOS L-band data used were only available until the end of November 2021 so the experiment using this data was truncated by 1 month compared to the other experiments. All the experiments and the control analyse soil moisture were obtained as described in Section 4.1.1.

To measure the impact of the dynamically updated daily LAI analysis on the NWP forecasts the LAI analysis was provided to the NWP system in place of the usual monthly climatology. Then several coupled IFS experiments were run in "forecast-only" mode where the atmospheric analysis is read in from an external source (here the operational system) and the land surface analysis is read in from the offline LDAS experiments detailed in Table 3 and consistent with the LAI analysis produced in each of these experiments. These analyses are then used to initialise the NWP forecasts which are run once per day. The details of these experiments can be seen in Table 4.

<table><caption>Table 4: Summary of IFS coupled forecast experiments</caption><thead><tr><td>Experiment ID</td><td>Land surface analysis and VOD source</td><td>Period</td></tr></thead><tbody><tr><td>IFS-CTRL</td><td>LDAS-CTRL</td><td>January 2018 to December 2021</td></tr><tr><td>IFS-L-VOD</td><td>LDAS-L-VOD</td><td>January 2018 to November 2021</td></tr><tr><td>IFS-C-VOD</td><td>LDAS-C-VOD</td><td>January 2018 to December 2021</td></tr><tr><td>IFS-X-VOD</td><td>LDAS-X-VOD</td><td>January 2018 to December 2021</td></tr></tbody></table>

## 4.2 Data

### 4.2.1 Observations

There are two main categories of satellite observations which have significant sensitivity to vegetation parameters. The first are optical sensors measuring at visible frequencies of the electromagnetic spectrum (e.g. Ocean and Land Cover Imager (OLCI) on the Sentinel-3 series of satellites). These instruments provide a direct measurement of LAI but the availability of usable observations is limited to daylight hours and by the presence of clouds. Therefore, it takes many days to get reliable global coverage. Observations from these sensors contribute to the monthly climatologies which are currently used in the NWP system.

The second are passive microwave (MW) sensors, such as for example the Advanced Microwave Scanning Radiometer-2 (AMSR-2) onboard the Global Change Observation Mission for Water (GCOM-W) satellite. The emitted microwave radiation is sensitive to the water content of the vegetation and therefore provides indirect information on the vegetation parameters themselves. At the lowest MW frequencies there is little to no sensitivity to the atmospheric state and no sensitivity to the time of day so these instruments provide near global coverage every 12 hours. This is the main reason to assimilate MW observations to produce a daily LAI analysis.

Observation operators to transform land-surface model parameters into simulated top of atmosphere MW radiances are not yet mature or accurate enough to consider directly assimilating level 1 MW radiances to analyse vegetation parameters (or other land surface model variables). This is due to the heterogeneity of the land-surface conditions and the complex relationship between the surface characteristics and the MW emissivity and
Note that the SMOS L-band data used were only available until the end of November 2021 so the experiment using this data was truncated by 1 month compared to the other experiments. All the experiments and the control analyse soil moisture were obtained as described in Section 4.1.1.

To measure the impact of the dynamically updated daily LAI analysis on the NWP forecasts the LAI analysis was provided to the NWP system in place of the usual monthly climatology. Then several coupled IFS experiments were run in "forecast-only" mode where the atmospheric analysis is read in from an external source (here the operational system) and the land surface analysis is read in from the offline LDAS experiments detailed in Table 3 and consistent with the LAI analysis produced in each of these experiments. These analyses are then used to initialise the NWP forecasts which are run once per day. The details of these experiments can be seen in Table 4.

<table><caption>Table 4: Summary of IFS coupled forecast experiments</caption><thead><tr><td>Experiment ID</td><td>Land surface analysis and VOD source</td><td>Period</td></tr></thead><tbody><tr><td>IFS-CTRL</td><td>LDAS-CTRL</td><td>January 2018 to December 2021</td></tr><tr><td>IFS-L-VOD</td><td>LDAS-L-VOD</td><td>January 2018 to November 2021</td></tr><tr><td>IFS-C-VOD</td><td>LDAS-C-VOD</td><td>January 2018 to December 2021</td></tr><tr><td>IFS-X-VOD</td><td>LDAS-X-VOD</td><td>January 2018 to December 2021</td></tr></tbody></table>

## 4.2 Data

### 4.2.1 Observations

There are two main categories of satellite observations which have significant sensitivity to vegetation parameters. The first are optical sensors measuring at visible frequencies of the electromagnetic spectrum (e.g. Ocean and Land Cover Imager (OLCI) on the Sentinel-3 series of satellites). These instruments provide a direct measurement of LAI but the availability of usable observations is limited to daylight hours and by the presence of clouds. Therefore, it takes many days to get reliable global coverage. Observations from these sensors contribute to the monthly climatologies which are currently used in the NWP system.

The second are passive microwave (MW) sensors, such as for example the Advanced Microwave Scanning Radiometer-2 (AMSR-2) onboard the Global Change Observation Mission for Water (GCOM-W) satellite. The emitted microwave radiation is sensitive to the water content of the vegetation and therefore provides indirect information on the vegetation parameters themselves. At the lowest MW frequencies there is little to no sensitivity to the atmospheric state and no sensitivity to the time of day so these instruments provide near global coverage every 12 hours. This is the main reason to assimilate MW observations to produce a daily LAI analysis.

Observation operators to transform land-surface model parameters into simulated top of atmosphere MW radiances are not yet mature or accurate enough to consider directly assimilating level 1 MW radiances to analyse vegetation parameters (or other land surface model variables). This is due to the heterogeneity of the land-surface conditions and the complex relationship between the surface characteristics and the MW emissivity and
penetration depth. Therefore, here a derived level 2 vegetation product is used instead of the raw level 1 radiances. This level 2 product is the vegetation optical depth (VOD) which is retrieved from the measured MW radiances and is a measure of how much the vegetation attenuates the MW radiation emitted by the surface. VOD is sensitive to the water content of the above ground vegetation in addition to the vegetation type and density.

VOD at different frequencies is sensitive to slightly different characteristics of the vegetation. For example at C-band (4-8GHz) and X-band (8-12GHz) the VOD is more sensitive to the vegetation cover and therefore strongly correlated with the LAI. At lower frequencies such as L-band (1-2GHz) the VOD is more sensitive and better correlated with above ground biomass (AGB) (Rodríguez-Fernández et al, 2018), but the measured value does not saturate as quickly as at the slightly higher frequencies, see Figure 8 which shows that for higher values of X-band VOD there remains relatively large variability in the L-band VOD.

<figure><img src="image_73.png" /><figcaption>Figure 8: Correlation between collocated SMOS L-band VOD and AMSR2 X-band VOD observations, for 1st January 2016</figcaption></figure>

In this study, retrieved VOD from three different frequencies and two different instruments are used. L-band VOD (1.41GHz) is used from the Microwave Imaging Radiometer using Aperture Synthesis (MIRAS) instrument onboard the ESA Soil Moisture Ocean Salinity (SMOS) satellite. In addition, C-band VOD (6.9GHz) and X-band VOD (10.65GHz) are used from the AMSR-2 instrument onboard the GCOM-W satellite. Both these satellites are in sun-synchronous polar orbits and therefore provide near complete global coverage every 24 hours, see Figure 9.
penetration depth. Therefore, here a derived level 2 vegetation product is used instead of the raw level 1 radiances. This level 2 product is the vegetation optical depth (VOD) which is retrieved from the measured MW radiances and is a measure of how much the vegetation attenuates the MW radiation emitted by the surface. VOD is sensitive to the water content of the above ground vegetation in addition to the vegetation type and density.

VOD at different frequencies is sensitive to slightly different characteristics of the vegetation. For example at C-band (4-8GHz) and X-band (8-12GHz) the VOD is more sensitive to the vegetation cover and therefore strongly correlated with the LAI. At lower frequencies such as L-band (1-2GHz) the VOD is more sensitive and better correlated with above ground biomass (AGB) (Rodríguez-Fernández et al, 2018), but the measured value does not saturate as quickly as at the slightly higher frequencies, see Figure 8 which shows that for higher values of X-band VOD there remains relatively large variability in the L-band VOD.

<figure><img src="image_74.png" /><figcaption>Figure 8: Correlation between collocated SMOS L-band VOD and AMSR2 X-band VOD observations, for 1st January 2016</figcaption></figure>

In this study, retrieved VOD from three different frequencies and two different instruments are used. L-band VOD (1.41GHz) is used from the Microwave Imaging Radiometer using Aperture Synthesis (MIRAS) instrument onboard the ESA Soil Moisture Ocean Salinity (SMOS) satellite. In addition, C-band VOD (6.9GHz) and X-band VOD (10.65GHz) are used from the AMSR-2 instrument onboard the GCOM-W satellite. Both these satellites are in sun-synchronous polar orbits and therefore provide near complete global coverage every 24 hours, see Figure 9.
<figure><img src="image_74.png" /><figcaption>Figure 9: Map of daily coverage of AMSR2 6.9GHz VOD observations on 1st July 2016</figcaption></figure>

4.2.2 Quality control

Prior to the assimilation stringent quality control procedures must be applied to the observational data in order to avoid gross errors in the data leading to a degraded LAI analysis. The observational data itself contains various flags to indicate where the data may not be of good quality. One of these flags relates to radio frequency interference (RFI). Passive microwave instruments generally measure the passively emitted MW radiation in spectral frequency bands that are protected from man-made emissions and therefore the measured radiation should only contain natural signals related to geophysical effects (e.g. Weston and de Rosnay 2021). However, in some protected bands it is possible to have out-of-band emissions from man-made sources in neighbouring bands which then cause RFI in the passive measurements. It is important to screen out these contaminated observations because they could lead to anomalous analysis increments and unphysical conditions. Figure 10 shows the global distribution of AMSR2 observations affected by RFI as indicated by the RFI flags which accompany the data. The higher the value of the flag, the more severe the contamination. There are VOD observations derived from the 6.9GHz and 7.3GHz channels on AMSR2 but here only the 6.9GHz channel is used because the 7.3GHz channel is more affected by RFI as shown in figure 10.
<figure><img src="image_75.png" /><figcaption>Figure 9: Map of daily coverage of AMSR2 6.9GHz VOD observations on 1st July 2016</figcaption></figure>

4.2.2 Quality control

Prior to the assimilation stringent quality control procedures must be applied to the observational data in order to avoid gross errors in the data leading to a degraded LAI analysis. The observational data itself contains various flags to indicate where the data may not be of good quality. One of these flags relates to radio frequency interference (RFI). Passive microwave instruments generally measure the passively emitted MW radiation in spectral frequency bands that are protected from man-made emissions and therefore the measured radiation should only contain natural signals related to geophysical effects (e.g. Weston and de Rosnay 2021). However, in some protected bands it is possible to have out-of-band emissions from man-made sources in neighbouring bands which then cause RFI in the passive measurements. It is important to screen out these contaminated observations because they could lead to anomalous analysis increments and unphysical conditions. Figure 10 shows the global distribution of AMSR2 observations affected by RFI as indicated by the RFI flags which accompany the data. The higher the value of the flag, the more severe the contamination. There are VOD observations derived from the 6.9GHz and 7.3GHz channels on AMSR2 but here only the 6.9GHz channel is used because the 7.3GHz channel is more affected by RFI as shown in figure 10.
<figure><img src="image_75.png" /><figcaption>Figure 10: Maps of AMSR2 observations contaminated by RFI for the 6.9GHz channel (top) and 7.3GHz channel (bottom), data is from 1st July 2016</figcaption></figure>

Observations also need to be filtered out where the VOD retrieval is outside of its physical bounds such as when the retrieved value is less than 0. In addition, frozen surfaces are screened out by a check on the surface temperature being less than 273K and where there is snow or ice on the surface.

Within the assimilation system itself there are additional quality control procedures to reject observations when the observed value and the model value are too different. Here, the rescaled VOD observations are rejected if the first-guess departure is larger than 1m²/m² in LAI space.

### 4.2.3 CDF-matching

The methodology for the CDF-matching procedure is described in section 4.1.2. The left panels of Figure 11 shows that the VOD observations are reasonably well correlated with the CONFESS LAI albeit with some significant spread around the diagonal and that the dynamic range of the VOD and LAI values are quite different.
<figure><img src="image_76.png" /><figcaption>Figure 10: Maps of AMSR2 observations contaminated by RFI for the 6.9GHz channel (top) and 7.3GHz channel (bottom), data is from 1st July 2016</figcaption></figure>

Observations also need to be filtered out where the VOD retrieval is outside of its physical bounds such as when the retrieved value is less than 0. In addition, frozen surfaces are screened out by a check on the surface temperature being less than 273K and where there is snow or ice on the surface.

Within the assimilation system itself there are additional quality control procedures to reject observations when the observed value and the model value are too different. Here, the rescaled VOD observations are rejected if the first-guess departure is larger than 1m²/m² in LAI space.

### 4.2.3 CDF-matching

The methodology for the CDF-matching procedure is described in section 4.1.2. The left panels of Figure 11 shows that the VOD observations are reasonably well correlated with the CONFESS LAI albeit with some significant spread around the diagonal and that the dynamic range of the VOD and LAI values are quite different.
<figure><img src="image_76.png" /><figcaption>Figure 11: Correlations between VOD observations and CONFESS LAI for high vegetation types for July 2018. The different panels show: upper left SMOS L-band VOD; upper right rescaled SMOS-L band VOD; middle left AMSR2 C-band VOD; middle right rescaled AMSR2 C-band VOD; lower right AMSR2 X-band VOD; lower left rescaled AMSR2 X-band VOD</figcaption></figure>

This is one of the reasons why, in the absence of an accurate physically-based observation operator, that a rescaling is required. Interestingly it appears that, conversely to previous findings, the L-band VOD is better correlated with the LAI than the C- or X-band VOD. The right panels of Figure 11 show that the correlation is much enhanced once the CDF-matching is applied to the data, as expected.

After the rescaling the L-band VOD is now not as well correlated as the C- and X-band VOD. The CDF-matching is designed to correct the bias between the VOD observations and LAI
<figure><img src="image_77.png" /><figcaption>Figure 11: Correlations between VOD observations and CONFESS LAI for high vegetation types for July 2018. The different panels show: upper left SMOS L-band VOD; upper right rescaled SMOS-L band VOD; middle left AMSR2 C-band VOD; middle right rescaled AMSR2 C-band VOD; lower right AMSR2 X-band VOD; lower left rescaled AMSR2 X-band VOD</figcaption></figure>

This is one of the reasons why, in the absence of an accurate physically-based observation operator, that a rescaling is required. Interestingly it appears that, conversely to previous findings, the L-band VOD is better correlated with the LAI than the C- or X-band VOD. The right panels of Figure 11 show that the correlation is much enhanced once the CDF-matching is applied to the data, as expected.

After the rescaling the L-band VOD is now not as well correlated as the C- and X-band VOD. The CDF-matching is designed to correct the bias between the VOD observations and LAI
which Figure 11 illustrates is being successfully achieved. The random and temporal variability of the observations should be retained, which is what is aimed to be exploited within the assimilation.

## 4.3 Results

To check the performance of the offline LDAS experiments the LAI analysis increments can be analysed. Figures 12 and 13 show that there are significant seasonal variations in the global distribution of LAI analysis increments.

<figure><img src="image_77.png" /><figcaption>**Figure 12:** Monthly mean analysis increments of leaf area index for high vegetation types for January 2018 (upper left), April 2018 (upper right), July 2018 (lower left), October 2018 (lower right)</figcaption></figure>

The most obvious differences are the large areas of zero analysis increments in the Northern hemisphere winter due to snow and frozen ground where the VOD observations are screened out and not assimilated. For the LAI for high vegetation types the mean analysis increments appear to be skewed to the negative with the analysis reducing the LAI in many more areas than increasing the LAI. This may be due to an extreme value check within the analysis where the analysed LAI (high) is prevented from becoming larger than a physically realistic maximum value of 8m²/m². Figure 13 shows that the increments for LAI (low) are more symmetric suggesting that this check (with a maximum value of 6.5m²/m²) is triggered on fewer points than for the LAI (high).
which Figure 11 illustrates is being successfully achieved. The random and temporal variability of the observations should be retained, which is what is aimed to be exploited within the assimilation.

## 4.3 Results

To check the performance of the offline LDAS experiments the LAI analysis increments can be analysed. Figures 12 and 13 show that there are significant seasonal variations in the global distribution of LAI analysis increments.

<figure><img src="image_78.png" /><figcaption>**Figure 12:** Monthly mean analysis increments of leaf area index for high vegetation types for January 2018 (upper left), April 2018 (upper right), July 2018 (lower left), October 2018 (lower right)</figcaption></figure>

The most obvious differences are the large areas of zero analysis increments in the Northern hemisphere winter due to snow and frozen ground where the VOD observations are screened out and not assimilated. For the LAI for high vegetation types the mean analysis increments appear to be skewed to the negative with the analysis reducing the LAI in many more areas than increasing the LAI. This may be due to an extreme value check within the analysis where the analysed LAI (high) is prevented from becoming larger than a physically realistic maximum value of 8m²/m². Figure 13 shows that the increments for LAI (low) are more symmetric suggesting that this check (with a maximum value of 6.5m²/m²) is triggered on fewer points than for the LAI (high).
<figure><img src="image_78.png" /><figcaption>Figure 13: As Figure 12 but for leaf area index for low vegetation types</figcaption></figure>

Figures 12 and 13 also show that the seasonal cycle in the tropics has an impact on the sign of the increments with positive increments generally in the inter-tropical convergence zone (ITCZ) and negative increments to the North and South of the ITCZ. The increments over the mid-latitudes are more mixed and don't seem to follow any obvious systematic pattern, possibly responding to inter-annual differences rather than intra-annual signals.

Figure 14 provides further evidence that the mid-latitude increments are responding to inter-annual differences caused by meteorological events.

<figure><img src="image_79.png" /><figcaption>Figure 14: Monthly mean analysis increments of leaf area index for low vegetation types for July 2021 from LDAS-X-VOD (left) and LDAS-L-VOD (right)</figcaption></figure>

Over western Europe in July 2021 there is a significant positive mean analysis increment which corresponds to the area affected by very heavy and sustained precipitation which caused devastating floods centred over North-West Germany and Belgium. Figure 14 also shows that the analysis increments are smaller for the experiment assimilating the L-band VOD compared to the X-band VOD which is a signal which is repeated for most other monthly
<figure><img src="image_79.png" /><figcaption>Figure 13: As Figure 12 but for leaf area index for low vegetation types</figcaption></figure>

Figures 12 and 13 also show that the seasonal cycle in the tropics has an impact on the sign of the increments with positive increments generally in the inter-tropical convergence zone (ITCZ) and negative increments to the North and South of the ITCZ. The increments over the mid-latitudes are more mixed and don't seem to follow any obvious systematic pattern, possibly responding to inter-annual differences rather than intra-annual signals.

Figure 14 provides further evidence that the mid-latitude increments are responding to inter-annual differences caused by meteorological events.

<figure><img src="image_80.png" /><figcaption>Figure 14: Monthly mean analysis increments of leaf area index for low vegetation types for July 2021 from LDAS-X-VOD (left) and LDAS-L-VOD (right)</figcaption></figure>

Over western Europe in July 2021 there is a significant positive mean analysis increment which corresponds to the area affected by very heavy and sustained precipitation which caused devastating floods centred over North-West Germany and Belgium. Figure 14 also shows that the analysis increments are smaller for the experiment assimilating the L-band VOD compared to the X-band VOD which is a signal which is repeated for most other monthly
mean increments. This is due to the smaller number of observations available for the SMOS L-band VOD compared to the AMSR2 X-band VOD. The difference in the analysis increments for soil moisture in the different experiments has also been analysed (not shown) but the differences are extremely small. This is due to the simple experimental setup not allowing for background error correlations between LAI and soil moisture, as well as not allowing the VOD observations to modify the soil moisture through cross correlations in the Jacobian matrix.

Overall, the mean LAI analysis increments look reasonable and are not making unrealistically large changes to the background climatological LAI which is encouraging and shows that the offline LDAS is performing in a satisfactory manner to produce realistic analysed LAI fields.

## 4.4 Evaluation

### 4.4.1 Impact on NWP

To measure the impact of a dynamically updated LAI analysis on NWP forecasts several experiments were run as shown in Table 4. To assess this impact the forecasts at different lead times from the experiments and control are compared by verifying against the ECMWF operational analysis at the same validity time. For example a 48 hour forecast from 00UTC on 1st July 2018 is compared to the operational analysis from 00UTC on 3rd July 2018. This is done for both the experiments and control and the differences to the verifying analysis are then compared to assess whether the experiment results in a better or worse match between the forecasts and the verifying analysis. Figure 15 shows that the IFS-X-VOD experiment results in smaller root mean squared errors (RMSE) for 2 metre temperature forecasts than the IFS-CTRL over land where the largest impacts are expected.

<figure><img src="image_80.png" /><figcaption>Figure 15: Global change in RMSE of 2 metre temperature forecasts from T+48 to T+120 hours between IFS-X-VOD and IFS-CTRL. Statistics cover January 2018 to December 2021</figcaption></figure>

The largest impacts (blue colours) are over the Amazon rainforest in South America (with up to a 5% reduction in RMSE), tropical Africa and India. Over the extra-tropics the impacts are
mean increments. This is due to the smaller number of observations available for the SMOS L-band VOD compared to the AMSR2 X-band VOD. The difference in the analysis increments for soil moisture in the different experiments has also been analysed (not shown) but the differences are extremely small. This is due to the simple experimental setup not allowing for background error correlations between LAI and soil moisture, as well as not allowing the VOD observations to modify the soil moisture through cross correlations in the Jacobian matrix.

Overall, the mean LAI analysis increments look reasonable and are not making unrealistically large changes to the background climatological LAI which is encouraging and shows that the offline LDAS is performing in a satisfactory manner to produce realistic analysed LAI fields.

## 4.4 Evaluation

### 4.4.1 Impact on NWP

To measure the impact of a dynamically updated LAI analysis on NWP forecasts several experiments were run as shown in Table 4. To assess this impact the forecasts at different lead times from the experiments and control are compared by verifying against the ECMWF operational analysis at the same validity time. For example a 48 hour forecast from 00UTC on 1st July 2018 is compared to the operational analysis from 00UTC on 3rd July 2018. This is done for both the experiments and control and the differences to the verifying analysis are then compared to assess whether the experiment results in a better or worse match between the forecasts and the verifying analysis. Figure 15 shows that the IFS-X-VOD experiment results in smaller root mean squared errors (RMSE) for 2 metre temperature forecasts than the IFS-CTRL over land where the largest impacts are expected.

<figure><img src="image_81.png" /><figcaption>Figure 15: Global change in RMSE of 2 metre temperature forecasts from T+48 to T+120 hours between IFS-X-VOD and IFS-CTRL. Statistics cover January 2018 to December 2021</figcaption></figure>

The largest impacts (blue colours) are over the Amazon rainforest in South America (with up to a 5% reduction in RMSE), tropical Africa and India. Over the extra-tropics the impacts are
more neutral. This is due to the screening of VOD observations over mid- to high-latitudes in winter when the land surface is either frozen or snow-covered. Figures 16 and 17 show the impact of the LAI analysis on the lower atmosphere.

<figure><img src="image_81.png" /><figcaption>Figure 16: As figure 15 but for 850hPa temperature forecasts</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>Figure 17: As figure 15 but for 850hPa relative humidity forecasts</figcaption></figure>
more neutral. This is due to the screening of VOD observations over mid- to high-latitudes in winter when the land surface is either frozen or snow-covered. Figures 16 and 17 show the impact of the LAI analysis on the lower atmosphere.

<figure><img src="image_82.png" /><figcaption>Figure 16: As figure 15 but for 850hPa temperature forecasts</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>Figure 17: As figure 15 but for 850hPa relative humidity forecasts</figcaption></figure>
Here the impact is more mixed with consistent improvements to temperature and relative humidity forecasts over the Amazon, Africa and Western Europe. However, there are some degradations over India.

The negative impacts over India are mostly coming from the Northern hemisphere winter (DJF) which is a period where the analysis increments are predominantly negative as shown in Figures 12 and 13 so it is possible that the reduction in LAI is being overdone.

Figure 18 shows the impact on the mean errors for temperature and relative humidity in the lower atmosphere is significant.

<figure><img src="image_83.png" /><figcaption>Figure 18: Mean temperature (left) and relative humidity (right) forecast errors at 850hPa (upper) and mean 2 metre temperature forecast errors (lower) for IFS-X-VOD (black) and IFS-CTRL (red). Statistics cover January 2018 to December 2021</figcaption></figure>

For relative humidity the existing moist biases in the forecasts are reduced by 5-10% at 850hPa. For temperature the existing negative bias is reduced at 850hPa and there is a similar signal for 2 metre temperature. This is consistent with the general reduction in LAI seen in Figures 12 and 13 and suggests that the use of the dynamic LAI analysis is having a significant impact on correcting systematic forecast errors compared to the use of a fixed LAI climatology in IFS-CTRL.

So far only results from the IFS-X-VOD experiment have been shown. Figure 19 shows a comparison between IFS-X-VOD, IFS-C-VOD and IFS-L-VOD showing that the results are extremely similar for all three experiments.
Here the impact is more mixed with consistent improvements to temperature and relative humidity forecasts over the Amazon, Africa and Western Europe. However, there are some degradations over India.

The negative impacts over India are mostly coming from the Northern hemisphere winter (DJF) which is a period where the analysis increments are predominantly negative as shown in Figures 12 and 13 so it is possible that the reduction in LAI is being overdone.

Figure 18 shows the impact on the mean errors for temperature and relative humidity in the lower atmosphere is significant.

<figure><img src="image_84.png" /><figcaption>Figure 18: Mean temperature (left) and relative humidity (right) forecast errors at 850hPa (upper) and mean 2 metre temperature forecast errors (lower) for IFS-X-VOD (black) and IFS-CTRL (red). Statistics cover January 2018 to December 2021</figcaption></figure>

For relative humidity the existing moist biases in the forecasts are reduced by 5-10% at 850hPa. For temperature the existing negative bias is reduced at 850hPa and there is a similar signal for 2 metre temperature. This is consistent with the general reduction in LAI seen in Figures 12 and 13 and suggests that the use of the dynamic LAI analysis is having a significant impact on correcting systematic forecast errors compared to the use of a fixed LAI climatology in IFS-CTRL.

So far only results from the IFS-X-VOD experiment have been shown. Figure 19 shows a comparison between IFS-X-VOD, IFS-C-VOD and IFS-L-VOD showing that the results are extremely similar for all three experiments.
<figure><img src="image_84.png" /><figcaption>**Figure 19:** Change in RMSE for 2 metre temperature forecasts between IFS-X-VOD (left), IFS-C-VOD (centre), IFS-L-VOD (right) and IFS-CTRL. Statistics cover January 2018 to December 2021</figcaption></figure>

There are some very small differences from assimilating the different VOD bands, but the overall results are very similar. It is expected that the L-band VOD has slightly different characteristics to the X- and C-band VODs. However, in the methodology used here the VOD observations from the different bands are all CDF-matched to the same LAI dataset, which will temper the differences between them and is the most obvious explanation for the similar results seen in Figure 19.

Figure 20 compares the results for the different seasons.

<figure><img src="image_85.png" /></figure>
<figure><img src="image_85.png" /><figcaption>**Figure 19:** Change in RMSE for 2 metre temperature forecasts between IFS-X-VOD (left), IFS-C-VOD (centre), IFS-L-VOD (right) and IFS-CTRL. Statistics cover January 2018 to December 2021</figcaption></figure>

There are some very small differences from assimilating the different VOD bands, but the overall results are very similar. It is expected that the L-band VOD has slightly different characteristics to the X- and C-band VODs. However, in the methodology used here the VOD observations from the different bands are all CDF-matched to the same LAI dataset, which will temper the differences between them and is the most obvious explanation for the similar results seen in Figure 19.

Figure 20 compares the results for the different seasons.

<figure><img src="image_86.png" /></figure>
<figure><img src="image_86.png" /><figcaption>**Figure 20:** Change in RMSE of 2 metre temperature forecasts between IFS-L-VOD and IFS-CTRL for December-January-February (top), March-April-May (2nd), June-July-August (3rd), September-October-November (bottom). Statistics cover January 2018 to December 2021</figcaption></figure>

In the extra-tropics the largest impacts are in the summer for that hemisphere (e.g. DJF for the Southern hemisphere; JJA for the Northern hemisphere). This is partially because land-atmosphere coupling processes are most active in summer and LAI inter-annual differences will be large in the mid- to high-latitudes so this is where the largest impacts of dynamically updating LAI would be expected to be. There is very little impact in the extra-tropics in the winter for that hemisphere due to the screening of the VOD observations over frozen and snow-covered ground and limited evapotranspiration. The impact in the tropics is fairly similar all year round where seasonal differences in meteorology and land-surface states are smaller. However, there is a slightly stronger positive impact in JJA and SON which is driven by larger improvements to 2 metre temperature forecasts over the Amazon during these seasons. One slightly surprising result is the lack of impact in the extra-tropics in the transition seasons MAM and SON where inter-annual differences are expected to be at their largest. However, Alessandri & van Oorschot (2022) found that the largest inter-annual variability in the CONFESS time-varying vegetation dataset occurs in the tropics, particularly South America, where the largest NWP impacts are found here.

Overall, the impact on NWP of assimilating satellite data to dynamically updated vegetation is mixed although there are some promising signs, especially the impact on 2 metre temperature forecasts and the improved mean forecast errors for near-surface relative humidity. This is encouraging and shows that microwave radiances can have an impact on NWP via improved analysis of the vegetation parameters. Future work as part of the CORSO project will aim at further improving the vegetation parameters in the model by developing an observation operator to transform model variables, such as soil moisture and LAI, into simulated level 1 satellite observations. This should enable the assimilation of level 1 satellite observations, such as passive microwave brightness temperatures and active backscatter observations, to analyse the soil moisture and vegetation parameters simultaneously and consistently.
<figure><img src="image_87.png" /><figcaption>**Figure 20:** Change in RMSE of 2 metre temperature forecasts between IFS-L-VOD and IFS-CTRL for December-January-February (top), March-April-May (2nd), June-July-August (3rd), September-October-November (bottom). Statistics cover January 2018 to December 2021</figcaption></figure>

In the extra-tropics the largest impacts are in the summer for that hemisphere (e.g. DJF for the Southern hemisphere; JJA for the Northern hemisphere). This is partially because land-atmosphere coupling processes are most active in summer and LAI inter-annual differences will be large in the mid- to high-latitudes so this is where the largest impacts of dynamically updating LAI would be expected to be. There is very little impact in the extra-tropics in the winter for that hemisphere due to the screening of the VOD observations over frozen and snow-covered ground and limited evapotranspiration. The impact in the tropics is fairly similar all year round where seasonal differences in meteorology and land-surface states are smaller. However, there is a slightly stronger positive impact in JJA and SON which is driven by larger improvements to 2 metre temperature forecasts over the Amazon during these seasons. One slightly surprising result is the lack of impact in the extra-tropics in the transition seasons MAM and SON where inter-annual differences are expected to be at their largest. However, Alessandri & van Oorschot (2022) found that the largest inter-annual variability in the CONFESS time-varying vegetation dataset occurs in the tropics, particularly South America, where the largest NWP impacts are found here.

Overall, the impact on NWP of assimilating satellite data to dynamically updated vegetation is mixed although there are some promising signs, especially the impact on 2 metre temperature forecasts and the improved mean forecast errors for near-surface relative humidity. This is encouraging and shows that microwave radiances can have an impact on NWP via improved analysis of the vegetation parameters. Future work as part of the CORSO project will aim at further improving the vegetation parameters in the model by developing an observation operator to transform model variables, such as soil moisture and LAI, into simulated level 1 satellite observations. This should enable the assimilation of level 1 satellite observations, such as passive microwave brightness temperatures and active backscatter observations, to analyse the soil moisture and vegetation parameters simultaneously and consistently.
4.4.2 Impact on carbon cycle

Figure 21 shows that in Northern hemisphere summer 2018 the absolute biases of accumulated CO2 gross primary production (GPP) against FLUXCOM are increased in the experiment using the analysed LAI compared to using the climatological LAI.

<figure><img src="image_87.png" /><figcaption>Figure 21: Differences in accumulated CO2 gross primary production absolute bias between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</figcaption></figure>

The increased biases are visible almost globally except for a few small reductions over Canada. However, in Northern hemisphere spring 2021 the results are more balanced with areas of increased bias in India, Asia and Western Europe but some reductions in bias over tropical Africa and the Eastern parts of both North and South America. These different results for different years and seasons potentially show the sensitivity the LAI analysis to inter-annual variability. These more positive results are also visible in summer 2021 (not shown) but for most other seasons and years the results are less positive and are more closely aligned with the 2018 results.

Figure 22 shows that in the Northern hemisphere summer 2018 the correlations between the model accumulated CO2 GPP and the FLUXCOM data are reduced in the experiment using the analysed LAI compared to the control using the climatological LAI.

<figure><img src="image_88.png" /><figcaption>Figure 22: Differences in accumulated CO2 gross primary production correlation between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</figcaption></figure>
4.4.2 Impact on carbon cycle

Figure 21 shows that in Northern hemisphere summer 2018 the absolute biases of accumulated CO2 gross primary production (GPP) against FLUXCOM are increased in the experiment using the analysed LAI compared to using the climatological LAI.

<figure><img src="image_88.png" /><figcaption>**Figure 21:** Differences in accumulated CO2 gross primary production absolute bias between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</figcaption></figure>

The increased biases are visible almost globally except for a few small reductions over Canada. However, in Northern hemisphere spring 2021 the results are more balanced with areas of increased bias in India, Asia and Western Europe but some reductions in bias over tropical Africa and the Eastern parts of both North and South America. These different results for different years and seasons potentially show the sensitivity the LAI analysis to inter-annual variability. These more positive results are also visible in summer 2021 (not shown) but for most other seasons and years the results are less positive and are more closely aligned with the 2018 results.

Figure 22 shows that in the Northern hemisphere summer 2018 the correlations between the model accumulated CO2 GPP and the FLUXCOM data are reduced in the experiment using the analysed LAI compared to the control using the climatological LAI.

<figure><img src="image_89.png" /><figcaption>**Figure 22:** Differences in accumulated CO2 gross primary production correlation between IFS-C-VOD and IFS-CTRL measured against FLUXCOM data for July-August-September 2018 (left) and April-May-June 2021 (right)</figcaption></figure>
The largest correlation degradations are over Southern Africa. In Northern hemisphere spring 2021 there are still areas of degraded correlations particularly over the Sahel but there are some improved correlations too over tropical Africa. The general pattern is that the correlations are degraded though. The results for the differences in root mean squared error (RMSE) mirror those for the biases and correlation (not shown). In addition to the maps the overall global budgets of GPP were calculated and compared to the FLUXCOM dataset. There are very small changes in the global budgets and any changes generally go in the wrong direction as they further increase a negative bias between the IFS modelled GPP and the FLUXCOM GPP. This is not a huge surprise given the negatively skewed LAI increments in Figure 12 and, as previously mentioned, the possible influence of the extreme value check for large values of LAI.

The very small changes in global budgets are also to be expected given that a CDF-matching approach is used which effectively will match the global biases in the VOD observations to the existing climatological LAI biases. In the absence of a physically based observation operator to transform model LAI into simulated VOD this was the best option available and has been used by previous studies too e.g. Mucia et al (2022). This limited impact on the carbon fluxes could also be due to the very small changes in soil moisture by assimilating the VOD observations, as discussed briefly in section 4.3. This was due to the design of the system with the VOD observations only able to directly modify the LAI analysis due to the lack of correlations in the background error covariance matrices and Jacobian matrices between the LAI and the soil moisture. In future this could be addressed by including non-zero correlations between soil moisture and LAI in a similar system.

Overall, it is interesting that the dynamically updated LAI results in improvements in the meteorological scores while apparently degrading the carbon fluxes when verified against observations. It is possible there are compensating biases within the NWP system which contribute to these conflicting results. As mentioned at the end of section 4.4.1, the development of an observation operator for vegetation and soil moisture should allow the assimilation of level 1 measurements to simultaneously analyse LAI and soil moisture. This approach using a coupled soil-vegetation-atmosphere assimilation will enhance the consistency between the meteorological and carbon flux analyses in the IFS.
The largest correlation degradations are over Southern Africa. In Northern hemisphere spring 2021 there are still areas of degraded correlations particularly over the Sahel but there are some improved correlations too over tropical Africa. The general pattern is that the correlations are degraded though. The results for the differences in root mean squared error (RMSE) mirror those for the biases and correlation (not shown). In addition to the maps the overall global budgets of GPP were calculated and compared to the FLUXCOM dataset. There are very small changes in the global budgets and any changes generally go in the wrong direction as they further increase a negative bias between the IFS modelled GPP and the FLUXCOM GPP. This is not a huge surprise given the negatively skewed LAI increments in Figure 12 and, as previously mentioned, the possible influence of the extreme value check for large values of LAI.

The very small changes in global budgets are also to be expected given that a CDF-matching approach is used which effectively will match the global biases in the VOD observations to the existing climatological LAI biases. In the absence of a physically based observation operator to transform model LAI into simulated VOD this was the best option available and has been used by previous studies too e.g. Mucia et al (2022). This limited impact on the carbon fluxes could also be due to the very small changes in soil moisture by assimilating the VOD observations, as discussed briefly in section 4.3. This was due to the design of the system with the VOD observations only able to directly modify the LAI analysis due to the lack of correlations in the background error covariance matrices and Jacobian matrices between the LAI and the soil moisture. In future this could be addressed by including non-zero correlations between soil moisture and LAI in a similar system.

Overall, it is interesting that the dynamically updated LAI results in improvements in the meteorological scores while apparently degrading the carbon fluxes when verified against observations. It is possible there are compensating biases within the NWP system which contribute to these conflicting results. As mentioned at the end of section 4.4.1, the development of an observation operator for vegetation and soil moisture should allow the assimilation of level 1 measurements to simultaneously analyse LAI and soil moisture. This approach using a coupled soil-vegetation-atmosphere assimilation will enhance the consistency between the meteorological and carbon flux analyses in the IFS.
# 5 ORCHIDEE

We provide here an overview of the calibration of ORCHIDEE parameters following the assimilation of remote sensing SIF data and *in situ* GPP data. The optimized parameters are then used to perform a simulation of GPPs at global scale, which are compared to those obtained with the standard parameters. The difference highlights the combined constraint brought by GPP and SIF (from Sentinel-5p observations) data to improve the model prediction. Note that the use of GPP data in addition to SIF data is crucial as it helps constraining photosynthesis model parameters and avoids overfitting issues that could else degrade the simulated GPP. Such work has been conducted mainly in the context of the Sen4GPP ESA project with a contribution of CoCO₂ resources and it has mainly been done by Cedric Bacour, Vladislav Bastrikov and Fabienne Maignan.

## 5.1 Methods

We first provide an overview of the calibration of ORCHIDEE parameters with the assimilation of SIF and GPP data. The optimized parameters are then used to perform a simulation of GPPs at global scale, which are compared to those obtained with the standard parameters. The difference highlights the combined constraint brought by GPP and SIF (from Sentinel-5p observations) data to improve the model prediction.

### 5.1.1 ORCHIDEE model and optimized parameters

ORCHIDEE is a mechanistic LSM that simulates the exchanges of carbon, water, and energy between biosphere and atmosphere (Krinner et al., 2005). It is the land surface component of the Earth System Model of Institut Pierre-Simon Laplace IPSL-CM. Photosynthesis and all components of the surface energy and water budgets are calculated at a half-hourly resolution while the dynamics of the carbon storage (including carbon allocation in plant reservoirs, soil carbon dynamics, and litter decomposition) are resolved on a daily basis. Photosynthesis depends on light availability and CO₂ concentration, soil moisture and temperature, and it is parameterized based on Farquhar et al. (1980) and Collatz et al. (1992) for C3 and C4 plants, respectively. The formulations of Yin and Struik (2009) are used to describe the main photosynthesis processes. As in most LSMs, the spatial distribution of vegetation is represented using fractions of plant functional types (PFTs) for each grid point. Except for the phenology, processes are described with the same governing equations for all PFTs, but usually with different parameter values.

The list of parameters that are optimized is provided in the table below.

<table><caption>Table 5: List of the ORCHIDEE optimized parameters.</caption><thead><tr><th>Name</th><th>Description</th></tr></thead><tbody><tr><td colspan="2">Photosynthesis</td></tr><tr><td>A1</td><td>Empirical factor involved in the calculation of leaf-to-air vapour difference effect on the coupling between CO<sub>2</sub> assimilation and stomatal conductance</td></tr><tr><td>B1</td><td>Empirical factor involved in the calculation of leaf-to-air vapour difference effect</td></tr><tr><td>V<sub>cmax25</sub></td><td>Maximum rate of Rubisco activity-limited carboxylation at 25°C (µmol[CO<sub>2</sub>]·m<sup>-2</sup>·s<sup>-1</sup>)</td></tr><tr><td>SLA</td><td>Specific leaf area (m<sup>2</sup>·g<sup>-1</sup>)</td></tr><tr><td>ALA</td><td>Average leaf angle (°)</td></tr><tr><td>Clumping</td><td>Leaf clumping index (-)</td></tr></tbody></table>
# 5 ORCHIDEE

We provide here an overview of the calibration of ORCHIDEE parameters following the assimilation of remote sensing SIF data and *in situ* GPP data. The optimized parameters are then used to perform a simulation of GPPs at global scale, which are compared to those obtained with the standard parameters. The difference highlights the combined constraint brought by GPP and SIF (from Sentinel-5p observations) data to improve the model prediction. Note that the use of GPP data in addition to SIF data is crucial as it helps constraining photosynthesis model parameters and avoids overfitting issues that could else degrade the simulated GPP. Such work has been conducted mainly in the context of the Sen4GPP ESA project with a contribution of CoCO₂ resources and it has mainly been done by Cedric Bacour, Vladislav Bastrikov and Fabienne Maignan.

## 5.1 Methods

We first provide an overview of the calibration of ORCHIDEE parameters with the assimilation of SIF and GPP data. The optimized parameters are then used to perform a simulation of GPPs at global scale, which are compared to those obtained with the standard parameters. The difference highlights the combined constraint brought by GPP and SIF (from Sentinel-5p observations) data to improve the model prediction.

### 5.1.1 ORCHIDEE model and optimized parameters

ORCHIDEE is a mechanistic LSM that simulates the exchanges of carbon, water, and energy between biosphere and atmosphere (Krinner et al., 2005). It is the land surface component of the Earth System Model of Institut Pierre-Simon Laplace IPSL-CM. Photosynthesis and all components of the surface energy and water budgets are calculated at a half-hourly resolution while the dynamics of the carbon storage (including carbon allocation in plant reservoirs, soil carbon dynamics, and litter decomposition) are resolved on a daily basis. Photosynthesis depends on light availability and CO₂ concentration, soil moisture and temperature, and it is parameterized based on Farquhar et al. (1980) and Collatz et al. (1992) for C3 and C4 plants, respectively. The formulations of Yin and Struik (2009) are used to describe the main photosynthesis processes. As in most LSMs, the spatial distribution of vegetation is represented using fractions of plant functional types (PFTs) for each grid point. Except for the phenology, processes are described with the same governing equations for all PFTs, but usually with different parameter values.

The list of parameters that are optimized is provided in the table below.

<table><caption>Table 5: List of the ORCHIDEE optimized parameters.</caption><thead><tr><th>Name</th><th>Description</th></tr></thead><tbody><tr><td colspan="2">Photosynthesis</td></tr><tr><td>A1</td><td>Empirical factor involved in the calculation of leaf-to-air vapour difference effect on the coupling between CO<sub>2</sub> assimilation and stomatal conductance</td></tr><tr><td>B1</td><td>Empirical factor involved in the calculation of leaf-to-air vapour difference effect</td></tr><tr><td>V<sub>cmax25</sub></td><td>Maximum rate of Rubisco activity-limited carboxylation at 25°C (µmol[CO<sub>2</sub>]·m<sup>-2</sup>·s<sup>-1</sup>)</td></tr><tr><td>SLA</td><td>Specific leaf area (m<sup>2</sup>·g<sup>-1</sup>)</td></tr><tr><td>ALA</td><td>Average leaf angle (°)</td></tr><tr><td>Clumping</td><td>Leaf clumping index (-)</td></tr></tbody></table>
<table><tr><td>ASV</td><td>Offset of the linear temperature acclimation relationship for the entropy parameter of the Vcmax (maximum carbonxylation rate) temperature-dependence function, following Kattge and Knorr (2007) <b></b> <b>erreur ! Source du renvoi introuvable.</b></td></tr><tr><td>ASJ</td><td>Offset of the linear temperature acclimation relationship for the entropy parameter of the Jmax (maximum electron transport rate) temperature-dependence function, following Kattge and Knorr (2007) <b>erreur ! Source du renvoi introuvable.</b></td></tr><tr><td>ARJV</td><td>Offset of the linear temperature acclimation relationship for the ratio Jmax/Vcmax, following Kattge and Knorr (2007)</td></tr><tr><td>GB<sub>REF</sub></td><td>Conductance for the leaf boundary layer (mmol.m<sup>-2</sup>.s<sup>-1</sup>)</td></tr><tr><td colspan="2"><b>Phenology</b></td></tr><tr><td>LAI<sub>MAX</sub></td><td>Maximum LAI value</td></tr><tr><td>T<sub>sen</sub></td><td>Temperature threshold for senescence (°C)</td></tr><tr><td>L<sub>age,crit</sub></td><td>Average critical leaves age (days)</td></tr><tr><td colspan="2"><b>Hydrology</b></td></tr><tr><td>Hum<sub>est</sub></td><td>Root profile (m)</td></tr><tr><td colspan="2"><b>SIF module</b></td></tr><tr><td>k<sub>F</sub></td><td>Relative rate constant of fluorescence</td></tr></table>

All parameters depend on PFT except for $GB_{ref}$. Most parameters constrain photosynthesis and phenology processes; one parameter is related to hydrology and the $k_F$ parameter only impacts the intensity of SIF simulated by ORCHIDEE (it therefore acts like a bias absorber).

The variation range of the parameters is set as ±15% from their a priori value, except for *Clumping*, $V_{cmax25}$ and *ARJV* (±25%), and ALA and $k_F$ (±50%) due to a higher uncertainty associated with their standard values. The prior error in the **B** matrix is determined to correspond to 15% of the variation range following Bacour et al (2023).

## 5.1.2 Data assimilation system

We used the ORCHIDAS Data Assimilation tool specifically developed around the ORCHIDEE LSM (https://orchidas.lsce.ipsl.fr/) (MacBean et al., 2022; Bacour et al., 2023).

The assimilation relies on a Bayesian framework that optimizes ORCHIDEE parameters gathered in a vector **x**, by finding the minimum of a global misfit function $J(\mathbf{x})$. $J(\mathbf{x})$ is a linear combination of the misfit functions associated with each data-stream (SIF and GPP data), also accounting for the difference between the values of the optimized parameters with some prior information $\mathbf{x}_b$.

$$J(\mathbf{x}) = \frac{1}{2} [\sum_o (H(\mathbf{x}) - \mathbf{y}_o)^T \cdot \mathbf{R}_o^{-1} \cdot (H(\mathbf{x}) - \mathbf{y}_o) + (\mathbf{x} - \mathbf{x}_b)^T \cdot \mathbf{B}^{-1} \cdot (\mathbf{x} - \mathbf{x}_b)] \quad (5.1)$$

where $\mathbf{y}_o$ are the observation vectors (with $o$ = SIF or GPP); $H$ is the ORCHIDEE model; $\mathbf{R}_o$ are error covariance matrices characterizing the observation error with respect to the model (therefore including the error in the model structure) associated with data-stream $o$; and **B** is the a priori error covariance matrix on model parameters **x**. Errors on observations and on the model parameters are assumed to be Gaussian and data-streams are considered independent from each other (i.e., co-variances between data streams are neglected).

The minimization of $J(\mathbf{x})$ is done iteratively, using a Genetic Algorithm (Goldberg, 1989), which, unlike gradient-descent approaches, limits the risk of getting stuck into a local minimum of the misfit function (Santaren et al., 2014; Bastrikov et al., 2018). The optimization also accounts for prescribed bounds in the parameter variations.
<table><tr><td>ASV</td><td>Offset of the linear temperature acclimation relationship for the entropy parameter of the Vcmax (maximum carbonxylation rate) temperature-dependence function, following Kattge and Knorr (2007) <b></b> <b>erreur ! Source du renvoi introuvable.</b></td></tr><tr><td>ASJ</td><td>Offset of the linear temperature acclimation relationship for the entropy parameter of the Jmax (maximum electron transport rate) temperature-dependence function, following Kattge and Knorr (2007) <b>erreur ! Source du renvoi introuvable.</b></td></tr><tr><td>ARJV</td><td>Offset of the linear temperature acclimation relationship for the ratio Jmax/Vcmax, following Kattge and Knorr (2007)</td></tr><tr><td>GB<sub>REF</sub></td><td>Conductance for the leaf boundary layer (mmol.m<sup>-2</sup>.s<sup>-1</sup>)</td></tr><tr><td colspan="2"><b>Phenology</b></td></tr><tr><td>LAI<sub>MAX</sub></td><td>Maximum LAI value</td></tr><tr><td>T<sub>sen</sub></td><td>Temperature threshold for senescence (°C)</td></tr><tr><td>L<sub>age,crit</sub></td><td>Average critical leaves age (days)</td></tr><tr><td colspan="2"><b>Hydrology</b></td></tr><tr><td>Hum<sub>est</sub></td><td>Root profile (m)</td></tr><tr><td colspan="2"><b>SIF module</b></td></tr><tr><td>k<sub>F</sub></td><td>Relative rate constant of fluorescence</td></tr></table>

All parameters depend on PFT except for $GB_{ref}$. Most parameters constrain photosynthesis and phenology processes; one parameter is related to hydrology and the $k_F$ parameter only impacts the intensity of SIF simulated by ORCHIDEE (it therefore acts like a bias absorber).

The variation range of the parameters is set as ±15% from their a priori value, except for Clumping, $V_{cmax25}$ and ARJV (±25%), and ALA and kF (±50%) due to a higher uncertainty associated with their standard values. The prior error in the B matrix is determined to correspond to 15% of the variation range following Bacour et al (2023).

## 5.1.2 Data assimilation system

We used the ORCHIDAS Data Assimilation tool specifically developed around the ORCHIDEE LSM (https://orchidas.lsce.ipsl.fr/) (MacBean et al., 2022; Bacour et al., 2023).

The assimilation relies on a Bayesian framework that optimizes ORCHIDEE parameters gathered in a vector $x$, by finding the minimum of a global misfit function $J(x)$. $J(x)$ is a linear combination of the misfit functions associated with each data-stream (SIF and GPP data), also accounting for the difference between the values of the optimized parameters with some prior information $x_b$.

$$J(\mathbf{x}) = \frac{1}{2} [\sum_{o} (H(\mathbf{x}) - \mathbf{y}_{o})^{T} \cdot \mathbf{R}_{o}^{-1} \cdot (H(\mathbf{x}) - \mathbf{y}_{o}) + (\mathbf{x} - \mathbf{x}_{b})^{T} \cdot \mathbf{B}^{-1} \cdot (\mathbf{x} - \mathbf{x}_{b})] \quad (5.1)$$

where $\mathbf{y}_o$ are the observation vectors (with $o$ = SIF or GPP); $H$ is the ORCHIDEE model; $\mathbf{R}_o$ are error covariance matrices characterizing the observation error with respect to the model (therefore including the error in the model structure) associated with data-stream $o$; and $\mathbf{B}$ is the a priori error covariance matrix on model parameters $x$. Errors on observations and on the model parameters are assumed to be Gaussian and data-streams are considered independent from each other (i.e., co-variances between data streams are neglected).

The minimization of $J(x)$ is done iteratively, using a Genetic Algorithm (Goldberg, 1989), which, unlike gradient-descent approaches, limits the risk of getting stuck into a local minimum of the misfit function (Santaren et al., 2014; Bastrikov et al., 2018). The optimization also accounts for prescribed bounds in the parameter variations.
**Model data errors:**

The diagonal of the error covariance matrix on observations R (accounting for the error in model and data assimilated) is populated by the root mean square difference (RMSD) between the observations and the corresponding model simulations using the prior standard parameter values (as in MacBean et al., 2022; Bacour et al., 2023). Given the uneven number of assimilated samples between GPP and SIF data (more GPP data assimilated, in particular due to the use of *in situ* observations at a daily time step), we decreased the prior error on SIF observations so that the misfit functions associated with SIF and GPP are similar at the first iteration using the prior parameter values.

5.1.3 Assimilated data

The assimilations are conducted on a PFT-basis, i.e., by optimizing the model parameters associated with the targeted PFT against GPP data (site scale estimates or satellite-derived data when no *in situ* data are available) and TROPOMI SIF retrievals (Guanter et al., 2019) for a collection of selected pixels. The co-assimilation of these two variables limits parameter overfitting.

SIF data

We use the daily averaged SIF retrievals of the TROPOSIF product (estimates from the 743-758 nm fitting window) (Guanter et al., 2019). Only observations passing the quality flag and associated with view zenith angles below 40° and cloud fraction below 0.2 were considered. The data were aggregated at 8-day/0.25° resolution. Note that the spatial resolution (0.25°) is driven by the resolution of the ERA5 meteorological forcing available at global scale, which is used as input to the ORCHIDEE simulations. We considered retrievals for the 2018-2020 period (hence two years and half).

We selected 15 pixels for each of the 14 vegetation PFTs of ORCHIDEE, with the highest thematic homogeneity possible at 0.25° (based on the PFT map for year 2018; see https://orchidas.lsce.ipsl.fr/dev/lccci/orchidee_pfts.php). One other selection criterion was the spatial distribution of the selected pixels: in the process, we favoured pixels away from each other in order to avoid sampling too much specific regions.

The spatial homogeneity of the selected pixels varies with the PFTs: a very high homogeneity is possible for TrEBF (above 95%), and BoC3GRA (mostly above 80%). The homogeneity is usually above 70% for the selected pixels of TeEBF, TeC3GRA, BoC3GRA and C4GRA. Lower values can be obtained for the other PFTs (but still higher than 55%).

GPP data

For most PFTs, we assimilated daily *in situ* GPP estimates from FLUXNET (La Thuile (Baldocchi et al., 2001) and FLUXNET2015 (Pastorello et al., 2020)), based on the night-time method to partition net ecosystem exchange measurements between respiration and GPP components (Reichstein et al., 2005). The number of sites and their temporal cover strongly depend on the PFT. For instance, only four sites are considered for C4GRA and C4CRO PFTs while twenty-eight sites are used for TeDBF and TeENF.

Only for three PFTs (TrDBF, BoDNF, and TrC3GRA), there was no *in situ* GPP available. We therefore used FLUXSAT-GPP data (Joiner et al., 2018). Similarly to FLUXCOM, FLUXSAT also derives from a data-driven approach relying on FLUXNET measurements and MODIS reflectances in seven spectral bands and calibrated against FLUXNET measurements. One noteworthy difference is that FLUXSAT does not use any meteorological forcing. We used the data for the same 15 pixels as those considered for SIF, and data binned at 0.25°/weekly resolution (also similarly to SIF). Data from 2017-2020 were considered
**Model data errors:**

The diagonal of the error covariance matrix on observations **R** (accounting for the error in model and data assimilated) is populated by the root mean square difference (RMSD) between the observations and the corresponding model simulations using the prior standard parameter values (as in MacBean et al., 2022; Bacour et al., 2023). Given the uneven number of assimilated samples between GPP and SIF data (more GPP data assimilated, in particular due to the use of *in situ* observations at a daily time step), we decreased the prior error on SIF observations so that the misfit functions associated with SIF and GPP are similar at the first iteration using the prior parameter values.

5.1.3 Assimilated data

The assimilations are conducted on a PFT-basis, i.e., by optimizing the model parameters associated with the targeted PFT against GPP data (site scale estimates or satellite-derived data when no *in situ* data are available) and TROPOMI SIF retrievals (Guanter et al., 2019) for a collection of selected pixels. The co-assimilation of these two variables limits parameter overfitting.

SIF data

We use the daily averaged SIF retrievals of the TROPOSIF product (estimates from the 743-758 nm fitting window) (Guanter et al., 2019). Only observations passing the quality flag and associated with view zenith angles below 40° and cloud fraction below 0.2 were considered. The data were aggregated at 8-day/0.25° resolution. Note that the spatial resolution (0.25°) is driven by the resolution of the ERA5 meteorological forcing available at global scale, which is used as input to the ORCHIDEE simulations. We considered retrievals for the 2018-2020 period (hence two years and half).

We selected 15 pixels for each of the 14 vegetation PFTs of ORCHIDEE, with the highest thematic homogeneity possible at 0.25° (based on the PFT map for year 2018; see https://orchidas.lsce.ipsl.fr/dev/lccci/orchidee_pfts.php). One other selection criterion was the spatial distribution of the selected pixels: in the process, we favoured pixels away from each other in order to avoid sampling too much specific regions.

The spatial homogeneity of the selected pixels varies with the PFTs: a very high homogeneity is possible for TrEBF (above 95%), and BoC3GRA (mostly above 80%). The homogeneity is usually above 70% for the selected pixels of TeEBF, TeC3GRA, BoC3GRA and C4GRA. Lower values can be obtained for the other PFTs (but still higher than 55%).

GPP data

For most PFTs, we assimilated daily *in situ* GPP estimates from FLUXNET (La Thuile (Baldocchi et al., 2001) and FLUXNET2015 (Pastorello et al., 2020)), based on the night-time method to partition net ecosystem exchange measurements between respiration and GPP components (Reichstein et al., 2005). The number of sites and their temporal cover strongly depend on the PFT. For instance, only four sites are considered for C4GRA and C4CRO PFTs while twenty-eight sites are used for TeDBF and TeENF.

Only for three PFTs (TrDBF, BoDNF, and TrC3GRA), there was no *in situ* GPP available. We therefore used FLUXSAT-GPP data (Joiner et al., 2018). Similarly to FLUXCOM, FLUXSAT also derives from a data-driven approach relying on FLUXNET measurements and MODIS reflectances in seven spectral bands and calibrated against FLUXNET measurements. One noteworthy difference is that FLUXSAT does not use any meteorological forcing. We used the data for the same 15 pixels as those considered for SIF, and data binned at 0.25°/weekly resolution (also similarly to SIF). Data from 2017-2020 were considered
5.1.4 Other data for evaluation and comparison

FluxCOM

The FLUXCOM GPP is estimated by upscaling data from eddy covariance sites globally via machine learning methods using satellite data such as MODIS and meteorological forcings (Tramontana et al., 2016d). We are using here a recent FLUXCOM dataset at 0.05°/monthly resolutions: label as v2 in the Figures of the product, available here: https://daac.ornl.gov/VEGETATION/guides/FluxSat_GPP_FPAR.html#datadescraccess.

FluxSat

Similarly to FLUXCOM, FLUXSAT (Joiner et al., 2018) also derives from a data- driven approach relying on FLUXNET measurements and MODIS reflectances in seven spectral bands and calibrated against FLUXNET measurements. One noteworthy difference is that FLUXSAT does not use any meteorological forcing.

TRENDY models

Global land surface models (LSMs), through the TRENDY model inter-comparison (Trends in net land carbon exchange; https://sites.exeter.ac.uk/trendy), provide independent global scale GPP estimates at 0.5° to 2° spatial resolutions and monthly temporal resolution. The GPP simulations from the most recent Global Carbon Budget - GCB (Friedlingstein et al., 2022) provide state-of-the-art information on the spatial and temporal GPP dynamics in response to climate drivers, land management, CO2 increase in the atmosphere, etc. The different TRENDY model simulations do not match perfectly the atmospheric CO2 growth rate but they provide a reasonable global land carbon sink. Friedlingstein et al. (2022) show that the TRENDY models provide over 2012-2021 a mean total land flux carbon sink of 1.5 +/- 0.5 GtC/yr. The land budget estimated as the residual from fossil and ocean fluxes is 1.5 +/- 0.6 GtC/yr. The TRENDY ensemble mean is thus similar to the one deduced from atmospheric CO2 but the spread of the TRENDY models (0.5 GtC/yr) indicates that some models strongly deviate from the mean.

5.2 Results

5.2.1 Model – data improvement

Figures 23 and 24 below illustrate the improvement of the model prediction after assimilation
with respect to GPP and SIF respectively by comparing the observed vs simulated mean
seasonal cycles, over all PFTs and pixels considered for the parameter optimization. Except
for grasses and crops, the prior GPPs simulated by ORCHIDEE agree well with the *in situ*
data, with an error (Root Mean Square Difference - RMSD) typically lower than 3 gC/m²/day.
Grasses and crops show a larger inter-pixels variability and a higher model-data mismatch. A
model improvement with respect to GPP following assimilation is observed for all PFTs (from
~1% only for TeENF for which the prior simulations were already in very good agreement with
the observation GPP data, to ~64% for BoDNF), except for boreal C3 grasses (BoC3GRA,
with a degradation of -0.3%). Comparatively, the prior SIF simulations largely overestimate
the TROPOMI SIF data. This initial model-data bias is largely decreased after assimilation
(RMSD reduction ranging from ~26% to ~80%), except for boreal deciduous broadleaf forest
(BoDBF). This PFT, which shows a slight model degradation (RMSD reduction of -0.9%), is
the only PFT for which the prior simulations underestimate the intensity of SIF as compared
to TROPOSIF data. Albeit a smaller error is associated with SIF data compared to GPP, we
can observe that the posterior model-data mismatch remains higher for SIF. In particular, the
optimized SIF simulations now underestimate the peak of the growing season for several PFTs
(TeENF, TeDBF, BoDBF, etc.) when the modelled GPP temporal dynamics are more
consistent with the observations. These different results obtained for SIF and GPP highlight
5.1.4 Other data for evaluation and comparison

FluxCOM

The FLUXCOM GPP is estimated by upscaling data from eddy covariance sites globally via machine learning methods using satellite data such as MODIS and meteorological forcings (Tramontana et al., 2016d). We are using here a recent FLUXCOM dataset at 0.05°/monthly resolutions: label as v2 in the Figures of the product, available here: https://daac.ornl.gov/VEGETATION/guides/FluxSat_GPP_FPAR.html#datadescraccess.

FluxSat

Similarly to FLUXCOM, FLUXSAT (Joiner et al., 2018) also derives from a data- driven approach relying on FLUXNET measurements and MODIS reflectances in seven spectral bands and calibrated against FLUXNET measurements. One noteworthy difference is that FLUXSAT does not use any meteorological forcing.

TRENDY models

Global land surface models (LSMs), through the TRENDY model inter-comparison (Trends in
net land carbon exchange; https://sites.exeter.ac.uk/trendy), provide independent global scale
GPP estimates at 0.5° to 2° spatial resolutions and monthly temporal resolution. The GPP
simulations from the most recent Global Carbon Budget - GCB (Friedlingstein et al., 2022)
provide state-of-the-art information on the spatial and temporal GPP dynamics in response to
climate drivers, land management, CO2 increase in the atmosphere, etc. The different
TRENDY model simulations do not match perfectly the atmospheric CO2 growth rate but they
provide a reasonable global land carbon sink. Friedlingstein et al. (2022) show that the
TRENDY models provide over 2012-2021 a mean total land flux carbon sink of 1.5 +/- 0.5
GtC/yr. The land budget estimated as the residual from fossil and ocean fluxes is 1.5 +/- 0.6
GtC/yr. The TRENDY ensemble mean is thus similar to the one deduced from atmospheric
CO2 but the spread of the TRENDY models (0.5 GtC/yr) indicates that some models strongly
deviate from the mean.

5.2 Results

5.2.1 Model – data improvement

Figures 23 and 24 below illustrate the improvement of the model prediction after assimilation
with respect to GPP and SIF respectively by comparing the observed vs simulated mean
seasonal cycles, over all PFTs and pixels considered for the parameter optimization. Except
for grasses and crops, the prior GPPs simulated by ORCHIDEE agree well with the *in situ*
data, with an error (Root Mean Square Difference - RMSD) typically lower than 3 gC/m²/day.
Grasses and crops show a larger inter-pixels variability and a higher model-data mismatch. A
model improvement with respect to GPP following assimilation is observed for all PFTs (from
~1% only for TeENF for which the prior simulations were already in very good agreement with
the observation GPP data, to ~64% for BoDNF), except for boreal C3 grasses (BoC3GRA,
with a degradation of -0.3%). Comparatively, the prior SIF simulations largely overestimate
the TROPOMI SIF data. This initial model-data bias is largely decreased after assimilation
(RMSD reduction ranging from ~26% to ~80%), except for boreal deciduous broadleaf forest
(BoDBF). This PFT, which shows a slight model degradation (RMSD reduction of -0.9%), is
the only PFT for which the prior simulations underestimate the intensity of SIF as compared
to TROPOSIF data. Albeit a smaller error is associated with SIF data compared to GPP, we
can observe that the posterior model-data mismatch remains higher for SIF. In particular, the
optimized SIF simulations now underestimate the peak of the growing season for several PFTs
(TeENF, TeDBF, BoDBF, etc.) when the modelled GPP temporal dynamics are more
consistent with the observations. These different results obtained for SIF and GPP highlight
some inconsistencies in our modelling framework that the data assimilation system fails to solve.

<figure><img src="image_89.png" /><figcaption>**Figure 23:** Comparison of the mean seasonal cycles of GPP for each PFT, between *in situ* data (black) and ORCHIDEE simulations (prior in blue and posterior in orange). Plain lines correspond to the mean values over the sites considered, and the shaded areas represent the corresponding standard deviation. Units are in gC/m²/day. The RMSD of the fit before and after assimilation are provided for each vegetation type, as well as the reduction in RMSD (in %). The PFTs indicated in green correspond to those for which FLUXSAT GPP data are used instead of *in situ* flux observations.</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>**Figure 24:** Same legend as figure 23 for SIF data (over 2018-2020).</figcaption></figure>
some inconsistencies in our modelling framework that the data assimilation system fails to solve.

<figure><img src="image_90.png" /><figcaption>Figure 23: Comparison of the mean seasonal cycles of GPP for each PFT, between *in situ* data (black) and ORCHIDEE simulations (prior in blue and posterior in orange). Plain lines correspond to the mean values over the sites considered, and the shaded areas represent the corresponding standard deviation. Units are in gC/m²/day. The RMSD of the fit before and after assimilation are provided for each vegetation type, as well as the reduction in RMSD (in %). The PFTs indicated in green correspond to those for which FLUXSAT GPP data are used instead of *in situ* flux observations.</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>Figure 24: Same legend as figure 23 for SIF data (over 2018-2020).</figcaption></figure>
Figure 25 summarizes the improvement of the model – data fit following the optimization.

<figure><img src="image_91.png" /><figcaption>Figure 25: Reduction in RMSD (in %) between ORCHIDEE simulations and SIF (left) and GPP (right) data following the data assimilation.</figcaption></figure>

As detailed above the improvement is greater for SIF that GPP but almost positive for all sites, except a few and mainly for the GPP.

5.2.2 Spatial distribution of GPP

Figure 26 provides the spatial GPP distribution for ORCHIDEE simulations, using prior and posterior to the data assimilation parameter values.

<figure><img src="image_92.png" /><figcaption>Figure 26: Yearly mean map for the simulations performed with the ORCHIDEE land surface model prior (left) and posterior (right) to data assimilation over the period 2018 (from May) - 2020. The global minimum, maximum, and mean values are provided (gC/m²/day), as well as the global budget (in PgC/m²/day).</figcaption></figure>

The simulations performed with the standard parameter values result in a mean global budget of 178 GtC/year, which is at the upper end of the estimates by other LSMs and data-driven products. The co-assimilation of SIF and GPP data decreases the global budget by about 35 GtC/year, resulting in a closer agreement with that of FLUXSAT or FLUXCOM-V2. The spatial distribution of the optimized GPP over the tropics is closer to that of FLUXSAT, than the one represented in FLUXCOM-V2. This can be partly explained by the fact that the constraint on the optimized model parameters relied on FLUXSAT estimates for the TrDBF (tropical deciduous broadleaf forest) PFT because no *in situ* data were available. However, in ORCHIDEE, this PFT is mostly dominant in the Northern and Southern parts of the African tropical forest, as well as in Northern Australia, and not over the Amazon basin (mostly tropical evergreen broadleaf forests) where large differences between FLUXSAT and FLUXCOM-V2 are observed.
Figure 25 summarizes the improvement of the model – data fit following the optimization.

<figure><img src="image_92.png" /><figcaption>Figure 25: Reduction in RMSD (in %) between ORCHIDEE simulations and SIF (left) and GPP (right) data following the data assimilation.</figcaption></figure>

As detailed above the improvement is greater for SIF that GPP but almost positive for all sites, except a few and mainly for the GPP.

5.2.2 Spatial distribution of GPP

Figure 26 provides the spatial GPP distribution for ORCHIDEE simulations, using prior and posterior to the data assimilation parameter values.

<figure><img src="image_93.png" /><figcaption>Figure 26: Yearly mean map for the simulations performed with the ORCHIDEE land surface model prior (left) and posterior (right) to data assimilation over the period 2018 (from May) - 2020. The global minimum, maximum, and mean values are provided (gC/m²/day), as well as the global budget (in PgC/m²/day).</figcaption></figure>

The simulations performed with the standard parameter values result in a mean global budget of 178 GtC/year, which is at the upper end of the estimates by other LSMs and data-driven products. The co-assimilation of SIF and GPP data decreases the global budget by about 35 GtC/year, resulting in a closer agreement with that of FLUXSAT or FLUXCOM-V2. The spatial distribution of the optimized GPP over the tropics is closer to that of FLUXSAT, than the one represented in FLUXCOM-V2. This can be partly explained by the fact that the constraint on the optimized model parameters relied on FLUXSAT estimates for the TrDBF (tropical deciduous broadleaf forest) PFT because no *in situ* data were available. However, in ORCHIDEE, this PFT is mostly dominant in the Northern and Southern parts of the African tropical forest, as well as in Northern Australia, and not over the Amazon basin (mostly tropical evergreen broadleaf forests) where large differences between FLUXSAT and FLUXCOM-V2 are observed.
## 5.3 Evaluation using other GPP products

Global GPP estimates are shown in Table 6 for the globe for ORCHIDEE prior and posterior to the optimization as well as for the FLUXCOM-V2 from 2002 to 2019.

<table><caption>Table 6: Global GPP estimations (in PgC/year) from the ORCHIDEE land surface model from 2002 to 2019, and from FLUXCOM-V2.</caption><thead><tr><th>Year</th><th>ORC prior</th><th>ORC posterior</th><th>FLUXCOM-V2</th></tr></thead><tbody><tr><td>2002</td><td>161.2</td><td>128.0</td><td>130.6</td></tr><tr><td>2003</td><td>163.4</td><td>129.4</td><td>130.6</td></tr><tr><td>2004</td><td>166.1</td><td>131.6</td><td>132.5</td></tr><tr><td>2005</td><td>164.0</td><td>129.9</td><td>132.2</td></tr><tr><td>2006</td><td>167.1</td><td>132.4</td><td>133.5</td></tr><tr><td>2007</td><td>168.0</td><td>133.1</td><td>133.3</td></tr><tr><td>2008</td><td>167.6</td><td>133.1</td><td>133.3</td></tr><tr><td>2009</td><td>166.5</td><td>132.1</td><td>133.9</td></tr><tr><td>2010</td><td>169.4</td><td>134.1</td><td>134.3</td></tr><tr><td>2011</td><td>171.7</td><td>136.5</td><td>135.3</td></tr><tr><td>2012</td><td>168.6</td><td>134.6</td><td>133.9</td></tr><tr><td>2013</td><td>169.3</td><td>134.8</td><td>134.7</td></tr><tr><td>2014</td><td>171.2</td><td>136.4</td><td>136.2</td></tr><tr><td>2015</td><td>170.2</td><td>136.0</td><td>134.6</td></tr><tr><td>2016</td><td>172.4</td><td>137.5</td><td>136.2</td></tr><tr><td>2017</td><td>173.2</td><td>138.5</td><td>136.8</td></tr><tr><td>2018</td><td>172.1</td><td>137.8</td><td>137.1</td></tr><tr><td>2019</td><td>171.6</td><td>137.6</td><td>137.3</td></tr><tr><td>Mean values</td><td>168.6</td><td>134.1</td><td>134.2</td></tr></tbody></table>

Model GPP estimates show a rapid increase of GPP during this period, of about 0.58 PgC/yr in ORCHIDEE prior and 0.53 PgC/yr after the optimization. However, it is significantly lower in the FLUXCOM-V2 product with 0.37 PgC/yr on average.

The assimilation of SIF data significantly reduces the GPP simulated by ORCHIDEE on average by 34 PgC/year (and slightly the trend). The mean GPP value following the optimisation, from 2002 to 2019, is 134.08 PgC/year, which is relatively close to the FLUXCOM value (136.27 PgC/yr). The reduction of the ORCHIDEE GPP occur over all regions but more significantly over the Tropics as shown by the spatial distribution of the ORCHIDEE GPP prior and posterior in Figure 26.

The optimization does not change too much the year-to-year GPP variations with a slight improvement of the correlation with FLUXCOM (see Table 6). It is however difficult with only the comparison to FLUXCOM to evaluate these interannual variations.
## 5.3 Evaluation using other GPP products

Global GPP estimates are shown in Table 6 for the globe for ORCHIDEE prior and posterior to the optimization as well as for the FLUXCOM-V2 from 2002 to 2019.

<table><caption>Table 6: Global GPP estimations (in PgC/year) from the ORCHIDEE land surface model from 2002 to 2019, and from FLUXCOM-V2.</caption><thead><tr><th>Year</th><th>ORC prior</th><th>ORC posterior</th><th>FLUXCOM-V2</th></tr></thead><tbody><tr><td>2002</td><td>161.2</td><td>128.0</td><td>130.6</td></tr><tr><td>2003</td><td>163.4</td><td>129.4</td><td>130.6</td></tr><tr><td>2004</td><td>166.1</td><td>131.6</td><td>132.5</td></tr><tr><td>2005</td><td>164.0</td><td>129.9</td><td>132.2</td></tr><tr><td>2006</td><td>167.1</td><td>132.4</td><td>133.5</td></tr><tr><td>2007</td><td>168.0</td><td>133.1</td><td>133.3</td></tr><tr><td>2008</td><td>167.6</td><td>133.1</td><td>133.3</td></tr><tr><td>2009</td><td>166.5</td><td>132.1</td><td>133.9</td></tr><tr><td>2010</td><td>169.4</td><td>134.1</td><td>134.3</td></tr><tr><td>2011</td><td>171.7</td><td>136.5</td><td>135.3</td></tr><tr><td>2012</td><td>168.6</td><td>134.6</td><td>133.9</td></tr><tr><td>2013</td><td>169.3</td><td>134.8</td><td>134.7</td></tr><tr><td>2014</td><td>171.2</td><td>136.4</td><td>136.2</td></tr><tr><td>2015</td><td>170.2</td><td>136.0</td><td>134.6</td></tr><tr><td>2016</td><td>172.4</td><td>137.5</td><td>136.2</td></tr><tr><td>2017</td><td>173.2</td><td>138.5</td><td>136.8</td></tr><tr><td>2018</td><td>172.1</td><td>137.8</td><td>137.1</td></tr><tr><td>2019</td><td>171.6</td><td>137.6</td><td>137.3</td></tr><tr><td>Mean values</td><td>168.6</td><td>134.1</td><td>134.2</td></tr></tbody></table>

Model GPP estimates show a rapid increase of GPP during this period, of about 0.58 PgC/yr in ORCHIDEE prior and 0.53 PgC/yr after the optimization. However, it is significantly lower in the FLUXCOM-V2 product with 0.37 PgC/yr on average.

The assimilation of SIF data significantly reduces the GPP simulated by ORCHIDEE on average by 34 PgC/year (and slightly the trend). The mean GPP value following the optimisation, from 2002 to 2019, is 134.08 PgC/year, which is relatively close to the FLUXCOM value (136.27 PgC/yr). The reduction of the ORCHIDEE GPP occur over all regions but more significantly over the Tropics as shown by the spatial distribution of the ORCHIDEE GPP prior and posterior in Figure 26.

The optimization does not change too much the year-to-year GPP variations with a slight improvement of the correlation with FLUXCOM (see Table 6). It is however difficult with only the comparison to FLUXCOM to evaluate these interannual variations.
The comparison is now performed for 3 latitudinal bands (Figures 27, 28 and 29) including the TRENDY ensemble of models and the two data-driven estimates, FLUXCOM and FLUXSAT.

<figure><img src="image_93.png" /><figcaption>Figure 27: Comparison of the average GPP time series for the simulations performed with ORCHIDEE (prior and optimized), the mean and standard deviation of the TRENDY simulations (black, grey shaded area) and the FluxCOM-V2 and FluxSat data products, for the northern hemisphere temperate latitudinal band (60°N / 30°N).</figcaption></figure>

<figure><img src="image_94.png" /><figcaption>Figure 28: Similar caption to Figure 27 for the northern tropics latitudinal band (30°N / 0°).</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Figure 29: Similar caption to Figure 27 for the southern tropics latitudinal band (0° / -30°N).</figcaption></figure>
The comparison is now performed for 3 latitudinal bands (Figures 27, 28 and 29) including the TRENDY ensemble of models and the two data-driven estimates, FLUXCOM and FLUXSAT.

<figure><img src="image_94.png" /><figcaption>Figure 27: Comparison of the average GPP time series for the simulations performed with ORCHIDEE (prior and optimized), the mean and standard deviation of the TRENDY simulations (black, grey shaded area) and the FluxCOM-V2 and FluxSat data products, for the northern hemisphere temperate latitudinal band (60°N / 30°N).</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Figure 28: Similar caption to Figure 27 for the northern tropics latitudinal band (30°N / 0°).</figcaption></figure>

<figure><img src="image_96.png" /><figcaption>Figure 29: Similar caption to Figure 27 for the southern tropics latitudinal band (0° / -30°N).</figcaption></figure>
5.4 Summary

The optimisation of the ORCHIDEE model with TROPOMI SIF and *in situ* GPP data, illustrates the observational constraint brought by the SIF data from Sentinel-5p. It provides large constraints on the ORCHIDEE model parameters. Overall, from this study and from an ensemble of tests not detailed in this report, we obtained:

* The parameter optimisation with TROPOMI SIF data and *in situ* GPP data leads to a large reduction of the simulated GPP by ORCHIDEE.

* The reduction of GPP concerns nearly all regions and it is associated to a strong reduction of the mean seasonal cycle. This improves the prior ORCHIDEE simulations which tend to overestimate the seasonal variability of GPP and NEE. Indeed prior ORCHIDEE simulations with many parameters prescribed from literature and trait databases are usually providing a too large seasonal cycle that is corrected through the assimilation of various data stream (including SIF data).

* Such GPP reduction brings the ORCHIDEE model close to the data-driven products FLUXCOM-V2 and FLUXSAT, especially in the northern temperate region.

* For the Tropics more substantial differences remain between the different GPP estimates.

* The use of *in situ* GPP data in the optimization procedure is important to stabilize the optimisation procedure and to avoid overfitting of SIF data and potential degradation of the GPP estimates (not shown here).

* Such an approach could be used in the context of the future CO2 monitoring system (or the Integrated Forecasting System of ECMWF) with near real time (NRT) data assimilation. In that case the assimilation of NRT SIF data would be complemented by the assimilation at the same time of a climatology of *in situ* GPP data (to stabilize the optimization procedure) with the target to optimize some state variable such as LAI and key parameters such as the maximum photosynthetic capacity (as well as some SIF related parameters)

* Such direction is currently being investigated within the Copernicus CAMS52a project directly with the CTESSEL model and the same SIF observation operator as the one used in ORCHIDEE.
5.4 Summary

The optimisation of the ORCHIDEE model with TROPOMI SIF and *in situ* GPP data, illustrates the observational constraint brought by the SIF data from Sentinel-5p. It provides large constraints on the ORCHIDEE model parameters. Overall, from this study and from an ensemble of tests not detailed in this report, we obtained:

* The parameter optimisation with TROPOMI SIF data and *in situ* GPP data leads to a large reduction of the simulated GPP by ORCHIDEE.

* The reduction of GPP concerns nearly all regions and it is associated to a strong reduction of the mean seasonal cycle. This improves the prior ORCHIDEE simulations which tend to overestimate the seasonal variability of GPP and NEE. Indeed prior ORCHIDEE simulations with many parameters prescribed from literature and trait databases are usually providing a too large seasonal cycle that is corrected through the assimilation of various data stream (including SIF data).

* Such GPP reduction brings the ORCHIDEE model close to the data-driven products FLUXCOM-V2 and FLUXSAT, especially in the northern temperate region.

* For the Tropics more substantial differences remain between the different GPP estimates.

* The use of *in situ* GPP data in the optimization procedure is important to stabilize the optimisation procedure and to avoid overfitting of SIF data and potential degradation of the GPP estimates (not shown here).

* Such an approach could be used in the context of the future CO2 monitoring system (or the Integrated Forecasting System of ECMWF) with near real time (NRT) data assimilation. In that case the assimilation of NRT SIF data would be complemented by the assimilation at the same time of a climatology of *in situ* GPP data (to stabilize the optimization procedure) with the target to optimize some state variable such as LAI and key parameters such as the maximum photosynthetic capacity (as well as some SIF related parameters)

* Such direction is currently being investigated within the Copernicus CAMS52a project directly with the CTESSEL model and the same SIF observation operator as the one used in ORCHIDEE.
# 6 SDBM

For the calculation of GPP with the Simple Diagnostic Biosphere Model (SDBM), ULUND and iLab use the Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) to constrain photosynthetic uptake of CO₂, instead of SIF as described in the previous sections. This option of using FAPAR instead of SIF has been explained before (e.g. in the WP3 section of the DMP, D9.5, as well as in the 2nd Periodic Report in form of a technical review report).

The reasons for doing so are the following:

* FAPAR products are more mature and the use of FAPAR within our modelling system directly yields a GPP estimate.

* FAPAR characterises radiation available for photosynthesis, in contrast to SIF, which quantifies the excess photon flux not used for photosynthesis.

* FAPAR has a clear model-independent physical definition as opposed to LAI, which in the model serves as a parameter, which does not necessarily coincide with a physical definition of LAI as leaf surface area.

* There are three EO data streams (VOD, SIF and LAI) listed in the CoCO₂ Description of Action. While VOD and SIF are used by other partners, to explore the constraint of all three EO data types (and hence the spread in GPP), we decided to use FAPAR as a substitute for LAI (see previous point).

We do not believe that this has a negative impacts on this task and the results, but rather the opposite, because within task 3.3 we can explore a larger observational database to constrain GPP (last point from the above list). The novelty then lies in the comparison of the GPP products derived from different EO datasets against the data-driven products and in situ measured fluxes.

## 6.1 Methods

The SDBM version used here is conceptually similar to the one used by Kaminski et al. (2017) and is based on Knorr and Heimann (1995). The spatial resolution of the model is 0.5° times 0.5° globally. The model classifies global terrestrial ecosystems into 8 biome types based on the MODIS land cover classification (Friedl et al., 2010): (1) cropland/urban/natural vegetation mosaic, (2) needleleaf forest, (3) broadleaf forest, (4) mixed forest, (5) shrubland, (6) savanna or grassland, (7) tundra, (8) barren or sparsely vegetated.

The terrestrial uptake of CO₂ by photosynthesis (gross primary productivity, GPP) is calculated using a light-use efficiency approach:

$$GPP(x, t) = \varepsilon \cdot \alpha(x, t) \cdot FAPAR(x, t) \cdot Srad(x, t) \qquad (6.1)$$

where ε denotes the biome-specific light use efficiency parameter, α a plant water stress factor, and Srad the incident solar radiation driving photosynthesis. Autotrophic respiration (Rₐ) is calculated from

$$R_a(x, t) = r \cdot GPP(x, t), \qquad (6.2)$$

with r being a constant with a value of r = 1 - 1/2.13 (Waring et al., 1998). Ecosystem heterotrophic respiration (Rₕ) is calculated following a Q₁₀ functional relationship with temperature T and is, as the photosynthesis calculation, modulated by the water stress factor:

$$R_h(x, t) = n(x) \cdot \alpha(x, t) \cdot Q_{10}^{T/10} \qquad (6.3)$$

The biome-specific parameter Q₁₀ expresses the ratio of respiration at ambient temperature T + 10 to that at T, with T given in °C. The spatially varying normalisation factor n(x) is the ratio of the temporal integrals of NPP(x,t) = GPP(x,t) - Rₐ(x,t) and Rₕ(x,t) computed for prior