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

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 acceptable 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 observations and long-term surface fluxes were used to calibrate 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 energy partitioning during winter compared with a single-layer model simulation, the increased model complexity does provide a more detailed description of the canopy micrometeorology 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 frequencies, and the consequences of biogenic volatile organic compound (BVOC) emissions from the terrestrial ecosystem.

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 future climates . In such a model, the land surface sub-model provides the surface fluxes to the atmospheric 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 .

Although present-day land surface models differ from each other in their formulation and details, their performances show similar deficiencies. For example, imposing the same land cover changes to seven land surface models resulted 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 . Difficulties in reproducing fluxes of sensible and latent heat for a wide range of vegetation types have been ascribed to the so-called “big-leaf” approach which treats the surface as a isothermal large leaf. Potentially, representing the vertical canopy structure in detail and simulating radiation partitioning and turbulent transport within the vegetation will result in an improved determination of sensible and latent heat flux estimates . For example, several multi-layer SVAT schemes have been proposed and validated with site-level observations . These studies demonstrated that both top-canopy and within-canopy fluxes and micrometeorological profiles could be captured by means of a sophisticated parametrization scheme to describe the vegetation dynamics and the coupling between the atmosphere and the canopy.

Symbolic notation used throughout the paper.

Symbol	Description	Unit
$a_{\mathrm{1}},a_{\mathrm{2}},a_{\mathrm{3}},a_{\mathrm{4}},a_{\mathrm{5}}$	tuning coefficients for $C_{\mathrm{Deff}}$	unitless
$a_{\mathrm{6}}$	factor ceiling of the slope	unitless
$a_{\mathrm{7}}$	critical friction velocity in the middle point of the S-shape function	unitless
$a_{\mathrm{8}}$	factor to constrain the S-shape function	unitless
$a_{\mathrm{9}}$	threshold for vegetation cover	unitless
$a_{10}$	linear weighting factor	unitless
$A$	assimilation rate	$\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$C_{\mathrm{Deff}}$	effective drag coefficient	unitless
$C_{\mathrm{S}}$	concentration of CO${}_{\mathrm{2}}$ at leaf surface	$\mathrm{ppm}$
$C_{\mathrm{D},i}$	vertically discretized estimate for canopy drag coefficient	unitless
$D_{\mathrm{h},\mathrm{air}}$	heat diffusivity of air	${\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$D_{\mathrm{h},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$	heat diffusivity of water vapour	${\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$d_{\mathrm{l}}$	characteristic leaf length	$\mathrm{m}$
$f_{P\mathrm{gap}}$	over-story gap probability from $P$ gap fraction	${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$
$G_{\mathrm{veg}}$	logic variable to indicate the growth status of the vegetation	unitless
$g_{\mathrm{0}}$	residual stomatal conductance if the irradiance approaches zero	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
$h_{\mathrm{s}}$	relative humidity at leaf surface	$\mathrm{\%}$
$h_{\mathrm{c}}$	canopy height	$\mathrm{m}$
$k_i$	diffusivity for level $i$	${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$k_i^{\ast }$	modified diffusivity for level $i$	${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$k_{\mathrm{surf}}$	conductance for the surface–atmosphere interface	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
LAI${}_i$	leaf area index at level $i$	${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$
Nu	Nusselt number	unitless
$P_{\mathrm{m},i}$	momentum shielding factor	unitless
PAI	plant area index	${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$
$R$	correlation coefficient between the simulation and the observation	unitless
$R_{\mathrm{0}}$	maximum correlation coefficient	unitless
$R_{\mathrm{b},i}$	boundary layer resistance at level $i$ for heat	$\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$
$R_{\mathrm{b},i}^{\prime }$	boundary layer resistance at level $i$ for water vapour	$\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$
$R_{\mathrm{s},i}$	stomatal resistance at level $i$	$\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$
Re	Reynold's number	unitless
SLA	specific leaf area	${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$
$S_{\mathrm{T}}$	Taylor skill score	unitless
$T_{\mathrm{week}}$	weekly mean air temperature	$\mathrm{K}$
$T_g$	temperature threshold for under-story phenology	$\mathrm{K}$
$T_{\mathrm{L}}$	Lagrangian timescale	$\mathrm{s}$
$u_{\ast }$	friction velocity	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
$u_i$	velocity at level $i$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
$V_{\mathrm{cmax}}$	carboxylation capacity	$\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
$W_{\mathrm{br}}$	weighting parameter for boundary layer resistance	unitless
$W_{\mathrm{nf}}$	near-field weighting factor	unitless
$W_{\mathrm{sf}}$	weighting parameter for soil–atmosphere conductance	unitless
$W_{\mathrm{sr}}$	linear reduction parameter for stomatal resistance	unitless
${\mathit{\beta }}_{\mathrm{3}}$	fraction of potential plant transpiration realized	unitless
${\mathit{\beta }}_{\mathrm{4}}$	fraction of soil evaporation realized	unitless
$\mathit{\mu }$	kinematic viscosity of air	${\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
${\widehat{\mathit{\sigma }}}_f$	ratio of the variance of the simulations over the variances of observations	unitless
${\mathit{\sigma }}_{\mathrm{w}}$	standard deviation in vertical velocity	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$

Because the standard version of ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems) makes use of a big-leaf approach , improved model capacity and performance were aimed for by implementation of a multi-layer energy budget scheme that was integrated with vertically discrete reflectivity, photosynthesis, stomatal resistance, and carbon allocation schemes. This new design resulted in a new version of ORCHIDEE named ORCHIDEE-CAN (ORCHIDEE-CANopy, revision 2290) . Despite its code including a multi-layer energy budget scheme , ORCHIDEE-CAN is currently applied using a single-layer energy budget, due to a lack of validated parameters for the multi-layer energy budget scheme.

In , 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${}_{\mathrm{2}}$ 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 demonstrate 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 doing 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.

Methodology

Multi-layer energy budget scheme

The multi-layer energy budget scheme used in this study was developed for global land surface models 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; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$) 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 the LAI profile. All symbols are explained in Table . 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 and . 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 neighbouring layers. The shortwave radiation transfer scheme, developed by , 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 .

• b.

  The wind profile and the vertical eddy diffusivity ($k$; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) are calculated using the one-dimensional second-order closure model of , 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_{\mathrm{b}}$; $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) is calculated based upon the leaf boundary-layer resistance, which is estimated according to . The stomatal resistance ($R_{\mathrm{s}}$; $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) is calculated using a Farquhar–von Caemmerer–Berry-type C3 and Collatz-type C4 photosynthesis model 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 to calculate layered stomatal resistances which depend on the ambient air temperature, humidity, within-canopy CO${}_{\mathrm{2}}$ concentration, vegetation-specific maximum carboxylation rate, and water supply from the roots to the stomata.

Readers are referred to 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; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$). The calculation of the vertical and horizontal distributions of the leaf mass, and thus the vegetation canopy, depends on plant phenology, intra-stand competition, forest management, and allometric relationships, and is detailed in .

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 , 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 ) 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.

Stand structure and data availability of the experimental sites. The maximum observed leaf area (LAI; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$) 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_{\mathrm{a}}$ denotes atmospheric temperature, and $q_{\mathrm{a}}$ denotes atmospheric humidity. LE, $H$, and $R_{\mathrm{n}}$ 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.

Site Code	FI-Hyy	FR-LBr	NL-Loo	DE-Bay	CA-Oas	AU-Tum	DE-Hai	BE-Vie
Species	Pinus sylvestris	Pinus pinaster	Pinus sylvestris	Picea abies	Populus sp.	Eucalyptus sp.	Fagus sylvatica	Fagus sylvatica${}^{\mathrm{a}}$
Leaf type	Needleleaved	Needleleaved	Needleleaved	Needleleaved	Broadleaved	Broadleaved	Broadleaved	Broadleaved
Growth form	Evergreen	Evergreen	Evergreen	Evergreen	Deciduous	Evergreen	Deciduous	Mixed
ORCHIDEE PFT	18	5	6	7	20	15	13	13
Over-story LAI	6.5	2.0	1.9	4.8	2.9	2.5	5.8	5.1
Under-story LAI	0.5	1.5	1.5	0.5	2.8	1.0	0.1	0.1
Height	17.0	23.0	15.0	15.0	22.0	50.0	30.0	25.0
$U$ profile	$+$	$-$	$+$	$+$	$+$	$+$	$+$	$+$
$T_{\mathrm{a}}$ profile	$+$	$+$	$+$	$+$	$+$	$+$	$+$	$+$
$q_{\mathrm{a}}$ profile	$+$	$+$	$+$	$+$	$+$	$+$	$-$	$+$
LE profile	$+$	$+$	$+$	$+$${}^{\mathrm{b}}$	$+$	$-$	$-$	$-$
$H$ profile	$+$	$+$	$+$	$+$	$+$	$+$	$+$	$-$
$R_{\mathrm{n}}$ profile	$-$	$+$	$+$	$+$${}^{\mathrm{b}}$	$-$	$-$	$-$	$-$
Reference		,	,	,		,		,

${}^{\mathrm{a}}$ This site is partially mixed with Pseudotsuga menziesii. ${}^{\mathrm{b}}$ The LE profile was available for the 2007 and 2008 period, but not 2011, and the $R_{\mathrm{n}}$ profile was partly available in 2007.

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 ). 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 ). Furthermore, our parametrization and validation set-up required that top-canopy observations had to be available for periods exceeding 1 year (Period IV; Table ). 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 . 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. ) and model experiments that aimed at testing the performance of only the multi-layer energy budget (see EXP1 and EXP3 in Sect. ) made use of the observed LAI profiles. For the remaining two model experiments (see EXP2 and EXP4 in Sect. ), ORCHIDEE-CAN v1.0 calculated the vertical LAI profiles following the carbon allocation and carbon turnover schemes, as described in .

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 , 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.

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 . 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.

Site code	FI-Hyy	FR-LBr	NL-Loo	DE-Bay	CA-Oas	AU-Tum	DE-Hai	BE-Vie
Period for short-term parameter optimiza- tion (Period I)	01/08/06 14/08/06	31/07/06 05/08/06	08/07/97 12/07/97	04/07/11 17/07/11	16/08/94 22/08/94	08/11/06 11/11/06	10/05/01 19/05/01	01/08/02 07/08/02
Closure gap (W m${}^{-\mathrm{2}}$)	43.34	41.56	10.48	18.97	19.82	18.40	29.89	28.19
Period for long-term parameter optimiza- tion (Period II)	01/01/02 31/12/02	01/01/03 31/12/03	01/01/02 31/12/02	01/01/97 31/12/97	01/01/05 31/12/05	01/06/01 31/06/02	01/01/05 31/12/05	01/01/97 31/12/97
Closure gap (W m${}^{-\mathrm{2}}$)	11.47	21.59	15.38	42.47	2.89	7.12	27.83	42.43
Period for single-year EXP1 and EXP3 validation (Period III)	01/01/05 31/12/05	01/01/06 31/12/06	01/01/97 31/12/97	01/01/99 31/12/99	01/01/04 31/12/04	01/06/04 31/06/05	01/01/01 31/12/01	01/01/02 31/12/02
Closure gap (W m${}^{-\mathrm{2}}$)	10.99	13.20	16.61	50.24	4.13	7.73	23.49	42.43
Period for multi-year EXP2 and EXP4 validation (Period IV)	01/01/02 31/12/06	01/01/03 31/12/06	01/01/02 31/12/06	01/01/97 31/12/99	01/01/04 31/12/05	01/06/01 31/06/05	01/01/00 31/12/06	01/01/97 31/12/06
Closure gap (W m${}^{-\mathrm{2}}$)	10.68	17.03	22.65	48.14${}^{\ast }$	3.51	9.40	23.69	33.77

${}^{\ast }$ The forest in 1997–1999 was strongly affected by forest decline; in 2011 the forest was again in a good state.

Description of parameters, code reference, initial values, and tuning ranges used in the multi-layer energy budget model in this work.

Parameter	Physical parameter	Empirical representation of	ORCHIDAS	Default	Tuning range
name			name	value
$a_{\mathrm{1}}$	Effective surface drag	Bending of tree branches	a_1	6.410	Use default
$a_{\mathrm{2}}$	Effective surface drag	Bending of tree branches	a_2	0.001	Use default
$a_{\mathrm{3}}$	Effective surface drag	Bending of tree branches	a_3	0.434	0.1 to 0.8
$a_{\mathrm{4}}$	Effective surface drag	Bending of tree branches	a_4	$-$0.751	$-$0.9 to $-$0.1
$a_{\mathrm{5}}$	Effective surface drag	Bending of tree branches	a_5	0.071	0.05 to 0.1
$a_{\mathrm{6}}$	Eddy diffusivity	Inner-canopy turbulent mixing	k_eddy_slope	5.0	1.0 to 20.0
$a_{\mathrm{7}}$	Eddy diffusivity	Inner-canopy turbulent mixing	k_eddy_ustar	0.3	0.0 to 0.6
$a_{\mathrm{8}}$	Surface–atmosphere conductance	Inner-canopy turbulent mixing	ks_slope	5.0	1.0 to 20.0
$a_{\mathrm{9}}$	Surface–atmosphere conductance	Under-story phenology	ks_veget	0.5	0.0 to 1.0
$a_{10}$	Surface–atmosphere conductance	Under-story phenology	ks_tune	1.0	0.5 to 1.5
$W_{\mathrm{br}}$	Layer boundary resistance	Upscaling the leaf coupling	br_fac	1.0	0.1 to 10.0
$W_{\mathrm{sr}}$	Layer stomatal resistance	Upscaling the leaf coupling	sr_fac	1.0	0.1 to 10.0

Effective drag coefficient $C_{\mathrm{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 that was used to estimate the vertical within-canopy wind profile. In this wind profile model , 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_{\mathrm{1}}$ and $a_{\mathrm{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_{\mathrm{D},i}$; unitless) and the momentum shielding factor ($P_{\mathrm{m},i}$; unitless) as follows: $$C_{\mathrm{Deff},i}=C_{\mathrm{D},i}/P_{\mathrm{m},i}.$$ Both the within-canopy drag and the momentum shielding were parametrized using a function of the cumulative leaf area index (LAI${}_{\mathrm{cum}}$; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$) from the top canopy layer to the bottom layer, which was modified from the original function as below: $$C_{\mathrm{Deff},i}=a_{\mathrm{1}}^{-{\mathrm{LAI}}_{\mathrm{cum},i}/a_{\mathrm{2}}}+a_{\mathrm{3}}^{-{\mathrm{LAI}}_{\mathrm{cum},i}/a_{\mathrm{4}}}+a_{\mathrm{5}},$$ where the subscript $i$ denotes the index of layering from the bottom layer ($i=\mathrm{1}$) to the top-canopy layer ($i=n$). $a_{\mathrm{1}}$ to $a_{\mathrm{5}}$ are tuning coefficients (unitless). The default parameter values for $a_{\mathrm{1}}$ to $a_{\mathrm{5}}$ are presented in Table .

Weighting functions for eddy diffusivity and surface conductance. (a) Weighting function for the eddy diffusivity ($k$) within the air column (Eq. ). The weighting is a function of the friction velocity ($u_{\ast }$) and was optimized by tuning the parameters $a_{\mathrm{6}}$ and $a_{\mathrm{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. ). This weighting function was optimized by tuning the parameters $a_{\mathrm{8}}$ to $a_{10}$. Two examples have the following parameter values: $a_{\mathrm{8}}=10.0$, $a_{\mathrm{9}}=0.5$, $a_{10}=1.0$, $T_{\mathrm{week}}\ge T_g$, and $T_{\mathrm{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_{\mathrm{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_{\mathrm{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.

[Figure omitted. See PDF]

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

After the vertical wind profile was derived from the one-dimensional second-order closure wind profile model, the friction velocity ($u_{\ast }$, $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$), the vertical wind velocity variance (${\mathit{\sigma }}_{\mathrm{w}}$; $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$), and the Lagrangian timescale ($T_{\mathrm{L}}$; $\mathrm{s}$) were calculated following the approach by . In this approach the vertical eddy diffusivity is a function of ${\mathit{\sigma }}_{\mathrm{w}}$ and $T_{\mathrm{L}}$. Subsequently, the vertical eddy diffusivity down the air column to the forest floor was calculated as follows: $$k_i={\mathit{\sigma }}_{\mathrm{w},i}^{\mathrm{2}}{T}_{\mathrm{L},i}.$$

Here we followed the approach proposed by for the Lagrangian timescale calculation. The Lagrangian timescale is thus calculated as $$T_{\mathrm{L},i}=0.66{\displaystyle \frac{\left(\mathrm{1},-,e^{-4.86(z/h_{\mathrm{c}})}\right)}{\left(\mathrm{1},-,e^{-4.86}\right)}}{\displaystyle \frac{h_{\mathrm{c}}}{u_{\ast }}}.$$

A previous effort to validate this model against in situ observations resulted in a bias of the air temperature profile within the canopy layer during nighttime . These issues have been well documented in the scientific literature . 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_{\ast }=u\left(h_{\mathrm{c}}\right)\times (0.32-0.264e^{-15.1\mathit{\zeta }\left(h_{\mathrm{c}}\right)})$, where $\mathit{\zeta }$ is the cumulative function of $C_{\mathrm{Deff}}$ and $h_{\mathrm{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_{\mathrm{nf}}$; unitless). Therefore the modified diffusivity for level $i$ ($k_i^{\ast }$; ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) was defined as $$k_i^{\ast }=W_{\mathrm{nf}}{\mathit{\sigma }}_{\mathrm{w},i}^{\mathrm{2}}{T}_{\mathrm{L},i},$$ where $W_{\mathrm{nf}}$ was calculated as $$W_{\mathrm{nf}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{(-a_{\mathrm{6}}(u_{\ast }-a_{\mathrm{7}}\left)\right)}}}.$$ This function has a sigmoidal shape, where $a_{\mathrm{6}}$ is the ceiling factor of the slope and $a_{\mathrm{7}}$ is the critical friction velocity at the inflection point of the sigmoid function (Fig. a). Consequently, atmospheric diffusivity is reduced if $u_{\ast }$ is low, which represents stable atmospheric conditions. Under turbulent atmospheric conditions, which are represented by a high $u_{\ast }$, $W_{\mathrm{nf}}$ is close to 1 and the simulated diffusivity will closely follow the relationship proposed by . The default parameter values for $a_{\mathrm{6}}$ and $a_{\mathrm{7}}$ are presented in Table . As an alternative to using $u_{\ast }$, 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 . 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).

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

Equation () describes the seasonality of the soil–atmosphere interface, which we believe is driven by the under-story and its phenology . 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_{\mathrm{surf}}$) or, in other words, the calculation of the water vapour exchange between the soil layer and the first air column (see the ${\mathit{\varphi }}_{\mathit{\lambda }E}$ and $K_{\mathrm{surf}}$ in Fig. 1 of , and the formal description of using $K_{\mathrm{surf}}$, which is given in the Supplement of , 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 . 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 ${\mathit{\beta }}_{\mathrm{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 ($\mathrm{1}-f_{P\mathrm{gap}}$), the mean air temperature during the previous week ($T_{\mathrm{week}}$), and a threshold temperature ($T_g$): $${\mathit{\beta }}_{\mathrm{0}}=\left\{\begin{array}{l}\frac{a_{10}}{\mathrm{1}+e^{(-a_{\mathrm{8}}((\mathrm{1}-f_{P\mathrm{gap}})-a_{\mathrm{9}}\left)\right)}},\\ \text{when }{G}_{\mathrm{veg}}=\text{true},\\ \frac{a_{10}}{\left(\mathrm{1},+,e^{(-a_{\mathrm{8}}((\mathrm{1}-f_{P\mathrm{gap}})-a_{\mathrm{9}}\left)\right)}\right)}\frac{T_g-T_{\mathrm{week}}}{T_g-273.15},\\ \text{ when }{G}_{\mathrm{veg}}=\text{false},\end{array}\right.$$ where $a_{\mathrm{8}}$ is a factor that constrains the slope of the function and $a_{\mathrm{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_{\mathrm{veg}}$ is a logic variable to indicate the growth status of the vegetation. $G_{\mathrm{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_{\mathrm{veg}}$ is set to true when the weekly averaged soil water content and temperature exceed the thresholds. $f_{P\mathrm{gap}}$ 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_{\mathrm{surf}}=\left(W_{\mathrm{sf}},{\mathit{\beta }}_{\mathrm{3}},+,(,\mathrm{1},-,W_{\mathrm{sf}},),{\mathit{\beta }}_{\mathrm{4}}\right)\left(u_{\mathrm{1}},C_{\mathrm{Deff},\mathrm{1}}\right),$$ where ${\mathit{\beta }}_{\mathrm{3}}$ and ${\mathit{\beta }}_{\mathrm{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 and . $u_{\mathrm{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_{\mathrm{Deff}}$ is the effective drag coefficient calculated according to Eq. (). $W_{\mathrm{sf}}$ is the weighting factor for the soil–atmosphere interface, which is described as the conditional function of over-story canopy cover fraction ($\mathrm{1}-f_{P\mathrm{gap}}$). $W_{\mathrm{sf}}={\mathit{\beta }}_{\mathrm{0}}$ when $(\mathrm{1}-f_{P\mathrm{gap}})>a_{\mathrm{9}}$, and $W_{\mathrm{sf}}=\mathrm{1}-{\mathit{\beta }}_{\mathrm{0}}$ when $(\mathrm{1}-f_{P\mathrm{gap}}$ $\le a_{\mathrm{9}}$ (see Fig. b). The default parameter values of $a_{\mathrm{8}}$, $a_{\mathrm{9}}$, $a_{10}$, and $W_{\mathrm{sf}}$ are presented in Table .

Boundary-layer resistance of the leaf surface $R_{\mathrm{b}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$)

The boundary-layer resistance of the leaf surface $R_{\mathrm{b},i}$ is described according to the expression from : $$R_{\mathrm{b}}=\left\{\begin{array}{l}{W}_{\mathrm{br}}\left(\frac{d_{\mathrm{l}}}{D_{\mathrm{h},\mathrm{air}}\mathit{\text{Nu}}}\right)\text{, for sensible heat},\\ W_{\mathrm{br}}\left(\frac{d_{\mathrm{l}}}{D_{\mathrm{h},{\mathrm{H}}_{\mathrm{2}}\mathrm{0}}\mathit{\text{Sh}}}\right)\text{, for latent heat},\end{array}\right.$$ where $W_{\mathrm{br}}$ accounts for the fact that the leaf length of the species under study differs from the characteristic leaf length (unitless), $d_{\mathrm{l}}$ is the characteristic leaf length (0.001 m was used as the default value), $D_{\mathrm{h},\mathrm{air}}$ is the heat diffusivity of still air (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), $D_{\mathrm{h},{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the heat diffusivity of water vapour (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), Sh is the Sherwood number (unitless), and Nu is the Nusselt number (unitless). The Sherwood number was calculated as $\mathit{\text{Sh}}=0.66{\mathit{\text{Re}}}^{0.5}{\mathit{\text{Sc}}}^{0.33}$ for laminar flow and $\mathit{\text{Sh}}=0.03{\mathit{\text{Re}}}^{0.8}{\mathit{\text{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 $\mathit{\text{Nu}}=0.66\mathit{\text{Re}}{\mathit{\text{Pr}}}^{0.33}$, where Pr is the Prandtl number (0.7 for air; unitless) and Re is the Reynolds number (unitless) which was calculated as $$\mathit{\text{Re}}={\displaystyle \frac{d_{\mathrm{l}}{u}_i}{\mathit{\mu }}},$$ where $u_i$ is the horizontal velocity at level $i$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and $\mathit{\mu }$ is the kinematic viscosity of air and was set to 0.0015 (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) . The default parameter value for $W_{\mathrm{br}}$ is provided in Table .

Stomatal resistance $R_{\mathrm{s}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$)

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 ). 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 ($R_{\mathrm{b}}$) and stomatal resistance ($R_{\mathrm{s}}$) (see Eqs. 9 and 13 in ). 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 $R_{\mathrm{s}}$ 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 $W_{\mathrm{sr}}$: $$R_{\mathrm{s},i}=W_{\mathrm{sr}}\left({\displaystyle \frac{\mathrm{1}}{\left(g_{\mathrm{0}},+,\left(\frac{A_i{h}_{\mathrm{s}}}{C_{\mathrm{s}}}\right)\right){\mathrm{LAI}}_i}}\right),$$ where $g_{\mathrm{0}}$ is the residual stomatal conductance if the solar irradiance approaches zero, $C_{\mathrm{s}}$ is the concentration of CO${}_{\mathrm{2}}$ at the leaf surface, and $h_{\mathrm{s}}$ is the relative humidity at the leaf surface. $A$ is the CO${}_{\mathrm{2}}$ assimilation rate which is solved analytically following . In Eq. () 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 $W_{\mathrm{sr}}$ is presented in Table .

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 $a_{\mathrm{3}}$ to $a_{10}$, $W_{\mathrm{br}}$, and $W_{\mathrm{sr}}$ 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 . The range assigned to each parameter is reported in Table . 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 . 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 ) were used to optimize $a_{\mathrm{3}}$ to $a_{\mathrm{7}}$, $W_{\mathrm{br}}$, and $W_{\mathrm{sr}}$. During this step, the parameters for the soil–atmosphere interface ($k_{\mathrm{surf}}$, i.e. $a_{\mathrm{8}}$ to $a_{10}$ and $W_{\mathrm{sf}}$) 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 ($C_{\mathrm{Deff}}$; $k$; $R_{\mathrm{b}}$; $R_{\mathrm{s}}$) were biased towards the summer. Secondly, the seasonal dynamics of $k_{\mathrm{surf}}$ was parametrized by trying to improve the correspondence between the simulated and observed top-canopy fluxes over 1 year (Period II in Table ). In this step, $a_{\mathrm{3}}$ to $a_{\mathrm{7}}$, $W_{\mathrm{br}}$, and $W_{\mathrm{sr}}$ were set to the values obtained from the first step of the optimization and $a_{\mathrm{8}}$ to $a_{10}$ and $W_{\mathrm{sf}}$ 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 ).

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.

Attribution of changes in model performance

The multi-layer energy budget scheme 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 . 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 canopyThe 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 ). In this experiment, the vertical LAI profile was only 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 canopyThe 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 canopyExperiment 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 canopyIn 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.

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 $\le$ 3.0).

[Figure omitted. See PDF]

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${}^{\circ }$ $\times$ 0.5${}^{\circ }$ 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 , 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 . 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. ). If the vertical LAI profile was not imposed, the LAI generated for the albedo calculation 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 .

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_{\mathrm{T}}$) .

$S_{\mathrm{T}}$ was calculated for the eight observational sites for the top-canopy fluxes of all four experiments making use of the simulated and observed half-hourly fluxes. The Taylor skill score was calculated as follows: $$S_{\mathrm{T}}={\displaystyle \frac{\mathrm{4}(\mathrm{1}+R)}{{\left({\widehat{\mathit{\sigma }}}_f,+,\mathrm{1},/,{\widehat{\mathit{\sigma }}}_f\right)}^{\mathrm{2}}\left(\mathrm{1},+,R_{\mathrm{0}}\right)}},$$ where $R$ is the correlation coefficient between the simulation and the observation, $R_{\mathrm{0}}$ is the maximum correlation coefficient and ${\widehat{\mathit{\sigma }}}_f$ is the ratio of the variance of the simulations to the variance of observations (${\widehat{\mathit{\sigma }}}_f=\mathit{\sigma }/{\mathit{\sigma }}_r$). Here, we set $R_{\mathrm{0}}$ to 1.0 for the maximum correlation between observation and model simulation. A value of 1.0 of $S_{\mathrm{T}}$ indicates that model simulations perfectly match the observations; values lower than 0.5 imply that the model has poor predictive ability.

Inter-annual and diurnal performance for both dense and sparse forest types, expressed as Taylor skill score ($S_{\mathrm{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_{\mathrm{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.

[Figure omitted. See PDF]

Results

Model parametrization

Using the default parameter set (i.e. $a_{\mathrm{1}}$ to $a_{\mathrm{5}}$) 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 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ after adjusting the parameters (see Table S3 in the Supplement). Tuning the conductance of the soil–atmosphere interface (i.e. $a_{\mathrm{8}}$ to $a_{10}$), 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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for LE and 123.2 to 50.3 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for $H$, respectively (see also Table S3 in the Supplement).

Inter-annual and diurnal performance for both dense and sparse forest types, expressed as Taylor skill score ($S_{\mathrm{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_{\mathrm{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.

[Figure omitted. See PDF]

At sites with sparse canopies, the net radiation at the forest floor was substantial, i.e. ranging nearly from 200 to 450 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Fig. S4 in the Supplement). Correctly simulating radiation transfer strongly contributed to correctly simulating 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 errors of net radiation estimation depended solely on the tree species. In sparse canopies, a positive air temperature gradient with higher temperatures at the forest floor compared to the top-canopy was also presented (Fig. S5 in the Supplement). 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_{\mathrm{3}}$ to $a_{10}$, $W_{\mathrm{br}}$, and $W_{\mathrm{sr}}$) had a large impact on the absolute values of the vertical profile in leaf temperature (Fig. S6 in the Supplement). Unfortunately, leaf temperature was not measured at any of the sites. Therefore, it remains to be assessed whether the model can concurrently reproduce observed energy 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 estimates was tested. Both the stomatal resistance and the boundary resistance weighting factors ($W_{\mathrm{sr}}$ and $W_{\mathrm{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_{\mathrm{6}}$ and $a_{\mathrm{7}}$, which tuned the eddy diffusivity, were largely unaffected by the initial conditions. Soil water content measurements would thus have helped to improve the parametrization, especially for the stomatal and leaf boundary-layer resistances.

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. a) was slightly higher and thus better than model performance at the dense forest sites (Fig. b). 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_{\mathrm{n}}$ were well simulated, displaying $S_{\mathrm{T}}$ scores continuously over 0.9 (Fig. b and d). For dense canopies, the $S_{\mathrm{T}}$ score of $R_{\mathrm{n}}$ dropped to 0.9 in winter, which might be attributed to an incorrect estimation of $R_{\mathrm{n}}$ during nighttime (Fig. c).

In general, the $S_{\mathrm{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_{\mathrm{T}}$ dropped below 0.5 for latent heat flux and 0.8 for sensible heat flux (Fig. a) 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. c and e). During these months nights are long and the inability of the model to simulate nighttime fluxes (Fig. c) 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. e).

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: $\mathit{\square }$, DE-Bay: $○$, DE-Hai: $\mathit{▵}$, and BE-Vie: $\mathit{◊}$; in the second column sub-plot b, d, f, and h are sparse forest sites; FR-LBr: $\mathit{\square }$, NL-Loo: $○$, CA-Oas: $\mathit{▵}$, and AU-Tum: $\mathit{◊}$). Simulation results are symbols filled with colours: EXP1 with grey colour and EXP3 with black colour.

[Figure omitted. See PDF]

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. a, b and e, f) due to small improvements in the $S_{\mathrm{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. c, c and e) were reflected in the improved partitioning of energy fluxes during wintertime (compare Figs. a and a). 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. f).

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. ; 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_{\mathrm{n}}$, $H$, and LE but increased for $G$ (see Fig. g and h, and ). Note that $G$ is an essential aspect in simulating the snow phenology . Therefore, improved simulations of the soil heat fluxes could have important indirect effects on climate simulations of regions with a pronounced snow season.

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: $\mathit{\square }$, DE-Bay: $○$, DE-Hai: $\mathit{▵}$, and BE-Vie: $\mathit{◊}$; in the second column sub-plot b, d, f, and h are sparse forest sites; FR-LBr: $\mathit{\square }$, NL-Loo: $○$, CA-Oas: $\mathit{▵}$, and AU-Tum: $\mathit{◊}$). Simulation results are symbols filled with colours: EXP1 with grey colour and EXP3 with dark-grey colour.

[Figure omitted. See PDF]

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. b 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 underneath 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 . However, the small loss (all fluxes except $G$) or gain (only for $G$) in model skill from introducing the multi-layer scheme can be strengthened (i.e. LE) or compensated for ($R_{\mathrm{n}}$, $H$, and $G$) by the small gain in model skill from the introduction of a prognostic vertical LAI profile.

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.

[Figure omitted. See PDF]

Discussion

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 , which applied the single-layer energy budget to simulate 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 , 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. ).

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.

[Figure omitted. See PDF]

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 ), 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.

• Within-canopy dragFor the inner-canopy drag parametrization, we modified an approach (Eq. ) that has previously only been tested and validated at grassland sites . 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 . 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.

• 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 . As previous studies have demonstrated, incorrect estimation on gradients may be accommodated to some extent by introducing a scaling factor (Eq. ) to constrain the within-canopy transport . 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 .

  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_{\mathrm{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 .

• 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 . 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 environmental conditions in the under-story.

  Starting from the point of view of the interaction between ecosystems and the climate, we introduced a weighting factor ($W_{\mathrm{sf}}$) as a function of a long-term average temperature and light conditions (gap fraction), a transpiration fractions described as ${\mathit{\beta }}_{\mathrm{3}}$ in the model code and a soil evaporation fraction (${\mathit{\beta }}_{\mathrm{4}}$) as environmental factors to parametrize surface conductance (Fig.  and consequently control the surface latent heat flux. This approach demonstrated the model's capability to simulate 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 , which is an environmental condition not accounted for in our approach. Furthermore, accounting for ecosystem-specific differences in root density profiles and aerial cover of the under-story might also help in the simulation of water and energy fluxes . From this perspective, detailed soil moisture profile observations would be very useful in developing a more advanced surface–atmosphere interface parametrization.

• Mismatch between low-resolution driver data and vertically resolved vegetation layersIn this study an apparent mismatch was present between the low resolution of the driver data that contain information derived from several different land cover types and the highly resolved vertical layering of the canopy. When low-resolution driver data are used, the benefit 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.

• The proposed parametrization approach and future workIn general, we provide a simple but useful parametrization approach for the multi-layer energy budget scheme in the ORCHIDEE-CAN v1.0 global land surface model. Comparing with other studies , our approach directly determines the energy and water fluxes and successfully avoids the iterative processes to meet the numerical requirement. In total, a set of 12 parameters need to be prescribed and calibrated regarding the empirical representation of surface drag, turbulent mixing, sub-canopy phenology, and leaf–atmosphere coupling processes. Our approach presents a good performance at all study sites, though we may have some deficits in wind speed estimation.

  In this study the model had been tested for several environmental conditions and demonstrated that the numerics 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, evapotranspiration. Only then will we be able to learn about the transferability of the parameters from the site level to the PFT level.

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.

• 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 .

• Natural disturbances, such as fires, pests, and windfall, can result in increases in leaf fall, individual tree mortality, or complete stand destruction , 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.

• Forest canopy structure plays an important role in regulating the provision of forest ecosystem services such as maintaining biodiversity or regulating streamflow . 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.

• 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 . To run at a regional or global scale, it is essential to first parametrize the model at the site level.

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.

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.

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 Bastrikov 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 Launiainen, 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