Version 16 (modified by Giersch, 9 years ago) (diff)

--

Boundary conditions

PALM offers a variety of boundary conditions. Dirichlet or Neumann boundary conditions can be chosen for u, v, θ, qv, and p at the bottom and top of the model. For the horizontal velocity components the choice of Neumann (Dirichlet) boundary conditions yields free-slip (no-slip) conditions. Neumann boundary conditions are also used for the SGS-TKE. Kinematic fluxes of heat and moisture can be prescribed at the surface instead (Neumann conditions) of temperature and humidity (Dirichlet conditions). At the top of the model, Dirichlet boundary conditions can be used with given values of the geostrophic wind. By default, the lowest grid level (k = 0) for the scalar quantities and horizontal velocity components is not staggered vertically and defined at the surface (z = 0). In case of free-slip boundary conditions at the bottom of the model, the lowest grid level is defined below the surface (z = - 0.5 Δz) instead. Vertical velocity is assumed to be zero at the surface and top boundaries, which implies using Neumann conditions for pressure.

Constant flux layer

Basics

Following Monin-Obukhov similarity theory (MOST) a constant flux layer can be assumed as boundary condition between the surface and the first grid level where scalars and horizontal velocities are defined (k = 1, zMO = 0.5 Δz). It is then required to provide the roughness lengths for momentum z0 and heat z0,h. Momentum and heat fluxes as well as the horizontal velocity components are calculated using the following framework. The formulation is theoretically only valid for horizontally-averaged quantities. In PALM we assume that MOST can be also applied locally and we therefore calculate local fluxes, velocities, and scaling parameters.

Following MOST, the vertical profile of the horizontal wind velocity

\begin{equation*}
u_\mathrm{h} = (u^2 + v^2)^{\frac{1}{2}} 
\end{equation*}

is given in the surface layer by

\begin{align}
  & \frac{\partial u_\mathrm{h}}{\partial z} =
  \frac{u_\ast}{\kappa z}\Phi_\mathrm{m}\left(\frac{z}{L}\right)\;,
\end{align}

where κ = 0.4 is the Von Kármán constant and Φm is the similarity function for momentum in the formulation of Businger-Dyer (see e.g. Panofsky and Dutton 1984)

\begin{align}
  & \Phi_\mathrm{m} =
  \begin{cases}
    1 + 5 \frac{z}{L} & \text{for~}  \frac{z}{L} \geq 0 \\
    \left(1 - 16 \frac{z}{L}\right)^{-\frac{1}{4}} & \text{for~}
    \frac{z}{L} < 0\;.
  \end{cases}
\end{align}

Here, L is the Obukhov length, calculated as

\begin{align}
  & 
L = \frac{\theta_\mathrm{v}(z) u_\ast^2}{\kappa g
    \left[\theta_\ast + 0.61 \theta(z) q_\ast + 0.61
      q_\mathrm{v}(z) \theta_\ast\right]}\;.
\end{align}

The scaling parameters θ* and q* are defined by MOST as:

\begin{align}
  & \theta_\ast = -
  \frac{\overline{w^{\prime\prime}\theta^{\prime\prime}}_0}{u_\ast},~q_\ast
  = -
  \frac{\overline{w^{\prime\prime}q_\mathrm{v}^{\prime\prime}}_0}{u_\ast}\;,
\end{align}

with the friction velocity u* defined as

\begin{align}
  & u_\ast =
  \left[\left(\overline{u^{\prime\prime} w^{\prime\prime}}_0\right)^2
    + \left(\overline{v^{\prime\prime} w^{\prime\prime}}_0\right)^2
  \right]^{\frac{1}{4}}\;.
\end{align}

In PALM, u* is calculated from uh at zMO by vertical integration over z from z0 to zMO.

From the equations above it is possible to derive a formulation for the horizontal wind components, viz.

\begin{align}
  & 
