%  LaTeX support: latex@mdpi.com 
%  In case you need support, please attach all files that are necessary for compiling as well as the log file, and specify the details of your LaTeX setup (which operating system and LaTeX version / tools you are using).

%=================================================================
\documentclass[journal,article,submit,moreauthors,pdftex]{Definitions/mdpi} 
\usepackage{float}
\usepackage{subfigure}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{adjustbox}
\usepackage{chngcntr}
\usepackage{booktabs}
% If you would like to post an early version of this manuscript as a preprint, you may use preprint as the journal and change 'submit' to 'accept'. The document class line would be, e.g., \documentclass[preprints,article,accept,moreauthors,pdftex]{mdpi}. This is especially recommended for submission to arXiv, where line numbers should be removed before posting. For preprints.org, the editorial staff will make this change immediately prior to posting.

%--------------------
% Class Options:
%--------------------
%----------
% journal
%----------
% Choose between the following MDPI journals:
% acoustics, actuators, addictions, admsci, aerospace, agriculture, agriengineering, agronomy, algorithms, animals, antibiotics, antibodies, antioxidants, applsci, arts, asc, asi, atmosphere, atoms, axioms, batteries, bdcc, behavsci , beverages, bioengineering, biology, biomedicines, biomimetics, biomolecules, biosensors, brainsci , buildings, cancers, carbon , catalysts, cells, ceramics, challenges, chemengineering, chemistry, chemosensors, children, cleantechnol, climate, clockssleep, cmd, coatings, colloids, computation, computers, condensedmatter, cosmetics, cryptography, crystals, dairy, data, dentistry, designs , diagnostics, diseases, diversity, drones, econometrics, economies, education, ejihpe, electrochem, electronics, energies, entropy, environments, epigenomes, est, fermentation, fibers, fire, fishes, fluids, foods, forecasting, forests, fractalfract, futureinternet, futurephys, galaxies, games, gastrointestdisord, gels, genealogy, genes, geohazards, geosciences, geriatrics, hazardousmatters, healthcare, heritage, highthroughput, horticulturae, humanities, hydrology, ijerph, ijfs, ijgi, ijms, ijns, ijtpp, informatics, information, infrastructures, inorganics, insects, instruments, inventions, iot, j, jcdd, jcm, jcp, jcs, jdb, jfb, jfmk, jimaging, jintelligence, jlpea, jmmp, jmse, jnt, jof, joitmc, jpm, jrfm, jsan, land, languages, laws, life, literature, logistics, lubricants, machines, magnetochemistry, make, marinedrugs, materials, mathematics, mca, medicina, medicines, medsci, membranes, metabolites, metals, microarrays, micromachines, microorganisms, minerals, modelling, molbank, molecules, mps, mti, nanomaterials, ncrna, neuroglia, nitrogen, notspecified, nutrients, ohbm, optics, particles, pathogens, pharmaceuticals, pharmaceutics, pharmacy, philosophies, photonics, physics, plants, plasma, polymers, polysaccharides, preprints , proceedings, processes, proteomes, psych, publications, quantumrep, quaternary, qubs, reactions, recycling, religions, remotesensing, reports, resources, risks, robotics, safety, sci, scipharm, sensors, separations, sexes, signals, sinusitis, smartcities, sna, societies, socsci, soilsystems, sports, standards, stats, surfaces, surgeries, sustainability, symmetry, systems, technologies, test, toxics, toxins, tropicalmed, universe, urbansci, vaccines, vehicles, vetsci, vibration, viruses, vision, water, wem, wevj

%---------
% article
%---------
% The default type of manuscript is "article", but can be replaced by: 
% abstract, addendum, article, benchmark, book, bookreview, briefreport, casereport, changes, comment, commentary, communication, conceptpaper, conferenceproceedings, correction, conferencereport, expressionofconcern, extendedabstract, meetingreport, creative, datadescriptor, discussion, editorial, essay, erratum, hypothesis, interestingimages, letter, meetingreport, newbookreceived, obituary, opinion, projectreport, reply, retraction, review, perspective, protocol, shortnote, supfile, technicalnote, viewpoint
% supfile = supplementary materials

%----------
% submit
%----------
% The class option "submit" will be changed to "accept" by the Editorial Office when the paper is accepted. This will only make changes to the frontpage (e.g., the logo of the journal will get visible), the headings, and the copyright information. Also, line numbering will be removed. Journal info and pagination for accepted papers will also be assigned by the Editorial Office.

%------------------
% moreauthors
%------------------
% If there is only one author the class option oneauthor should be used. Otherwise use the class option moreauthors.

%---------
% pdftex
%---------
% The option pdftex is for use with pdfLaTeX. If eps figures are used, remove the option pdftex and use LaTeX and dvi2pdf.

%=================================================================
\firstpage{1} 
\makeatletter 
\setcounter{page}{\@firstpage} 
\makeatother
\pubvolume{xx}
\issuenum{1}
\articlenumber{5}
\pubyear{2019}
\copyrightyear{2019}
%\externaleditor{Academic Editor: name}
\history{Received: date; Accepted: date; Published: date}
%\updates{yes} % If there is an update available, un-comment this line

%% MDPI internal command: uncomment if new journal that already uses continuous page numbers 
%\continuouspages{yes}

%------------------------------------------------------------------
% The following line should be uncommented if the LaTeX file is uploaded to arXiv.org
%\pdfoutput=1

%=================================================================
% Add packages and commands here. The following packages are loaded in our class file: fontenc, calc, indentfirst, fancyhdr, graphicx, lastpage, ifthen, lineno, float, amsmath, setspace, enumitem, mathpazo, booktabs, titlesec, etoolbox, amsthm, hyphenat, natbib, hyperref, footmisc, geometry, caption, url, mdframed, tabto, soul, multirow, microtype, tikz

%=================================================================
%% Please use the following mathematics environments: Theorem, Lemma, Corollary, Proposition, Characterization, Property, Problem, Example, ExamplesandDefinitions, Hypothesis, Remark, Definition, Notation, Assumption
%% For proofs, please use the proof environment (the amsthm package is loaded by the MDPI class).

%=================================================================
% Full title of the paper (Capitalized)
\Title{Supplementary Material for ECLand: The ECMWF Land surface modelling system}

% Author Orchid ID: enter ID or remove command
\newcommand{\orcidauthorA}{0000-0001-8646-8701} % Add \orcidA{} behind the author's name
\newcommand{\orcidauthorB}{0000-0002-1745-3634} % Add \orcidB{} behind the author's name
\newcommand{\orcidauthorC}{0000-0002-6564-1699}
\newcommand{\orcidauthorD}{0000-0002-0643-2643}
\newcommand{\orcidauthorE}{0000-0003-0783-0448}
%\newcommand{\orcidauthorF}{}
\newcommand{\orcidauthorG}{0000-0002-7495-7273}
\newcommand{\orcidauthorH}{0000-0003-2794-2631}
\newcommand{\orcidauthorI}{0000-0001-5324-2749}
\newcommand{\orcidauthorJ}{0000-0002-5997-290X}
\newcommand{\orcidauthorK}{0000-0002-7374-3820}
\newcommand{\orcidauthorL}{0000-0002-1215-3288}
\newcommand{\orcidauthorM}{0000-0002-9049-1952}
\newcommand{\orcidauthorN}{0000-0001-7582-6497}
%\newcommand{\orcidauthorO}{}
\newcommand{\orcidauthorP}{0000-0003-1722-2497}
\newcommand{\orcidauthorQ}{0000-0002-6478-1841}
\newcommand{\orcidauthorR}{0000-0002-7998-0130}
%\newcommand{\orcidauthorS}{}



