The interactions between soil-biosphere-atmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8 - Part 1: Model description

Abstract

Land surface models (LSMs) are pushing towards improved realism owing to an increasing number of observations at the local scale, constantly improving satellite data sets and the associated methodologies to best exploit such data, improved computing resources, and in response to the user community. As a part of the trend in LSM development, there have been ongoing efforts to improve the representation of the land surface processes in the interactions between the soil-biosphere-atmosphere (ISBA) LSM within the EXternalized SURFace (SURFEX) model platform. The force-restore approach in ISBA has been replaced in recent years by multi-layer explicit physically based options for sub-surface heat transfer, soil hydrological processes, and the composite snowpack. The representation of vegetation processes in SURFEX has also become much more sophisticated in recent years, including photosynthesis and respiration and biochemical processes. It became clear that the conceptual limits of the composite soil-vegetation scheme within ISBA had been reached and there was a need to explicitly separate the canopy vegetation from the soil surface. In response to this issue, a collaboration began in 2008 between the high-resolution limited area model (HIRLAM) consortium and Météo-France with the intention to develop an explicit representation of the vegetation in ISBA under the SURFEX platform. A new parameterization has been developed called the ISBA multi-energy balance (MEB) in order to address these issues. ISBA-MEB consists in a fully implicit numerical coupling between a multi-layer physically based snowpack model, a variable-layer soil scheme, an explicit litter layer, a bulk vegetation scheme, and the atmosphere. It also includes a feature that permits a coupling transition of the snowpack from the canopy air to the free atmosphere. It shares many of the routines and physics parameterizations with the standard version of ISBA. This paper is the first of two parts; in part one, the ISBA-MEB model equations, numerical schemes, and theoretical background are presented. In part two (Napoly et al., 2016), which is a separate companion paper, a local scale evaluation of the new scheme is presented along with a detailed description of the new forest litter scheme.

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/ doi:10.5194/gmd-10-843-2017 Author(s) 2017. CC Attribution 3.0 License.

Aaron Boone1, Patrick Samuelsson2, Stefan Gollvik2, Adrien Napoly1, Lionel Jarlan3, Eric Brun1, and

Bertrand Decharme1

1CNRM UMR 3589, Mto-France/CNRS, Toulouse, France

2Rossby Centre, SMHI, 601 76 Norrkping, Sweden

3CESBIO UMR 5126 UPS, CNRS, CNES, IRD, Toulouse, France

Correspondence to: Aaron Boone ([email protected])

Received: 14 October 2016 Discussion started: 27 October 2016

Revised: 23 January 2017 Accepted: 27 January 2017 Published: 21 February 2017

Abstract. Land surface models (LSMs) are pushing towards improved realism owing to an increasing number of observations at the local scale, constantly improving satellite data sets and the associated methodologies to best exploit such data, improved computing resources, and in response to the user community. As a part of the trend in LSM development, there have been ongoing efforts to improve the representation of the land surface processes in the interactions between the soilbiosphereatmosphere (ISBA) LSM within the EXternalized SURFace (SURFEX) model platform. The forcerestore approach in ISBA has been replaced in recent years by multi-layer explicit physically based options for sub-surface heat transfer, soil hydrological processes, and the composite snowpack. The representation of vegetation processes in SURFEX has also become much more sophisticated in recent years, including photosynthesis and respiration and biochemical processes. It became clear that the conceptual limits of the composite soilvegetation scheme within ISBA had been reached and there was a need to explicitly separate the canopy vegetation from the soil surface. In response to this issue, a collaboration began in 2008 between the high-resolution limited area model (HIRLAM) consortium and Mto-France with the intention to develop an explicit representation of the vegetation in ISBA under the SURFEX platform. A new parameterization has been developed called the ISBA multi-energy balance (MEB) in order to address these issues. ISBA-MEB consists in a fully implicit numerical coupling between a multi-layer physically based snow-

The interactions between soilbiosphereatmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8 Part 1: Model description

pack model, a variable-layer soil scheme, an explicit litter layer, a bulk vegetation scheme, and the atmosphere. It also includes a feature that permits a coupling transition of the snowpack from the canopy air to the free atmosphere. It shares many of the routines and physics parameterizations with the standard version of ISBA. This paper is the rst of two parts; in part one, the ISBA-MEB model equations, numerical schemes, and theoretical background are presented. In part two (Napoly et al., 2016), which is a separate companion paper, a local scale evaluation of the new scheme is presented along with a detailed description of the new forest litter scheme.

1 Introduction

Land surface models (LSMs) are based upon fundamental mathematical laws and physics applied within a theoretical framework. Certain processes are modeled explicitly while others use more conceptual approaches. They are designed to work across a large range of spatial scales, so that unresolved scale-dependent processes represented as a function of some grid-averaged state variable using empirical or statistical relationships. LSMs were originally implemented in numerical weather prediction (NWP) and global climate models (GCMs) in order to provide interactive lower boundary conditions for the atmospheric radiation and turbulence parameterization schemes over continental land surfaces. In the

Published by Copernicus Publications on behalf of the European Geosciences Union.

844 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

past 2 decades, LSMs have evolved considerably to include more biogeochemical and biogeophysical processes in order to meet the growing demands of both the research and the user communities (Pitman, 2003; van den Hurk et al., 2011).A growing number of state-of-the-art LSMs, which are used in coupled atmospheric models for operational numerical weather prediction (Ek et al., 2003; Boussetta et al., 2013), climate modeling (Oleson et al., 2010; Zhang et al., 2015), or both (Best et al., 2011; Masson et al., 2013), represent most or all of the following processes: photosynthesis and the associated carbon uxes, multi-layer soil water and heat transfer, vegetation phenology and dynamics (biomass evolution, net primary production), sub-grid lateral water transfer, river routing, atmospherelake exchanges, snowpack dynamics, and near-surface urban meteorology. Some LSMs also include processes describing the nitrogen cycle (Castillo et al., 2012), groundwater exchanges (Vergnes et al., 2014), aerosol surface emissions (Cakmur et al., 2004), isotopes (Braud et al., 2005), and the representation of human impacts on the hydrological cycle in terms of irrigation (de Rosnay et al., 2003) and ground water extraction (Pokhrel et al., 2015), to name a few.

As a part of the trend in LSM development, there have been ongoing efforts to improve the representation of the land surface processes in the Interactions between the soil biosphereatmosphere (ISBA) LSM within the EXternalized SURFace (SURFEX; Masson et al.,2013) model platform.The original two-layer ISBA forcerestore model (Noilhan and Planton, 1989) consists in a single bulk soil layer (generally having a thickness on the order of 50 cm to several meters) coupled to a supercially thin surface composite soilvegetationsnow layer. Thus, the model simulates fast processes that occur at sub-diurnal timescales, which are pertinent to short-term numerical weather prediction, and it provides a longer-term water storage reservoir, which provides a source for transpiration, a time lter for water reaching a hydro-graphic network, and a certain degree of soil moisture memory in the ground amenable to longer-term forecasts and climate modeling. Additional modications were made to this scheme over the last decade to include soil freezing (Boone et al., 2000; Giard and Bazile, 2000), which improved hydrological processes (Mahfouf and Noilhan, 1996; Boone et al., 1999; Decharme and Douville, 2006). This scheme was based on the pioneering work of Deardorff (1977) and it has proven its value for coupled landatmosphere research and applications since its inception. For example, it is currently used for research within the mesoscale non-hydrostatic research model (Meso-NH) (Lafore et al., 1998). It is also used within the operational high-resolution short-term numerical weather prediction at Mto-France within the limited area model AROME (Seity et al., 2011) and by HIRLAM countries within the ALADINHIRLAM system as the HARMONIEAROME model conguration (Bengtsson et al., 2017). Finally, it is used for climate research within the global climate model (GCM) Ac-

tion de Researche Petite Echelle Grande Echelle (ARPEGE-climat; Voldoire et al., 2013) and by HIRLAM countries within the ALADINHIRLAM system as HARMONIEAROME and HARMONIEALARO Climate congurations (Lind et al., 2016).

1.1 Rationale for improved vegetation processes

Currently, many LSMs are pushing towards improved realism owing to an increasing number of observations at the local scale, constantly improving satellite data sets and the associated methodologies to best exploit such data, improved computing resources, and in response to the user community via climate services (and seasonal forecasts, drought indexes, etc. . . ). In the SURFEX context, the forcerestore approach has been replaced in recent years by multi-layer explicit physically based options for sub-surface heat transfer (Boone et al., 2000; Decharme et al., 2016), soil hydrological processes (Boone et al., 2000;Decharme et al., 2011, 2016), and the composite snow-pack (Boone and Etchevers, 2001; Decharme et al., 2016).These new schemes have recently been implemented in the operational distributed hydro-meteorological hindcast system SAFRANISBAMODCOU (Habets et al., 2008), Meso-NH, and ARPEGE-climat and ALADINHIRLAM HARMONIEAROME and HARMONIEALARO Climate congurations. The representation of vegetation processes in SURFEX has also become much more sophisticated in recent years, including photosynthesis and respiration (Calvet et al., 1998), carbon allocation to biomass pools (Calvet and Soussana, 2001; Gibelin et al., 2006), and soil carbon cycling (Joetzjer et al., 2015). However, for a number of reasons it has also become clear that we have reached the conceptual limits of using of a composite soilvegetation scheme within ISBA and there is a need to explicitly separate the canopy vegetation from the soil surface:

in order to distinguish the soil, snow, and vegetation surface temperatures since they can have very different amplitudes and phases in terms of the diurnal cycle, and therefore accounting for this distinction facilitates (at least conceptually) incorporating remote sensing data, such as satellite-based thermal infrared temperatures (e.g., Anderson et al., 1997), into such models;

as it has become evident that the only way to simulate the snowpack beneath forests in a robust and a physically consistent manner (i.e., reducing the dependence of forest snow cover on highly empirical and poorly constrained snow fractional cover parameterizations) and including certain key processes (i.e., canopy interception and unloading of snow) is to include a forest canopy above or buried by the ground-based snow-pack (e.g., Rutter et al., 2009);

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 845

for accurately modeling canopy radiative transfer, within or below canopy turbulent uxes and soil heat uxes;

to make a more consistent photosynthesis and carbon allocation model (including explicit carbon stores for the vegetation, litter, and soil in a consistent manner);

to allow the explicit treatment of a ground litter layer, which has a signicant impact on ground heat uxes and soil temperatures (and freezing) and, by extension, the turbulent heat uxes.

In response to this issue, a collaboration began in 2008 between the high-resolution limited area model (HIRLAM) consortium and Mto-France with the intention to develop an explicit representation of the vegetation in ISBA under the SURFEX platform. A new parameterization has been developed called the ISBA multi-energy balance (MEB) in order to account for all of the above issues.

MEB is based on the classic two-source model for snow-free conditions, which considers explicit energy budgets (for computing uxes) for the soil and the vegetation, and it has been extended to a three-source model in order to include an explicit representation of snowpack processes and their interactions with the ground and the vegetation. The vegetation canopy is represented using the big-leaf method, which lumps the entire vegetation canopy into a single effective leaf for computing energy budgets and the associated uxes of heat, moisture, and momentum. One of the rst examples of a two-source model designed for atmospheric model studies is Deardorff (1978), and further renements to the vegetation canopy processes were added in the years that followed leading to fairly sophisticated schemes, which are similar to those used today (e.g., Sellers et al., 1986). The two-source big-leaf approach has been used extensively within coupled regional and global scale landatmosphere models (Xue et al., 1991;Sellers et al., 1996; Dickinson et al., 1998; Lawrence et al., 2011; Samuelsson et al., 2011). In addition, more recently multi-layer vegetation schemes have also been developed for application in GCMs (Bonan et al., 2014; Ryder et al., 2016).

ISBA-MEB has been developed taking the same strategy that has been used historically for ISBA: inclusion of the key rst-order processes while maintaining a system that has minimal input data requirements and computational cost while being consistent with other aspects of ISBA (with the ultimate goal of being used in coupled operational numerical weather forecast and climate models, and spatially distributed monitoring and hydrological modeling systems). In 2008, one of the HIRLAM partners, the Swedish Meteorological and Hydrological Institute (SMHI), had already developed and applied an explicit representation of the vegetation in the Rossby Centre Regional Climate Model (RCA3) used at SMHI (Samuelsson et al., 2006, 2011). This representation was introduced into the operational NWP HIRLAMv7.3 system, which became opera-

tional in 2010. In parallel, the dynamic vegetation model LJP-GUESS was coupled to RCA3 as RCAGUESS (Smith et al., 2011), making it possible to simulate complex biogeo-physical feedback mechanisms in climate scenarios. Since then RCAGUESS has been applied over Europe (Wramneby et al., 2010), Africa (Wu et al., 2016), and the Arctic (Zhang et al., 2014). The basic principles developed by SMHI have been the foundation since the explicit representation of the vegetation was introduced in ISBA and SURFEX, but now in a more general and consistent way. Implementation of canopy turbulence scheme, longwave radiation transmission function, and snow interception formulations in MEB largely follows the implementation used in RCA3 (Samuelsson et al., 2006, 2011). In addition, we have taken this opportunity to incorporate several new features into ISBA-MEB compared to the original SMHI scheme:

a snow fraction that can gradually bury the vegetation vertically thereby transitioning the turbulence coupling from the canopy air space directly to the atmosphere (using a fully implicit numerical scheme);

the use of the detailed solar radiation transfer scheme that is a multi-layer model that considers two spectral bands, direct and diffuse ux components, and the concept of sunlit and shaded leaves, which was primarily developed to improve the modeling of photosynthesis within ISBA (Carrer et al., 2013);

a more detailed treatment of canopy snow interception and unloading processes and a coupling with the ISBA physically based multi-layer snow scheme;

a reformulation of the turbulent exchange coefcients within the canopy air space for stable conditions, such as over a snowpack;

a fully implicit Jacobean matrix for the longwave uxes from multiple surfaces (snow, below-canopy snow-free ground surface, vegetation canopy);

all of the energy budgets are numerically implicitly coupled with each other and with the atmosphere using the coupling method adapted from Best et al. (2004), which was rst proposed by Polcher et al. (1998);

an explicit forest litter layer model (which also acts as the below-canopy surface energy budget when litter covers the soil).

This paper is the rst of two parts: in part one, the ISBAMEB model equations, numerical schemes, and theoretical background are presented. In part two, a local-scale evaluation of the new scheme is presented along with a detailed description of the new forest litter scheme (Napoly et al., 2016). An overview of the model is given in the next section, followed by conclusions.

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

846 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

2 Model description

SURFEX uses the tile approach for the surface, and separate physics modules are used to compute surfaceatmosphere exchange for oceans or seas, lakes, urbanized areas, and the natural land surface (Masson et al., 2013). The ISBA LSM is used for the latter tile, and the land surface is further split into upwards of 12 or 19 patches (refer to Table 1), which represent the various land cover and plant functional types.Currently, forests make up eight patches for the 19-class option, and three for the 12-class option. The ISBA-MEB (referred to hereafter simply as MEB) option can be activated for any number of the forest patches. By default, MEB is coupled to the multi-layer soil (ISBA-DF: explicit DiFfusion equation for heat and Richards equation for soil water ow; Boone et al.,2000; Decharme et al., 2011) and snow (ISBA-ES: multi-layer Explicit Snow processes with 12 layers by default; Boone and Etchevers, 2001, Decharme et al., 2016) schemes. These schemes have been recently updated (Decharme et al., 2016) to include improved physics and increased layering (14 soil layers by default). MEB can also be coupled to the simple three-layer soil forcerestore (3-L) option (Boone et al., 1999) in order to be compatible with certain applications, which have historically used 3-L, but by default, it is coupled with ISBA-DF since the objective is to move towards a less conceptual LSM.

A schematic diagram illustrating the various resistance pathways corresponding to the turbulent uxes for the three fully (implicitly) coupled surface energy budgets is shown in Fig. 1. The water budget prognostic variables are also indicated. Note that the subscripts, which are used to represent the different prognostic and diagnostic variables and the aerodynamic resistance pathways, are summarized in Table 2. The canopy bulk vegetation layer is represented using green, the canopy-intercepted snow and ground-based snow-pack are shaded using turquoise, and the ground layers are indicated using dark brown at the surface, which fade to white with increasing depth.

There are six aerodynamic resistance, Ra (s1), pathways dened as being between (i) the non-snow-buried vegetation canopy and the canopy air, Ravgc, (ii) the non-snow-

buried ground surface (soil or litter) and the canopy air, Ragc, (iii) the snow surface and the canopy air, Ranc,

(iv) the ground-based snow-covered part of the canopy and the canopy air, Ravnc, (v) the canopy air with the overly

ing atmosphere, Raca), and (vi) the ground-based snow sur

face (directly) with the overlying atmosphere, Rana. Pre

vious papers describing ISBA (Noilhan and Planton, 1989;

Mahfouf and Noilhan, 1991) expressed heat uxes using a dimensionless heat and mass exchange coefcient, CH; however, for the new MEB option, it is more convenient to express the different uxes using resistances (s m1), which are related to the exchange coefcient as Ra = 1/(Va CH), where

Va represents the wind speed at the atmospheric forcing level (indicated by using the subscript a) in m s1.

Figure 1. A schematic representation of the turbulent aerodynamic resistance, Ra, pathways for ISBA-MEB. The prognostic temperature, liquid water, and liquid water equivalent variables are shown.

The canopy air diagnostic variables are enclosed by the red-dashed circle. The ground-based snowpack is indicated using turquoise, the vegetation canopy is shaded green, and ground layers are colored dark brown at the surface, fading to white with depth. Atmospheric variables (lowest atmospheric model or observed reference level) are indicated using the a subscript. The ground snow fraction, png, and canopy-snow-cover fraction, pn , are indicated.

