A multi-layer land surface energy budget model for implicit coupling with global atmospheric simulations

In Earth system modelling, a description of the energy budget of the vegetated surface layer is fundamental as it determines the meteorological conditions in the planetary boundary layer and as such contributes to the atmospheric conditions and its circulation. The energy budget in most Earth system models has been based on a ‘big-leaf approach’, with averaging schemes that represent in-canopy processes. Furthermore, to be stable, that is to say, over large time steps and without large 5 iterations, a surface layer model should be capable of implicit coupling to the atmospheric model. Surface models with large time steps, however, have difficulties in reproducing consistently the energy balance in field observations. We here outline a newly developed numerical model for energy budget simulation, as a component of the land surface model ORCHIDEE-CAN (Organising Carbon and Hydrology In Dynamic Ecosystems CANopy). This new model implements techniques from 10 single-site canopy models in a practical way. It includes representation of in-canopy transport, a multilayer longwave radiation budget, height-specific calculation of aerodynamic and stomatal conductance, and interaction with the bare soil flux within the canopy space. Significantly, it avoids iterations over the height of the canopy and so maintains implicit coupling to the atmospheric model LMDz (Laboratoire de Météorologie Dynamique Zoomed model). As a first test, the model is eval15 uated against data from both an intensive measurement campaign and longer term eddy covariance measurements for the intensively studied Eucalyptus stand at Tumbarumba, Australia. The model performs well in replicating both diurnal and annual cycles of energy and water fluxes, as well as the vertical gradients of temperature and of sensible heat fluxes.


Introduction
Earth system models are the most advanced tools to predict future climate (Bonan, 2008).These models represent the interactions between the atmosphere and the surface beneath, with the surface formalised as a combination of open oceans, sea-ice, and land.For land, a description of the energy budget of the vegetated surface layer is fundamental as it determines the meteorological conditions in the planetary boundary layer and, as such, contributes to the atmospheric conditions and its circulation.
The vegetated surface layer of the Earth is subject to incoming and outgoing fluxes of energy, namely atmospheric sensible heat (H , Wm −2 ), latent heat (λE, Wm −2 ), shortwave radiation from the sun (R SW , Wm −2 ), long-wave radiation (R LW , Wm −2 ) emitted from other radiative sources such as clouds and atmospheric compounds, and soil heat exchange with the subsurface (J soil , Wm −2 ).The sum of these fluxes is equal to the amount of energy that is stored or released from the surface layer over a given time period t (s).So, for a surface of overall heat capacity C p (J K −1 m −2 ) the Published by Copernicus Publications on behalf of the European Geosciences Union.

J. Ryder et al.: A multi-layer land surface energy budget model
temperature change over time, T , is described as The sign convention used here makes all upward fluxes positive.So a positive sensible or latent heat flux from the surface cools the ground.Likewise a positive radiation flux towards the surface warms the ground.
One key concept in modelling the energy budget of the surface Eq. ( 1) is the way in which the surface layer is defined.In many cases the surface layer describes both the soil cover and the vegetation above it as a uniform block.Such an approach is known as a big leaf model, because the entirety of the volume of the trees or crops and the understorey, as well as the surface layer, are simulated in one entity, to produce fluxes parametrised from field measurements.In the model under study, ORCHIDEE-CAN (Organising Carbon and Hydrology In Dynamic Ecosystems -CANopy) (Naudts et al., 2015), the land surface is effectively simulated as an infinitesimal surface layer -a conceptual construct of zero thickness.As demonstrated in the original paper describing this model, such an approach, whilst reducing the canopy to simple components, was nevertheless able to simulate surface fluxes to an acceptable degree of accuracy for the sites that were evaluated as the original SECHIBA (Schematic of Hydrological Exchange at the Biosphere to Atmosphere Interface) model (Schulz et al., 2001) and later as a component of the original ORCHIDEE model (Krinner et al., 2005), the basis of ORCHIDEE-CAN (revision 2966).
The proof that existing land to surface simulations may now be inadequate comes from inter-comparison studies, such as Pitman et al. (2009), which evaluated the response of such models to land use change scenarios.That study found a marked lack of consistency between the models, an observation they attributed to a combination of the varying implementation of land cover change (LCC) maps, the representation of crop phenology, the parametrisation of albedo, and the representation of evapotranspiration for different land cover types.Regarding the latter two issues, the models they examined did not simulate in a transparent, comparable manner the changes in albedo and evapotranspiration as a result of changes in vegetation cover, such as from forest to cropland.It was not possible to provide a definitive description of the response of latent heat flux to land cover change across the seven models under study, because there was substantial difference in the mechanisms, which describe the evaporative response to the net radiation change across the conducted simulations.
Furthermore, the latent and sensible heat fluxes from offline land surface models were reported to depend very strongly on the process-based parametrisation, even when forced with the same micro-meteorological data (Jiménez et al., 2011).The structure of land surface models, it has been suggested (Schlosser and Gao, 2010), may be more important than the input data in simulating evapotranspiration.
Hence, improvements to the soil-surface-atmosphere interaction (Seneviratne et al., 2010), and to the hydrology (Balsamo et al., 2009), are considered essential for better simulating evapotranspiration.We can, therefore, assert that refinements in the numerical schemes of land surface models represent a logical approach to the further constraint of global energy and water budgets.
Large-scale validation has revealed that the big-leaf approach has difficulties in reproducing fluxes of sensible and latent heat (Jiménez et al., 2011;Pitman et al., 2009;de Noblet-Ducoudré et al., 2012) for a wide range of vegetated surfaces.This lack of modelling capability is thought to be due to the big-leaf approach not representing the vertical canopy structures in detail and thus not simulating factors such as radiation partition, separation of height classes, turbulent transport within the vegetation, and canopy-atmosphere interactions -all of which are crucial factors in the improved determination of sensible and latent heat-flux estimates (Baldocchi and Wilson, 2001;Ogée et al., 2003;Bonan et al., 2014), as well as the presence of an understorey, or mixed canopies, as is proposed by Dolman (1993).Furthermore, a model that is able to determine the temperatures of elements throughout the canopy profile will provide for a more useful comparison with remote sensing devices, for which the "remotely sensed surface temperature" (Zhao andQualls, 2005, 2006) also depends on the viewing angle.
This gap in modelling capability provides the motivation for developing and testing a new, multi-layer, version of the energy budget simulation based on Eq. (1).A multi-layer approach is expected to model more subtle but important differences in the energy budget in relation to multi-layer vegetation types such as forests, grasses and crops.Through the simulation of more than one canopy layer, the model could simulate the energy budget of different plant types in two or more layers such as that found in savannah, grassland, wood species, and agro-forestry systems (Verhoef and Allen, 2000;Saux-Picart et al., 2009) Where stand-alone surface models have few computational constraints, the typical applications of an Earth system model (ESM) require global simulations at a spatial resolution of 2 • × 2 • or a higher spatial resolution for century long timescales.Such applications come with a high computational demand that must be provided for by using a numerical scheme that can run stably over longer time steps (∼ 15 to 30 min), and that can solve a coupled or interdependent set of equations without iterations.In numerics, such a scheme is known as an implicit solution, and requires that all equations in the coupled systems are linearised.Given that ORCHIDEE is the land surface model of the IPSL (Institute Pierre Simon Laplace) ESM, the newly developed multi-layer model was specifically designed in a numerically implicit way.

