Home

Context Navigation

Land surface model

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 namelist parameters.

Overview

Since r1551 a full land surface model (LSM) is available in PALM (see 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 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 C₀ and T₀ are the heat capacity and radiative temperature of the surface skin layer, respectively. Note that C₀ 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

$\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 θ₀ and θ₁ 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:

$\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 surface_layer_fluxes_mod.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 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:

$\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₀ is then set to a non-zero value according to the material properties.

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, l_v = 2.5 * 10⁶ 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₀.

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

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

$\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}^N R_k\ max(m_\mathrm{soil,k}, m_\mathrm{wilt}) \end{equation*}$

where

$R_k$: Root fraction in layer $k$

and N is the number of soil layers.

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

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

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

$\begin{equation*} m_\mathrm{liq} / m_\mathrm{liq,max} \end{equation*}$

for vegetation (following the HTESSEL scheme, Viterbo & Beljaars, 1995) and

$\begin{equation*} (m_\mathrm{liq} / m_\mathrm{liq,max})^{0.67} \end{equation*}$

following 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.

Figure 12: Default layer distribution in the soil model.

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$: Heat conductivity of the soil layer$

The heat conductivity 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}$

Root extraction

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

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

$\begin{equation*} \sum \limits_{k = 1}^N R_\mathrm{fr}(k) = 1 . \end{equation*}$

The root extraction is then given by

$\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 et al. (1995) and 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₀, z_0,h, and z_0,q respectively) are hence calculated after Beljaars (1994) at each grid point as:

$\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 (pave_depth), a heat capacity (pave_heat_capacity), and a 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

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

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

$\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 vegetation_type) for which typical parameters are used to intialize the model. Moreover, a set of predefined soil types can be chosen (parameter 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 albedo_type.

A complete list of parameters and a detailed description can be found here.

References

Beljaars, ACM 1994. The parametrization of surface fluxes in large-scale models under free convection. Q. J. R. Met. Soc. 121: 255–270.

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, 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.

Masson, V. 2000. A physically-based scheme for the urban energy budget in atmospheric models, Boundary-layer Meteorol., 94, 357–397.

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.

Last modified 16 months ago Last modified on Aug 8, 2023 12:11:41 PM

Attachments (2)

2015_LSM.pdf (2.8 MB) - added by maronga 10 years ago. LSM introduction
wall_concept_type1.png (34.1 KB) - added by maronga 7 years ago. soil model (vegetation)

Download in other formats:

Plain Text