The surface energy budgets are formulated in terms of prognostic equations governing the evolutions of the bulk vegetation canopy, Tv, the snow-free ground surface (soil or litter), Tg, and the ground-based snowpack, Tn (K). The prognostic hydrological variables consist of the liquid soil water content, Wg, equivalent water content of ice, Wgf, snow water equivalent (SWE), Wn, vegetation canopy-intercepted liquid water, Wr, and intercepted snow, Wrn (kg m2). The diagnosed canopy air variables, which are determined implicitly during the simultaneous solution of the energy budgets, are enclosed within the red-dashed circle and represent the canopy air specic humidity, qc (kg kg1), air temperature

Tc, and wind speed Vc. The ground surface specic humidity is represented by qg. The surface snow cover fraction area is represented by png while the fraction of the canopy buried by the ground-based snowpack is dened as p n. The snowpack has Nn layers, while the number of soil layers is dened as

Ng where k is the vertical index (increasing from 1 at the surface downward). The ground and snowpack uppermost layer

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 847

Table 1. Description of the patches for the natural land surface sub-grid tile. The values for the 19-class option are shown in the leftmost three columns, and those for the 12-class option are shown in the rightmost three columns (the name and description are only given if they differ from the 19-class values). MEB can currently be activated for the forest classes: 46 (for both the 12- and 19-class options), and 1317.

Index Name Description Index Name Description

1 NO Bare soil 12 ROCK Rock 23 SNOW Permanent snow or ice 34 TEBD Temperate broad leaf 4 TREE Broad leaf5 BONE Boreal evergreen needle leaf 5 CONI Evergreen needle leaf6 TRBE Tropical evergreen broad leaf 6 EVER Evergreen broad leaf7 C3 C3 crops 78 C4 C4 crops 89 IRR Irrigated crops 910 GRAS Temperate grassland 1011 TROG Tropical grassland 1112 PARK Bog, park, garden 1213 TRBD Tropical broad leaf14 TEBE Temperate evergreen broad leaf15 TENE Temperate evergreen needle leaf16 BOBD Boreal broad leaf17 BOND Boreal needle leaf18 BOGR Boreal grassland19 SHRB Shrubs

temperatures correspond to those used for the surface energy budget (i.e., k = 1).

2.1 Snow fractions

Snow is known to have a signicant impact on heat conduction uxes, owing to its relatively high insulating properties.In addition, it can signicantly reduce turbulent transfer owing to reduced surface roughness, and it has a relatively large surface albedo thereby impacting the surface net radiation budget. Thus, the parameterization of its areal coverage turns out to be a critical aspect of LSM modeling of snowpack atmosphere interactions and sub-surface soil and hydrological processes. The fractional ground coverage by the snow-pack is dened as

png = Wn/Wn,crit 0 png 1

[parenrightbig]

[parenrightbig]

, (1)

where currently the default value is Wn,crit = 1 (kg m2).

Note that this is considerably lower than the previous value of 10 kg m2 used in ISBA (Douville et al., 1995), but this value has been shown to improve the ground soil temperatures, using an explicit snow scheme within ISBA (Brun et al., 2013).

The fraction of the vegetation canopy, which is buried by ground-based snow, is dened as

pn = Dn zhv,b

[parenrightbig]

/ zhv zhv,b

+ png Ggn + n,NnSWnet,n

[parenrightbig]

[parenrightbig]

(0 pn 1), (2) where Dn is the total ground-based snowpack depth (m) and zhvb represents the base of the vegetation canopy (m) (see Fig. 2), which is currently dened as

zhvb = ahv zhv zhv,min

[parenrightbig]

(zhvb 0), (3)

where ahv = 0.2 and the effective canopy base height is set to

zhv,min = 2 (m) for forests. The foliage distribution should be

reconsidered in further development since literature suggests, e.g., Massman (1982), that the foliage is not symmetrically distributed in the crown but skewed upward.

2.2 Energy budget

The coupled energy budget equations for a three-source model can be expressed for a single bulk canopy, a ground-based snowpack, and a underlying ground surface as

Cv @Tv@t =Rnv Hv LEv + Lf [Phi1]v, (4)

Cg,1 @Tg,1@t = 1 png

[parenrightbig][parenleftBigg]

Rng Hg LEg

Gg,1 + Lf [Phi1]g,1, (5)

Cn,1 @Tn,1@t =Rnn Hn LEn n,1SWnet,n + n,1

Gn,1 + Lf [Phi1]n,1, (6)

where Tg,1 is the uppermost ground (surface soil or litter layer) temperature, Tn,1 is the surface snow temperature, and

Tv is the bulk canopy temperature (K). Note that the subscript 1 indicates the uppermost layer or the base of the layer (for uxes) for the soil and snowpack. All of the following ux terms are expressed in W m2. The sensible heat uxes are dened between the canopy air space and the vegetation Hv, the snow-free ground Hg, and the ground-based snow-

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

848 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Table 2. Subscripts used to represent the prognostic and diagnostic variables. In addition, the symbols used to represent the aerodynamic resistance pathways (between the two elements separated by the dash) are also shown (refer also to Fig. 1). These symbols are used throughout the text.

Subscript Name Description

v Vegetation Bulk canopy layerg Ground Temperature or liquid water (for Ng layers)

gf Ground Frozen water (for Ng layers)a Atmosphere At the lowest atmospheric or forcing levelc Canopy air space Diagnosed variablesn Ground-based snowpack For Nn layersng Ground-based snowpack Fractional ground snow coverage n Ground-based snowpack Fractional vegetation snow coverager Interception reservoir Intercepted rain and snow meltwaterrn Interception reservoir Intercepted snow and frozen meltwater or rainvg-c Aerodynamic resistance Non-snow-buried vegetation canopy and canopy air g-c Aerodynamic resistance Non-snow-buried ground surface and canopy airn-c Aerodynamic resistance Snow surface and canopy airvn-c Aerodynamic resistance Ground-based snow-covered canopy and canopy air c-a Aerodynamic resistance Canopy air with overlying atmospheren-a Aerodynamic resistance Ground-based snow surface and overlying atmosphere

pack Hn. In an analogous fashion to the sensible heat ux, the latent heat uxes are dened for the vegetation canopy Ev, the snow-free ground Eg, and the ground-based snow-pack En. The net radiation uxes are dened for the vegetation canopy, ground, and snowpack as Rnv, Rng, and Rnn, respectively. Note that part of the incoming shortwave radiation is transmitted through the uppermost snow layer, and this energy loss is expressed as n,1 SWnet,n, where is the dimensionless transmission coefcient. The conduction uxes between the uppermost ground layer and the underlying soil and the analogue for the snowpack are dened as Gg,1 and

Gn,1, respectively. The conduction ux between the base of the snowpack and the ground surface is dened as Ggn. The last term on the right-hand side (RHS) of Eq. (6), n,1, represents the effective heating or cooling of a snowpack layer caused by exchanges in enthalpy between the surface and sub-surface model layers when the vertical grid is reset (the snow model grid-layer thicknesses vary in time).

The ground-based snow fraction is dened as png. Note that certain terms of Eq. (5) are multiplied by png to make them patch relative (or grid box relative in the case of single-patch mode) since the snow can potentially cover only part of the patch. Within the snow module itself, the notion of png is not used (the computations are snow relative). But note that when simultaneously solving the coupled equations Eqs. (4)(6), Eq. (6) must be multiplied by png since again, snow only covers a fraction of the area: further details are given in Appendices G and I. The formulation for png is described in

Sect. 2.1.

The phase change terms (freezing less melting: expressed in kg m2 s1) for the snow water equivalent intercepted by the vegetation canopy, the uppermost ground layer, and the

uppermost snowpack layer are represented by [Phi1]v, [Phi1]g,1, and [Phi1]n,1, respectively, and Lf represents the latent heat of fusion (J kg1). The computation of [Phi1]g,1 uses the Gibbs free-energy method (Decharme et al., 2016), [Phi1]n,1 is based on available liquid for freezing or cold content for freezing (Boone and Etchevers, 2001), and [Phi1]v is described herein (see Eq. 83).

Note that all of the phase change terms are computed as adjustments to the surface temperatures (after the uxes have been computed); therefore, only the energy storage terms are modied directly by phase changes for each model time step.

The surface ground, snow, and vegetation effective heat capacities, Cg,1, Cv, and Cn,1 (J m2 K1) are dened, re

spectively, as

Cg,1= [Delta1]zg,1 cg,1 (7)

Cv= Cvb + Ci Wr,n + Cw Wr, (8)

Cn,1= Dn,1 cn,1, (9)

where Ci and Cw are the specic heat capacities for solid (2.106 [notdef] 103 J kg1 K1) and liquid water (4.218 [notdef]

103 J kg1 K1), respectively. The uppermost ground-layer thickness is [Delta1]zg,1 (m), and the corresponding heat capacity of this layer is dened as cg 1 (J m3 K1). The uppermost soil layer ranges between 0.01 and 0.03 m for most applications so that the interactions between surface uxes and fast temperature changes in the surface soil layer can be represented. There are two options for modeling the thermal properties of the uppermost ground layer. First, they can be dened using the default ISBA conguration for a soil layer with parameters based on soil texture properties, which can also incorporate the thermal effects of soil organics (Decharme et al., 2016). The second option, which is the default when using MEB, is to model the uppermost

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 849

Figure 2. A schematic sketch illustrating the role of pn , the fraction of the vegetation canopy, which is buried by ground-based snow. In panel (a), the snow is well below the canopy base, zhvb, resulting in pn = 0, and the snow has no direct energy exchange with the

atmosphere. In panel (b), the canopy is partly buried by snow (0 < pn < 1) and the snow has energy exchanges with both the canopy air and the atmosphere. In panel (c), the canopy is fully buried by snow (pn = 1) and the snow has energy exchange only with the atmosphere,

whereas the soil and canopy only exchange with the canopy air space (png < 1). Finally, in panel (d), both png = 1 and pn = 1 so that the

only exchanges are between the snow and the atmosphere.

ground layer as forest litter. The ground surface in forest regions is generally covered by a litter layer consisting of dead leaves and or needles, branches, fruit, and other organic material. Some LSMs have introduced parameterizations for litter (Gonzalez-Sosa et al., 1999; Oge and Brunet, 2002;Wilson et al., 2012), but the approach can be very different from one case to another depending on their complexity. The main goal of this parameterization within MEB is to account for the generally accepted rst-order energetic and hydrological effects of litter; this layer is generally accepted to have a strong insulating effect owing to its particular thermal properties (leading to a relatively low thermal diffusivity), it causes a signicant reduction of ground evaporation (capillary rise into this layer is negligible), and it constitutes an interception reservoir for liquid water, which can also lose water by evaporation. See Napoly et al. (2016) for a detailed description of this scheme and its impact on the surface energy budget.

The canopy is characterized by low heat capacity, which means that its temperature responds fast to changes in uxes. Thus, to realistically simulate diurnal variations in 2 m temperature this effect must be accounted for. Sellers et al. (1986) dened the value as being the heat capacity of 0.2 kg m2 of water per unit leaf area index (m2 m2).

This results in values on the order of 1 [notdef] 104 J m2 K1 for

forest canopies in general. For local-scale simulations, Cvb can be dened based on observational data. In spatially distributed simulations (or when observational data is insufcient), Cvb = 0.2/CV where the vegetation thermal inertia

CV is dened as a function of vegetation class by the SURFEX default physiographic database ECOCLIMAP (Faroux et al., 2013). Note that CV has been determined for the composite soilvegetation scheme, and the factor 0.2 is used to reduce this value to be more representative of vegetation and on the order of the value discussed by Sellers et al. (1986).Numerical tests have shown that using this value, the canopy heat storage is on the order of 10 W m2 at mid-day for a typical mid-latitude summer day for a forest. The minimum veg-

etation heat capacity value is limited at 1 [notdef] 104 (J m2 K1)

in order to model, in a rather simple fashion, the thermal inertia of stems, branches, trunks, etc. The contributions from intercepted snow and rain are incorporated, where Wr,n and

Wr (kg m2) represent the equivalent liquid water content of intercepted canopy snow and liquid water, respectively.

The uppermost snow-layer thickness is Dn,1 (m), and the corresponding heat capacity is represented by cn,1 (Boone and Etchevers, 2001). Note that Dn,1 is limited to values no larger than several centimeters in order to model a reasonable thermal inertia (i.e., in order to represent the diurnal cycle) in a fashion analogous to the soil. For more details, see Decharme et al. (2016).

The numerical solution of the surface energy budget, sub-surface soil and snow temperatures, and the implicit numerical coupling with the atmosphere is described in Appendix I.

2.3 Turbulent uxes

In this section, the turbulent heat and water vapor uxes in Eqs. (4)(6) are described.

2.3.1 Sensible heat uxes

The MEB sensible heat uxes are dened as

Hv =a

(Tv Tc) Ravc

, (10)

Hg =a Tg Tc[parenrightbig]

Ragc

, (11)

Hn =a

(1 pn )

(Tn Tc)

Ranc +

pn (Tn Ta) Rana

, (12)

Hc =a

(Tc Ta) Raca

, (13)

H =a [bracketleftbigg][parenleftBigg]

1 pn png

[parenrightbig]

(Tc Ta)

Raca +

pn png (Tn Ta)

Rana

, (14)

where a represents the lowest atmospheric layer average air density (kg m3). The sensible heat uxes appear in the sur-

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

850 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

face energy budget equations (Eqs. 46). The sensible heat ux from the ground-based snowpack (Eq. 12) is partitioned by the fraction of the vegetation, which is buried by the ground-based snowpack pn , between an exchange between the canopy air space, and the overlying atmosphere (Eq. 2).The heat ux between the overlaying atmosphere and the canopy air space is represented by Hc, and it is equivalent to the sum of the uxes between the different energy budgets and the canopy air space. The total ux exchange between the overlying atmosphere and the surface (as seen by the atmosphere) is dened by H. It is comprised of two components: the heat exchange between the overlying atmosphere and the canopy air space and the part of the ground-based snowpack that is burying the vegetation. This method has been developed to model the covering of low vegetation canopies by a ground-based snowpack. Finally, the nal uxes for the given patch are aggregated using png and pn : the full expressions are given in Appendix C1.

The thermodynamic variable (T : J kg1) is linearly related

to temperature as

Tx

= Bx + Ax Tx, (15)

where x corresponds to one of the three surface temperatures (Tg, Tv, or Tn), canopy air temperature, Tc, or the overlying atmospheric temperature, Ta. The denitions of Ax and Bx

depend on the atmospheric variable in the turbulent diffusion scheme and are usually dened to cast T in the form

of dry static energy, or potential temperature and are determined by the atmospheric model in coupled mode (see Appendix A). The total canopy aerodynamic resistance is comprised of snow buried, Ravnc, and non-snow buried, Ravgc,

resistances from

Ravc =

"(1 pn ) png Ravnc +[parenleftBigg][parenleftBigg]1

for Tc and using Eq. (15) to determine Tc (see Appendix A

for details).

2.3.2 Water vapor uxes

The MEB water vapor uxes are expressed as

Ev =a hsv

(qsatv qc) Ravc

, (17)

Eg =a

qg qc

Ragc, (18)

En =a hsn [bracketleftbigg]

(1 pn )

(qsatin qc)

Ranc +

pn (qsatin qa)

Rana

,

(19)

Ec =a

(qc qa) Raca

, (20)

E =a

[bracketleftbigg][parenleftBigg]

1 pn png

[parenrightbig]

(qc qa) Raca

+pn png hsn (qsati

n

qa)

Rana

. (21)

The vapor ux between the canopy air and the overlying atmosphere is represented by Ec, and the total vapor ux exchanged with the overlying atmosphere is dened as E. The specic humidity (kg kg1) of the overlying atmosphere is represented by qa, whereas qsat and qsati represent the specic humidity at saturation over liquid water and ice, respectively. For the surface specic humidities at saturation, the convention qsatx = qsat (Tx) is used. The same holds true for

saturation over ice so that qsatin = qsati (Tn). The canopy air

specic humidity qc is diagnosed assuming that Ec is balanced by the vapor uxes between the canopy air and each of the three surfaces considered (the methodology for diagnosing the canopy air thermal properties is described in Appendix I, Sect. I3). The effective ground specic humidity is dened as

qg = hsg qsatg + (1 + ha)qc, (22)

where the humidity factors are dened as

hsg = g hug 1 pgf

[parenrightbig][parenleftbigg]

Lv

L

png

R a vgc

#1. (16)

The separation of the resistances is done to mainly account for differences in the roughness length between the buried and non-covered parts of the vegetation canopy; therefore, the primary effect of snow cover is to increase the resistance relative to a snow-free surface assuming the same temperature gradient owing to a lower surface roughness, and thus Ravnc Ravgc. The formulation also provides a con

tinuous transition to the case of vanishing canopy turbulent uxes as the canopy becomes entirely buried (as pn ! 1).

In this case, the energy budget equations collapse into a simple coupling between the snow surface and the overlying atmosphere, and the ground energy budget simply consists in heat conduction between the ground surface and the snow-pack base. The formulations of the resistances between the different surfaces and the canopy airspace and the overlying atmosphere are described in detail in Sect. 2.6. The canopy air temperature, which is needed by different physics routines, is diagnosed by combining Eqs. (10)(14) and solving

[parenrightbigg] +

gf hugf pgf

Ls L

, (23)

ha = g 1 pgf

[parenrightbig][parenleftbigg]

Lv

L

[parenrightbigg] +

. (24)

The latent heats of fusion and vaporization are dened as Ls and Lv (J kg1), respectively. The fraction of the surface layer that is frozen, pgf, is simply dened as the ratio of the liquid water equivalent ice content to the total water content. The average latent heat L is essentially a normalization factor that ranges between Ls and Lv as a function of snow cover and surface soil ice (see Appendix B). The soil coefcient g in Eqs. (23)(24) is dened as

g =

[parenleftbigg]

Ragc

Ragc + Rg

gf pgf

Ls L

gcor, (25)

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

[parenrightbigg]

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 851

where the soil resistance Rg is dened by Eq. (67). Note that the composite version of ISBA did not include an explicit soil resistance term, and therefore this also represents a new addition to the model. This term was found to further improve results for bare-soil evaporation within MEB, and its inclusion is consistent with other similar multi-source models (e.g., Xue et al., 1991). See Sect. 2.6 for further details.The delta function gcor is a numerical correction term that is required owing to the linearization of qsatg and is unity unless both hug qsatg < qc and qsatg > qc, in which case it is set to zero. The surface ground humidity factor is dened using the standard ISBA formulation from Noilhan and Planton (1989) as