\frac{\partial u}{\partial z} =
  \frac{-\overline{u^{\prime\prime} w^{\prime\prime}}_0}{u_\ast \kappa
    z}
  \Phi_\mathrm{m}\left(\frac{z}{L}\right)\,\text{and~}\,\frac{\partial
    v}{\partial z} = \frac{-\overline{v^{\prime\prime}
      w^{\prime\prime}}_0}{u_\ast \kappa z}
  \Phi_\mathrm{m}\left(\frac{z}{L}\right)\;.
\end{align}

Vertical integration of the above equation over z from z0 to zMO then yields the surface momentum fluxes

\begin{equation*}
\overline{u^{\prime\prime} w^{\prime\prime}}_0,\;\; \overline{v^{\prime\prime} w^{\prime\prime}}_0
\end{equation*}

The formulations above all require knowledge of the scaling parameters θ* and q*. These are deduced from vertical integration of

\begin{align}
  &
\frac{\partial \theta}{\partial z} =
  \frac{\theta_\ast}{\kappa z}
  \Phi_\mathrm{h}\left(\frac{z}{L}\right)~{\text{and}}~\frac{\partial
    q_\mathrm{v}}{\partial z} = \frac{q_\ast}{\kappa z}
  \Phi_\mathrm{h}\left(\frac{z}{L}\right)
\end{align}

over z from z0,h to zMO. The similarity function Φh is given by

\begin{align}
  & \Phi_\mathrm{h} =
  \begin{cases}
    1 + 5 \frac{z}{L} & \text{for~}  \frac{z}{L} \geq 0 \\
    \left(1 - 16 \frac{z}{L}\right)^{-1/2} & \text{for~} \frac{z}{L}
    < 0\;.
  \end{cases}
\end{align}

Implementation

Currently, there are three different options to calculate the Obukhov length and the surface fluxes which are steered via the NAMELIST parameter most_method.

Method 1: circular

The traditional implementation in PALM (most_method = 'circular') requires the use of data from the previous time step. The following steps are thus carried out in sequential order. First of all, θ* and q* are calculated by integration using the value of zMO/L from the previous time step. Second, the new value of zMO/L is derived using the new values of θ* and q* but using u* from the previous time step. Then, the new values of u*, and subsequently the momentum fluxes are calculated by integration, respectively. At last, the new surface fluxes are derived from θ* and q*, and u*. In the special case, when surface fluxes are prescribed instead of surface temperature and humidity, the first and last steps are omitted and θ* and q* are directly calculated from u* and the surface fluxes.

In summary, the following actions are performed in sequential order:

  1. calculate θ* and q*
  2. calculate uh
  3. determine Obukhov length
  4. calculate u*
  5. derive surface fluxes

Method 2: Newton iteration / lookup table

Alternatively, the Obukhov length can be calculated by solving an implicit equation relating the L to the bulk Richardson number. This can be achieved either by a Newton iteration algorithm (most_method = 'newton') or by using a lookup table (most_method = 'lookup'). Note that the latter is the new default in PALM as it is much faster than the Newton iteration method and the results are more precise compared to the circular method. However, it can only be used when the roughness lengths are homogeneously set on each processor.

Both methods require a different sequential order to derive the surface fluxes:

  1. calculate uh
  2. determine Obukhov length (Newton iteration or lookup table)
  3. calculate u*
  4. calculate θ* and q*
  5. derive surface fluxes

Depending on whether Neumann (prescribed fluxes) or Dirichlet boundary conditions are used for temperature and humidity, the bulk Richardson number is related to the Obukhov length via

\begin{equation*}
Ri_\mathrm{b,Di} = \dfrac{z}{L} \cdot \dfrac{[\phi_\mathrm{H}]}{[\phi_\mathrm{M}]^2} \;\;\;\;\textnormal{(Dirichlet conditions)}
\end{equation*}
\begin{equation*}
Ri_\mathrm{b,Ne} = \dfrac{z}{L} \cdot \dfrac{1}{[\phi_\mathrm{M}]^3} \;\;\;\;\textnormal{(Neumann conditions)}
\end{equation*}

where

