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

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

## Supplement of

# Evaluating the performance of land surface model ORCHIDEE-CAN v1.0 on water and energy flux estimation with a single- and multi-layer energy budget scheme

Yiying Chen et al.

Correspondence to: Yiying Chen (yiyingchen@gate.sinica.edu.tw)

The copyright of individual parts of the supplement might differ from the CC-BY 3.0 licence.
# Evaluating the performance of land surface model ORCHIDEE-CAN v1.0 on water and energy flux estimation with a single- and multi-layer energy budget scheme

Yiying Chen, James Ryder, Vladislav Bastrikov, Matthew J. Mcgrath, Kim Naudts, Juliane Otto, Catherine Ottle, Philippe Peylin, Jan Polcher, Aude Valade, et al.

► To cite this version:

Yiying Chen, James Ryder, Vladislav Bastrikov, Matthew J. Mcgrath, Kim Naudts, et al. Evaluating the performance of land surface model ORCHIDEE-CAN v1.0 on water and energy flux estimation with a single- and multi-layer energy budget scheme. Geoscientific Model Development, 2016, 9 (9), pp.2951 - 2972. 10.5194/gmd-9-2951-2016. hal-01587561

HAL Id: hal-01587561

https://hal.science/hal-01587561

Submitted on 27 May 2020

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

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

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

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

**Figure S1.** Model simulation and observation of the wind speed profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimized parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.

**Table S1.** Description of the experimental design. The model was forced either by the site-level observations (SITE) or the CRU-NCEP re-analysis (CRU) and was run with the single-layer energy budget scheme (SINGLE) or the multi-layer energy budget scheme (MULTI). The model could be forced to follow the observed *LAI* profiles (IMPOSE) or made use of the internal calculation of the seasonal dynamics and vertical profile of *LAI* (SIM). EXP denotes the experiment name, PERIOD refers to the periods for which the simulations were run as defined in Table 3.

<table><thead><tr><th rowspan="2">EXP</th><th colspan="2">FORCING</th><th colspan="2">ENERGY BUDGET</th><th colspan="2">LAI PROFILE</th><th rowspan="2">PERIOD</th></tr><tr><th>SITE</th><th>CRU</th><th>SINGLE</th><th>MULTI</th><th>IMPOSE</th><th>SIM</th></tr></thead><tbody><tr><td>SPINUP</td><td></td><td>+</td><td>+</td><td></td><td></td><td>+</td><td>20yrs</td></tr><tr><td>optimizE</td><td>+</td><td></td><td>+</td><td></td><td>+</td><td></td><td>I & II</td></tr><tr><td>EXP1</td><td>+</td><td></td><td>+</td><td></td><td>+</td><td></td><td>III</td></tr><tr><td>EXP2</td><td>+</td><td></td><td>+</td><td></td><td>+</td><td></td><td>IV</td></tr><tr><td>EXP3</td><td>+</td><td></td><td></td><td>+</td><td>+</td><td></td><td>III</td></tr><tr><td>EXP4</td><td>+</td><td></td><td></td><td>+</td><td>+</td><td></td><td>IV</td></tr></tbody></table>

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

# Evaluating the performance of land surface model ORCHIDEE-CAN v1.0 on water and energy flux estimation with a single- and multi-layer energy budget scheme

Yiying Chen<sup>1,a</sup>, James Ryder<sup>1</sup>, Vladislav Bastrikov<sup>1</sup>, Matthew J. McGrath<sup>1</sup>, Kim Naudts<sup>1,b</sup>, Juliane Otto<sup>1,c</sup>, Catherine Ottlé<sup>1</sup>, Philippe Peylin<sup>1</sup>, Jan Polcher<sup>2</sup>, Aude Valade<sup>3</sup>, Andrew Black<sup>4</sup>, Jan A. Elbers<sup>5</sup>, Eddy Moors<sup>5</sup>, Thomas Foken<sup>6</sup>, Eva van Gorsel<sup>7</sup>, Vanessa Haverd<sup>7</sup>, Bernard Heinesch<sup>8</sup>, Frank Tiedemann<sup>9</sup>, Alexander Knohl<sup>9</sup>, Samuli Launiainen<sup>10</sup>, Denis Loustau<sup>11</sup>, Jérôme Ogée<sup>11</sup>, Timo Vessala<sup>12,13</sup>, and Sebastiaan Luyssaert<sup>1,d</sup>

<sup>1</sup>Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France

<sup>2</sup>Laboratoire de Météorologie Dynamique (LMD, CNRS), Ecole Polytechnique, Palaiseau, France

<sup>3</sup>Institut Pierre Simon Laplace, Place Jussieu 4, 75010 Paris, France

<sup>4</sup>Land and Food Systems, University of British Columbia, Vancouver, BC, Canada

<sup>5</sup>Alterra, Wageningen UR, Wageningen, the Netherlands

<sup>6</sup>Department of Micrometeorology, University of Bayreuth, Bayreuth Center of Ecology and Environmental Research, Bayreuth, Germany

<sup>7</sup>CSIRO, Marine and Atmospheric Research, Canberra, Australia

<sup>8</sup>Dept. Biosystem Engineering (BIOSE), University of Liege, Gembloux, Belgium

<sup>9</sup>Dept. Bioclimatology, Georg-August University of Göttingen, Büsgenweg, Göttingen, Germany

<sup>10</sup>Natural Resources Institute Finland, Vantaa, Finland

<sup>11</sup>INRA UMR 1391 ISPA Centre de Bordeaux Aquitaine, Bordeaux, France

<sup>12</sup>Department of Physics, University of Helsinki, Helsinki, Finland

<sup>13</sup>Department of Forest Sciences, University of Helsinki, Helsinki, Finland

<sup>a</sup>now at: Research Center for Environmental Changes (RCEC), Academia Sinica, Taipei, Taiwan

<sup>b</sup>now at: Department of Land in the Earth System, Max Planck Institute for Meteorology, Hamburg, Germany

<sup>c</sup>now at: Climate Service Center Germany (GERICS), Helmholtz-Zentrum Geesthacht, Hamburg, Germany

<sup>d</sup>now at: Department of Ecological Sciences, VU University, Amsterdam, the Netherlands

*Correspondence to:* Yiying Chen (yiyingchen@gate.sinica.edu.tw)

Received: 2 February 2016 – Published in Geosci. Model Dev. Discuss.: 25 February 2016

Revised: 5 August 2016 – Accepted: 8 August 2016 – Published: 2 September 2016

**Abstract.** Canopy structure is one of the most important vegetation characteristics for land-atmosphere interactions, as it determines the energy and scalar exchanges between the land surface and the overlying air mass. In this study we evaluated the performance of a newly developed multi-layer energy budget in the ORCHIDEE-CAN v1.0 land surface model (Organising Carbon and Hydrology In Dynamic Ecosystems – CANopy), which simulates canopy structure and can be coupled to an atmospheric model using an implicit coupling procedure. We aim to provide a set of accept-

able parameter values for a range of forest types. Top-canopy and sub-canopy flux observations from eight sites were collected in order to conduct this evaluation. The sites crossed climate zones from temperate to boreal and the vegetation types included deciduous, evergreen broad-leaved and evergreen needle-leaved forest with a maximum leaf area index (LAI; all-sided) ranging from 3.5 to 7.0. The parametrization approach proposed in this study was based on three selected physical processes – namely the diffusion, advection, and turbulent mixing within the canopy. Short-term sub-canopy
<table><caption>**Table S2.** Optimized parameter values per site. The uncertainties (1 standard deviation) were derived from the sensitivity analysis for the soil water content at the end of the spin-up.</caption><thead><tr><th>Site</th><th>FL-Hyy</th><th>FR-LBr</th><th>NL-Loo</th><th>DE-Bay</th><th>CA-Oas</th><th>AU-Tum</th><th>DE-Hai</th><th>BE-Vie</th></tr></thead><tbody><tr><td>Code</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>a<sub>3</sub></td><td>0.420(±0.0038)</td><td>0.300(±0.0027)</td><td>0.302(±0.0027)</td><td>0.387(±0.0035)</td><td>0.234(±0.0021)</td><td>0.360(±0.0032)</td><td>0.301(±0.0027)</td><td>0.341(±0.0031)</td></tr><tr><td>a<sub>4</sub></td><td>-0.374(±0.0041)</td><td>-0.098(±0.0011)</td><td>-0.111(±0.0012)</td><td>-0.306(±0.0034)</td><td>-0.051(±0.0006)</td><td>-0.081(±0.0009)</td><td>-0.400(±0.0044)</td><td>-0.223(±0.0025)</td></tr><tr><td>a<sub>5</sub></td><td>0.050(±0.0010)</td><td>0.050(±0.0010)</td><td>0.085(±0.0017)</td><td>0.006(±0.0001)</td><td>0.079(±0.0016)</td><td>0.028(±0.0006)</td><td>0.059(±0.0012)</td><td>0.086(±0.0017)</td></tr><tr><td>a<sub>6</sub></td><td>16.82(±0.0841)</td><td>11.52(±0.0576)</td><td>11.29(±0.0565)</td><td>19.21(±0.0961)</td><td>10.56(±0.0528)</td><td>20.10(±0.1005)</td><td>10.01(±0.0501)</td><td>11.00(±0.0550)</td></tr><tr><td>a<sub>7</sub></td><td>0.06(±0.0005)</td><td>0.32(±0.0026)</td><td>0.18(±0.0014)</td><td>0.11(±0.0009)</td><td>0.21(±0.0017)</td><td>0.40(±0.0032)</td><td>0.13(±0.0010)</td><td>0.05(±0.0004)</td></tr><tr><td>a<sub>8</sub></td><td>4.57(±0.0914)</td><td>4.82(±0.0964)</td><td>1.71(±0.0342)</td><td>5.10(±0.1200)</td><td>6.53(±0.1360)</td><td>1.50(±0.0300)</td><td>5.20(±0.1400)</td><td>4.70(±0.0940)</td></tr><tr><td>a<sub>9</sub></td><td>0.52(±0.0026)</td><td>0.45(±0.0015)</td><td>0.77(±0.0022)</td><td>0.56(±0.0015)</td><td>0.57(±0.0029)</td><td>0.62(±0.0031)</td><td>0.46(±0.0023)</td><td>0.53(±0.0027)</td></tr><tr><td>a<sub>10</sub></td><td>0.99(±0.0198)</td><td>0.95(±0.0180)</td><td>0.52(±0.0102)</td><td>0.93(±0.0186)</td><td>0.95(±0.0190)</td><td>1.60(±0.0320)</td><td>0.97(±0.0194)</td><td>0.95(±0.0190)</td></tr><tr><td>W<sub>br</sub></td><td>0.81(±0.0353)</td><td>2.63(±0.1147)</td><td>1.83(±0.0798)</td><td>7.57(±0.3301)</td><td>3.20(±0.1395)</td><td>0.86(±0.0375)</td><td>7.56(±0.3296)</td><td>4.53(±0.1975)</td></tr><tr><td>W<sub>sr</sub></td><td>2.97(±0.0624)</td><td>1.88(±0.0395)</td><td>5.53(±0.1161)</td><td>2.87(±0.0603)</td><td>6.70(±0.1407)</td><td>2.43(±0.0510)</td><td>4.27(±0.0897)</td><td>4.35(±0.0914)</td></tr></tbody></table>
observations and long-term surface fluxes were used to cali-
brate the parameters in the sub-canopy radiation, turbulence,
and resistance modules with an automatic tuning process.
The multi-layer model was found to capture the dynamics
of sub-canopy turbulence, temperature, and energy fluxes.
The performance of the new multi-layer model was further
compared against the existing single-layer model. Although
the multi-layer model simulation results showed few or no
improvements to both the nighttime energy balance and en-
ergy partitioning during winter compared with a single-layer
model simulation, the increased model complexity does pro-
vide a more detailed description of the canopy micrometeo-
rology of various forest types. The multi-layer model links to
potential future environmental and ecological studies such as
the assessment of in-canopy species vulnerability to climate
change, the climate effects of disturbance intensities and fre-
quencies, and the consequences of biogenic volatile organic
compound (BVOC) emissions from the terrestrial ecosystem.

# 1 Introduction

Today's Earth system models (ESMs) integrate ocean, ice
sheet, atmosphere, and land surface in order to provide a
powerful tool to simulate the Earth's past, present, and fu-
ture climates (Drobinski et al., 2012). In such a model, the
land surface sub-model provides the surface fluxes to the at-
mospheric sub-model, affects the dynamics of the planetary
boundary layer, and exerts a strong influence on the climate.
The dynamics of the simulated surface fluxes rely on the land
surface sub-model that, over the past 40 years, has evolved
from a simple bucket model approach towards sophisticated
soil-vegetation-atmosphere transfer (SVAT) schemes (Pit-
man, 2003; Stöckli and Vidale, 2005).

Although present-day land surface models differ from
each other in their formulation and details, their perfor-
mances show similar deficiencies. For example, imposing the
same land cover changes to seven land surface models re-
sulted in diverging climate effects. Amongst other factors,
this divergence was due to the parametrization of albedo and
the representation of evapotranspiration for different land
cover types (Pitman et al., 2009). Difficulties in reproduc-
ing fluxes of sensible and latent heat for a wide range of
vegetation types have been ascribed to the so-called “big-
leaf” approach (Bonan, 1996; Sellers et al., 1996; Dickinson
et al., 1998; Jiménez et al., 2011) which treats the surface as
a isothermal large leaf. Potentially, representing the vertical
canopy structure in detail and simulating radiation partition-
ing and turbulent transport within the vegetation will result
in an improved determination of sensible and latent heat flux
estimates (Baldocchi and Wilson, 2001; Ogée et al., 2003;
Bonan et al., 2014). For example, several multi-layer SVAT
schemes have been proposed and validated with site-level
observations (Ogée et al., 2003; Staudt et al., 2011; Haverd

et al., 2012; Launiainen et al., 2015). These studies demon-
strated that both top-canopy and within-canopy fluxes and
micrometeorological profiles could be captured by means of
a sophisticated parametrization scheme to describe the veg-
etation dynamics and the coupling between the atmosphere
and the canopy.

Because the standard version of ORCHIDEE (Organis-
ing Carbon and Hydrology In Dynamic Ecosystems) makes
use of a big-leaf approach (Ducoudré et al., 1993; Krin-
ner et al., 2005), improved model capacity and performance
were aimed for by implementation of a multi-layer energy
budget scheme (Ryder et al., 2016) that was integrated with
vertically discrete reflectivity, photosynthesis, stomatal resis-
tance, and carbon allocation schemes. This new design re-
sulted in a new version of ORCHIDEE named ORCHIDEE-
CAN (ORCHIDEE-CANopy, revision 2290) (Naudts et al.,
2015). Despite its code including a multi-layer energy bud-
get scheme (Ryder et al., 2016), ORCHIDEE-CAN is cur-
rently applied using a single-layer energy budget, due to a
lack of validated parameters for the multi-layer energy bud-
get scheme.

In Ryder et al. (2016), the model was developed and
tested for a single site. In this study, we compiled a set of
within-canopy and above-canopy measurements of energy,
water, and CO₂ fluxes and used these data to parametrize
and validate the new multi-layer energy budget scheme in
the ORCHIDEE-CAN v1.0 (revision 2754) global-scale land
surface model. The data set allowed us to test the model
under diverse environmental conditions in order to demon-
strate that the numerics can deal with the variation that can
be found in global ecosystems. For this we granted ourselves
the freedom to derive a separate parameter set for each site.
Model performance of the new multi-layer parametrization
was compared against the existing single-layer model. By do-
ing so we learned about the strengths and weaknesses of the
model and its parameters. In subsequent studies, we will have
to derive a single parameter set for each plant functional type
(PFT) and test how well the model reproduces global patterns
in, for example, evapotranspiration.