Model requirements
Several alternative approaches to the big-leaf model have been developed.These alternatives share the search for a more detailed representation of some of the interactions between the heat and radiation fluxes and the surface layer.Following Baldocchi and Wilson (2001), the range and evolution of such models includes 1. the big-leaf model (e.g.Penman and Schofield, 1951); 2. the big-leaf with dual sources (e.g.Shuttleworth and Wallace, 1985); 3. two-layer models, which split the canopy from the soil layer (e.g.Dolman, 1993;Verhoef and Allen, 2000;Yamazaki et al., 1992); 4. three-layer models, which split the canopy from the soil layer, and simulate the canopy as a separate understorey and overstorey (e.g.Saux-Picart et al., 2009); 5. one-dimensional multi-layer models (e.g.Baldocchi and Wilson, 2001); 6. three-dimensional models that consist of an array of plants and canopy elements (e.g.Sinoquet et al., 2001).
For coupling to an atmospheric model (see below), and thus running at a global scale, simplicity, robustness, generality, and computational speed need to be balanced.We, therefore, propose a one-dimensional multi-layer model combined with a detailed description of the three-dimensional canopy characteristics.We aim for a multi-layer canopy model that -simulates processes that are sufficiently well understood at a canopy level such that they can be parametrised at the global scale through (semi-)mechanistic, rather than empirical, techniques, e.g. the description of stomatal conductance (Ball et al., 1987;Medlyn et al., 2011), and the partition of radiation in transmitted, reflected, and absorbed radiation at different canopy levels (Pinty et al., 2006;McGrath et al., 2016); -simulates the exposure of each section of the canopy, and the soil layer, to both shortwave and long-wave radiation, and simulate in-canopy gradients, separating between soil-surface-atmosphere and vegetationatmosphere interactions; -simulates non-standard canopy set-ups, for instance combining different species in the same vertical structure, e.g.herbaceous structures under trees, as explored by Dolman (1993); Verhoef and Allen (2000); Saux-Picart et al. (2009); -describes directly the interaction between the soil surface and the sub-canopy using an assigned soil resistance rather than a soil-canopy amalgamation; -is flexible, i.e. sufficiently stable to be run over fifty layers or over just two; -avoids introducing numerics that would require iterative solutions.
Where the first five requirements relate to the process description of the multi-layer model, the last requirement is imposed by the need to couple ORCHIDEE to an atmospheric model.Generally, coupling an implicit scheme will be more stable than an explicit scheme, which means that it can be run over longer time steps.Furthermore, the approach is robust: for example, if there is an instability in the land surface model, it will tend to be dampened in subsequent time steps, rather than diverge progressively.For this work, the model needs to be designed to be run over time steps as long as 30 min in order to match the time steps of the IPSL atmospheric model LMDz (Laboratoire de Météorologie Dynamique Zoomed model), to which it is coupled, and therefore to conserve processing time.However, the mathematics of an implicit scheme have to be linearised and is thus by necessity rigidly and carefully designed.As discussed in Polcher et al. (1998) and subsequently in (Best et al., 2004), the use of implicit coupling was widespread in models when the land surface was a simple bucket model, but as the land surface schemes have increased in complexity, explicit schemes have, for most models, been used instead, because complex explicit schemes are more straightforward to derive than implicit schemes.As they demonstrate, there is nevertheless a framework for simulating all land surface fluxes and processes (up to a height of, say, 50 m, so including abovecanopy physics) in a tiled non-bucket surface model coupled, using an implicit scheme, to an atmospheric model.