hug =

1

2

+ 1 png

[parenrightbig][parenleftbigg]

Ravc

Ravgc

. (28b)

The Halstead coefcients in Eq. (28a) are dened as

hvg =

[parenleftbigg]

Ravgc

Ravgc + Rs

(1 ) + , (29a)

,hvn =

[parenleftbigg]

Ravnc

Ravnc + Rsn

"1 cos

wg,1 w fc,1

[parenrightBigg][bracketrightBigg][parenleftBigg][parenleftBigg]0

(1 ) + . (29b)

The stomatal resistance Rs can be computed using either the Jarvis method (Jarvis, 1976) described by Noilhan and Plan-ton (1989) or a more physically based method that includes a representation of photosynthesis (Calvet et al., 1998). The stomatal resistance for the partially snow-buried portion dened as

Rsn =Rs/

1 min pn , 1 Rs/Rs,max [parenrightbig][bracketrightbig]

Rsn Rs,max

[parenrightbig]

(30)

hug 1

[parenrightbig]

. (26)

In the case of condensation (qsatg < qa), hug = 1 (see Mah

fouf and Noilhan, 1991, for details). The effective eld capacity w fc,1 is computed relative to the liquid water content of the uppermost soil layer (it is adjusted in the presence of soil ice compared to the default eld capacity). The analogous form holds for the humidity factor over the frozen part of the surface soil layer, hugf, with wg,1 and w fc,1 replaced by wgf,1 and w fcf,1 (m3 m3) in Eq. (26), respectively (Boone et al., 2000). Note that it would be more accurate to use qsati in place of qsat for the sublimation of the canopy-intercepted snow and the soil ice in Eqs. (17)(18), respectively, but this complicates the linearization and this has been neglected for now. The snow factor is dened as hsn = Ls/L. This factor

can be modied so that En includes both sublimation and evaporation (Boone and Etchevers, 2001), but the impact of including a liquid water ux has been found to be negligible; thus, for simplicity only sublimation is accounted for currently.

The leading coefcient for the canopy evapotranspiration is dened as

hsv = (1 pnv)hsvg (Lv/L) + pnv hsvn (Ls/L), (27)

where pnv is an evaporative efciency factor that is used to partition the canopy interception storage mass ux between evaporation of liquid water and sublimation (see Eq. 79). When part of the vegetation canopy is buried (i.e., pn > 0), a different roughness is felt by the canopy air space so that a new resistance is computed over the pn -covered part of the canopy as is done for sensible heat ux. This is accounted for by dening

hsvg =png (1 pn )

[parenleftbigg]

Ravc

Ravnc

so that the effect of coverage by the snowpack is to increase the canopy resistance. Note that when the canopy is not partially or fully buried by ground-based snowpack (pn = 0)

and does not contain any intercepted snow (pnv = 0), the

leading coefcient for the canopy evapotranspiration simplies to the Halstead coefcient from the composite version of ISBA (Mahfouf and Noilhan, 1991)

hsv =

[parenleftbigg]

Ravgc

Ravgc + Rs

(1 ) +

(pn = 0 and pnv = 0). (31)

The fraction of the vegetation covered by water is and is described in Sect. 2.8.2.

The evapotranspiration from the vegetation canopy, Ev, is comprised of three components:

Ev = Etr + Er + Ern, (32)

where the transpiration, evaporation from the canopy liquid water interception store, and sublimation from the canopy snow interception store are represented by Etr, Er, and Ern, respectively. The expressions for these uxes are given in Appendix C.

2.4 Radiative uxes

The Rn terms in Eqs. (4)(6) represent the surface net radiation terms (longwave and shortwave components)

Rnx = SWnet,x + LWnet,x, (33)

where x = n, g, or v. The total net radiation of the surface is Rn =Rnn + Rng + Rnv

= SW # SW " +LW # LW ", (34)

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

[parenrightbigg]

hvn

+ 1 png

[parenrightbig][parenleftbigg]

Ravc

Ravgc

[parenrightbigg]

hvg, (28a)

hsvn =png (1 pn )

[parenleftbigg]

Ravc

Ravnc

[parenrightbigg]

852 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

where the total down-welling solar (shortwave) and atmospheric (longwave) radiative uxes (W m2) at the top of the canopy or snow surface (in the case snow is burying the vegetation) are represented by SW # and LW #, respec

tively. The total upwelling (towards the atmosphere) short-wave and longwave radiative uxes SW " and LW ", respec

tively, are simply dened as the downward components less the total surface net radiative uxes (summed over the three surfaces). The effective total surface albedo and surface radiative temperature (and emissivity) can then be diagnosed (see the Sect. 2.4.2) for coupling with the host atmospheric model. The n is dened as the solar radiation transmission at the base of a snowpack layer, so for a sufciently thin snow-pack, solar energy penetrating the snow to the underlying ground surface is expressed as n,NnSWnetn, where Nn represents the number of modeled snowpack layers (for a deep snowpack, this term becomes negligible).

2.4.1 Shortwave radiative uxes

The total land surface shortwave energy budget can be shown to satisfy

SW #= SWnetg + SWnetv + SWnetn + SW ", (35)

where SWnetg, SWnetv, and SWnetn represent the net short-wave terms for the ground, vegetation canopy, and the ground-based snowpack. The effective surface albedo (which may be required by the atmospheric radiation scheme or for comparison with satellite-based data) is diagnosed as

s = SW " /SW # . (36)

The multi-level transmission computations for direct and diffuse radiation are from Carrer et al. (2013). The distinction between the visible (VIS) and near-infrared (NIR) radiation components is important in terms of interactions with the vegetation canopy. Here, we take into account two spectral bands for the soil and the vegetation, where visible wavelengths range from approximately 0.3 to 0.7 [notdef]106 m,

and NIR wavelengths range from approximately 0.7 to1.4 [notdef]106 m. The spectral values for the soil and the veg

etation are provided by ECOCLIMAP (Faroux et al., 2013) as a function of vegetation type and climate.

The effective all-wavelength ground (below-canopy) albedo is dened as

gn = png n + 1 png

[parenrightbig]

1/4, (39)

where is the StefanBoltzmann constant, and [epsilon1]s represents the effective surface emissivity. In Eq. (39), there are two that are known (LW uxes) and two that are unknown (Trad and [epsilon1]s). Here we opt to pre-dene [epsilon1]s in a manner that is consistent with the various surface contributions as

[epsilon1]s = png [epsilon1]sn + 1 png

[parenrightbig]

g, (37)

where g represents the ground albedo.

The ground-based snow albedo, n, is prognostic and depends on the snow grain size. It currently includes up to three spectral bands (Decharme et al., 2016); however, when coupled to MEB, only the two aforementioned spectral bands are currently considered for consistency with the vegetation and soil.

The effective canopy albedo, v, represents the combined canopy vegetation, v, and intercepted snow albedos. Currently, however, we assume that v = v, which is based

on recommendations by Pomeroy and Dion (1996). They showed that multiple reections and scattering of light from patches of intercepted snow together with a high probability of reected light reaching the underside of an overlying branch implied that trees actually act like light traps. Thus, they concluded that intercepted snow had no signicant inuence on the shortwave albedo or the net radiative exchange of boreal conifer canopies.

In addition to baseline albedo values required by the radiative transfer model for each spectral band, the model requires the direct and diffusive downwelling solar components. The diffuse fraction can be provided by observations (ofine mode) or a host atmospheric model. For the case when no diffuse information is provided to the surface model, the diffuse fraction is computed using the method proposed by Erbs et al. (1982).

2.4.2 Longwave radiative uxes

The longwave radiation scheme is based on a representation of the vegetation canopy as a plane-parallel surface. The model considers one reection with three reecting surfaces (ground, ground-based snowpack, and the vegetation canopy; a schematic is shown in Appendix E). The total land surface longwave energy budget can be shown to satisfy

LW #= LWnetg + LWnetv + LWnetn + LW ", (38)

where LWnetg, LWnetv, and LWnetn represent the net long-wave terms for the ground, vegetation canopy, and the ground-based snowpack. The effective surface radiative temperature (which may be required by the atmospheric radiation scheme or for comparison with satellite-based data) is diagnosed as

Trad =

LW " LW # (1 [epsilon1]s) [epsilon1]s

[epsilon1]sg. (40)

The canopy-absorption-weighted effective snow and ground emissivities are dened, respectively, as

[epsilon1]sn =n LW [epsilon1]v + (1 n LW) [epsilon1]n, (41) [epsilon1]sg =g LW [epsilon1]v + 1 g LW

[parenrightbig]

[epsilon1]g, (42)

where [epsilon1]v, [epsilon1]g, and [epsilon1]n represent the emissivities of the vegetation, snow-free ground, and the ground-based snowpack, respectively. The ground and vegetation emissivities are given

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 853

by ECOCLIMAP are vary primarily as a function of vegetation class for spatially distributed simulations, or they can be prescribed for local scale studies. The snow emissivity is currently dened as [epsilon1]n = 0.99. The effect of longwave ab

sorption through the non-snow-buried part of the vegetation canopy is included as

n LW =

1 png pn 1 png[parenrightbig][bracketrightbig]

LW +

png

+pn 1 png

[parenrightbig][bracketrightbig]

fLW, (43)

g LW =

LW + png (1 pn )fLW, (44)

where the canopy absorption is dened as

LW = 1 exp(LW LAI) = 1 [notdef]v (45)

and LW represents a longwave radiation transmission factor that can be species (or land classication) dependent, [notdef]v is dened as a vegetation view factor, and LAI represents the leaf area index (m2 m2). The absorption over the understory snow-covered fraction of the grid box is modeled quite simply from Eq. (45)

fLW = 1 exp

LW LAI(1 pn )

[bracketrightbig]

= 1 exp[LW LAIn] (46)

so that transmission is unity (no absorption or reection by

the canopy: LW = fLW = 0) when pn = 1 (i.e., when the

canopy has been buried by snow); LAIn is used to represent the LAI, which has been reduced owing to burial by the snowpack. From Eqs. (40)(44), it can be seen that when there is no snowpack (i.e., png = 0 and pn = 0), then the ef

fective surface emissivity is simply an absorption-weighted soilvegetation value dened as [epsilon1]s = LW [epsilon1]v + (1 LW) [epsilon1]g.

See Appendix E for the derivation of the net longwave radiation terms in Eq. (38).

2.5 Heat conduction uxes

The sub-surface snow and ground heat conduction uxes are modeled using Fouriers law (G = @T /@z). The heat con

duction uxes in Eqs. (5)(6) are written in discrete form as

Gg,1=

2 Tg,1 Tg,2

1 png (1 pn )

[bracketrightbig]

1. (50)

The parameterization of the bulk canopy aerodynamic conductance gav between the canopy and the canopy air is based on Choudhury and Monteith (1988). It is dened as

gav =

2LAIaav

[prime]v

[1 exp([prime]v/2)], (51)

where uhv represents the wind speed at the top of the canopy (m s1),and lw represents the leaf width (m: see Table 3). The remaining parameters and their values are dened in Table 3. The conductance accounting for the free convection correction from Sellers et al. (1986) is expressed as

g av = [bracketleftBigg]

LAI

890

uhv lw

1/2

(Tv Tc). (52)

Note that this correction is only used for unstable conditions.

The effect of snow burying the vegetation impacts the aerodynamic resistance of the canopy is simply modeled by modifying the LAI using

LAIn = LAI(1 pn ). (53)

The LAIn is then used in Eq. (50) to compute Ravnc, and

this resistance is limited to 5000 s m1 as LAIn ! 0.

2.6.2 Aerodynamic resistance between the ground and the canopy air

The resistance between the ground and the canopy air space is dened as

Ragc = Ragn/ H , (54)

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

Tv Tc lw

1/4

[bracketrightBigg]

[Delta1]zg,1/ g,1

[parenrightbig]

+ [Delta1]zg,2/ g,2

[parenrightbig]

= [Lambda1]g,1 Tg,1 Tg,2

,

(47)

Gn,1=

2 Tn,1 Tn,2

Dn,1/ n,1

[parenrightbig]

+ Dn,2/ n,2

[parenrightbig]

= [Lambda1]n,1 Tn,1 Tn,2

, (48)

Ggn=

2 Tn,Nn Tg,1

= [Lambda1]g,n Tn,Nn Tg,1

, (49)

where Ggn represents the snowground inter-facial heat ux which denes the snow scheme lower boundary condition. All of the internal heat conduction uxes (k = 2,N 1) use

the same form as in Eq. (48) for the snow (Boone and Etchevers, 2001) and Eq. (47) for the soil (Boone et al., 2000;

Decharme et al., 2011). The heat capacities and thermal conductivities g for the ground depend on the soil texture, organic content (Decharme et al., 2016), and potentially on the thermal properties of the forest litter in the uppermost layer (Napoly et al., 2016); all of the aforementioned properties depend on the water content. The snow thermal property parameterization is described in Decharme et al. (2016).

2.6 Aerodynamic resistances

The resistances between the surface and the overlying atmosphere, Rana and Raca, are based on Louis (1979) modi

ed by Mascart et al. (1995) to account for different roughness length values for heat and momentum as in ISBA: the full expressions are given in Noilhan and Mahfouf (1996).

2.6.1 Aerodynamic resistance between the bulk vegetation layer and the canopy air

The aerodynamic resistance between the vegetation canopy and the surrounding airspace can be dened as

Ravgc = gav + g av

Dn,Nn/ n,Nn

[parenrightbig]

+ [Delta1]zg,1/ g,1

[parenrightbig]

854 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Table 3. Surface vegetation canopy turbulence parameters that are constant.

Symbol Denition Unit Value Reference Comment

aav Canopy conductance scale factor m s1/2 0.01 Choudhury and Monteith (1988) Eq. (26) [prime]

v Attenuation coeff. for wind 3 Choudhury and Monteith (1988) p. 386 lw Leaf width m 0.02v Attenuation coeff. for mom. 2 Choudhury and Monteith (1988) p. 386 z0g Roughness of soil surface m 0.007[notdef]L RossGoudriaan leaf angle dist. 0.12 Monteith (1975) p. 26ul Typical local wind speed m s1 1 Sellers et al. (1996) Eq. (B7)

Kinematic viscos. of air m2 s1 0.15 [notdef] 10

4

The friction velocity at the top of the vegetation canopy is dened as

u hv =

k uhv

ln

where Ragn is the default resistance value for neutral conditions. The stability correction term H depends on the canopy structural parameters, wind speed, and temperature gradient between the surface and the canopy air. The aerodynamic resistance is also based on Choudhury and Monteith (1988). It is assumed that the eddy diffusivity K (m2 s1) in the vegetation layer follows an exponential prole:

K (z) = K (zhv) exp

v

1

z

zhv[parenrightbigg][bracketrightbigg]

, (55)

where zhv represents the canopy height. Integrating the reciprocal of the diffusivity dened in Eq. (55) from z0g to d +z0v

yields

Ragn =

zhv

v K (zhv)

(zhv d)/z0v[bracketrightbig]

, (61) where the wind speed at the top of the canopy is

uhv = fhv Va (62)

and Va represents the wind speed at the reference height za above the canopy. The canopy height is dened based on vegetation class and climate within ECOCLIMAP as a primary parameter. It can also be dened using an external dataset, such as from a satellite-derived product (as a function of space and time). The vegetation roughness length for momentum is then computed as a secondary parameter as a function of the vegetation canopy height. The factor fhv ( 1) is a stability-dependent adjustment factor (see Ap

pendix D).The dimensionless height scaling factor is dened as

z =

exp

v

1

z0g zhv [parenrightbigg][bracketrightbigg]

exp

v

1

d + z0v zhv

[parenrightbigg][bracketrightbigg]

. (56)

The diffusivity at the canopy top is dened as

K (zhv) = k u hv (zhv d). (57)

The von Karman constant k has a value of 0.4. The displacement height is dened as Choudhury and Monteith (1988)

d = 1.1zhv ln

h1 + (cd LAI)1/4[bracketrightBig]

, (58)

where the leaf drag coefcient cd,is dened from Sellers et al. (1996)

cd = 1.328

[bracketleftbigg]

2 Re1/2

(zhv d)

zr (z 1). (63)

The reference height is dened as zr = za d for simulations

where the reference height is sufciently above the top of the vegetation canopy. This is usually the case for local scale studies using observation data. When MEB is coupled to an atmospheric model, however, the lowest model level can be below the canopy height; therefore, for coupled model simulations zr = max(za, zhv d + zmin) where zmin = 2 (m).

Finally, the stability correction factor from Eq. (54) is dened as

H= (1 ahv Ri)1/2 (Ri 0), (64a)

1

= 1 + b Ri(1 + c Ri)1/2 [bracketleftbigg]

[bracketrightbigg] +

0.45

[bracketleftbigg]

1 (1 [notdef]L)

1.6, (59)

1 +

[parenleftbigg]

Ri Ri,crit

(fz0 1) [bracketrightbigg]

Ri > 0 and Ri Ri,crit

[parenrightbig]

, (64b)

where [notdef]L represents the RossGoudriaan leaf angle distribution function, which has been estimated according to Monteith (1975) (see Table 3), and Re is the Reynolds number dened as

Re =

ul lw

. (60)

= 1 + b Ri(1 + c Ri)1/2[parenleftBigg][parenleftBigg]R

fz0

i > Ri,crit

[parenrightbig]

, (64c)

where the Richardson number is dened as

Ri =

g zhv (Ts Tc)

Ts uhv2 . (65)

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 855

Note that strictly speaking, the temperature factor in the denominator should be dened as (Ts + Tc)/2, but this has only

a minor impact for our purposes. The critical Richardson number, Ri,crit, is set to 0.2. This parameter has been dened assuming that some turbulent exchange is likely always present (even if intermittent), but it is recognized that eventually a more robust approach should be developed for very stable surface layers (Galperin et al., 2007). The expression for unstable conditions (Eq. 64a) is from Sellers et al. (1996) where the structural parameter is dened as ahv = 9.