## 2 Methodology

### 2.1 Multi-layer energy budget scheme

The multi-layer energy budget scheme used in this study was developed for global land surface models (Ryder et al., 2016) and the calculations differ from the more common big-leaf energy budget scheme in three aspects. The new scheme calculates the following: (a) within-canopy longwave and shortwave radiation based on a vertical leaf area index (LAI; m² m⁻²) profile, (b) a within-canopy and below-canopy wind profile based on the vertical LAI profile, and (c) the dependency of stomatal resistance and aerodynamic resistance based on the microclimatological conditions along
<table><caption>Table S3. Calibration results during observation Period I and II for each site.</caption><thead><tr><th rowspan="2">Site Code</th><th colspan="3">Period I</th><th colspan="2">Period II</th></tr><tr><th>optimized variable</th><th>RMSE prior(default)</th><th>RMSE optimized</th><th>RMSE prior(default)</th><th>RMSE optimized</th></tr></thead><tbody><tr><td rowspan="6">AU-Tum</td><td>R<sub>n</sub></td><td>51.4</td><td>51.9</td><td></td><td></td></tr><tr><td>LE</td><td>86.6</td><td>38.9</td><td>39.5</td><td>42.4</td></tr><tr><td>H</td><td>150.9</td><td>33.1</td><td>46.3</td><td>36.4</td></tr><tr><td>U</td><td>0.15</td><td>0.07</td><td></td><td></td></tr><tr><td>T<sub>a</sub></td><td>0.48</td><td>0.35</td><td></td><td></td></tr><tr><td>q<sub>a</sub></td><td>0.00030</td><td>0.00027</td><td></td><td></td></tr><tr><td rowspan="6">BE-Vie</td><td>R<sub>n</sub></td><td>32.9</td><td>39.6</td><td></td><td></td></tr><tr><td>LE</td><td>102.6</td><td>38.1</td><td>125.8</td><td>23.5</td></tr><tr><td>H</td><td>97.3</td><td>44.8</td><td>127.7</td><td>31.3</td></tr><tr><td>U</td><td>0.64</td><td>0.64</td><td></td><td></td></tr><tr><td>T<sub>a</sub></td><td>0.61</td><td>0.86</td><td></td><td></td></tr><tr><td>q<sub>a</sub></td><td>0.00087</td><td>0.00083</td><td></td><td></td></tr><tr><td rowspan="6">CA-Oas</td><td>R<sub>n</sub></td><td>35.1</td><td>34.1</td><td></td><td></td></tr><tr><td>LE</td><td>54.0</td><td>34.7</td><td>150.9</td><td>62.4</td></tr><tr><td>H</td><td>73.9</td><td>50.2</td><td>155.3</td><td>72.2</td></tr><tr><td>U</td><td>0.25</td><td>0.21</td><td></td><td></td></tr><tr><td>T<sub>a</sub></td><td>1.27</td><td>1.24</td><td></td><td></td></tr><tr><td>q<sub>a</sub></td><td>n.a.</td><td>n.a.</td><td></td><td></td></tr><tr><td rowspan="6">DE-Bay</td><td>R<sub>n</sub></td><td>33.3</td><td>33.3</td><td></td><td></td></tr><tr><td>LE</td><td>76.3</td><td>74.7</td><td>128.1</td><td>23.4</td></tr><tr><td>H</td><td>60.7</td><td>30.2</td><td>136.6</td><td>34.2</td></tr><tr><td>U</td><td>0.62</td><td>0.21</td><td></td><td></td></tr><tr><td>T<sub>a</sub></td><td>0.82</td><td>0.64</td><td></td><td></td></tr><tr><td>q<sub>a</sub></td><td>n.a.</td><td>n.a.</td><td></td><td></td></tr><tr><td rowspan="6">DE-Hai</td><td>R<sub>n</sub></td><td>21.0</td><td>24.7</td><td></td><td></td></tr><tr><td>LE</td><td>138.6</td><td>35.7</td><td>87.4</td><td>32.3</td></tr><tr><td>H</td><td>148.9</td><td>48.9</td><td>88.2</td><td>43.5</td></tr><tr><td>U</td><td>2.05</td><td>1.21</td><td></td><td></td></tr><tr><td>T<sub>a</sub></td><td>0.78</td><td>0.79</td><td></td><td></td></tr><tr><td>q<sub>a</sub></td><td>n.a.</td><td>n.a.</td><td></td><td></td></tr></tbody></table>
<table><caption>**Table 1.** Symbolic notation used throughout the paper.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>a<sub>1</sub>, a<sub>2</sub>, a<sub>3</sub>, a<sub>4</sub>, a<sub>5</sub></td><td>tuning coefficients for C<sub>Deff</sub></td><td>unitless</td></tr><tr><td>a<sub>6</sub></td><td>factor ceiling of the slope</td><td>unitless</td></tr><tr><td>a<sub>7</sub></td><td>critical friction velocity in the middle point of the S-shape function</td><td>unitless</td></tr><tr><td>a<sub>8</sub></td><td>factor to constrain the S-shape function</td><td>unitless</td></tr><tr><td>a<sub>9</sub></td><td>threshold for vegetation cover</td><td>unitless</td></tr><tr><td>a<sub>10</sub></td><td>linear weighting factor</td><td>unitless</td></tr><tr><td>A</td><td>assimilation rate</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>C<sub>Deff</sub></td><td>effective drag coefficient</td><td>unitless</td></tr><tr><td>C<sub>S</sub></td><td>concentration of CO<sub>2</sub> at leaf surface</td><td>ppm</td></tr><tr><td>C<sub>D,i</sub></td><td>vertically discretized estimate for canopy drag coefficient</td><td>unitless</td></tr><tr><td>D<sub>h,air</sub></td><td>heat diffusivity of air</td><td>cm<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>D<sub>h,H<sub>2</sub>O</sub></td><td>heat diffusivity of water vapour</td><td>cm<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>d<sub>l</sub></td><td>characteristic leaf length</td><td>m</td></tr><tr><td>f<sub>P gap</sub></td><td>over-story gap probability from P gap fraction</td><td>m<sup>2</sup> m<sup>-2</sup></td></tr><tr><td>G<sub>veg</sub></td><td>logic variable to indicate the growth status of the vegetation</td><td>unitless</td></tr><tr><td>g<sub>0</sub></td><td>residual stomatal conductance if the irradiance approaches zero</td><td>m s<sup>-1</sup></td></tr><tr><td>h<sub>s</sub></td><td>relative humidity at leaf surface</td><td>%</td></tr><tr><td>h<sub>c</sub></td><td>canopy height</td><td>m</td></tr><tr><td>k<sub>i</sub></td><td>diffusivity for level i</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>k<sup>*</sup><sub>i</sub></td><td>modified diffusivity for level i</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>k<sub>surf</sub></td><td>conductance for the surface-atmosphere interface</td><td>m s<sup>-1</sup></td></tr><tr><td>LAI<sub>i</sub></td><td>leaf area index at level i</td><td>m<sup>2</sup> m<sup>-2</sup></td></tr><tr><td>Nu</td><td>Nusselt number</td><td>unitless</td></tr><tr><td>P<sub>m,i</sub></td><td>momentum shielding factor</td><td>unitless</td></tr><tr><td>PAI</td><td>plant area index</td><td>m<sup>2</sup> m<sup>-2</sup></td></tr><tr><td>R</td><td>correlation coefficient between the simulation and the observation</td><td>unitless</td></tr><tr><td>R<sub>0</sub></td><td>maximum correlation coefficient</td><td>unitless</td></tr><tr><td>R<sub>b,i</sub></td><td>boundary layer resistance at level i for heat</td><td>s m<sup>-1</sup></td></tr><tr><td>R'<sub>b,i</sub></td><td>boundary layer resistance at level i for water vapour</td><td>s m<sup>-1</sup></td></tr><tr><td>R<sub>s,i</sub></td><td>stomatal resistance at level i</td><td>s m<sup>-1</sup></td></tr><tr><td>Re</td><td>Reynold's number</td><td>unitless</td></tr><tr><td>SLA</td><td>specific leaf area</td><td>m<sup>2</sup> g<sup>-1</sup></td></tr><tr><td>S<sub>T</sub></td><td>Taylor skill score</td><td>unitless</td></tr><tr><td>T<sub>week</sub></td><td>weekly mean air temperature</td><td>K</td></tr><tr><td>T<sub>g</sub></td><td>temperature threshold for under-story phenology</td><td>K</td></tr><tr><td>T<sub>L</sub></td><td>Lagrangian timescale</td><td>s</td></tr><tr><td>u<sub>*</sub></td><td>friction velocity</td><td>m s<sup>-1</sup></td></tr><tr><td>u<sub>i</sub></td><td>velocity at level i</td><td>m s<sup>-1</sup></td></tr><tr><td>V<sub>cmax</sub></td><td>carboxylation capacity</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>W<sub>br</sub></td><td>weighting parameter for boundary layer resistance</td><td>unitless</td></tr><tr><td>W<sub>nf</sub></td><td>near-field weighting factor</td><td>unitless</td></tr><tr><td>W<sub>sf</sub></td><td>weighting parameter for soil-atmosphere conductance</td><td>unitless</td></tr><tr><td>W<sub>sr</sub></td><td>linear reduction parameter for stomatal resistance</td><td>unitless</td></tr><tr><td>β<sub>3</sub></td><td>fraction of potential plant transpiration realized</td><td>unitless</td></tr><tr><td>β<sub>4</sub></td><td>fraction of soil evaporation realized</td><td>unitless</td></tr><tr><td>μ</td><td>kinematic viscosity of air</td><td>cm<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>σ<sub>f</sub></td><td>ratio of the variance of the simulations over the variances of observations</td><td>unitless</td></tr><tr><td>σ<sub>w</sub></td><td>standard deviation in vertical velocity</td><td>m s<sup>-1</sup></td></tr></tbody></table>

the LAI profile. All symbols are explained in Table 1. In the following paragraphs these calculations are further described.

a. The multi-layer energy budget scheme makes use of the longwave radiation transfer scheme proposed by Gao

et al. (1989) and Gu et al. (1999). The scheme simulates longwave radiation transport, as well as scattering and absorption, along a vertically layered leaf area distribution. The simulated longwave radiation within a layer depends on the emitted longwave radiation by all of its
<table><caption>Table S3. Continuation of Table ??</caption><thead><tr><th rowspan="2">Site Code</th><th colspan="3">Period I</th><th colspan="2">Period II</th></tr><tr><th>optimize variable</th><th>RMSE prior(default)</th><th>RMSE optimized</th><th>RMSE prior(default)</th><th>RMSE optimized</th></tr></thead><tbody><tr><td rowspan="6">FI-Hyy</td><td><i>R</i><sub>n</sub></td><td>33.5</td><td>33.0</td><td></td><td></td></tr><tr><td><i>LE</i></td><td>157.9</td><td>49.3</td><td>44.5</td><td>21.2</td></tr><tr><td><i>H</i></td><td>155.5</td><td>52.5</td><td>46.9</td><td>32.3</td></tr><tr><td><i>U</i></td><td>0.23</td><td>0.15</td><td></td><td></td></tr><tr><td><i>T</i><sub>a</sub></td><td>1.15</td><td>1.14</td><td></td><td></td></tr><tr><td><i>q</i><sub>a</sub></td><td>0.00024</td><td>0.00015</td><td></td><td></td></tr><tr><td rowspan="6">FR-LBr</td><td><i>R</i><sub>n</sub></td><td>27.4</td><td>25.6</td><td></td><td></td></tr><tr><td><i>LE</i></td><td>89.4</td><td>49.5</td><td>44.5</td><td>40.4</td></tr><tr><td><i>H</i></td><td>73.4</td><td>47.3</td><td>51.7</td><td>32.8</td></tr><tr><td><i>U</i></td><td>0.17</td><td>0.15</td><td></td><td></td></tr><tr><td><i>T</i><sub>a</sub></td><td>1.46</td><td>1.46</td><td></td><td></td></tr><tr><td><i>q</i><sub>a</sub></td><td>0.00037</td><td>0.00038</td><td></td><td></td></tr><tr><td rowspan="6">NL-Loo</td><td><i>R</i><sub>n</sub></td><td>33.6</td><td>33.4</td><td></td><td></td></tr><tr><td><i>LE</i></td><td>71.2</td><td>47.9</td><td>63.2</td><td>22.1</td></tr><tr><td><i>H</i></td><td>122.4</td><td>56.9</td><td>63.9</td><td>33.3</td></tr><tr><td><i>U</i></td><td>0.88</td><td>0.75</td><td></td><td></td></tr><tr><td><i>T</i><sub>a</sub></td><td>0.81</td><td>0.78</td><td></td><td></td></tr><tr><td><i>q</i><sub>a</sub></td><td>0.00072</td><td>0.00067</td><td></td><td></td></tr><tr><td rowspan="6">All Sites</td><td><i>R</i><sub>n</sub></td><td>33.5</td><td>34.5</td><td></td><td></td></tr><tr><td><i>LE</i></td><td>91.2</td><td>46.1</td><td>85.5</td><td>38.2</td></tr><tr><td><i>H</i></td><td>123.2</td><td>50.3</td><td>89.6</td><td>40.4</td></tr><tr><td><i>U</i></td><td>0.62</td><td>0.42</td><td></td><td></td></tr><tr><td><i>T</i><sub>a</sub></td><td>0.92</td><td>0.93</td><td></td><td></td></tr><tr><td><i>q</i><sub>a</sub></td><td>0.00047</td><td>0.00043</td><td></td><td></td></tr></tbody></table>