Model description
We here summarise the key components of the new implicit multi-layer energy budget model.The important innovation, compared to existing multi-layer canopy models that work at the local scale (e.g.Baldocchi, 1988;Ogée et al., 2003), is that we will solve the problems implicitly; i.e. all variables are described in terms of the next time step.The notation used here is listed in full in Table 1, and is chosen to complement the description of the LMDz coupling scheme, as is described in Polcher et al. (1998).A complete version of the derivation of the numerical scheme is provided in the Supplement.
We propose to regard the canopy as a network of potentials and resistances, as shown in Fig. 1, a variation of which was first proposed in Waggoner et al. (1969).At each level in the network, we have the state variable potentials: the temperature of the atmosphere at that level, the atmospheric humidity, and the leaf level temperature.We include in the network fluxes of latent heat and sensible heat between the Symbol Description a 1 , a 2 , . ..a 5 coefficients for canopy turbulence A q,i , B q,i , C q,i , D q,i components for substituted Eq. (ii) A T ,i , B T ,i , C T ,i , D T ,i components for substituted Eq. (i) A i , B i , C i , D i , E i matrix substitutions (from the alternative derivation) C air p specific heat capacity of air (J (kg K) Reynold's number (−) q a i atmospheric-specific humidity at level i (kg kg −1 ) q leaf,i leaf-specific humidity at level i (kg kg −1 ) q T leaf sat saturated-specific humidity of leaf at level at i (kg kg −1 ) q t specific humidity (kg kg leaves at each level and the atmosphere, and vertically between each canopy level.The soil surface interacts with the lowest canopy level, and uppermost canopy level interacts with the atmosphere.We also consider the absorption and reflection of radiation by each vegetation layer and by the surface (SW and LW) and emission of radiation (LW only).This represents the classic multi-layer canopy model formulation, with a network of resistances that simulate the connection between the soil surface temperature and humidity, and fluxes passing through the canopy to the atmosphere.The analogy is the circuit diagram approach, for which T a and q a represent the atmospheric potentials of temperature and specific humidity at different heights, and H and λE are the sensible and latent heat fluxes that act as currents for these potentials.At each level within the vegetation, T a and q a interact with the leaf level temperature and humidity T L and q L through the resistances R i (for resistance to sensible heat flux) and R i (for resistance to latent heat flux).The change in leaf level temperature is determined by the energy balance at each level.The modelling approach formalises the following constraints and assumptions.

Leaf vapour pressure assumption
We assume that the air within leaf level cavities is completely saturated.This means that the vapour pressure of the leaf can be calculated as the saturated vapour pressure at that leaf temperature (Monteith and Unsworth, 2008).Therefore, the change in pressure within the leaf is assumed proportional to the difference in temperature between the present time step and the next one, multiplied by the rate of change in saturated pressure against temperature.The symbolic notation of the subsequent equations is fully explained in Table 1.
where α i and β i are regarded as constants for each particular level and time step; therefore, and .
To find a solution we still need to find an expression for the terms q T t leaf,i sat and ∂q sat ∂T | T t leaf,i in α i and β i above.Using the empirical approximation of Tetens (e.g.Monteith and Unsworth, 2008, Sect. 2.1) and the specific humidity vapour pressure relationship, we can describe the saturation vapour pressure to within 1 Pa up to a temperature of about 35 • C.
So the specific humidity of the leaf follows a relationship to the leaf temperature that is described by a saturation curve.

Derivation of the leaf layer resistances (R i and R i )
The variables R i and R i represent, in our circuit diagram analogue, resistances to the sensible and latent heat flux, respectively.The resistance to the sensible heat flux is equal to the boundary-layer resistance, R b,i , of the leaf surface: For sensible heat flux, R b,i is calculated as (e.g.Monteith and Unsworth, 2008) for which D h,air is the thermal diffusivity of air and d l is the characteristic leaf length.The Nusselt number, N u, is calculated as in Grace and Wilson (1976), for which Nu = 0.66Re 0.5 P r 0.33 , where P r is the Prandtl number (which is 0.70 for air), and Re is the Reynolds number, for which where µ is the kinematic viscosity of air (i.e.0.15 cm 2 s −1 ), d l is again the characteristic dimension of the leaf, and u is the wind speed at the level i in question.
The resistance to latent heat flux is calculated as the sum of the boundary-layer resistance (which is calculated slightly differently) and the leaf stomatal resistance: In this case we use the following expression (e.g.Monteith and Unsworth, 2008) is the molecular diffusivity of water vapour and Sh is the Sherwood number, which for laminar flow is and for turbulent flow is for which Sc is the Schmidt number (which is 0.63 for water; Grace and Wilson, 1976).The transition from laminar to turbulent flow takes place in the model when the Reynolds number exceeds a value of 8000 (Baldocchi, 1988).
The stomatal conductance, g s,i is calculated according to the Ball-Berry approximation, per level i.In summary where g 0 is the residual stomata conductance, A the assimilation rate, h s the relative humidity at the leaf surface and C s the concentration of CO 2 at the leaf surface.
The description here is related to that of the standard OR-CHIDEE model (e.g.LSCE/IPSL, 2012, Sect.2.1), for which the g s that is used to determine the energy budget is calculated as an amalgamated value, over the sum of all levels i.However, in this new energy budget description we keep separate the g s for each level i, and use the inverse of this conductance value to determine the resistance that is R s,i .Furthermore, the amount of water that is supplied to the plant and transported through the plant is calculated (Naudts et al., 2015).In times of drought, the water supply term may be lower than the theoretical latent flux that can be emitted for a certain g s , using Eq. ( 29).In these cases, the g s term at leaf level is restricted to that corresponding to the supply term limited latent heat flux at the level in question.

Leaf interaction with precipitation
Both soil interactions and leaf level evaporation components are parametrised using the same interception and evaporation coefficients as are used in the existing ORCHIDEE model (Krinner et al., 2005; LSCE/IPSL, 2012), extended by ORCHIDEE-CAN.Notably, ORCHIDEE-CAN assumes horizontal clumping of plant species, and hence canopy gaps, as opposed to the uniform medium that is applied in the original ORCHIDEE.A portion of rainfall is intercepted by the vegetation (i.e. a canopy interception reservoir), as determined by the total canopy leaf area index (LAI) and by the plant functional type (PFT), where it will be subject to evaporation as standing water.The rest falls on the soil surface, and is treated in the same way as that for soil water in the existing model.

The leaf energy balance equation for each layer
For vegetation, we assume the energy balance is satisfied for each layer.We extend Eq. ( 1) in order to describe a vegetation layer of volume V i , area A i , and thickness h i : All terms are defined in Table 1.The specific heat of each vegetation layer (θ i ) is assumed equal to that of water, and is modulated according to the leaf area density (m 2 m −3 ) at that level.Since the fluxes in the model are described per square metre, A i may be represented by the plant area density (PAD; m 2 m −3 ) for that layer, where plant denotes leaves, stems, grasses or any other vegetation included in optical LAI measurements.Note that LAI, which has units of m 2 m −2 , is a value that describes the integration over the whole of the canopy profile of PAD (which is applied per metre of height, hence the dimension m 2 m −3 ).Canopy layers that do not contain foliage may be accounted for at a level by assigning that R i = R i = ∞ for that level (i.e. an open circuit).
Rewriting Eq. ( 14) in terms of the state variables and resistances that are shown in Fig. 1 means that R i is the resistance to sensible heat flux and R i the resistance to latent heat flux.Dividing both sides of the equation by V i , the volume of the vegetation layer (equal to h i multiplied by A i ), expresses the sensible and latent heat fluxes between the leaf and the atmosphere as Note that this is the first of three key equations that are labelled (a), (b) or (c) on the left hand side, throughout.

Vertical transport within a column
The transport equation between each of the vegetation layer segments may be described as where div is the operator that calculates the divergence of the vector field, χ is the property under question, ρ is the fluid density, u is the horizontal wind speed vector, S χ is the concentration for the property in question, and is a parameter that will in this case be the diffusion coefficient k(z).
To derive from this expression the conservation of scalars equation, as might be applied to vertical air columns, we proceed according to the finite volume method, as used in the FRAME (Fine Resolution Atmospheric Multi-pollutant Exchange; Singles et al., 1998) model and as outlined in Vieno (2006) and derived from Press (1992).The final equation is specific to a one-dimensional model, and therefore does not include a term of the influence of horizontal wind.The resulting expression is sufficiently flexible to allow for variation in the height of each layer, but we preserve vegetation layers of equal height here for simplicity: where F is the vertical flux density, and z represents coordinates in the vertical and x coordinates in the stream-wise direction.χ may represent the concentration of any constituent that may include water vapour or heat, but also gas or aerosol phase concentration of particular species.S represents the source density of that constituent (in this case the fluxes of latent and sensible heat from the vegetation layer), and the transport k(z) term represents the vertical transport between each layer.In the equation above, we substitute the flux-gradient relationship according to the expression: This approach allows future applications to include a supplementary term to simulate emissions or deposition of gas or aerosol-based species using the same technique.
The transport terms, per level i in the vertically discretised form, 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.More details are outlined in that paper, but the in-canopy wind speed is dependent on C Deff , the effective phyto-element canopy drag coefficient.This is defined according to Wohlfahrt and Cernusca (2002): where LAD is the leaf area density and a 1 , a 3 , and a 5 are parameters to be defined.This second-order closure model also provides profiles of σ w , the standard deviation in vertical velocity and T L , the Lagrangian timescale within the canopy.The term T L is defined as in the model of Raupach (1989a) and represents the time, since emission, at which a flux transitions from the near field (emitted equally in all directions, and not subject to eddy diffusivity), and the far field (which is subject to normal eddy diffusivity and gradient influences).The eddy diffusivity k i (z) is then derived in the far field using the expressions from Raupach (1989b): However, the simulation of near-field transport requires ideally a Lagrangian solution (Raupach, 1989a).As that is not directly possible in this implicit solution, we instead adopt a method developed by Makar et al. (1999) (andlater Stroud et al., 2005 andWolfe andThornton, 2011) for the transport of chemistry species in canopies for which a nearfield correction factor R nf is introduced to the far-field solution, which is based on the ratio between the Lagrangian timescale T L and τ , which represents the time since emission for a theoretical near-field diffusing cloud of a canopy source as defined in Raupach (1989a), which (unlike for the far-field) acts as point source travelling uniformly in all directions.R nf appears to depend on canopy structure and on venting (Stroud et al., 2005), but has yet to be adequately described.
To obtain a satisfactory match for the profiles here, we apply a weighting factor W nf , based on the friction velocity u * (Chen et al., 2016).Resistance analogue for a multi-layer canopy approximation of n levels, to which the energy balance applies at each level.Refer to Table 1 for interpretation of the symbols.
There is thus a modified expression for k i , with R nf acting effectively as a tuning coefficient for the near-field transport: The necessity to account for the near-field transport effect in canopies, remains a question under discussion (Mc-Naughton and van den Hurk, 1995;Wolfe and Thornton, 2011).

Fluxes of sensible and latent heat between the canopy layers
We re-write the expression for scalar conservation (Eq.16, above), as applied to canopies, as a pair of expressions for the fluxes of sensible and latent heat; therefore, comparing with Eq. 16, χ ≡ T or q, F ≡ H or λE, and S ≡ (the source sensible or latent heat flux at each vegetation layer).Neither the sensible or latent heat-flux profile is constant over the height of the canopy.The rate of change of T a,i (the temperature of the atmosphere surrounding the leaf at level i) and q a,i (the specific humidity of the atmosphere surrounding the leaf at level i) are proportional to the rate of change of the respective fluxes with height and the source of heat fluxes from the leaf at that level: We now assume the flux-gradient relation and therefore write Eq. ( 19) according to sensible heat flux at level i: which is substituted in Eq. ( 23) and following the same approach for the expression for latent heat flux at level i, λE a,i : which is, again, substituted in Eq. ( 23): We have now defined the three key equations in the model.
-Eq.(a) balances the energy budget at each canopy air level.
-Eq.(b) balances heat fluxes vertically between each vegetation level and horizontally between each vegetation level and the surrounding air.
-Eq.(c) balances humidity fluxes in the same sense as for Eq.(b).
The equations must be solved simultaneously, whilst at the same time satisfying the constraints of an implicit scheme.

Write equations in implicit format
The difference between explicit and implicit schemes is that an explicit scheme will calculate each value of the variable (i.e.temperature and humidity) at the next time step entirely in terms of values from the present time step.An implicit scheme requires the solution of equations that couple together values at the next time step.The basic differencing scheme for implicit equations is described by Richtmyer and Morton (1967).In that work, they introduce the method with an example equation where B denotes a linear finite difference operator, t, x, y are increments in the respective co-ordinates and u t , u t+1 are the solutions at respectively steps t and t + 1.
It is therefore assumed that B depends on the size of the increments t, x, y and that, once known, it may be used to derive u t+1 from u t .So if B can be determined we can use this relationship to calculate the next value in the temporal sequence.However, we necessarily need to know the initial value in the sequence (i.e.u 0 ).This means that it is an initial value problem.Now, the equivalent of Eq. ( 30), in the context of a column model, such as LMDz, takes the form: This describes the state variable X (for example temperature) at level i, in relation to the value at level i − 1. C X i and D X i are coupling coefficients that are derived in that scheme.In this particular example, the value of W i at time t is defined in terms of X i−1 at the same time step.
To maintain the implicit coupling between the atmospheric model (i.e.LMDz) and the land surface model (i.e.OR-CHIDEE), we need to express the relationships that are outlined above in terms of a linear relationship between the present time step t and the next time step t + 1.We therefore re-write equations (a), (b), and (c) in implicit form (i.e. in terms of the next time step, which is t + 1), as explained in the following subsections.

Implicit form of the energy balance equation
We substitute the expressions for leaf level vapour pressure Eq. ( 4) to the energy balance equation ( 15), which we rewrite in implicit form:

Implicit form of the sensible heat-flux equation
We differentiate Eq. ( 25) according to the finite volume method Eq. ( 17), and divide by V i :

Implicit form of the latent heat-flux equation
We differentiate Eq. ( 29) according to the finite volume method Eq. ( 17), and divide by V i :

Solution by induction
These equations are solved by deducing a solution based on the form of the variables in Eqs. ( 32), (34), and (36) above.The coefficients within this solution can then be determined,with respect to the boundary conditions, by substitution.This is solution by induction.With respect to Eq. ( 34), we wish to express T t+1 a,i in terms of values further down the column, to allow the equation to be solved by moving up the column, as in Richtmyer and Morton (1967).There is also an alternative method to solve these equations also derived from that text, which we describe in the Supplement.In order to solve by implicit means, we make the assumption (later to be proved by induction) that We then also re-write these expressions in terms of the values of the next level: (ii) q t+1 a,i+1 =A q,i+1 q t+1 a,i + B q,i+1 where A T ,i , B T ,i , C T ,i , D T ,i , A q,i , B q,i , C q,i , and D q,i are constants for that particular level and time step and are, as yet, unknown, but will be derived.We thus substitute Eqs. ( 38) and (40) into Eq.( 34) to eliminate T t+1 .Symmetrically, we substitute Eqs. ( 39) and (41) into Eq.( 36) to eliminate q t+1 .For the vegetation layer, we conduct a similar procedure, in which the leaf level temperature is described as follows (where E i , F i , and G i are known assumed constants for the level and time step in question): Now the coefficients A T ,i , B T ,i , C T ,i , D T ,i , A q,i , B q,i , C q,i , and D q,i can be described in terms of the coefficients from the level above and the potentials (i.e.T and q) at the previous time step, which we can in turn determine by means of the boundary conditions.So we have a set of coefficients that may be determined for each time step, and we have the means to determine T S (and q S via the saturation assumption).We thus have a process to calculate the temperature and humidity profiles for each time step by systematically calculating each of the coefficients from the top of the column (the downwards sweep) then calculating the initial value (the surface temperature and humidity) and finally calculating each T a , q a , and T leaf by working up the column (the upwards sweep).The term T t+1 leaf,i can also be described in terms of the variables at the level below by T t+1 leaf,i+1 using Eq.(iii) and its terms E i , F i , and G i .
Table 2. Input coefficients at the top layer of the model, where A T ,n , B T ,n . . ., etc. are the respective coefficients at the top of the surface model and A T ,atmos , B T ,atmos are the coefficients at the lowest level of the atmospheric model.

Stand-alone model Coupled model
A q,n = 0 A q,n = A q,atmos B q,n = B q,input B q,n = B q,atmos C q,n = 0 C q,n = 0 D q,n = 0 D q,n = 0

The upper boundary conditions
In stand-alone simulations, the top level variables A T ,n , C T ,n , D T ,n and A q,n , C q,n , D q,n are set to zero, and B T ,n and B q,n are set to the input temperature and specific humidity, respectively, for the relevant time step (as in Best et al., 2004).In coupled simulations, A T ,n , B T ,n and B q,n , C q,n are taken from the respective values at the lowest level of the atmospheric model.Table 2 summarises the boundary conditions for both the coupled and un-coupled simulations.

The lower boundary condition
We need to solve the lowest level transport equations separately, using an approach that accounts for the additional effects of radiation emitted, absorbed, and reflected from the vegetation layers: where η 1,S , η 2,S , and η 3,S are components of the radiation scheme, and ξ 1 , ξ 2 , ξ 3 , and ξ 4 are components of the surface flux (where φ H = ξ 1 + ξ 2 T t+1 S and φ LE = ξ 3 + ξ 4 T t+1 S ; refer to Sect.S3.2 of the Supplement).
The interaction with the soil temperature is by means of the soil flux term J soil .Beneath the soil surface layer, there is a seven layer thermal soil model simulating energy transport and storage in the soil (Hourdin, 1992), which is unchanged from the standard version of ORCHIDEE.

Radiation scheme
The radiation approach is the application of the long-wave radiation transfer matrix (LRTM) (Gu, 1988;Gu et al., 1999), as applied in Ogée et al. (2003).This approach separates the calculation of the radiation distribution completely from the Geosci.Model Dev., 9, 223-245, 2016 www.geosci-model-dev.net/9/223/2016/implicit expression.Instead, a single source term for longwave radiation is added at each level.This means that the distribution of radiation is no longer completely explicit (i.e. an explicit scheme makes use of information only from the present and not the next time step).However, an advantage of the approach is that it accounts for a higher order of reflections from adjacent levels than the single order that is assumed in the process above.
The components for long-wave radiation are abbreviated as The shortwave radiation component is abbreviated as where η 1,i , η 2,i , and η 3,i are components of the radiation scheme.η 1,i accounts for the components relating to emission and absorption of LW radiation from the vegetation at level i (i.e. the implicit parts of the long-wave scheme relating to the level i) and η 2,i the components relating to radiation from vegetation at all other levels incident on the vegetation at level i (i.e. the non-implicit part of the long-wave scheme), as well as explicit parts from the level i. η 3,i is the component of the SW radiation scheme -it describes the fraction of the total downwelling shortwave light that is absorbed at each layer, including over multiple forward-and back-reflections, as simulated by the multilayer albedo scheme (McGrath et al., 2016).The fraction of original downwelling SW radiation that is ultimately reflected from the surface and from the vegetation cover back to the canopy can then be calculated using this information.

Selected site and observations
Given the desired capability of the multi-layer model to simulate complex within-canopy interactions, we selected a test site with an open canopy.This is because open canopies may be expected to be more complex in terms of their interactions with the overlying atmosphere.In addition, long-term data measurements of the atmospheric fluxes had to be available in order to validate the performance of the model across years and seasons, and within-canopy measurements were required in order to validate the capacity of the model to simulate within-canopy fluxes.One site that fulfilled these requirements was the long-term measurement site at Tumbarumba in south-eastern inland Australia (35.6 • S, 148.2 • E, elevation ∼ 1200 m) which is included in the global Fluxnet network (Baldocchi et al., 2001).The measurement site is a Eucalyptus delegatensis canopy, a tall, temperate evergreen species (∼ 40 m).With an LAI of ∼ 2.4, the canopy is described as moderately open (Ozflux, 2013).

Forcing and model comparison data
As a test of stability over a long-term run, the model was forced (i.e.run offline, independently from the atmospheric model) using above-canopy measurements.The forcing data that were used in this simulation was derived from the long-term Fluxnet measurements for the years 2002 to 2007, specifically above-canopy measurements of long-wave and shortwave radiation, temperature, humidity, wind speed, rainfall, and snowfall.The first 4 years of data, from 2002 to 2005, were used as a spin-up to charge the soil to its typical water content for the main simulation.The biomass from the spin-up was overwritten by the observed leaf biomass to impose the observed LAI profile.Soil carbon is not required in this study, which justify the short spin-up time.The years 2006 and 2007 were then used as the main part of the run.Although the shortwave radiation was recorded at the field site in upwelling and downwelling components (using a set of directional radiometers), the long-wave radiation was not.As a consequence, the outgoing long wave was calculated using the recorded above-canopy temperature by the Stefan-Boltzmann law with an emissivity factor of 0.96 (a standard technique for estimating this variable; e.g. Park et al., 2008).This value is then subtracted from the net radiation, together with the two shortwave components, to obtain an estimation of the downwelling long-wave radiation with which to force the model.
For the validation of the within-canopy processes, more detailed measurement data were required.For the same site there exists data from an intensive campaign of measurements made during November 2006 (Austral summer), described by Haverd et al. (2009).Within the canopy, profiles of temperature and potential temperature were recorded over the 30-day period and, for a number of days (7-14 November), sonic anemometers were used to measure wind speed and sensible heat flux in the vertical profile at eight heights as well.Measurements were also made over the 30-day period of the soil heat flux and the soil water content.These within-canopy data were used for validation of the modelled output, but the same above-canopy long-term data (i.e. the Fluxnet data) were used in the forcing file in all cases.No further measurements were collected specifically for this publication.The measurement data (i.e. the data both from the 1-month intensive campaign and the long-term Fluxnet measurements at the same site; Ozflux, 2013) were prepared as an ORCHIDEE forcing file, according to the criteria for gapfilling missing data (Vuichard and Papale, 2015).

Model set-up
The multi-layer module that is described in this paper only calculates the energy budget.Its code was therefore integrated in the enhanced model ORCHIDEE-CAN (revision 2966), and relies on that larger model for input-output operations of drivers and simulations, as well as the calcula- tion of soil hydrology, soil heat fluxes, and photosynthesis (see Table 3 for other input).A more detailed description of how these processes are implemented in ORCHIDEE-CAN is provided in Naudts et al. (2015).
The ORCHIDEE-CAN model is capable of simulating the canopy vegetation structure prognostically, and these prognostic vegetation stands have now been linked to the multilayer energy budget calculations.In these tests, a vegetation profile was prescribed, in order to obtain a simulation as close as possible to the observed conditions.That is to say, the stand height to canopy radius ratios of the trees across several size classes in ORCHIDEE-CAN were prescribed over the course of the spin-up phase to an approximation of the Tumbarumba LAD profile.The assigned height to radius profiles are provided in Table 3. LAD is an estimate of the sum of the surface area of all leaves growing on a given land area (e.g. per m 2 ) over a metre of height.It is effectively LAI (m 2 per m 2 ) per canopy levels, and thus has units of m 2 per level of the canopy.As there were no LAD profiles available for the field site at the time of measurement, data from Lovell et al. (2012) for the Tumbatower profile, as depicted in Fig. 3 of that publication, were used as a template.The profile was scaled according to the measured site LAI of 2.4, resulting in the profile shown in Fig. 2. As no gap-forming or stand replacement disturbances have been recorded at the site, the vertical distribution of foliage was assumed unchanged over the period between the different measurement campaigns.
Several tuning coefficients were applied to constrain the model, which are listed in Table 3.A combination of manual and automated tuning was used to tune the model as closely as possible to the measurement data.The key tuning coefficients were R b, fac , a tuning coefficient for the leaf boundarylayer resistance, R s, fac , a coefficient for the stomatal resistance, and R nf , the near-field correction factor to the modified eddy diffusivity coefficient K * i , the coefficients a 1 , a 3 , and a 5 , corresponding to the definition C deff from Eq. ( 20), and , a correction factor for the total LAI to allow for canopy gaps.A more detailed guide to the model tuning is provided in Chen et al. (2016).

Results
Although the aim of this study is to check the performance of our multi-layer energy budget model against site-level observations, it should be noted that site-level energy fluxes come with their own limitations that result in a so-called closure gap.The closure gap is reflected in a mismatch between the net radiation and the fluxes of latent, sensible and soil heat.For the observations used in this study, the closure gap was ∼ 37 W m −2 (7.5 % of total fluxes) during the day and 4 W m −2 (4.6 %) during the night.Underestimation of the data and mismatches exceeding the closure gap likely indicate a shortcoming in the model.At a fundamental level, energy budget models distribute the net radiation between sensible, latent, and soil heat fluxes.Evaluation of these component fluxes becomes only meaningful when the model reproduces the net radiation (Fig. 3).Note that through its dependency on leaf temperature the calculation of the long-wave component of net radiation depends on the sensible, latent, and soil heat fluxes.Taken as a whole, there is a very good correlation between the observation-driven and model-driven net long-wave radiation (r 2 = 0.96).However, when the data are separated into night-time and daytime, as shown, a clear cycle is revealed, for which the model overestimates daytime radiation and underestimates radiation at night.This discrepancy is likely a result of actual daytime heat storage in the soil being underestimated in the model.A portion of the upwelling long-wave radiation is sourced from temperature changes in fluxes from the soil model, and the rest from vegetation.So if the daytime surface layer temperature is underestimated by the model, we expect reduced net long-wave predicted radiation, compared to that which is measured, and vice versa for the night-time scenario.The use of above-canopy air temperature, instead of radiative temperature (which was not measured) may also contribute to inaccuracies in the predicted long-wave radiation.In terms of the current parametrisation, and for the site under study, the annual cycles for both sensible and latent heat are well simulated (Fig. 4a and c).In addition, no clear systematic bias was observed between summer and winter (Fig. 4b and d).But, as shown, there is an overall systematic bias of +12.7 W m −2 for sensible heat and +10.7 W m −2 for latent heat flux, when averaged over the whole year.Such a bias represents ∼ 23 % of sensible heat and ∼ 15 % of latent heat fluxes.
The analysis proceeded by further increasing the temporal resolution and testing the capacity of the model to reproduce diurnal flux cycles.The model overestimates the diurnal peak in sensible heat flux, whilst the latent heat flux is underestimated by a slightly larger magnitude (Fig. 5b).The diurnal pattern of the model biases persists in all four seasons (Fig. S1a-d in the Supplement).We see that the maximum mean discrepancy between measured and modelled sensible heat-flux ranges from +95 to −84 W m −2 (Fig. 5b) and latent heat flux from −49 to 43 W m −2 (Fig. 5d).Over the course of the year, the difference is largest in the autumn and smallest in the summer (Fig. S2a-d).However, from the net radiation (i.e. the sum of downwelling minus upwelling for long wave and shortwave) we can see that there is a discrepancy between measured and modelled that acts to offset in part the discrepancy observed in the flux plots (Fig. 5a-f).
Long-term measurements from above the forest and data from a short intensive field campaign were jointly used to evaluate model performance at different levels within the canopy.As was the case for the annual cycle, the sinusoidal cycles resulting from the diurnal pattern in solar angle are well matched (Fig. 6a-d).Sensible heat flux was measured below and above the canopy and the model was able to simulate the trend of the gradient, though inexactly (Fig. 6a,  b).Latent heat flux at an equivalent height of 2 m was not recorded (Fig. 6c).However, the match in magnitude of the measured data is not accurately simulated hour by hour (Fig. 6e).
Using the current parameters, there is a discrepancy between the measured and the modelled temperature gradients within the canopy (Fig. 7).It should be noted that the mean values are strongly determined by a few extreme hours (not shown).As such the model is capable of simulating the majority of the time steps but fails to reproduce the more extreme observations (not shown).During the daytime, the strong negative gradient in the measured output is only partly reflected in the modelled slopes.At night-time, there is a clear positive gradient for the measured data, which is matched by the model.These profiles demonstrate that the trend of in-canopy gradients can be replicated by the parametrisation of the model.
The version of the model used in these tests so far is composed of 30 levels, with 10 levels in the understorey, 10 in the canopy vegetation profile, and 10 between the top of the canopy and the lowest layer of the boundary-layer model, in order to provide a high-resolution simulation and a test of the stability of the scheme.However, a canopy simulation of such detail might be overly complex for a canopy model that is to be coupled to an atmospheric simulation, in terms of additional run time required.To provide an evaluation of the difference in fluxes that were predicted by a model of lower resolution, the same tests were conducted with the model composed of 30, 15, 8, and 5 levels, overall, within which there are an equal number of empty levels below and above sections containing 10, 5, 2, and a single vegetation layer, respectively (note that in all cases the vegetation levels are simulated separately from the surface soil level, treated separately in each case, and represent a separate layer) (cf.Dolman, 1993).
Tests were conducted for both hourly mean (Fig. 8) and daily mean (Fig. 9), both calculated over the course of a year, and for a moving average.These plots show the root mean square (RMS) error between the original set-up and the modified number of levels.Figure 8 shows that there is already a significant difference between the calculated hourly sensible heat flux for the version of the model with 10 canopy layers (30 total profile layers) and 5 canopy layers (15 total profile layers), that reaches a peak of 28 W m −2 , but that the discrepancy is substantially larger for the 2 canopy layer (8 total levels) and a single canopy layer (5 total levels) cases.In the case of the latent heat flux (Fig. 8b), the discrepancy is most marked for the single canopy layer case, with a peak difference of 60 W m −2 .Considering the daily averages, for sensible heat flux the difference between the different model set-ups is always below 25 W m −2 (Fig. 9a) in all cases.For latent heat flux, there is more considerable divergence, up to 42 W m −2 , for the single canopy later set-up (Fig. 9b).

Discussion
The proposed model is able to simulate fluxes of sensible and latent heat above the canopy over a long-term period, as has been shown by simulation of conditions at a Fluxnet site on a long-term, annual scale (Figs. 4 and 5), and over a concentrated, week-long period (Fig. 6).Although these figures show a discrepancy between measured and modelled fluxes, we see from Fig. 5 that the modelled overestimate of sensible heat flux is offset by an underestimation of latent heat flux and of net radiation.In the study of land-atmosphere interactions, the multi-layer model functions to a standard comparable to single-layer models, but calculates in more detail in-canopy transport, the sources and sinks of heat, and profiles of temperature and specific humidity profiles.An earlier iterative model applied to the same site (Haverd et al., 2009) found differences of the order of 50 W m −2 at maximum for the mean daily average latent and sensible above canopy heat fluxes.The innovation of this model 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.Uniquely for a canopy model, this is achieved without iterations, as the mathematics has been derived to use the same implicit coupling technique as the existing surface-atmosphere coupling applied in ORCHIDEE-LMDz (Polcher et al., 1998;Best et al., 2004), but now over the height of the canopy.This also means that the model is scalable without impacting heavily on runtimes.For largescale applications, performance within the canopy must be further constrained through comparison with intensive incanopy field campaigns from diverse ecosystems.

Simulation of aerodynamic resistance
In this study, the aerodynamic coefficient that is used in single-layer models was replaced by an eddy diffusivity pro- file, the purpose of which is twofold: (i) to develop a transport coefficient that is based on the vertical canopy profile, and (ii) to more accurately represent the in-canopy gradients of temperature and specific humidity.In this way, it was hoped to contribute to a model that can better allow for such features as vertical canopy gaps (i.e.trunk space between a well separated under and overstorey), horizontal gaps, transport and chemistry between different sections of the canopy, tree growth and the mix of different kinds of vegetation in the same surface layer simulation (e.g.Dolman, 1993).To be able to do this, a height-based transport closure model was used to simulate within-canopy transport.
The transport closure model used here can be compared to the previous single-layer approach within ORCHIDEE.In that approach, aerodynamic interaction between the land surface and the atmosphere is parametrised by the atmospheric resistance R a and the architectural resistance R 0 .R a is typically calculated through consideration of the roughness height of the canopy (i.e.small for flat surfaces, large for uneven tall surfaces), which in turn is parametrised in surface layer models by canopy height (e.g.LSCE/IPSL, 2012); however, LAI can display a better correlation with roughness length (a critical parameter) than it does to canopy height (Beringer et al., 2005).In parametrising the roughness length in terms of canopy height alone, no account is made for the clumping of trees, the density of the forest, or the phenological changes in stand profile (other than the height) as the stand grows.Some of these changes are compensated for in R 0 , the structural coefficient that is unique to each PFT grouping, but does not allow for more subtle effects.To be able to satisfactorily explore such results in a modelling study requires an accurate parametrisation of within-canopy transport.
In this study, canopy transport is parametrised by Ktheory, applying the closure model of Massman and Weil  8 for a set-up with 10 canopy layers (out of 30 total profile levels), 5 canopy layers (of 15 total profile levels), 2 canopy levels (of 8 total profile levels), and 1 canopy levels (of 5 total profile levels), for a run over the course of a year, expressed as a daily average (a) 20-day moving average of the root mean square difference of the daily mean of sensible heat flux; (b) as (a) for latent heat flux (1999) to derive the in-canopy turbulence statistics, based both on the LAI profile and the canopy height.The simulation produces a good estimation of above-canopy fluxes, but the differences between daytime and night-time profiles are not well described using the original parametrisation (Fig. 7).This means that the model overestimates the nighttime canopy transport, as compared to the daytime simulation.Taking a broader view, studies of chemical species transport have demonstrated that K-theory, sometimes constrained by a scaling factor, remains a reasonable approximation for above-canopy fluxes, even if the within-canopy gradients are not entirely correct (Gao et al., 1989;Dolman and Wallace, 1991;Makar et al., 1999;Wolfe and Thornton, 2011).The justification for such a scaling factor seems to vary in terms of the form of the canopy structure, likely related to canopy openness (McNaughton and van den Hurk, 1995;Stroud et al., 2005).Here, too, we find that a scaling factor is necessary to match the gradient fluxes though the scaling factor required varies according to the time of day.We now also parametrise the effective phyto-element canopy drag coefficient, C Deff , in order to obtain a more accurate simulation.For a completely satisfactory resolution of this issue, it will be necessary to derive a method to reformulate the method of Raupach (1989a, b) in an implicit form, which lies outside the scope of this paper.

Simulation of energy partition throughout canopy and soil surface
Trees in a spruce forest have been reported to account for 50-60 % of the latent heat flux; moisture in the soil itself would have a reduced impact due to soil shading (Baldocchi et al., 2000).Another study found that the fraction of radiation that reaches the soil ranges from 0.05 (forest) to 0.12 (tundra) (Beringer et al., 2005).The same study found that the latent heat flux correlates most closely with the leaf-level vapour pressure deficit -that is to say the difference between the leaf level saturation vapour pressure and the actual vapour pressure of the outside air, rather than between air water vapour pressure and the saturation vapour pressure at the soil level.Since a single-layer canopy model regards both the canopy and soil surface as the same entity, the aforementioned subtleties will inevitably be lost in the modelling.Although, the partition of energy between soil surface and vegetation is site dependent -a well-hydrated site would behave differently to one in an arid region -it is effects such as these that a more realistic energy budget scheme should be able to simulate.Being able to simulate separately the vegetation allows for the partitioning of fluxes between the vegetation and the soil.For example, from the measurements (Fig. 6a and b), we see that approximately 50 % of the sensible heat that is measured above the canopy is sourced not from the soil surface, but from the overlying vegetation, as this is the difference between the measured flux at 1 m and that above the canopy.The model simulated approximately 70 % of the sensible heat being sourced from the soil.Hence, the model overestimates the measured contribution from the soil.There is no equivalent measurement at 1m for the latent heat flux, but the model calculates that approximately 50 % of the latent heat flux is sourced from the vegetation rather than the soil surface.
This model also simulates leaf temperature that may be verified by leaf level measurements, where such measurements exist (Helliker and Richter, 2008).Such a comparison would require additional developments (as is discussed in the following section) because leaf temperature measurements strongly depend on the approach that is used.

Outlook
This document lays out the framework for the model design, but it allows for the further implementation of many features in site-level to global-scale scenarios: -As the method calculates leaf temperature and incanopy radiation, it will be possible to simulate the explicit emission by leaves of certain common biogenic volatile organic compounds (BVOCs), such as isoprene and monoterpene (Guenther et al., 1995(Guenther et al., , 2006)).As the method calculates in-canopy gradients of temperature, specific humidity, and radiation, it is possible to simulate more accurately chemical reactions that depend on these factors such as the NO x and O 3 cycle within and above canopies (Walton et al., 1997) and the formation and size distribution of aerosol interactions (Atkinson and Arey, 2003;Nemitz et al., 2004a, b;Ehn et al., 2014), which may act as cloud condensation nuclei and thus again feedback into radiation absorption interactions at the atmospheric component of a coupled model such as LMDz-ORCHIDEE.
-Separate computation of vegetation and soil temperatures, which can be very different, and then estimate the directional whole canopy temperature as measured by remote sensing sensors.It may then be possible to assimilate this variable from remote sensing products in order to better constrain the energy budget.
-Recent research in ecology demonstrates further the need to better understand canopy microclimates, and in particular gradients of state variables such as temperature and specific humidity, and radiation penetration.For example, temperature gradients in the rainforest exert a key influence on the habitat choices of frogs, and changes to such a microclimate threaten their survival (Scheffers et al., 2013).In a similar vein, microclimate affects in canopies can act as a buffer to changes in the climate overall (i.e. the macro-climate) in terms of the survival of species in the sub-canopy (Defraeye et al., 2014).Therefore structural forest changes, such as forest thinning, will reduce buffer lag effect, but it is only with well-designed canopy models that an informed prediction of the long-term consequences of land management policies can be made.

Conclusions
A new numerical model for ORCHIDEE-CAN (revision 2966) has been developed that enables the simulation of vertical canopy profiles of temperature and moisture using a non-iterative implicit scheme.This means that the new model may also be used when coupled to an atmospheric model, without compromising computer run-time.
Initial tests demonstrated that the model runs stably, balances the energy budget at all levels, and provides a good simulation of the measured field data, both on short timescales of a few days, and over the course of a year.As demonstrated, the model structure allows for coupling/linking to a more physical-based albedo scheme (Pinty et al., 2006;Mc-Grath et al., 2016;Naudts et al., 2015) and the implementation of a vertically discretised stomatal conductance scheme.
Reducing the vertical discretisation of the canopy from 10 layers to 5, 2, and 1 layer increased the root mean square error (RMSE) between the model and the observations for λE and H and thus demonstrates the overall benefits of introducing a multi-layer energy budget scheme.The multi-layer energy budget model component outlined here may be used to simulate canopies in more detail and variety.It also offers the potential to integrate with other parts of ORCHIDEE for enhanced simulation of CO 2 transport, emission of volatile organic compounds (VOCs) and leaf-scale plant hydraulics.
The Supplement related to this article is available online at doi:10.5194/gmd-9-223-2016-supplement.
Figure 1.Resistance analogue for a multi-layer canopy approximation of n levels, to which the energy balance applies at each level.Refer to Table1for interpretation of the symbols.

Figure 2 .
Figure 2. Simulated leaf area density profile corresponding to the forest canopy at the Tumbarumba site used in this study.

MFigure 3 .
Figure 3. Correlation of observed upwelling long-wave radiation (derived from the measured above-canopy temperature) and the upwelling long-wave radiation that is simulated by ORCHIDEE-CAN.Night-time data (corresponding to a downwelling shortwave radiation of < 10 Wm −2 ) are plotted in black, and daytime data are plotted in orange.

Figure 4 .
Figure 4. Daily mean for measured (circles) and modelled (triangles) over a year-long run for (a) sensible heat flux; (b) difference between measured and modelled sensible heat flux; (c) latent heat flux; and (d) difference between measured and modelled latent heat flux; 1 in every 5 data points is shown, for clarity.Thick lines show the respective 20-day moving average for each data set.Graphs (b) and (d) also show the overall mean of individual data points.

Figure 5 .
Figure 5. Hourly means for measured (circles) and modelled (triangles) for (a) sensible heat flux; (b) difference between measured and modelled sensible heat flux, as calculated over the period of 2006; 1 in every 10 days is plotted for clarity; (c) and (d): as above for latent heat flux; (e) and (f): as above for net radiation.Continuous lines show the overall mean.

Figure 6 .
Figure 6.Short-term simulated and observed energy fluxes: (a) sensible heat fluxes at a height of 50 m; (b) as (a) for latent heat flux; (c) measured and modelled sensible heat flux at 2 m above the ground; (d) modelled latent heat flux at 2 m above the ground (measurements not available); (e) difference in measured and modelled sensible and latent heat flux at a height of 50 m.

Figure 8 .
Figure8.Comparison of model performance for a set-up with 10 canopy layers (out of 30 total profile levels), 5 canopy layers (of 15 total profile levels), 2 canopy levels (of 8 total profile levels), and 1 canopy levels (of 5 total profile levels), expressed as an hourly average for a year-long run.(a) Root mean square difference between the different model runs for sensible heat flux; (b) as (a) for latent heat flux

Figure 9 .
Figure9.As for Fig.8for a set-up with 10 canopy layers (out of 30 total profile levels), 5 canopy layers (of 15 total profile levels), 2 canopy levels (of 8 total profile levels), and 1 canopy levels (of 5 total profile levels), for a run over the course of a year, expressed as a daily average (a) 20-day moving average of the root mean square difference of the daily mean of sensible heat flux; (b) as (a) for latent heat flux

Table 3 .
Parameters and tuning coefficients used in the model for simulation described in this work.