It is generally accepted that there is a need to improve the parameterization of the exchange coefcient for extremely stable conditions typically encountered over snow (Niu and Yang, 2004; Andreadis et al., 2009). Since the goal here is not to develop a new parameterization, we simply modify the expression for stable conditions by using the standard function from ISBA. The standard ISBA stability correction for stable conditions is given by Eq. (64c), where b = 15 and c = 5

(Noilhan and Mahfouf, 1996). The factor that takes into account differing roughness lengths for heat and momentum is dened as

fz0 =

ln zhv/z0g

ln zhv/z0gh

, (66)

where z0gh is the ground roughness length for scalars. The weighting function (i.e., ratio of Ri to Ri,crit) in Eq. (64b) is used in order to avoid a discontinuity at Ri = 0 (the rough

ness length factor effect vanishes at Ri = 0) in Eq. (64c). An

example of Eq. 64c is shown in Fig. 3 using the z0g from Table 3, and for z0gh/z0g of 0.1 and 1.0. Finally, the resistance between the ground-based snowpack Ranc and the canopy

air use the same expressions as for the aerodynamic resistance between the ground and the canopy air outlined herein, but with the surface properties of the snowpack (namely the roughness length and snow surface temperature).

2.6.3 Ground resistance

The soil resistance term is dened based on Sellers et al. (1992) as

Rg = exp

aRg bRg

wg/wsat

Figure 3. Stability correction term is shown using the Sellers formulation for Ri 0 while the function for stable conditions adapted

from ISBA (Ri > 0) for two ratios of z0g/z0gh. The ground surface roughness length is dened in Table 3.

currently these values are used, in part, since the litter formulation is the default conguration for MEB for forests as it generally gives better surface uxes (Napoly et al., 2016).

2.7 Water budget

The governing equations for (water) mass for the bulk canopy, and surface snow and ground layers are written as

@Wr

@t = Prv + max(0, Etr) Er Drv [Phi1]v, (68) @Wrn

@t = In Un Ern + [Phi1]v (69)

png @Wn,1

@t = Ps In + Un + png

Pr Prv + Drv Fnl,1 En + [Phi1]n,1 + nl,1

[parenrightbig]

, (70)

w [Delta1]zg,1 @wg,1

@t = Pr Prv + Drv Eg

[parenrightbig][parenleftBigg]1

png

[parenrightbig]

[parenrightbig][bracketrightbig]

+png Fnl,Nn R0 Fg,1 [Phi1]g,1, (71)

w [Delta1]zg,1 @wgf,1

@t = [Phi1]g,1 Egf 1 png

[parenrightbig]

. (67)

The coefcients are aRg = 8.206 and bRg = 4.255, and the

vertically averaged volumetric water content and saturated volumetric water content are given by wg and wsat, respectively. The averaging is done from one to several upper layers. Indeed, the inclusion of an explicit ground surface energy budget makes it more conceptually straightforward to include a ground resistance compared to the original composite soil vegetation surface. The ground resistance is often used as a surrogate for an additional resistance arising due to a forest litter layer, therefore the soil resistance is set to zero when the litter-layer option is activated. Finally, the coefcients aRg and bRg were determined from a case study for a specic location, and could possibly be location dependent. But

, (72)

where Wr and Wrn represent the vegetation canopy water stores (intercepted water) and the intercepted snow and frozen water (all in kg m2), respectively. Wn,1 represents the snow liquid water equivalent (SWE) for the uppermost snow layer of the multi-layer scheme. The soil liquid water content and water content equivalent of frozen water are dened as wg and wgf, respectively (m3 m3).

The interception reservoir Wr is modeled as single-layer bucket, with losses represented by evaporation Er, and canopy drip Drv of liquid water that exceeds a maximum holding capacity (see Sect. 2.8.2 for details). Sources include condensation (negative Er and Etr) and Prv, which repre-

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

856 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

sents the intercepted precipitation. The positive part of Etr is extracted from the sub-surface soil layers as a function of soil moisture and a prescribed vertical root zone distribution (Decharme et al., 2016). This equation is the same as that used in ISBA, except for the addition of the phase change term, [Phi1]v (kg m2 s1). This term has been introduced owing to the introduction of an explicit canopy snow interception reservoir Wrn; the canopy snow and liquid water reservoirs can exchange mass via this term which is modeled as melt less freezing. The remaining rainfall (Pr Prv) is parti

tioned between the snow-free and snow-covered ground surface, where Pr represents the total grid cell rainfall rate. The canopy snow interception is more complex, and represents certain baseline processes such as snow interception In and unloading Un; see Sect. 2.8.1 for details.

The soil water and snow liquid water vertical uxes at the base of the surface ground and snow are represented, respectively, by Fg,1 using Darcys Law and by Fnl,1 using a tipping-bucket scheme (kg m2 s1). The liquid water ux at the base of the snowpack Fnl,Nn is directed downward into the soil and consists in the liquid water in excess of the lowest model liquid water-holding capacity. A description of the snow and soil schemes are given in (Boone and Etchevers, 2001) and (Decharme et al., 2011), respectively. R0 is the surface runoff. It accounts for sub-grid heterogeneity of precipitation, soil moisture, and for when potential inltration exceeds a maximum rate (Decharme and Douville, 2006).The soil liquid water equivalent ice content can have some losses owing to sublimation in the uppermost soil layer Egf but it mainly evolves owing to phase changes from soil water freezethaw [Phi1]g. The remaining symbols in Eqs. (68)(69)

are dened and described in Sect. 2.8.2 and 2.8.1.

2.8 Precipitation interception

2.8.1 Canopy snow interception

The intercepted snow mass budget is described by Eq. (69), while the energy budget is included as a part of the bulk canopy prognostic equation (Eq. 4). The positive mass contributions acting to increase intercepted snow on canopy are snowfall interception In, water on canopy that freezes [Phi1]v < 0, and sublimation of water vapor to ice Ern < 0. Unloading Un, sublimation Ern > 0, and snowmelt [Phi1]v > 0, are the sinks. All of the terms are in kg m2 s1. It is assumed that intercepted rain and snow can co-exist on the canopy. The intercepted snow is assumed to have the same temperature as the canopy Tv; thus, there is no advective heat exchange with the atmosphere that simplies the equations. For simplicity, when intercepted water on the canopy freezes, it is assumed to become part of the intercepted snow.

The parameterization of interception efciency is based upon Hedstrom and Pomeroy (1998). It determines how much snow is intercepted during the time step and is dened

as

In,v,0 = W rn Wrn

[parenrightbig][bracketleftbig]1

exp kn,v Ps [Delta1]t

[parenrightbig][bracketrightbig]

, (73)

where Wrn is the maximum snow load allowed, Ps the frozen precipitation rate, and kn,v a proportionality factor.

kn,v is a function of Wrn and the maximum plan area of the snowleaf contact area per unit area of ground Cn,vp:

kn,v =

Cn,vp

Wrn . (74)

For a closed canopy, Cn,vp would be equal to one, but for a partly open canopy it is described by the relationship

Cn,vp =

Cn,vc

1 Cc uhv zhv/(wn Jn)

, (75)

where Cn,vc is the canopy coverage per unit area of ground which can be expressed as 1 [notdef]v where [notdef]v is the sky-view

factor (see Eq. 45), and uhv represents the mean horizontal wind speed at the canopy top (Eq. 62), which corresponds to the height zhv (m). The characteristic vertical snow-ake velocity, wn, is set to 0.8 m s1 (Isymov, 1971). Jn is set to 103 m, which is assumed to represent the typical size of the mean forested down wind distance.

For calm conditions and completely vertically falling snowakes, Cn,vp = Cc. For any existing wind, snow could

be intercepted by the surrounding trees so that high wind speed increases interception efciency. Generally for open boreal conifer canopies, Cn,vc < Cn,vp < 1. Under normal wind speed conditions (i.e., wind speeds larger than 1 m s1),

Cn,vc (and Cn,vp) values are usually close to unity.

The maximum allowed canopy snow load, Wrn , is a function of the maximum snow load per unit branch area, Sn,v (kg m2), and the leaf area index:

Wrn = Sn,v LAI (76)

where Sn,v is dened as

Sn,v = Sn,v

0.27 +46 n,v

. (77)

Based on measurements, Schmidt and Gluns (1991) estimated average values for Sn,v of 6.6 for pine and 5.9 kg m2 for spruce trees. Because the average value for this parameter only varies by about 10 % across these two fairly common tree species, and ECOCLIMAP does not currently make a clear distinction between these two forest classes, we currently use 6.3 as the default value for all forest classes. n,v is the canopy snow density (kg m3) dened by the relationship

n,v = 67.92 + 51.25exp

(Tc Tf)/2.59[bracketrightbig]

(Tc Tcmax), (78)

where Tc is the canopy air temperature and Tcmax is the temperature corresponding to the maximum snow density. Assuming a maximum snow density of 750 kg m3 and solving

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 857

Eq. (78) for canopy temperature yields Tcmax = 279.854 K.

This gives values of Sn,v in the range 46 kg m2.

The water vapor ux between the intercepted canopy snow and the canopy air Ern (Eq. C6), includes the evaporative efciency pnv. This effect was rst described by Nakai et al. (1999). In the ISBA-MEB parameterization, the formulation is slightly modied so that it approaches zero when there is no intercepted snow load:

pnv =

0.89Snv0.3

1 + exp[4.7(Snv 0.45)]

, (79)

where Snv is the ratio of snow-covered area on the canopy to the total canopy area

Snv =

Wrn

Wrn (0 Snv 1). (80)

A numerical test is performed to determine if the canopy snow becomes less than zero within one time step due to sublimation. If this is true, then the required mass is removed from the underlying snowpack so that the intercepted snow becomes exactly zero during the time step to ensure a high degree of mass conservation. Note that this adjustment is generally negligible.

The intercepted snow unloading, due to processes such as wind and branch bending, has to be estimated. Hedstrom and Pomeroy (1998) suggested an experimentally veried exponential decay in load over time t, which is used in the parameterization:

Un,v = In,v,0 exp(UnLt) = In,v,0 cnL, (81) where UnL is an unloading rate coefcient (s1) and cnL the dimensionless unloading coefcient. Hedstrom and Pomeroy (1998) found that cnL = 0.678 was a good approximation

that, with a time step of 15 min, gives UnL = 4.498 [notdef]

106 s1. A tuned value for the RCA-LSM from the Snow Model Intercomparison Project phase 2 (SnowMIP2) experiments (Rutter et al., 2009) is UnL = 3.4254 [notdef] 106 s1,

which has been adopted for MEB for now. All unloaded snow is assumed to fall to the ground where it is added to the snow storage on forest ground. Further, corrections to compensate for changes in the original LSM due to this new parameterization have been made for heat capacity, latent heat of vaporization, evapotranspiration, snow storages, and uxes of latent heat.

Finally, canopy snow will partly melt if the temperature rises above the melting point and become intercepted water, where the intercepted (liquid and frozen) water phase change is simply proportional to the temperature:

[Phi1]v =

Ci Wrn

, (84)

where [notdef]v is a view factor indicating how much of the precipitation that should fall directly to the ground (see Eq. 45). The overstory canopy drip rate Drv is dened simply as the value of water in the reservoir which exceeds the maximum holding capacity:

Drv = max 0, Wrv Wrv,max

[parenrightbig]

/[Delta1]t, (85)

where the maximum liquid water-holding capacity is dened simply as

Wrv,max = cwrv LAI. (86)

Generally speaking, cwrv = 0.2 (Dickinson, 1984), al

though it can be modied slightly for certain vegetation cover. Note that Eq. (68) is rst evaluated with Drv = 0, and

then the canopy drip is computed as a residual. Thus, the nal water amount is corrected by removing the canopy drip or throughfall. This water can then become a liquid water source for the soil and the ground-based snowpack.

The fraction of the vegetation covered with water is dened as

v = (1 !rv)

[parenleftbigg]

Wr Wr,max

Ci Snv W rn

Lf [Phi1] (Tf Tv), (82)

where [Phi1]v < 0 signies melting. Tf represents the melting point temperature (273.15 K) and the characteristic phase

change timescale is [Phi1] (s). If it is assumed that the available heating during the time step for phase change is proportional to canopy biomass via the LAI then Eq. (82) can be written (for both melt and refreezing) as

[Phi1]v = Snv k[Phi1]v (Tf Tv). (83)

Note that if energy is available for melting, the phase change rate is limited by the amount of intercepted snow, and likewise freezing is limited by the amount of intercepted liquid water. The melting of intercepted snow within the canopy can be quite complex, thus currently the simple approach in Eq. (83) adopted herein. The phase change coefcient was tuned to a value of k[Phi1]v = 5.56 [notdef] 106 kg m2 s1 K1 for

the SNOWMIP2 experiments with the RCA-LSM. Currently, this value is the default for ISBA-MEB.

2.8.2 Canopy rain interception

The rain intercepted by the vegetation is available for potential evaporation, which means that it has a strong inuence on the uxes of heat and consequently also on the surface temperature. The rate of change of intercepted water on vegetation canopy is described by Eq. (68). The rate that water is intercepted by the overstory (which is not buried by the ground-based snow) is dened as

Prv = Pr (1 [notdef]v) 1 pngp n

[parenrightbig]

Lf [Phi1] (Tf Tv) =

2/3

. (87)

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

!rvWr

+ (1 + arv LAI)Wr,max arv Wr

feature that permits a coupling transition of the snowpack from the canopy air to the free atmosphere as a function of snow depth and canopy height using a fully implicit numerical scheme. MEB has been developed in order to meet the criteria associated with computational efciency, high coding standards (especially in terms of modularity), conservation (of mass, energy, and momentum), numerical stability for large (time step) scale applications, and state-of-the-art representation of the key land surface processes required for current hydrological and meteorological modeling research and operational applications at Mto-France and within the international community as a part of the HIRLAM consortium. This includes regional scale real-time hindcast hydro-meteorological modeling, coupling within both research and operational non-hydrostatic models, regional climate models, and a global climate model, not to mention being used for ongoing ofine land surface reanalysis projects and fundamental research applications.

The simple composite soilvegetation surface energy budget approach of ISBA has proven its ability to provide solid scientic results and realistic boundary conditions for hydrological and meteorological models since its creation over 2 decades ago. However, owing to the ever increasing demands of the user community, it was decided to improve the representation of the vegetation processes as a priority.The key motivation of the MEB development was to move away from the composite scheme in order to address certain known issues (such as excessive bare-soil evaporation in forested areas, the neglect of canopy snow interception processes), to improve consistency in terms of the representation of the carbon cycle (by modeling explicit vegetation energy and carbon exchanges), to add new key explicit processes (forest litter, the gradual covering of vegetation by ground-based snow cover), and to open the door to potential improvements in land data assimilation (by representing distinct surface temperatures for soil and vegetation). Finally, note that while some LSMs intended for GCMs now use multiple-vegetation layers, a single bulk vegetation layer is currently used in MEB since it has been considered as a reasonable rst increase in complexity level from the composite soilvegetation scheme. However, MEB has been designed such that the addition of more canopy layers could be added if deemed necessary in the future.

This is part one of two companion papers describing the model formulation of ISBA-MEB. Part two describes the model evaluation at the local scale for several contrasting well-instrumented sites in France, and for over 42 sites encompassing a wide range of climate conditions for several different forest classes over multiple annual cycles (Napoly et al., 2016). This two-part series of papers will be followed by a series of papers in upcoming years that will present the evaluation and analysis of ISBA-MEB with a specic focus (coupling with snow processes, regional to global scale hydrology, and nally fully coupled runs in a climate model).

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

858 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Delire et al. (1997) used the rst term on the RHS of Eq. (87) for relatively low vegetation (Deardorff, 1978) and the second term for tall vegetation (Manzi and Planton, 1994). Currently in ISBA, a weighting function is used which introduces the vegetation height dependence using the roughness length as a proxy from

!rv = 2z0v 1 (0 !rv 1), (88)

where the current value for the dimensionless coefcient is arv = 2.

2.8.3 Halstead coefcient

In the case of wet vegetation, the total plant evapotranspiration is partitioned between the evaporation of intercepted water, and transpiration via stomata by the Halstead coefcient.In MEB, two such coefcients are used for the non-snow-buried and buried parts of the vegetation canopy, hvg and hvn (Eqs. 29a and 29b, respectively). In MEB, the general form of the Halstead coefcient, as dened in Noilhan and Planton (1989), is modied by introducing the factor kv to take into account the fact that saturated vegetation can transpire, i.e., when v = 1 (Bringfelt et al., 2001). Thus, for MEB we de

ne = kv v. The intercepted water forms full spheres just

touching the vegetation surface when kv = 0, which allows

for full transpiration from the whole leaf surface. In contrast, kv = 1 would represent a situation where a water lm covers

the vegetation completely and no transpiration is allowed. To adhere to the interception model as described above, where the intercepted water exists as droplets, we set the value of kv to 0.25. Note that in the case of condensation, i.e., E < 0, hv = 1.

Without a limitation of hvg and hvn, the evaporative demand could exceed the available intercepted water during a time step, especially for the canopy vegetation which experiences a relatively low aerodynamic resistance. To avoid such a situation, a maximum value of the Halstead coefcient is imposed by calculating a maximum value of the v. See Appendix F for details.

3 Conclusions

This paper presents the description of a new multi-energy balance (MEB) scheme for representing tall vegetation in the ISBA land surface model component of the SURFEX land atmosphere coupling and driving platform. This effort is part of the ongoing effort within the international scientic community to continually improve the representation of land surface processes for hydrological and meteorological research and applications.

MEB consists in a fully implicit numerical coupling between a multi-layer physically based snowpack model, a variable-layer soil scheme, an explicit litter layer, a bulk vegetation scheme, and the atmosphere. It also includes a

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 859

4 Code availability

The MEB code is a part of the ISBA LSM and is available as open source via the surface modeling platform called SURFEX, which can be downloaded at http://www.cnrm-game-meteo.fr/surfex/

Web End =http://www. http://www.cnrm-game-meteo.fr/surfex/