<table><caption>Table S4. Evaluation of the model performance, Taylor score ($S_T$), correlation coefficient (R) and root mean square error (RMSE) for four experiments and changes in performance.</caption><thead><tr><th>Experiment</th><th>EXP1</th><th>EXP2</th><th>|EXP1-EXP2|</th><th>EXP3</th><th>EXP4</th><th>|EXP3-EXP4|</th></tr></thead><tbody><tr><td><strong>Rn</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><i>S</i><sub>T</sub> (0 – 1)</td><td>0.961</td><td>0.931</td><td>0.030</td><td>0.893</td><td>0.924</td><td>0.031</td></tr><tr><td><i>R</i> (0 – 1)</td><td>0.986</td><td>0.874</td><td></td><td>0.763</td><td>0.903</td><td></td></tr><tr><td>RMSE (Wm<sup>-2</sup>)</td><td>33.21</td><td>87.30</td><td></td><td>113.1</td><td>64.31</td><td></td></tr><tr><td><strong>H</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><i>S</i><sub>T</sub> (0 – 1)</td><td>0.863</td><td>0.828</td><td>0.035</td><td>0.810</td><td>0.865</td><td>0.054</td></tr><tr><td><i>R</i> (0 – 1)</td><td>0.777</td><td>0.689</td><td></td><td>0.774</td><td>0.788</td><td></td></tr><tr><td>RMSE (Wm<sup>-2</sup>)</td><td>59.64</td><td>71.51</td><td></td><td>45.88</td><td>42.15</td><td></td></tr><tr><td><strong>LE</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><i>S</i><sub>T</sub> (0 – 1)</td><td>0.822</td><td>0.778</td><td>0.044</td><td>0.786</td><td>0.745</td><td>0.041</td></tr><tr><td><i>R</i> (0 – 1)</td><td>0.804</td><td>0.710</td><td></td><td>0.649</td><td>0.645</td><td></td></tr><tr><td>RMSE (Wm<sup>-2</sup>)</td><td>48.06</td><td>56.44</td><td></td><td>51.64</td><td>41.01</td><td></td></tr><tr><td><strong>G</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><i>S</i><sub>T</sub> (0 – 1)</td><td>0.234</td><td>0.275</td><td>0.041</td><td>0.410</td><td>0.454</td><td>0.044</td></tr><tr><td><i>R</i> (0 – 1)</td><td>0.544</td><td>0.451</td><td></td><td>0.424</td><td>0.507</td><td></td></tr><tr><td>RMSE (Wm<sup>-2</sup>)</td><td>23.64</td><td>24.83</td><td></td><td>20.04</td><td>19.14</td><td></td></tr></tbody></table>
neighbouring layers. The shortwave radiation transfer scheme, developed by Pinty et al. (2006), was applied to the albedo calculation. The scheme computes the absorption, transmission, and reflection of incoming radiation by vegetation canopies, which depends on the solar zenith angle, the type of illumination (direct or diffuse), the vegetation type, and the vegetation structure. This scheme considers shortwave radiation both from visible and near-infrared bands and was originally developed for single-layer canopies, but has since been extended for use with layered canopies (McGrath et al., 2016).

b. The wind profile and the vertical eddy diffusivity ($k$; $m^2 s^{-1}$) are calculated using the one-dimensional second-order closure model of Massman and Weil (1999), which makes use of the LAI profile of the stand. It calculates wind profile and vertical eddy diffusivity based on Lagrangian theory.

c. The aerodynamic resistance ($R_b$; $sm^{-1}$) is calculated based upon the leaf boundary-layer resistance, which is estimated according to Baldocchi (1988). The stomatal resistance ($R_s$; $sm^{-1}$) is calculated using a Farquhar–von Caemmerer–Berry-type C3 (Farquhar et al., 1980) and Collatz-type C4 photosynthesis model (Collatz et al., 1992) which simultaneously solves carbon assimilation and stomatal conductance at the leaf level but excludes mesophyll conductance calculation. ORCHIDEE-CAN v1.0 uses an analytical approach as described by Yin and Struik (2009) to calculate layered stomatal resistances which depend on the ambient air temperature, humidity, within-canopy $CO_2$ concentration, vegetation-specific maximum carboxylation rate, and water supply from the roots to the stomata.

Readers are referred to Ryder et al. (2016) for a comprehensive description of the multi-layer energy budget, its assumptions, mathematical details, and a proof of concept. Note that in ORCHIDEE-CAN v1.0 LAI is calculated from a prognostic leaf mass by making use of a vegetation-specific specific leaf area (SLA; $m^2 g^{-1}$). The calculation of the vertical and horizontal distributions of the leaf mass, and thus the vegetation canopy, depends on plant phenology, intrastand competition, forest management, and allometric relationships, and is detailed in Naudts et al. (2015).

## 2.2 Observational data

For this study forest sites were retained if the following data were available: (a) short but intensive campaigns taking flux and profile measurements within and/or below the tree canopy, and (b) multi-year monitoring of top-canopy fluxes. Through numerous regional projects such as CARBOEUROPE, AMERIFLUX, Fluxnet Canada, OZFLUX, ICOS, and NEON, and efforts such as FLUXNET (Baldocchi and Wilson, 2001), multiple year-long time series are

now commonly available, especially for the temperate and boreal zones in Europe, Japan, Australia, and North America. Site selection was thus mostly limited by the availability of within-canopy and below-canopy measurements.

Eight flux observation sites (Table 2) met the aforementioned criteria, and represented various climates from the Mediterranean to the boreal zone and different vegetation types including broad-leaved summer green, broad-leaved evergreen and needle-leaved evergreen. Data were thus missing from needle-leaved summer green vegetation such as larch (*Larix* sp.) and tropical vegetation, so it was not possible to cover all of the forest types that are considered in ORCHIDEE-CAN.

The short intensive campaigns, taking measurements within-canopy and below-canopy, usually extended for periods ranging from several days to a few weeks (Period I; Table 3). During intensive campaigns, vertical profile measurements of wind speed, temperature, and atmospheric humidity were typically conducted. Such measurements were sometimes complemented with profile measurements of sensible and latent heat fluxes, as well as sub-canopy radiation measurements (Period II and III; Table 3). Furthermore, our parametrization and validation set-up required that top-canopy observations had to be available for periods exceeding 1 year (Period IV; Table 3). A typical long-term set-up measured sensible and latent heat fluxes, longwave and shortwave incoming radiation, wind speed, atmospheric temperature, and humidity.

Parametrization and validation utilize the ORCHIDEE-CAN v1.0 model simulations, and so climate forcing data were required to drive the simulations. Site-level weather observation, i.e. shortwave incoming radiation, longwave incoming radiation, two-dimensional wind speed, precipitation, snow, near-surface air pressure, and specific humidity were reformatted and gap-filled using the method proposed by Vuichard and Papale (2015). Weather observations are an integral part of both intensive campaigns and multi-year top-canopy flux monitoring. Hence, within a measurement site, flux, profile, and weather data were usually available at the same temporal resolution and over the same time periods.

Finally, the forcing files were completed with the observed vertical LAI profiles. However, the temporal resolution of LAI was much lower than the resolution of the meteorological variables. When the total LAI was measured at a higher time resolution than its vertical profile, the observed total LAI was vertically distributed according to the observed relative vertical LAI distribution. Model parametrization (Sect. 2.3) and model experiments that aimed at testing the performance of only the multi-layer energy budget (see EXP1 and EXP3 in Sect. 2.5) made use of the observed LAI profiles. For the remaining two model experiments (see EXP2 and EXP4 in Sect. 2.5), ORCHIDEE-CAN v1.0 calculated the vertical LAI profiles following the carbon allocation and carbon turnover schemes, as described in Naudts et al. (2015).
<figure><img src="image_4.png" /></figure>

**Figure S2.** Model simulation and observation of the sensible heat flux profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimized parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.
<table><caption>**Table 2.** Stand structure and data availability of the experimental sites. The maximum observed leaf area (LAI; m² m⁻²) of the over-story and under-story LAI (all-sided) are reported separately. The height of the over-story is expressed in metres. *U* denotes wind speed, *T*<sub>a</sub> denotes atmospheric temperature, and *q*<sub>a</sub> denotes atmospheric humidity. *LE*, *H*, and *R*<sub>n</sub> denote the latent heat flux, the sensible heat flux, and the net radiation, respectively. + indicates that profile measurements were available. — indicates that no profile measurements were available.</caption><thead><tr><th>Site Code</th><th>Fl-Hyy</th><th>FR-LBr</th><th>NL-Loo</th><th>DE-Bay</th><th>CA-Oas</th><th>AU-Tum</th><th>DE-Hai</th><th>BE-Vie</th></tr></thead><tbody><tr><td>Species</td><td><i>Pinus sylvestris</i></td><td><i>Pinus pinaster</i></td><td><i>Pinus sylvestris</i></td><td><i>Picea abies</i></td><td><i>Populus</i> sp.</td><td><i>Eucalyptus</i> sp.</td><td><i>Fagus sylvatica</i></td><td><i>Fagus sylvatica</i><sup>a</sup></td></tr><tr><td>Leaf type</td><td>Needleleaved</td><td>Needleleaved</td><td>Needleleaved</td><td>Needleleaved</td><td>Broadleaved</td><td>Broadleaved</td><td>Broadleaved</td><td>Broadleaved</td></tr><tr><td>Growth form</td><td>Evergreen</td><td>Evergreen</td><td>Evergreen</td><td>Evergreen</td><td>Deciduous</td><td>Evergreen</td><td>Deciduous</td><td>Mixed</td></tr><tr><td>ORCHIDEE PFT</td><td>18</td><td>5</td><td>6</td><td>7</td><td>20</td><td>15</td><td>13</td><td>13</td></tr><tr><td>Over-story LAI</td><td>6.5</td><td>2.0</td><td>1.9</td><td>4.8</td><td>2.9</td><td>2.5</td><td>5.8</td><td>5.1</td></tr><tr><td>Under-story LAI</td><td>0.5</td><td>1.5</td><td>1.5</td><td>0.5</td><td>2.8</td><td>1.0</td><td>0.1</td><td>0.1</td></tr><tr><td>Height</td><td>17.0</td><td>23.0</td><td>15.0</td><td>15.0</td><td>22.0</td><td>50.0</td><td>30.0</td><td>25.0</td></tr><tr><td>U profile</td><td>+</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>T<sub>a</sub> profile</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td></tr><tr><td>q<sub>a</sub> profile</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td></tr><tr><td>LE profile</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>-</td><td>-</td></tr><tr><td>H profile</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>-</td></tr><tr><td>R<sub>n</sub> profile</td><td>-</td><td>+</td><td>+</td><td>+</td><td>+</td><td>-</td><td>+</td><td>-</td></tr><tr><td>Reference</td><td>Launiainen<br>(2007)</td><td>et al.<br>Ogée et al. (2003),<br>Porte et al. (2000)</td><td>Dolman et al. (2002),<br>Moors (2012)</td><td>Foken et al. (2012),<br>Staudt et al. (2011)</td><td>Barr et al. (2004)</td><td>Haverd<br>(2012),<br>Lovell et al. (2012)</td><td>et al.<br>Knohl et al. (2003)</td><td>Aubinet et al. (2001),<br>Laitat et al. (1998)</td></tr></tbody></table>

## 2.3 Model parametrization

At the start of this study the multi-layer energy budget did not yet have a working set of parameters for ORCHIDEE-CAN v1.0. Therefore, we refrained from performing a sensitivity analysis prior to optimizing the model parameters (Kuppel et al., 2014; MacBean et al., 2015), but instead selected three processes described by a total of 10 parameters for optimization. The selected processes were related to the physical processes within the canopy — that is to say, diffusion, advection, and turbulent mixing.

### 2.3.1 Effective drag coefficient $C_{\text{Deff}}$ (unitless)

The canopy structure is a very important characteristic for the land–atmosphere interaction, which can now be simulated by the ORCHIDEE-CAN v1.0 land surface model. We assumed that the drag coefficient is scalar independent and can be parametrized by the canopy structure. The effective drag coefficient is used in the one-dimensional second-order closure wind profile model (Massman and Weil, 1999) that was used to estimate the vertical within-canopy wind profile. In this wind profile model (Massman and Weil, 1999), the drag coefficient is assumed to be a constant throughout the canopy layer, but it can also be treated as a function of the vertical canopy structure.

In this study, we made use of a prototype parametrization approach proposed by Wohlfahrt and Cernusca (2002). Wohlfahrt and Cernusca (2002) provided the basic idea for considering the effective drag coefficient in grasslands that can be varied due to changes in canopy structure, such as bending effects. Thus, we adapted this parametrization to our model; however, we left the first two tuning coefficients ($a_1$ and $a_2$) as constants. This modification allows the effective drag to decrease from a large value to a constant while moving from the top of the canopy to the soil surface layer. Hence, we applied the ideas derived in grassland research to a forest canopy. This approach requires an effective drag coefficient, which relates to the vertically discretized estimate of the canopy drag coefficient ($C_{D,i}$; unitless) and the momentum shielding factor ($P_{m,i}$; unitless) as follows:

$$C_{\text{Deff},i} = C_{D,i} / P_{m,i}. \qquad (1)$$

Both the within-canopy drag and the momentum shielding were parametrized using a function of the cumulative leaf area index (LAI<sub>cum</sub>; m² m⁻²) from the top canopy layer to the bottom layer, which was modified from the original function (Wohlfahrt and Cernusca, 2002) as below:

$$C_{\text{Deff},i} = a_1^{-\text{LAI}_{\text{cum},i}/a_2} + a_3^{-\text{LAI}_{\text{cum},i}/a_4} + a_5, \qquad (2)$$

where the subscript *i* denotes the index of layering from the bottom layer (*i* = 1) to the top-canopy layer (*i* = *n*). $a_1$ to $a_5$ are tuning coefficients (unitless). The default parameter values for $a_1$ to $a_5$ are presented in Table 4.

<footnote><sup>a</sup> This site is partially mixed with *Pseudotsuga menziesii*. <sup>b</sup> The LE profile was available for the 2007 and 2008 period, but not 2011, and the R<sub>n</sub> profile was partly available in 2007.</footnote>
<figure><img src="image_5.png" /></figure>

**Figure S3.** Model simulation and observation of the latent heat flux profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimize parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.
<table><caption>**Table 3.** Observation periods for the different data uses in this study. Date format: dd/mm/yy. The information of the energy closure gap for each site over different selected periods was also calculated based on Chen and Li (2012). EXP1: single-layer scheme with a prescribed LAI profile; EXP2: single-layer scheme with a long-term simulated LAI profile; EXP3: multi-layer scheme with a prescribed LAI profile; EXP4: multi-layer scheme with a simulated LAI profile.</caption><thead><tr><th>Site code</th><th>FI-Hyy</th><th>FR-LBr</th><th>NL-Loo</th><th>DE-Bay</th><th>CA-Oas</th><th>AU-Tum</th><th>DE-Hai</th><th>BE-Vie</th></tr></thead><tbody><tr><td>Period for short-term parameter optimization (Period I)</td><td>01/08/06</td><td>31/07/06</td><td>08/07/97</td><td>04/07/11</td><td>16/08/94</td><td>08/11/06</td><td>10/05/01</td><td>01/08/02</td></tr><tr><td>Closure gap (W m<sup>-2</sup>)</td><td>14/08/06</td><td>05/08/06</td><td>12/07/97</td><td>17/07/11</td><td>22/08/94</td><td>11/11/06</td><td>19/05/01</td><td>07/08/02</td></tr><tr><td>Period for long-term parameter optimization (Period II)</td><td>01/01/02</td><td>01/01/03</td><td>01/01/02</td><td>01/01/97</td><td>01/01/05</td><td>01/06/01</td><td>01/01/05</td><td>01/01/97</td></tr><tr><td>Closure gap (W m<sup>-2</sup>)</td><td>31/12/02</td><td>31/12/03</td><td>31/12/02</td><td>31/12/97</td><td>31/12/05</td><td>31/06/02</td><td>31/12/05</td><td>31/12/97</td></tr><tr><td>Period for single-year EXP1 and EXP3 validation (Period III)</td><td>01/01/05</td><td>01/01/06</td><td>01/01/97</td><td>01/01/99</td><td>01/01/04</td><td>01/06/04</td><td>01/01/01</td><td>01/01/02</td></tr><tr><td>Closure gap (W m<sup>-2</sup>)</td><td>31/12/05</td><td>31/12/06</td><td>31/12/97</td><td>31/12/99</td><td>31/12/04</td><td>31/06/05</td><td>31/12/01</td><td>31/12/02</td></tr><tr><td>Period for multi-year EXP2 and EXP4 validation (Period IV)</td><td>01/01/02</td><td>01/01/03</td><td>01/01/02</td><td>01/01/97</td><td>01/01/04</td><td>01/06/01</td><td>01/01/00</td><td>01/01/97</td></tr><tr><td>Closure gap (W m<sup>-2</sup>)</td><td>31/12/06</td><td>31/12/06</td><td>31/12/06</td><td>31/12/99</td><td>31/12/05</td><td>31/06/05</td><td>31/12/06</td><td>31/12/06</td></tr><tr><td></td><td>10.68</td><td>17.03</td><td>22.65</td><td>48.14*</td><td>3.51</td><td>9.40</td><td>23.69</td><td>33.77</td></tr></tbody></table>

