= Land surface model =

[[NoteBox(note,This page is part of the **Land Surface Model** (LSM) documentation. \\ It describes the physical and numerical realization of the LSM. \\ Please also see the **[wiki:doc/app/land_surface_parameters namelist parameters]**., 450px)]]

== Overview ==
Since r1551 a full land surface model (LSM) is available in PALM (see [source:palm/trunk/SOURCE/land_surface_model_mod.f90 land_surface_model_mod.f90]). It consists of a multi layer soil model, predicting soil temperature and moisture content, and a solver for the energy balance, predicting the temperature of the surface or the skin layer (depending on land use classification). Moreover, a liquid water reservoir accounts for the presence of liquid water on plants, pavement, and bare soil due to precipitation. The implementation is based on the ECMWF-IFS land surface parametrization (H-TESSEL) and its adaptation in the DALES model ([#heus Heus et al. 2010]).

Note that the use of the LSM requires using some kind of [wiki:doc/tec/radiation radiation model] to provide radiative fluxes at the surface.

== Energy balance solver ==
The energy balance of the Earth's surface reads
{{{
#!Latex
\begin{equation*}
   C_0 \dfrac{dT_0}{dt} = R_\mathrm{n} - H - LE - G
\end{equation*}
}}}
where ''C'',,0,, and ''T'',,0,, are the heat capacity and radiative temperature of the surface skin layer, respectively. Note that  ''C'',,0,, is usually zero as it is assumed that the skin layer has no vertical extent and thus does not have a heat capacity (see also below). ''R'',,n,,, ''H'', ''LE'', and ''G'' are the net radiation, sensible heat flux, latent heat flux, and ground (soil) heat flux at the surface, respectively. 

=== Parameterization of ''H'' ===
''H'' is calculated as
{{{
#!Latex
\begin{equation*}
  H = - \rho\ c_\mathrm{p}\ \dfrac{1}{r_\mathrm{a}} ( \theta_1 - \theta_0 )
\end{equation*}
}}}
where ''ρ'' is the density of the air, ''c'',,p,, = 1005 J kg^-1^ K^-1^ is the specific heat at constant pressure, ''r'',,a,, is the aerodynamic resistance, and ''θ'',,0,, and ''θ'',,1,, are the potential temperature at the surface and at the first grid level above the surface, respectively. ''r'',,a,, is calculated via Monin-Obukhov similarity theory, based on roughness lengths for heat and momentum and the assumption of a constant flux layer between the surface and the first grid level:
{{{
#!Latex
\begin{equation*}
r_\mathrm{a} = \dfrac{\theta_1 - \theta_0}{u_*\ \theta_*}
\end{equation*}
}}}
where ''u'',,*,, and ''θ'',,*,, are the friction velocity and the characteristic temperature scale according to Monin-Obukhov similarity scaling (these are calculated in [source:palm/trunk/SOURCE/surface_layer_fluxes_mod.f90 surface_layer_fluxes_mod.f90]).

=== Parameterization of ''G'' ===
''G'' is parametrized as (see [#dynkerke Duynkerke 1999])
{{{
#!Latex
\begin{equation*}
  G = \Lambda ( T_0 - T_{\mathrm{soil},1} )
\end{equation*}
}}}
with ''Λ'' being the total thermal conductivity between the skin and the center of the uppermost soil layer, and ''T'',,soil,1,, being the temperature of the uppermost soil layer. ''Λ''  is calculated as a combination of the conductivity between the canopy and the soil-top (constant value) and the conductivity of the top half of the uppermost soil layer:
{{{
#!Latex
\begin{equation*}
\Lambda = \dfrac{\Lambda_\mathrm{skin}\Lambda_\mathrm{soil}}{\Lambda_\mathrm{skin}+\Lambda_\mathrm{soil}}
\end{equation*}
}}}
When no skin layer is used (i.e. in case of bare soil and pavements), ''Λ'' reduces to the heat conductivity of the uppermost soil layer (divided by the layer depth). In that case, it is assumed that the soil temperature is constant within the uppermost 1/4 of the top soil layer and equals the radiative temperature. ''C'',,0,,  is then set to a non-zero value according to the material properties.

=== Parameterization of ''LE'' ===
The latent heat flux is calculated as 
{{{
#!Latex
\begin{equation*}
  LE = - \rho\ l_\mathrm{v}\ \dfrac{1}{r_\mathrm{a} + r_\mathrm{s}} ( q_{\mathrm{v},1} - q_{\mathrm{v,sat}}(T_0) )\;.
\end{equation*}
}}}
Here, ''l'',,v,, = 2.5 * 10^6^ J kg^-1^ is the latent heat of vaporisation, ''r'',,s,, is the surface resistance, ''q'',,v,1,, is the water vapor mixing ratio at first grid level, and ''q'',,v,sat,, is the water vapor mixing ratio at temperature ''T'',,0,,. 

All equations above are solved locally for each surface element of the LES grid. Each element for the surface type 'vegetation' can consist of patches of bare soil, vegetation, and a liquid water reservoir, which is the interception water stored on plants and soil from precipitation. Therefore, an additional equation is solved for the liquid water reservoir. A liquid water reservoir is also available when the surface type is set to 'pavement'. ''LE'' is then calculated for each of the three components (bare soil, vegetation, liquid water). The resistances are calculated separately for bare soil and vegetation following Jarvis (1976). The canopy resistance is calculated as
{{{
#!Latex
\begin{equation*}
r_\mathrm{c} = \dfrac{r_\mathrm{c,min}}{LAI}\ f_1(R_\mathrm{sw,in})\ f_2(\widetilde m)\ f_3(e_\mathrm{def})
\end{equation*}
}}}
with
{{{
#!Latex
$r_\mathrm{c,min}$: Minimum stomatal resistance\\
$LAI$: Leaf area index\\
$f_i$: Correction functions ($f_i \geq 1$)\\
$R_\mathrm{sw,in}$: Incoming shortwave radiation\\
$\widetilde m$: Layer-averaged soil moisture\\
$e_\mathrm{def}$: Water-vapor pressure deficit
}}}
The correction functions read
{{{
#!Latex
\begin{equation*}
\dfrac{1}{f_1(R_\mathrm{sw,in})} = min\left(1, \dfrac{0.004\ R_\mathrm{sw,in}}{0.81(0.004\ R_\mathrm{sw,in} + 1)}\right)
\end{equation*}
}}}
which accounts for the reaction of plants to shortwave radiation (opening/closing stomata),
{{{
#!Latex
\begin{equation*}
\dfrac{1}{f_3(e_\mathrm{def})} = exp(g_\mathrm{D}\ e_\mathrm{def})
\end{equation*}
}}}
where ''g'',,D,, is a correction factor (for high vegetation only, otherwise zero). Moreover, the reaction of plants to water availability in the soil is considered:
{{{
#!Latex
\begin{equation*}
\dfrac{1}{f_2(\tilde m)} = 
\begin{cases}
0 & \tilde m < m_\mathrm{wilt}\\
\dfrac{\tilde m - m_\mathrm{wilt}}{m_\mathrm{fc} - m_\mathrm{wilt}} & m_\mathrm{wilt} \leq \tilde m \leq m_\mathrm{fc}\\
1 & \tilde m > m_\mathrm{fc}
\end{cases}
\end{equation*}
}}}
where
{{{
#!Latex
$m_\mathrm{wilt}$: Permanent wilting point\\
$m_\mathrm{fc}$: Field capacity
}}}
and
{{{
#!Latex
\begin{equation*}
 \tilde m = \sum\limits_{k = 1}^N R_k\ max(m_\mathrm{soil,k}, m_\mathrm{wilt})
\end{equation*}
}}}
where
{{{
#!Latex
$R_k$: Root fraction in layer $k$
}}}
and N is the number of soil layers.

The bare soil resistance is given by
{{{
#!Latex
\begin{equation*}
r_\mathrm{soil} = r_\mathrm{soil,min}\ f_{2b}(m_\mathrm{soil,1})
\end{equation*}
}}}
with
{{{
#!Latex
$r_\mathrm{soil,min}$: Minimum soil resistance\\
$f_{2b}$: Correction function ($f_i \geq 1$)\\
$m_\mathrm{soil,1}$: Soil moisture of the uppermost layer)
}}}
The total evapotranspiration is then calculated as
{{{
#!Latex
\begin{equation*}
LE = c_\mathrm{veg} (1 - c_\mathrm{liq})\ LE_\mathrm{veg} + c_\mathrm{liq}\ LE_\mathrm{liq} + (1 - c_\mathrm{liq}) (1 - c_\mathrm{veg}) \ LE_\mathrm{soil}
\end{equation*}
}}}
where ''c'',,veg,,, and ''c'',,liq,, is the surface fraction covered with vegetation and liquid water, respectively.


=== Prognostic equation for the liquid water reservoir ===
The prognostic equation for the liquid water stored on plants and bare soil ''m'',,liq,, reads
{{{
#!Latex
\begin{equation*}
\dfrac{d m_\mathrm{liq}}{dt} = \dfrac{LE_\mathrm{liq}}{\rho_\mathrm{l}\ l_\mathrm{v}}
\end{equation*}
}}}
The maximum amount of water that can be stored on plants is calculated via
{{{
#!Latex
\begin{equation*}
m_\mathrm{liq,max} = \min\left(1, m_\mathrm{liq,ref}\ c_\mathrm{veg}\ LAI + (1 - c_\mathrm{veg})\right)\;,
\end{equation*}
}}}
where ''m'',,,liq,ref,,, = 0.2 mm is the reference liquid water column on a single leaf or bare soil. Exceeding liquid water is directly removed from the surface and infiltrated in the underlying soil. For positive values of ''LE'',,liq,,, liquid water is evaporating from the surface, while negative values indicate precipitation (rain, dew). The values of ''c,,,liq,,, are calculated either as the ratio
{{{
#!Latex
\begin{equation*}
m_\mathrm{liq} / m_\mathrm{liq,max}
\end{equation*}
}}}
for vegetation (following the HTESSEL scheme, [#viterbo Viterbo & Beljaars, 1995]) and 
{{{
#!Latex
\begin{equation*}
(m_\mathrm{liq} / m_\mathrm{liq,max})^{0.67}
\end{equation*}
}}}
following [#masson Masson (2000)].

== Soil model ==
The soil model consists of prognostic equations for the soil temperature and the volumetric soil moisture which are solved for multiple layers. The soil model only takes into account vertical transport within the soil and no ice phase is considered. By default, the soil model consists of eight layers (see Fig. 12 below), in which the vertical heat and water transport is modelled using the Fourier law of diffusion and Richards' equation, respectively. Also, root fractions can be assigned to each soil layer to account for the explicit water withdrawal of plants used for transpiration from the respective soil layer. 

[[Image(wall_concept_type1.png,600px, border=1)]]

Figure 12: Default layer distribution in the soil model. 


=== Soil heat transport ===
The prognostic equation for soil temperature reads
{{{
#!Latex
\begin{equation*}
(\rho C)_\mathrm{soil} \dfrac{\partial T_\mathrm{soil}}{\partial t} = \dfrac{\partial}{\partial z} \left( \lambda_T \dfrac{\partial T_\mathrm{soil}}{\partial z}\right)
\end{equation*}
}}}
with
{{{
#!Latex
$(\rho C)_\mathrm{soil}$: Volumetric heat capacity of the soil layer\\
$\lambda_T$: Heat conductivity of the soil layer
}}}
The heat conductivity is calculated as
{{{
#!Latex
\begin{equation*}
   \lambda_T = Ke\ (\lambda_{T,\mathrm{sat}} - \lambda_{T,\mathrm{dry}}) + \lambda_{T,\mathrm{dry}}
\end{equation*}
}}}
with
{{{
#!Latex
$\lambda_{T,\mathrm{sat}}$: Heat conductivity of saturated soil\\
$\lambda_{T,\mathrm{dry}}$: Heat conductivity of dry soil
$Ke$: Kersten number
}}}
The heat conductivity of saturated soil is given by
{{{
#!Latex
\begin{equation*}
\lambda_{T,\mathrm{sat}} = \lambda_{T,\mathrm{sm}}^{1-m_\mathrm{soil,sat}}\ \lambda_\mathrm{m}
\end{equation*}
}}}
with
{{{
#!Latex
$\lambda_{T,\mathrm{sm}}$: Heat conductivity of the soil matrix\\
$\lambda_\mathrm{m}$: Heat conductivity of water\\
}}}
The Kersten number is calculated as
{{{
#!Latex
\begin{equation*}
Ke = \log_{10} \left[max\left(0.1, \dfrac{m_\mathrm{soil}}{m_\mathrm{sat}}\right) \right] + 1
\end{equation*}
}}}


=== Soil moisture transport ===
The prognostic equation for volumetric soil moisture reads
{{{
#!Latex
\begin{equation*}
\dfrac{\partial m_\mathrm{soil}}{\partial t} = \dfrac{\partial}{\partial z} \left( \lambda_m \dfrac{\partial m_\mathrm{soil}}{\partial z} - \gamma\right) + S_m
\end{equation*}
}}}
with
{{{
#!Latex
$\lambda_m$: Hydraulic diffusion coefficient\\
$\gamma$: Hydraulic conductivity\\
$S_m$: Sink term due to root extraction 
}}}

The hydraulic diffusion coefficient is calculated after [#clapp Clapp and Hornberger (1978)] as
{{{
#!Latex
\begin{equation*}
\lambda_m = \dfrac{b \gamma_\mathrm{sat}(-\Psi_\mathrm{sat})}{m_\mathrm{sat}}\left(\dfrac{m_\mathrm{soil}}{m_\mathrm{sat}}\right)^{b+2}
\end{equation*}
}}}
with
{{{
#!Latex
$b = 6.04$: exponent\\
$\gamma_\mathrm{sat}$: Hydraulic conductivity at saturation\\
$\Psi_\mathrm{sat} = -338 m$: Matrix potential at saturation
}}}

The hydraulic conductivity is calculated either after [#vangenuchten Van Genuchten (1980)] (as in H-TESSEL):
{{{
#!Latex
\begin{equation*}
 \gamma = \gamma_\mathrm{sat} \dfrac{\left[(1 + (\alpha h)^n)^{1 - 1/n} - (\alpha h)^{n-1}\right]^2}{(1+ (\alpha h)^n)^{(1 - 1/n)(l + 2)}}
\end{equation*}
}}}
where
{{{
#!Latex
$\alpha$: Van Genuchten coefficient\\
$h$: Pressure head\\
$n$: Van Genuchten coefficient\\
$l$: Van Genuchten coefficient
}}}
and
{{{
#!Latex
\begin{equation}
m_\mathrm{soil}(h) = m_\mathrm{res} + \dfrac{m_\mathrm{sat} - m_\mathrm{res}}{(1 + (\alpha h)^n)^{1 - 1/n}}\;.
\end{equation}
}}}
=== Root extraction ===
The root extraction of water from the respective soil layer ''S,,m,,'' is calculated as follows:
{{{
#!Latex
\begin{equation*}
m_\mathrm{total} = \sum \limits_{k = 1}^N R_\mathrm{fr}(k)\ m_\mathrm{soil}(k)
\end{equation*}
}}}
where ''m'',,total,, is the total water content of the soil and ''R'',,fr,, is the root fraction in soil layer ''k''. Only those layer are summed up which have a soil moisture above wilting point (plants do not extract water from such layers). The root profiles must prescribed in such a way that
{{{
#!Latex
\begin{equation*}
\sum \limits_{k = 1}^N R_\mathrm{fr}(k) = 1 .
\end{equation*}
}}}
The root extraction is then given by
{{{
#!Latex
\begin{equation*}
S_{m,k} = \dfrac{LE_\mathrm{veg}}{\rho_l\ l_\mathrm{v}} \dfrac{R_\mathrm{fr}(k)}{dzs(k)}\ \dfrac{m_\mathrm{soil}(k)}{m_\mathrm{total}} .
\end{equation*}
}}}
Here, ''dzs(k)'' is the thickness of layer ''k''. Again, only layers where the soil moisture is above its value at wilting point are used for root extraction (which is zero otherwise).


=== Boundary conditions ===
Neumann boundary conditions are used for the transport of heat and moisture at the upper boundary (surface). The values are given by the energy balance in terms of ''G'' for heat and ''LE'',,soil,,, for moisture. At the bottom boundary a deep soil temperature is prescribed (Dirichlet conditions), whereas two options are available for soil moisture. The underlying surface can be set to either bedrock (no moisture flux at the bottom, water conserving) or to open bottom (implies non-conservation of water). 


For more details, see also [#viterbo Viterbo et al. (1995)] and [#balsamo Balsamo et al. (2009)].

== Exception: water surfaces ==
When prescribing grid points with water surface (i.e. {{{surface_type = 'water'}}}), the energy balance is solved as for land surface without evaporation from vegetation canopy and bare soil. Moreover, the water temperature is constant in time and is derived from the value of {{{pt_surface}}} at model start. In order to obtain realistic results, it is required to set a zero heat capacity for the ocean skin layer (e.g. {{{c_surface = 0.0}}}).

A special feature of the treatment of water surfaces is the treatment of surface roughness chances due to water waves. The roughness lengths for momentum, heat, and moisture (''z'',,0,,, ''z'',,0,h,,, and ''z'',,0,q,, respectively) are hence calculated after [#beljaars Beljaars (1994)] at each grid point as:
{{{
#!Latex
\begin{align*}
z_0 &= \dfrac{0.11 \nu}{u_\ast} + \alpha_\mathrm{Ch} \dfrac{u_\ast^2}{g}\\
z_\mathrm{0,h} &= \dfrac{0.4 \nu}{u_\ast}\\
z_\mathrm{0,q} &= \dfrac{0.62 \nu}{u_\ast}
\end{align*}
}}}
Here, ''ν'' is the molecular viscosity, and ''α'',,Ch,, = 0.018 is the Charnock constant taken from the ECMWF-IFS model.

Please note that this parameterization was developed for large-scale models where waves are a pure subgrid-scale phenomenon. When using large-eddy simulations, however, this might no longer be the case. It has not been studied whether this parameterization is appropriate in such cases and should only be used with caution.

== Exception: pavement ==
It is possible to account for urban land surfaces such as roads by adding a pavement layer to the soil model. This is realized by setting the parameter {{{surface_type = 'pavement'}}}. The pavement is steered via a depth ([wiki:doc/app/land_surface_parameters#pave_depth pave_depth]), a heat capacity ([wiki:doc/app/land_surface_parameters#pave_heat_capacity pave_heat_capacity]), and a heat conductivity ([wiki:doc/app/land_surface_parameters#pave_heat_conductivity pave_heat_conductivity]). The pavement then replaces the upper soil layers up to a depth of {{{pave_depth}}}. In case that {{{pave_depth}}} is between two soil layers, the respective heat conductivity and heat capacities are linearly interpolated between the soil value and the pavement value to the respective grid point.

Note that the pavement layer must be at least have the same depth as the uppermost soil layer. If the prescribed pavement depth is too small, it is automatically set to this minimum depth.

The pavement is able to hold a maximum liquid water column of 1 mm from precipitation, which can also evaporate. The soil below the pavement is assumed to be completely dry.


== Technical details ==
The discretized and linearized energy budget equation in PALM reads
{{{
#!Latex
\begin{equation*}
T_\mathrm{0,p} = \dfrac{A \Delta t + C_\mathrm{sk} T_0}{C_\mathrm{sk} + B \Delta t} 
\end{equation*}
}}}
with
{{{
#!Latex
\begin{equation*}
A = R_\mathrm{n} + 3 \sigma T_0^4 + f_H \theta_1 + f_{LE} \left( q_1 - q_s + \dfrac{d q_s}{d T} T_0 \right) + \Lambda T_\mathrm{soil}
\end{equation*}
}}}
and
{{{
#!Latex
\begin{equation*}
B = \Lambda + f_{LE} \dfrac{d q_s}{d T} + \dfrac{f_H}{\Pi} + 4 \sigma T_0^3
\end{equation*}
}}}
with (in order of occurence):
{{{
#!Latex
$T_\mathrm{0,p}$:  prognostic skin temperature\\
$\Delta t$: time step\\
$C_\mathrm{sk}$: skin heat capacity\\
$T_\mathrm{0}$: skin temperature before time step\\
$R_\mathrm{n}$: net radiation at the surface\\
$\sigma$: Stefan-Boltzmann constant\\
$f_{H} = \dfrac{\rho c_\mathrm{p}}{r_\mathrm{a}}$\\
$\Pi$: conversion factor from actual temperature to potential temperature\\
$T_1$: temperature at first grid level\\
$f_{LE} = \dfrac{\rho l_\mathrm{v}}{r_\mathrm{a} + r_\mathrm{s}}$\\
$q_1$: water vapor mixing ratio at first grid level\\
$q_s$: saturation water vapor mixing ratio at the surface\\
$\Lambda$: heat conductivity of the skin layer\\
$T_\mathrm{soil}$: temperature of the uppermost soil layer
}}}
Time stepping is the same as in the atmospheric part of the model (default: 3rd-order Runge-Kutta). 

Note that for ''C'',,sk,, = 0, the prognostic equation for ''T'',,0,p,, reduces to a diagnostic equation:
{{{
#!Latex
\begin{equation*}
T_\mathrm{0,p} = \dfrac{A}{B} 
\end{equation*}
}}}


== Job preparation ==
The lsm is activated by specifying a {{{&land_surface_parameters}}} namelist in the _p3d file, e.g.:
{{{
 &land_surface_parameters
          surface_type = 'vegetation',
          vegetation_type = 2, 
          soil_type = 3,
          conserve_water_content = .T.,
          c_surface = 0.0,
          dz_soil = 0.01, 0.02, 0.04, 0.07, 0.15, 0.21, 0.72, 1.89,
          root_fraction = 0.1, 0.2, 0.3, 0.2, 0.1, 0.05, 0.05, 0.0,
          soil_temperature = 290.0, 289.0, 288.0, 286.0, 285.0, 285.0, 285.0, 285.0,
          deep_soil_temperature = 280.0,
/
}}}
In particular, PALM provides a set of predefined land surface types (parameter [wiki:doc/app/land_surface_parameters#vegetation_type vegetation_type]) for which typical parameters are used to intialize the model. Moreover, a set of predefined soil types can be chosen (parameter [wiki:doc/app/land_surface_parameters#soil_type soil_type]). You can overwrite the default parameters by explicitly prescribing them in the namelist.

Note that the surface albedo is part of the radiation scheme and is steered via the parameter [wiki:doc/app/radiation_parameters#albedo_type albedo_type].

A complete list of parameters and a detailed description can be found [wiki:doc/app/land_surface_parameters here].

== References ==
* [=#beljaars]'''Beljaars, ACM''' 1994. The parametrization of surface fluxes in large-scale models under free convection. Q. J. R. Met. Soc. 121: 255–270.

* [=#balsamo]'''Balsamo G, Vitebo P, Beljaars A, van den Hurk B, Hirschi M, Betts AK, Scipal K.''' 2009. A revised hydrology for the ECMWF model: Verification from field site to terrestrial water storage and impact in the integrated forecast system. J. Hydrometeorol. 10: 623–643.

* [=#clapp]'''Clapp RB, Hornberger GM.''' 1978. Empirical Equations for Some Soil Hydraulic Properties. Water Res. Res. 14: 601-604.
 
* [=#duynkerke]'''Duynkerke PG.''' 1999. Turbulence, radiation and fog in Dutch stable boundary layers. Boundary-Layer Meteorol. 90: 447–477, doi:10.1023/A:1026441904734.

* [=#heus]'''Heus T, Van Heerwaarden CC, Jonker HJJ, Siebesma AP, Axelsen S, Dries K, Geoffroy O, Moene AF, Pino D, De Roode SR, Vil`a-Guerau de Arellano J.''' 2010. Formulation of the dutch atmospheric large-eddy simulation (dales) and overview of its applications. Geosci. Model Dev. 3: 415–444.

* [=#jarvis]'''Jarvis PG.''' 1976. The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field. Philos. Trans. Roy. Soc. London 273B: 593–610.

* [=#masson]'''Masson, V.''' 2000. A physically-based scheme for the urban energy budget in atmospheric models, Boundary-layer Meteorol., 94, 357–397.

* [=#vangenuchten]'''van Genuchten M.''' 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Amer. J. 44: 892–898.

* [=#viterbo]'''Viterbo P, Beljaars ACM.''' 1995. An Improved Land Surface Parameterization Scheme in the ECMWF Model and Its Validation. J. Climate 8: 2716–2748.