Web End =cnrm-game-meteo.fr/surfex/ http://www.cnrm-game-meteo.fr/surfex/

Web End = . SURFEX is updated at a relatively low frequency (every 3 to 6 months) and the developments presented in this paper are available starting with SURFEX version 8.0. If more frequent updates are needed, or if what is required is not in Open-SURFEX (DrHOOK, FA/LFI formats, GAUSSIAN grid), you are invited to follow the procedure to get a SVN account and to access real-time modications of the code (see the instructions at the previous link).

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

860 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Appendix A: Thermodynamic coupling variable

If potential temperature is used as the thermodynamic variable in the coupled model diffusion scheme, then the thermodynamic variable T (J kg1; see Eqs. 1014) coefcients are

dened as

Bx =0 (x = v,g,n,c,a), (A1)

Ax =Cp/[Pi1]s (x = v,g,n,c), (A2)

Aa =Cp/[Pi1]a, (A3)

where [Pi1] is the non-dimensional Exner function and Cp is

the heat capacity of dry air (J kg1 K1). If the atmospheric variable being diffused is dry static energy then

Bx =0 (x = v,g,n,c), (A4)

Ba =g za, (A5)

Ax =Cp (x = v,g,n,c,a), (A6) where za is the height (m) of the simulated or observed over-lying atmospheric temperature, Ta and g is the gravitational constant. The choice of the atmospheric thermodynamic variable is transparent to ISBA-MEB (it is made within the surfaceatmosphere coupler). The default (in ofine mode and in inline mode with certain atmospheric models) is using Eqs. (A1)(A3). Note that the method can be extended to use the actual air heat capacity (including water vapor) if a linearization of the heat capacity is used.

Appendix B: Latent heat normalization factor

The L is a normalization factor (Lv L Ls), which could

be determined in a number of ways. This coefcient ensures conservation of mass between the different surfaces and the atmosphere.

One possible method is to diagnose it by inverting the equation for LEc (multiplying Eq. 20 by L thereby eliminating it from the RHS of this equation, and then solving for L), but the resulting equation is difcult to apply since the terms can be either positive or negative, and division by a small number is possible. Here, a more smooth (in time) function is proposed, which accounts for each of the surfaces weighted by its respective fraction:

L =

aLs Ls + aLv Lv aLs + aLv

Appendix C: Turbulent ux expressions

The turbulent uxes of heat and water vapor can be further decomposed into different components, which are required for computing different diagnostics and coupling with the water budgets. They are presented herein.

C1 Sensible heat ux

It is convenient to split Hn into two components since one governs the coupling between the canopy air space and the snow surface, while the other modulates the exchanges with the overlying atmosphere (as the canopy layer becomes buried).

The ground-based snowpack heat ux, Hn (Eq. 12), can be split into a part that modulates the heat exchange with the canopy air space, Hnc and the other part which controls

the exchanges directly with the overlying atmosphere, Hna,

dened as

Hnc =a

(Tn Tc) Ranc

, (C1)

. (C2)

Tc is diagnosed by imposing conservation of the heat uxes between the surface and the canopy air (as described in Appendix I). Using the denition in Eq. (C2), the total sensible heat ux exchange with the atmosphere (Eq. 14) can also be written in more compact form as

H = a[bracketleftbig][parenleftBigg]1

png pn

Hna =a

(Tn Ta) Rana

[parenrightbig]

Hc + png pn Hna

[bracketrightbig]

. (C3)

C2 Water vapor ux

The various water vapor ux terms must be broken into different components for use within the different water balance equations for the vegetation, soil, and snowpack. Using the denitions in Eqs. (27)(29b), the components of the canopy evapotranspiration Ev can be expressed as

Etr = a

Lv L

(qsatv qc)

[bracketleftbigg]

png (1 pn ) Ravnc + Rsn +

1 png Ravgc + Rs

[bracketrightbigg]

(1 pnv) (1 ),

(C4)

, (B1)

where

aLv=

f (1 pnv) +

Er = a

Lv L

1 png

[parenrightbig][parenleftBigg]1

pgf

1 png[parenrightbig]

pgf + png

[bracketrightbig][parenleftBigg]1

[parenrightbig][bracketrightbig][parenleftBigg]

1 pngpn

, (B2)

(qsatv qc)

png (1 pn )

Ravnc +1 png Ravgc

aLs=

f pnv +

pngpn

[parenrightbig]

[bracketrightbigg]

(1 pnv) , (C5)

Ern = a

Ls L

(qsatv qc)

png (1 pn )

Ravnc +1 pngRavgc[bracketrightbigg]

pnv. (C6)

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

+ pngpn . (B3)

In the limit as the snow totally buries the canopy vegetation, L ! Ls. In contrast, for snow and surface ice-free condi

tions, L = Lv.

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 861

The complex resistances (bracketed terms in Eqs. C4C6) arise owing to the inclusion of the effects of burying the snow canopy by the ground-based snowpack. If the ground-based snowpack is not sufciently deep to bury any of the canopy (pn = 0), then the bracketed term in Eq. (C4) simplies to

1/ Ravgc + Rs

[parenrightbig]

(note that Ravgc = Ravc when pn = 0

from Eq. 16), and likewise the bracketed terms in Eqs. (C5) (C6) simplify to 1/Ravgc. Finally, the partitioning between

the vapor uxes from intercepted snow and the snow-free canopy reservoir and transpiration is done using pnv, which represents the fraction of the snow interception reservoir that is lled (see Eq. 79).

Using the denitions of qg from Eq. (22) together with those for the humidity factors, hsg and ha (Eqs. 23 and 24, respectively) and the soil coefcient, g (Eq. 25), the bare-soil evaporation Eg components can be expressed as

Egl= a [parenleftbigg]

Lv

L

[parenrightbigg][parenleftBigg][parenleftBigg]h

ug qsatg qc

[parenrightbig][parenleftbigg]

gcor

Rag + Rg

[parenrightbigg][parenleftBigg][parenleftBigg]1

pgf

, (C7)

[parenrightbig][parenleftbigg]

gfcor

Rag + Rgf

[parenrightbigg]

Figure E1. Simple schematic for longwave radiation transfer for one reection and up to three emitting surfaces (in addition to the down-welling atmospheric ux). Hollow arrows indicate uxes after one reection.

where the Richardson number Ri is dened in Eq. (65). The coefcients are dened as

Cv,N =ln

1 + z

exp[parenleftbigg]

k

pCDN [parenrightbigg] 1

pgf, (C8)

where Eg = Egl + Egf. The delta function, gfcor, is a nu

merical correction term, which is required owing to the linearization of qsatg and is unity unless both hugf qsatg < qc and qsatg > qc, in which case it is set to zero. Note that the ground resistances, Rg and Rgf, are set to zero if the forest litter option is active (the default for forests).

The ground-based snowpack sublimation, En (Eq. 19), can be partitioned into a vapor exchange with the canopy air space, Enc and the overlying atmosphere, Ena, as

Enc = a

Ls L

Egf = a

LsL

[parenrightbigg][parenleftBigg][parenleftBigg]hugf qsatg qc

[bracketrightbigg]

, (D2)

Cv,S = z

[parenleftbigg]

k pCDN

k pCD [parenrightbigg]

, (D3)

exp[parenleftbigg]

k

pCDN

k

pCD [parenrightbigg] 1[bracketrightbigg]

, (D4)

where the drag coefcient CD and the drag coefcient for neutral conditions CDN are computed between the canopy air space and the free atmosphere above using the standard ISBA surface-layer transfer functions (Noilhan and Mahfouf, 1996).

Appendix E: Longwave radiative ux expressions

The complete expression for the vegetation canopy net long-wave radiation with an innite number of reections can be expressed as a series expansion (e.g., Braud, 2000) as a function of the temperatures of the emitting surfaces (Tv,

Tg,1, Tn,1), their respective emissivities ([epsilon1]v, [epsilon1]g and [epsilon1]n) and the canopy longwave absorption function, LW (Eq. 45). The MEB expressions are derived by explicitly expanding the series and assuming one reection from each emitting source, which is a good approximation since emissivities are generally close to unity (uxes from a single reection are proportional to 1 [epsilon1]x where x represents g, v, or n, and [epsilon1] is close

to unity for most natural surfaces).

Snow is considered to be intercepted by the vegetation canopy and to accumulate on the ground below. The corresponding schematic of the radiative transfer is shown in Fig. E1. The canopy-intercepted snow is treated using a

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

Cv,U = ln

1 + z

[parenrightbigg] [parenleftbigg]

qsatin qc Ranc

, (C9)

Ena = a

Ls L

[parenrightbigg] [parenleftbigg]

qsatin qa Rana

. (C10)

The corresponding latent heat uxes can be determined by simply multiplying Eqs. (C4)(C8) by L. Finally, using the denition in Eq. (C10), the total vapor exchange with the atmosphere (Eq. 21) can also be written in more compact form as

E = a[bracketleftbig][parenleftBigg]1

png pn

[parenrightbig]

Ec + png pn Ena

[bracketrightbig]

. (C11)

Appendix D: Canopy-top wind stability factor

The expressions for the stability factor fhv (Eq. 62), which is used to compute the wind at the top of the vegetation canopy uhv, are taken from Samuelsson et al. (2006, 2011). They are dened as

fhv = Cv,N + Cv,S

[parenrightbig][radicalbig]C

D /k (Ri > 0), (D1a)

= Cv,N + Cv,U

[parenrightbig][radicalbig]C

D /k (Ri 0), (D1b)

862 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

composite approach so that the canopy temperature Tv represents the effective temperature of the canopy-intercepted snow composite. The canopy emissivity is therefore simply dened as

[epsilon1]v = (1 pnv)[epsilon1]v + pnv [epsilon1]n. (E1)

In order to facilitate the use of a distinct multi-layer snow-process scheme, we split the uxes between those interacting with the snowpack and the snow-free ground. The expressions for the snow-free surface are

Ag = LW # 1 png [parenrightbig]

dened as

[prime]LW = [parenleftBig]

1 p[prime]ng[parenrightBig]

, (E2a)

Bg = Ag LW (1 [epsilon1]v), (E2b)

Cg = Ag (1 LW), (E2c)

Dg = Cg 1 [epsilon1]g

[parenrightbig]

, (E2d)

Eg = Dg 1 [prime]LW

, (E2e)

Fg = [prime]LW [epsilon1]v T 4v 1 png

[parenrightbig]

, (E2f)

Gg = Fg 1 [epsilon1]g

[parenrightbig]

LW + p[prime]ng fLW. (E4)

The factor, fLW, over the understory snow-covered fraction of the grid box is modeled quite simply from Eq. (46). The net longwave radiation for the understory, snowpack, and vegetation canopy are therefore dened, respectively, as

LWnetg =Cg + Fg + Jg + Jn Dg Gg Ig, (E5a) LWnetn =Cn + Fn + Kn + Kg Dn Gn In, (E5b)

LWnetv =Ag + Dg + Gg + Ig + An + Dn + Gn

+ In Bg Cg Eg Hg 2Fg

Jg Lg Kg Bn Cn En

Hn 2Fn Jn Ln Kn, (E5c)

where the upwelling longwave radiation is computed from

LW "= LW # LWnetg LWnetn LWnetv. (E6)

The inclusion of the snow-buried canopy fraction in Eqs. (E2m) and (E3m) causes all of the vegetation transmission and below canopy uxes to vanish as png and pn ! 0

so that the only longwave radiative exchanges occur between

the atmosphere and the snowpack in this limit.

E1 Net longwave radiation ux derivatives

The rst-order derivatives of the net longwave radiation terms are needed in order to solve the system of linearized surface energy budget equations (Eqs. I1I3). The Taylor series expansion (neglecting higher-order terms) is expressed as

LW+neti =LWneti +

Nseb

[parenrightbig]

, (E2g)

Hg = Gg 1 [prime]LW

[parenrightbig]

, (E2h)

Ig = [epsilon1]g T 4g 1 png

Jg = Ig [prime]LW (1 [epsilon1]v) [parenleftBig]

[parenrightbig]

, (E2i)

1 p[prime]ng[parenrightBig]

, (E2j)

Kg = Ig [prime]LW (1 [epsilon1]v) p[prime]ng, (E2k)

Lg = Ig 1 [prime]LW

[parenrightbig]

, (E2l) p[prime]ng = png (1 pn ), (E2m)

and the equations for the snow-covered understory fraction are

An = LW # png, (E3a)

Bn = An fLW (1 [epsilon1]v), (E3b)

Cn = An (1 fLW), (E3c)

Dn = Cn (1 [epsilon1]n), (E3d)

En = Dn 1 [prime]LW

[parenrightbig]

, (E3e)

Fn = f LW [epsilon1]v T 4v png, (E3f)

Gn = Fn (1 [epsilon1]n), (E3g)

Hn = Gn 1 [prime]LW

[parenrightbig]

Xj=1 @Lneti @Tj

T +j Tj [parenrightBig]

(i = 1,Nseb), (E7)

where Nseb represents the number of surface energy budgets, and i and j represent the indexes for each energy budget. The superscript + represents the variable at time t +[Delta1]t, while by

default, no superscript represents the value at time t. Equation (E7) therefore results in a Nseb [notdef] Nseb Jacobian matrix

(3 [notdef] 3 for MEB). The matrix coefcients are expressed as @LWnetv

@Tv =

@Gg

@Tv

, (E3h)

In = [epsilon1]n T 4n png, (E3i)

Jn = In [prime]LW (1 [epsilon1]v) [parenleftBig]

@Hg

@Tv 2

@Fg

@Tv +

@Gn

@Tv

1 p[prime][prime]ng[parenrightBig]

@Hn

@Tv

, (E3j)

Kn = In [prime]LW (1 [epsilon1]v) p[prime][prime]ng, (E3k)

Ln = In 1 [prime]LW

[parenrightbig]

2

@Fn

, (E3l)

p[prime][prime]ng = png + pn 1 png

[parenrightbig]

, (E3m)

where the different terms are indicated in Fig. E1. In MEB, the ground-based snowpack depth can increase to the point that it buries the canopy; thus, for both the snow-covered and snow-free understory fractions a modied snow fraction is

@Tv , (E8a) @LWnetv

@Tg =

@Ig @Tg

@Jg @Tg

@Kg

@Tg

@Lg

@Tg , (E8b)

@LWnetv

@In @Tn

@Jn @Tn

@Kn

@Tn

@Ln

@Tn , (E8c)

@LWnetg

@Tn =

@Tv =

@Fg

@Tv

@Gg

@Tv , (E8d)

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 863

@LWnetg

@Tg =

@Jg @Tg

@Ig@Tg , (E8e)

@LWnetg

(1 [notdef]v)

1 pngp n

[parenrightbig]P

r

+ (Wrv/[Delta1]t)

[bracketrightbig](L/L

v)

a (1 pnv)kv

[braceleftBig][bracketleftbig]

png (1 p n)/Ravnc

[bracketrightbig]

+

[bracketleftbig][parenleftBigg]1

png

[parenrightbig]/R

avgc

[bracketrightbig][bracerightBig]

(qsatv qc)

.

@Tn =

@Jn@Tn , (E8f)

@LWnetn

Equation (F2) is an approximation since all of the variables on the RHS use conditions from the start of the time step; however, this method has proven to greatly reduce the risk for occasional numerical artifacts (jumps) and the associated need for mass corrections (if net losses in mass exceed the updated test value for interception storage).

Appendix G: Energy and mass conservation

G1 Energy conservation

The soil and snowpack prognostic temperature equations can be written in ux form for k = 1,Ng soil layers and k = 1,Nn

snow layers as

Cg,k @Tg,k

@t = Gg,k1 Gg,k + Lf [Phi1]g,k, (G1)

Cn,k @Tn,k

@t = Gn,k1 Gn,k + Lf [Phi1]n,k + n,k1

n,k + SWnet,n n,k1 n,k

[parenrightbig]

@Tv =

@Fn

@Tv

@Gn

@Tv , (E8g)

@LWnetn

@Kg

@Tg =

@Tg , (E8h) @LWnetn

@Tn =

@Jn @Tn

@In@Tn . (E8i)

Using Eq. (E5) to evaluate the derivatives we have

@LWnetv

@Tv =

4Tv Gg Hg 2Fg + Gn Hn 2Fn

[parenrightbig]

, (E9a)

@LWnetv

@Tg =

4Tg Ig Jg Kg Lg

[parenrightbig]

, (E9b)

@LWnetv

@Tn =

4Tn (In Jn Kn Ln), (E9c)

@LWnetg

@Tv =

4Tv Fg Gg

[parenrightbig]

, (E9d)

@LWnetg

. (G2)

The total energy balance of the vegetation canopysoil snowpack system is conserved at each time step [Delta1]t and can be obtained by summing the discrete time forms of Eqs. (4), (G1), and (G2) for the vegetation and all soil and snow layers, respectively, yielding

Cv[Delta1]Tv +

Ng

Pk=1 Cg,k

[Delta1]Tg,k + pngNn

Pk=1 Cn,k

[Delta1]Tn,k =

[Delta1]t

@Tg =

4Tg Jg Ig

[parenrightbig]

, (E9e)

@LWnetg

@Tn =

4Tn Jn, (E9f)

@LWnetn

@Tv =

4Tv (Fn Gn), (E9g)

@LWnetn

@Tg =

4Tg Kg, (E9h)

@LWnetn

4Tn (Jn In), (E9i)

and therefore from a coding perspective, the computation of the derivatives is trivial (using already computed quantities).

Appendix F: Halstead coefcient maximum

A maximum Halstead coefcient is imposed by estimating which value of v that is needed to just evaporate any existing intercepted water Wrv given the conditions at the beginning of the time step. Assuming that phase changes are small, and neglecting canopy drip and any condensation from transpiration, the time-differenced prognostic equation for intercepted water on canopy vegetation (Eq. 68) can be approximated as

Wrv+ Wrv

[Delta1]t = (1 [notdef]v)(1 pngp n)Pr Er. (F1) Assuming that all existing water evaporates in one time step(i.e., W+rv = 0), and substituting the full expression for Er

(Eq. C5) into Eq. (F1), the maximum value of v can be determined as