% AUTHORS GpComment (ideally CP-Team plus ...
\Author{\orcidA{}Souhail Boussetta $^{1}$, \orcidB{}Gianpaolo Balsamo $^{1}$, \orcidC{}Gabriele Arduini$^{1}$, \orcidD{}Emanuel Dutra$^{2,3}$, \orcidE{}Joe McNorton$^{1}$,  Margarita Choulga$^{1}$,\orcidG{}Anna Agust\'{\i}-Panareda$^{1}$, \orcidH{}Anton Beljaars$^{1}$, \orcidI{}Nils Wedi$^{1}$, \orcidJ{}Joaqu\'{\i}n Munõz-Sabater$^{1}$, \orcidK{}Patricia De Rosnay$^{1}$, \orcidL{}Irina Sandu$^{1}$, \orcidM{}Ioan Hadade$^{1}$, \orcidN{}Glenn Carver$^{1}$,  Cinzia Mazzetti$^{1}$, \orcidP{}Christel Prudhomme$^{1}$, \orcidQ{}Dai Yamazaki$^{4}$, \orcidR{}Ervin Zsoter $^{1}$}


% Authors, for metadata in PDF
\AuthorNames{Firstname Lastname, Firstname Lastname and Firstname Lastname}

% Affiliations / Addresses (Add [1] after \address if there is only one affiliation.)
\address{%
$^{1}$ \quad European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG2 9AX, UK 1; Souhail.Boussetta@ecmwf.int (S.B.)\\
$^{2}$ \quad Instituto Português do Mar e da Atmosfera, 1749-077 Lisbon, Portugal emanuel.dutra@ipma.pt (E.D.)\\
$^{3}$ \quad Instituto Dom Luiz, IDL, Faculty of Sciences, University of Lisbon, 1749-016 Lisbon, Portugal; \\
$^{4}$ \quad Institute of Industrial Science, The University of Tokyo, Tokyo, Japan; yamadai@rainbow.iis.u-tokyo.ac.jp (D.Y.)\\
}

% Contact information of the corresponding author
\corres{Correspondence: souhail.boussetta@ecmwf.int (S.B.); gianpaolo.balsamo@ecmwf.int (G.B.)}
% The commands \thirdnote{} till \eighthnote{} are available for further notes

%\simplesumm{} % Simple summary

%\conference{} % An extended version of a conference paper


%(list three to ten pertinent keywords specific to the article, yet reasonably common within the subject discipline.)

% The fields PACS, MSC, and JEL may be left empty or commented out if not applicable
%\PACS{J0101}
%\MSC{}
%\JEL{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Diversity
%\LSID{\url{http://}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Applied Sciences:
%\featuredapplication{Authors are encouraged to provide a concise description of the specific application or a potential application of the work. This section is not mandatory.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Data:
%\dataset{DOI number or link to the deposited data set in cases where the data set is published or set to be published separately. If the data set is submitted and will be published as a supplement to this paper in the journal Data, this field will be filled by the editors of the journal. In this case, please make sure to submit the data set as a supplement when entering your manuscript into our manuscript editorial system.}

%\datasetlicense{license under which the data set is made available (CC0, CC-BY, CC-BY-SA, CC-BY-NC, etc.)}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Only for the journal Toxins
%\keycontribution{The breakthroughs or highlights of the manuscript. Authors can write one or two sentences to describe the most important part of the paper.}
%\setcounter{secnumdepth}{4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{ {./Figures/} }

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setcounter{section}{-1} %% Remove this when starting to work on the template.
%\section{How to Use this Template}
%The template details the sections that can be used in a manuscript. Note that the order and names of article sections may differ from the requirements of the journal (e.g., the positioning of the Materials and Methods section). Please check the instructions for authors page of the journal to verify the correct order and names. For any questions, please contact the editorial office of the journal or support@mdpi.com. For LaTeX related questions please contact latex@mdpi.com.
%The order of the section titles is: Introduction, Materials and Methods, Results, Discussion, Conclusions for these journals: aerospace,algorithms,antibodies,antioxidants,atmosphere,axioms,biomedicines,carbon,crystals,designs,diagnostics,environments,fermentation,fluids,forests,fractalfract,informatics,information,inventions,jfmk,jrfm,lubricants,neonatalscreening,neuroglia,particles,pharmaceutics,polymers,processes,technologies,viruses,vision

%ED - The introduction reads fine, and has the mains components. The coupled part (e.g. GLACE2) could be further explored, and it also worth mentioning that the model is used in the EC-EARTH Earth System model 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







\counterwithin{table}{section}

\section{Current model description}
\label{appendixA}
\subsection{Tiles and surface fluxes}
\subsubsection{Tile and vegetation characteristics}\label{P4S7.2.1}

Grid-box surface fluxes are calculated separately for the
different subgrid surface fractions (or ``tiles''), leading to a
separate solution of the surface energy balance equation and skin
temperature for each of these tiles. This is an analogue of the
``mosaic'' approach of~\citet{KosterandSuarez1992}. Note that the
tiles at the interface soil-atmosphere are in energy and
hydrological contact with one single atmospheric profile above and
one single soil profile below. Each grid box is divided into eight
fractions: two vegetated fractions (high and low vegetation
without snow), one bare soil fraction, three snow/ice fractions
(snow on bare ground/low vegetation, high vegetation with snow
beneath, and sea-ice, respectively), and two water fractions
(interception reservoir, ocean/lakes). The tile for ``high
vegetation with snow beneath'' is a combined tile with a separate
energy balance and evaporation model for the high vegetation and
the underlying snow. A mixture of land and water (ocean/inland
water) tiles is not allowed, i.e. a grid box is either 100\% land
or 100\% sea.

In each grid box two vegetation types are present: a high and a
low vegetation type. An external climate database, based on the
Global Land Cover Characteristics (GLCC) data has been derived using one year of Advanced Very High Resolution Radiometer(AVHRR) data and ancillary information
 \citep{Lovelandetal2000}. The nominal resolution is 1~km. The data used provides each pixel with a biome classification based on the Biosphere-Atmosphere Transfer Scheme (BATS)
model~\citep{Dickinsonetal1993}, and four parameters have been
derived for each grid box: dominant vegetation type, $T_{\rm H}$
and $T_{\rm L}$, and the area fraction, $A_{\rm H}$ and $A_{\rm
L}$, for each of the high- and low-vegetation components,
respectively.

The coverage $C_{i}$ for the tile $i$ depends on the type and
relative area of low and high vegetation, and the presence of snow
and intercepted water. In the absence of snow and interception,
the vegetation coverage of high ($c_{\rm H}$) and low ($c_{\rm
L}$) vegetation are calculated as $A_{\rm H}c_{\rm veg}(T_{\rm
H})$ and $A_{\rm L}c_{\rm veg}(T_{\rm L})$, respectively, with
$c_{\rm veg}$ a vegetation type dependent density (see Table 
\ref{P4T7.1}). The bare ground fraction $c_{\rm B}$ is the
residual.
\begin{equation}%eqA1
\begin{split}
c_{\rm H} &= A_{\rm H} c_{\rm veg}(T_{\rm H}) \\
c_{\rm L} &= A_{\rm L} c_{\rm veg}(T_{\rm L}) \\
c_{\rm B} &= (1- c_{\rm H} - c_{\rm L})
\end{split}
 \label{eqA1}
\end{equation}
Each vegetation type is characterized by a series of (fixed)
parameters as detailed in Table \ref{P4T7.1}.
\begin{enumerate}
\item[(i)] A minimum canopy resistance, $r_{\rm s,min}$.

\item[(ii)] A leaf area index, {\it LAI}.

\item[(iii)] A vegetation coverage, $c_{\rm veg}$.

\item[(iv)] A coefficient, $g_{\rm D}$, for the dependence of the
canopy resistance, $r_{\rm c}$, on water vapour pressure deficit.

\item[(v)] The root distribution over the soil layers, specified by
an exponential profile involving attenuation coefficients, $a_{\rm
r}$, and $b_{\rm r}$.
\end{enumerate}

%T7.1
\begin{table}
 \centering
\caption{Vegetation types and parameter values (see text). H/L
refer to the distinction between high and low
vegetation.}
\label{P4T7.1}
\small
\begin{tabular}{clcccccc}
\toprule
 &  &  & $\boldsymbol{r_{\rm s,min}}$ &  & $\boldsymbol{g_{\rm D}}$ & &  \\
{\bf Index} & \multicolumn{1}{c}{\bf Vegetation type} & {\bf H/L} & $\boldsymbol{({\bf sm}^{-1})}$ & $\boldsymbol{c_{\rm veg}}$ & $\boldsymbol{({\rm hPa}^{-1})}$ &$\boldsymbol{a_{\rm r}}$ & $\boldsymbol{b_{\rm r}}$ \\
%\midrule
\hline \hline
\phantom{0}1 & Crops, mixed farming & L & 100 & 0.90
& 0\phantom{00.} & \phantom{0}5.558 & 2.614 \\
\phantom{0}2 & Short grass & L & 100 & 0.85 & 0\phantom{00.} & 10.739 & 2.608 \\
\phantom{0}3& Evergreen needleleaf trees & H & 250 & 0.90 & 0.03 & \phantom{0}6.706 & 2.175 \\
\phantom{0}4 & Deciduous needleleaf trees & H & 250 & 0.90 & 0.03 & \phantom{0}7.066 & 1.953 \\
\phantom{0}5 & Deciduous broadleaf trees & H & 175 & 0.90 & 0.03 & \phantom{0}5.990 & 1.955 \\
\phantom{0}6 & Evergreen broadleaf trees & H & 240 & 0.99 & 0.03 & \phantom{0}7.344 & 1.303 \\
\phantom{0}7 & Tall grass & L & 100 & 0.70 & 0\phantom{00.} & \phantom{0}8.235 & 1.627 \\
\phantom{0}8 & Desert & -- & 250 & 0\phantom{00.} & 0\phantom{00.} & \phantom{0}4.372 & 0.978 \\
\phantom{0}9 & Tundra & L & \phantom{0}80 & 0.50 & 0\phantom{00.} & \phantom{0}8.992 & 8.992 \\
10 & Irrigated crops & L & 180 & 0.90 & 0\phantom{00.} & \phantom{0}5.558 & 2.614 \\
11 & Semidesert & L & 150 & 0.10 & 0\phantom{00.} & \phantom{0}4.372 & 0.978 \\
12 & Ice caps and glaciers & -- & -- & -- & -- & -- & -- \\
13 & Bogs and marshes & L & 240 & 0.60 & 0\phantom{00.} & \phantom{0}7.344 & 1.303 \\
14 & Inland water & -- & -- & -- & -- & -- & -- \\
15 & Ocean & -- & -- & -- & -- & -- & -- \\
16 & Evergreen shrubs & L & 225 & 0.50 & 0\phantom{00.} & \phantom{0}6.326 & 1.567 \\
17 & Deciduous shrubs & L & 225 & 0.50 & 0\phantom{00.} & \phantom{0}6.326 & 1.567 \\
18 & Mixed forest/woodland & H & 250 & 0.90 & 0.03 & \phantom{0}4.453 & 1.631 \\
19 & Interrupted forest & H & 175 & 0.90 & 0.03 & \phantom{0}4.453 & 1.631 \\
20 & Water and land mixtures & L & 150 & 0.60 & 0\phantom{00.} & -- & -- \\
\bottomrule
\end{tabular}
\end{table}


The numerical values for the parameters of Table \ref{P4T7.1} are
based both on experiments conducted as described
in \citet{vandenHurketal2000} and on literature review, in
particular~\citet{Mahfoufetal1995},~\citet{ManziandPlanton1994},~\citet{GiardandBazile2000},~\citet{DormanandSellers1989},~\citet{Bonan1994},~\citet{Pitmanetal1991},
and~\citet{Zengetal1998}. The values for minimum canopy resistance have been revised
for crops, needle-leaf forest and short-grass in \citet{Boussettaetal2013}.

The presence of snow and intercepted water dynamically modifies the
coverage fractions. The coverage of snow, $c_{\rm sn}$, has been revised to be linearly
dependent from the snow depth $D_{\rm sn}$ (units m),
therefore taking into account both the snow mass $S$ (units $\text{kg m}^{-2}$)  and
the snow density $\rho_{\rm sn}$ (units $\text{kg m}^{-3}$), to account for the hysteresis 
characterizing snow accumulation and melting processes. 

The interception reservoir fraction, $c_{1}$, is given by $W_{1}/W_{\rm
1m}$, with $W_{\rm 1m}$, the maximum value for the
intercepted water in the grid box, defined from the leaf area index
contributions from the high and low vegetation tiles. The water
contents of the interception reservoir, $W_{1}$ (units m), and
$S$ are prognostic quantities in the model. Snow cover is assumed to
be overlying vegetation and bare ground with the same fraction. The
interception reservoir occupies an identical fraction of all
snow-free tiles.
\begin{equation}%eq7.2
\begin{split}
c_{\rm sn} &= \min \biggl(1, \frac{ S / \rho_{\rm sn}}{D_{\rm cr}}\biggr) \\
W_{\rm 1m} &= W_{\rm lmax} [c_{\rm B} + c_{\rm H} \cdot {\it
LAI}(T_{\rm H}) + c_{\rm L} \cdot {\it LAI} (T_{\rm L})] \\
c_{1} &= \min \biggl(1, \frac{W_{\rm l}}{W_{\rm lm}}\biggr) \\
\end{split}
 \label{P4eq7.2}
\end{equation}
In the expressions above the minimum snow depth that ensures
complete coverage of the grid box is $D_{\rm cr} = 0.10~{\rm m}$
and the maximum water over a single layer of leaves or over bare
ground is $W_{1\max}=0.0002~{\rm m}$. The leaf area index
%{\it LAI}, is specified in \hyperTable{P4T7.1} as a function of surface type.
{\it LAI}, is prescribed from a MODIS satellite-based dataset
as detailed in \citet{Boussettaetal2013}.
The full set of fractional tile coverages is given by
(\ref{P4eq7.3}) and (\ref{P4eq7.4}), where the indexing of the
tiles is detailed in Table \ref{P4T7.2}. Since a mixture of land
and ocean tiles is not allowed, a grid box is either 100\% water
(open water and ice, with ice fraction $c_{\rm i}$):
\begin{equation}%eq7.3
\begin{split}
C_{1} &= 1 - c_{i} \\
C_{2} &= c_{i} \\
C_{i} &= 0, i \in [3, N_{\rm T}] \\
\end{split}
 \label{P4eq7.3}
\end{equation}
or 100\% land (tiles 3 to $N_{\rm T}$, where $N_{\rm T}=8$ is the number of tiles):
\begin{equation}%eq7.4
\begin{split}
C_{1} &= C_{2} = 0 \\
C_{3} &= (1 - c_{\rm sn}) \cdot c_{1} \\
C_{4} &= (1 - c_{\rm sn}) \cdot (1 - c_{1}) \cdot c_{\rm L} \\
C_{5} &= c_{\rm sn} \cdot (1 - c_{\rm H}) \\
C_{6} &= (1 - c_{\rm sn}) \cdot (1 - c_{1}) \cdot c_{\rm H} \\
C_{7} &= c_{\rm sn} \cdot c_{\rm H} \\
C_{8} &= (1 - c_{\rm sn}) \cdot (1 - c_{1}) \cdot (1 - c_{\rm L} -
c_{\rm H})
\end{split}
 \label{P4eq7.4}
\end{equation}
Apart from the fractional gridbox coverage, each tile has a couple
of additional parameters (see Table \ref{P4T7.2}).
\begin{enumerate}
\item[(i)] The skin conductivity, $\Lambda_{\rm sk}$, provides the
thermal connection between the skin level and the soil or snow deck.
For high vegetation, $\Lambda_{\rm sk}$, is different for a stable and
unstable stratification of the temperature gradient between the skin
level and the upper soil or snow layer. This difference is considered
to represent the asymmetric coupling between the ground surface and
the tree canopy layer: an effective convective transport within the
tree trunk space for unstable conditions, and a limited turbulent
exchange for stable stratification \citep{Bosveldetal1999}.

\item[(ii)] A small fraction $f_{R_{\rm s}}$ of net short-wave radiation
that is transmitted directly to the top soil or snow layer. The
remaining fraction of the short-wave radiation ($1 - f_{R_{\rm s}}$) is absorbed
by the skin layer.
\end{enumerate}
Finally, the surface albedo, $\alpha_{i}$, is similar for all land
tiles within a grid box except for those covered with snow (see
the snow scheme description below). The climate database provides
the snow-free background albedo on a monthly basis. Long-wave
emissivity, $\varepsilon$, outside the window region is equal to
0.99 for all tiles; emissivity in the window region is tile
dependent and varies between 0.93 and 0.98. The remaining surface
characteristics (roughness length for momentum, $z_{\rm 0m}$, and
heat, $z_{\rm 0h}$) are derived from vegetation type, vegetation cover, and snow cover through Table \ref{table_z0}.

%T7.2
%\begin{table}
% \centering
%\caption{Tile specific values.\label{P4T7.2}}
%\small
%\begin{tabular}{p{4pc}p{6pc}p{5.5pc}p{5.5pc}p{4pc}p{5pc}}
%\toprule
%\hbox{}\par{\bf Index} & \hbox{}\par{\bf Tile} & $\boldsymbol{\Lambda_{\rm sk}}$ {\bf unstable}\par $\boldsymbol{({\bf Wm}^{-2}{\rm K}^{-1})}$ & $\boldsymbol{\Lambda_{\rm sk}}$ {\bf stable}\par $\boldsymbol{({\bf W m}^{-2}{\rm K}^{-1})}$ & \hbox{}\par$\boldsymbol{f_{R_{\rm s}}}$ & {\bf Resistance} \par {\bf scheme} \\
%\midrule
%1 & Open water & $\infty$ & $\infty$ & 0\phantom{.00} & Potential \\
%2 & Ice water & 58 & 58 & 0\phantom{.00} & Potential \\
%3 & Interception \par reservoir & 10 & 10 & 0.05 & Potential \\
%4 & Low vegetation & 10 & 10 & 0.05 & Resistance \\
%5 & Snow on low \par vegetation/bare \par
%ground & 7 & 7 & 0\phantom{.00} & Potential \\
%6 & High vegetation & $\Lambda_{\rm a,u}+5$ & $\Lambda_{\rm a,s}+5$ & 0.03 & Resistance \\
%7 & High vegetation \par with snow \par beneath & $\Lambda_{\rm a,u}+5$ & $\Lambda_{\rm a,s}+5$ & 0.03 & Canopy
%and \par snow resistance \\
%8 & Bare ground & 15 & 15 & 0\phantom{.00} & Resistance \\
%\midrule
%\multicolumn{6}{p{35pc}}{The resistance scheme describes the way
%of coupling with the atmosphere: {\it Potential} denotes
%atmospheric resistance only; {\it Resistance} denotes aerodynamic
%resistance in series with a canopy or soil resistance; {\it Canopy
%and snow resistance} denotes a canopy resistance for the
%vegetation and an extra aerodynamic coupling to the snow surface
%(see \hyperFigs{P4F7.1}--\ref{P4F7.2} and
%\hyperSubsection{P4S7.2.2}). For tiles~6 and~7, $\Lambda_{\rm a,u} =
%15~\text{W m}^{-2}{\rm K}^{-1}$ and $\Lambda_{\rm a,s} = 10~\text{W
%m}^{-2}{\rm K}^{-1}$ represent the aerodynamic coupling between
%the canopy and the soil in the unstable and stable cases,
%respectively, and the factor 5 represents the long-wave radiative
%exchanges. {\it Unstable/stable} refers to the temperature
%gradient between the skin layer and the top soil or snow layer.} \\
%\bottomrule
%\end{tabular}
%\end{table}

%TNew7.2
\begin{table}
\centering
\caption{The resistance scheme describes the way of coupling with the atmosphere: {\it P=Potential} denotes
atmospheric resistance only; {\it R=Resistance} denotes aerodynamic
resistance in series with a canopy or soil resistance; 
In presence of snow on high vegetation {\it $R_s$=Canopy and snow resistance} 
a canopy resistance for the vegetation and an extra aerodynamic coupling to the snow surface are used (see figure \ref{P4F7.2} and section \ref{P4S7.2.2}). }
\small
\begin{tabular}{clcccc}
\toprule
%\begin{tabular}{|p{0.25in}|p{1.2in}|p{1.0in}|p{1.0in}|p{0.5in}|p{0.5in}|} \hline 
 & & $\boldsymbol{\Lambda_{\rm sk}}$ {\bf unstable}\par  & 
 $\boldsymbol{\Lambda_{\rm sk}}$ {\bf stable}\par  & 
 $\boldsymbol{f_{R_{\rm s}}}$ & 
 {\bf Scheme} 
 \\  
%
\bf{\textit{Type}} & \bf{\textit{Vegetation type}} & 
$\boldsymbol{({\bf Wm}^{-2}{\rm K}^{-1})}$ & 
$\boldsymbol{({\bf W m}^{-2}{\rm K}^{-1})}$ & 
$-$ & 
$-$ 
\\  \hline \hline 
%\midrule
%\textbf{\textit{Type code}} & \textbf{\textit{Vegetation type}} & \textbf{\textit{R${}_{0}$\newline }}[mgCO${}_{2}$ m${}^{-2}$s${}^{-1}$]\textbf{\textit{}} & 
%\textbf{\textit{g${}_{m}$*${\left(25 \right)}$\newline }}[mm s${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{g${}_{c}$\newline }}[mm s${}^{-1}$]\textbf{\textit{}} & 
%\textbf{\textit{Dmax${}^{*}$\newline }}[kg kg${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{Am,max${\left(25 \right)}$
%\newline}}[mgCO${}_{2}$         m${}^{-2}$ s${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{f${}_{0}$${}^{*}$\newline }}[-] & 
%\textbf{\textit{${\Gamma}{\left(25 \right)}$\newline }}[ppm]\textbf{\textit{}} \\  \hline \hline 
1 & Crops, mixed farming & 10.0 & 10.0 & 0.05 & $R$ \\ 
2 & Short grass & 10.0 & 10.0 & 0.05 & $R$ \\ 
3 & Evergreen needleleaf  & 10.0 & 10.0 & 0.03 & $R$,$R_s$ \\
4 & Deciduous needleleaf  & 10.0 & 10.0 & 0.03 & $R$,$R_s$ \\ 
5 & Deciduous broadleaf  & 10.0 & 10.0 & 0.03 & $R$,$R_s$ \\ 
6 & Evergreen broadleaf  & 10.0 & 10.0 & 0.035 & $R$,$R_s$ \\
7 & Tall grass & 10.0 & 10.0 & 0.05 & $R$ \\ 
8 & Desert & 15.0 & 15.0 & 0. & $R$  \\
9 & Tundra & 10.0 & 10.0 & 0.05 & $R$ \\ 
10 & Irrigated crops & 10.0 & 10.0 & 0.05 & $R$ \\ 
11 & Semidesert & 10.0 & 10.0 & 0.05 & $R$ \\ 
12 & Ice caps and glaciers & 58.0 & 58. & 0. & $P$ \\
13 & Bogs and marshes & 10.0 & 10.0 & 0.05& $R$ \\ 
14 & Inland water & $\infty$ & $\infty$ & 0. & $P$ \\ 
15 & Ocean & $\infty$ & $\infty$ & 0. & $P$ \\ 
16 & Evergreen shrubs & 10.0 & 10.0 & 0.05 & $R$ \\ 
17 & Deciduous shrubs & 10.0 & 10.0 & 0.05 & $R$ \\ 
18 & Mixed forest- Wood & 10.0 & 10.0 & 0.03 & $R$,$R_s$ \\ 
19 & Interrupted forest & 10.0 & 10.0 & 0.03 & $R$,$R_s$ \\ 
20 & Water -land mixtures & $\infty$ & $\infty$ & 0. & $P$ \\ %\hline 
\bottomrule
\label{P4T7.2}
\end{tabular}
\end{table}

\subsubsection{Surface momentum flux}\label{P4S7.2.1new}
The roughness lengths for momentum and heat are prescribed as constant values depending on each vegetation type following \citet{Mahfoufetal1995} since IFS cycle 31r1 as used in ERA-Interim re-analysis. 
In Cy37r3 some of these values were revised based on a study by \citet{Sandu-2011}. 
This study emphasized that the model generally tends to overestimate the 10-m winds compared to routine observations and 
that the forecast errors depend on the vegetation type because the representation of the 10-m wind speed is mainly controlled 
by the value of the momentum roughness length associated with each vegetation type.
This parameter is difficult to determine from observations and the values used so far in the ECMWF model were empirically chosen
as 20\% and 10\% of the vegetation height for forests and respectively for the other vegetation types \citep{Mahfoufetal1995}. 
The overestimation of the near-surface wind speed for most of the vegetation types suggests that the values used for the
momentum roughness length were too low. They were therefore revised based on theoretical considerations and synop
observations of wind speed at 10 m. The basic idea was to search, for each vegetation type, for a new value
of the momentum roughness length for which the mean 10-m wind speed forecast error with respect to synop
observations drops to zero. This calibration showed that the momentum roughness length values should be increased for
nine and decreased for one of the 18 vegetation types characterizing land areas. The newly derived values for ten of the
vegetation types were introduced in CY37r3 (Table \ref{table_z0}). As the roughness length
for momentum was on average increased, the roughness length for heat was decreased in order to account for terrain heterogeneity (Table \ref{table_z0}).

The wet skin tile obtains values that are weighted between low and high vegetation according to their fractional cover. 
The exposed snow tile is set to table entry 12 (ice caps and glaciers), the snow under high vegetation tile has the roughness length 
of the high vegetation and the bare soil tile is set to the roughness length of table entry 8 (desert).
The model uses the roughness lengths for the individual tiles and aggregates the fluxes.
For postprocessing also an aggregated roughness length field is computed by tile averaging the neutral transfer coefficients
and backing out the aggregated roughness lengths $z_{oma}$ and $z_{oha}$.
%
\begin{eqnarray}
\frac{1}{(ln10/z_{oma})^2}= \Sigma_i \frac{Fr_i}{(ln10/z_{omi})^2}  \label{eqn_z0ma} \\
\frac{1}{(ln10/z_{oha})^2}= \Sigma_i \frac{Fr_i}{(ln10/z_{ohi})^2}  \nonumber
\end{eqnarray}
%

% 
\begin{table}
 \centering
\caption{Roughness lengths for momentum and heat associated with high and low vegetation types.
\label{table_z0}}
\small
\begin{tabular}{clccc}
\toprule
{\bf Index} & \multicolumn{1}{c}{\bf Vegetation type} &
{\bf H/L veg} & $\boldsymbol{z_{0m}}$ & $\boldsymbol{z_{0h}}$\\
%\midrule
\hline \hline
\phantom{0}1 & Crops, mixed farming & L & 0.25 & 0.25 10$^{-2}$\\
\phantom{0}2 & Short grass & L & 0.1 & 0.1 10$^{-2}$ \\
\phantom{0}3& Evergreen needleleaf trees & H & 2.0 & 2.0 \\
\phantom{0}4 & Deciduous needleleaf trees & H & 2.0 & 2.0 \\
\phantom{0}5 & Deciduous broadleaf trees & H & 2.0 & 2.0 \\
\phantom{0}6 & Evergreen broadleaf trees & H & 2.0 & 2.0 \\
\phantom{0}7 & Tall grass & L & 0.47 & 0.47 10$^{-2}$ \\
\phantom{0}8 & Desert & -- & 0.013 & 0.013 10$^{-2}$ \\
\phantom{0}9 & Tundra & L & 0.034 & 0.034 10$^{-2}$ \\
10 & Irrigated crops & L & 0.5 & 0.5 10$^{-2}$ \\
11 & Semidesert & L & 0.17 & 0.17 10$^{-2}$ \\
12 & Ice caps and glaciers & -- & 1.3 10$^{-3}$ & 1.3 10$^{-4}$ \\
13 & Bogs and marshes & L & 0.5 & 0.5 10$^{-2}$ \\
14 & Inland water & -- & -- & -- \\
15 & Ocean & -- & -- & -- \\
16 & Evergreen shrubs & L & 0.100 & 0.1 10$^{-2}$ \\
17 & Deciduous shrubs & L & 0.25 & 0.25 10$^{-2}$ \\
18 & Mixed forest/woodland & H & 2.0 & 2.0 \\
19 & Interrupted forest & H & 1.1 & 1.1 \\
20 & Water and land mixtures & L & -- & -- \\
\bottomrule
\end{tabular}
\end{table}

\subsubsection{Surface energy balance and fluxes}\label{P4S7.2.2}

The surface energy balance is satisfied independently for each tile by
calculating its skin temperature. The skin layer represents the
vegetation layer, the top layer of the bare soil, or the top layer of
the snow pack. It has no heat capacity and therefore responds
instantaneously to changes in, e.g.  radiative forcing. In order to
calculate the skin temperature, the surface energy balance equation is
linearized for each tile leading to an expression for the skin
temperature. This procedure is equivalent to the Penmann--Monteith
approach which can be derived by eliminating the skin temperature from
the surface energy balance equation assuming that the net radiation
minus ground heat flux is known (e.g. \citet{Brutsaert1982}). The
approach followed here is an extension to the Penmann--Monteith
derivation in the sense that it allows for coupling with the
underlying soil (or snow, ice). Because of the short time scale
associated with the skin layer, the equation for its temperature is
solved implicitly together with the vertical turbulent transport in
the boundary layer. These are solved in a fully implicit approach where the skin
temperatures depend on each other following the coupler approach suggested by
\citet{Bestetal2004}.

The following general discussion applies to each tile but the
parameters are tile dependent. The surface energy balance equation can be written as:
%
\begin{equation}
R_{SW} + R_{LW} +H + { LJ}_q = \Lambda_{skin}(T_{sk}-T_{ s}) 
\label{sfc_energy_balance}
\end{equation}
%
where $R_{SW}$ and $R_{LW}$ are the net
short-wave and long-wave radiation fluxes at the surface and the
right-hand side represents the ground heat flux through coupling
with the underlying soil, snow or ice with temperature $T_s$.

A resistance parametrization is used to calculate the turbulent
fluxes. Momentum exchange is parametrized with the same roughness
length for all tiles, but with a different stability correction
for each tile. The resistance scheme for water vapour and heat
exchanges is different for different tiles (see figure \ref{P4F7.2}). For ocean, sea ice and snow on low vegetation,
the turbulent fluxes of heat and water vapour are given by
\begin{gather}%eq7.5
H_{i} = \rho_{\rm a} c_{p} |U_{\rm L}|C_{{\rm H},i} (T_{\rm L} +{\it
gz}_{\rm L} / c_{\rm p} - T_{{\rm sk},i})
 \label{P4eq7.5} \\%7.6
E_{i} = \rho_{\rm a} |U_{\rm L}|C_{{\rm H},i} [q_{\rm L} -q_{\rm sat}
(T_{{\rm sk},i})]
 \label{P4eq7.6}
\end{gather}
with $\rho_{\rm a}$ the air density, $c_{\rm p}$ the heat capacity of
moist air, $g$ the acceleration of gravity, $|U_{\rm L}|$, $T_{\rm
L}$, $q_{\rm L}$, $z_{\rm L}$ the wind speed, temperature,
humidity and height of the lowest atmospheric model level, and
$C_{{\rm H},i}$ the turbulent exchange coefficient, that varies
from tile to tile because of different atmospheric stabilities.
For the exchange coefficients, a Monin--Obukhov formulation is adopted for the stability dependence and the different roughness lengths for heat and momentum are assumed accordingly.

%F7.2
\begin{figure}
 \centering
\includegraphics[width=\textwidth]{Figures/figA1.eps}
\caption{Resistance scheme for three categories of coupling. {\it
Potential} refers to ocean, sea ice and snow on low vegetation;
{\it (Canopy) resistance} to dry low and dry high vegetation, bare
soil, and interception reservoir when potential evaporation
exceeds the maximum reservoir content; {\it Resistance to snow} to
snow under high vegetation.\label{P4F7.2}}
\end{figure}

For high and low vegetation, an additional canopy resistance $r_{\rm
c}$ is added with
\begin{equation}%eq7.7
E_{i} = \frac{\rho_{\rm a}}{r_{\rm a}+r_{\rm c}}[q_{\rm L} - q_{\rm
sat}(T_{{\rm sk},i})]
 \label{P4eq7.7}
\end{equation}
with $r_{\rm a} = (|U_{\rm L}|C_{{\rm H},i})^{-1}$ and $i$ indicating
the high or low vegetation tiles. $r_{\rm c}$ is a function of
downward short-wave radiation $R_{\rm s}$, leaf area index ${\it
LAI}$, average unfrozen root soil water $\bar{\theta}$, atmospheric
water vapour deficit $D_{\rm a}$ and a minimum stomatal resistance
$r_{\rm s,min}$, following~\citet{Jarvis1976} given by
\begin{equation}%eq7.8
r_{\rm c} = \frac{r_{\rm S,min}}{\it LAI}f_{1} (R_{\rm s})f_{2}
(\bar{\theta})f_{3}(D_{\rm a})
 \label{P4eq7.8}
\end{equation}
$f_{1}$ is a hyperbolic function of downward short-wave radiation only so that
\begin{equation}%eq7.9
\frac{1}{f_{1}(R_{\rm s})} = \min \biggl[1, \frac{{\it bR}_{\rm
s}+c}{a({\it bR}_{\rm s}+1)}\biggr]
 \label{P4eq7.9}
\end{equation}
where $a = 0.81$, $b =0.004\ {\rm W}^{-1} {\rm m}^{2}$ and
$c=0.05$.

Function $f_{2}$ is defined as
\begin{equation}%eq7.10
\frac{1}{f_{2}(\bar{\theta)}} =
\begin{cases}
0 & \bar{\theta}< \theta_{\rm pwp} \\
\dfrac{\bar{\theta} -\theta_{\rm pwp}}{\theta_{\rm cap} - \theta_{\rm
pwp}} & \theta_{\rm pwp}\leq \bar{\theta}\leq \theta_{\rm cap} \\
1 & \bar{\theta}>\theta_{\rm cap}
\end{cases}
 \label{P4eq7.10}
\end{equation}
where the soil moisture at permanent wilting point and at field
capacity, $\theta_{\rm pwp}$ and $\theta_{\rm cap}$, respectively,
are defined in Table \ref{P4T7.7}. $\bar{\theta}$ is a weighted
average of the unfrozen soil water given by
\begin{equation}%eq7.11
\bar{\theta} = \sum^{4}_{k=1} R_{k}\max [f_{{\rm liq},k}\theta_{k},
\theta_{\rm pwp}]
 \label{P4eq7.11}
\end{equation}
where $R_{k}$ is the fraction of roots in layer $k$ and the
fraction of unfrozen soil water, $f_{{\rm liq},k} = 1- f_{\rm
fr}(T_{k})$, is a parametrized function of the soil temperature
of layer $k$, $T_{k}$, as specified in \ref{P4S7.5.2}.
Table \ref{P4T7.1} lists the coefficients $a_{\rm r}$ and $b_{\rm
r}$ which are used to calculate the root fraction $R_{k}$
according to~\citet{Zengetal1998}:
\begin{equation}%eq7.12
R_{k} = 0.5 [\exp(-a_{\rm r}z_{k-1/2}) + \exp(-b_{\rm r}z_{k-1/2}) -
\exp(-a_{\rm r}z_{k+1/2}) - \exp(-b_{\rm r}z_{k+1/2})]
 \label{P4eq7.12}
\end{equation}
where $z_{k+1/2}$ is the depth of the bottom of layer $k$ (in m;
$z_{1/2} = 0$ m). Contributions from levels exceeding the column
depth are added to the deepest soil layer in order to ensure that
$\sum R_{k}=1$. Table \ref{P4T7.3} lists the distribution of the
roots over the four soil layers.

%T7.3
\begin{table}
\centering
\caption{Root distribution per vegetation type (in \%) over the
four layers. Vegetation indexes refer to Table \ref{P4T7.1}.\label{P4T7.3}}
\begin{tabular}{lcccccccccccccccc}
\toprule
{\bf Vegetation}\\
{\bf index} & {\bf 1} & {\bf 2} & {\bf 3} & {\bf 4} & {\bf 5} & {\bf 6} & {\bf 7} & {\bf 8} & {\bf 9} & {\bf 10} & {\bf 11} & {\bf 13} & {\bf 16} & {\bf 17} & {\bf 18} & {\bf 19} \\
%\midrule
\hline \hline
Layer 1 & 24 & 35 & 26 & 26 & 24 & 25 & 27 & 100 & 47 & 24 & 17 & 25 & 23 & 23 & 19 & 19 \\
Layer 2 & 41 & 38 & 39 & 38 & 38 & 34 & 27 & \phantom{00}0 & 45 & 41 & 31 & 34 & 36 & 36 & 35 & 35 \\
Layer 3 & 31 & 23 & 29 & 29 & 31 & 27 & 27 & \phantom{00}0 & \phantom{0}8 & 31 & 33 & 27 & 30 & 30 & 36 & 36 \\
Layer 4 & \phantom{0}4 & \phantom{0}4 & \phantom{0}6 & \phantom{0}7 &
\phantom{0}7 & 14 & \phantom{0}9 & \phantom{00}0 & \phantom{0}0 & \phantom{0}4 & 19 & 11 & 11 & 11 & 10 & 10 \\
\bottomrule
\end{tabular}
\end{table}

A dependence on atmospheric humidity deficit ($D_{\rm a}=e_{\rm
sat}(T_{\rm L})-e_{\rm L}$, with e the vapour pressure) is included
according to
\begin{equation}%eq7.13
\frac{1}{f_{3}(D_{\rm a})} = \exp (-g_{\rm D}D_{\rm a})
 \label{P4eq7.13}
\end{equation}
where $g_{\rm D}$ depends on the vegetation type
(Table \ref{P4T7.1}), and is non-zero for high vegetation only.

Evaporation from the interception reservoir is given by
(\ref{P4eq7.6}) only when the amount of water in the interception
reservoir, $W_{1}$, is sufficient to sustain potential evaporation
during the entire time step $\Delta t$. If $W_{1}$ is limited, an
additional resistance $r_{1}$, analogue to $r_{\rm c}$ in
(\ref{P4eq7.7}), is introduced. $r_{1}$ is calculated from the
potential evaporation of the previous time step. Note that this
type of flux-limiter is a time-step dependent feature of the model
numerics.

The evaporation from non-vegetated areas responds to a different physical mechanism
compared to densely vegetated areas.
Over bare soil the vaporisation of water in the soil pores takes place in a thin layer close
to the surface-atmosphere interface as a direct effect of incoming solar radiation providing
the latent heat requirements 2.5106~J/kg. Atmospheric conditions (e.g. air temperature,
humidity, wind velocity and radiation) as well as soil conditions
(e.g. water content, roughness length) all play a role in modulating the evaporation processes.

The relationship between soil moisture and bare soil evaporation is generally parametrized
in land surface models and several formulations have been compared in the study by~\citep{MahfoufandNoilhan1991}.
They found that in all cases the bare soil evaporation formulation stops when the soil is completely dry
(for soil moisture close to zero).
In the TESSEL scheme~\citep{vandenHurketal2000,ViterboandBeljaars1995,Viterboetal1999}
the soil-moisture-evaporation link was assumed to be varying linearly between the wilting point and
the field capacity values for soil moisture. With the introduction of a tiling approach,
the same stress function was applied to both vegetated and non-vegetated tiles,
neglecting the fact that wilting point is a soil moisture threshold that applies uniquely to vegetated areas.
A justification for this approach is that the soil underneath the surface tiles is still a unique column
for both vegetated and non-vegetated fractions.
In the most recent version of the HTESSEL scheme
\citep{Balsamoetal2011}, included in IFS Cy36r4 and operational since 10 November 2010),
the formulation of the bare soil evaporation has been revisited to allow a smooth transition
between vegetated and non-vegetated areas and to realign the formulation of bare ground evaporation
with studies in the literature.

The soil evaporation resistance, $r_{\rm soil}$, is
\begin{equation}%eq7.14
r_{\rm soil} = r_{\rm soil,min} f_{2b} (f_{\rm liq} \theta_{1})
 \label{P4eq7.14}
\end{equation}	
 	
with $r_{\rm soil,min} = 50~\text{s m}^{-1}$.
By this parametrization, evaporation from bare ground is treated similar to a single leaved canopy
with a minimum resistance $r_{\rm soil,min}$, extracting water from the upper soil layer only,
and not experiencing any additional stress due to limited radiation or dry air.				
In case of a vegetated tile the canopy resistance is formulated according to \citet{Jarvis1976}
and the soil moisture inhibition function (see (\ref{P4eq7.7})), which depends on the root-zone soil wetness,
normalized between the wilting point and the field capacity, for non-vegetated tiles
is uniquely dependent on the soil moisture of the first soil layer and a minimum soil resistance.
              
In the new formulation the $f_2b$ for bare ground is calculated as
\begin{equation}%eq7.14b
\frac{1}{f_{2b}({\theta_{1})}} =
\dfrac{{\theta_{1}} -\theta_{\rm min}}{\theta_{\rm cap} - \theta_{\rm min}} 
\label{P4eq7.14b}
\end{equation}

where ${\theta_{1}}$ is the unfrozen soil water in the first soil layer and
$\theta_{\rm min}$ is a vegetation-weighted average of the wilting point and a residual soil moisture
content ($\theta_{\rm res}$ which is a function of the soil textural class, as defined in Table \ref{P4T7.5}).
The weights are given by the vegetation cover $veg$ \citep{vandenHurketal2000} over the gridbox, so that 
\begin{equation}%eq7.14c
\theta_{\rm min}={veg} \theta_{\rm pwp} + (1-{veg}) \theta_{\rm res} 
\label{P4eq7.14c}
\end{equation}	
 
This modification of the soil moisture inhibition function allows a
much lower level of soil moisture to be reached in dry areas as a
consequence of direct bare ground evaporation under strong insulation.
Also results are more consistent with the soil moisture values observed, for instance, over deserts. 

%Bare-soil evaporation uses a resistance approach, an analogue to
%the canopy transpiration formulation (see (\ref{P4eq7.7})). The
%soil evaporation resistance, $r_{\rm soil}$, is
%\begin{equation}%eq7.14
%r_{\rm soil} = r_{\rm soil,min} f_{2\'} (f_{\rm liq} \theta_{1})
% \label{P4eq7.14}
%\end{equation}
%with $f_{2}$ given by (\ref{P4eq7.10}), and $r_{\rm soil,min} = 50~\text{s m}^{-1}$.
%By this parametrization, evaporation from bare
%ground is treated similar to a single leaved canopy with a minimum
%resistance $r_{\rm soil,min}$, extracting water from the upper
%soil layer only, and not experiencing any additional stress due to
%limited radiation or dry air.
%Equation (\ref{P4eq7.14}) shuts off
%evaporation when the top soil moisture reaches permanent wilting
%point. When compared to observations over semi-arid areas, an
%alternative relative humidity
%formulation~\citep{MahfoufandNoilhan1991,ViterboandBeljaars1995},
%that does not have a similar limitation, gave excessive
%evaporation~\citep{vandenHurketal2000}.

A special treatment is included in the calculation of evaporation
over high vegetation with snow underneath (see figure \ref{P4F7.2}).
Evaporation takes place from both the canopy component in the tile
($E_{\rm veg,7}$) and from the snow lying under the vegetation.
The canopy evaporation uses a canopy resistance and saturation
specific humidity at the canopy skin temperature, while the snow
evaporation $E_{\rm sn,7}$ is parametrized with an additional
constant aerodynamic resistance $r_{\rm a,sn}$ and saturation
specific humidity at snow temperature $T_{\rm sn}$. The
evaporation from tile~7 is the combination of the canopy
transpiration and the snow evaporation so that
\begin{equation}%eq7.15
E_{7} = \rho_{\rm a}\frac{q_{\rm L}-q_*}{r_{\rm a}} = \rho_{\rm
a}\frac{q_*-q_{\rm sat}(T_{\rm sn})}{r_{\rm a,sn}} + \rho_{\rm
a}\frac{q_*-q_{\rm sat}(T_{\rm sk})}{r_{\rm c}}
 \label{P4eq7.15}
\end{equation}
where $q_{*}$ is the humidity at the connection point of the three
resistances (figure \ref{P4F7.2}). After elimination of $q_{*}$,
$E_{7}$ can be rewritten as
\begin{equation}%eq7.16
E_{7} = \rho_{\rm a} \frac{q_{\rm L}-q_{\rm sat}(T_{\rm
sk})}{r_{\rm a}+r_{\rm c}+r_{\rm c}\frac{r_{\rm a}}{r_{\rm a,sn}}}
+ \rho_{\rm a} \frac{q_{\rm L}-q_{\rm sat}(T_{\rm sn})}{r_{\rm
a}+r_{\rm a,sn}+r_{\rm a,sn}\frac{r_{\rm a}}{r_{\rm c}}}
 \label{P4eq7.16}
\end{equation}
The first term in the equation above is interpreted as $E_{\rm
veg,7}$ and is treated in the standard way (i.e., implicit in the
tile skin temperature). The second term is interpreted as
evaporation from snow ($E_{\rm sn,7}$) and is handled explicitly.
The values of $r_{\rm a,sn}$ depend on the stability of the
subcanopy layer and are functions of $\Lambda_{\rm a,u}$ and
$\Lambda_{\rm a,s}$ (see Table \ref{P4T7.2}); $r_{\rm a,sn} = 67\
{\rm sm}^{-1}$ and $r_{\rm a,sn} = 220\ \text{s m}^{-1}$ for an
unstable and stable subcanopy layer, respectively. In spring, the
latent heat flux of that tile, $L_{\rm v}E_{\rm veg,7}+L_{\rm
s}E_{\rm sn,7}$ will be dominated by snow evaporation since the
frozen soil under the snow deck will lead to very large values of
$r_{\rm c}$.

The grid box total sensible and latent heat fluxes are expressed as
an area weighted average:
\begin{gather}%eq7.17
H = \sum^{8}_{i=1}C_{i}H_{i}
 \label{P4eq7.17} \\%7.18
E = \sum^{8}_{i=1}C_{i}E_{i}
 \label{P4eq7.18}
\end{gather}
with $H_{i}$ given by (\ref{P4eq7.5}), and $E_{i}$ by
(\ref{P4eq7.6}) for ocean, sea-ice and snow on low vegetation,
(\ref{P4eq7.7}) for dry high and low vegetation, the interception
reservoir (with $r_{\rm c}$ replaced by $r_{1}$) and for bare soil
(with $r_{\rm c}$ replaced by $r_{\rm soil}$) and (\ref{P4eq7.16})
for high vegetation with underlying snow.

\subsubsection{Land coupling}\label{P4S7.3}

A skin temperature $T_{\rm sk}$ forms the interface between the
soil and the atmosphere. It is calculated for each grid box tiles separately, by scanning the surface energy balance solver over the eight tiles, assuming a
complete coverage of the specific tile. For a single tile, this
procedure is very similar to the derivation of the Penman--Monteith
equation in which the skin temperature is eliminated from the
surface energy balance equation. The numerical approach used in
HTESSEL has the advantage that the feedback of skin temperature on
net radiation and ground heat flux is included. The input radiation and reference
atmospheric temperature ($T_{\rm L}$), specific humidity ($q_{\rm
L}$) and wind speed ($U_{\rm L}$) are identical for each tile. The
surface fluxes ``seen'' by the atmosphere are calculated as an
area-weighted average over the tiles (see (\ref{P4eq7.17}) and
(\ref{P4eq7.18})). For the high vegetation with snow underneath,
the skin temperature is that of the high vegetation; the
temperature of the underlying snow is calculated separately.

The energy balance equation solved for each tile takes into
account partial absorption of net short-wave radiation, $1-f_{{\rm
Rs},i}$, in the skin layer (see Table \ref{P4T7.2}). The remaining
energy is directly passed to the soil or snow so that
\begin{equation}%eq7.19
(1-f_{{\rm Rs},i})(1-\alpha_{i})R_{\rm s}+\varepsilon (R_{\rm T} -
\sigma T^{4}_{{\rm sk},i}) + H_{i} + L_{\rm v,s} E_{i} = \Lambda_{{\rm
sk},i}(T_{{\rm sk},i}-T_{1})
 \label{P4eq7.19}
\end{equation}
where $i$ denotes the tile index, $R_{\rm s}$ and $R_{\rm T}$ are
downward short-wave radiation and long-wave radiation,
respectively, $\sigma$ is the Stefan--Bolzman constant, $T_{1}$ the
temperature of the upper soil or snow layer, $H_{i}$ the sensible
heat flux, and $L_{\rm v,s}E_{i}$ the latent heat flux from the
skin layer, and $\Lambda_{{\rm sk},i}$, the skin conductivity for
tile~$i$. Latent heat of evaporation, $L_{\rm v}$, is used for all
evaporation terms except snow evaporation, while $L_{\rm s}$, the
latent heat of sublimation, is used for evaporation of snow (i.e.,
tile 5 and the contribution $E_{\rm sn,7}$ from tile 7, defined by
(\ref{P4eq7.16})).

The tiled surface is thermally coupled to the snow deck, when
present, and to a single soil profile. The net flux into the soil is
a weighted average of the flux from each tile.

The solution of (\ref{P4eq7.19}) is performed inside the code for
turbulent exchanges in the atmosphere. The atmospheric vertical diffusion equations yield a tridiagonal system of equations, with the coupling to the skin temperature
given by the matrix row corresponding to the lowest model level.
The first step for the solution of the system of equations, an LU decomposition, is followed by the solution of (\ref{P4eq7.19}) before back-substitution.

Due to the linearisation performed in (\ref{P4eq7.19}) the net long-wave radiation $R_{\rm LW}=\varepsilon (R_{\rm T} - \sigma T^{4}_{{\rm sk}})$ for each tile is re-computed at the end of the energy balance by linearising it around the new skin
temperature so that

\begin{equation}
\begin{split}
R_{{\rm LW},i}=R_{\rm LW}-\varepsilon 4 \sigma T^{3}_{{\rm sk}}( T_{{\rm sk},i}-T_{\rm sk}) \,     \\
T_{\rm sk} = \sum_i Fr^i T_{{\rm sk},i} \,
\label{eq_lwtil}
\end{split}
\end{equation}

where $T_{\rm sk}$ is the mean grid-box skin temperature and $R_{{\rm LW},i}$ the 
tiled net long-wave radiation used to compute the net energy passed to the 
remaining surface processes. 

\subsection{Snow}\label{P4S7.4}

The snow scheme represents an additional ``layer'' on top of the
upper soil layer, with an independent, prognostic, thermal and
mass contents. The snow pack is represented by a single snow
temperature, $T_{\rm sn}$ and the snow mass per unit area (snow
mass for short) $S$. The net energy flux at the top of the snow
pack, $G^{\rm T}_{\rm sn}$, is the residual of the skin energy
balance from the snow covered tiles and the snow evaporation from
the tile with high vegetation over snow (see (\ref{P4eq7.15})).
The basal heat flux, $G^{\rm B}_{\rm sn}$, is given by a
resistance formulation modified in case of melting. The absorbed
energy is used to change the snow temperature or melt the snow,
when $T_{\rm sn}$ exceeds the melting point. Liquid water within the
snow pack is diagnosed from snow temperature, snow mass and snow density. Liquid water coexists in the snowpack leading to internal phase changes (freezing/melting), and rainfall can be intercepted.

The heat capacity of the snow deck is a function of its depth,
and snow density, and is modified in the presence of liquid water.
The snow thermal conductivity
changes with changing snow density. Snow density changes due to overburden and thermal meta-morphisms \citep{Anderson1976,BooneAndEtchevers2001} and compaction related to
melt water retained in the snowpack, adapted from \citet{Lynchstieglitz1994}.
The snow albedo changes exponentially with snow age. For snow on low vegetation it ranges
between 0.50 for old snow and 0.85 for fresh snow. The
albedo for high vegetation with snow underneath depends on
vegetation type adapted from \citet{Moodyetal2007}.
The details of the snow scheme and its validation can be found
in \citet{Dutraetal2009} and \citet{Dutraetal2010}.

\subsubsection{Snow mass and energy budget}\label{P4S7.4.1}

The snow energy budget in the presence of snow liquid water reads as
\begin{equation}%eq7.23a
\begin{split}
(\rho C)_{\rm sn} D_{\rm sn} \frac{\partial T_{\rm sn}}{\partial t} &=
R^{\rm N}_{\rm sn}+L_{\rm s}E_{\rm sn}+H_{\rm sn}-G^{\rm B}_{\rm sn}-L_{\rm f}M_{\rm sn}-
Q^{\rm INT}_{\rm sn} \\
Q^{\rm INT}_{\rm sn} &= L_{\rm f}M^{\rm INT}_{\rm sn}=
L_{\rm f}\frac{\partial S_{\rm l}}{\partial t}\\
\end{split}
\label{P4eq7.23a}
\end{equation}
where $(\rho C)_{\rm sn}$ is the snow volumetric heat capacity (J ${\rm m}^{-3}{\rm K}^{-1}$), $D_{\rm sn}$ is the snowpack depth (m), $T_{\rm sn}$ is the snow temperature (K), and the energy fluxes $R^{\rm N}_{\rm sn}$, $H_{\rm sn}$ and $G^{\rm B}_{\rm sn}$ are the net radiation (shortwave and longwave), sensible heat flux and basal heat flux (W ${\rm m}^{-2}$), respectively. The mass fluxes $E_{\rm sn}$ and $M_{\rm sn}$ are the snow sublimation and melting (kg ${\rm m}^{-2}{\rm s}^{-1}$), respectively, that are associated with the latent heat of sublimation $L_{\rm s}$ and fusion $L_{\rm f}~({\rm J~kg^{-1}})$. The superscript $INT$ denotes internal phase changes, where $Q^{\rm INT}_{\rm sn}$ is the heat change associated with internal phase changes, and $S_{\rm l}$ the snow liquid water content (SLW) $({\rm kg~m^{-2}})$. Without loss of generality, it can be assumed that for the grid box characteristic of NWP models, the following expression is valid
\begin{equation}%eq7.23b
S_{\rm l}=S_{\rm l}(T_{\rm sn},S)\approx f(T_{\rm sn})S^c_l(S,\rho_{\rm sn})
\label{P4eq7.23b}
\end{equation}
where $S^{\rm c}_{\rm l}~({\rm {kg~m^{-2}}})$ is the snow liquid water capacity, $S$ is the sum of snow and water in the snowpack (also referred as SWE) and $\rho_{\rm sn}$ is the snow density $\rm{(kg~m^{\rm -3}}$). The snow temperature function is prescribed in an analytical form -- following a similar approach described by \citet{Viterboetal1999} for soil phase changes
\begin{equation}%eq7.23c
f(T_{\rm sn}) =\left\{
\begin{array}{l l}
0  & \quad  T_{\rm sn} < T_{\rm f}-d/2 \\
1+\sin \biggl\{ \dfrac{\pi(T_{\rm sn}-T_{\rm f})}{d}  \biggr\}   & \quad  T_{\rm sn} \geq T_{\rm f}-d/2 \\
\end{array} \right. \}
\label{P4eq7.23c}
\end{equation}
where $T_{\rm f}$ is the triple-point temperature for water and $d$ is a characteristic temperature difference, with respect to $T_{\rm f}$, limiting the phase change regime. In the numerical implementation, $d=4~{\rm K}$ was chosen. Snow liquid water capacity is approximated as a function of $S$ and $\rho_{\rm sn}$, following \citet{Anderson1976} 
\begin{equation}%eq7.23d
S^c_{\rm l}=S[r_{\rm l,min} + ( r_{\rm l,max} - r_{\rm l,min} ) 
      \max(0,\rho_{\rm sn,l}-\rho_{\rm sn} )/\rho_{\rm sn,l} ]
\label{P4eq7.23d}
\end{equation}
with the constants $r_{\rm l,min}=$0.03, $r_{\rm l,max}=$0.1 and $\rho_{\rm sn,l}=200~{\rm kg m^{-3}}$. This equation is a simple parametrization of a very complex phenomenon and has been used recently in other snow schemes, for example, \citet{BooneAndEtchevers2001}.

Combining (\ref{P4eq7.23a}) and (\ref{P4eq7.23b}) results in a modified snow energy budget equation,
\begin{equation}%eq7.23
\biggl[ (\rho C)_{\rm sn} D_{\rm sn} + 
         L_{\rm f}S^{\rm c_l}\frac{\partial f(T_{\rm sn})}{\partial T_{\rm sn}} \biggr]
         \frac{\partial T_{\rm sn}}{\partial t} =  
         R^{\rm N}_{\rm sn}+L_sE_{\rm sn}+H_{\rm sn}-G^{\rm B}_{\rm sn}-L_{\rm f}M_{\rm sn}
\label{P4eq7.23}
\end{equation}

with one extra term in the lhs of the equation, that can be interpreted as an additional snow heat capacity -- or heat capacity barrier. In compacted snowpacks, the representation of SLW as a diagnostic increases the snow heat capacity by a factor of up to five. This increase acts as a heat barrier near $T_{\rm f}$, representing the increased snow temperature inertial due to freeze-melt events. The snow heat capacity in (\ref{P4eq7.23a}),(\ref{P4eq7.23}) is approximated as
\begin{equation}%eq7.24
(\rho C)_{\rm sn} \approx \frac{(\rho C)_{\rm i}}{\rho_{\rm i}}{\rho_{\rm sn}}
\label{P4eq7.24}
\end{equation}

This diagnostic approach for SLW also allows the representation of rainfall interception. The snow mass balance reads as
\begin{equation}%eq7.20
\frac{\partial S}{\partial t} = F + c_{\rm sn}F_{\rm l} +
   c_{\rm sn}E_{\rm sn} - R_{\rm sn},   
 \label{P4eq7.20}
\end{equation}
where $F$, $F_{\rm l}$, and $R_{\rm sn}$ are the mass fluxes of snowfall, rainfall and runoff $({\rm kg~m^{-2}~s^{-1}})$ and $c_{\rm sn}$ is the snow cover fraction. Rainfall is considered to reach the snowpack at $T_{\rm f}$, and latent heat released by the freezing of the intercepted rainfall, if $T_{\rm sn}<T_{\rm f}$, is also accounted in the energy-balance solution. Runoff is defined as the rate at which liquid water leaves the snowpack and parametrized as follows:
\begin{equation}%eq7.25a
R_{\rm sn}=c_{\rm sn}M_{\rm sn}+\max \biggl( c_{\rm sn}F_{\rm l} - \frac{S^{\rm c}_{\rm l}(1-f(T_{\rm sn}))}{\Delta \it{t}},0 \biggr)
\label{P4eq7.25a}
\end{equation}

Liquid water is generated by melting ($M_{\rm sn}$) and by rainfall interception ($F_{\rm l}$). When snow liquid water content exceeds the snow liquid water capacity (defined in (\ref{P4eq7.23c})) runoff is generated. 

In this section, all fluxes are per unit area and apply only to the snow area (i.e. tile 5 and 7). The snow mass $S$ applies to the entire grid square, as well as the snowfall flux from the atmospheric model. As a general rule, all quantities with subscript $sn$ refer to the snow area. Snow evaporation (\ref{P4eq7.23}, \ref{P4eq7.20}) is defined as 
\begin{equation}%eq7.21
c_{\rm sn} E_{\rm sn} = c_{5}E_{5} + c_{7} E_{{\rm sn},7}
\label{P4eq7.21}
\end{equation}

Snow mass and snow depth are related by
\begin{equation}%eq7.22
D_{\rm sn} = \frac{1}{\rho_{\rm sn}}\frac{S}{c_{\rm sn}}
 \label{P4eq7.22}
\end{equation}
where $D_{\rm sn}$ is snow depth for the snow-covered area (units m;
note that $D_{\rm sn}$ is not a grid-averaged quantity) and $\rho_{\rm sn}$ is
the snow density (units kg ${\rm m}^{-3}$).

\subsubsection{Prognostic snow density and albedo}\label{P4S7.4.2}

Snow density is assumed to be constant with depth, and the rate of density change is parametrized as
\begin{equation}%eq7.26a
\frac{1}{\rho_{\rm sn}}\frac{\partial \rho_{\rm sn}}{\partial t}=
\frac{\sigma_{\rm sn}}{\eta_{\rm sn}(T_{\rm sn},\rho_{\rm sn})}+
\xi_{\rm sn}(T_{\rm sn},\rho_{\rm sn})+
\frac{\max(0,Q^{\rm INT}_{\rm sn})}{L_f(S-S_l)}
\label{P4eq7.26a}
\end{equation}
where the first two terms in (\ref{P4eq7.26a}) represent overburden and thermal metamorphism \citep{Anderson1976,BooneAndEtchevers2001}, respectively, and the last term represents the compaction related to meltwater retained in the snowpack, adapted from \citet{Lynchstieglitz1994}. In the overburden term (first term on the rhs of (\ref{P4eq7.26a})), $\sigma_{\rm sn}$ and $\eta_{\rm sn}$ are the pressure of the overlying snow (Pa) and snow viscosity (Pa~s), respectively. Melted water retained in the snowpack leads to a decrease of snow depth while keeping $S$ constant. In snowfall conditions snow density is updated as the ratio between the new snow mass and new snow depth using the snowfall mass ($\Delta {\it tF}$) and density ($\rho_{\rm new}$)
\begin{equation}%eq7.26
\rho^{*}_{\rm sn} = 
\frac{S+\Delta {\it tF}}
{\frac{S}{\rho^{t}_{\rm sn}}+\frac{\Delta {\it tF}}{\rho_{\rm new}}}
\label{P4eq7.26}
\end{equation}
where $\rho^{*}_{\rm sn}$ is an updated snow density. Snowfall density $\rho_{\rm new}$ is given by an expression from CROCUS \citep{Brunetal1989,Brunetal1992} were fresh snow density is a function of near surface air temperature and wind speed
\begin{equation}%eq7.26b
\rho_{\rm new}=a_{\rm sn}+b_{\rm sn}(T_{\rm air}-T_{\rm f})+
  c_{\rm sn}(V_{\rm a})^{1/2}
\label{P4eq7.26b}
\end{equation}
where $T_{\rm air}$ and $V_{\rm a}$ are the near surface air temperature (K) and wind speed $({\rm m s^{-1}})$, respectively. The coefficients are $a_{\rm sn}=109~{\rm {kg~m^{-3}}}$, $b_{\rm sn}=6~{\rm {kg~m^{-3} K^{-1}}}$, and $c_{\rm sn}=26~{\rm {kg~m^{-7/2} s^{1/2}}}$. Snow density is constrained to be between 50 to 450 ${\rm {kg~m^{-3}}}$. 

The snow viscosity in (\ref{P4eq7.26a}) is formulated following \citet{Anderson1976} 
\begin{equation}%eq7.26c
\eta_{\rm sn}=\eta_0\exp\biggl(a_{\rm \eta}(T_{\rm f}-T_{\rm sn})+
       b_{\rm \eta}\rho_{\rm sn} \biggr)
\label{P4eq7.26c}
\end{equation}
where $\eta_0={\rm {3.7\times10^7~(Pa~s)}}$, $a_{\rm \eta}={\rm{8.1\times10^{-2}~{K^{-1}}}}$ and $b_{\rm \eta}={\rm{1.8\times10^{-2}~{m^3 kg^{-1}}}}$. The pressure of the overlying snow is given by $\sigma_{\rm sn}=\frac{1}{2}Sg$, where $g$ is the standard gravity $\rm{(m^{2}s^{-2})}$.

The thermal metamorphism (second term in the rhs of (\ref{P4eq7.26a})) is parametrized as
\begin{equation}%eq7.26d
\xi_{\rm sn}=a_{\rm \xi}\exp\biggl[ -b_{\rm \xi}(T_f-T_{\rm sn}) - 
     c_{\rm \xi}\max(0,\rho_{\rm sn}-\rho_{\rm \xi})  \biggr]
\label{P4eq7.26d}
\end{equation}
using the constant values of \citet{Anderson1976}: $a_{\rm \xi}={\rm{2.8\times10^{-6}~s^{-1}}}$, $b_{\rm \xi}={\rm {4.2\times10^{-2}}}$, $c_{\rm \xi}={\rm {460~m^3~kg^{-1}}}$, and $\rho_{\rm \xi}={\rm {150~kg~m^{-3}}}$.

Snow albedo in exposed areas evolves according to the formulation
of~\citet{Bakeretal1990},~\citet{Verseghy1991},~\citet{Douvilleetal1995} and \citet{Dutraetal2010}. For melting and non melting-conditions:
\begin{equation}%eq7.28
\alpha^{\rm t+1}_{\rm sn} = \left\{
\begin{array}{l l}
\alpha^{\rm t}_{\rm sn} -\tau_a \Delta t / \tau_1 
& \quad  M_{\rm sn}=0 \\
(\alpha^{\rm t}_{\rm sn}-\alpha_{\rm min})\exp(-\tau_f \Delta t / \tau_1)+\alpha_{\rm min}
& \quad  M_{\rm sn}>0 {\rm or} T_{\rm sn} \geq T_f-d/2 \\
\end{array} \right.
\label{P4eq7.28}
\end{equation}
where $\tau_{a}=0.008$, which will decrease the albedo by 0.1 in
12.5~days, $\alpha_{ min}=0.5$ and $\alpha_{ max}=0.85$. The timescales $\tau_{1}=86400$ s, and $\tau_{ f} = 0.24$ corresponding to an e-folding time of about 4 days (see Table \ref{P4T7.4}). 

A continuous reset is applied to snow albedo after snowfall events:
\begin{equation}%eq7.29
\alpha^{\rm t+1}_{\rm sn}=
\alpha^{\rm t}_{\rm sn}+\min\biggl( 1,\frac{F\Delta t}{10} \biggr)
(\alpha_{\rm max}-\alpha^{\rm t}_{\rm sn})
 \label{P4eq7.29}
\end{equation}
This formulation assumes that $10~{\rm {kg~m^{-2}}}$ of fresh snowfall are needed to reset the snow albedo to its maximum value ($\alpha_{\max}=0.85$).

The above formulae are inadequate to describe the evolution of the
surface albedo of snow cover with high vegetation. Observations
suggest a dependence on forest type but, by and large, the albedo
changes from a value around 0.3 just after a heavy snowfall to a
value around 0.2 after a few days (see~\citealp{BettsandBall1997} and
the discussion in~\citealp{ViterboandBetts1999}). This change reflects
the disappearance of intercepted snow, due to melt (for sufficiently
warm temperatures) or wind drift (for cold temperatures). Ways of
describing those two mechanisms would involve either a separate
albedo variable for the snow in the presence of high vegetation, or
the introduction of an interception reservoir for snow. In the
absence of any of the two, a vegetation-type-dependent albedo adapted from \citet{Moodyetal2007} is used (see Table \ref{P4T7.4a}). \citet{Moodyetal2007} provide a five-year (2000--2004) climatological statistic of northern hemisphere broadband (0.--5.0 $\mu$m) white-sky albedo for the 16 IGBP ecosystem classes when accompanied by the presence of snow on the ground. 

%T7.4
\begin{table}
 \centering
\caption{Snow-related parameters.\label{P4T7.4}}
\begin{tabular}{lp{20pc}l}
\toprule
\multicolumn{1}{c}{\bf Symbol} & \multicolumn{1}{c}{\bf Parameter} & \multicolumn{1}{c}{\bf Value} \\
%\midrule
\hline \hline
$D^{\max}_{\rm sn}$ & Maximum snow thermal depth & 1.00~m \\
$D_{\rm cr}$ & Threshold value for grid box coverage of snow & 0.1~m \\
$\alpha_{\min}$ & Minimum albedo of exposed snow & 0.50 \\
$\alpha_{\max}$ & Maximum albedo of exposed snow & 0.85 \\
$\lambda_{\rm i}$ & Ice heat conductivity & 2.2 $\text{W~m}^{-1}{\rm
K}^{-1}$ \\
$\rho_{\rm i}$ & Ice density & 920 $\text{kg~m}^{-3}$ \\
$(\rho C)_{\rm i}$ & Ice volumetric heat capacity & $2.05 \times10^{6}~\text{J~m}^{-3}{\rm K}^{-1}$ \\
$\tau_{\rm a}$ & Linear coefficient for decrease of albedo of non-melting snow & 0.008 \\
$\tau_{\rm f}$ & Coefficient for exponential decrease of snow density and melting snow albedo & 0.24 \\
$\tau_{1}$ & Length of day & 86400 s \\
\bottomrule
\end{tabular}
\end{table}

%T7.4a
\begin{table}
\centering
\caption{
Mean values of northern hemisphere broadband surface albedo (in presence of snow) aggregated by high vegetation type.
\label{P4T7.4a}}
\small
\begin{tabular}{llc}
\toprule
{\bf Index} & {\bf Vegetation type} & {\bf Albedo} \\
%\midrule
\hline \hline
3 & Evergreen needle leaf trees & 0.27 \\
4 & Deciduous needle leaf trees & 0.33 \\
5 & Deciduous broad leaf trees  & 0.31 \\
6 & Evergreen broad leaf trees  & 0.38 \\
18 & Mixed forest / woodland & 0.29 \\
19 & Interrupted forest & 0.29 \\
\bottomrule
\end{tabular}
\end{table}


\subsubsection{Additional details}\label{P4S7.4.3}

\textbf{Limiting of snow depth in the snow energy equation}\\
Initial experimentation with the snow model revealed that the time
evolution of snow temperature was very slow over Antarctica. The
reason is rather obvious; the snow depth over Antarctica is set to
a climatological value of 10 m which can respond only very slowly
to the atmospheric forcing due to its large thermal inertia. In
previous model versions, the properties of layer 1 were replaced
by snow properties when snow was present, which kept the timescale
short. A limit is put on the depth of the snow layer in the thermal budget,
$D^{\max}_{\rm sn} = 1.0 \text{ m}$. The energy equation reads
\begin{equation}%eq7.30
\begin{split}
\biggl[ (\rho C)_{\rm sn}D^{*}_{\rm sn} + 
         L_fS^c_l\frac{\partial f(T_{\rm sn})}{\partial T_{\rm sn}} \biggr]
         \frac{\partial T_{\rm sn}}{\partial t} &=  
         R^{\rm N}_{\rm sn}+L_sE_{\rm sn}+H_{\rm sn}-G^{\rm B}_{\rm sn}-L_fM_{\rm sn} \\
D^{*}_{\rm sn} &= \min (D_{\rm sn}, D^{\max}_{\rm sn})
\end{split}
 \label{P4eq7.30}
\end{equation}
A physical solution will be explored in a future scheme with a
multilayer snow model, with for example four layers to represent
timescales from one day to a full annual cycle.

\textbf{Basal heat flux and thermal coefficients}

The heat flux at the bottom of the snow pack is written as a finite difference using
\begin{equation}%eq7.31
G^{\rm B}_{\rm sn} = \frac{T_{\rm sn}-T_{1}}{r_{\rm sn}}
 \label{P4eq7.31}
\end{equation}
where $r_{\rm sn}$ is the resistance between the middle of the snow
pack and the middle of soil layer 1, with two components: the
resistance of the lower part of the snow pack and the resistance of
the top half of soil layer 1. Therefore
\begin{equation}%eq7.32
r_{\rm sn} = 0.5 \frac{D^{*}_{\rm sn}}{\lambda_{\rm
sn}}+\frac{1}{\Lambda_{\rm sk,8}}
 \label{P4eq7.32}
\end{equation}
where the second term is the skin layer conductivity for bare soil
(tile 8), which can be seen as an approximation of
$0.5(D_{1}/\lambda_{\rm T})$. The snow thermal conductivity, is
related to the ice thermal conductivity according
to~\citet{Douvilleetal1995} given by
\begin{equation}%eq7.33
\lambda_{\rm sn} = \lambda_{\rm i}\biggl(\frac{\rho_{\rm
sn}}{\rho_{\rm i}}\biggr)^{\!1.88}
 \label{P4eq7.33}
\end{equation}
Table \ref{P4T7.4} contains the numerical values of the ice
density and ice heat conductivity.

\textbf{Numerical solution for non-melting situations}

The net heat flux that goes into the top of the snow deck is an output of the vertical diffusion scheme
\begin{equation}%eq7.34
H^{\rm N}_{\rm sn} = R^{\rm N}_{\rm sn} - L_{\rm s} E_{\rm
sn}-H_{\rm sn}
 \label{P4eq7.34}
\end{equation}
In the absence of melting, the solution of (\ref{P4eq7.30}) is
done implicitly. The preliminary snow temperature, prior to the
checking for melting conditions, $T^{*}_{\rm sn}$, is given by
\begin{equation}%eq7.35
A_{1} \frac{T^{*}_{\rm sn}-T^{t}_{\rm sn}}{\Delta t} =H^{\rm N}_{\rm
sn} - \frac{T^{*}_{\rm sn}-T_{1}}{r_{\rm sn}}
 \label{P4eq7.35}
\end{equation}
\begin{equation}%eq7.36
\begin{split}
A_{1} &= \min\biggl[\frac{(\rho C)_{\rm i}}{\rho_{\rm
i}c_{\rm sn}}S, A^{\max}_{1}\biggr]+
{\left. {L_{\rm f}S^{\rm c}_{\rm l}\frac{\partial f(T_{\rm sn})}{\partial T_{\rm sn}}} \right|_{T^t_{\rm sn}}} \\
A^{\max}_{1} &= \frac{(\rho C)_{\rm i}}{\rho_{\rm i}}\rho_{\rm sn}D^{\max}_{\rm sn}
 \label{P4eq7.36}
\end{split}
\end{equation}
where superscript $t$ refers to the current time step and superscript
$^*$ to the preliminary value at the next time step. The solution for
$T^{*}_{\rm sn}$ is obtained from
\begin{equation}%eq7.37
T^{*}_{\rm sn} \biggl(1+\frac{\Delta t}{r_{\rm sn}A_{1}}\biggr) =
T^{t}_{\rm sn} + \frac{\Delta t}{A_{1}}\biggl(H^{\rm N}_{\rm sn} +
\frac{T^{t}_{1}}{r_{\rm sn}}\biggr)
 \label{P4eq7.37}
\end{equation}
The basal snow heat flux to be used as input for the thermal budget
of the soil (in the snow covered fraction only) is
\begin{equation}%eq7.38
G^{\rm B}_{\rm sn} = \frac{T^{*}_{\rm sn}-T_{1}}{r_{\rm sn}}
 \label{P4eq7.38}
\end{equation}
Finally, a preliminary new value for the snow mass, $S^{*}$, is
computed from snow fall and snow evaporation using
\begin{equation}%eq7.39
\frac{S^{*}-S^{t}}{\Delta t} = F+c_{\rm sn}E_{\rm sn}
 \label{P4eq7.39}
\end{equation}

\subsubsection{Treatment of melting}\label{P4S7.4.4}

\subsubsubsection{No melting occurs}

If $T^{*}_{\rm sn}<T_{0}$ no melting occurs and the preliminary
values $T^{*}_{\rm sn}$ and $S^{*}$ become the $t+1$ values, while
the basal heat flux is given by (\ref{P4eq7.38}).

\subsubsubsection{Melting conditions}

If $T^{*}_{\rm sn}>T_{0}$, snow melting occurs and the time step  is
divided in two fractions, $\Delta t = \Delta_{1}t+\Delta_{2}t$, where
the first fraction, $\Delta_{1}t$ brings the temperature to $T_{0}$
with no melting so that
\begin{equation}%eq7.40
\Delta_{1}t = \frac{A_{1}(T_{0}-T^{t}_{\rm sn})}{H^{\rm N}_{\rm sn}-
(T_{0}-T^{t}_{1})/r_{\rm sn}}
 \label{P4eq7.40}
\end{equation}
while, during the second fraction, $\Delta_{2}t$, melting occurs with
no resultant warming of the snow so that
\begin{equation}%eq7.41
\begin{split}
&T^{t+1} = T_{0} \\
&Q_{\rm sn} = H^{\rm N}_{\rm sn} - G^{\rm B}_{\rm sn} \\
&\frac{S^{t+1}-S^{*}}{\Delta_{2}t} = -c_{\rm sn} M_{\rm
sn} = -c_{\rm sn} \frac{Q_{\rm sn}}{L_{\rm f}} = -c_{\rm
sn}\frac{H^{N}_{\rm sn}-G^{\rm B}_{\rm sn}}{L_{\rm f}}
\end{split}
 \label{P4eq7.41}
\end{equation}
{\it If not all the snow melts}, i.e., if $S^{t+1} >0$, the heat flux
passed to the soil is
\begin{equation}%eq7.42
G^{\rm B}_{\rm sn} = \frac{T_{0}-T_{1}}{r_{\rm sn}}
 \label{P4eq7.42}
\end{equation}
{\it When all the snow melts}, i.e., if $S^{t+1} <0$, the melting time step is redefined as
\begin{equation}%eq7.43
\begin{split}
S^{t+1} & = 0 \\
\Delta_{2} t & = L_{\rm f}\frac{S^{*}}{c_{\rm
sn}(H^{\rm N}_{\rm sn}-G^{\rm B}_{\rm sn})} \\
\Delta_3 t &= 1-(\Delta_{1}t + \Delta_{2} t)
\end{split}
 \label{P4eq7.43}
\end{equation}
and the basal heat flux is redefined as
\begin{equation}%eq7.44
G^{\rm B}_{\rm sn} = \frac{\Delta_{1}t + \Delta_{2}t}{\Delta
t}\frac{T_{0}-T^{t}_{1}}{r_{\rm sn}}+ \frac{\Delta_{3}t}{\Delta
t}H^{\rm N}_{\rm sn}
 \label{P4eq7.44}
\end{equation}

\subsubsection{Numerical stability of the snow scheme}
Considering the snow scheme is adopting an explicit formulation occasional numerical instability can occur with long timesteps. 
This is particularly the case since the conductivity between the forest canopy and the underlying snow (skin layer conductivity for tile 7) 
has increased (since the IFS cycle 39r1). The snow tiles are represented in figure \ref{P4F7.3n}.

%F7.3 new snow scheme representation
\begin{figure}[h!]
\centering
\includegraphics[width=120mm]{Figures/figA2.pdf}
\caption{Schematic representation of the snow tiles. \label{P4F7.3n}}
\end{figure}

A standard solution for stability problems would imply the coupled formulation to be fully implicit. 
However a simpler solution that does not require the heat flux $G$ and the snow temperature $T_{sn}^{n+1}$ at time level $n+1$ was found. 
This involves an assumed snow temperature where $T_{sn}$ respond to heat flux $G$ and it consider a simplified version of the snow energy equation \ref {P4eq7.30} with the linear relation
%
\begin{eqnarray}
   \frac{\partial T_{sn}}{\partial t} = A_{sn} + B_{sn} G \;.\label{eqn_ABsn}
\end{eqnarray}

In this case the snow heat flux is approximated as
%
\begin{eqnarray}
   G &=& \Lambda_{sk} [T_{sk}^{n+1} - T_{sn}^{n+1}]   \nonumber  \\
     &=& \Lambda_{sk} [T_{sk}^{n+1} - T_{sn}^{n} 
                        -(A_{sn}+B_{sn} G)\Delta t] \;,  \nonumber  \\
   G &=& \frac{\Lambda_{sk}}{1+\Lambda_{sk} B_{sn} \Delta t} 
            [T_{sk}^{n+1} - T_{sn}^{n} - A_{sn} \Delta t] \;. \label{eqn_G2}
\end{eqnarray}
%
The coefficients $A_{sn}$ and $B_{sn}$ can only be calculated correctly by the snow scheme, so a call to the snow scheme would be necessary before the surface
energy balance is called involving a substantial redesign effort. A simplified estimate of $A_{sn}$ and $B_{sn}$ is used instead with $A_{sn}=0$ 
and that $B_{sn}=(D(\rho C)_{sn})^{-1}$ as it was found efficient to prevent oscillatory evolution of the snow temperature.
A further security is introduced in the snow temperature absolute values which are not allowed to be below $180 K$.

For introducing the numerical stability term, it is sufficient to replace $\Lambda_{{\rm sk},i}$ in the surface energy balance equation \ref{P4eq7.19} of tile 5 and tile 7 
with the $\Lambda_{{\rm sks},i}$ given by the formulation below.
%
\begin{eqnarray}
  \Lambda_{sks} \longrightarrow 
     \frac{\Lambda_{sk}}{1+\Lambda_{sk} B_{sn} \Delta t} 
\label{eqn_L2}
\end{eqnarray}
%
Equation (\ref{eqn_L2}) enhances the original $\Lambda_{sk}$ to guarantee numerical stability and it mimics a fully implicit formulation provided the approximated estimate of $A_{sn}$ and $B_{sn}$.

\subsection{Soil heat transfer}\label{P4S7.5}

In the absence of internal phase changes, the soil heat transfer is
assumed to obey the following Fourier law of diffusion
\begin{equation}%eq7.45
(\rho C)_{\rm soil} \frac{\partial T}{\partial t} =
\frac{\partial}{\partial z}\biggl[\lambda_{T}\frac{\partial
T}{\partial z}\biggr]
 \label{P4eq7.45}
\end{equation}
where $(\rho C)_{\rm soil}$ is the volumetric soil heat capacity
$(\text{J m}^{-3}{\rm K}^{-1})$, $T$ is the soil temperature (units
K), $z$ is the vertical coordinate (i.e. the distance from the surface,
positive downwards -- units m), and $\lambda_{\rm T}$ is the thermal
conductivity $(\text{W m}^{-1}{\rm K}^{-1})$. The above equation
assumes that heat fluxes are predominantly in the vertical direction,
that the effects of phase changes in the soil and the heat transfer
associated with the vertical movement of water in the soil can be
neglected~\citep{deVries1975}, and that the effects of hysteresis can
be neglected~\citep{Milly1982}.

The boundary condition at the bottom, no heat flux of energy, is
an acceptable approximation provided that the total soil depth is
large enough for the time-scales represented by the model or, in
other words, the bottom of the soil is specified at a depth where
the amplitude of the soil heat wave is a negligible fraction of
its surface amplitude (see~\citet{deVries1975} and next
\ref{P4S7.6}).

\subsubsection{Discretization and choice of
parameters}\label{P4S7.5.1}


%T7.5
\begin{table}
 \centering
\caption{Parameters in the land-surface scheme. See
\hyperTable{P4T7.4} for snow-related parameters.\label{P4T7.5}}
\begin{tabular}{lll}
\toprule
\multicolumn{1}{c}{\bf Symbol} & \multicolumn{1}{c}{\bf Parameter} & \multicolumn{1}{c}{\bf Value} \\
%\midrule
\hline \hline
$b_I$ & Interception efficiency & 0.25 \\
$D_{1}$ & Depth of soil layer 1 & 0.07 m \\
$D_{2}$ & Depth of soil layer 2 & 0.21 m \\
$D_{3}$ & Depth of soil layer 3 & 0.72 m \\
$D_{4}$ & Depth of soil layer 4 & 1.89 m \\
$F_{\rm cv}$ & Fraction of gridbox covered by convective rainfall & 0.5 \\
$k$ & Heterogeneity factor for convective precipitation & 0.5 \\
$T_{\rm f1}$ & Highest temperature for existence of ice water &
$T_{0}+1$ \\
$T_{\rm f2}$ & Lowest temperature for existence of liquid water & $T_{0}-3$ \\
$w_{1\max}$ & Maximum water amount on single leaf & 0.0002~m \\
%$\gamma_{\rm sat}$ & Hydraulic conductivity at saturation  & $4.57
%\times 10^{-4} \text{ m s}^{-1}$ \\
%$b$ & Clapp and Hornberger soil parameter & 6.04 \\
%$\lambda_{\rm dry}$ & Heat conductivity of dry soil & 0.190~W ${\rm m}^{-1}{\rm K}^{-1}$ \\
%$\lambda_{\rm sm}$ & Heat conductivity of soil matrix & 3.44~W ${\rm m}^{-1}{\rm K}^{-1}$ \\
%$\lambda_{\rm w}$ & Heat conductivity of liquid water & 0.57~W ${\rm m}^{-1}{\rm K}^{-1}$ \\
%$(\rho C)_{\rm soil}$ & Volumetric soil heat capacity & $2.19 \times
%10^{6} \text{ J m}^{-3}{\rm K}^{-1}$ \\
%$\theta_{\rm sat}$ & Soil moisture at saturation  & 0.472~${\rm m}^{3}{\rm m}^{-3}$ \\
%$\theta_{\rm cap}$ & Soil moisture at field capacity & 0.323~${\rm m}^{3}{\rm m}^{-3}$ \\
%$\theta_{\rm pwp}$ & Soil moisture at permanent wilting point & 0.171~${\rm m}^{3}{\rm m}^{-3}$ \\
%$\psi_{\rm sat}$ & Matric potential at saturation & $-0.338$~m \\
\bottomrule
\end{tabular}
\end{table}

For the solution of (\ref{P4eq7.45}) the soil is discretized in
four layers, of depths $D_{k}$ $(k = 1,2,3,4)$, the temperatures
are defined at full layers $(T_{k})$, and the heat fluxes, at half
layers ($\hat{G}_{k+1/2}$ is the heat flux, positive downwards,
units $\text{W m}^{-1}$, at the interface between layer $k$ and
$k+1$). An energy-conserving implicit algorithm is used, leading
to a tridiagonal system of equations with solution detailed in
\hyperSection{P4S7.9}.

The boundary condition at the bottom is
\begin{equation}%eq7.46
G_{4+1/2} = 0
 \label{P4eq7.46}
\end{equation}
At the top, the boundary condition is the soil heat flux at the
surface, computed as a weighted average over the tiles. For the
snow free tiles, the flux into the soil consists of two parts.
Apart from the diffusion of heat governed by $\Lambda_{{\rm
sk},i}(T_{{\rm sk},i}-T_{1})$ (see (\ref{P4eq7.19})), the net
short-wave radiation not absorbed by the skin layer ($f_{{\rm
Rs},i}$) provides energy to the soil. Table \ref{P4T7.2} lists the
values of $\Lambda_{{\rm sk},i}$ and $f_{{\rm Rs},i}$ for each of
the tiles. For the snow tiles, the heat flux into the soil is the
snow basal flux, calculated using a resistance formulation and
modified in the case of partial melting (see (\ref{P4eq7.31}),
(\ref{P4eq7.38}), (\ref{P4eq7.42}), and (\ref{P4eq7.44})).

The net heat flux into the soil is given by
\begin{equation}%eq7.47
G_{1/2} = \sum_{i}C_{i}[\Lambda_{{\rm sk},i}(T_{{\rm
sk},i}-T_{1})+f_{{\rm Rs},i}(1-\alpha_{i})R_{\rm s}]+c_{\rm
s}G^{\rm B}_{\rm sn} \label{P4eq7.47}
\end{equation}
where the summation scans all snow free tiles.

The volumetric soil heat capacity is assumed constant, with value
$2.19 \times 10^6 \text{ J m}^{-3}{\rm K}^{-1}$ (see
\hyperTable{P4T7.5} for a list of constants used by the model). The
heat conductivity, $\lambda$, depends on the soil-water content
following~\citet{PetersLidardetal1998}~\citep[see
also][]{Farouki1986,Johansen1975} and is given by a combination of
dry $\lambda_{\rm dry}$ and saturated $\lambda_{\rm sat}$ values,
weighted by a factor known as the Kersten number, $K_{\rm e}$, so
that
\begin{equation}%eq7.48
\lambda = K_{\rm e}(\lambda_{\rm sat} - \lambda_{\rm
dry})+\lambda_{\rm dry} \label{P4eq7.48}
\end{equation}
where $\lambda_{\rm dry} = 0.190~{\rm W}~\text{m}^{-1}~{\rm K}^{-1}$ and
\begin{equation}%eq7.49
\lambda_{\rm sat} = \lambda^{1-\theta_{\rm sat}}_{\rm sm}
\lambda_{\rm w}^{\theta} \label{P4eq7.49}
\end{equation}
where the heat conductivity of the soil matrix, $\lambda_{\rm sm}
= 3.44~{\rm W}~{\rm m}^{-1}~{\rm K}^{-1}$ and the thermal
conductivity of water is $\lambda_{\rm w} = 0.57~{\rm W}~{\rm
m}^{-1}~{\rm K}^{-1}$. Equation (\ref{P4eq7.49}) represents a
simplification of Peters-Lidard formulation, neglecting the
changes in conductivity due to ice water and assuming the quartz
content typical of a loamy soil. Finally, the Kersten number for
fine soils was selected in~\citet{PetersLidardetal1998}
as
\begin{equation}%eq7.50
K_{e} = \log_{10}\biggl[\max\biggl(0.1,\frac{\theta}{\theta_{\rm
sat}}\biggr)\biggr]+1
 \label{P4eq7.50}
\end{equation}
The depths of the soil layers are chosen in an approximate
geometric relation (see Table \ref{P4T7.5}), as suggested in
\cite{Deardorff1978}. \citet{Warrilowetal1986} have shown that four
layers are enough for representing correctly all timescales from
one day to one year.  Using the numerical values of the heat
capacity and soil depths defined in Table \ref{P4T7.5}, the
amplitude and phase response of the numerical solution of
(\ref{P4eq7.45}) were analysed by~\citet{ViterboandBeljaars1995}
for typical values of soil moisture in (\ref{P4eq7.48}), and for
harmonic forcings at the surface with periods ranging from half a
day to two years. The analysis points to an error in the numerical
solution of less than 20\% in amplitude and 5\% in phase for
forcing periods between one day and one year.

\subsubsection{Soil-water phase changes}\label{P4S7.5.2}

At high and mid latitudes the phase changes of water in the soil have
an important effect on the water and energy transfer in the soil. A
proper consideration of the solid phase of soil water requires
modifications including, in order of importance, the following.
\begin{enumerate}
\item[(i)] The thermal effects related to the latent heat of
fusion/freezing~\citep[e.g.][]{Rouse1984}.

\item[(ii)] Changes in the soil thermal conductivity due to the
presence of ice (e.g.~\citealp{Penner1970}, not included in TESSEL as
mentioned in the previous section).

\item[(iii)] Suppression of transpiration in the presence of frozen
ground~\citep[e.g.][]{Bettsetal1998} and already described in
(\ref{P4eq7.11}).

\item[(iv)] Soil water transfer dependent on a soil water potential
including the effects of frozen water~\citep[e.g.][]{Lundin1989},
represented in a proxy way by (\ref{P4eq7.66}).
\end{enumerate}
The latent-heat effects are described in the following. The main
impact will be to delay the soil cooling in the beginning of the cold
period, and to delay the soil warming in spring, although the latter
effect is less important because it occurs when the solar forcing is
significant. Both effects make the soil temperatures less responsive
to the atmospheric forcing and damp the amplitude of the annual soil
temperature cycle. More details on the soil-freezing scheme and its
impact on forecasts and the model climate are described
in~\citet{Viterboetal1999}.

The soil energy equation, (\ref{P4eq7.45}), is modified in the
presence of soil water phase changes as
\begin{equation}%eq7.51
(\rho C)_{\rm soil}\frac{\partial T}{\partial t} =
\frac{\partial}{\partial z}\biggl[\lambda_{T}\frac{\partial
T}{\partial z}\biggr]+L_{\rm fus}\rho_{\rm w}\frac{\partial \theta_{\rm
I}}{\partial t}
 \label{P4eq7.51}
\end{equation}
where $\theta_{\rm I}$ is the volumetric ice-water content. Without
loss of generality, for the grid squares characteristic of NWP models
it can be assumed that
\begin{equation}%eq7.52
\theta_{\rm I} = \theta_{\rm I} (\theta, T) = f(T)\theta
 \label{P4eq7.52}
\end{equation}
where $\theta$ is the total soil-water content (liquid $+$ ice), and
\begin{equation}%eq7.53
\begin{array}{ll}
f_{\rm fr}(T)=0 & T>T_{\rm f1} \\
0<f_{\rm fr}(T)<1 & T_{f1}\leq T\leq T_{\rm f2} \\
f_{\rm fr}(T) =1 & T< T_{\rm f2}
\end{array}
 \label{P4eq7.53}
\end{equation}
where $T_{\rm f1}$ and $T_{\rm f2}$ are characteristic
temperatures limiting the phase change regime. In reality, the
values of $T_{\rm f1}$ and $T_{\rm f2}$ and the function $f_{\rm
fr}(T)$ have complicated dependencies on soil texture and
composition~\citep[see e.g.][]{WilliamsandSmith1989}, but here
they are approximated in a simple way. For an idealized
homogeneous, one-component soil, $f_{\rm fr}(T)$ would be a
step-function. According to~\cite{WilliamsandSmith1989} physical
reasons for having an interval over which melting/freezing is
active, rather than a threshold temperature, include the
following.
\begin{enumerate}
\item[(i)] Adsorption, resulting from forces between the mineral
parts of the soil and the water.

\item[(ii)] Capillarity, related to the fact that the water-free
surface is not plane.

\item[(iii)] Depression of the freezing point due to the effect of
dissolved salts.

\item[(iv)]      Soil heterogeneity.
\end{enumerate}
To avoid an undesirable coupling between the temperature and water
equations in the soil, (\ref{P4eq7.52}) is simplified to
\begin{equation}%eq7.54
\theta_{\rm I} = f_{\rm fr}(T)\theta_{\rm f}
 \label{P4eq7.54}
\end{equation}
where $\theta_{\rm f}$ is a constant, representing the amount of
soil water that can be frozen (thawed). For simplicity,
$\theta_{\rm f}= (c_{\rm H}+c_{\rm L})\theta_{\rm cap}$. The
scaling with the vegetated fractions is the simplest way of
distinguishing between dry (vegetation-sparse areas, e.g. deserts)
and wet (vegetated) areas. Combining (\ref{P4eq7.54}) with
(\ref{P4eq7.51}) results~in
\begin{equation}%eq7.55
\biggl[(\rho C)_{\rm soil} -L_{\rm fus}\rho_{\rm w}\frac{\partial
f_{\rm fr}}{\partial T}\biggr]\frac{\partial T}{\partial t} =
\frac{\partial}{\partial z}\biggl[\lambda_{T}\frac{\partial
T}{\partial z}\biggr]
 \label{P4eq7.55}
\end{equation}
showing that the effect of freezing can be interpreted as an
additional soil heat capacity, sometimes referred in the literature
as the `heat-capacity barrier' around freezing; not considering the
process of soil water freezing/melting can lead to very large
artificial temperature changes that do not
occur in nature when sufficient soil water is available.

Finally, function $f_{\rm fr}(T)$, is given by
\begin{equation}%eq7.56
f_{\rm fr}(T) =
\begin{cases}
0  & \quad  T > T_{{\rm f}1}\\
0.5 \biggl\{ 1-\sin \biggl[\dfrac{\pi (T-0.5 T_{{\rm f}1} -0.5
T_{{\rm f}2})}{T_{{\rm f}1} -T_{{\rm f}2}}\biggr]\biggr\}    & \quad  T_{{\rm f}2} \geq T \geq T_{{\rm f}1}\\
1  & \quad  T < T_{{\rm f}2}\\
\label{P4eq7.56}
\end{cases}
\end{equation}
with $T_{{\rm f}1} = T_0 + 1, T_{{\rm f}2} = T_0 - 3$.

\subsection{Soil-water budget}\label{P4S7.6}

The vertical movement of water in the unsaturated zone of the soil
matrix obeys the following equation
(see~\citet{Richards1931},~\citet{Philip1957},~\citet{Hillel1982}
and~\citet{Milly1982} for the conditions under which
(\ref{P4eq7.57}) and~(\ref{P4eq7.58}) are valid) for the
volumetric water content $\theta$:
\begin{equation}%eq7.57
\rho_{\rm w}\frac{\partial \theta}{\partial t} = -\frac{\partial
F_{\rm w}}{\partial z} +\rho_{\rm w}S_{\theta} \label{P4eq7.57}
\end{equation}
$\rho_{\rm w}$ is the water density (${\rm kg\ m}^{-3}), F_{\rm w}$
is the water flux in the soil (positive downwards, ${\rm kg\ m}^{-2}
{\rm s}^{-1}$), and $S_{\theta}$ is a volumetric sink term (${\rm
m}^3 {\rm m}^{-3} {\rm s}^{-1}$), corresponding to root extraction. Using
Darcy's law, $F_{\rm w}$ can be specified as
\begin{equation}%eq7.58
F_{\rm w} =\rho_{\rm w} \biggl(\lambda \frac{\partial
\theta}{\partial z} -\gamma \biggr) \label{P4eq7.58}
\end{equation}
$\lambda$ (m$^2$ s$^{-1}$) and $\gamma$ (m s$^{-1}$) are the
hydraulic diffusivity and hydraulic conductivity, respectively.

Replacing~(\ref{P4eq7.58}) in~(\ref{P4eq7.57}), specifying
$S_{\theta} =S_{\theta}(\theta, z)$, and defining parametric
relations for $\lambda$ and $\gamma$ as functions of soil water, a
partial differential equation for $\theta$ is obtained; it can be
numerically integrated if the top boundary condition is
precipitation minus evaporation minus surface runoff. The bottom
boundary condition assumes free
drainage.~\citet{Abramopoulosetal1988} specified free drainage or
no drainage, depending on a comparison of a specified geographical
distribution of bedrock depth, with a model-derived water-table
depth.  For the sake of simplicity the assumption of no bedrock
everywhere has been adopted.

\sububsection{Interception}\label{P4S7.6.1}

The interception reservoir is a thin layer on top of the
soil/vegetation, collecting liquid water by the interception of rain
and the collection of dew, and evaporating at the potential rate.
The water in the interception reservoir, $W_1$, obeys
\begin{equation}%eq7.59
\rho_{\rm w}\frac{\partial W_l}{\partial t} =c_1E_1 +D+I
\label{P4eq7.59}
\end{equation}
where $C_1E_1$ is the water evaporated by the interception
reservoir (or dew collection, depending on its sign), $D$
represents the dew deposition from other tiles, and $I$
(kg~m$^{-2}$~s$^{-1}$) is the interception---the fraction of
precipitation that is collected by the interception reservoir and
is later available for potential evaporation. Because the
interception reservoir has a very small capacity (a maximum of the
order of 1~mm, see (\ref{P4eq7.2})), it can fill up or evaporate
completely in one time step; special care has to be taken in order
to avoid numerical problems when integrating (\ref{P4eq7.59}). In
addition, since $E_1$ is defined in the vertical diffusion code,
it might impose a rate of evaporation that depletes entirely the
interception layer in one time step. In order to conserve water in
the atmosphere-intercepted water--soil continuum, the mismatch of
evaporation of tile 3 plus dew deposition from the other tiles
(which is not explicitly dealt with by the vertical diffusion) as
seen by the vertical diffusion and the intercepted water has to be
fed into the soil.

The equation is solved in three fractional steps: evaporation, dew
deposition, and rainfall interception. The solver provides the
following as outputs.
\begin{enumerate}
\item[(i)] Interception layer contents at time step $n +1, W_{\rm i}^{n +1}$.

\item[(ii)] Throughfall (i.e. rainfall minus intercepted water).

\item[(iii)] The evaporation effectively seen by the intercepted layer in each tile $i$.
\end{enumerate}
First, the upward evaporation ($E_1 <0$) contribution is
considered; because $C_1E_1$ depends linearly on $W_1$
(see~(\ref{P4eq7.2})), an implicit version of the evaporating part
of~(\ref{P4eq7.59}) is obtained by linearizing $C_1(W_1)E_1$
giving
\begin{equation}%eq7.60
\rho_{\rm w}\frac{W_1^* -W^t_1}{\Delta t} =C_1(W_1^t)E_1
+\frac{E_1}{W_{\rm 1m}}(W_1^* -W_1^t) \label{P4eq7.60}
\end{equation}
where $W_1^*$ is the new value of interception-reservoir content
after the evaporation process has been taken into account. After
solving for $W_1^*$, a non-negative value of evaporation is obtained
and the evaporation seen by this fractional time step is calculated
\begin{equation}%eq7.61
\begin{split}
W^1_1 &= \max (0, W_1^*)\\
E^1   &= \rho_{\rm w}\frac{W^1_1 -W^t_1}{\Delta t}
\label{P4eq7.61}
\end{split}
\end{equation}
The dew deposition is dealt with explicitly for each non-snow tile
in succession, for tiles 3, 4, 6, 7, 8, where tile 7 is also
considered because in the exposed snow tile, the canopy is in
direct evaporative contact with the atmosphere. When the
evaporative flux is downwards ($E_1 > 0$)
\begin{equation}%eq7.62
\begin{split}
W^2  &=W^1  +\min \biggl(W_{\rm 1m} -W^1, \frac{\Delta t}{\rho_{\rm
w}}c_iD_i\biggr)\\
D_i  &= \rho_{\rm w}\frac{W^2_{1, i}-W^1_1}{\Delta t}
\label{P4eq7.62}
\end{split}
\end{equation}
where superscript 2 denotes the final value at the end of the this
fractional time step.

The interception of rainfall is considered by applying the following
set of equations to large-scale and convective rainfall
\begin{equation}%eq7.63
\begin{split}
W^3_1  & = W^2_1 +\min \biggl(W_{\rm 1m} -W^2_1, \frac{\Delta
t}{\rho_{\rm w}}b_{\rm I} (c_{\rm H}+c_{\rm L})R_{\rm 1s}\biggr)\\
T_{\rm 1s}  &= R_{\rm 1s} -\rho_{\rm w}\frac{W^3_1 -W^2_1}{\Delta
t}\\
W^{t+1}_1  & = W^3_1 +\min \biggl(W_{\rm 1m} -W^3_1, \frac{\Delta
t}{\rho_{\rm w}}b_{\rm I} (c_{\rm H} + c_{\rm L})\frac{R_{\rm cv}}{F_{\rm cv}}\biggr)\\
T_{\rm cv}  &= R_{\rm cv} -\rho_{\rm w}\frac{W^{t+1}_{1}- W^3_1}{\Delta t}\\
\label{P4eq7.63}
\end{split}
\end{equation}
$R_{\rm cv}/F_{\rm cv}$ is a modified convective rainfall flux,
computed by applying the heterogeneity assumption that convective
rainfall only covers a fraction $F_{\rm cv}=0.5 $ of the grid box,
$b_{\rm I} =0.5$ is a coefficient of efficiency of interception of rain.
The total evaporation seen by the interception reservoir is $D_i$
for tiles 4, 6, 7, and 8 and $c_{\rm l}E_{\rm l} +D_i$ for tile 3.

The interception reservoir model described in this section is
probably the simplest water-conserving formulation based on
Rutter's original
proposition~\citep{Rutteretal1972,Rutteretal1975}. For more
complicated formulations still based on the Rutter concept see,
for
instance,~\citet{MahfoufandJacquemin1989},~\citet{DolmanandGregory1992}
and~\citet{Ridder2001}.

\sububsection{Soil properties}\label{P4S7.6.2}

Integration of (\ref{P4eq7.57}) and~(\ref{P4eq7.58}) requires the
specification of hydraulic conductivity and diffusivity as a
function of soil-water content.  
In the previous TESSEL version the parametric relations of~\citet{ClappandHornberger1978}
(see also~\citealp{Cosbyetal1984}) were adopted (still available as option). 
These are given by
\begin{equation}
\begin{split}
\gamma &=\gamma_{\rm sat}\biggl(\frac{\theta}{\theta_{\rm
sat}}\biggr)^{\!2b+3}\\
\lambda &=\frac{b \gamma_{\rm sat}(-\psi_{\rm sat})}{\theta_{\rm
sat}}\biggl(\frac{\theta}{\theta_{\rm sat}}\biggr)^{\!b+2}\\
\label{P4eq7.64}
\end{split}
\end{equation}
where $b$ is a non-dimensional exponent, $\gamma_{\rm sat}$ and
$\psi_{\rm sat}$ are the values of the hydraulic conductivity and
matric potential at saturation, respectively.  A minimum value is
assumed for $\lambda$ and $\gamma$ corresponding to permanent
wilting-point water content.

\citet{Cosbyetal1984} tabulate best estimates of $b, \gamma_{\rm
sat}, \Psi_{\rm sat}$  and $\theta_{\rm sat}$, for the 11 soil
classes of the US Department of Agriculture (USDA) soil
classification, based on measurements over large samples. Since the
model described here specifies only one soil type everywhere, and
because the determination of the above constants is not independent
of the values of $\theta_{\rm cap}$ and $\theta_{\rm pwp}$, the
following procedure is adopted.

A comprehensive review of measurements of $\theta_{\rm cap}$ and
$\theta_{\rm pwp}$ may be found in~\citet{Patterson1990}. Starting
from Patterson's estimates of $\theta_{\rm cap}$ and $\theta_{\rm
pwp}$ for the 11 USDA classes, a mean of the numbers corresponding
to the medium-texture soils (classes 4, 5, 7, and 8, corresponding
to silt loam, loam, silty clay loam and clay loam, respectively) is
taken. The resulting numbers are $\theta_{\rm cap} =0.3\dot{2}3~{\rm
m}^3~{\rm m}^{-3}$ and $\theta_{\rm pwp} =0.171~{\rm m}^3~{\rm
m}^{-3}$. Averaging the values of~\citet{Cosbyetal1984} for soil
moisture and soil-water conductivity at saturation for the same
classes gives the numerical values $\gamma_{\rm sat} =0.57 \times
10^{-6}~{\rm m~s}^{-1}$  and  $\theta_{\rm sat} =0.472~{\rm m}^3~{\rm
m}^{-3}$ . The Clapp and Hornberger expression for the matric
potential is
\begin{equation}%7.65
\psi  =\psi_{\rm sat}\biggl(\frac{\theta}{\theta_{\rm
sat}}\biggr)^{\!-b} \label{P4eq7.65}
\end{equation}
is used with  $\psi(\theta_{\rm pwp}) =-153$~m ($-15$ bar) and
$\psi(\theta_{\rm cap}) =-3.37$~m ($-0.33$ bar)
(see~\citealp{Hillel1982,JacqueminandNoilhan1990}) to find the
remaining constants $b$ and $\psi_{\rm sat}$. The results are $b
=6.04$ and $\Psi_{\rm sat} =-0.338$~m. The above process ensures a
soil that has an availability corresponding to the average value
of medium-texture soils, and yields a quantitative definite
hydraulic meaning to $\theta_{\rm cap}$ and $\theta_{\rm pwp}$
compatible with the Van Genuchten relations (see Table
\ref{P4T7.6} for a summary of the soil constants).

%---------------------------------------------------------------------------

The \cite{vg80} formulation used in HTESSEL provides a closed-form analytical expression for the conductivity, given as a function of the
pressure head, $h$, as
\begin{equation}
\gamma =\gamma_{\rm sat} \frac {[(1+(\alpha h)^{n})^{1-1/n}-(\alpha h)^{n-1}]^2}{(1+(\alpha h)^{n})^{(1-1/n)(l+2)}}
\end{equation}
where $\alpha$, $n$ and $l$ are soil-texture dependent parameters. Pressure head $h$ is linked to the soil moisture by the expression
\begin{equation}
\theta(h) = \theta_r + \frac{\theta_{sat}-\theta_{r}}{(1+(\alpha h)^{n})^{1-1/n}}
\end{equation}
The VG scheme is recognized among soil physicists as capable of reproducing both the soil water retention and the
hydraulic conductivity, and has shown good agreement with observations in intercomparison studies \citep{shao99}.
%T7.6
\begin{table}
\centering
\caption{
Van Genuchten soil parameters.
\label{P4T7.6}}
\small
\begin{tabular}{lcccc} \hline
\bf{Texture} & $\alpha$ & $l$ & $n$ & $\gamma_{\rm sat}$ \\
Units & $m^{-1}$ & - & - & $10^{-6}{\rm m/s}$ \\ \hline \hline
Coarse      & 3.83 &  1.250  & 1.38 & 6.94\\
Medium      & 3.14 & -2.342  & 1.18 & 1.16 \\
Medium-Fine & 0.83 & -0.588  & 1.25 & 0.26 \\
Fine        & 3.67 & -1.977  & 1.10 & 2.87 \\
Very Fine   & 2.65 &  2.500  & 1.10 & 1.74 \\
Ext.Trop. Organic     & 1.30 &  0.400  & 1.20 & 0.93 \\
Trop. Organic      & 3.14 & -2.342  & 1.18 & 1.16 \\ \hline
\end{tabular}
\end{table}
Table \ref{P4T7.6} lists parameter values for six soil textures for the VG scheme. HTESSEL
uses the dominant soil texture class for each gridpoint. This information is taken from
the FAO \citep{fao03} dataset.
The permanent wilting point and the soil field capacity are obtained by a specified matric potential of
$\psi(\theta_{\rm pwp})=-15~$bar and $\psi(\theta_{\rm cap}) =-0.10~$bar, respectively.
%T7.7
\begin{table}
 \centering
\caption{
Values for the volumetric soil moisture in Van Genuchten and Clapp-Hornberger (CH, loamy; bottom row),
at saturation, $\theta_{\rm sat}$, field capacity, $\theta_{\rm cap}$, permanent wilting point, $\theta_{\rm pwp}$, and
$\theta_{\rm res}$, residual moisture. 
Last column reports the plant available soil moisture. Units are $[ m^3m^{-3}]$.
\label{P4T7.7}}
\small
\begin{tabular}{lccccc} \hline
\bf{Texture} & $\theta_{\rm sat}$ & $\theta_{\rm cap}$ & $\theta_{\rm pwp}$ & $\theta_{\rm res}$ & $\theta_{\rm cap}-\theta_{\rm pwp}$ \\ \hline \hline
Coarse & 0.403 & 0.244 & 0.059 & 0.025 & 0.185 \\
Medium & 0.439 & 0.347 & 0.151 & 0.010 & 0.196 \\
Medium-Fine & 0.430 & 0.383 & 0.133 & 0.010 & 0.251 \\
Fine & 0.520 & 0.448 & 0.279 & 0.010 & 0.170 \\
Very Fine & 0.614 & 0.541 & 0.335 & 0.010 & 0.207 \\
Ext.Trop. Organic & 0.766 & 0.663 & 0.267 & 0.010 & 0.396 \\
Trop. Organic & 0.439 & 0.347 & 0.151 & 0.010 & 0.196 \\ \hline
Old-Loamy (CH) & 0.472 & 0.323 & 0.171 & - & 0.151 \\ \hline
\end{tabular}
\end{table}
In Table \ref{P4T7.7} the volumetric soil moistures associated with each soil class are shown for saturation, field capacity, wilting point, and residual moisture.
Also shown is the plant available water content and the percentage of land points in each class.
The last row shows the corresponding values for the single loamy soil used in the CH formulation in TESSEL.
Note that the plant available soil water is greater for all the new soil classes in HTESSEL.
%########################## FIGURE 2
\begin{figure}[h!]
\centering \includegraphics[width=130mm]{Figures/figA3.eps} 
\caption{Hydraulic properties of TESSEL and HTESSEL: (a) Diffusivity  and (b) conductivity. The (+) symbols
on the curves highlight (from high to low values) saturation, field capacity permanent wilting point.
}
\label{F2}
\end{figure}
Figure \ref{F2} shows the soil hydraulic diffusivity and conductivity for the TESSEL CH formulation
and the six VG soil texture classes in HTESSEL.
In TESSEL those were not allowed to fall below their wilting point values.
At saturation, TESSEL has the highest diffusivity and conductivity.
The reduced values for fine soils in HTESSEL reduces the infiltration of water and consequently the baseflow.

%---------------------------------------------------------------------------


\sububsection{Runoff}\label{P4S7.6.3}
In general when the water flux at the surface exceeds the 
maximum infiltration rate, the excess water
is put into surface runoff. 
A general formulation of surface runoff can be written as:
\begin{equation}
R =T + M - I_{\rm max}
\label{Rinf}
\end{equation}
where $I_{\rm max}$ is the maximum infiltration rate,
$T$ the throughfall precipitation and $M$ the snow melting. 
Different runoff schemes differ in the formulation of the infiltration.
%########################## FIGURE 3
\begin{figure}[h!]
\centering \includegraphics[width=65mm]{Figures/figA4.eps}
\caption{Surface runoff generation (rate $mm/h$) as a function of the b parameter 
(accounting for sub-grid effects of orography), when exposed to a precipitation rate of $10~mm/h$.
\label{F3}}
\end{figure}
In TESSEL a maximum infiltration rate at the surface, 
defined by the maximum downward diffusion from a saturated surface 
was used to define the runoff term. 
The maximum infiltration rate $I_{\rm max}$ is calculated as
\begin{equation}
I_{\rm max} = \rho_{\rm w} \biggl( \frac{b_c\gamma_{\rm sat}(-\psi_{\rm sat})}{\theta_{\rm sat}} \frac{\theta_{\rm sat}-\theta_{1}}{z_1/2}
+\gamma_{\rm sat}\biggr)
\label{Rinf2}
\end{equation}
where $\rho_{\rm w}$ is the water density, and $z_1$ is the depth of the first soil model layer (7 cm).
At typical NWP model resolutions this scheme is active only in the presence of frozen soil,
when downward soil water transfer is inhibited, otherwise it hardly ever produces runoff, as shown in \cite{boone04}.

In HTESSEL A variable infiltration rate, first introduced in the so-called Arno scheme by \cite{dumenil92},
accounts for the sub-grid variability related to orography and considers that the runoff can
(for any precipitation amount and soil condition) occur on a fraction $s$ of the grid-point area $S$.
\begin{equation}
\frac {s}{S} = 1 - \biggl( 1 - \frac {W} {W_{\rm sat}} \biggr) ^b\\ \hspace{2cm}
b = \frac {\sigma_{\rm or} - \sigma_{\rm min}} {\sigma_{\rm or} + \sigma_{\rm max}}
\label{S}
\end{equation}
where $W$ and $W_{\rm sat}$ are vertically integrated soil water contents ($\theta$ and $\theta_{\rm sat}$) over the first 50 cm of
soil defined as an effective depth for surface runoff.
Parameter $b$ is spatially variable, depends on standard deviation of orography ($\sigma_{\rm or}$), and is allowed to vary
between $0.01$ and $0.5$. The parameters $\sigma_{\rm min}$ and $\sigma_{\rm max}$ are set to 100 m and 1000 m respectively
as in \cite{vhpv03}.

The surface runoff is obtained by the Hortonian runoff formulation by integrating (\ref{S}) over the gridbox.
\begin{equation}
I_{\rm max} = (W_{\rm sat} - W ) + {\rm max} \biggl[ 0, W_{\rm sat} \biggl[ \biggl( 1 - \frac {W}{W_{\rm sat}} \biggr)^{\frac{1}{b+1}}
-\biggl( \frac {T+M}{(b+1)W_{\rm sat}}\biggr) \biggr]^{b+1} \biggr]
\label{vicr}
\end{equation}
Whenever rain or snow melt occurs, a fraction of the water is removed as surface runoff.
The ratio runoff/precipitation scales with the standard deviation of orography,
and therefore depends on the complexity represented in the gridbox, as well as on
soil texture and soil water content via $W$ and $W_{\rm sat}$.

In Figure \ref{F3} the response to a $10~{\rm mm/h}$ rain rate for the six VG soil
types and for the CH case in TESSEL is shown as a function of the $b$ parameter.
At field capacity, the surface runoff may vary from roughly $1\%$ to $50\%$ 
of the rainfall (snow melting) rate, generally increasing with finer textures 
and orographic complexity.

\subsubsection{Water transport in frozen soil}\label{P4S7.6.4}

Finally, the water transport is limited in the case
of a partially frozen soil, by considering the effective hydraulic
conductivity and diffusivity to be a weighted average of the
values for total soil water and a very small value (for
convenience, taken as the value of (\ref{P4eq7.64}) at the
permanent wilting point) for frozen water.

\subsubsection{Discretization and the root profile}\label{P4S7.6.5}

A common soil discretization is chosen for the thermal and water
soil balance for ease of interpretation of the results, proper
accounting of the energy involved in freezing/melting soil water,
and simplicity of the code. Equations (\ref{P4eq7.57})
and~(\ref{P4eq7.58})  are discretized in space in a similar way to
the temperature equations, i.e., soil water and root extraction
defined at full layers, $\theta_k$ and $\rho_{\rm w}S_{\theta,
k}$, and $F_{k+1/2}$ the flux of water at the interface between
layer $k$ and $k+1$. The resulting system of equations represents
an implicit, water-conserving method.

For improved accuracy, the hydraulic diffusivity and conductivity
are taken as (see~\citealp{MahrtandPan1984})
\begin{equation}%eq7.66
\begin{split}
\lambda_{k +1/2}  &= (1 -f_{\rm fr}^*)\lambda [\max (\theta_k^n,
\theta^n_{k+1})] +f_{\rm fr}^* \lambda (\theta_{\rm pwp})\\
\lambda_{k +1/2}  &= (1 -f_{\rm fr}^*)\gamma [\max (\theta_k^n,
\theta^n_{k+1})] +f_{\rm fr}^* \gamma (\theta_{\rm pwp})
\label{P4eq7.66}
\end{split}
\end{equation}
where $f_{\rm fr}^*= \min [f_{\rm fr}(\theta_k),f_{\rm fr}
(\theta_{k+1})]$. The boundary conditions are given by
\begin{equation}%eq7.67
\begin{split}
F_{4+1/2}  &= \rho_{\rm w}\gamma_4\\
F_{1/2}   &= T +M_{\rm sn} -y_{\rm sfc} +E_{1/2} \label{P4eq7.67}
\end{split}
\end{equation}
The difference between throughfall $T$ and surface runoff $Y_{\rm
sfc}$ (kg~m$^{-2}$ s$^{-1}$) is the soil infiltration at the
surface:
\begin{align}%eq7.68
T  &= T_{\rm 1s} +T_{c{\rm v}} \nonumber\\
y_{\it sfc}  &= \max (0, T_{\rm 1s} +M_{\rm sn} -I_{\rm f, mx})
+\frac{\max(0, F_{\rm cv} +T_{\rm cv} -I_{\rm f, mx})}{F_{\rm cv}}\\
I_{\rm f, mx}  &= \rho_{\rm w}
\biggl[\lambda_{1/2}\frac{(\theta_{\rm sat}- \theta_1)}{0.5D_1}
+\gamma_{1,2}\biggr] \nonumber \label{P4eq7.68}
\end{align}
and $\lambda_{1/2} =f_{\rm fr}^* \lambda (\theta_{\rm pwp})+(1
-f_{\rm fr}^* )\lambda (\theta_{\rm sat})$, with a similar equation
for $\gamma_{1/2}$. The evaporation at the top of the soil layer,
$E_{1/2}$, is computed as the sum of the evaporations of tile 8 plus
the contributions necessary to conserve water with the solver of the
interception layer.
\begin{enumerate}
\item[(i)] Tile 3 mismatch (after the evaporated water used by the
interception reservoir for the given tile is subtracted).

\item[(ii)] When the evaporative fluxes are downward (i.e., dew
deposition), the evaporation for tiles 4, 6 and the canopy
evaporation of tile~7.
\end{enumerate}
Root extraction is computed as sum of the transpiration contributions
(expressed in $kg m^{-2}$) over tiles $i$ and soil layers $j$ divided
by the layer depth, such as 
\begin{equation}%eq7.69
\rho_{\rm w} S_{\theta,k} =\sum_i C_i \frac{E_i}{D_k}\frac{R_k
\theta_k}{\sum_j R_j \theta_j} \label{P4eq7.69}
\end{equation}
where the sum over tiles $i$ is done for tiles 4, 6, and 7.
In case of dew deposition (i.e., tile downward evaporative flux), $S_{\theta, k} =0$.

\subsection{The ecosystem exchanges parametrization}\label{P4S7.7}

In order to obtain the net exchange of CO${}_{2 }$ between the land surface and atmosphere by a NWP model, the Gross Primary Production (GPP) 
and the ecosystem respiration ($R_{\rm eco} $) needs to be represented. 
In CHTESSEL  \citep{Boussettaetal2013} the GPP [kg CO${}_{2}$ m${}^{-2}$ s${}^{-1}$] is parametrized following the so-called A-gs approach while $R_{\rm eco} $ [kg CO${}_{2}$ m${}^{-2}$ s${}^{-1}$] 
is parametrized as a function of soil temperature, soil moisture, snow depth and vegetation type rather than prognostically computed. Schemes relying 
on prognostic land carbon "pools" (CLM-DGVM, \citealp{Levisetal2004}; JULES, \citealp{Clarketal2011}; ORCHIDEE, \citealp{Krinneretal2005}) 
are less practical for NWP purposes because they are difficult to initialize without a very long spin-up.
The Net Ecosystem Exchange (NEE) of CO${}_{2 }$ [kg CO${}_{2}$ m${}^{-2}$ s${}^{-1}$] between the surface and atmosphere is then given by
\begin{equation} 
NEE=GPP-R_{\rm eco}   \label{P4eq7.70}        
\end{equation} 
This quantity is the main observed variable in field site experiments.

\subsubsection{The photosynthesis scheme (A-gs)}\label{P4S7.7.1}

The A-gs approach is based on plant physiological considerations and describes the plant photosynthesis 
process and its dependence on CO${}_{2}$, temperature and soil moisture (\citealp{Jacobs1994}; \citealp{Jacobsetal1996}; \citealp{Calvetetal1998}; \citealp{Calvetetal2004}). 
The stomatal conductance is active for regulating both water vapor and CO${}_{2}$ fluxes. The gross CO${}_{2}$ assimilation 
by the canopy $A_{\rm g} $ [kg CO${}_{2}$ m${}^{-2}$s${}^{-1}$] is calculated using a photosynthesis module following \citet{Goudriaanetal1985}. 
The net assimilation $A_{\rm n} $ [kg CO${}_{2}$ m${}^{-2}$s${}^{-1}$] (i.e. the net flow of CO${}_{2}$ through 
the stomata) is $A_{\rm g} $ minus the dark respiration $R_{\rm d} $ [kg CO${}_{2}$ m${}^{-2}$s${}^{-1}$]. 
Once the net assimilation is known, the stomatal conductance for CO${}_{2}$, $g_{\rm sc} $ [m s${}^{-1}$] can be derived by Kirchhoff's resistance/conductance analogy. 
It is defined as the net flow of CO${}_{2}$ through the stomata divided by the difference between the CO${}_{2}$ concentration 
outside the leaves $C_{\rm s} $ [kg CO${}_{2}$ m${}^{-3}$] and the concentration in the intercellular cavities $C_{\rm i} $ [kg CO${}_{2}$ m${}^{-3}$] (see figure \ref{P4F7.5}):
\begin{equation}  
g_{\rm sc}^{} =\frac{A_{\rm n} (Env)}{C_{\rm s} -C_{\rm i} (Env)}   \label{P4eq7.71}       
\end{equation} 
where the functions ($Env$) indicate the dependence of $A_{\rm n} $ and $C_{\rm i} $ on the different environment factors. 
In a future version, it is planned that the stomatal conductance from the carbon module would be used for water vapor 
transpiration, which leads to a full coupling between carbon and water vapor fluxes. 

\begin{figure}
\centering
\includegraphics[width=120mm]{Figures/figA5.jpg}
\caption{Schematic representation of a leaf with the resistance analogies for carbon and water vapor, where $g_{\rm cu} $ is the cuticular conductance, 
$g_{\rm sc} $ the stomatal conductance to CO${}_{2}$, $g_{\rm s} $ the stomatal conductance to water vapor, $g_{\rm m} $ the mesophyll conductance, 
$C_{\rm s} $ and $C_{\rm i} $ the CO${}_{2}$ concentration at the canopy surface and inside the leaf cavity  respectively, 
$q_{\rm s} $ and $q_{\rm sat} (T_{\rm s}) $ the humidity at the canopy surface and the saturated humidity at 
the temperature  $T_{\rm s} $ of the canopy surface respectively.
\label{P4F7.5}}
\end{figure}

The basics of the A-gs model are described by \citet{Jacobs1994} and \citet{Jacobsetal1996}. 
Further details of the current extended formulation are given in publications by \citet{Calvetetal1998} and \citet{Calvetetal2004}. 
For the description of the canopy conductance to CO${}_{2 }$, $g_{\rm sc} $, a stepwise approach is followed with: (i) the definition of the temperature dependent parameters, 
(ii), the radiation response, (iii) the calculation of the ratio between internal and external CO${}_{2}$ concentration, 
(iv) the computation of stomatal conductance, (v) inclusion of the soil moisture response and (vi) the vertical integration over the canopy.

{\it  (i) Temperature responses}

There are several parameters in the photosynthesis model that are temperature dependent, namely the compensation point, the mesophyll conductance 
and the maximum photosynthetic capacity. The compensation point $\Gamma $ is defined as the CO${}_{2 }$ concentration at which the net CO${}_{2}$ assimilation 
of a fully lit leaf becomes zero. It can be measured in a laboratory by exposing plants to a variable CO${}_{2 }$ concentration. 
The mesophyll conductance $g_{\rm m} $ describes the transport of CO${}_{2}$ from the substomatal cavities to the mesophyll cells where the carbon is fixed. 
It includes the representation of physical and chemical processes. The maximum photosynthetic capacity $A_{\rm m,\max }$ is specified as an absolute upper 
limit to the photosynthesis rate in full sunlight and non-limiting CO${}_{2}$ concentration. 
The temperature dependence is described with so-called Q${}_{10 }$ functions, where Q${}_{10}$ represents the proportional increase 
of a parameter for a $10^{\rm o}$ increase in temperature \citep{Berryetal1982}. 

For the compensation point the formulation is:
\begin{equation} 
\Gamma (T_{\rm s} )=\Gamma (25^{\rm o} )Q_{10\Gamma }^{(T_{\rm s} -25)/10} \label{P4eq7.72} 
\end{equation} 
 where $\Gamma (25^{\rm o} )$ is the compensation point at 25${}^{\rm o}$ C, $Q_{10\Gamma } $ is the $Q_{10} $-constant and $T_{\rm s} $is the leaf surface temperature. For $g_{\rm m} $ and $A_{\rm m,\max } $, the temperature dependence is further adjusted by the inhibition functions after \citet{Collatzetal1992}:
\begin{equation} 
g_{\rm m} (T_{\rm s} )=\frac{g_{\rm m} (25^{\rm o} )Q_{\rm 10gm}^{(T_{\rm s} -25)/10} }{(1+e^{0.3(T_{\rm 1gm} -T_{\rm s} )} )(1+e^{0.3(T_{\rm s} -T_{\rm 2gm} )} )} \label{P4eq7.73}
\end{equation}
\begin{equation}
A_{\rm m,\max } (T_{\rm s} )=\frac{A_{\rm m,\max } (25^{\rm o} )Q_{\rm 10Am,\max }^{(T_{\rm s} -25)/10} }{(1+e^{0.3(T_{\rm 1am,\max } -T_{\rm s} )} )(1+e^{0.3(T_{\rm s} -T_{2am,\max } )} )} \label{P4eq7.74}
\end{equation}
where $Q_{\rm 10gm} ,\, Q_{\rm 10Am,\max } ,\, T_{\rm 1gm} ,\, T_{\rm 2gm} ,\, T_{\rm 1Am,\max } ,\, T_{\rm 2Am,\max } $, 
are constants affecting the sensitivity to the plant surface temperature $T_{\rm s} $. Parameter $g_{\rm m} (25^{\rm o} )$
depends on soil moisture stress and will be further described in the following sections. 
Its unstressed value $g_{\rm m}^{*} (25^{\rm o} )$ is optimized here with the help of observations. The constants in these functions are vegetation type dependent and 
are listed in Tables \ref{P4T7.CO2_1} and \ref{P4T7.CO2_2}.

\begin{table}
\centering
\caption{Parameter values as specified in the optimized CTESSEL model. The vegetation types are from the 
GLCC database and used in the same way as in TESSEL \citet{vandenHurketal2000}
\label{P4T7.CO2_1}}
\small
%\begin{tabular}{clccccccc}
%\begin{tabular}{|p{0.25in}|p{1.2in}|p{0.5in}|p{0.55in}|p{0.55in}|p{0.55in}|p{0.7in}|p{0.4in}|p{0.3in}|} \hline 
\begin{tabular}{p{0.25in} p{1.2in} p{0.5in} p{0.55in} p{0.55in} p{0.55in} p{0.7in} p{0.4in} p{0.3in}}
\toprule
 & & \bf{\textit{R${}_{0}$}} & \bf{\textit{g${}_{\rm m}$*${\left(25 \right)}$}} & \bf{\textit{g${}_{\rm c}$ }} & \bf{\textit{Dmax${}^{*}$ }}& 
\bf{\textit{Am,max${\left(25 \right)}$}} &\bf{\textit{f${}_{0}$${}^{*}$ }}& \bf{\textit{${\Gamma}{\left(25 \right)}$ }}\\  
%
\bf{\textit{Type}} & \bf{\textit{Vegetation type}} & [mgCO${}_{2}$ m${}^{-2}$s${}^{-1}$] & 
[mm s${}^{-1}$] & [mm s${}^{-1}$] & [kg kg${}^{-1}$] & [mgCO${}_{2}$ m${}^{-2}$ s${}^{-1}$] & [-] & 
[ppm] \\  \hline \hline 
%\midrule
%\textbf{\textit{Type code}} & \textbf{\textit{Vegetation type}} & \textbf{\textit{R${}_{0}$\newline }}[mgCO${}_{2}$ m${}^{-2}$s${}^{-1}$]\textbf{\textit{}} & 
%\textbf{\textit{g${}_{m}$*${\left(25 \right)}$\newline }}[mm s${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{g${}_{c}$\newline }}[mm s${}^{-1}$]\textbf{\textit{}} & 
%\textbf{\textit{Dmax${}^{*}$\newline }}[kg kg${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{Am,max${\left(25 \right)}$
%\newline}}[mgCO${}_{2}$         m${}^{-2}$ s${}^{-1}$]\textbf{\textit{}} & \textbf{\textit{f${}_{0}$${}^{*}$\newline }}[-] & 
%\textbf{\textit{${\Gamma}{\left(25 \right)}$\newline }}[ppm]\textbf{\textit{}} \\  \hline \hline 
1 & Crops, mixed farming & 0.100 & 1.3 & 0.15 & Eq.\ref{P4eq7.89} & 2.20 & 0.85 & 42 \\ 
2 & Short grass & 0.080 & 1.3 & 0.20 & Eq.\ref{P4eq7.89} & 3.00 & 0.65 & 42 \\ 
3 & Evergreen needleleaf  & 0.360 & 0.8 & 0.20 & 0.124 & 2.20 & Eq.\ref{P4eq7.93} & 42 \\
4 & Deciduous needleleaf  & 0.330 & 0.8 & 0.20 & 0.124 & 2.20 & Eq.\ref{P4eq7.93} & 42 \\ 
5 & Deciduous broadleaf  & 0.280 & 1.4 & 0.00 & 0.109 & 1.83 & Eq.\ref{P4eq7.93} & 42 \\ 
6 & Evergreen broadleaf  & 0.270 & 1.1 & 0.25 & 0.124 & 1.83 & Eq.\ref{P4eq7.93} & 42 \\
7 & Tall grass & 0.150 & 2.3 & 0.20 & Eq.\ref{P4eq7.89} & 1.83 & 0.70 & 2.6 \\ 
8 & Desert & - & - & - & - & - & - & - \\
9 & Tundra & 0.360 & 2.0 & 0.25 & Eq.\ref{P4eq7.89} & 3.00 & 0.95 & 42 \\ 
10 & Irrigated crops & 0.096 & 1.4 & 0.25 & Eq.\ref{P4eq7.89} & 1.83 & 0.92 & 42 \\ 
11 & Semidesert & 0.019 & 1.0 & 0.25 & Eq.\ref{P4eq7.89} & 1.83 & 0.80 & 42 \\ 
12 & Ice caps and glaciers & - & - & - & - & - & - & - \\
13 & Bogs and marshes & 0.270 & 0.5 & 0.25 & Eq.\ref{P4eq7.89} & 1.83 & 0.96 & 42 \\ 
14 & Inland water & - & - & - & - & - & - & - \\ 
15 & Ocean & - & - & - & - & - & - & - \\ 
16 & Evergreen shrubs & 0.110 & 0.9 & 0.15 & Eq.\ref{P4eq7.89} & 1.83 & 0.72 & 2.6 \\ 
17 & Deciduous shrubs & 0.080 & 1.9 & 0.20 & Eq.\ref{P4eq7.89} & 1.83 & 0.96 & 42 \\ 
18 & Mixed forest- Wood & 0.420 & 1.0 & 0.00 & 0.124 & 2.20 & Eq.\ref{P4eq7.93} & 42 \\ 
19 & Interrupted forest & 0.160 & 0.8 & 0.10 & 0.124 & 2.20 & Eq.\ref{P4eq7.93} & 42 \\ 
20 & Water -land mixtures & 0.270 & 1.0 & 0.25 & Eq.\ref{P4eq7.89} & 1.83 & 0.95 & 42 \\ %\hline 
\bottomrule
\end{tabular}
\end{table}

%\noindent 

%\noindent \eject 

\begin{table}
\centering
\caption{
CTESSEL temperature response, quantum use efficiency and soil moisture stress parameters
\label{P4T7.CO2_2}}
%\begin{tabular}{|p{0.8in}|p{0.5in}|p{1.7in}|} \hline
\small
\begin{tabular} {lcccc}
\toprule
\bf{Parameter} & \bf{Value} & \bf{Notes} \\
 \hline \hline
$Q_{10}$${}_{\Gamma }$ & 2 &  \\
$Q_{\rm 10gm}$ & 2 &  \\ 
$Q_{\rm 10Am,max}$ & 2 &  \\ 
$T_{\rm 1gm}$ (${}^{\rm o}$C) & 5 & (13 for veg. types 7 and 16) \\ 
$T_{\rm 1Am,max }$(${}^{\rm o}$C) & 8 & (13 for veg. types 7 and 16) \\ 
$T_{\rm 2gm }$(${}^{\rm o}$C) & 36 &  \\ 
$T_{\rm 2Am,max }$(${}^{\rm o}$C) & 38 &  \\ 
\textit{a} & 2.381 & (5.323 for veg.types 7 and 16 \\
\textit{b} & 0.6103 & (0.8923 for veg.types 7 and 16) \\
$D_{\max }^{\rm X} $(g/kg) & 300 &  \\ 
$\varepsilon $${}_{\rm o}$ (mg CO${}_{2}$/J PAR) & 0.0142 & (0.0117 for vegetation types 7 and 16) \\ \hline 
\end{tabular}
\end{table}

{\it  (ii) Radiation and CO${}_{2}$ response}

For $A_{n} $, two regimes are distinguished: the radiation limiting regime and the CO${}_{2}$ limiting regime (\citealp{Goudriaanetal1985}; \citealp{Jacobs1994}). 
In the radiation limiting regime with sufficient CO${}_{2}$, $A_{\rm n}$ is controlled by the amount of photosynthetically active radiation (PAR) $I_{\rm a} $
\begin{equation}
A_{\rm n} =\varepsilon {\kern 1pt} I_{\rm a} -R_{\rm d} \label{P4eq7.75}
\end{equation}
where $R_{\rm d}$ is the dark respiration and where $\varepsilon$ is the quantum efficiency expressed according to 
\begin{equation}
\varepsilon {\kern 1pt} =\varepsilon _{\rm o} \frac{C_{\rm i} -\Gamma }{C_{\rm i} +2\Gamma } \label{P4eq7.76}
\end{equation}
Parameter $\varepsilon _{\rm o} $ is the maximum quantum use efficiency and $C_{\rm s}$ is the ambient CO${}_{2}$ concentration at the leaf surface. 

At high radiation intensities, net assimilation saturates at a level $A_{\rm m}$ and becomes CO${}_{2}$ limited according 
to $A_{\rm m} =\left(C_{\rm i} -\Gamma \right)\, g_{\rm m} $ with $C_{\rm i}$ the CO${}_{2}$ concentration inside the leaf cavities (see below). 
An absolute limit to account for the maximum photosynthetic capacity of the leaves is further applied as follows 
\begin{equation}
A_{\rm m} =A_{\rm m,\max } \left[1-\, \exp \{ -g_{\rm m} \left(C_{\rm i} -\Gamma \right)\, /A_{\rm m,\max } \} \right]\label{P4eq7.77}
\end{equation}

The radiation and CO${}_{2}$ limiting regimes are combined via a smooth exponential transition function
\begin{equation}
A_{\rm n} =(A_{\rm m} +R_{\rm d} )\left[1-\exp \left(\frac{-\varepsilon I_{\rm a} }{A_{\rm m} +R_{\rm d} } \right)\right]-R_{\rm d} \label{P4eq7.78}
\end{equation}
The autotrophic dark respiration is simply parametrized according to \citet{VanHeemst1986} and includes the respiration from the leaves only
\begin{equation}
R_{\rm d} =A_{\rm m} /9 \label{P4eq7.79}
\end{equation} The respiration from other parts of the vegetation is included in the parametrized heterotrophic respiration.

{\it  (iii) The $C_{\rm i} /C_{\rm s} $ ratio}

The CO${}_{2}$ concentration inside the leaf cavities \textit{$C_{\rm i} $} needs to be known in order to derive the stomatal conductance from the net 
assimilation. Observations indicate that the ratio $C_{\rm i} /C_{\rm s} $ is a rather conservative quantity for moist atmospheric conditions and that increasing 
humidity deficit exerts a strong stomatal control affecting this ratio. Therefore $C_{\rm i} /C_{\rm s}$ is specified as a function of atmospheric moisture 
deficit \textit{$D_{\rm s}$} at the leaf surface.
\begin{equation}
\frac{C_{\rm i} }{C_{\rm s} } =f+(1-f)\frac{\Gamma }{C_{\rm s} } \label{P4eq7.80}  
\end{equation}
where $f$ is the coupling factor defined by:
\begin{equation}
f=f_{\rm o} \left(1-\frac{D_{\rm s} }{D_{\max } } \right)+f_{\min } \frac{D_{\rm s} }{D_{\max } } \label{P4eq7.81}
\end{equation}
and $f_{\rm o}$ is the value of $f$ at $D_{\rm s}$, $D_{\max }$ is the maximum saturation deficit and 
\begin{equation}
f_{\min } =\frac{g_{\rm cu} }{g_{\rm cu} +g_{\rm m} } \label{P4eq7.82}
\end{equation}
The transport of CO${}_{2}$ is maintained in the situation where $f_{} =f_{\min } $ through the leaf cuticle or because of imperfect closure of the stomata. 
This process is represented by the cuticule with conductance  $g_{\rm cu}^{} $.


{\it  (iv) Stomatal conductance }

The first computation of the stomatal conductance for CO${}_{2}$, $g_{\rm sc}^{1}$ is achieved by dividing net assimilation by the difference 
between CO${}_{2}$ concentration in and outside the leaves. It is modified here to account for the limiting cases of very dry air and dark respiration:
\begin{equation}
g_{\rm sc}^{1} =\frac{A_{\rm n} -A_{\min } \left(\frac{D_{\rm s} }{D_{\max } } \frac{A_{\rm n} +R_{\rm d} }{A_{\rm m} +R_{\rm d} } 
\right)+R_{\rm d} \left(1-\frac{A_{\rm n} +R_{\rm d} }{A_{\rm m} +R_{\rm d} } \right)}{C_{\rm s} -C_{\rm i} } \label{P4eq7.83}
\end{equation}
where $A_{\min } $ represents the residual photosynthesis rate (at full light intensity) associated with cuticular transfers when the stomata are closed because 
of a high specific humidity deficit:
\begin{equation}
A_{\min } =g_{\rm m} (C_{\min } -\Gamma )\label{P4eq7.84}
\end{equation}
In this equation, $C_{\min }$ is the value of $C_{\rm i}$ at maximum specific humidity deficit:
\begin{equation}
C_{\min } =\frac{g_{\rm cu} C_{\rm s} +g_{\rm m} \Gamma }{g_{\rm cu} +g_{\rm m} } \label{P4eq7.85}
\end{equation}
The diffusion of CO${}_{2}$ through the stomatal openings interacts with that of water vapour and therefore stomatal conductance to CO${}_{2}$ is corrected for 
this interaction by an iterative refinement:
\begin{equation}
g_{\rm sc} =g_{\rm sc}^{1} +E\frac{M_{\rm a} }{\rho _{\rm a} M_{\rm v} } \frac{C_{\rm s} +C_{\rm i} }{2(C_{\rm s} -C_{\rm i} )} \label{P4eq7.86}
\end{equation}
where $M_{\rm v}$ and $M_{\rm a}$ are molecular masses of water vapour and air respectively, $\rho _{\rm a} $ is the air density and $E$ is the leaf transpiration based on the 
previous guess of the stomatal conductance: 
\begin{equation}
E=(1.6g_{\rm sc}^{1} )D_{\rm s} \rho _{\rm a} \label{P4eq7.87}
\end{equation}
Finally, the stomatal conductance to water vapour $g_{\rm s} $ is given by:
\begin{equation}
g_{\rm s} =1.6g_{\rm sc} \label{P4eq7.88}
\end{equation}
The total conductance used by the transpiration scheme is $g_{\rm s} +g_{\rm cu} $, where $g_{\rm cu} $ is the vegetation type dependent cuticular conductance (Table \ref{P4T7.CO2_1}).


{\it  (v) Soil moisture stress response}

Among other possible A-gs formulations for which the soil moisture stress response is directly applied to the gross assimilation $A_{\rm g}$ \citep{Rondaetal2001} 
or the net assimilation $A_{\rm n} $ \citep{Salaetal1996}, \citet{Calvet2000} found that the soil moisture stress response is driven in a complex 
way through the mesophyll conductance $g_{\rm m} $, the maximum specific humidity deficit tolerated by the vegetation $D_{\max } $, and the ratio $C_{\rm i} /C_{\rm s}$ controlled 
by $f$. The soil moisture response behaves differently for high and low vegetation. In CHTESSEL the adopted soil moisture stress response follow the function 
described in \citet{Calvet2000} and \citet{Calvetetal2004} and based on a meta-analysis of several herbaceous and woody vegetation types. 


\noindent \textit{Low vegetation formulation:}

\citet{Calvet2000} found that the mesophyll conductance $g_{\rm m} $ and the maximum atmospheric moisture deficit $D_{\max } $ vary with soil moisture but that they remain 
correlated according to:
\begin{equation}
\ln \left(g_{\rm m} (25{}^\circ )\right)=a-b\, \, \ln (D_{\max }^{} )\label{P4eq7.89}
\end{equation}
where $g_{\rm m} (25{}^\circ )$ in [mm s${}^{-1}$] and $D_{\max }$  in [g kg${}^{-1}$] 

Therefore this equation is used to derive $D_{\max }^{*} $ (maximum saturation deficit without soil moisture stress) from the tabulated $g_{\rm m}^{*} (25{}^\circ )$, 
where superscript * indicates optimal soil moisture conditions and $a$ and $b$ are tabulated empirical coefficients (Table \ref{P4T7.CO2_2}). Then the soil moisture stress 
index $f_{2} $ (see \ref{P4eq7.10}) is applied to find $D_{\max }^{} $in stressed conditions according to a bilinear function with a breakpoint at a critical 
soil moisture stress index $f_{2c} $:
\begin{equation}
\begin{array}{ll}
D_{\max } =D_{\max }^{\rm X} \frac{f_{2} }{f_{\rm 2c} } & for\; f_{2} <f_{\rm 2c}\\    
D_{\max } =D_{\max }^{\rm X} +(D_{\max }^{*} -D_{\max }^{\rm X} )\frac{f_{2} -f_{\rm 2c} }{1-f_{\rm 2c} } & for\; f_{2} \ge f_{\rm 2c}
\end{array}  \label{P4eq7.90}   
\end{equation}
where $D_{\max }^{\rm X} $ is the maximum value of $D_{\max } $ corresponding to $f_{\rm 2c} $. The resulting $D_{\max }^{} $ is substituted in (\ref{P4eq7.89}) to 
find $g_{\rm m} (25{}^\circ )$. 


\noindent \textit{High vegetation formulation:}

Observations show that $g_{\rm m} (25{}^\circ )$ is well correlated with the coupling factor $f_{0} $ according to the following empirical expression \citep{Calvetetal2004}:
\begin{equation}
\ln \left(g_{\rm m}^{*} (25{}^\circ )\right)=4.7-7.0f_{0}^{*} \label{P4eq7.91} 
\end{equation}
where $g_{\rm m} (25{}^\circ )$ is in [mm s${}^{-1}$]. 
In this case (\ref{P4eq7.91}) is used to derive $g_{\rm m}^{*} (25{}^\circ )$ from the value $f_{0}^{*}$ as tabulated according to vegetation type ( Table \ref{P4T7.CO2_1}). 
Subsequently a soil moisture stress function is applied to find $g_{\rm m}^{} (25{}^\circ )$:
\begin{equation}
\begin{array}{ll}
g_{\rm m} (25{}^\circ )=g_{\rm m}^{\rm N} \frac{f_{2} }{f_{\rm 2c} } & for\; f_{2} <f_{\rm 2c} \\
g_{\rm m} (25{}^\circ )=g_{\rm m}^{\rm N} +(g_{\rm m}^{*} (25{}^\circ )-g_{\rm m}^{\rm N} )\frac{f_{2} -f_{\rm 2c} }{1-f_{\rm 2c} } & for\; f_{2} \ge f_{\rm 2c} 
\end{array}  \label{P4eq7.92}
\end{equation}
where $g_{\rm m}^{\rm N} $ is the stressed value of $g_{\rm m} $ derived from \citet{Calvetetal2004} meta-analysis with the following empirical function:
\begin{equation}
\ln \left(g_{\rm m}^{\rm N} \right)=2.8-7.0f_{0}^{*} \label{P4eq7.93}
\end{equation}
$g_{\rm m}^{\rm N}$ is in [mm s${}^{-1}$].

After computing $g_{\rm m} (25{}^\circ $ according to (\ref{P4eq7.92}), the stressed value for $f_{0}^{}$ is derived with
\begin{equation}
\ln \left(g_{\rm m}^{} (25{}^\circ )\right)=2.8-7.0f_{0}^{} \label{P4eq7.94}
\end{equation}
$g_{\rm m} (25{}^\circ )$ is in [mm s${}^{-1}$].

Further details on the soil stress parametrization can be found in \citet{Calvet2000}, \citet{Calvetetal2004} and  \citet{Voogtetal2006}.


{\it  (vi) Vertical integration from leaf to canopy}

The net CO${}_{2}$ assimilation calculated at the leaf scale is upscaled to the canopy scale assuming that leaf parameters do not vary with height in the canopy, 
and that the attenuation of the incoming shortwave radiation in the canopy can be computed using a simple radiative extinction model. The incoming PAR above 
the vegetation ($I_{\rm a} (h)$, with $h$ the canopy height) is assumed to be 48\% of the incoming shortwave radiation and then further attenuated in the canopy. 
The dependence of PAR on height $z$ within the canopy is described by \citet{Roujean1996} according to:
\begin{equation}
I_{\rm a} (z)=I_{\rm a} (h)(1-K(z)) \label{P4eq7.95}
\end{equation}
where $K$ is the extinction function given by:
\begin{equation}
K(z)=\delta (\mu _{\rm s} )K_{\rm df} (z)+(1-\delta (\mu _{\rm s} )K_{\rm dr} (z)) \label{P4eq7.96}
\end{equation}
$K_{\rm df} (z)$ and $K_{\rm dr} (z)$ are the extinction coefficients of diffuse and direct light, respectively:
\begin{equation}
K_{\rm df} (z)=1-e^{\left(\frac{-0.8bLAI(h-z)}{h} \right)}, \label{P4eq7.97}
\end{equation}
\begin{equation}
K_{\rm dr} (z)=1-e^{\left(-\frac{G_l}{\cos (\mu _{\rm s} )} \frac{bLAI(h-z)}{h} \right)}  \label{P4eq7.98}
\end{equation}
where $\mu _{\rm s} $ is the solar zenith angle and $G_l$ is a parameter that describes the distribution of leaves (a spherical angular distribution is assumed with $G_l=0.5$), 
$\delta $ is the ratio of diffuse to total downward shortwave radiation at the top of the canopy, $LAI(h-z)$ is the cumulative leaf area index above height $z$ and 
$b$ is the foliage scattering coefficient given by:
\begin{equation}
b=1-\frac{1-\sqrt{1-\omega } }{1+\sqrt{1-\omega } } \label{P4eq7.99} 
\end{equation}
based on the leaf single scattering albedo $\omega $ (=0.2) for the solar spectrum corresponding to the PAR. Parameter $\delta $ is given by:
\begin{equation}
\delta (\mu _{\rm s} )=\frac{0.25}{0.25+\cos (\mu _{\rm s} )} \label{P4eq7.100} 
\end{equation}

Assuming a homogeneous leaf vertical distribution, the integrated canopy net CO${}_{2}$ assimilation, dark respiration and conductance can be written as:
\begin{equation}
A_{\rm nI} =LAI\int _{0}^{1}A_{\rm n} d(z/h)  \label{P4eq7.101} 
\end{equation}
\begin{equation}
R_{\rm dI} =LAI\int _{0}^{1}R_{\rm d} d(z/h)  \label{P4eq7.102} 
\end{equation}
\begin{equation}
g_{\rm c} =LAI\int _{0}^{1}g_{\rm s} d(z/h)  \label{P4eq7.103} 
\end{equation}

In the above equations, $LAI$ is defined as the ratio of leaf area covering a unit of ground area (m${}^{2}$ m${}^{-2}$).

The integrations are parametrized with a three-point Gaussian quadrature method following  \citet{Goudriaan1986}:
\begin{equation}
A_{\rm nI} =LAI\sum _{i=1}^{3}W_{i} A_{\rm n} (z_{i} )  \label{P4eq7.104} 
\end{equation}
\begin{equation}
R_{\rm dI} =LAI\sum _{i=1}^{3}W_{i} R_{\rm d} (z_{i} )  \label{P4eq7.105} 
\end{equation}
\begin{equation}
g_{\rm c} =LAI\sum _{i=1}^{3}W_{i} g_{\rm s} (z_{i} )  \label{P4eq7.106} 
\end{equation}
 where $W_{i} $ and $z_{i}$ are the Gauss weights and levels, respectively.

\noindent  
\ 
\subsubsection{Ecosystem respiration and Gross Primary Production}\label{P4S7.7.2}

The ecosystem respiration $R_{\rm eco} $ is given by two terms: the autotrophic dark respiration $R_{\rm dI} $ (\ref{P4eq7.79}), 
and $R_{\rm soilstr} $ which represents both heterotrophic respiration from the soil and autotrophic respiration from the above and below ground structural biomass. 
It is parametrized following a modified formulation of \citet{Normanetal1992} as a function of soil temperature, soil moisture, snow depth and vegetation type:
\begin{equation}
R_{\rm soilstr} =R_{0}{\left(25 \right)}Q_{\rm 10Ro}^{\left(\frac{T_{\rm soil} -25}{10} \right)} f_{\rm sm} f_{\rm sn}  \label{P4eq7.107} 
\end{equation}
In this equation $f_{\rm sn} $ and $f_{\rm sm} $ are snow and soil moisture attenuation functions respectively defined as:
\begin{equation}
f_{\rm sn} =1-C_{\rm vs} (1-e^{-\alpha z_{\rm snow} } ) \label{P4eq7.108} 
\end{equation}
$C_{\rm vs} $ is the surface fraction covered by snow, $\alpha$ is a constant expressing the attenuation of the soil carbon emission within the snow pack and 
$z_{\rm snow} $ is the snow depth. The soil moisture stress function for soil respiration is defined following a study by \citet{Albergeletal2010} as: 
\begin{equation}
f_{\rm sm} =\frac{\bar{\theta }}{\theta _{\rm cap} }  \label{P4eq7.109}
\end{equation}
$Q_{\rm 10Ro} $ represents the proportional increase of a parameter for a $10^{\rm o}$ increase in temperature  \citep{Berryetal1982}, 
this case, given its variability with climate regimes, $Q_{\rm 10Ro}$ is defined as a function of soil temperature after  \citet{McGuireetal1992}. 
The vegetation types are affecting the ecosystem respiration through a reference respiration at 25${}^\circ$C, $R_{\rm 0}{\left(25 \right)}$, estimated 
by minimizing the root mean square errors between simulated and observed $R_{\rm eco} $ for each vegetation type. 

Finally, the gross primary production $GPP$ and the ecosystem respiration are given by:
\begin{equation}
GPP=A_{\rm nI} +R_{\rm dI}  \label{P4eq7.110}
\end{equation}
\begin{equation}
R_{\rm eco} =R_{\rm dI} +R_{\rm soilstr}.  \label{P4eq7.111}
\end{equation}
These two quantities are often derived from flux-tower NEE observations and used for evaluation of the two main CO${}_{2}$ processes.

\subsubsection{Vegetation description}\label{P4S7.7.3}

The state of vegetation is given by the Leaf Area Index; crucial for deriving the plant assimilation and transpiration activity. 
Currently a satellite observation-based climatology was considered for the representation of $LAI$. The satellite product (MOD15A2) is 
derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on board of TERRA. It is produced daily for the land surfaces 
at 1 km spatial resolution from MODIS spectral reflectance with a global coverage, and synthesized on an 8-day time interval based on simultaneously 
retrieved maximum Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) in order to remove the atmospheric noise \citep{Mynenietal2002}.
The collection 5 of the product (released in 2008 available from February 2000 to present) is used. To derive the climatological time series, 9 years of data 
(2000--2008) were re-projected from the sinusoidal to a geographic regular latitude/longitude projection, spatially averaged to 1/12th degree resolution, 
then temporally smoothed, monthly averaged \citep{Jarlanetal2008} and finally interpolated to the IFS reduced Gaussian grid. 
The MODIS LAI products were analyzed and validated in previous studies (\citealp{Garriguesetal2008}; MODIS Land team http://landval.gsfc.nasa.gov/, 
\citealp{Weissetal2012}). After a positive assessment within the IFS system \citep{Boussettaetal2013}, this product was adopted by ECMWF for operational use. 
This climatology is used within CHTESSEL to drive the surface carbon exchange module. 

The land use classification follows from the Global Land Cover Characteristics (GLCC) data \citep{Lovelandetal2000} according to 
the Biosphere-Atmosphere Transfer Scheme (BATS) classification to assign dominant high and low vegetation types and associated parameters 
within each grid box (as detailed in Tables \ref{P4T7.CO2_1}, and \ref{P4T7.CO2_2}).


\subsection{Lakes and other inland water bodies}\label{P4S7.8}
The representation of inland water bodies (lakes, reservoirs, rivers and coastal waters) is important in order to account for 
the thermal inertia effects and specific albedo and roughness characteristics associated to open water and to account for phase change during freezing/melting.
This is simulated by the Fresh-water Lake model FLake \citep{Mironovetal2010}, which was chosen to be included in the IFS for its intermediate complexity,
particularly adapted for numerical weather prediction and climate applications. 
Moreover FLake benefits from a large research community efforts, contributing to validation and development (see {\sf http://www.flake.igb-berlin.de/papers.shtml}).

The FLake implementation in the HTESSEL surface scheme has been introduced in subsequent steps, with ad-hoc simplifications that are necessary 
at global scale due to the lack of detailed knowledge of the lake conditions, i.e. turbidity or presence of lake bottom sediments.
A freezing parametrisation and specific ice prognostic variables with their energy and mass balance allows FLake to be used in all climate conditions. 
At present, there is no account for a water balance equation and lake depth and surface area (or fractional cover) are kept constant in time.

The offline implementation and surface testing at global scale is presented in \citet{Dutraetal2010b}. 
\citet{Balsamoetal2010} performed a sensitivity study on the lake depth parameter, which is most directly related to lake thermal inertia.
The forecast sensitivity of FLake in the IFS is evaluated in \citet{Balsamoetal2012},
while the preparation of a retrospective lake reanalysis for initial conditions 
together with the coupled data assimilation impact assessment is provide in \citet{Balsamo2013}.

Previously to the introduction of FLAKE only grid-scale lakes (lake fraction larger than $50\%$) were considered in the IFS and treated similarly to ocean-points.
In these grid points, the surface temperature was prescribed, and evolved daily during integration according to climatological temperatures (derived from ERA-40),
without diurnal changes of surface temperature and inter-annual variability.
With the activation of FLAKE both partial lake cover within a grid-box and fully resolved lakes are modelled prognostically and evolve with the same 
timestep and spatial grid resolution used to the atmospheric and land surface processes in the IFS.

A new tile (tile 9) was created in HTESSEL to represent the lakes interface with the atmosphere.
The FLake model has been extensively modified removing components that were already part of
HTESSEL scheme such as the turbulent fluxes calculation in the tiling scheme,
and the omission of parametrisation for the treatment of lake bottom sediments for which global dataset are not available.
The treatment of snow cover is omitted although the HTESSEL snow scheme can be coupled to the lake tile over resolved lakes
allowing a consistent representation of snow cover over frozen lakes \citep{Dutraetal2010b}.

%########################## FIGURE 3
\begin{figure}[h!]
\centering \includegraphics[width=120mm]{Figures/figA6.pdf}
\caption{Schematic representation of the lake model.
\label{F6}}
\end{figure}

\subsubsection{Lake prognostic and ancillary variables}
FLake meets the requirements for physical parametrisations in weather forecast systems, in terms of both physical realism and numerical efficiency.
\cite{Mironovetal2010} presented a brief description of FLake, in its implementation into a limited-area atmospheric model.
FLake adopts a similar concept to force-restore \citep{Deardorff1978} in land surface models,
since it has intermediate features between purely bulk models and multi-layer one-dimensional models.

The vertical temperature profile specified in FLake consists of a top mixed layer, with uniform distribution of temperature, 
and a thermocline with its upper boundary located at the mixed layer bottom, and the lower boundary at the lake bottom. 
The profile in the thermocline is parameterized according to the self-similarity concept, 
which has been supported by many theoretical and observational studies (see \citealp{Mironov2008} for a discussion). 

In the IFS implementation of FLake the surface fluxes of heat, moisture and momentum are computed by the HTESSEL routines.
A new lake tile (tile 9) has been introduced to enable Flake to run on any single surface grid-point in the surface scheme,
independently from the lake fractional cover.
This is desirable to ease the interoperability of lake output at diverse spatial resolutions where the ratio of resolved/unresolved lakes and coastal water varies.

The prognostic variables included in FLake are: mixed-layer temperature $T_{\rm ML}$, mixed-layer depth $d_{\rm ML}$,
bottom temperature $T_{\rm BL}$, mean temperature of the water column $T_{\rm MN}$, shape factor ${\xi}_{\rm ML}$
(with respect to the temperature profile in the thermocline), temperature at the ice upper surface $T_{\rm Ice}$, and ice thickness $d_{\rm Ice}$,
which correspond to the schematic representation of figure \ref{F6}.


\subsubsection{Lake temperature and self-similarity}
The self-similarity assumption leads to a significant simplification in comparison with multi-layer one-dimensional models, 
while still retaining a more realistic physics than bulk models.
In particular, the unrealistic assumption of fully mixed temperature profile in a lake, applied in many bulk models, can be avoided.
The mixed layer depth is calculated taking into account both wind and thermally driven mixing \citep{Mironov2008}, that follows 
from the volumetric heating of the water column by solar radiation. 
The structure of the stratified layer between the upper mixed-layer and the lake bottom, the thermocline,
reposes on the concept of self-similarity (assumed shape) of a temperature-depth curve. 
The same concept is used to describe the temperature structure of the lake ice. 
Considering fresh water density as

%(1)
\begin{equation}
\rho_{\rm w} =\rho_{\rm r} \left[ 1 - \frac{1}{2} \alpha_{\rm tl} (T_{\rm MN} - T_{\rm MNr})^2 \right]   \label{P4LakeEq1}
\end{equation}

where $\rho_{\rm w}$ is the water density, $\rho_{\rm r}=999.98 kg m^{-3}$ is the maximum density of fresh water at the temperature $T_{\rm MNr} = 277.13 K$, 
and $\alpha_{\rm tl} = 1.6509~K^{-2}$ is an empirical coefficient. According to \ref{P4LakeEq1},
the thermal expansion coefficient $\alpha_{\rm tl}$ and a buoyancy parameter $\beta$ depend on the water temperature according to:

%(2)
\begin{equation}
\beta (T_{l})=g \alpha_{T_l} = g~ a_{T} (T_{l}(z)-T_{r})
\label{P4LakeEq2}
\end{equation}
where $g$ is the acceleration due to gravity.
The following two-layer parametric representation of the evolving temperature profile $T_{l}(z) $is adopted

%(3)
\begin{equation}
T_{l}(z) = \left\{ \begin{array}{rl}
    T_{ML} &\mbox{ if $0 < z \le d_{ML}$} \\
T_ {ML}-(T_{ML}-T_{BL})\Phi_{T_l}(\zeta) &\mbox{ if $d_{ML} < z\le d_{l}$}
       \end{array} \right.
       \label{P4LakeEq3}
\end{equation}

as a function of depth $z$, where $T_{ML}$ is the temperature of the upper mixed layer of depth $d_{ML}$. 
$T_{BL}$ is the temperature at the bottom of the water body ($z = d_l$), and
$\Phi_{T}=[T_{ML}-T_{l}(z,t)]/[T_{ML}-T_{BL}]$ is a dimensionless function of dimensionless depth
$\zeta=[z / d_{ML}(t)]/[d_{l} / d_{ML}(t)]$ that satisfies the boundary conditions $\Phi_{T}(0) = 0$ and $\Phi_{T}(1) = 1$.

The mean temperature of the water column $T_{MN}$ is obtained integrating \ref{P4LakeEq3}, and relates $d_{ML}$ , $d_{l}$, $T_{ML}$, $T_{BL}$ as

%(4) %(5)
\begin{equation}
T_{MN}=T_{ML}-\xi_{ML} \Big( 1- \frac{d_{ML}}{d_{l}} \Big) (T_{ML}-T_{BL})
\label{P4LakeEq5}
\end{equation}

where $\xi_{ML}$ is the shape factor with respect to the temperature profile in the thermocline.

The equation for the mean temperature of the water column
is obtained by integrating the one-dimensional heat transfer equation over $z$ from $0$ to $d_l$, and it reads as

\begin{equation}
d_{l} \frac{d T_{MN}}{dt}=\frac{1}{\rho_w c_w} \Big[ G_s + I_s - G_{BL} - I(d_l) \Big]
\label{P4LakeEq7}
\end{equation}

where $c_{w} = 4.2 ~10^{3} J~kg^{-1} K^{-1}$ is the specific heat of water, $G_s$ and $I_s$ are the values at the surface of the vertical heat flux 
and of the heat flux due to solar radiation, respectively, and $G_{BL}$ is the heat flux through the lake bottom.
The radiation heat flux $I$ that penetrates into the water is the surface value of the incident solar radiation flux from the atmosphere
multiplied by $1-\alpha_w$, where $\alpha_w$ is the albedo of the water surface with respect to solar radiation.
The surface flux $G_{s}$ is a sum of the sensible and latent heat fluxes and the net heat flux due to long-wave radiation at the air-water interface.
The equation of heat budget of the mixed layer reads as

%(8)
\begin{equation}
d_{ML} \frac{d T_{ML}}{dt}=\frac{1}{\rho_w c_w} \Big[ G_s + I_s - G_{ML} - I(d_{ML}) \Big]
\label{P4LakeEq8}
\end{equation}

where $G_{ML}$ is the heat flux at the bottom of the mixed layer, and assumes two variants for stationary/retreating or deepening mixed layers as detailed in \citet{Mironovetal2010} .
In the case of the mixed-layer being stationary or retreating, $d~d_{ML}/dt \le 0$, the bottom temperature is

%(9)
\begin{equation}
\frac{d T_{BL}}{dt}=0
\label{P4LakeEq9}
\end{equation}

In the case of the mixed-layer deepening, $d~d_{ML}/dt > 0$, the following equation is used:
%(10)
\begin{eqnarray}
&\frac{1}{2}&(d_{l}-d_{ML})^2 \frac{d T_{ML}}{dt} - \frac{d}{dt} \Big[C_P (d_{l}-d_{ML})^2 (T_{ML}-T_{BL})\Big] = \nonumber \\
&=& \frac{1}{\rho_w c_w} \Big[C_G (d_{l}-d_{ML})(G_{ML}-G_{BL}) + (d_{l}-d_{ML}) I_{ML}-I_{ML_BL} \Big]
\label{P4LakeEq10}
\end{eqnarray}

Where $C_P$ and $C_G$ are dimensionless parameters. If $d_{ML}=d_l$ both $T_{ML}$ and $T_{BL}$ are equal to $T_{MN}$ and the water column profile is isothermal.
The depth of a stably or neutrally stratified wind-mixed layer is determined from:

%(15)
\begin{equation}
\frac{d~d_{ML}}{dt} = \frac{d_{EQ}-d_{ML}}{d_{EQ}} C_{rh} u^*
\label{P4LakeEq15}
\end{equation}

where $u^*$ is the surface friction velocity, $C_{rh}=0.03$ is a dimensionless constant, and $d_{EQ}$ is a mixed-layer depth in radiative equilibrium \citep{Mironovetal2010}.

\subsubsection{Lake freezing and melting}
The parameterisation of lake freezing is essential for treating inland water bodies in global models as mid-latitude and high-altitude lakes are 
often subject to extensive frozen periods especially when characterized by a shallow depth.

In the current implementation the snow accumulation over ice is not allowed and the existence of fractional ice condition on a given grid-point is not contemplated.
The approach used to describe the temperature structure of the lake ice is conceptually similar to the approach used to describe the temperature structure of the lake thermocline.
The following parametric representation of the evolving temperature profile within the ice is adopted (cf. \ref{P4LakeEq3}):

%(20)
\begin{equation}
T_{ice}=T_f-(T_f-T_{sk,i})\Phi_{ice}(\zeta_{i}) \mbox{ if $-d_{ice} < z \le 0$}
\label{P4LakeEq20}
\end{equation}
where $z$ is the vertical co-ordinate (positive downward) with the origin at the ice-water interface, $d_{ice}(t)$ is the ice thickness,
$T_f = 273.15~K$ is the fresh-water freezing point, and $T_{sk,i}(t)$ is the temperature at the ice upper surface. 
The dimensionless universal function $\Phi_{ice}=[T_{f}-T_{ice}(z,t)]/[T_{f}-T_{sk,i}]$
of dimensionless depth $\zeta_i=z/d_{ice}(t)$ satisfies the boundary conditions $\Phi_{ice}(0) = 0$ and $\Phi_{ice}(1) = 1$.
The equation of the heat budget of the ice layer is 

%(21)
\begin{eqnarray}
&\frac{d}{dt}& \left\{ \rho_{ice} c_i d_{ice} [T_f - C_i (T_f-T_{ice})] \right\}- \rho_{ice} c_i T_{ice} \frac{d~d_{ice}}{dt} = \nonumber \\
&=& G_{s} +I_{s} -I(0)+ \kappa_i \frac{T_f-T_{ice}}{d_{ice}}\Phi_{ice}^{'}(0)
\label{P4LakeEq21}
\end{eqnarray}

where $\rho_{ice} =9.1~10^{2} kg~m^{-3}$, $c_i =2.1~10^{3} J~kg^{-1}~K^{-1}$ and
$\kappa_{i} =2.29~J~m^{-1} s^{-1} K^{-1}$ are the ice mean property for density, specific heat and heat conductivity, respectively.
$G_s$ and $I_s$ are the values of $G$ and $I$, respectively, at the ice upper surface $z = d_{ice}(t)$, and $\Phi_{ice}^{'}(0) = d\Phi_{ice}/d\zeta_i$ at $\zeta_i = 0$. 
The radiation heat flux $I_s$ that penetrates into the ice medium is the surface value of the incident solar radiation flux from the atmosphere
multiplied by $1-\alpha_i$, where $\alpha_{i}$ is the ice surface albedo with respect to solar radiation and
$C_i$  dimensionless parameter representing the shape factor with respect to the temperature profile within the ice.
%
%(22)
%\begin{equation}
%\label{P4LakeEq22}
%\end{equation}

Equation \ref{P4LakeEq21} serves to determine $T_{ice}$ when this temperature is below the freezing point $T_f$,
i.e. when no melting at the ice upper surface takes place. During the ice melting from above, $T_{ice}$ remains equal to $T_f$.
During the ice growth or ice melting from below (these occur as $T_{ice} < T_f$), the ice thickness is computed as:

%(23)
\begin{equation}
L_f \frac{d (\rho_{ice} d_{ice})}{dt} = G_w + \kappa_{ice} \frac{T_f-T_{ice}}{d_{ice}} \Phi_{ice}^{'}(0) 
\label{P4LakeEq23}
\end{equation}
where $L_f = 3.3~10^{5}~J~kg^{-1}$ is the latent heat of fusion, and $G_w$ is the heat flux in the near-surface water layer just beneath the ice. 
If the right-hand side of \ref{P4LakeEq23} is negative (this may occur due to a negative $G_w$), ice ablation takes place.
During the ice melting from above, the following equation is used:

%(24)
\begin{equation}
L_f \frac{d(\rho_{ice} d_{ice})}{dt} = G_w - (G_s+I_s) + I(0)
\label{P4LakeEq24}
\end{equation}
that holds as the atmosphere heats the ice upper surface and $T_{ice}$ is equal to $T_f$.
The evolution of the temperature profile beneath the ice is described as follows.
The temperature at the ice-water interface is fixed at the freezing point, $T_{ML}=T_f$. 
The mean temperature of the water column is computed from \ref{P4LakeEq23},
where $G_s$ and $I_s$ are replaced with $G_w$ and $I(0)$, respectively. 
If the bottom temperature is less than the temperature of maximum density, $\rho$,
the mixed-layer depth and the shape factor with respect to the temperature profile in the thermocline are kept unchanged,
$d~d_{ML}/dt = 0$ and $\xi/dt = 0$, and the bottom temperature is computed from \ref{P4LakeEq5}.
If the entire water column appears to be mixed at the moment of freezing, i.e. $d_{ML} = d_l$ and $T_{ML} = T_{BL}$,
the mixed layer depth is reset to zero, $d_{ML} = 0$, and the shape factor is reset to its minimum value, 
$\xi_{ML} = \xi_{min}$. 

As the bottom temperature reaches the temperature of maximum density, its further increase is prevented and is kept constant equal to $\rho_r$. 
If $d_{ML} > 0$, the shape factor $\xi_{ML}$ is kept unchanged, and the mixed-layer depth is computed from \ref{P4LakeEq5}. 
As the mixed-layer depth approaches zero, \ref{P4LakeEq5} is used to compute the shape factor
$\xi_{ML}$ that in this regime increases towards its maximum value,
$\xi_{ML} = \xi_{max}$ (the values of $\xi_{min}=0.5$ and $\xi_{max}=0.8$ are used).
If $d_{ML}=0$, the heat flux from water to ice is estimated from

%(25)
\begin{equation}
G_w=- \kappa_w \frac{T_b-T_{ice}}{d_{l}} {max}[1,\Phi_{w}(0)]
\label{P4LakeEq25}
\end{equation}

where $\kappa_w = 5.46~10^{-1}J~m^{-1}~s^{-1}~K^{-1}$ is the molecular heat conductivity of water, and
$\Phi_{w}^{'}(0) = d\Phi_{w}/d\zeta$ at $\zeta = 0$. 
Note that as soon as a mixed layer is formed ($d_{ML} > 0$) the above water flux is neglected ($G_w = 0$).

The shape factor with respect to the temperature profile in the thermocline is computed from

%(26)
\begin{equation}
\frac{d{\xi}_{\rm ML}}{dt}={sign}(d~d_{ML}/dt) \frac{\xi_{max}-\xi_{min}}{t_{rc}}
\label{P4LakeEq26}
\end{equation}
where $sign$ is the sign function $[sign(x)=x>0]$. The shape factor $\xi_{\rm ML}$ evolves towards its maximum value during the mixed-layer deepening,
and towards its minimum value during the mixed-layer stationary state or retreat. 
The adjustment occurs on a relaxation time scale $t_{rc}$ estimated using the buoyancy frequency in the thermocline \citep{Mironovetal2010}.

%%(27)
%\begin{equation}
%\label{P4LakeEq27}
%\end{equation}
%where $Crc = 0.003$ is a dimensionless constant, and the mean-square buoyancy frequency in the thermocline is given by \ref{P4LakeEq18}. 
%Eqs. 26 and 27 are used during the period of open water. During the period of ice cover, a different procedure is used as outlined above.
%The dimensionless parameter $C$ defined through \ref{P4LakeEq11} is given by

%(28)
%\begin{equation}
%\label{P4LakeEq28}
%\end{equation}
%and the quantity $b8(0)$ that enters \ref{P4LakeEq25} is given by

%(29)
%\begin{equation}
%\label{P4LakeEq29}
%\end{equation}
%The shape factor with respect to the temperature profile within the ice is computed from

%(30)
%\begin{equation}
%\label{P4LakeEq30}
%\end{equation}
%where $himax = 3 m$ and $i = 2$. The quantity $i(0)$ that enters Eqs. 21 and 23 is given by

%(31)
%\begin{equation}
%\label{P4LakeEq31}
%\end{equation}
The exponential approximation of the decay law for the flux of solar radiation is used:

%(32)
\begin{equation}
I=I_s\Sigma_{k=1}^n a_k e^{-\gamma_k(z+d_{ice})}
\label{P4LakeEq32}
\end{equation}
where $I_s$ is the surface value of the incident solar radiation heat flux multiplied by $1-\alpha$ with $\alpha$ 
being the albedo of the water surface or of the ice surface with respect to solar radiation, $n$ is the number of wavelength bands,
$\alpha_k$ are fractions of the total radiation flux for different wavelength bands, and $a_k$ are attenuation coefficients for different bands.
The attenuation coefficients are piece-wise constant functions of $z$, i.e. they have different values for water and ice but remain constant within each media.
The following parametrisation of the ice surface albedo with respect to solar radiation is adopted:

%(33)
\begin{equation}
\alpha_i={\alpha_i}^{max}-({\alpha_i}^{max}-{\alpha_i}^{min}) e^{-C_{\alpha}(T_f-T_{ice})/T_f}
\label{P4LakeEq33}
\end{equation}
where ${\alpha_i}^{max} = 0.7$ and ${\alpha_i}^{min}= 0.4$ are maximum and minimum values of the ice albedo, respectively, and $C_{\alpha}= 95.6$ is a fitting coefficient. 
The presence of snow over lake ice and its seasonal changes are parametrized in \ref{P4LakeEq33} as a function of $T_{ice}$. 
During the melting season, the ice surface temperature is close to the fresh-water freezing point.
The presence of wet snow, puddles, melt-water ponds is again parametrized implicitly and results in a decrease of the area-averaged surface albedo. 
The water surface albedo with respect to solar radiation, $\alpha_w=0.07$, is assumed to be constant.
\
\unskip

\section{Vegetation growth model}
\label{appendixB}

The vegetation growth is the result of environmental factors affecting the accumulation of carbon by assimilating atmospheric $CO_{\rm2}$ and senescence which results from deficit of photosynthesis. The LAI is then derived from the leaf biomass $B$ as a result of the net assimilation of the canopy $A_{nI}$ (equation \ref{P4eq7.101}) and is obtained from $B$ as:

\begin{equation}
LAI=\frac{B}{\alpha_{B}},
\end{equation}

where ${\alpha_{B}}$ is the inverse of the specific Leaf area (SLA) being the ratio of the green leaves biomass to LAI (equation \ref{eqLAI_SLA}) and the leaf biomass $B$ is obtained from the following differential equation:

\begin{equation}
dB=\frac{M_{\rm c}}{P_{\rm c}M_{\rm CO_{\rm2}}}{A_{nI}}{dt}-Bd(t/\tau) 
 \label{eqBio}
\end{equation}

$P_{\rm c}$ here is the proportion of carbon in the dry plant biomass, and $M_{\rm c}$ and $M_{\rm CO_{\rm2}}$ are the molecular weights of carbon and $CO_{\rm2}$ , respectively. The mortality increment term of equation (\ref{eqBio}) represents an exponential extinction of $B$ characterized by a time-dependent effective life expectancy (expressed in units of days):

\begin{equation}
\tau(t)=\tau_{M} \frac{A_{n f m}(t)}{A_{n, \max }}
\end{equation}

This development relies on having a simple allocation scheme for the total nonwoody aboveground biomass $B_{T}$ which is considered as the sum of the leaf biomass $B$ and the structural biomass $B_{s}$. In the nitrogen decline (N-decline) model, $B_{T}$ is related to the the leaf biomass $B$ by the allometric equation following \citet{CalvetSousana2001}:

\begin{equation}
B_{T}=\left(\frac{B}{c}\right)^{1 /(1-a)},
\label{eqNdecline}
\end{equation}

where $c$ and $a$ are constant parameters ($c = 0.754 , a = 0.38$).

The B-decline term of equation (\ref{eqBio}) is split into a mortality and a storage term $(M_{ B}$ and $S_{ B}$ , respectively):
\begin{equation}
B d(t / \tau)=M_{B}+S_{B}
\end{equation}

In equation (\ref{eqBio}) when net assimilation exceeds the B-decline term a growing phase is happening with positive values of $dB$. In this case the plant N-decline model is triggered and the aboveground biomass $B_{ T}$ is derived from $B$ using equation \ref{eqNdecline}. While the structural biomass $B_{s}$ is expressed as the difference between the total aboveground biomass $B_{ T}$ and the leaf biomass $B$. The mortality of the structural biomass $M_{B s}$ is then computed by:

\begin{equation}
M_{B s}=B_{s} d\left(t / \tau_{M}\right)
\label{eqStrucMort}
\end{equation}
In addition, a loss of carbon is simulated through the structural biomass respiration with: 

\begin{equation}
R_{B s}=\eta_{R} B_{S} Q_{10}^{(T-25) / 10} d t
\label{eqStrucResp}
\end{equation}

where $dt$ is one day step, $T$ the surface temperature (in $^{\circ}C$), $R$ a respiration coefficient of 1\% of $B_{ s}$ per day, and $Q_{ 10} = 2$.

The storage $S_{ B}$ term is calculated as the residual of the structural biomass budget:

\begin{equation}
S_{B}=dB_{s} +M_{B_{s}}+R_{B_{s}} 
\label{eqStorgB}
\end{equation}

An alternative formulation of B-decline is used when the storage exceeds the mortality of active biomass $M_{B}$. $B_{ T}$ is recalculated assuming that there is no loss of leaf biomass outside the plant system, during the considered time step, $M_{B} = 0$ and $dB_{ T}$ is the difference between the daily net assimilation and the mortality and respiration losses of structural biomass following equations \ref{eqStrucMort} and \ref{eqStrucResp}.
The leaf biomass $B$ is derived from $B_{ T}$ using the N-decline equation, and $B_{ s}$ is the difference between the two terms, an updated value of $S_{B}$ is then computed by equation (\ref{eqStorgB}. When $db$ becomes negative,the vegetation is senescent and the N-decline equation is not valid anymore; In such configuration, the structural biomass reservoir $B_{ s}$ becomes independent from the leaf biomass reservoir $B$ to which a nil value term is assigned and the mortality and respiration losses are applied to $B_{ s}$ using equations \ref{eqStrucMort} and \ref{eqStrucResp}.

\begin{table}
\centering
\caption{Vegetation growth model parameters values as specified in ECLand model. The vegetation types foloows the BATS classification
\label{TabGrowth}}
\small

\begin{tabular}{p{0.25in} p{1.2in} p{0.8in} p{0.55in} p{0.55in}}
\toprule
 & & \bf{\textit{e} & \bf{\textit{f} & \bf{\textit{N${}_{\rm l}$ }}  }}\\  
%
\bf{\textit{Type}} & \bf{\textit{Vegetation type}} & [m${}^{2}$ kg${}^{-1}\%^{-1}$] & 
[m${}^{2}$kg${}^{-1}$] & [\%]  \\  \hline \hline 

1 & Crops, mixed farming & 3.79 & 9.84 & 1.3  \\ 
2 & Short grass & 5.56 & 6.73 & 1.3  \\ 
3 & Evergreen needleleaf  & 4.85 & -0.24 & 2.8  \\
4 & Deciduous needleleaf  & 4.85 & -0.24 & 2.8  \\ 
5 & Deciduous broadleaf  & 4.83 & 2.53 & 2.0  \\ 
6 & Evergreen broadleaf  & 4.83 & 2.53 & 2.5 \\
7 & Tall grass & 7.68 & -4.33 & 1.3 \\ 
8 & Desert & - & - & -  \\
9 & Tundra & 5.56 & 6.73 & 1.3  \\ 
10 & Irrigated crops & 3.79 & 9.84 & 1.3  \\ 
11 & Semidesert & 5.56 & 6.73 & 1.3  \\ 
12 & Ice caps and glaciers & - & - & -  \\
13 & Bogs and marshes & 5.56 & 6.73 & 1.3  \\ 
14 & Inland water & - & - & -  \\ 
15 & Ocean & - & - & -  \\ 
16 & Evergreen shrubs & 7.68 & -4.33 & 1.3  \\ 
17 & Deciduous shrubs & 5.56 & 6.73 & 1.3  \\ 
18 & Mixed forest- Wood & 4.85 & -0.24 & 2.8  \\ 
19 & Interrupted forest & 4.85 & -0.24 & 2.8  \\ 
20 & Water -land mixtures & 5.56 & 6.73 & 1.3  \\ %\hline 
\bottomrule
\end{tabular}
\end{table}


\

\unskip
%\setcounter{table}{0}
%\setcounter{figure}{0}



\section{Ancillary data}
\label{climfield_apdx}
The land surface model uses several ancillary climate fields, that do not depend on initial condition or forecast step. These fields are the land-sea-mask, the lake cover and depth, soil type, surface albedo, Leaf area index and  vegetation type and vegetation cover. Most fields are constant in time but surface albedo and LAI are monthly specified to describe the seasonal cycle. A list of these fields is provided in Table \ref{tab_clim}. 

\begin{table}[htbp!]
\caption{Climatological fields of the land surface model}
\centering
\begin{adjustbox}{width=1\textwidth}
\small
\begin{tabular}{llll}
\hline
\textbf{Field} & \textbf{variability} & \textbf{Data} & \textbf{reference} \\ \hline
Land sea mask & static & Global Land Cover service (GLOBCOVER 2006) &  \citet{Arinoetal2007} \\ 
Orography & static & Shuttle Radar Topography Mission  (SRTM 30) & \citet{Farretal2007} \\ 
Land Use/Land Use & static & Global Land Cover Characteristics (GLCCv1.2) &  \citet{Lovelandetal2000} \\ 
Leaf Area Index & Monthly & Moderate Resolution Imaging Spectroradiometer (MODIS collection 5) &  \citet{Mynenietal2002} \\ 
Albedo & Monthly & Moderate Resolution Imaging Spectroradiometer (MODIS ) & \citet{ Schaafetal2002} \\ 
Soil Type & static & Digital Soil Map of the World (DSMW) & \citet{fao03} \\ 
Lake depth & static &  global lake database (GLDB) & \citet{Kourzeneva2010}, \citet{Balsamoetal2012} \\ \hline
\end{tabular}
\end{adjustbox}
\label{tab_clim}
\end{table}


\section{New surface water data}

 \textbf{Glacier mask} \\
 The glacier fractional mask is based on:
\begin{itemize}
\item Svalbard glacier data from the Norwegian Institute \citep{SvalbardGlacier2014}, shape file nominal resolution 1:100 000;
\item Iceland glacier data from the Icelandic Meteorological Office (Palmason, personal communication 2019), data for the year 2017, nominal resolution 100 m;
\item Greenland glacier data from The Greenland Ice Mapping Project (GIMP) (Howat et al., 2014), data based on years 1999-2001, nominal resolution 15 m,
\item Antarctica glacier data is based on observations recorded by the CryoSat-2 satellite (Slater et al., 2018), data based on years 2010-2016, nominal resolution 1 km (30”);
\item Global Land Ice Measurements from Space (GLIMS) project (GLIMS and NSIDC, 2005, updated 2013), data for the year 2017, polygonal boundary dataset is used for the rest of the globe glaciers.
\end{itemize}

 \textbf{land-water and lake fraction mask} \\
In addition to the glacier mask data (above), the upgraded land-water mask main source of data is the Joint Research Centre (JRC) Global Surface Water 30 m (1”) resolution dataset. An expert system (Pekel et al., 2016) to separate open water from other surfaces over the period from March 1984 to 2018 was used to create this data based on satellite thermal images. This relies on the contrasting spectral properties of water and other features (e.g. snow) from the Landsat 5, 7 and 8 sensors' six visible, near and shortwave infrared channels.
Monthly Water History for years 2010-2018 (Pekel et al., 2016) was also used to confirm permanent water pixels. The global digital elevation map based on 1987-2017 data with 3” nominal resolution (~90m) MERIT DEM (Yamazaki et al., 2017) was also used to fill large permanent open water areas that are away from coastal points (e.g. middle of Caspian Sea, middle of the ocean with no islands, etc.) since Global Surface Water data is masked in these areas. These data are merged to create a consistent land-water mask at 1km resolution. The orographic data generated at 1km and used in the current IFS (\citet{IFSdocCY47R1}, Ch11) is then used to fill additional points/islands that are missing from the previous steps and to ensure consistency with the new mask. The lake fraction mask is then generated by a flooding algorithm (Choulga et al., 2019) at the same resolution. \\
 \textbf{Lake depth} \\
The lake depth is specified for all land points mainly following the Global Lake DataBase v3, (GLDBv3; Kourzeneva, 2010 and Choulga et al., 2014).  GLDBv3 is based on a lake mask extracted from Ecoclimap2 ecosystem dataset (Champeaux et al., 2004) and the depths are set according to:
\begin{itemize} 
\item	mean depth from in situ data of about 14500 lakes, 
\item	6 large lakes bathymetry data, 
\item	all navigable and most of non-navigable Finish lakes bathymetry data, 
\item	indirect mean depth estimates based on lake geological origin (Choulga et al., 2014), 
\item	default depths set to 5 m for saline lakes, 10 m for natural freshwater lakes, artificial reservoirs and unidentified waterbodies, 50 m for crater/caldera lakes.
\end{itemize}
In addition, the Great Lakes bathymetry is set from the Global Relief Model ETOPO1 (Amante and Eakins, 2009). The Caspian Sea, the Azov Sea and the ocean bathymetry are obtained from the General Bathymetric Chart of the Oceans (GEBCO) data (Weatherall et al., 2015).
These data are merged into a 1km regular lat/lon file and the interpolation to a target grid with coarser resolution is based on mode (most common) value, land water mask (ocean or inland water) and data source (availability of in situ measurements, indirect estimates or default values) types. This new approach helps to enhance in situ value contribution at model high resolutions where the lake signal is most evident. Further details are specified in Choulga et al. (2019).
Figure \ref{figlakeD} shows the difference between the new and the current operational lake depth 

\begin{figure}[hbt!]
    \centering
    \centering {{\includegraphics[width=15cm, trim={0 0 0 25cm},clip]{Figures/Diff_dl_1279_4.png}}}

    \caption{Water depth difference: GLDBv3 based - Current operational map in [m]}
    \label{figlakeD}
\end{figure}


\subsubsection{Data merging and consistency}

To preserve The JRC high resolution information all calculations were done at 30 m and then aggregated to 1 km. Due to high volume of data needed to be processed, the Google Earth Engine web-based code editor was used (https://earthengine.google.com/, \citet{gorelick2017google} ).
In the case of the glacier mask all data were interpolated on a 30 m latitude/longitude regular grid. Grid-cell values vary from 0 (no glacier at all) to 1 (fully covered by glacier) depending on what size of grid-cell is covered by initial data. Data from several regional sources for Svalbard, Iceland, Greenland and Antarctica are merged together and completed globally with GLIMS project data. In addition, Antarctica region was checked for areas with no glacier further from the coast; these grid-cells were forced to be fully covered by glacier, as they might cause unrealistic values in the atmospheric model (e.g. ERA5).
The land-water mask was then created according to the following steps:
\begin{enumerate}
\item from JRC Global Surface Water Mapping Layers, v1.1, the permanent water layer (PWL) is created according to the change classification flags;
\item from JRC Global Surface Water Metadata, v1.1, a water detection layer (WDL) is created based on  a ratio of 0.7 or more of the number of valid water observations to the total number of observation;
\item from JRC Monthly Water History, v1.1, 108 monthly maps (January 2010 to December 2018) a  water land layer (WLL) is created by classifying each grid-cell as (i) “land” if grid-cell was defined as NotWater more than 0 times and as Water 0 times, (ii) “water” if grid-cell was defined as NotWater 0 times and as Water 0 or more times, and (iii) "nodata” for the remaining cases;
\item from MERIT DEM a global template layer (GTL) is created where grid-cells are classified only as fully land or fully no land;
\item in addition, a Black sea regional mask (BM) is created;
\item all previously generated layers including the glacier mask are merged at 30 m resolution according to Figure \ref{figlsmProc} and combined 30m layer is created;
\item  The combined 30m layer is aggregated to 1 km by simple averaging and a fractional combined 1km layer is created;
\item The combined 1km layer is then checked against the global template (GTL) to delete erroneous points or to reshape existing islands (e.g. satellite blobs in the middle of the ocean where number of observations is very low due to geographical location or meteorological conditions) and finally, the upgraded fractional land-water mask is created.
\end{enumerate}

\begin{figure}[H]
    \centering
    \centering {{\includegraphics[width=5cm]{Figures/LSM_creation_scheme_new.jpg}}}

    \caption{Land-Water Mask processing chart}
    \label{figlsmProc}
\end{figure}


\section{Urban Parametrisation}

 The complete list of variables required for the urban scheme are based on previous similar schemes  (\citep{masson2000physically}, \citep{oleson2008urban}, \citep{porson2010implementation}) and listed in Table \ref{urbTab}


\begin{table}[htbp!]
\caption{Urban scheme parameters}
\begin{tabular}{|l|c|c|c|}
\hline
Street-to-Building Ratio & N & Y & 0.1-1.0 \\ \hline
\textbf{Radiation Variables} & \multicolumn{1}{l|}{\textbf{}} & \multicolumn{1}{l|}{\textbf{}} & \multicolumn{1}{l|}{\textbf{}} \\ \hline
Wall Albedo & N & Y or N & 0.04-0.4 \\ \hline
Roof Albedo & N & Y or N & 0.04-0.4 \\ \hline
Street Albedo & N & Y or N & 0.02-0.2 \\ \hline
Wall Emissivity & N & Y or N & 0.80-0.95 \\ \hline
Roof Emissivity & N & Y or N & 0.80-0.95 \\ \hline
Street Emissivity & N & Y or N & 0.90-0.99 \\ \hline
Fraction of Radiation Scattered in Canyon & N & N & 0.3 \\ \hline
\multicolumn{ 4}{|l|}{Canyon Albedo} \\ 
Canyon Emissivity & Y & N & - \\ \hline
Effective Urban Albedo & Y & N & - \\ \hline
Effective Urban Emissivity & Y & N & - \\ \hline
\textbf{Roughness Variables } & \multicolumn{1}{l|}{\textbf{}} & \multicolumn{1}{l|}{\textbf{}} & \multicolumn{1}{l|}{} \\ \hline
Roughness Length for Momentum of Building Materials & N & N & 1x10-4-1x10-6 \\ \hline
Drag Coefficient for Canyon & N & N & 1.2 \\ \hline
\multicolumn{ 4}{|l|}{Alpha Coefficient for Canyon} \\ 
Displacement Height & N & N & 1.4 \\ \hline
Roughness Length for Momentum of Urban Tile & Y & N & - \\ \hline
Roughness Length for Heat of Urban Tile & Y & N & - \\ \hline
\textbf{Thermal Variables} &  &  &  \\ \hline
Volumetric Heat Capacity of Wall (Jm-3K-1) & N & Y or N & 5.0x105-5.0x107 \\ \hline
Volumetric Heat Capacity of Roof (Jm-3K-1) & N & Y or N & 3.0x105-3.0x107 \\ \hline
Volumetric Heat Capacity of Street (Jm-3K-1) & N & Y or N & 5.0x105-5.0x107 \\ \hline
Thermal Conductivity of Wall (Wm-1K-1) & N & Y or N & 0.2-2.0 \\ \hline
Thermal Conductivity of Roof (Wm-1K-1) & N & Y or N & 0.1-1.0 \\ \hline
Thermal Conductivity of Street (Wm-1K-1) & N & Y or N & 0.2-2.0 \\ \hline
\multicolumn{ 4}{|l|}{Canyon Heat Capacity} \\ 
Canyon Thermal Conductivity & Y & N & - \\ \hline
Effective Urban Heat Capacity & Y & N & - \\ \hline
Effective Urban Conductivity & Y & N & - \\ \hline
\multicolumn{1}{|c|}{\textbf{Hydrology Variables}} &  &  &  \\ \hline
Urban Run-off Rate & N & Y or N & 0.8 \\ \hline
Alpha-factor (Van Genuchten) & N & N & 0.765 \\ \hline
Lambda-factor (Van Genuchten) & N & N & 35.2 \\ \hline
Saturation limit & N & N & 0.01-1.0 \\ \hline
\end{tabular}
\label{urbTab}
\end{table}



\subsubsection{Urban Radiation}

The initial urban radiation developments, described here, will be superseded by a more detailed description of
urban-atmosphere interactions when the required mapped parameters become available \citep{hogan2019flexible}. Our scheme
here is based upon previous schemes developed for the Climate Land Model - Urban (CLMU, \citep{oleson2008urban}) and the
Met Office - Reading Urban Surface Exchange Scheme (MORUSES, \citep{porson2010implementation}). The solar-radiation scheme for the
roof fraction of the tile follows a surface albedo simple incoming/reflected exchange. For the canyon fraction of the
tile the impact of shadowing based on the solar zenith, building height and aspect ratio is considered along with the
albedo of the wall and street. The canyon bulk albedo,  $\alpha _{\mathit{can}}$, which is the same for both direct and diffuse solar radiation, is derived using

\begin{equation}
\alpha _{\mathit{can}}=1-\left[\left(\left\{1-\chi_f\right\}\chi_r+Q_{S,r}\right)+\frac{2H} W\left(\left\{1-\chi_f\right\}\chi_w+Q_{S,w}\right)\right]
\label{equrbrad}
\end{equation}

where  $\chi _f$ \ is the fraction of solar radiation scattered by the sky (estimated at 0.3),  $\chi _r$ \ and  $\chi
_w$ \ are the shadowing factors of street and walls following \citet{masson2000physically}, Q is the net radiation exchange following \citet{harman2004radiative} and $\frac H W$ \ is the building height to street width ratio. The solution is reached by first calculating the shadowing factors, taken from equations 13 and 14 of  \citet{masson2000physically}. Next, to satisfy equation 11 of \citet{porson2010implementation}, a matrix of exchange coefficients is calculated from equations 2 to 5 of \citet{harman2004radiative}. Finally, the incoming/outgoing exchange-diffuse radiation is calculated using equations 12-14 of \citet{porson2010implementation}.

For bulk emissivity ($\varepsilon $) a similar approach is adopted where net solar radiation is replaced with long-wave
radiation using material emissivity. The bulk emissivity for the canyon is influenced by the building height to canyon
width ratio but in a different way to albedo. The code has been written to calculate bulk emissivity based on  $\frac H
W$ \ of the canyon taken from the appendix of \citet{porson2010implementation} and such that it would satisfy the net long-wave
radiation,  $R_{L,c}$ \ using

\begin{equation}
\centering
 R_{L,c}=\varepsilon (\mathit{LW}_{\downarrow }-\sigma T_c^4),
\label{equrbrad2}
\end{equation}

where  $\varepsilon $ \ is the canyon bulk emissivity,  $\mathit{LW}_{\downarrow }$ \ is the downward long-wave
radiation,  $\sigma $ \ is the Stefan-Boltzmann constant and  $T_c$ \ is the canyon temperature, which is assumed to
equal the surface temperature. The long and shortwave radiation routine used here assumes building azimuthal
orientation is random and accounts for all possible directions equally. Here, as a simplification, albedo and
emissivity are assumed to be spectrally constant in both long and shortwave,  additionally gas absorption/emission in the
longwave is ignored. Future work will seek to address this based on the work of \citet{hogan2019flexible}.

\subsubsection{Urban Roughness}

The roughness length for momentum and heat exchange with the atmosphere is calculated based on the \citet{macdonald1998improved} formulation, where roughness is a function of  $\frac H W$ \ and the ratio of street width to building width (
$\frac R B$). The estimated displacement height is calculated based on a calibrated coefficient for staggered arrays (A
= 4.43) and drag coefficient,  $\beta $\ = 1, taken from \citet{macdonald1998improved}. These formulations are


\begin{equation}
\begin{split}
 D=1-\frac R B (A^{\left(\frac R B-1\right)}), \\
 z=(\beta \ast \left(1-D\right)\frac{\frac{\frac H W}{2(2\ast \tan ^{-1}\left(1\right))}\ast R}
B/k^2), \\
 \mathit{Zm}=H\ast \left(1-D\right)\ast e^{-z}, \\
 \end{split}
\label{equrbrad2}
\end{equation}


where  $k$ \ is the Von Karman constant (0.4) and  $\mathit{Zm}$ \ is the canyon roughness length for momentum.


\bigskip

The approach to estimate roughness length for heat uses morphometric methods relating to aerodynamic parameters based on
empirical relations. We use previously defined flow regimes of isolated roughness, wake interference and skimming \citep{harman2004scalar}. These regimes separate the infinite canyon into recirculation and ventilated regions. Full
details of the roughness length calculation can be found in \citet{harman2004scalar}. As with the radiation, a random
orientation assumption is made for the infinite length canyon model used. This formulation is dependent on bulk
resistance of both the canyon and roof (see equation 18 and 19 of \citet{porson2010implementation} ). This results in a roughness length for heat which is several order magnitudes smaller than the roughness length for momentum in agreement with
similar studies (e.g. \citep{kanda2005simple}).

\subsubsection{Urban Thermal Properties}

Urban thermal properties are implemented following a similar methodology to the radiation, whereby individual properties
are defined for the urban wall, roof and street structure. These combine to provide bulk properties weighted by the
fraction of the urban tile covered by each. These bulk variables, thermal conductivity and volumetric heat content, are
used to compute the transfer of heat from the atmosphere and between upper most and second soil layer. Further updates
aim to replace the four soil layers with four defined urban properties (wall, roof, street and internal) including
interactions between each medium. For accurate canyon representation additional properties are likely to be required
representing canyon temperature for heat exchange.

\subsubsection{Urban Hydrology}

As with previous properties, the impervious nature of the urban environment is represented as a weighted fraction of the
total tile. Currently this is only implemented for the saturation point of the grid cell. An increased urban fraction
results in a low saturation point of the top soil level, resulting in a high run-off rate. Various adaptations to the
tile average hydraulic properties are made by adjusting the Van Genuchten variables \citet{vg80}, for example the alpha- and lambda-
factors. Further developments are required to detail urban hydrology, including drainage, however these first steps
offer a parameterised solution to those more complex processes. Currently the urban tile affects run-off and infiltration in the top surface layer and the remaining sub-surface models are kept as soil.

\subsubsection{Urban Snow Interaction}

The urban properties relevant to snow accumulation are incorporated as part of the high vegetation snow, as a result the
treatment of snow is kept the same as in the non-urban version with only updated underlining thermodynamic and
hydrological properties to account for the urban environment. The albedo over the snow-covered urban environment is
expected to be lower than the surrounding bare soil for several reasons including increased snowmelt due to
anthropogenic heating, surface orography, anthropogenic snow clearing activities and darkening caused by aerosol
deposition. To account for the reduced urban snow albedo, the albedo value is calculated as the average between the urban and snow albedos, using the snow albedo description of \citet{Dutraetal2010}. The revised albedo of snow under high-vegetation albedo is given by equation \ref{equrbsnow}, which is simplistic and will require further consideration in future versions of the scheme.
\begin{equation}
 \alpha _{\mathit{snow}_{\mathit{new}}}=\alpha _{\mathit{snow}}\ast \frac{\left(C_H\right)}{(C_H+C_U)}+\frac 1 2(\alpha _{\mathit{snow}}+\alpha _{\mathit{urban}})\ast \frac{(C_U)}{(C_H+C_U)}
\label{equrbsnow}
\end{equation}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 \newpage

% References, variant B: external bibliography
%=====================================
\reftitle{References}
\externalbibliography{yes}
\bibliography{EcLand_ref.bib}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Citations and References in Supplementary files are permitted provided that they also appear in the reference list here. 

%=====================================
% References, variant A: internal bibliography
%=====================================
%\reftitle{References}
%\begin{thebibliography}{999}
% Reference 1
%\bibitem[Author1(year)]{ref-journal}
%Author1, T. The title of the cited article. {\em Journal Abbreviation} {\bf 2008}, {\em 10}, 142--149.
% Reference 2
%\bibitem[Author2(year)]{ref-book}
%Author2, L. The title of the cited contribution. In {\em The Book Title}; Editor1, F., Editor2, A., Eds.; Publishing House: City, Country, 2007; pp. 32--58.
%\end{thebibliography}

% The following MDPI journals use author-date citation: Arts, Econometrics, Economies, Genealogy, Humanities, IJFS, JRFM, Laws, Religions, Risks, Social Sciences. For those journals, please follow the formatting guidelines on http://www.mdpi.com/authors/references
% To cite two works by the same author: \citeauthor{ref-journal-1a} (\citeyear{ref-journal-1a}, \citeyear{ref-journal-1b}). This produces: Whittaker (1967, 1975)
% To cite two works by the same author with specific pages: \citeauthor{ref-journal-3a} (\citeyear{ref-journal-3a}, p. 328; \citeyear{ref-journal-3b}, p.475). This produces: Wong (1999, p. 328; 2000, p. 475)

%=====================================
% References, variant B: external bibliography
%=====================================
%\reftitle{References}
%\externalbibliography{yes}
%\bibliography{EcLand_ref.bib}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% optional
%\sampleavailability{Samples of the compounds ...... are available from the authors.}

%% for journal Sci
%\reviewreports{\\
%Reviewer 1 comments and authors’ response\\
%Reviewer 2 comments and authors’ response\\
%Reviewer 3 comments and authors’ response
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}