\begin{equation*}
[\phi_\mathrm{H}] = \log\left(\dfrac{z_\mathrm{MO}}{z_\mathrm{0,h}}\right) - \Phi_\mathrm{H}\left(\dfrac{z_\mathrm{MO}}{L}\right) + \Phi_\mathrm{H}\left(\dfrac{z_\mathrm{0,h}}{L}\right)\;,
\end{equation*}
\begin{equation*}
[\phi_\mathrm{M}] = \log\left(\dfrac{z_\mathrm{MO}}{z_\mathrm{0}}\right) - \Phi_\mathrm{M}\left(\dfrac{z_\mathrm{MO}}{L}\right) + \Phi_\mathrm{M}\left(\dfrac{z_\mathrm{0,h}}{L}\right)\;,
\end{equation*}

are the integrated universal profile stability functions of Φm and Φh (see Paulson 1970, Holtslag and Bruin 1988), and

\begin{equation*}
Ri_\mathrm{b,Di} = \dfrac{g z_\mathrm{MO} \left(\theta_\mathrm{v,1} - \theta_\mathrm{v,0}\right)}{u_\ast^2 \theta_v}\;,
\end{equation*}
\begin{equation*}
Ri_\mathrm{b,Ne} = - \dfrac{g z_\mathrm{MO} \overline{v^{\prime\prime} \theta_\mathrm{v}^{\prime\prime}}_0}{\kappa^2 u_\mathrm{h}^3 \theta_v}\;.
\end{equation*}

In case of most_method = 'newton', the above equations are solved for L by Newton iteration, i.e. finding the root of the equation

\begin{equation*}
f = Ri_\mathrm{b} - \dfrac{z}{L} \cdot \dfrac{[\phi_\mathrm{H}]^x}{[\phi_\mathrm{M}]^y}
\end{equation*}

where x and y depend on the chosen boundary conditions (see above). The solution is given by iteration of

\begin{equation*}
L^{t+1} = L^t - \dfrac{f(L^t)}{f'(L^t)}
\end{equation*}

with

\begin{equation*}
f'(L) = \dfrac{df}{dL}
\end{equation*}

until L meets a convergence criterion.

If most_method = 'lookup' is used, a table of Rib against zMO/L is created at model start, based on the prescribed values of the roughness lengths. During the model run, Rib is calculated and the respectively value of zMO/L is retrieved from the lookup table and using linear interpolation between the discrete values in the table. In order to speed up this method, the value of Rib from the previous time step is used as initial value. Due to the fact that the lookup table is created at model start, it is essential that the roughness lengths are 1) homogeneous on each processor (limiting this method to homogeneous surface configurations) and should not vary during the simulation (should thus not be used when using dynamic roughness length over water surfaces).

For more details, see surface_layer_fluxes.f90.

References

  • Holtslag AAM, Bruin HARD. 1988. Applied modelling of the night-time surface energy balance over land. J. Appl. Meteorol., 27, 689–704.
  • Panofsky HA, Dutton JA. 1984. Atmospheric Turbulence, Models and Methods for Engineering Applications, John Wiley & Sons, New York.
  • Paulson CA 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J. Appl. Meteorol., 9, 857–861.

Attachments (6)

  • 03.png (465.4 KB) - added by Giersch 9 years ago. Turbulence recycling method
  • 04.png (91.3 KB) - added by Giersch 9 years ago. Sketch of the 2.5-D implementation of topography
  • 05.png (5.0 MB) - added by Giersch 9 years ago. Snapshot of turbulence structures induced by a densely built-up artificial island of the coast of Macau, China
  • grid_lr.png (89.4 KB) - added by Giersch 9 years ago. Grid structure at the lateral boundaries with non-cyclic lateral boundary conditions along the left-right direction
  • grid_ns.png (74.7 KB) - added by Giersch 9 years ago. Figure 2: Grid structure of the lateral boundaries with non-cyclic lateral boundary conditions along the north-south direction.
  • mask_method.png (19.6 KB) - added by suehring 8 years ago.