v,max = (F2)

@Tn =

[bracketleftbigg][parenleftBigg]1

png

[parenrightbig]

Gg,0 + +png Ggn + n,NnSWnet,n + Gn,0

[parenrightbig]

,

Ng

+ Rnv Hv LEv + Lf [Phi1]v +

Pk=1 [Phi1]g,k

+png

Nn

Pk=1 [Phi1]n,k

[parenrightBigg][bracketrightbigg]

,

(G3)

where [Delta1]Tx = Tx(t + [Delta1]t) Tx(t). Note that Eq. (G2) must

rst be multiplied by png in order to make it patch or grid cell relative when it is combined with the soil and vegetation budget equations. The surface boundary conditions for Eqs. (4) and (6) are, respectively,

Gg,0 = Rng Hg LEg, (G4)

Gn,0 = Rnn Hn LEn, (G5)

n,0 = 1, (G6) n,0 = 0. (G7)

Equation (G6) signies that the net shortwave radiation at the surface enters the snowpack, and Eq. (G7) represents the fact that energy changes owing to the time-evolving snow grid

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

864 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

can only arise in the surface layer owing to exchanges with the sub-surface layer. Snowfall is assumed to have the same temperature as the snowpack; thus, a corresponding cooling/heating term does not appear in Eq. (G5), although the corresponding mass increase must appear in the snow water budget equation (see Sect. 2.7).

The lower boundary conditions for Eqs. (G1) and (G2) are, respectively,

Gg,Ng = 0, (G8) n,Nn = 0. (G9)

The appearance of the same discrete form for [Phi1] in both the energy and mass budget equations ensures enthalpy conservation. Owing to Eqs. (G7) and (G9), the total effective heating of the snowpack owing to grid adjustments is

DNn

[integraldisplay]0

where Eq. (G11) has been multiplied by png to make it patch or grid box relative (as was done for energy conservation in Sect. G1). R0 can simply be a diagnostic or coupled with a river routing scheme (Habets et al., 2008; Decharme et al., 2012; Getirana et al., 2015). The soil water lower boundary condition Fg,Ngw represents the base ow or drainage leaving the lowest hydrological layer, which can then be transferred as input to a river routing scheme (see references above) or to a ground water scheme. In such instances, it can be negative if an option to permit a ground water inow is activated (Vergnes et al., 2014). The soil liquid water and equivalent frozen water equivalent volumetric water content extend down to layer Ngw, where Ngw Ng. Note that the verti

cal soil water transfer or evolution is not computed below zg k = Ngw

, whereas heat transfer can be. In order to compute the thermal properties for deep soil temperature (thermal conductivity and heat capacity for example), soil moisture estimates are needed: values from the soil are extrapolated downward assuming hydrostatic equilibrium A detailed description of the soil model is given by Decharme et al. (2011) and Decharme et al. (2013).

Note that Eq. (G11) is snow relative; therefore, this equation must be multiplied by the ground-based snow fraction png to be grid box relative for coupling with the soil and vegetation water storage terms. The lower boundary condition for liquid water ow Fnl,Nn is dened as the liquid water exceeding the lowest maximum snow-layer liquid water-holding capacity. nl represents the internal mass changes of a snowpack layer when the vertical grid is reset. When integrated over the entire snowpack depth, this term vanishes (analogous to Eq. (G10) for the snowpack temperature equation). See Boone and Etchevers (2001) and Decharme et al. (2016) for details on the snow model processes.

The equations describing ooding are not described in detail here as this parameterization is independent of MEB, and it is described in detail by Decharme et al. (2012). The coupling of MEB with the interactive ooding scheme will be the subject of a future paper.

Appendix H: Implicit numerical coupling with the atmosphere

The landatmosphere coupling is accomplished through the atmospheric model vertical diffusion (heat, mass, momentum, chemical species, aerosols, etc.) and radiative schemes.Owing to the potential for relatively large diffusivity, especially in the lower atmosphere near the surface, fairly strict time step constraints must be applied. In this section, a fully implicit time scheme (with an option for explicit coupling) is described. There are two reasons for using this approach;(i) an implicit coupling is more numerically stable not only for time steps typical of GCM applications but also for some NWP models, and (ii) the methodology permits code modularity in that the land surface model routines can be indepen-

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

n dDn = 0, (G10)

where DNn represents the total snow depth. Thus, this term only represents a contribution from contiguous snow layers, not from a source external to the snowpack. The energy storage of the snowsoilvegetation system is balanced by the net surface radiative and turbulent uxes and internal phase changes (solid and liquid phases of water substance).

G2 Mass conservation

The soil and snowpack prognostic mass equations can be written in ux form for k = 2,Ngw soil layers and k = 1,Nn

snow layers as

@Wn,k

@t = Fnl,k1 Fnl,k [Phi1]n,k

+nl,k nl,k1 (k = 2,Nn), (G11)

w [Delta1]zg,1 @wg,k

@t = Fg,k1 Fg,k [Phi1]g,k

F2,k max(0, Etr) k = 2,Ngw [parenrightbig]

, (G12)

w [Delta1]zg,1 @wgf,k

@t = [Phi1]g,k k = 2,Ngw

[parenrightbig]

. (G13)

The total grid box water budget at each time step is obtained by summing the budget equations for the surface layers (Eqs. 6872) together with those for the sub-surface layers (Eqs. G11G13) to have

[Delta1]Wr + [Delta1]Wrn + png

Nn

Xk=1[Delta1]Wn,k + w

Ngw

Xk=1[Delta1]zg,k wgk + wgfk

[parenrightbig]

= [Delta1]t

hPr + Ps R0 Fg,Ngw 1 png

[parenrightbig]

Eg Ev

Ng

png En [Phi1]v

Xk=1[Phi1]g,k pngNn

Xk=1 [Phi1]n,k

i, (G14)

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 865

dent of the atmospheric model code and they can be called using a standard interface, which is the philosophy of SURFEX (Masson et al., 2013). The coupling follows the methodology rst proposed by Polcher et al. (1998), which was further generalized by Best et al. (2004).

The atmospheric turbulence scheme is generally expressed as a second-order diffusion equation in the vertical (which is assumed herein) and it is discretized using the backward difference time scheme. Note that a semi-implicit scheme, such as the CrankNicolson (Crank and Nicolson, 1947), could also be used within this framework. Thus, the equations can be cast as a tri-diagonal matrix. Assuming a xed for zero (the general case) upper boundary condition at the top of the atmosphere, the diffusion equations for the generic variable can be cast as a linear function of the variable in the layer below (Richtmeyer and Morton, 1967) as

+k = B,k + A,k +k+1 (k = 1,Na 1), (H1) where Na represents the number of atmospheric model layers, k = 1 represents the uppermost layer with k increasing

with decreasing height above the surface, and the superscript

+ indicates the value of at time t + [Delta1]t (at the end of the

time step). The coefcients A,k and B,k are computed in a downward sweep within the turbulence scheme and thus consist in atmospheric prognostic variables, diffusivity, heat capacities, and additional source terms from layer k and above evaluated at time level t (Polcher et al., 1998). As shown by Best et al. (2004), the equation for the lowest atmospheric model layer can be expressed using a ux lower boundary condition as

+Na = B,Na + A,Na F +,Na+1, (H2)

where F +,Na+1 is the implicit surface ux from one or mul

tiple surface energy budgets. Technically, only the B,Na and A,Na coefcients are needed by the LSM in order to compute the updated land surface uxes and temperatures, which are fully implicitly coupled with the atmosphere. Once F +,Na+1 has been computed by the LSM, it can be returned

to the atmospheric turbulence scheme, which can then solve for +k from k = Na to k = 1 (i.e., the upward sweep). For

explicit landatmosphere coupling or ofine land-only applications, the coupling coefcients can be set to A,Na = 0

and B,Na = Na in the driving code.

Appendix I: Numerical solution of the surface energy budgets

I1 Discretization of surface energy budgets

The surface energy budget equations (Eqs. 46) are integrated in time using the implicit backward difference scheme. They can be written in discretized form as

Cv

T +v Tv

@LWnetv

@Tg,1

T +g,1 Tg,1[parenrightBig]

[parenrightbig][parenleftBigg]+ T +n,1 Tn,1[parenrightBig] +

SWnetv + LWnetv

+ 'v Av T +v Ac T +c

[parenrightbig]

+ hsv 'v L

qsatv +

@qsatv@Tv T +v Tv

[parenrightbig]

+

@LWnetv

@Tn,1

q+c[bracketrightbigg]

, (I1)

T +g,1 Tg,1[parenrightBig]

[Delta1]t =

Cg,1

@LWnetg @Tv

T +v Tv

@LWnetg

@Tg,1

T +g,1 Tg,1[parenrightBig]

@LWnetg

@Tn,1

[parenrightbig][parenleftBigg]+ T +n,1 Tn,1[parenrightBig] +

SWnetg + LWnetg

+ 'g

+ Ag T +g Ac T +c[parenrightBig]

+ 'g L

hsg

qsatg + @qsatg @Tg

T +g Tg[parenrightBig][bracketrightbigg]

ha q+c

[bracketrightbigg][parenleftBigg][parenleftBigg]1

png

[parenrightbig]

+ png [Lambda1]g,n

T n,Nn T +g,1[parenrightBig]

,

[Lambda1]g,1

T +g,1 T +g,2[parenrightBig]

, (I2)

T +n,1 Tn,1[parenrightBig]

[Delta1]t =

png Cn,1

@LWnetn@Tv T +v Tv

[parenrightbig]

+

@LWnetn

@Tg,1

T +g,1 Tg,1[parenrightBig]

T +n,1 Tn,1[parenrightBig] +

SWnetn + LWnetn

+ (1 pn ) 'nc An T +n Ac T +c

[parenrightbig]

+ pn 'na Bn Ba + An T +n Aa T +a

[parenrightbig]

+ (1 pn )'nc Ls

qsatin +

@qsatin@Tn T +n T +c

+

@LWnetn

@Tn,1

[parenrightbig]

q+c[bracketrightbigg]

+ pn 'na Ls

qsatin +

@qsatin@Tn T +n T +a

[parenrightbig]

q+a[bracketrightbigg]

T +n,1 T +n,2[parenrightBig]

png. (I3)

Note that Eq. (I3) has been multiplied by png since the snow-pack must be made patch relative when solving the coupled equations. The q+satx and longwave radiation terms have been linearized with respect to Tx (the longwave radiation derivatives are given by Eq. E9). The superscript + corresponds

to the values of variables at time t + [Delta1]t, while the absence

of a superscript indicates variables evaluated at time t. Note that we have dened 'x = a/Rax (kg m2 s1) for simplic

ity. The thermodynamic variable, Tx, in the sensible heat ux

terms have been expressed as a function of Tx using Eq. (15).

Several of the Bx terms have canceled out in the sensible heat

ux terms in Eqs. (I1)(I3) since they are dened such that

Bc= Bv = Bg = Bn. Note that compared to Eqs. (4)(6), the

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

[Lambda1]g,1

T +v Tv

[Delta1]t = @LWnetv @Tv

866 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

phase change terms ([Phi1]x) do not appear in Eqs. (I1)(I3). This is because they are evaluated as an adjustment after the energy budget and the uxes have been computed.

In Eq. (I2), T n,Nn represents a test temperature for the lowest snowpack layer. It is rst computed using an implicit calculation of the combined snowsoil layers to get a rst estimate of the snowground heat conduction inter-facial ux when simultaneously solving the surface energy budgets.The nal snow temperature in this layer, T +n,Nn, is computed afterwards within the snow scheme; any difference between the resulting conduction ux and the test ux in Eq. (I2) is added to the soil as a correction at the end of the time step in order to conserve energy. In practice, this correction is generally small, especially since the snow fraction goes to unity very rapidly (i.e., for a fairly thin snowpack when using MEB; see Eq. 1). Thus, in this general case, the difference between the test ux and the nal ux arise only owing to updates to snow properties within the snow scheme during the time step. Since T n,Nn is computed using an implicit solution method for the entire soilsnow continuum, it is also quite numerically stable. The use of a test ux permits a modular coupling between the snow scheme and the soilvegetation parts of ISBA-MEB.

In order to solve Eqs. (I1)(I3) for the three unknown surface energy budget temperatures, T +v, T +g,1, and T +n,1, equations for the six additional unknown surface energy budget temperatures, T +a, T +c, q+a, q+c, T +g,2, and T +n,2, must be dened. They can be expressed as linear equations in terms of

T +v, T +g,1, and T +n,1, and their derivations are presented in the remaining sections of this Appendix.

I2 Atmospheric temperature and specic humidity

The rst step in solving the surface energy budget is to eliminate the lowest atmospheric energy and water vapor variables from the snow surface energy budget equation. They will also be used to diagnose the nal ux exchanges between the canopy air space and overlying atmosphere.

From Eq. (H2), the thermodynamic variable of the lowest atmospheric model variable at time t + [Delta1]t is dened as

T +Na = BT ,Na + AT ,Na H+. (I4)

Note that using Eq. (15), we can rewrite Eq. (I4) in terms of air temperature as

Ta+ = BTa + ATa H+, (I5)

where BTa = BT ,Na Ba

[parenrightbig]

eATa = ATa 'ca Ac 1 png p n

/C, (I7b)

eBTa =

nBTa Ba + ATa [bracketleftBig][parenleftBigg]1

png p n

[parenrightbig]

'ca (Bc Ba)+

png p n 'na (Bc Ba)

[bracketrightBig][bracerightBig]

[parenrightbig]

/C, (I7c)

eCTa = ATa png p n 'na Ac/C. (I7d) In analogous fashion to determining the air temperature,

the specic humidity of the lowest atmospheric model variable at time t + [Delta1]t is dened from Eq. (H2) asq+a = Bq,a + Aq,a E+, (I8)

where again the subscript q, a represents the values of the coefcients A and B for the lowest atmospheric model layer (k = Na). Substitution of Eq. (21) for E in Eq. (I8) and solv

ing for T +a yields

q+a =

eBq,a + eAq,a q+c +

eCq,a q+satin, (I9) where the coefcients are dened as

C = 1 + Aq,a

[bracketleftbig][parenleftBigg]1

png p n

[parenrightbig]

'ca + 'na hsn p n png

[bracketrightbig]

, (I10a)

eAq,a = Aq,a 'ca 1 png p n

[parenrightbig]

/C, (I10b)

eBq,a = Bq,a/C, (I10c)

eCq,a = Aq,a 'na hsn p n png/C. (I10d) I3 Canopy air temperature and specic humidity

In order to close the energy budgets, T +c and q+c must be determined.

Assuming conservation of the heat ux between the different surfaces and the canopy air space, we have

1 pngpn

[parenrightbig]

H+c = png (1 pn ) H+nc

+ 1 png

[parenrightbig]

H+g + H+v, (I11)

which can be expanded as

'ca 1 png p n

[parenrightbig]

[notdef]

Bc+ Ac T +c Ba Aa T +a

=

Ac

h'g

T +g T +c[parenrightBig][parenleftBigg]1

[parenrightbig]

png

[parenrightbig]

i. (I12)

Note that the above conservation equation does not include the part of the snow sensible heat ux, which is in direct contact with the atmosphere (Hna), since it was already ac

counted for in the expression for T +a via Eq. (I5). Eliminating T +a using Eq. I6 and solving for T +c yields

T +c = aTc + bTc T +v + cTc T +g + dTc T +n (I13)

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

/Aa, ATa = AT ,Na/Aa, and Ta is

shorthand for T (k = Na). Substitution of Eq. (14) for H in

Eq. (I5) and solving for T +a yields

T +a =

eBTa + eATa T +c +

eCTa T +n, (I6)

where

+ 'v T +v T +c

'nc T +n T +c

png (1 p n)

C = Aa

n1 + ATa

'ca 1 png p n

[parenrightbig]

+png p n 'na

[bracketrightbig][bracerightBig]

,

(I7a)

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 867

with the coefcients

C ='ca 1 png p n

[parenrightbig][parenleftBigg]A

c

Aa

eAT a

[parenrightbig]

+

Ac

'v + 'g

1 png [parenrightbig]

+'nc png (1 p n)

[bracketrightbig]

, (I14a)

I4 Sub-surface temperatures

The sub-surface conduction heat uxes (Eqs. 4749) can be expressed in compact form as

G+x,k = [Lambda1]x,k [parenleftBig]

T +x,k T +x,k+1[parenrightBig]

aTc =

h'ca 1 png p n

[parenrightbig][parenleftBigg]B

a

Bc + Aa

eBT a

[parenrightbig][bracketrightBig]

/C, (I14b)

bTc =Ac 'v/C, (I14c)

cTc =Ac 'g 1 png

[parenrightbig]

/C, (I14d)

dTc =

Ac 'nc png (1 p n)

+Aa

eCT a 'ca 1 png p n

[parenrightbig][bracketrightbig]

/C. (I14e)

In an analogous fashion for canopy air temperature determination, assuming conservation of the vapor ux between the different surfaces and the canopy air space,

1 pngpn

[parenrightbig]

E+c = png (1 pn ) E+nc

+ 1 png

[parenrightbig]

E+g + E+v, (I15)

which can be expanded using the denitions of the evaporative uxes Ex from Eqs. (17)(I15) together with the denitions of qg from Eq. (22) and q+a from Eq. (I9) as

'ca 1 png p n

[parenrightbig]

, (I19)

where [Lambda1]x,k represents the ratio of the inter-facial thermal conductivity to the thickness between the mid-points of contiguous layers (k and k + 1). Using the methodology de

scribed in Appendix H for the atmospheric diffusion scheme, the soil and snow heat diffusion equation (both using the form of Eq. G1) can be dened in an analogous fashion as

T +g,k = Bg,k + Ag,k T +g,k1 k = 2,Ng

[parenrightbig]

, (I20)

where the coefcients Bg,k and Ag,k are determined during the upward sweep (rst step of the tri-diagonal solution) from the base of the soil to the sub-surface soil and snow layers as described by Richtmeyer and Morton (1967). The resulting coefcients for the soil are dened as

C = Cgk/[Delta1]t