<table><caption>**Table 4.** Description of parameters, code reference, initial values, and tuning ranges used in the multi-layer energy budget model in this work.</caption><thead><tr><th>Parameter name</th><th>Physical parameter</th><th>Empirical representation of</th><th>ORCHIDAS name</th><th>Default value</th><th>Tuning range</th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>Effective surface drag</td><td>Bending of tree branches</td><td>a_1</td><td>6.410</td><td>Use default</td></tr><tr><td>a<sub>2</sub></td><td>Effective surface drag</td><td>Bending of tree branches</td><td>a_2</td><td>0.001</td><td>Use default</td></tr><tr><td>a<sub>3</sub></td><td>Effective surface drag</td><td>Bending of tree branches</td><td>a_3</td><td>0.434</td><td>0.1 to 0.8</td></tr><tr><td>a<sub>4</sub></td><td>Effective surface drag</td><td>Bending of tree branches</td><td>a_4</td><td>-0.751</td><td>-0.9 to -0.1</td></tr><tr><td>a<sub>5</sub></td><td>Effective surface drag</td><td>Bending of tree branches</td><td>a_5</td><td>0.071</td><td>0.05 to 0.1</td></tr><tr><td>a<sub>6</sub></td><td>Eddy diffusivity</td><td>Inner-canopy turbulent mixing</td><td>k_eddy_slope</td><td>5.0</td><td>1.0 to 20.0</td></tr><tr><td>a<sub>7</sub></td><td>Eddy diffusivity</td><td>Inner-canopy turbulent mixing</td><td>k_eddy_ustar</td><td>0.3</td><td>0.0 to 0.6</td></tr><tr><td>a<sub>8</sub></td><td>Surface-atmosphere conductance</td><td>Inner-canopy turbulent mixing</td><td>ks_slope</td><td>5.0</td><td>1.0 to 20.0</td></tr><tr><td>a<sub>9</sub></td><td>Surface-atmosphere conductance</td><td>Under-story phenology</td><td>ks_veget</td><td>0.5</td><td>0.0 to 1.0</td></tr><tr><td>a<sub>10</sub></td><td>Surface-atmosphere conductance</td><td>Under-story phenology</td><td>ks_tune</td><td>1.0</td><td>0.5 to 1.5</td></tr><tr><td>W<sub>br</sub></td><td>Layer boundary resistance</td><td>Upscaling the leaf coupling</td><td>br_fac</td><td>1.0</td><td>0.1 to 10.0</td></tr><tr><td>W<sub>sr</sub></td><td>Layer stomatal resistance</td><td>Upscaling the leaf coupling</td><td>sr_fac</td><td>1.0</td><td>0.1 to 10.0</td></tr></tbody></table>

### 2.3.2 Eddy diffusivity for vertical energy and water transport $k$ (m$^2$s$^{-1}$)

After the vertical wind profile was derived from the one-dimensional second-order closure wind profile model, the friction velocity ($u_*$, m s$^{-1}$), the vertical wind velocity variance ($\sigma_w$; m s$^{-1}$), and the Lagrangian timescale ($T_L$; s) were calculated following the approach by Raupach (1989). In this approach the vertical eddy diffusivity is a function of $\sigma_w$ and $T_L$. Subsequently, the vertical eddy diffusivity down the air

column to the forest floor was calculated as follows:

$$k_i = \sigma_{w,i}^2 T_{L,i} \quad (3)$$

Here we followed the approach proposed by Haverd et al. (2009) for the Lagrangian timescale calculation. The Lagrangian timescale is thus calculated as

$$T_{L,i} = 0.66 \frac{(1 - e^{-4.86(z/h_c)})}{(1 - e^{-4.86})} \frac{h_c}{u_*} \quad (4)$$

A previous effort to validate this model against in situ observations resulted in a bias of the air temperature pro-

<footnote>* The forest in 1997–1999 was strongly affected by forest decline; in 2011 the forest was again in a good state.</footnote>
<figure><img src="image_6.png" /></figure>

**Figure S4.** Model simulation and observation of the net radiation profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimized parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Weighting functions for eddy diffusivity and surface conductance. (a) Weighting function for the eddy diffusivity (k) within the air column (Eq. 3). The weighting is a function of the friction velocity ($u_*$) and was optimized by tuning the parameters $a_6$ and $a_7$. Three different parameter sets show the response of the weighting function to different parameter values. (b) The weighting function for the surface conductance is a function of the vegetation cover and air temperature (Eq. 7). This weighting function was optimized by tuning the parameters $a_8$ to $a_{10}$. Two examples have the following parameter values: $a_8 = 10.0$, $a_9 = 0.5$, $a_{10} = 1.0$, $T_{\text{week}} \ge T_g$, and $T_{\text{week}} = 278.15$. Both of the two cases demonstrate the seasonal cycle of the weighting which will be used to scale the value of $k_{\text{surf}}$. Values to the left of the deflection point show the effect of an increasing/decreasing over-story cover with an increasing/decreasing temperature in spring/autumn. In spring and autumn, under-story growth, and thus its contribution to evapotranspiration, was assumed to be temperature limited. Values right of the deflection point ($a_9 = 0.5$) show the dependency of the evapotranspiration on the soil surface layer on the over-story canopy cover when air temperature is no longer limiting under-story growth.</figcaption></figure>

file within the canopy layer during nighttime (Ryder et al., 2016). These issues have been well documented in the scientific literature (Gao et al., 1989; Dolman and Wallace, 1991; Makar et al., 1999; Wolfe et al., 2011). One possible, although empirical, solution is to apply a different scaling for $k_i$, according to the time of the day. Here we build on a similar approach but, rather than using the time of the day, we used the calculated friction velocity ($u_* = u(h_c) \times (0.32 - 0.264e^{-15.1\zeta(h_c)})$, where $\zeta$ is the cumulative function of $C_{\text{Deff}}$ and $h_c$ is the canopy height) to account for the observed differences in vertical transport within the canopy between daytime and nighttime by applying a weighting factor ($W_{\text{nf}}$; unitless). Therefore the modified diffusivity for level i ($k_i^*$; $m^2 s^{-1}$) was defined as

$$k_i^* = W_{\text{nf}} \sigma_{w,i}^2 T_{\text{L},i}, \quad (5)$$

where $W_{\text{nf}}$ was calculated as

$$W_{\text{nf}} = \frac{1}{1 + e^{(-a_6(u_* - a_7))}}. \quad (6)$$

This function has a sigmoidal shape, where $a_6$ is the ceiling factor of the slope and $a_7$ is the critical friction velocity at the inflection point of the sigmoid function (Fig. 1a). Consequently, atmospheric diffusivity is reduced if $u_*$ is low, which represents stable atmospheric conditions. Under turbulent atmospheric conditions, which are represented by a high $u_*$, $W_{\text{nf}}$ is close to 1 and the simulated diffusivity will

closely follow the relationship proposed by Raupach (1991). The default parameter values for $a_6$ and $a_7$ are presented in Table 4. As an alternative to using $u_*$, it has been proposed to use a mixing length scale to classify flow regimes in order to give a better description of the coupling process below and above the forest canopy (Thomas and Foken, 2007; Staudt et al., 2011; Foken et al., 2012). The numerical scheme of this approach relies on iterations. Since ORCHIDEE-CAN v1.0 is designed to be coupled to regional or global atmospheric models, its numerics has been designed to avoid iterations in order to run efficiently.

Future studies may focus on replacing this empirical solution by a more mechanistic solution. In the context of ORCHIDEE and its coupling to the atmospheric model, this implies that we will have to search for an implicit solution of the near-field far-field theory by Raupach (1989).

### 2.3.3 Conductance for the soil–atmosphere interface $k_{\text{surf}}$ (m s$^{-1}$)

Equation (7) describes the seasonality of the soil–atmosphere interface, which we believe is driven by the under-story and its phenology (Launiainen et al., 2015). Currently, the model does not simulate the production or the phenology of the under-story. As a substitute for this rather complex process, we made use of a weighting coefficient for the conductance of the soil–atmosphere interface ($k_{\text{surf}}$) or, in other words, the calculation of the water vapour exchange between the soil
<figure><img src="image_7.png" /></figure>

**Figure S5.** Model simulation and observation of the air temperature profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimized parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.
layer and the first air column (see the $\phi_{\lambda,E}$ and $K_{surf}$ in Fig. 1 of Ryder et al., 2016, and the formal description of using $K_{surf}$, which is given in the Supplement of Ryder et al., 2016, in Eqs. S4.30 and S4.31).

A relationship between under-story phenology and the conductance for the soil-atmosphere interface has been observed in boreal forest (Launiainen et al., 2015). In winter, when the under-story is senescent, the characteristics in terms of the evapotranspiration at the interface will closely resemble the evapotranspiration of a bare soil. In summer, however, an under-story will be present and its density relates to the gap fraction of the over-story canopy. Hence, the summertime evapotranspiration of the interface will be closer to the evapotranspiration of a vegetation canopy. Therefore, we introduced $\beta_0$ (unitless) as a weighting function ranging from zero to unity, in order to scale the surface conductivity as a function of over-story phenology. Under-story phenology was described as a function of the over-story canopy coverage ($1 - f_{Pgap}$), the mean air temperature during the previous week ($T_{week}$), and a threshold temperature ($T_g$):

$$ \beta_0 = \begin{cases} \frac{a_{10}}{1+e^{(-a_8((1-f_{Pgap})-a_9))}}, & \text{when } G_{\text{veg}} = \text{true,} \\ \frac{a_{10}}{\left(1+e^{(-a_8((1-f_{Pgap})-a_9))}\right) \frac{T_g-T_{\text{week}}}{T_g-273.15}}, & \text{when } G_{\text{veg}} = \text{false,} \end{cases} \quad (7) $$

where $a_8$ is a factor that constrains the slope of the function and $a_9$ is a vegetation cover threshold. $a_{10}$ is a linear weighting factor. $T_g$ is a temperature threshold set to 283.15 K. $G_{veg}$ is a logic variable to indicate the growth status of the vegetation. $G_{veg}$ is an existing variable in ORCHIDEE-CAN v1.0 and depends on a threshold for soil water content and temperature $T_g$. Growth can be expected and therefore $G_{veg}$ is set to true when the weekly averaged soil water content and temperature exceed the thresholds. $f_{Pgap}$ is calculated in ORCHIDEE-CAN v1.0 and describes the over-story gap probability, which is a function of the canopy structure of the vegetation and the solar zenith angle and is calculated in ORCHIDEE-CAN v1.0.

For the lowest layer in the air column, i.e. the layer adjacent to the surface, the surface conductance is then calculated as

$$ k_{\text{surf}} = (W_{\text{sf}}\beta_3 + (1 - W_{\text{sf}})\beta_4) (u_1 C_{\text{Deff},1}), \quad (8) $$

where $\beta_3$ and $\beta_4$ are coefficients describing respectively the fractions of the potential plant transpiration and soil evaporation that are realized. The definition of these coefficients and the numerical approaches are presented in Ryder et al. (2016) and Dufresne and Ghattas (2009). $u_1$ is the wind speed at the lowest canopy layer thus close to the forest floor and is derived from the one-dimensional second-order closure model. $C_{Deff}$ is the effective drag coefficient calculated according to Eq. (2). $W_{sf}$ is the weighting factor for the soil-atmosphere

interface, which is described as the conditional function of over-story canopy cover fraction ($1 - f_{Pgap}$). $W_{sf} = \beta_0$ when ($1 - f_{Pgap}$) > $a_9$, and $W_{sf} = 1 - \beta_0$ when ($1 - f_{Pgap}$) $\le$ $a_9$ (see Fig. 1b). The default parameter values of $a_8$, $a_9$, $a_{10}$, and $W_{sf}$ are presented in Table 4.

### 2.3.4 Boundary-layer resistance of the leaf surface $R_b$ (sm⁻¹)

The boundary-layer resistance of the leaf surface $R_{b,i}$ is described according to the expression from Baldocchi (1988):

$$ R_b = \begin{cases} W_{br} \left( \frac{d_1}{D_{h,air}Nu} \right), & \text{for sensible heat,} \\ W_{br} \left( \frac{d_1}{D_{h,H_20Sh}} \right), & \text{for latent heat,} \end{cases} \quad (9) $$

where $W_{br}$ accounts for the fact that the leaf length of the species under study differs from the characteristic leaf length (unitless), $d_1$ is the characteristic leaf length (0.001 m was used as the default value), $D_{h,air}$ is the heat diffusivity of still air (m² s⁻¹), $D_{h,H_2O}$ is the heat diffusivity of water vapour (m² s⁻¹), $Sh$ is the Sherwood number (unitless), and $Nu$ is the Nusselt number (unitless). The Sherwood number was calculated as $Sh = 0.66 Re^{0.5} Sc^{0.33}$ for laminar flow and $Sh = 0.03 Re^{0.8} Sc^{0.33}$ for turbulent flow, where $Sc$ is the Schmidt number (0.63 for water vapour; unitless). The transition from laminar to turbulent flow takes place in the model when the Reynolds number exceeds a value of 8000. The Nusselt number was calculated as $Nu = 0.66 Re Pr^{0.33}$, where $Pr$ is the Prandtl number (0.7 for air; unitless) (Grace, 1978) and $Re$ is the Reynolds number (unitless) which was calculated as

$$ Re = \frac{d_1 u_i}{\mu}, \quad (10) $$

where $u_i$ is the horizontal velocity at level $i$ (m s⁻¹) and $\mu$ is the kinematic viscosity of air and was set to 0.0015 (m² s⁻¹) (Garratt, 1992). The default parameter value for $W_{br}$ is provided in Table 4.

### 2.3.5 Stomatal resistance $R_s$ (sm⁻¹)

The stomatal resistance of the leaves was calculated for each canopy layer based on the parameters within the layer under consideration. Two stomatal resistances were calculated with the concurrent assimilation rate: (a) the stomatal resistance assuming unlimited soil water availability (the atmospheric demand) and (b) the stomatal resistance that exactly satisfies the amount of water the plant can transport from its roots to its stomata (the plant supply). ORCHIDEE-CAN v1.0 calculates the plant supply of the water available for transpiration as the pressure difference between the soil and the leaves divided by the sum of hydraulic resistances of fine roots, sapwood, and leaves (see Eq. 20 in Naudts et al., 2015). The atmospheric demand for water for transpiration is calculated as the vapour pressure difference between the leaves and atmosphere divided by the sum of boundary-layer resistance
<figure><img src="image_8.png" /></figure>

**Figure S6.** Model simulation and observation of the leaf temperature profile at eight forest sites during the short-term campaign (Period I). All the dashed lines indicate the prior simulation with default parameter values and the solid lines present the optimized simulation with optimized parameter values. The filled circles are the observation means and the bars are stand deviations over the simulation period at 13:00. The gray bars in the background indicate the measured maximum LAI at each level in the reference year.
(Rb) and stomatal resistance (Rs) (see Eqs. 9 and 13 in Ryder et al., 2016). When the supply can satisfy the demand, there is no water stress and photosynthesis (A) is calculated. When the demand is limited by the supply term, A and Rs are recalculated such that they satisfy the supply. Water stress thus enters Eq. (11) in the value of A. ORCHIDEE-CAN v1.0 scales stomatal resistance to account for the part of the canopy that is coupled to the atmosphere and thus contributes to the latent heat flux. In this study, this weighting was formalized through a linear parameter Wsr:

