Version 26 (modified by maronga, 8 years ago) (diff)

--

Land surface model

Overview

Since r1551 a full land surface model (LSM) is available in PALM (see land_surface_model_mod.f90). It consists of a four layer soil model, predicting soil temperature and moisture content, and a solver for the energy balance, predicting the temperature of the skin layer. Moreover, a liquid water reservoir accounts for the presence of liquid water on plants and 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 et al. 2010).

Note that the use of the LSM requires using some kind of radiation model to provide radiative fluxes at the surface.

Energy balance solver

The energy balance of the Earth's surface reads

\begin{equation*}
   C_0 \dfrac{dT_0}{dt} = R_\mathrm{n} - H - LE - G
\end{equation*}

where C0 and T0 are the heat capacity and radiative temperature of the surface skin layer, respectively. Rn, 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

\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, cp = 1005 J kg-1 K-1 is the specific heat at constant pressure, ra 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. ra 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:

\begin{equation*}
r_\mathrm{a} = \dfrac{u_*\ \theta_*}{\theta_1 - \theta_0}
\end{equation*}

where u* and θ* are the friction velocity and the characteristic temperature scale according to Monin-Obukhov similarity scaling (these are calculated in surface_layer_fluxes.f90).

Parameterization of G

G is parametrized as (see Duynkerke 1999)

\begin{equation*}
  G = \Lambda ( T_0 - T_{\mathrm{soil},1} )
\end{equation*}

with Λ being the heat conductivity between skin layer and the soil, and Tsoil,1 being the temperature of the uppermost soil layer.

Parameterization of LE

The latent heat flux is calculated as

\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, lv = 2.5 * 106 J kg-1 is the latent heat of vaporisation, rs is the surface resistance, qv,1 is the specific humidity at first grid level, and qv,sat is the saturation specific humidity at temperature T0.

All equations above are solved locally for each surface element of the LES grid. Each element can consist of both 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. 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

\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

$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

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

\begin{equation*}
\dfrac{1}{f_3(e_\mathrm{def})} = exp(g_\mathrm{D}\ e_\mathrm{def})
\end{equation*}

where gD is a correction factor (for high vegetation only, otherwise zero). Moreover, the reaction of plants to water availability in the soil is considered:

\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

$m_\mathrm{wilt}$: Permanent wilting point\\
$m_\mathrm{fc}$: Field capacity

and

\begin{equation*}
 \tilde m = \sum\limits_{k = 1}^4 R_k\ max(m_\mathrm{soil,k}, m_\mathrm{wilt})
\end{equation*}

where

$R_k$: Root fraction in layer $k$

The bare soil resistance is given by

\begin{equation*}
r_\mathrm{soil} = r_\mathrm{soil,min}\ f_{2b}(m_\mathrm{soil,1})
\end{equation*}

with

$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

\begin{equation*}
LE = c_\mathrm{veg} (1 - c_\mathrm{liq})\ LE_\mathrm{veg} + c_\mathrm{liq}\ c_\mathrm{veg}\ LE_\mathrm{liq} + (1 - c_\mathrm{veg}) \ LE_\mathrm{soil}
\end{equation*}

where cveg, and cliq 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 mliq reads

\begin{equation*}
\dfrac{d m_\mathrm{liq}}{dt} = \dfrac{LE_\mathrm{liq}}{\rho_\mathrm{l}\ l_\mathrm{v}}
\end{equation*}

For positive values of LEliq, liquid water is evaporating from the surface, while negative values indicate precipitation (rain, dew).

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 four layers (see figure 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. No image "lsm_simple.png" attached to doc/tec/lsm

Soil heat transport

The prognostic equation for soil temperature reads

\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

$(\rho C)_\mathrm{soil}$: Volumetric heat capacity of the soil layer\\
$\lambda_T$: Thermal heat conductivity of the soil layer

The volumetric heat capacity is calculated as

\begin{equation*}
   \lambda_T = Ke\ (\lambda_{T,\mathrm{sat}} - \lambda_{T,\mathrm{dry}}) + \lambda_{T,\mathrm{dry}}
\end{equation*}

with

$\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

\begin{equation*}
\lambda_{T,\mathrm{sat}} = \lambda_{T,\mathrm{sm}}^{1-m_\mathrm{soil,sat}}\ \lambda_\mathrm{m}
\end{equation*}

with

$\lambda_{T,\mathrm{sm}}$: Heat conductivity of the soil matrix\\
$\lambda_\mathrm{m}$: Heat conductivity of water\\

The Kersten number is calculated as

\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

\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

$\lambda_m$: Hydraulic diffusion coefficient\\
$\gamma$: Hydraulic conductivity\\
$S_m$: Sink term due to root extraction

The hydraulic diffusion coefficient is calculated after Clapp and Hornberger (1978) as

\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

$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 Van Genuchten (1980) (as in H-TESSEL):

\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

$alpha$: Van Genuchten coefficient\\
$h$: Pressure head\\
$n$: Van Genuchten coefficient\\
$l$: Van Genuchten coefficient

and

\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}

or after Clapp and Hornberger (1978):

\begin{equation}
\gamma = \gamma_\mathrm{sat} \left(\dfrac{m_\mathrm{soil}}{m_\mathrm{sat}}\right)^{2b + 3}
\end{equation}

Root extraction

The root extraction of water from the respective soil layer Sm is calculated as follows:

\begin{equation*}
m_\mathrm{total} = \sum \limits_{k = 1}^4 R_\mathrm{fr}(k)\ m_\mathrm{soil}(k)
\end{equation*}

where mtotal is the total water content of the soil and Rfr 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}4 R_\mathrm{fr}(k) = 1 \end{equation*} }}} The root extraction is then given by

\begin{equation*}
S_m = \dfrac{LE_\mathrm{veg}}{\rho_l\ l_\mathrm{v}} \sum \limits_{k = 0}^3 R_\mathrm{fr}(k)\ \dfrac{m_\mathrm{soil}(k)}{m_\mathrm{total}(k)} 
\end{equation*}

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 LEsoil, 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 et al. (1995) and Balsamo et al. (2009).

Technical details

The discretized and linearized energy budget equation in PALM reads

\begin{equation*}
T_\mathrm{0,p} = \dfrac{A \Delta t + C_\mathrm{sk} T_0}{C_\mathrm{sk} + B \Delta t} 
\end{equation*}

with

\begin{equation*}
A = R_\mathrm{n} + 3 \sigma T_0^4 + \dfrac{f_H}{\Pi} T_1 + f_{LE} \left( q_1 - q_s + \dfrac{d q_s}{d T} T_0 \right) + \Lambda T_\mathrm{soil}
\end{equation*}

and

\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):

$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$: specific humidity at first grid level\\
$q_s$: saturation specific humidity 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 Csk = 0, the prognostic equation for T0,p reduces to a diagnostic equation:

\begin{equation*}
T_\mathrm{0,p} = \dfrac{A}{B} 
\end{equation*}

Usage

References

  • 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, RB and Hornberger GM. 1978. Empirical Equations for Some Soil Hydraulic Properties. Water Res. Res., 14: 601-604.

  • Duynkerke PG. 1999. Turbulence, radiation and fog in Dutch stable boundary layers. Boundary-Layer Meteorol. 90: 447–477, doi:10.1023/A:1026441904734.
  • 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 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.
  • 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 P, Beljaars ACM. 1995. An Improved Land Surface Parameterization Scheme in the ECMWF Model and Its Validation. J. Climate 8: 2716–2748.

Attachments (2)