[parenrightbig]

+ [Lambda1]gk1 + [Lambda1]gk 1 Agk+1

[parenrightbig]

, (I21a)

Bgi =

[bracketleftbig][parenleftBigg]C

g k/[Delta1]t

[parenrightbig]

q+c 1

eAq,a

, (I21b)

Agk =[Lambda1]gk1/C. (I21c)

The same form holds for the snow layers. The upward sweep is performed before the evaluation of the energy budget; thus, Eq. (I20) is used to eliminate T +g,2 and T +n,2 from Eqs. (I2) and (I3), respectively. To do this, the sub-surface implicit uxes in Eqs. (5) and (6) can be expressed, respectively, as

G+g,1 =[Lambda1]g,1 [bracketleftBig]

T +g,1 1 Ag,2

hsg q+satg ha q+c[parenrightBig][parenleftBigg]1

[parenrightbig]

[notdef] eBq,a

eCq,a q+satin

=

T gk + [Lambda1]gk Bgk+1

/C

2 k Ng 1

[parenrightbig]

h'g

png

[parenrightbig]

+ 'v hsv q+satv q+c

[parenrightbig]

'nc hsn q+satin q+c

[parenrightbig]

i. (I16)

Owing to the linearization of the qsatx, terms about Tx, Eq. (I16) can be solved for q+c as a function of the surface energy budget temperatures as

q+c = aqc + bqc T +v + cqc T +g + dqc T +n, (I17)

where the coefcients are dened as

C ='ca 1 png pn

[parenrightbig][parenleftBigg]1

png (1 p n)

+ Bg,2

i, (I22a)

G+n,1 =[Lambda1]n,1 [bracketleftBig]

[parenrightbig]

T +n,1 1 An,2

[parenrightbig]

+ Bn,2

i. (I22b)

I5 Surface stresses

Using the same surfaceatmosphere coupling methodology as for temperature and specic humidity, the u wind component in the lowest atmospheric model layer can be expressed as

u+a = Bua + Aua +x. (I23)

The surface u component momentum exchange with the atmosphere is expressed as

+x = u+a

[bracketleftbig][parenleftBigg]1

png pn

eAq,a

[parenrightbig]

+ 'g hN 1 png

[parenrightbig]

+ 'v hsv + 'nc hsn png (1 pn ), (I18a)

aqc =

[braceleftBig][parenleftBigg]1

png pn

[parenrightbig]

'ca

eBq,a + 'v hsv

qsatv

@qsatv@Tv Tv[parenrightbigg]

+ 'g hsg

qsatg

@qsatg@Tg Tg[parenrightbigg][parenleftBigg][parenleftBigg]1

png

[parenrightbig]

+ 'nc hsn

qsati n

@qsati n@Tn Tn[parenrightbigg]

(I18b)

[parenrightbig]

png (1 pn )

o/C, (I18c)

bqc =hsv 'v

@qsatv

@Tv /C, (I18d)

cqc =hsg 'g

@qsatg

@Tg 1 png

'Dca + png pn 'Dna

[bracketrightbig]

, (I24)

where it includes stresses from the snow-buried and non-snow-buried portions of the surface consistent with the uxes of heat and water vapor. For simplicity, we have dened

'Dx = a Va CDx (I25)

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

[parenrightbig]

/C, (I18e)

dqc =hsn 'nc

@qsati n

@Tn png (1 pn )/C. (I18f)

868 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

and CD is the surface drag coefcient, which is dened following Noilhan and Mahfouf (1996). Eliminating +x from

Eq. (I24) using Eq. (I25) gives

u+a =

Bua

1 + Aua 'Dc

, (I26)

where for convenience we have dened the average drag coefcient as

'Dc = 1 png pn

[parenrightbig]

'Dca + png pn 'Dna. (I27)

The net u-momentum ux from the surface to the canopy air space is expressed as

+x =

Bua 'Dc

1 + Aua 'Dc

. (I28)

Finally, the scalar friction velocity can be computed from

u = [parenleftbigg]

'Dc V +a

a

[parenrightbig]

1/2, (I29)

where V +a is the updated wind speed (computed from u+a and v+a). Note that v+a and +y are computed in the same manner, but using Bva from the atmosphere (note that Ava = Aua).

I6 Summary: nal solution of the implicitly coupled equations

The fully implicit solution of the surface and atmospheric variables proceeds for each model time step as follows:

1. Within the atmospheric model, perform the downward sweep of the tri-diagonal matrix within the turbulent diffusion scheme of the atmospheric model to obtain the A,k and B,k coefcients for each diffused variable ( = T , q, u, and v) for each layer of the atmo

sphere (k = 1,Na). Update Aa and Ba, then pass these

values along with the aforementioned coupling coefcients at the lowest atmospheric model layer (i.e., AT,a,

BT,a, Aq,a, Bq,a, Au,a, Bu,a, and Bv,a) to the land surface model. These coefcients are then used to eliminate T +a and q+a from the implicit surface energy budget equations (Eqs. I1I3).

2. Within the land surface model, perform the upward sweep of the tri-diagonal matrix within the soil and snow layers to determine the An,k, Bn,k, Ag,k, and Bg,k, coefcients for the soil and snow layers (from soil-layer Ng to layer 2, and again from soil-layer Ng to layer 2 of the snow scheme). Note that coefcients for layer 1 of the snow and soil schemes are not needed since they correspond to the linearized surface energy budgets (next step).

3. Within the land surface model, the expressions for T +a (Eq. I6), q+a (Eq. I9), T +c (Eq. I13), q+c (Eq. I17), T +g,2

(Eq. I22a)and T +n,2 (Eq. I22b) can now be substituted into the energy budget equations (Eqs. I1I3), which can then be readily solved for T +v, T +g,1, and T +n,1.

4. Within the land surface model, perform back substitution (using T +g,1 as the upper boundary condition) to obtain T +g,k for soil layers k = 2,Ng using Eq. (I20).

5. Within the land surface model, call the explicit snow-process scheme to update the snow scheme temperature, T +n,k, and the snow mass variables for snow layers k = 2,Nn. The implicit snow surface uxes, R+n,n, H+n

and E+n, are used as the upper boundary condition along with the implicit soil temperature, T +g,1, to compute the updated lower snowpack boundary condition (i.e., the snowsoil inter-facial ux, Ggn).

6. Within the land surface model, compute V +a (see Sect. I5). Diagnose T +a, Tc+, q+a and q+c (again, using

the equations mentioned in step 3) in order to compute the updated (implicit) uxes. The updated evapotranspi-ration (Eqs. C4C8) and snowmelt water mass uxes are used within the hydrology schemes to update the different water storage variables for the soil and vegetation canopy (Eqs. 6872).

7. Within the atmospheric model, perform back substitution (using H+, E+, +x and +y as the lower boundary conditions: Eq. H2) to obtain updated proles (or turbulent tendencies, depending on the setup of the atmospheric model) of Tk, qk, uk and vk for atmospheric

layers k = 1,Na. Finally, the updated upwelling short-

wave, SW ", and implicit longwave ux, LW"+ (or

equivalently, the effective emissivity and implicit long-wave radiative temperature, T +rad) are returned to the atmospheric model as lower boundary conditions for the respective radiative schemes.

Alternately, in ofine mode, A,a = 0 and B,a = a in the

driving routine in step 1, and the solution procedure ends at step 6. Finally, if multiple patches and/or tiles are being used within the grid call of interest, the corresponding fractional-area-weighted uxes are passed to the atmospheric model in step 7.

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 869

Acknowledgements. This work was initiated within the international HIRLAM consortium as part of the ongoing collective effort to improve the SURFEX platform for research and operational hydrological and meteorological applications. We wish to make a posthumous acknowledgement of the contribution to this work by Jol Noilhan, who was one of the original supporters of this project and helped initiate this endeavor; his scientic vision was essential at the early stages of this work. We wish to thank other contributors to this development in terms of discussions and evaluation, such asG. Boulet, E. Martin, J.-C. Calvet, P. Le Moigne, C. Canac, and G.

Aouad. The technical support of S. Faroux of the SURFEX team is also greatly appreciated. We also wish to thank E. Lebas for preparing the MEB schematic. Part of this work was supported by a grant from Mto-France.

Edited by: S. ValckeReviewed by: two anonymous referees

References

Anderson, M. C., Norman, J. M., Diak, G. R., Kustas, W. P., and Mecikalski, J. R.: A two-source time-integrated model for estimating surface uxes using thermal infrared remote sensing, Remote Sens. Environ., 60, 195216, 1997.

Andreadis, K. M., Storck, P., and Lettenmaier, D. P.: Modeling snow accumulation and ablation processes in forested environments, Water Resour. Res., 45, w05429, doi:http://dx.doi.org/10.1029/2008WR007042

Web End =10.1029/2008WR007042 http://dx.doi.org/10.1029/2008WR007042

Web End = , 2009.

Bengtsson, L., Andrae, U., Aspelien, T., Batrak, Y., de Rooy,W., Gleeson, E., Hansen-Sass, B., Hortal, M., Ivarsson, K.-I., Lenderink, G., Niemel, S., Nielsen, K. P., Onvlee, J., Rontu,L., Samulesson, P., Subias, A., Tijm, S., Toll, V., Yang, X., and Kltzow, M. .: The HARMONIE-AROME model conguration in the ALADIN-HIRLAM NWP system, Mon. Weather Rev., doi:http://dx.doi.org/10.1175/MWR-D-16-0417.1

Web End =10.1175/MWR-D-16-0417.1 http://dx.doi.org/10.1175/MWR-D-16-0417.1

Web End = , in press, 2017.

Best, M. J., Beljaars, A., Polcher, J., and Viterbo, P.: A proposed structure for coupling tiled surfaces with the planetary boundary layer, J. Hydrometeor., 5, 12711278, 2004.

Best, M. J., Pryor, M., Clark, D. B., Rooney, G. G., Essery, R. L.H., Mnard, C. B., Edwards, J. M., Hendry, M. A., Porson, A., Gedney, N., Mercado, L. M., Sitch, S., Blyth, E., Boucher, O., Cox, P. M., Grimmond, C. S. B., and Harding, R. J.: The Joint UK Land Environment Simulator (JULES), model description Part 1: Energy and water uxes, Geosci. Model Dev., 4, 677699, doi:http://dx.doi.org/10.5194/gmd-4-677-2011

Web End =10.5194/gmd-4-677-2011 http://dx.doi.org/10.5194/gmd-4-677-2011

Web End = , 2011.

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

Web End =10.5194/gmd-7-2193-2014 http://dx.doi.org/10.5194/gmd-7-2193-2014

Web End = , 2014.

Boone, A. and Etchevers, P.: An intercomparison of three snow schemes of varying complexity coupled to the same land-surface model: Local scale evaluation at an Alpine site, J. Hydrometeor., 2, 374394, 2001.

Boone, A., Calvet, J.-C., and Noilhan, J.: Inclusion of a Third Soil Layer in a Land-Surface Scheme using the Force-Restore method, J. Appl. Meteor., 38, 16111630, 1999.

Boone, A., Masson, V., Meyers, T., and Noilhan, J.: The inuence of the inclusion of soil freezing on simulations by a soil-vegetation-atmosphere transfer scheme, J. Appl. Meteor., 9, 15441569, 2000.

Boussetta, S., Balsamo, G., Beljaars, A., Agusti-Panareda, A., Calvet, J.-C., Jacobs, C., van den Hurk, B., Viterbo, P., Lafont, S., Dutra, E., Jarlan, L., Balzarolo, M., Papale, D., and van der Werf,G.: Natural land carbon dioxide exchanges in the ECMWF Integrated Forecasting System: Implementation and ofine validation, J. Geophys. Res., 18, 59235946, doi:http://dx.doi.org/10.1002/jgrd.50488

Web End =10.1002/jgrd.50488 http://dx.doi.org/10.1002/jgrd.50488

Web End = , 2013.

Braud, I.: SiSPAT Users Manual, Model Documentation, September 2000 3.0, LTHE, LTHE, BP 53, 38041 Grenoble Cdex 9, France, 2000.

Braud, I., Bariac, T., Gaudeta, J. P., and Vauclin, M.: SiSPAT-Isotope, a coupled heat, water and stable isotope (HDO and H18

2 O) transport model for bare soil. Part I. Model description and rst verications, J. Hydrol., 309, 277300, 2005.

Bringfelt, B., Risnen, J., Gollvik, S., Lindstrm, G., Graham,L. P., and Ullerstig, A.: The land surface treatment for the Rossby Centre Regional Atmosphere Climate Model version 2, Reports of Meteorology and Climatology 98, SMHI, SE-601 76 Norrkping, Sweden, 2001.

Brun, E., Vionnet, V., Decharme, B., Peings, Y., Valette, R., Boone,A., Karbou, F., and Morin, S.: Simulation of northern Eurasian local snow depth, mass and density using a detailed snowpack model and meteorological reanalyses, J. Hydrometeor., 14, 203 219, 2013.

Cakmur, R. V., Miller, R. L., and Torres, O.: Incorporating the effect of small-scale circulations upon dust emission in an atmospheric general circulation model, J. Geophys. Res., 109, D07201, doi:http://dx.doi.org/10.1029/2003JD004067

Web End =10.1029/2003JD004067 http://dx.doi.org/10.1029/2003JD004067

Web End = , 2004.

Calvet, J. C. and Soussana, J.-F.: Modelling CO2-enrichment effects using an interactive vegetation SVAT scheme, Agr. Forest

Meteorol., 108, 129152, 2001.

Calvet, J.-C., Noilhan, J., Roujean, J.-L., Bessemoulin, P., Cabelguenne, M., Olioso, A., and Wigneron, J.-P.: An interactive vegetation SVAT model tested against data from six contrasting sites, Agr. Forest Meteorol., 92, 7395, 1998.

Carrer, D., Roujean, J.-L., Lafont, S., Calvet, J.-C., Boone, A., Decharme, B., Delire, C., and Gastellu-Etchegorry, J.-P.: A canopy radiative transfer scheme with explicit FAPAR for the interactive vegetation model ISBA-A-gs: Impact on carbon uxes,J. Geophys. Res., 188, 888903, 2013.

Castillo, G., Levis, C. K. S., and Thornton, P.: Evaluation of the New CNDV Option of the Community Land Model: Effects of Dynamic Vegetation and Interactive Nitrogen on CLM4 Means and Variability, J. Climate, 25, 37023714, 2012.

Choudhury, B. J. and Monteith, J. L.: A four-layer model for the heat budget of homogeneous land surfaces, Q. J. Roy. Meteor. Soc., 114, 373398, 1988.

Crank, J. and Nicolson, P.: A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type, Proc. Camb. Philol. Soc., 43, 5067, 1947. Deardorff, J. W.: A parameterization of ground-surface moisture content for use in atmosphere prediction models, J. Appl. Meteor., 16, 11821185, 1977.

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

870 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Deardorff, J. W.: Efcient prediction of ground surface temperature and moisture, with inclusion of a layer of vegetation, J. Geophys.Res., 83, 18891903, 1978.

Decharme, B. and Douville, H.: Introduction of a sub-grid hydrology in the ISBA land surface model, Clim. Dynam., 26, 6578, 2006.

Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Local evaluation of the ISBA soil multilayer diffusion scheme using four pedotransfer functions, J. Geophys. Res., 116, 2012620155, doi:http://dx.doi.org/10.1029/2011JD016002

Web End =10.1029/2011JD016002 http://dx.doi.org/10.1029/2011JD016002

Web End = , 2011.

Decharme, B., Alkama, R., Papa, F., and Faroux, S.: Global ofine evaluation of the ISBA-TRIP ood model, Clim. Dynam., 38, 13891412, doi:http://dx.doi.org/10.1007/s00382-011-1054-9

Web End =10.1007/s00382-011-1054-9 http://dx.doi.org/10.1007/s00382-011-1054-9

Web End = , 2012.

Decharme, B., Martin, E., and Faroux, S.: Reconciling soil thermal and hydrological lower boundary conditions in land surface models, J. Geophys. Res., 118, 78197834, doi:http://dx.doi.org/10.1002/jgrd.50631.

Web End =10.1002/jgrd.50631. http://dx.doi.org/10.1002/jgrd.50631.

Web End = , 2013.

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature proles simulated by the ISBA land surface model, The Cryosphere, 10, 853877, doi:http://dx.doi.org/10.5194/tc-10-853-2016

Web End =10.5194/tc-10-853-2016 http://dx.doi.org/10.5194/tc-10-853-2016

Web End = , 2016.

Delire, C., Calvet, J.-C., Noilhan, J., Wright, I., Manzi, A., and

Nobre, C.: Physical properties of Amazonian soils: A modeling study using the Anglo-Brazilian Amazonian Climate Observation Study data, J. Geophys. Res., 102, 3011930133, doi:http://dx.doi.org/10.1029/97JD01836

Web End =10.1029/97JD01836 http://dx.doi.org/10.1029/97JD01836

Web End = , 1997.de Rosnay, P., Polcher, J., Laval, K., and Sabre, M.: Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over the Indian Peninsula, Geophys. Res.Lett., 30, 1986, doi:http://dx.doi.org/10.1029/2003GL018024

Web End =10.1029/2003GL018024 http://dx.doi.org/10.1029/2003GL018024

Web End = , 2003.Dickinson, R. E.: Modelling evapotranspiration for three-dimensional global models, in: Climate Processes and Climate Sensitivity, Geophys. Monogr., 29, Maurice Ewing, edited by: Hansen, J. E. and Takahashi, T., 1984.

Dickinson, R. E., Shaikh, M., Bryant, R., and Graumlich, L.: Interactive canopies for a climate model, J. Climate, 11, 28232836, 1998.

Douville, H., Royer, J.-F., and Mahfouf, J.-F.: A new snow parameterization for the Mto-France climate model, Part I : Validation in stand-alone experiments, Clim. Dynam., 12, 2135, 1995.Ek, M., Mitchell, K. E., Lin, Y., Rogers, E., Grunmann, P.,