$$R_{s,i} = W_{sr} \left( \frac{1}{\left( g_0 + \left( \frac{A_i h_s}{C_s} \right) \right) \text{LAI}_i} \right), \quad (11)$$

where g0 is the residual stomatal conductance if the solar irradiance approaches zero, C_s is the concentration of CO2 at the leaf surface, and hs is the relative humidity at the leaf surface. A is the CO2 assimilation rate which is solved analytically following Yin and Struik (2009). In Eq. (11) the relative humidity used is the top-canopy forcing instead of a layered relative humidity in order to avoid an iterative process. The default parameter value for Wsr is presented in Table 4.

## 2.4 Model optimization

### Optimization procedure

Parametrizing the scaling coefficients and weighting factors enabled us to simultaneously improve the match between the simulated and observed sub-canopy micrometeorology, including temperature and specific humidity when available, and between the simulated and observed top-canopy heat fluxes (LE and H). Within-canopy fluxes were also simulated but are not usually measured. The parametrization made use of an in-house optimization package called ORCHIDAS (ORCHIDEE Data Assimilation Systems; http://orchidas.lsce.ipsl.fr/). ORCHIDAS provides a range of numerical approaches for assimilating multiple data streams in ORCHIDEE.

We used the maximum gradient approach to tune the parameters a3 to a10, Wbr, and Wsr for each study site independently. Over the course of several iterations, the optimization approach minimized the mismatch between the model output and the observations, using a gradient-based algorithm called L-BFGS-B (Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with Bound constraints), which provides the possibility to prescribe boundaries for each parameter (Byrd et al., 1995). The range assigned to each parameter is reported in Table 4. Furthermore, this approach allowed for measurement uncertainties in the eddy covariance LE measurement by reducing its weight in the cost function from 1.0 to 0.66. This value of 0.66 was set based on the outcome of a paired tower experiment to estimate the random errors of the eddy covariance measurements (Richardson et al., 2006).

For the optimization, the LAI in ORCHIDEE-CAN v1.0 was set to match the observed vertical LAI profile.

A three-step optimization procedure was carried out in this study. Firstly, the within-canopy and below-canopy observations from the short-term intensive measurement campaigns (Period I in Table 3) were used to optimize a3 to a7, Wbr, and Wsr. During this step, the parameters for the soil-atmosphere interface (ksurf, i.e. a8 to a10 and Wsf) were set to their default values. Since these campaigns took place during summer, parameters related to the within-canopy effective drag profiles, eddy diffusivity, boundary layer resistance, and stomatal resistance (CDeff; k; Rb; Rs) were biased towards the summer. Secondly, the seasonal dynamics of ksurf was parametrized by trying to improve the correspondence between the simulated and observed top-canopy fluxes over 1 year (Period II in Table 4). In this step, a3 to a7, Wbr, and Wsr were set to the values obtained from the first step of the optimization and a8 to a10 and Wsf were tuned. Finally, the performance of the calibrated model was evaluated based on a second single year of top-canopy observations (Period III in Table 3).

Although the spin-up was stopped on 30 June (Table S1 in the Supplement) and all simulations thus used the 30 June soil water content as their initial condition, this approach does not guarantee that this typical summer soil water content matches the soil water content in the year of the intensive measurement campaign. The effect of this possible mismatch was quantified by running a sensitivity analysis in which the whole parametrization approach, which was repeated for seven different initial soil water contents, varied from -30 to 30% in increments of 10% of the 30 June value.

## 2.5 Attribution of changes in model performance

The multi-layer energy budget scheme (Ryder et al., 2016) that was parametrized and tested in this study required realistic spatial and temporal soil water content and a value for the ground heat flux from surface level as initial conditions. This need was satisfied by implementing this scheme within the newly enhanced ORCHIDEE-CAN v1.0 land surface model (Naudts et al., 2015). Integration of the multi-layer energy budget into ORCHIDEE-CAN v1.0, however, complicated the design of the validation study, as it was now necessary to separate, as much as possible, the performance of the multi-layer energy budget scheme from the performance of the rest of the model. To this aim, four experiments were designed in order to better understand the performance of the new scheme (Table S1 in the Supplement).

### - Experiment 1 (EXP1): single-layer scheme with a prescribed canopy

The first experiment was run at the site level and made use of the default single-layer energy budget scheme. The energy budget scheme was driven by the observed climate forcing and the observed total LAI (Table 2). In this experiment, the vertical LAI profile was only
<figure><img src="image_9.png" /><figcaption>**Figure S7.** Sensitivity test of using default *k*<sub>surf</sub> value with different initial soil moisture conditions to determine optimized parameter values for short term period at FR-LBr site. (A) parameters from *a*<sub>3</sub> to *a*<sub>5</sub> to determine the effective surface drag coefficient, *C*<sub>Deff</sub> (B) parameters *a*<sub>6</sub> and *a*<sub>7</sub> to determine the weighting factor for eddy diffusivity, *W*<sub>nf</sub> (C) parameter from *a*<sub>8</sub> to *a*<sub>10</sub> to determine the weighting factor for surface-air interface conductance, *W*<sub>sf</sub> (D) weighting factor for stomatal resistance *W*<sub>sr</sub> and boundary layer resistance *W*<sub>br</sub>, respectively.</figcaption></figure>
used for the photosynthesis module in ORCHIDEE-CAN v1.0. Note that vertical LAI profiles cannot be used by the single-layer scheme and the results are therefore limited to the top-canopy fluxes. This experiment was used as the reference simulation to document the performance of the single-layer approach.

- **Experiment 2 (EXP2): single-layer scheme with a simulated canopy**

The second experiment was identical to the first experiment, except that the LAI was now simulated by ORCHIDEE-CAN v1.0, rather than using the observed LAI. Given that these experiments make use of observed climate drivers and LAI, changes in model performance between experiments 1 and 2 are derived by the introduction of a dynamic and prognostic vertical LAI profile. A large decrease in performance between experiments 1 and 2 would suggest that ORCHIDEE-CAN v1.0 does a poor job in simulating the vertical LAI profile.

- **Experiment 3 (EXP3): multi-layer scheme with a prescribed canopy**

Experiment 3 differs from EXP1 through the use of the multi-layer energy budget scheme, rather than the single-layer scheme. As a consequence, the observed vertical LAI profiles rather than the observed total LAI are now applied to drive the simulations with a multi-layer energy budget. This experiment was used for quantifying the change in performance when switching from the single-layer approach to the multi-layer approach. Although these simulations calculate the turbulent fluxes for each canopy level, the change in performance was based on a comparison of experiments 1 and 3, and as such the analysis had to be limited to the top-canopy fluxes, as within-canopy fluxes cannot be calculated by the single-layer approach used in the first experiment. A large decrease in performance between experiments 1 and 3 would suggest that the multi-layer energy budget in ORCHIDEE-CAN v1.0 does not help to better simulate the top-canopy fluxes.

- **Experiment 4 (EXP4): multi-layer scheme with a simulated canopy**

In Experiment 4 the vertical LAI profile was calculated by ORCHIDEE-CAN v1.0. Thus, this experiment made use of the full functionality of ORCHIDEE-CAN v1.0 and the multi-layer energy budget. As such, albedo, photosynthesis, and the energy budget calculations were fully consistent. Comparing the performance of experiments 2 and 4 quantifies the actual change in performance for a prognostic LAI profile and its interactions in ORCHIDEE-CAN v1.0. A large decrease in performance between experiments 2 and 4 would therefore suggest that the multi-layer energy budget in ORCHIDEE-CAN v1.0 does not help to better simulate

the top-canopy fluxes. Furthermore, a large decrease in performance between experiments 3 and 4 would indicate that ORCHIDEE-CAN v1.0 does a poor job in simulating the vertical LAI profile.

All four experiments were started from 20 years of spin-up simulations, which were driven by CRU-NCEP climate re-analysis from 1991 to 2010 with a spatial resolution of 0.5° × 0.5° (Maignan et al., 2011) at selected study sites. These spin-up simulations allow the model to build up a realistic soil water pool at the start of each simulation. The climate forcing to spin up the model can be obtained from local high-resolution climate observations for a usually very limited time period or low-resolution regional re-analysis for a much longer time period. Using the local high-resolution data would have the advantage that local information is used, but due to the fact that some time series are only 2 to 4 years long (Table 3, Period IV), the spin-up would have to cycle 5 to 10 times over the same data. Although local data could then still have been used, cycling gives a lot of weight to the climatic events in the time series and may as such result in a biased spin-up. The alternative is to use 20 years of a climate re-analysis; these data represent the inter-annual variability better than cycling over the same 2 or 4 years of data, but this has the disadvantage that the data are less likely to represent the local conditions (especially in mountainous regions). Given the fact that we did not have access to soil water content data, we could not evaluate which method is better to spin up the soil water content in the model. For this reason, we performed a sensitivity analysis of the parametrization of the initial soil water content at one of the driest sites used in this study (see Sect. 3.1).

A 10-layer LAI profile was applied for each site – the number of layers chosen follows the approach from a previous study (Ryder et al., 2016). If the vertical LAI profile was prescribed, the total LAI was re-scaled within these 10 layers to follow the observed vertical LAI profile at each site (Fig. 2). If the vertical LAI profile was not imposed, the LAI generated for the albedo calculation (McGrath et al., 2016) was used instead. Note that, contrary to previous versions of ORCHIDEE, ORCHIDEE-CAN v1.0 no longer places a constraint on the maximum LAI. In ORCHIDEE-CAN v1.0, the total LAI is the outcome of carbon allocation to the canopy through a pipe model and carbon removal from the canopy through leaf turnover (Naudts et al., 2015).

## 2.6 Model performance

The change in model performance due to the use of the multi-layer rather than single-layer scheme for a prescribed LAI profile (EXP1 vs. EXP3), and a simulated LAI profile (EXP2 vs. EXP4), were quantified by comparing the Taylor skill score ($S_T$) (Taylor, 2001).

$S_T$ was calculated for the eight observational sites for the top-canopy fluxes of all four experiments making use of the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Vertical LAI profile for maximal total LAI. The LAI was discretized in 10 evenly spaced layers and the canopy height was normalized. The canopies of FI-Hyy, DE-Bay, DE-Hai, and BE-Vie were considered dense (over-story LAI > 3.0), whereas the canopies of FR-LBr, NL-Loo, CA-Oas, and AU-Tum were considered sparse (over-story LAI ≤ 3.0).</figcaption></figure>

simulated and observed half-hourly fluxes. The Taylor skill score was calculated as follows:

$$S_T = \frac{4(1 + R)}{(\hat{\sigma}_f + 1/\hat{\sigma}_f)^2 (1 + R_0)}, \quad (12)$$

where R is the correlation coefficient between the simulation and the observation, R₀ is the maximum correlation coefficient and $\hat{\sigma}_f$ is the ratio of the variance of the simulations to the variance of observations ($\hat{\sigma}_f = \sigma/\sigma_r$). Here, we set R₀ to 1.0 for the maximum correlation between observation and model simulation. A value of 1.0 of S_T indicates that model simulations perfectly match the observations; values lower than 0.5 imply that the model has poor predictive ability.

# 3 Results

## 3.1 Model parametrization

Using the default parameter set (i.e. a₁ to a₅) resulted in an underestimation of the wind speed in the lower canopy level at all study sites. Optimized parameters could be roughly grouped according to canopy structure (see Table S1 in the Supplement). For forest sites with a dense canopy (see the second row of Fig. S1 in the Supplement), the parameters had to be adjusted to simulate a low wind speed in the lower

canopy. For forest sites with a sparse canopy, the parameters had to be adjusted to simulate relatively high wind speeds at the bottom of the canopy. At these sites, flux observations showed a substantial contribution from the forest floor to the sensible and latent heat fluxes at the top of the canopy. The average model error of wind profile estimation, in terms of root mean square error (RMSE), was reduced from 0.62 to 0.42 m s⁻¹ after adjusting the parameters (see Table S3 in the Supplement). Tuning the conductance of the soil-atmosphere interface (i.e. a₈ to a₁₀), rather than tuning the stomatal conductance and leaf boundary-layer resistances, enabled a closer match between the simulations and observations (Figs. S2 and S3 in the Supplement).

At sites with dense canopies, however, tuning the weightings of stomatal resistance and weighting the boundary layer resistance improved the match between the simulated and observed inner-canopy and top-canopy fluxes of sensible and latent heat (Figs. S2 and S3 in the Supplement). The model errors of heat and water flux estimations were reduced substantially from 91.2 to 46.1 W m⁻² for LE and 123.2 to 50.3 W m⁻² for H, respectively (see also Table S3 in the Supplement).

At sites with sparse canopies, the net radiation at the forest floor was substantial, i.e. ranging nearly from 200 to 450 W m⁻² (Fig. S4 in the Supplement). Correctly simulating
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Inter-annual and diurnal performance for both dense and sparse forest types, expressed as Taylor skill score ($S_T$), of the single-layer energy budget scheme. Taylor skill score was calculated for each component in the energy budget. Simulations made use of the single-layer energy budget scheme in ORCHIDEE-CAN v1.0 according to the settings described for Experiment 1 (EXP1). Taylor skill scores were aggregated according to canopy density (dense vs. sparse). A value of 1.0 of $S_T$ indicates that model simulations perfectly match the observations; values lower than 0.5 imply that the model has poor predictive ability. FI-Hyy, DE-Bay, DE-Hai, and BE-Vie are dense forest sites, and FR-LBr, NL-Loo, CA-Oas, and AU-Tum are sparse forest sites.</figcaption></figure>

ing radiation transfer strongly contributed to correctly simu-
lating the within-canopy flux profiles and top-canopy latent
and sensible heat fluxes. Nevertheless, radiation transfer was
not re-parametrized in this study and, hence, the model er-
rors of net radiation estimation depended solely on the tree
species. In sparse canopies, a positive air temperature gra-
dient with higher temperatures at the forest floor compared
to the top-canopy was also presented (Fig. S5 in the Supple-
ment). Using default parameter values for all factors resulted
in a good simulation of the air temperature gradient for all
eight sites. However, optimizing the parameters (i.e. $a_3$ to
$a_{10}$, $W_{br}$, and $W_{sr}$) had a large impact on the absolute val-
ues of the vertical profile in leaf temperature (Fig. S6 in the
Supplement). Unfortunately, leaf temperature was not mea-
sured at any of the sites. Therefore, it remains to be assessed
whether the model can concurrently reproduce observed en-
ergy fluxes and soil water content.

At one site with an open canopy (FR-LBr) the effect of
the initial soil water content on the optimized parameter esti-
mates was tested. Both the stomatal resistance and the bound-
ary resistance weighting factors ($W_{sr}$ and $W_{br}$) were found to
be very sensitive to the optimization procedure, with changes
in their values exceeding 5% (Fig. S7 in the Supplement).

After parameter adjustment the sensitivity to the initial soil
water content was 5% less than that using the originally
optimized values. Changes in parameters $a_6$ and $a_7$, which
tuned the eddy diffusivity, were largely unaffected by the ini-
tial conditions. Soil water content measurements would thus
have helped to improve the parametrization, especially for
the stomatal and leaf boundary-layer resistances.

**3.2 Performance of the single-layer scheme**

Model performance of the single-layer model was evaluated
by making use of EXP1. Overall model performance for
sparse canopies (Fig. 3a) was slightly higher and thus better
than model performance at the dense forest sites (Fig. 3b).
Moreover, model performance at the forests with sparse
canopies showed less variability within a year than model
performance at sites with a dense canopy.

At the sparse canopy sites, both the intra-annual and diurnal variations in net radiation $R_n$ were well simulated, displaying $S_T$ scores continuously over 0.9 (Fig. 3b and d). For dense canopies, the $S_T$ score of $R_n$ dropped to 0.9 in winter, which might be attributed to an incorrect estimation of $R_n$ during nighttime (Fig. 3c).
<figure><img src="image_6.png" /></figure>

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

**Figure 4.** Inter-annual and diurnal performance for both dense and sparse forest types, expressed as Taylor skill score ($S_T$), of the multi-layer energy budget scheme. Taylor skill score was calculated for each component in the energy budget. Simulations made use of the multi-layer energy budget scheme in ORCHIDEE-CAN v1.0 according to the settings described for Experiment 3 (EXP3). Taylor skill scores were aggregated according to canopy density (dense vs. sparse). A value of 1.0 of $S_T$ indicates that model simulations perfectly match the observations; values lower than 0.5 imply that the model has poor predictive ability. FI-Hyy, DE-Bay, DE-Hai, and BE-Vie are dense forest sites, and FR-LBr, NL-Loo, CA-Oas, and AU-Tum are sparse forest sites.

In general, the $S_T$ for the single-layer or big-leaf model for the sensible heat flux was higher than for the latent heat flux, both at the annual and daily resolution. The $S_T$ dropped below 0.5 for latent heat flux and 0.8 for sensible heat flux (Fig. 3a) from December to February (or June to August at Au-Tum), indicating that the single-layer model incorrectly partitioned energy during the cold season (Fig. 5c and e). During these months nights are long and the inability of the model to simulate nighttime fluxes (Fig. 3c) may well be the cause of the observed model deficiencies during the winter months. The low model performance on latent heat flux estimation was due to the model overestimation during these months (see Fig. 5e).

### 3.3 Performance of the multi-layer scheme

Model performance of the multi-layer model was evaluated by making use of EXP3. By introducing the multi-layer energy budget scheme, model performance for sparse and dense canopies became more comparable (Figs. 4a, b and 5e, f) due to small improvements in the $S_T$ for simulation of dense canopies and small losses in the skill to simulate the energy budget of sparse canopies. Improved simulations of nighttime fluxes under dense canopies (Figs. 4c, 6c and e) were

reflected in the improved partitioning of energy fluxes during wintertime (compare Figs. 3a and 4a). The multi-layer energy budget model gains some skills compared to the single-layer model in the simulation of the latent heat flux from sparse canopies between December and April (see Fig. 5f).

Overall, the introduction of the multi-layer energy budget and its integration into ORCHIDEE-CAN v1.0 resulted in a small decrease in model skill (Fig. 7; Table S4 in the Supplement). When moving from the single-layer scheme with a prescribed LAI (EXP1) to the multi-layer scheme with a prescribed LAI profile (EXP3), the model skill decreased for $R_n$, $H$, and LE but increased for $G$ (see Fig. 5g and h, and 7). Note that $G$ is an essential aspect in simulating the snow phenology (Wang et al., 2015). Therefore, improved simulations of the soil heat fluxes could have important indirect effects on climate simulations of regions with a pronounced snow season.

Despite this improvement, the overall model performance in the ground heat flux estimation at all eight forest sites was still very low < 0.5 (Fig. 4b and c; Table S4 in the Supplement). The low performance may be due to either deficiencies in the model or the inability of point measurements to represent the large variation in ground heat fluxes under-
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Inter-annual variation of measured and simulated energy fluxes. The lines indicate mean values of selected sites (dense or sparse forests). The observed mean is shown as a solid line, and the simulations of the single-layer energy budget scheme (EXP1) and the multi-layer energy budget schemes (EXP3) are shown as a dashed and dotted line, respectively. The symbols represent the monthly averaged values of energy fluxes at one site (in the first column sub-plot **a**, **c**, **e**, and **g** are dense forest sites; FI-Hyy: □, DE-Bay: ○, DE-Hai: △, and BE-Vie: ◊; in the second column sub-plot **b**, **d**, **f**, and **h** are sparse forest sites; FR-LBr: □, NL-Loo: ○, CA-Oas: △, and AU-Tum: ◊). Simulation results are symbols filled with colours: EXP1 with grey colour and EXP3 with black colour.</figcaption></figure>

neath a canopy or the errors made in estimating the rate of
heat storage change in the layer of soil between the soil heat
flux plates and the soil surface (Mayocchi and Bristow, 1995;
Kustas et al., 2000). However, the small loss (all fluxes ex-
cept G) or gain (only for G) in model skill from introducing

the multi-layer scheme can be strengthened (i.e. LE) or com-
pensated for (R<sub>n</sub>, H, and G) by the small gain in model skill
from the introduction of a prognostic vertical LAI profile.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Mean diurnal cycle of measured and simulated energy fluxes. The lines indicate mean values of selected sites (dense or sparse forests). The observed mean is shown as a solid line, and the simulations of the single-layer energy budget scheme (EXP1) and the multi-layer energy budget scheme (EXP3) are shown as dashed line and dotted line, respectively. The symbols represent the annual hourly averaged values of energy fluxes at one site (in the first column sub-plot **a**, **c**, **e**, and **g** are dense forest sites; FI-Hyy: □, DE-Bay: ○, DE-Hai: △, and BE-Vie: ◊; in the second column sub-plot **b**, **d**, **f**, and **h** are sparse forest sites; FR-LBr: □, NL-Loo: ○, CA-Oas: △, and AU-Tum: ◊). Simulation results are symbols filled with colours: EXP1 with grey colour and EXP3 with dark-grey colour.</figcaption></figure>

# 4 Discussion

## 4.1 Single-layer vs. multi-layer energy budget

Three major deficiencies of the single-layer energy budget scheme have been identified: (1) poor model performance in the net radiation estimation during nighttime in dense canopy

forests; (2) incorrect energy partitioning during winter seasons at dense forest sites; and (3) incorrect simulation of soil heat flux for all forest sites. These site-level findings are consistent with previous large-scale validation work (Pitman et al., 2009; Jiménez et al., 2011; de Noblet-Ducoudré et al., 2012), which applied the single-layer energy budget to sim-
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Change in model performance, expressed as Taylor skill score, with increasing experimental complexity for both the single-layer and multi-layer energy budget schemes for all eight study sites. EXP1: single-layer scheme with a prescribed LAI profile; EXP2: single-layer scheme with a simulated LAI profile; EXP3: multi-layer scheme with a prescribed LAI profile; EXP4: multi-layer scheme with a simulated LAI profile.</figcaption></figure>

ulate land surface fluxes dynamically and demonstrated that this approach has difficulties in the reproduction of surface energy fluxes.

In this study, we tried to overcome these difficulties by implementing a multi-layer energy budget scheme. The multi-layer energy and water calculations make use of a vertically resolved radiation transfer scheme for shortwave and longwave radiation (replacing prescribed shortwave reflection values), a within-canopy wind velocity profile (replacing empirical formulations for roughness length), a vertical prognostic LAI profile (replacing a prescribed LAI value), within-canopy leaf boundary-layer resistance profiles for energy and water transport, a within-canopy stomatal resistance profile, a vertical discrete eddy diffusivity profile, and a soil-atmosphere layer conductivity.

This approach resulted in small improvements in simulating energy partitioning during nighttime for dense canopies, small losses in model performance in terms of energy partitioning for sparse canopies, and year-round gains in model performance for simulation of the ground heat flux. As such, the multi-layer energy and water vapour flux scheme did not solve the long-standing issues related to simulating nighttime energy partitioning (Jordan and Smith, 1994; Prihodko et al., 2008; Wild, 2009; He et al., 2011), but it succeeded in obtaining a similar model performance, while much of the empiricism of the big-leaf approach was replaced by a more realistic process description. A more realistic model description opens new avenues of research (see Sect. 4.3).

## 4.2 Parametrization approach

Despite the direction of the land surface model community towards the development of more mechanistic models, all large-scale land surface models contain an important level of empiricism. When the model is carefully developed and validated, the empirical parameters mimic an overly complex (for the purpose of the model) or incompletely understood process. As we tried to follow this philosophy, we believe that our parameters have a plausible natural background (Table 4), but this does not overcome the issue of equifinality of the model. Ideally, future developments should aim at replacing such parameters by a more mechanistic approach if the empirical module represents a process that is at the core of the objectives of the model. In this study, the parametrization of the new scheme and its underlying processes revealed strengths and weaknesses of the model as well as avenues for future experimental work.

### 1. Within-canopy drag

For the inner-canopy drag parametrization, we modified an approach (Eq. 2) that has previously only been tested and validated at grassland sites (Wohlfahrt and Cernusca, 2002). In that study, LAI was treated as equal to the plant area index (PAI), which is a separate measure that accounts not only for leaves, but also for other vegetation material such as stems and seed heads. In forests, however, the difference between LAI and PAI is made up of the branches and trunks and becomes especially important in winter in deciduous stands as canopy drag still exists. As a first parametrization this simplification allowed a better comparison with the observations and with the single-layer model. We applied a formulation that makes use of LAI and, by doing so, some model errors might have been introduced, especially for the deciduous forest sites. ORCHIDEE-CAN v1.0 now simulates both LAI and PAI, and so this enhanced approach could be adopted. Results confirmed that substituting PAI by LAI is acceptable during the leaf-on seasons.

Alternative approaches have been proposed by Cescatti and Marcolla (2004). For example, the inner-canopy drag could also be modelled as the function of the percentage of horizontal gaps in the forest canopy – a canopy characteristic that is presently simulated in ORCHIDEE-CAN v1.0. Measurement sites such as DE-Bay or AU-Tum have detailed wind and vertical LAI profile observations and could thus be used in a pilot study for developing a suitable parametrization approach linking inner-canopy drag and shielding to the canopy gaps. Such a development would also meet the requirements for calculating drag and shielding following small-scale mortality from forest management, fires, wind damages, and pests.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Effect of under-story phenology on the vertical profile of the latent and sensible heat fluxes at the FR-LBr site. (a) Simulated latent heat flux assuming that the interface between the soil and the lowest atmospheric layer behaves as a bare soil (dotted line), a fully vegetated surface (dashed line), or a partly vegetated, partly bare surface where the ratio between bare soil and vegetated soil depends on the under-story phenology (full line). The observed profile is shown as black dots where the error bars denote the 5-day temporal variance. (b) Simulated sensible heat flux assuming that the interface between the soil and the lowest atmospheric layer behaves as a bare soil (dotted line) or a fully vegetated surface (dashed line), or depends on the under-story phenology (full line). The observed profile is shown as black dots where the error bars denote the 5-day temporal variance.</figcaption></figure>

## 2. Within-canopy transport

In this study, within-canopy transport was parametrized by K-theory. A one-dimensional second-order closure model was applied to derive the within-canopy turbulence statistics, based both on the LAI profile and the canopy height. This approach has been reported to produce a reasonable approximation of above-canopy fluxes estimation, even if the within-canopy temperature and humidity gradients are not always well captured (Raupach, 1989). As previous studies have demonstrated, incorrect estimation on gradients may be accommodated to some extent by introducing a scaling factor (Eq. 6) to constrain the within-canopy transport (Makar et al., 1999; Wolfe et al., 2011; Ryder et al., 2016). Alternatively, such a scaling factor might vary in terms of the form of the canopy structure or openness, though the determination of the factor has yet to be adequately described due to a restricted range of measurements (McNaughton and Van Den Hurk, 1995; Stroud et al., 2005).

At sparse forest sites, the temperature measurements showed a general positive gradient during the daytime (Fig. S5 in the Supplement) and a negative gradient during the nighttime (not shown). For the sparse forests, the temperature gradient is even more complex, having

a negative or reversed gradient throughout the vertical profiles. By using the current parametrization approach, most of the sparse forest sites required a higher shear stress (a stronger threshold friction velocity $a_7$) for the within-canopy mixing, compared to dense forest sites (Table S2 in the Supplement), in order to replicate the measurement results. This observation relates to a general difficulty in being able to simulate canopy transport based on limited general measurements (Stroud et al., 2005).

## 3. Sub-canopy and surface-atmosphere conditions

In this study, we treated the under-story and over-story as the same species to construct the vertical LAI profile based on the observed LAI profile. This treatment only allowed the under-story growth to follow over-story canopy phenology. In fact, the forest floor is often occupied by plants with very different traits, of which one of the most obvious is the difference in leaf onset and/or leaf fall (Barr et al., 2004). Given the aforementioned model formulation, simulation of the under-story phenology and traits could be further improved in the future. For example, over-story and under-story vegetation could be simulated as different plant functional types or plant species within the same energy budget
column. Also, the microclimate created by the over-
story could be used as an input to simulate the envi-
ronmental conditions in the under-story.

Starting from the point of view of the interaction be-
tween ecosystems and the climate, we introduced a
weighting factor (W<sub>sf</sub>) as a function of a long-term av-
erage temperature and light conditions (gap fraction), a
transpiration fractions described as β<sub>3</sub> in the model code
and a soil evaporation fraction (β<sub>4</sub>) as environmental
factors to parametrize surface conductance (Fig. 8 and
consequently control the surface latent heat flux. This
approach demonstrated the model's capability to simu-
late the flux profile in agreement with observations. It
may, however, not be valid for the savanna ecosystem
because the under-story phenology of this ecosystem
relies on water availability in the top soil layer (Bal-
docchi and Wilson, 2001; Hutley et al., 2000), which is
an environmental condition not accounted for in our ap-
proach. Furthermore, accounting for ecosystem-specific
differences in root density profiles and aerial cover of
the under-story might also help in the simulation of wa-
ter and energy fluxes (El Masri et al., 2015; Launiainen
et al., 2015). From this perspective, detailed soil moist-
ture profile observations would be very useful in de-
veloping a more advanced surface-atmosphere interface
parametrization.

4. Mismatch between low-resolution driver data and vertically resolved vegetation layers

In this study an apparent mismatch was present between
the low resolution of the driver data that contain infor-
mation derived from several different land cover types
and the highly resolved vertical layering of the canopy.
When low-resolution driver data are used, the bene-
fit from replacing the big-leaf approach in favour of a
multi-layer approach becomes questionable.

In this study the spin-up of the soil water content made use of low-resolution driver data, but the simulations themselves were driven by spatially and temporally high-resolution site observations. Nevertheless, the apparent mismatch touches upon an interesting issue: how does one account for the average surface fluxes from the contribution of different subgrid-scale land cover types? The present ORCHIDEE single-layer model calculates a weighted average of different PFTs across a grid square to calculate a total representative flux. An alternative approach, and one that we are investigating using this multi-layer model in ORCHIDEE-CAN v1.0, is to calculate the heat fluxes of each vegetation type separately (sub-grid-scale modelling) so that the mixing occurs above the canopy.

5. The proposed parametrization approach and future work
In general, we provide a simple but useful parametriza-
tion approach for the multi-layer energy budget scheme

in the ORCHIDEE-CAN v1.0 global land surface
model. Comparing with other studies (Ogée et al., 2003;
Staudt et al., 2011; Launiainen et al., 2015), our ap-
proach directly determines the energy and water fluxes
and successfully avoids the iterative processes to meet
the numerical requirement. In total, a set of 12 param-
eters need to be prescribed and calibrated regarding the
empirical representation of surface drag, turbulent mix-
ing, sub-canopy phenology, and leaf-atmosphere cou-
pling processes. Our approach presents a good perfor-
mance at all study sites, though we may have some
deficits in wind speed estimation.

In this study the model had been tested for several en-
vironmental conditions and demonstrated that the nu-
 merics can deal with the variation that can be found in
global ecosystems. A separate parameter set for each
site has been provided. Next, we will have to derive a
single parameter set for each PFT and test how well the
model reproduces global patterns in, for example, evap-
otranspiration. Only then will we be able to learn about
the transferability of the parameters from the site level
to the PFT level.

4.3 Increased model capacity

The innovation of the multi-layer energy and water scheme is the capacity to simulate the behaviour of fluxes within the canopy and the separation of the soil-level temperature from the temperature of the vegetation levels. The multi-layer scheme helps to address how forest management such as thinning or shelterwood cutting may alter the forest-atmosphere coupling and resulting fluxes. It also paves the way for the consideration of mixed forests where different plant species or functional types can be in a different microclimatic environment to that of the high canopy. This capacity is essential for the following types of future potential applications.

1. The simulation of emission of biogenic volatile organic compounds (BVOCs) from plants, linking climate change, atmospheric chemistry, and the terrestrial biosphere. The implemented multi-layer energy and water budget calculates the leaf temperature and within-canopy radiation, and therefore would potentially allow us to simulate the emission of BVOCs, such as isoprene or monoterpene from plants (Guenther et al., 1995, 2006).

2. Natural disturbances, such as fires, pests, and windfall, can result in increases in leaf fall, individual tree mortality, or complete stand destruction (Lugo, 2008; Seidl et al., 2011; Yue et al., 2014), which in turn determine the vertical LAI profile. The implemented multi-layer energy and water budget scheme calculates the vertical eddy diffusivity and effective drag coefficient as a function of the vertical LAI profile; hence, the new scheme
allows the study of effects of changes in disturbance intensity on the energy budget and thus the climate system.

3. Forest canopy structure plays an important role in regulating the provision of forest ecosystem services such as maintaining biodiversity (Scheffers et al., 2013; Defraeye et al., 2014) or regulating streamflow (Jackson, 2005). Therefore, structural changes to the forest canopy, through, for example, forest thinning or species changes, will reduce the buffering effect of the canopy. It is only with models including a multi-layer energy budget that an informed prediction of the long-term consequences of land-management policies can be made.

4. This work takes the first step in exploring the use of vertical canopy profiles in coupled vegetation/atmospheric models, particularly in relation to the calculation of GPP, which is sensitive to the vertical profiles of light, water and nitrogen (Bonan et al., 2012, 2014). To run at a regional or global scale, it is essential to first parametrize the model at the site level.

## 5 Conclusion

Although the first parametrization of a multi-layer energy and water budget scheme did not greatly improve the model performance over the use of the so-called big-leaf approach for energy and water calculations, it provides a more detailed description of the within-canopy micrometeorology of various forest types. A more detailed process description is essential when linking climate change to studies addressing, for example, species vulnerability to climate change, the climate feedbacks from different disturbance intensities, changes in under-story habitat following management changes, and BVOCs as a result of climate change.

In this study, multiple-site calibration and optimization were performed in order to better understand the functionality of the newly implemented multi-layer energy budget in ORCHIDEE-CAN v1.0. Developing the multi-layer energy budget requires accurate field measurements for model calibration and validation. Here we were able to collect and make use of many of the few data sets that exist for intensive in-canopy profile time series measurements. We suggest that more intensive field campaigns, with soil water content observations, especially during the winter season, would help in the development of a more reliable parametrization scheme for the within-canopy eddy diffusivity and soil-atmosphere interface conductance. For future model developments, adding an extra soil-atmosphere interface representation such as moss or herbs on the forest floor would be beneficial for a more complete multi-layer energy budget with the objective of describing the surface-atmosphere interface gas and water vapour exchanges.

## 6 Code and data availability

The code and the run environment are open source. Nevertheless, readers interested in running ORCHIDEE-CAN v1.0 (revision 2754) are encouraged to contact the corresponding author for full details and the latest bug fixes. The ORCHIDEE-CAN branch is available via the follow web link: https://forge.ipsl.jussieu.fr/orchidee/browser/branches/ORCHIDEE-DOFOCO/ORCHIDEE.

**The Supplement related to this article is available online at doi:10.5194/gmd-9-2951-2016-supplement.**

*Author contributions.* Yiying Chen, James Ryder, and Sebastiaan Luyssaert developed the parametrization scheme. Yiying Chen, Sebastiaan Luyssaert, and Philippe Peylin designed the study and Yiying Chen wrote the manuscript with contributions from all co-authors. James Ryder, Matthew J. McGrath, Juliane Otto, Kim Naudts, Sebastiaan Luyssaert, and Aude Valade helped Yiying Chen with integrating the parametrization scheme for the multi-layer energy budget in ORCHIDEE-CAN v1.0. Vladislav Bas-trikov and Philippe Peylin provided the optimization tools and helped with the configuration of these tools. Eva van Gorsel, Vanessa Haverd, Bernard Heinesch, Alexander Knohl, Samuli Lau-niainen, Denis Loustau, Eddy Moors, Jérôme Ogée, Thomas Foken, and Timo Vessala provided field observations for all study sites.

*Acknowledgements.* Yiying Chen, James Ryder, Matthew J. McGrath, Juliane Otto, Kim Naudts and Sebastiaan Luyssaert were funded through ERC starting grant 242564 (DOFOCO), and Aude Valade was funded through ADEME (BiCaFF). We thank the two anonymous reviewers for their very helpful observations and suggestions.

Edited by: J. Kala

Reviewed by: two anonymous referees

## References

Aubinet, M., Chermanne, B., Vandenhaute, M., Longdoz, B., Yernaux, M., and Laitat, E.: Long term carbon dioxide exchange above a mixed forest in the Belgian Ardennes, Agr. Forest Meteorol., 108, 293–315, doi:10.1016/S0168-1923(01)00244-1, 2001.

Baldocchi, D.: A Multi-layer model for estimating sulfur dioxide deposition to a deciduous oak forest canopy, Atmos. Environ., 22, 869–884, doi:10.1016/0004-6981(88)90264-8, 1988.

Baldocchi, D. D. and Wilson, K. B.: Modeling CO₂ and water vapor exchange of a temperate broadleaved forest across hourly to decadal time scales, Ecol. Model., 142, 155–184, doi:10.1016/S0304-3800(01)00287-3, 2001.

Barr, A. G., Black, T. A., Hogg, E. H., Kljun, N., Morgenstern, K., and Nesic, Z.: Inter-annual variability in the leaf
area index of a boreal aspen-hazelnut forest in relation to net ecosystem production, Agr. Forest Meteorol., 126, 237–255, doi:10.1029/2002JD003011, 2004.

Bonan, G. B.: A land surface model (LSM version 1.0) for ecological, hydrological, and atmospheric studies, Technical description and user's guide, NCAR Tech. Note NCAR/TN-417+STR, Tech. rep., Boulder, Colorado, USA, 1996.

Bonan, G. B., Oleson, K. W., Fisher, R. A., Lasslop, G., and Reichstein, M.: Reconciling leaf physiological traits and canopy flux data: Use of the TRY and FLUXNET databases in the Community Land Model version 4, J. Geophys. Res.-Biogeo., 117, G02026, doi:10.1029/2011JG001913, 2012.

Bonan, G. B., Williams, M., Fisher, R. A., and Oleson, K. W.: Modeling stomatal conductance in the earth system: linking leaf water-use efficiency and water transport along the soil–plant–atmosphere continuum, Geosci. Model Dev., 7, 2193–2222, doi:10.5194/gmd-7-2193-2014, 2014.

Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C.: A Limited Memory Algorithm for Bound Constrained Optimization, SIAM J. Sci. Comput., 16, 1190–1208, doi:10.1137/0916069, 1995.

Cescatti, A. and Marcolla, B.: Drag coefficient and turbulence intensity in conifer canopies, Agr. Forest Meteorol., 121, 197–206, doi:10.1016/j.agrformet.2003.08.028, 2004.

Chen, Y.-Y. and Li, M.-H.: Determining Adequate Averaging Periods and Reference Coordinates for Eddy Covariance Measurements of Surface Heat and Water Vapor Fluxes over Mountainous Terrain, Terr. Atmos. Ocean. Sci., 23, 685–701, doi:10.3319/TAO.2012.05.02.01(Hy), 2012.

Collatz, G., Ribas-Carbo, M., and Berry, J.: Coupled Photosynthesis-Stomatal Conductance Model for Leaves of C₄ Plants, Aust. J. Plant Physiol., 19, 519–538, doi:10.1071/PP9920519, 1992.

Defraeye, T., Derome, D., Verboven, P., Carmeliet, J., and Nicolai, B.: Cross-scale modelling of transpiration from stomata via the leaf boundary layer, Ann. Bot., 114, 711–723, doi:10.1093/aob/mct313, 2014.

de Noblet-Ducoudré, N., Boisier, J.-P. P., Pitman, A., Bonan, G. B., Brovkin, V., Cruz, F., Delire, C., Gayler, V., van den Hurk, B. J. J. M., Lawrence, P. J., van der Molen, M. K., Müller, C., Reick, C. H., Strengers, B. J., and Voldoire, A.: Determining robust impacts of land-use-induced land cover changes on surface climate over North America and Eurasia: Results from the first set of LU-CID experiments, J. Climate, 25, 3261–3281, doi:10.1175/JCLI-D-11-00338.1, 2012.

Dickinson, R. E., Shaikh, M., Bryant, R., and Graumlich, L.: Interactive Canopies for a Climate Model, J. Climate, 11, 2823–2836, doi:10.1175/1520-0442(1998)011<2823:ICFACM>2.0.CO;2, 1998.

Dolman, A. J. and Wallace, J. S.: Lagrangian and K-theory approaches in modelling evaporation from sparse canopies, Q. J. Roy. Meteorol. Soc., 117, 1325–1340, doi:10.1002/qj.49711750210, 1991.

Dolman, A. J., Moors, E. J., and Elbers, J. A.: The carbon uptake of a mid latitude pine forest growing on sandy soil, Agr. Forest Meteorol., 111, 157–170, doi:10.1016/S0168-1923(02)00024-2, 2002.

Drobinski, P., Anav, A., Lebeaupin Brossier, C., Samson, G., Stéfanon, M., Bastin, S., Baklouti, M., Béranger, K., Beuvier, J., Bourdallé-Badie, R., Coquart, L., D'Andrea, F., de Noblet-Ducoudré, N. A., Diaz, F., Dutay, J. C., Ethe, C., Foujols, M. A., Khvorostyanov, D., Madec, G., Mancip, M., Mason, S., Menut, L., Palmieri, J., Polcher, J., Turquety, S., Valcke, S., and Viovy, N.: Model of the Regional Coupled Earth system (MORCE): Application to process and climate studies in vulnerable regions, Environ. Modell. Softw., 35, 1–18, doi:10.1016/j.envsoft.2012.01.017, 2012.

Ducoudré, N. I., Laval, K., and Perrier, A.: SECHIBA, a New Set of Parameterizations of the Hydrologic Exchanges at the Land-Atmosphere Interface within the LMD Atmospheric General Circulation Model, J. Climate, 6, 248–273, doi:10.1175/1520-0442(1993)006<0248:SANSOP>2.0.CO;2, 1993.

Dufresne, J.-L. and Ghattas, J.: Description du schema de la couche limite turbulente et l'interface avec la surface planetaire dans LMDZ, Tech. rep., LMDZ, Paris, France, 2009.

El Masri, B., Shu, S., and Jain, A. K.: Implementation of a dynamic rooting depth and phenology into a land surface model: Evaluation of carbon, water, and energy fluxes in the high latitude ecosystems, Agr. Forest Meteorol., 211–212, 85–99, doi:10.1016/j.agrformet.2015.06.002, 2015.

Farquhar, G. D., von Caemmerer, S., and Berry, J. A.: A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species, Planta, 90, 78–90, doi:10.1007/BF00386231, 1980.

Foken, T., Meixner, F. X., Falge, E., Zettsch, C., Serafimovich, A., Bargsten, A., Behrendt, T., Biermann, T., Breuninger, C., Dix, S., Gerken, T., Hunner, M., Lehmann-Pape, L., Hens, K., Jocher, G., Kesselmeier, J., Lüers, J., Mayer, J.-C., Moravek, A., Blake, D., Riederer, M., Rütz, F., Scheibe, M., Siebicke, L., Sörgel, M., Staudt, K., Trebs, I., Tsokankunku, A., Welling, M., Wolff, V., and Zhu, Z.: Coupling processes and exchange of energy and reactive and non-reactive trace gases at a forest site – results of the EGER experiment, Atmos. Chem. Phys., 12, 1923–1950, doi:10.5194/acp-12-1923-2012, 2012.

Gao, W., Shaw, R. H., and Paw U, K. T.: Observation of organized structure in turbulent flow within and above a forest canopy, Bound.-Lay. Meteorol., 47, 349–377, doi:10.1007/BF00122339, 1989.

Garratt, J. R.: The Atmospheric Boundary Layer, Cambridge University Press, New York, USA, 1992.

Grace, J.: The turbulent boundary layer over a flapping Populus leaf, Plant Cell Environ., 1, 35–38, doi:10.1111/j.1365-3040.1978.tb00743.x, 1978.

Gu, L., Shugart, H. H., Fuentes, J. D., Black, T., and Shewchuk, S. R.: Micrometeorology, biophysical exchanges and NEE decomposition in a two-story boreal forest – development and test of an integrated model, Agr. Forest Meteorol., 94, 123–148, doi:10.1016/S0168-1923(99)00006-4, 1999.

Guenther, A., Hewitt, C. N., Erickson, D., Fall, R., Geron, C., Graedel, T., Harley, P., Klinger, L., Lerdau, M., Mckay, W. A., Pierce, T., Scholes, B., Steinbrecher, R., Tallamraju, R., Taylor, J., and Zimmerman, P.: A global model of natural volatile organic compound emissions, J. Geophys. Res., 100, 8873, doi:10.1029/94JD02950, 1995.

Guenther, A., Karl, T., Harley, P., Wiedinmyer, C., Palmer, P. I., and Geron, C.: Estimates of global terrestrial isoprene emissions using MEGAN (Model of Emissions of Gases and Aerosols from Nature), Atmos. Chem. Phys., 6, 3181–3210, doi:10.5194/acp-6-3181-2006, 2006.
Haverd, V., Leuning, R., Griffith, D., Gorsel, E. V., and Cuntz, M.: The Turbulent Lagrangian Time Scale in Forest Canopies Constrained by Fluxes, Concentrations and Source Distributions, Bound.-Lay. Meteorol., 130, 209–228, doi:10.1007/s10546-008-9344-4, 2009.

Haverd, V., Lovell, J. L., Cuntz, M., Jupp, D. L. B., Newnham, G. J., and Sea, W.: The Canopy Semi-analytic P gap And Radiative Transfer (CanSPART) model: Formula- tion and application, Agr. Forest Meteorol., 160, 14–35, doi:10.1016/j.agrformet.2012.01.018, 2012.

He, Y., De Wekker, S. F., Fuentes, J. D., and D'Odorico, P.: Coupled land-atmosphere modeling of the effects of shrub encroachment on nighttime temperatures, Agr. Forest Meteorol., 151, 1690– 1697, doi:10.1016/j.agrformet.2011.07.005, 2011.

Hutley, L. B., O'Grady, A. P., and Eamus, D.: Evapotranspiration from eucalypt open-forest savanna of northern australia, Funct. Ecol., 14, 183–194, doi:10.1046/j.1365-2435.2000.00416.x, 2000.

Jackson, R. B.: Trading Water for Carbon with Biolog- ical Carbon Sequestration, Science, 310, 1944–1947, doi:10.1126/science.1119282, 2005.

Jiménez, C., Prigent, C., Mueller, B., Seneviratne, S. I., Mc- Cabe, M. F., Wood, E. F., Rossow, W. B., Balsamo, G., Betts, A. K., Dirmeyer, P. A., Fisher, J. B., Jung, M., Kanamitsu, M., Reichle, R. H., Reichstein, M., Rodell, M., Sheffield, J., Tu, K., and Wang, K.: Global intercomparison of 12 land sur- face heat flux estimates, J. Geophys. Res.-Atmos., 116, D02102, doi:10.1029/2010JD014545, 2011.

Jordan, D. and Smith, W.: Energy balance analysis of nighttime leaf temperatures and frost formation in a subalpine environ- ment, Agr. Forest Meteorol., 71, 359–372, doi:10.1016/0168- 1923(94)90020-5, 1994.

Knohl, A., Schulze, E. D., Kolle, O., and Buchmann, N.: Large carbon uptake by an unmanaged 250-year-old deciduous for- est in Central Germany, Agr. Forest Meteorol., 118, 151–167, doi:10.1016/S0168-1923(03)00115-1, 2003.

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C.: A dynamic global vegetation model for studies of the cou- pled atmosphere-biosphere system, Global Biogeochem. Cy., 19, GB1015, doi:10.1029/2003GB002199, 2005.

Kuppel, S., Peylin, P., Maignan, F., Chevallier, F., Kiely, G., Montagnani, L., and Cescatti, A.: Model–data fusion across ecosystems: from multisite optimizations to global simulations, Geosci. Model Dev., 7, 2581–2597, doi:10.5194/gmd-7-2581- 2014, 2014.

Kustas, W. P., Prueger, J. H., Hatfield, J. L., Ramalingam, K., and Hipps, L. E.: Variability in soil heat flux from a mesquite dune site, Agr. Forest Meteorol., 103, 249–264, doi:10.1016/S0168- 1923(00)00131-3, 2000.

Laitat, E., Chermanne, B., and Portier, B.: Biomass, carbon and ni- trogen allocation in open top chambers under ambient and ele- vated CO₂ and in a mixed forest stand A tentative approach for scaling up from the experiments of Vielsalm, in: Forest Ecosys- tem Modelling, Upscaling and Remote Sensing, 33–59, Aca- demic Publishing, The Hague, the Netherlands, 1998.

Launiainen, S., Vesala, T., Mölder, M., Mammarella, I., Smolander, S., Rannik, Ü., Kolari, P., Hari, P., Lindroth, A., and Katul, G. G.: Vertical variability and effect of stability on turbulence character-

istics down to the floor of a pine forest, Tellus B, 59, 919–936, doi:10.1111/j.1600-0889.2007.00313.x, 2007.

Launiainen, S., Katul, G. G., Lauren, A., and Kolari, P.: Coupling boreal forest CO₂, H₂O and energy flows by a vertically structured forest canopy – Soil model with separate bryophyte layer, Ecol. Model., 312, 385–405, doi:10.1016/j.ecolmodel.2015.06.007, 2015.

Lovell, J., Haverd, V., Jupp, D., and Newnham, G.: The Canopy Semi-analytic Pgap And Radiative Transfer (CanSPART) model: Validation using ground based lidar, Agr. Forest Meteorol., 158– 159, 1–12, doi:10.1016/j.agrformet.2012.01.020, 2012.

Lugo, A. E.: Visible and invisible effects of hurricanes on forest ecosystems: an international review, Austral Ecol., 33, 368–398, doi:10.1111/j.1442-9993.2008.01894.x, 2008.

MacBean, N., Maignan, F., Peylin, P., Bacour, C., Bréon, F.-M., and Ciais, P.: Using satellite data to improve the leaf phenology of a global terrestrial biosphere model, Biogeosciences, 12, 7185– 7208, doi:10.5194/bg-12-7185-2015, 2015.

Maignan, F., Bréon, F.-M., Chevallier, F., Viovy, N., Ciais, P., Gar- rec, C., Trules, J., and Mancip, M.: Evaluation of a Global Veg- etation Model using time series of satellite vegetation indices, Geosci. Model Dev., 4, 1103–1114, doi:10.5194/gmd-4-1103- 2011, 2011.

Makar, P. A., Fuentes, J. D., Wang, D., Staebler, R. M., and Wiebe, H. A.: Chemical processing of biogenic hydrocarbons within and above a temperate deciduous forest, J. Geophys. Res., 104, 3581, doi:10.1029/1998JD00065, 1999.

Massman, W. J. and Weil, J. C.: An analytical one-dimensional second-order closure model of turbulence statistics and the Lagrangian time scale within and above plant canopies of arbitrary structure, Bound.-Lay. Meteorol., 91, 81–107, doi:10.1023/A:1001810204560, 1999.

Mayocchi, C. and Bristow, K.: Soil surface heat flux: some general questions and comments on measurements, Agr. Forest Meteo- rol., 75, 43–50, doi:10.1016/0168-1923(94)02198-S, 1995.

McGrath, M. J., Pinty, B., Ryder, J., Otto, J., and Luyssaert, S.: A multilevel canopy radiative transfer scheme based on a do- mainaveraged structure factor, in preparation, 2016.

McNaughton, K. G. and Van Den Hurk, B. J. J. M.: A “Lagrangian” revision of the resistors in the two-layer model for calculating the energy budget of a plant canopy, Bound.-Lay. Meteorol., 74, 261–288, doi:10.1007/BF00712121, 1995.

Moors, E. J.: Water Use of Forests in the Netherlands, PhD thesis, Wageningen, the Netherlands, 2012.

Naudts, K., Ryder, J., McGrath, M. J., Otto, J., Chen, Y., Valade, A., Bellasen, V., Berhongaray, G., Bönisch, G., Campioli, M., Ghattas, J., De Groote, T., Haverd, V., Kattge, J., MacBean, N., Maignan, F., Merilä, P., Penuelas, J., Peylin, P., Pinty, B., Pretzsch, H., Schulze, E. D., Solyga, D., Vuichard, N., Yan, Y., and Luyssaert, S.: A vertically discretised canopy description for ORCHIDEE (SVN r2290) and the modifications to the energy, water and carbon fluxes, Geosci. Model Dev., 8, 2035–2065, doi:10.5194/gmd-8-2035-2015, 2015.

Ogée, J., Brunet, Y., Loustau, D., Berbigier, P., and Delzon, S.: MuSICA, a CO₂, water and energy multilayer, multi- leaf pine forest model: Evaluation from hourly to yearly time scales and sensitivity analysis, Glob. Change Biol., 9, 697–717, doi:10.1046/j.1365-2486.2003.00628.x, 2003.
Pinty, B., Lavergne, T., Dickinson, R. E., Widlowski, J. L., Gobron, N., and Verstraete, M. M.: Simplifying the interaction of land surfaces with radiation for relating remote sensing products to climate models, J. Geophys. Res.-Atmos., 111, D02116, doi:10.1029/2005JD005952, 2006.

Pitman, A. J.: The evolution of, and revolution in, land surface schemes designed for climate models, Int. J. Climatol., 23, 479–510, doi:10.1002/joc.893, 2003.

Pitman, A. J., De Noblet-Ducoudré, N., Cruz, F. T., Davin, E. L., Bonan, G. B., Brovkin, V., Claussen, M., Delire, C., Ganzeveld, L., Gayler, V., Van Den Hurk, B. J. J. M., Lawrence, P. J., Van Der Molen, M. K., Müller, C., Reick, C. H., Seneviratne, S. I., Strengen, B. J., and Voldoire, A.: Uncertainties in climate responses to past land cover change: First results from the LU-CID intercomparison study, Geophys. Res. Lett., 36, L14814, doi:10.1029/2009GL039076, 2009.

Porte, A., Bosc, A., Champion, I., and Loustau, D.: Estimating the foliage area of Maritime pine (Pinus pinaster Ait.) branches and crowns with application to modelling the foliage area distribution in the crown, Ann. For. Sci., 57, 73–86, doi:10.1051/forest:2000110, 2000.

Prihodko, L., Denning, A., Hanan, N., Baker, I., and Davis, K.: Sensitivity, uncertainty and time dependence of parameters in a complex land surface model, Agr. Forest Meteorol., 148, 268–287, doi:10.1016/j.agrformet.2007.08.006, 2008.

Raupach, M. R.: Applying Lagrangian Fluid-Mechanics To Infer Scalar Source Distribution Concentration Profiles in Plant Canopies, Agr. Forest Meteorol., 47, 85–108, 1989.

Raupach, M. R.: Vegetation-atmosphere interaction in homogeneous and heterogeneous terrain: some implications of mixed-layer dynamics, Vegetatio, 91, 105–120, doi:10.1007/BF00036051, 1991.

Richardson, A. D., Hollinger, D. Y., Burba, G. G., Davis, K. J., Flanagan, L. B., Katul, G. G., Munger, J. W., Ricciuto, D. M., Stoy, P. C., Suyker, A. E., Verma, S. B., and Wofsy, S. C.: A multi-site analysis of random error in tower-based measurements of carbon and energy fluxes, Agr. Forest Meteorol., 136, 1–18, doi:10.1016/j.agrformet.2006.01.007, 2006.

Ryder, J., Polcher, J., Peylin, P., Ottlé, C., Chen, Y., van Gorsel, E., Haverd, V., McGrath, M. J., Naudts, K., Otto, J., Valade, A., and Luyssaert, S.: A multi-layer land surface energy budget model for implicit coupling with global atmospheric simulations, Geosci. Model Dev., 9, 223–245, doi:10.5194/gmd-9-223-2016, 2016.

Scheffers, B. R., Phillips, B. L., Laurance, W. F., Sodhi, N. S., Diesmos, A., and Williams, S. E.: Increasing arboreality with altitude: a novel biogeographic dimension, P. Roy. Soc. B Biol. Sci., 280, 1–9, doi:10.1098/rspb.2013.1581, 2013.

Seidl, R., Fernandes, P. M., Fonseca, T. F., Gillet, F., Jönsson, A. M., Merganičová, K., Netherer, S., Arpaci, A., Bontemps, J.-D., Bugmann, H., González-Olabarria, J. R., Lasch, P., Meredieu, C., Moreira, F., Schelhaas, M.-J., and Mohren, F.: Modelling natural disturbances in forest ecosystems: a review, Ecol. Model., 222, 903–924, doi:10.1016/j.ecolmodel.2010.09.040, 2011.

Sellers, P. J., Los, S. O., Tucker, C. J., Justice, C. O., Dazlich, D. A., Collatz, G. J., and Randall, D. A.: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part II: The generation of global fields of terrestrial biophysical parameters

from satellite data, J. Climate, 9, 706–737, doi:10.1175/1520-0442(1996)009<0706:ARLSPF>2.0.CO;2, 1996.

Staudt, K., Serafimovich, A., Siebicke, L., Pyles, R. D., and Falge, E.: Vertical structure of evapotranspiration at a forest site (a case study), Agr. Forest Meteorol., 151, 709–729, doi:10.1016/j.agrformet.2010.10.009, 2011.

Stöckli, R. and Vidale, P. L.: Modeling diurnal to seasonal water and heat exchanges at European Fluxnet sites, Theor. Appl. Climato., 80, 229–243, doi:10.1007/s00704-004-0102-3, 2005.

Stroud, C., Makar, P., Karl, T., Guenther, A., Geron, C., Turnipseed, A., Nemitz, E., Baker, B., Potosnak, M., and Fuentes, J. D.: Role of canopy-scale photochemistry in modifying biogenic-atmosphere exchange of reactive terpene species: Results from the CELTIC field study, J. Geophys. Res.-Atmos., 110, 149–162, doi:10.1029/2005JD005775, 2005.

Taylor, K. E.: Summarizing multiple aspects of model performance in a single diagram, J. Geophys. Res., 106, 7183, doi:10.1029/2000JD900719, 2001.

Thomas, C. and Foken, T.: Flux contribution of coherent structures and its implications for the exchange of energy and matter in a tall spruce canopy, Bound.-Lay. Meteorol., 123, 317–337, doi:10.1007/s10546-006-9144-7, 2007.

Vuichard, N. and Papale, D.: Filling the gaps in meteorological continuous data measured at FLUXNET sites with ERA-Interim reanalysis, Earth Syst. Sci. Data, 7, 157–171, doi:10.5194/essd-7-157-2015, 2015.

Wang, T., Peng, S., Krinner, G., Ryder, J., Li, Y., Dantec-Nédélec, S., and Ottlé, C.: Impacts of Satellite-Based Snow Albedo Assimilation on Offline and Coupled Land Surface Model Simulations, PLOS ONE, 10, e0137275, doi:10.1371/journal.pone.0137275, 2015.

Wild, M.: How weil do IPCC-AR4/CMIP3 climate models simulate global dimming/brightening and twentieth-century daytime and nighttime warming?, J. Geophys. Res.-Atmos., 114, 1–10, doi:10.1029/2008JD011372, 2009.

Wohlfahrt, G. and Cernusca, A.: Momentum transfer by a mountain meadow canopy: A simulation analysis based on Massman's (1997) model, Bound.-Lay. Meteorol., 103, 391–407, doi:10.1023/A:1014960912763, 2002.

Wolfe, G. M., Thornton, J. A., Bouvier-Brown, N. C., Goldstein, A. H., Park, J.-H., McKay, M., Matross, D. M., Mao, J., Brune, W. H., LaFranchi, B. W., Browne, E. C., Min, K.-E., Wooldridge, P. J., Cohen, R. C., Crounse, J. D., Faloona, I. C., Gilman, J. B., Kuster, W. C., de Gouw, J. A., Huisman, A., and Keutsch, F. N.: The Chemistry of Atmosphere-Forest Exchange (CAFE) Model – Part 2: Application to BEARPEX-2007 observations, Atmos. Chem. Phys., 11, 1269–1294, doi:10.5194/acp-11-1269-2011, 2011.

Yin, X. and Struik, P. C.: $C_3$ and $C_4$ photosynthesis models: An overview from the perspective of crop modelling, NJAS – Wagen. J. Life Sc., 57, 27–38, doi:10.1016/j.njas.2009.07.001, 2009.

Yue, C., Ciais, P., Cadule, P., Thonicke, K., Archibald, S., Poulter, B., Hao, W. M., Hantson, S., Mouillot, F., Friedlingstein, P., Maignan, F., and Viovy, N.: Modelling the role of fires in the terrestrial carbon balance by incorporating SPITFIRE into the global vegetation model ORCHIDEE – Part 1: simulating historical global burned area and fire regimes, Geosci. Model Dev., 7, 2747–2767, doi:10.5194/gmd-7-2747-2014, 2014.