Koren, V., Gayno, G., and Tarpley, J. D.: Implementation of Noah land-surface model advances in the NCEP operational mesoscale Eta model, J. Geophys. Res., 108, 8851, doi:http://dx.doi.org/10.1029/2002JD003296

Web End =10.1029/2002JD003296 http://dx.doi.org/10.1029/2002JD003296

Web End = , 2003.

Erbs, D. G., Klein, S. A., and Dufe, J. A.: Estimation of the diffuse radiation fraction for hourly, daily and monthly average global radiation, Sol. Energy, 28, 293304, 1982.

Faroux, S., Kaptu Tchuent, A. T., Roujean, J.-L., Masson, V., Martin, E., and Le Moigne, P.: ECOCLIMAP-II/Europe: a twofold database of ecosystems and surface parameters at 1 km resolution based on satellite information for use in land surface, meteorological and climate models, Geosci. Model Dev., 6, 563 582, doi:http://dx.doi.org/10.5194/gmd-6-563-2013

Web End =10.5194/gmd-6-563-2013 http://dx.doi.org/10.5194/gmd-6-563-2013

Web End = , 2013.

Galperin, B., Sukoriansky, S., and Anderson, P. S.: On the critical Richardson number in stably stratied turbulence, Atmos. Sci.

Lett., 8, 6569, doi:http://dx.doi.org/10.1002/asl.153

Web End =10.1002/asl.153 http://dx.doi.org/10.1002/asl.153

Web End = , 2007.

Getirana, A. C. V., Boone, A., and Peugeot, C.: Evaluating

LSM-based water budgets over a West African basin assisted with a river routing scheme, J. Hydrometeor., 15, 23312346, doi:http://dx.doi.org/10.1175/JHM-D-14-0012.1

Web End =10.1175/JHM-D-14-0012.1 http://dx.doi.org/10.1175/JHM-D-14-0012.1

Web End = , 2015.

Giard, D. and Bazile, E.: Implementation of a new assimilation scheme for soil and surface variables in a global NWP model, Mon. Weather Rev., 128, 9971015, 2000.

Gibelin, A.-L., Calvet, J. C., Roujean, J.-L., Jarlan, L., and

Los, S. O.: Ability of the land surface model ISBA-A-gs to simulate leaf area index at the global scale: Comparison with satellites products, J. Geophys. Res., 111, d18102, doi:http://dx.doi.org/10.1029/2005JD006691

Web End =10.1029/2005JD006691 http://dx.doi.org/10.1029/2005JD006691

Web End = , 2006.

Gonzalez-Sosa, E., Braud, I., Thony, J.-L., Vauclin, M., Bessemoulin, P., and Calvet, J.-C.: Modelling heat and water exchanges of fallow land covered with plant-residue mulch, Agr.Forest Meteorol., 97, 151169, 1999.

Habets, F., Boone, A., Champeaux, J., Etchevers, P., Leblois, E.,

Ledoux, E., Moigne, P. L., Martin, E., Morel, S., Segui, Q., Rousset-Regimbeau, F., and Viennot, P.: The SAFRAN-ISBAMODCOU hydrometeorological model applied over France, J.Geophys. Res., 113, D06113, doi:http://dx.doi.org/10.1029/2007JD008548

Web End =10.1029/2007JD008548 http://dx.doi.org/10.1029/2007JD008548

Web End = , 2008.Hedstrom, N. R. and Pomeroy, J. W.: Measurements and modelling of snow interception in the boreal forest, Hydrol. Process., 12, 16111625, 1998.

Isymov, N.: An approach to the prediction of snow loads, Tech. rep.,

Faulty of Engineering Science, The Univ. of Western Ontario, London, 442 pp., 1971.

Jarvis, P. G.: The interpretation of variations in leaf water potential and stomatal conductance found in canopies in the eld, Philos.T. R. Soc. Lond., 273, 593610, 1976.

Joetzjer, E., Delire, C., Douville, H., Ciais, P., Decharme, B., Carrer, D., Verbeeck, H., De Weirdt, M., and Bonal, D.: Improving the ISBACC land surface model simulation of water and carbon uxes and stocks over the Amazon forest, Geosci. Model Dev., 8, 17091727, doi:http://dx.doi.org/10.5194/gmd-8-1709-2015

Web End =10.5194/gmd-8-1709-2015 http://dx.doi.org/10.5194/gmd-8-1709-2015

Web End = , 2015.

Lafore, J. P., Stein, J., Asencio, N., Bougeault, P., Ducrocq, V.,

Duron, J., Fischer, C., Hreil, P., Mascart, P., Masson, V., Pinty, J.P., Redelsperger, J. L., Richard, E., and Vil-Guerau de Arellano,J.: The Meso-NH Atmospheric Simulation System. Part I: adiabatic formulation and control simulations, Ann. Geophys., 16, 90109, doi:http://dx.doi.org/10.1007/s00585-997-0090-6

Web End =10.1007/s00585-997-0090-6 http://dx.doi.org/10.1007/s00585-997-0090-6

Web End = , 1998.

Lawrence, D. M., Oleson, K. W., Flanner, M. G., Thornton, P. E.,

Swenson, S. C., Lawrence, P. J., Zeng, X., Yang, Z.-L., Levis,S., Sakaguchi, K., Bonan, G. B., and Slater, A. G.: Parameterization Improvements and Functional and Structural Advances in Version 4 of the Community Land Model, J. Adv. Model. Earth Sys., 3, M03001, doi:http://dx.doi.org/10.1029/2011MS00045

Web End =10.1029/2011MS00045 http://dx.doi.org/10.1029/2011MS00045

Web End = , 2011.

Lind, P., Lindstedt, D., Kjellstrm, E., and Jones, C.: Spatial and

Temporal Characteristics of Summer Precipitation over Central Europe in a Suite of High-Resolution Climate Models, J. Climate, 29, 35013518, 2016.

Louis, J.-F.: A parametric model of vertical eddy uxes in the atmosphere, Bound.-Lay. Meteorol., 17, 187202, 1979.

Mahfouf, J.-F. and Noilhan, J.: Comparative study of various formulations of evaporation from bare soil using in situ data., J. Appl.Meteorol., 9, 351362, 1991.

Mahfouf, J.-F. and Noilhan, J.: Inclusion of gravitational drainage in a land surface scheme based on the force-restore method, J.Appl. Meteor., 35, 987992, 1996.

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/

A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8 871

Manzi, A. O. and Planton, S.: Implementation of the ISBA parameterization scheme for land surface processes in a GCM-an annual cycle experiment, J. Hydrol., 155, 353387, 1994.

Mascart, P., Noilhan, J., and Giordani, H.: A modied parameterization of ux-prole relationships in the surface layer using different roughness length values for heat and momentum, Bound.-Lay. Meteorol., 72, 331344, 1995.

Massman, W.: Foliage distribution in old-growth coniferous tree canopies, Can. J. Forest Res., 12, 1017, 1982.

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse,M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf,J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon,B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and ocean surface platform for coupled or ofine simulation of earth surface variables and uxes, Geosci. Model Dev., 6, 929960, doi:http://dx.doi.org/10.5194/gmd-6-929-2013

Web End =10.5194/gmd-6-929-2013 http://dx.doi.org/10.5194/gmd-6-929-2013

Web End = , 2013.

Monteith, J. L. (Ed.): Vegetation and the atmosphere, vol. 1, Academic Press, 1975.

Nakai, Y., Sakamoto, T., Terajima, T., Kitamura, K., and Shirai, T.: The effect of canopy-snow on the energy balance above a coniferous forest, Hydrol. Process., 13, 23712382, 1999.

Napoly, A., Boone, A., Samuelsson, P., Gollvik, S., Martin, E., Seferian, R., Carrer, D., Decharme, B., and Jarlan, L.: The Interactions between Soil-Biosphere-Atmosphere (ISBA) land surface model Multi-Energy Balance (MEB) option in SURFEX Part 2: Model evaluation for local scale forest sites, Geosci. Model Dev. Discuss., doi:http://dx.doi.org/10.5194/gmd-2016-270, in review

Web End =10.5194/gmd-2016-270, in review , 2016.Niu, G.-Y. and Yang, Z.-L.: Effects of vegetation canopy processes on snow surface energy and mass balances, J. Geophys. Res., 109, d23111, doi:http://dx.doi.org/10.1029/2004JD004884

Web End =10.1029/2004JD004884 http://dx.doi.org/10.1029/2004JD004884

Web End = , 2004.

Noilhan, J. and Mahfouf, J.-F.: The ISBA land surface parameterisation scheme, Global Planet. Change, 13, 145159, 1996.Noilhan, J. and Planton, S.: A simple parameterization of land surface processes for meteorological models, Mon. Weather Rev., 117, 536549, 1989.

Oge, J. and Brunet, Y.: A forest oor model for heat and moisture including a litter layer, J. Hydrol., 255, 212233, 2002.

Oleson, K. W., Lawrence, D. M., Bonan, G. B., Flanner, M. G., Kluzek, E., Lawrence, P. J., Levis, S., Swenson, S. C.,P. E. Thornton, A. D., A., Decker, M., Dickinson, R., Feddema,J., C. L. Heald, F. H., Lamarque, J. F., Mahowald, N., Niu, G. Y., Qian, T., Randerson, J., Running, S., Sakaguchi, K., A. Slater,R. S., Wang, A., Yang, Z. L., and Zeng, X.: Technical Description of version 4.0 of the Community Land Model (CLM), NCAR Technical Note TN-478+STR, NCAR, NCAR, P.O. Box 3000, Boulder, CO, USA, 80307-3000, 2010.

Pitman, A. J.: The evolution of, and revolution in, land surface schemes designed for climate models, Int. J. Clim., 23, 479510, 2003.

Pokhrel, Y. N., Koirala, S., Yeh, P. J.-F., Hanasaki, N., Longuevergne, S. K., and Oki, T.: Incorporation of groundwater pumping in a global land surface model with the representation of human impacts, Water Resour. Res., 51, 7896, 2015.

Polcher, J., McAvaney, B., Viterbo, P., Gaertner, M.-A., Hahamann,A., Mahfouf, J.-F., Noilhan, J., Phillips, T., Pitman, A., Schlosser,C., Schulz, J.-P., Timbal, B., Verseghy, D., and Xue, Y.: A proposal for a general interface between land-surface schemes and general circulation models, Global Planet. Change, 19, 263278, 1998.

Pomeroy, J. W. and Dion, K.: Winter radiation extinction and reection in a boreal pine canopy: measurements and modelling, Hydrol. Process., 10, 15911608, 1996.

Richtmeyer, R. and Morton, K.: Difference method for initial values problems, Interscience Publishers, 2nd Edn., 1967.

Rutter, N., Essery, R., Pomeroy, J., Altimir, N., Andreadis, K., Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Douville, H., Dutra, E., Elder, K., Ellis, C., Feng, X., Gelfan, A., Goodbody,A., Gusev, Y., Gustafsson, D., Hellstrm, R., Hirabayashi, Y., Hirota, T., Jonas, T., Koren, V., Kuragina, A., Lettenmaier, D., Li,W.-P., Luce, C., Martin, E., Nasonova, O., Pumpanen, J., Pyles,R. D., Samuelsson, P., Sandells, M., Schdler, G., Shmakin,A., Smirnova, T. G., Sthli, M., Stckli, R., Strasser, U., Su,H., Suzuki, K., Takata, K., Tanaka, K., Thompson, E., Vesala,T., Viterbo, P., Wiltshire, A., Xia, K., Xue, Y., and Yamazaki,T.: Evaluation of forest snow processes models (SnowMIP2), J.

Geophys. Res., 114, doi:http://dx.doi.org/10.1029/2008JD011063

Web End =10.1029/2008JD011063 http://dx.doi.org/10.1029/2008JD011063

Web End = , 2009.

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

Web End =10.5194/gmd-9-223-2016 http://dx.doi.org/10.5194/gmd-9-223-2016

Web End = , 2016.

Samuelsson, P., Gollvik, S., and Ullerstig, A.: The land-surface scheme of the Rossby Centre regional atmospheric climate model (RCA3)., Report in Meteorology 122, SMHI, SE-60176 Norrkping, Sweden, 2006.

Samuelsson, P., Jones, C., Willn, U., Ullerstig, A., Gollvik,S., Hansson, U., Jansson, C., Kjellstrm, E., Nikulin, G., and Wyser, K.: The Rossby Centre Regional Climate Model RCA3: Model description and performance, Tellus A, 63, 13, doi:http://dx.doi.org/10.1111/j.1600-0870.2010.00478.x

Web End =10.1111/j.1600-0870.2010.00478.x http://dx.doi.org/10.1111/j.1600-0870.2010.00478.x

Web End = , 2011.

Schmidt, R. A. and Gluns, D. R.: Snowfall interception on branches of three conifer species, Can. J. Forest Res., 21, 12621269, 1991.

Seity, Y., Brousseau, P., Malardel, S., Hello, G., Bnard, P., Bouttier,F., Lac, C., and Masson, V.: The AROME-France Convective-Scale Operational Model, Mon. Weather Rev., 131, 976991, 2011.

Sellers, P. J., Mintz, Y., Sud, Y. C., and Dalcher, A.: A Simple Biosphere Model (SiB) for use within General Ciculation Models, J.Atmos. Sci., 43, 505531, 1986.

Sellers, P. J., Heiser, M. D., and Hall, F. G.: Relations between surface conductance and spectral vegetation indices at intermediate (100 m2 to 15 km2) length scales, J. Geophys. Res., 97, 19033 19059, doi:http://dx.doi.org/10.1029/92JD01096

Web End =10.1029/92JD01096 http://dx.doi.org/10.1029/92JD01096

Web End = , 1992.

Sellers, P. J., Randall, D. A., Collatz, G. J., Berry, J. A., Field,C. B., Dazlich, D. A., Zhang, C., Collelo, G. D., and Bounoua, L.: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model formulation, J. Climate, 9, 676705, 1996.Smith, B., Samuelsson, P., Wramneby, A., and Rummukainen, M.:

A model of the coupled dynamics of climate, vegetation and terrestrial ecosystem biogeochemistry for regional applications,

www.geosci-model-dev.net/10/843/2017/ Geosci. Model Dev., 10, 843872, 2017

872 A. Boone et al.: The interactions between ISBA-MEB in SURFEXv8

Tellus A, 63, 87106, doi:http://dx.doi.org/10.1111/j.1600-0870.2010.00477.x

Web End =10.1111/j.1600-0870.2010.00477.x http://dx.doi.org/10.1111/j.1600-0870.2010.00477.x

Web End = , 2011.van den Hurk, B., Best, M., Dirmeyer, P., Pitman, A., Polcher, J., and Santanello, J.: Acceleration of Land Surface Model Development over a Decade of Glass, B. Am. Meteorol. Soc., 92, 1591 1608, 2011.

Vergnes, J.-P., Decharme, B., and Habets, F.: Introduction of groundwater capillary rises using subgrid spatial variability of topography into the ISBA land surface model, J. Geophys. Res., 119, 1106511086, doi:http://dx.doi.org/10.1002/2014JD021573

Web End =10.1002/2014JD021573 http://dx.doi.org/10.1002/2014JD021573

Web End = , 2014.Voldoire, A., Sanchez-Gomez, E., y Mlia, D. S., B. Decharme,C. C., Snsi, S., Valcke, S., Beau, I., Alias, A., Chevallier, M., Dqu, M., Deshayes, J., Douville, H., E. Fernandez, G. M., Maisonnave, E., Moine, M. P., Planton, S., Saint-Martin, D.S. S., Tyteca, S., Alkama, R., Belamari, S., A. Braun, L. C., and Chauvin, F.: The CNRM-CM5.1 global climate model: description and basic evaluation, Clim. Dynam., 40, 20912121, doi:http://dx.doi.org/10.1007/s00382-011-1259-y

Web End =10.1007/s00382-011-1259-y http://dx.doi.org/10.1007/s00382-011-1259-y

Web End = , 2013.

Wilson, T., Meyers, T., Kochendorfer, J., Anderson, M., and Heuer,M.: The effect of soil surface litter residue on energy and carbon uxes in a deciduous forest, Agr. Forest Meteorol., 161, 134 147, 2012.

Wramneby, A., Smith, B., and Samuelsson, P.: Hot spots of vegetation-climate feedbacks under future greenhouse forcing in Europe, J. Geophys. Res.-Atmos., 115, D21119, doi:http://dx.doi.org/10.1029/2010JD014307

Web End =10.1029/2010JD014307 http://dx.doi.org/10.1029/2010JD014307

Web End = , 2010.

Wu, M., Schurgers, G., Rummukainen, M., Smith, B., Samuelsson,P., Jansson, C., Siltberg, J., and May, W.: Vegetationclimate feedbacks modulate rainfall patterns in Africa under future climate change, Earth Syst. Dynam., 7, 627647, doi:http://dx.doi.org/10.5194/esd-7-627-2016

Web End =10.5194/esd- http://dx.doi.org/10.5194/esd-7-627-2016

Web End =7-627-2016 , 2016.

Xue, Y., Sellers, P. J., Kinter, J. L., and Shukla, J.: A simplied Biosphere Model for Global Climate Studies, J. Climate, 4, 345 364, 1991.

Zhang, W., Jansson, C., Miller, P. A., Smith, B., and Samuelsson,P.: Biogeophysical feedbacks enhance the Arctic terrestrial carbon sink in regional Earth system dynamics, Biogeosciences, 11, 55035519, doi:http://dx.doi.org/10.5194/bg-11-5503-2014

Web End =10.5194/bg-11-5503-2014 http://dx.doi.org/10.5194/bg-11-5503-2014

Web End = , 2014.

Zhang, Z., Xue, Y., MacDonald, G., Cox, P. M., and Collatz, G. J.: Investigation of North American vegetation variability under recent climate A study using the SSiB4/TRIFFID biophysical/dynamic vegetation model, J. Geophys. Res., 120, 1300 1321, doi:http://dx.doi.org/10.1002/2014JD021963

Web End =10.1002/2014JD021963 http://dx.doi.org/10.1002/2014JD021963

Web End = , 2015.

Geosci. Model Dev., 10, 843872, 2017 www.geosci-model-dev.net/10/843/2017/