[[NoteBox(warn,This site is currently under construction!)]] = Boundary conditions = == Basics == PALM offers a variety of boundary conditions. Dirichlet or Neumann boundary conditions can be chosen for ''u'', ''v'', ''θ'', ''q'',,v,,, 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. 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, ''z'',,MO,, = 0.5 ''Δz''). It is then required to provide the roughness lengths for momentum ''z'',,0,, and heat ''z'',,0,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 {{{ #!Latex \begin{equation*} u_\mathrm{h} = (u^2 + v^2)^{\frac{1}{2}} \end{equation*} }}} is given in the surface layer by {{{ #!Latex \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 Panofsky and Dutton 1984]) {{{ #!Latex \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 {{{ #!Latex \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: {{{ #!Latex \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 {{{ #!Latex \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 ''u'',,h,, at ''z'',,MO,, by vertical integration over ''z'' from ''z'',,0,, to ''z'',,MO,,. From the equations above it is possible to derive a formulation for the horizontal wind components, viz. {{{ #!Latex \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 over ''z'' from ''z'',,0,, to ''z'',,MO,, of the equation above then yields the surface momentum fluxes {{{ #!Latex $\overline{u^{\prime\prime} w^{\prime\prime}}_0,\;\; \overline{v^{\prime\prime} w^{\prime\prime}}_0$ }}} The formulations above all require knowledge of the scaling parameters ''θ'',,*,, and ''q'',,*,,. These are deduced from vertical integration of {{{ #!Latex \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 ''z'',,0,h,, to ''z'',,MO,,. The similarity function Φ,,h,, is given by {{{ #!Latex \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} }}} Currently, there are three different options to calculate the Obukhov length and the surface fluxes which are steered via the NAMELIST parameter [wiki:doc/app/inipar#most_method 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 of the corresponding vertical derivation functions mentioned above using the value of ''z'',,MO,,/L from the previous time step. Second, the new value of ''z'',,MO,,/L is derived from the equation for ''L'' 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 above equations for the scaling parameters are employed to calculate the new surface fluxes by using ''θ'',,*,, 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 ''u'',,h,, 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 ''u'',,h,, 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 {{{ #!Latex \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 {{{ #!Latex \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 Paulson 1970], [#holtslag Holtslag and Bruin 1988]), and {{{ #!Latex \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 {{{ #!Latex \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 {{{ #!Latex \begin{equation*} L^{t+1} = L^t - \dfrac{f(L^t)}{f'(L^t)} \end{equation*} }}} with {{{ #!Latex \begin{equation*} f'(L) = \dfrac{df}{dL} \end{equation*} }}} until ''L'' meets a convergence criterion. If {{{most_method = 'lookup'}}} is used, a table of ''Ri'',,b,, against ''z'',,MO,,/''L'' is created at model start, based on the prescribed values of the roughness lengths. During the model run, ''Ri,,b,,'' is calculated and the respectively value of ''z'',,MO,,/''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 ''Ri,,b,,'' 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 [source:palm/trunk/SOURCE/surface_layer_fluxes.f90 surface_layer_fluxes.f90]. Furthermore, the flat bottom of the model can be replaced by a Cartesian topography (see Sect. [wiki:doc/tec/bc#Topography Topography]). By default, lateral boundary conditions are set to be cyclic in both directions. Alternatively, it is possible to opt for non-cyclic conditions in one direction, i.e., a laminar or turbulent inflow boundary (see Sect. [wiki:doc/tec/bc#Laminarandturbulentinflowboundaryconditions Laminar and turbulent inflow boundary conditions]) and an open outflow boundary on the opposite site (see Sect. [wiki:doc/tec/bc#Openoutflowboundaryconditions Open outflow boundary conditions]). The boundary conditions for the other direction have to remain cyclic. In order to prevent gravity waves from being reflected at the top boundary, a sponge layer (Rayleigh damping) can be applied to all prognostic variables in the upper part of the model domain ([#klemp1978 Klemp and Lilly, 1978]). Such a sponge layer should be applied only within the free atmosphere, where no turbulence is present. The model is initialized by horizontally homogeneous vertical profiles of potential temperature, specific humidity (or a passive scalar), and the horizontal wind velocities. The latter can be also provided from a 1-D precursor run (see Sect.[wiki:doc/tec/1dmodel 1-D model for precursor runs]). Uniformly distributed random perturbations with a user-defined amplitude can be imposed to the fields of the horizontal velocities components to initiate turbulence. == Laminar and turbulent inflow boundary conditions == In case of laminar inflow, Dirichlet boundary conditions are used for all quantities, except for the SGS-TKE ''e'' and perturbation pressure ''π^∗^'' for which Neumann boundary conditions are used. Vertical profiles, as taken for the initialization of the simulation, are used for the Dirichlet boundary conditions. In order to allow for a fast turbulence development, random perturbations can be imposed on the velocity fields within a certain area behind the inflow boundary (inlet). These perturbations may persist for the entire simulation. For the purpose of preventing gravity waves from being reflected at the inlet, a relaxation area can be defined after [#davies1976 Davies (1976)]. So far, it was found to be sufficient to implement this method for temperature only. This is hence realized by an additional term in the prognostic equation for ''θ'' (see third equation in Sect. [wiki:/doc/tec/gov governing equation]): {{{ #!Latex \begin{align*} \frac{\partial \theta}{\partial t} = \ldots - C_{\text{relax}} \left(\theta - \theta_{\text{inlet}}\,\right). \end{align*} }}} Here, ''θ'',,inlet,, is the stationary inflow profile of ''θ'', and ''C'',,relax,, is a relaxation coefficient, depending on the distance ''d'' from the inlet, viz. {{{ #!Latex \begin{align*} &C_{\text{relax}}(d) = \begin{cases} F_{\text{inlet}} \cdot \sin^2 \left(\frac{\pi}{2} \frac{D - d}{D} \right) & \text{for~} d < D,\\ 0 & \text{for~} d \ge D,\\ \end{cases} \end{align*} }}} with ''D'' being the length of the relaxation region and ''F'',,inlet,, being a damping factor. == Turbulence recycling == If non-cyclic horizontal boundary conditions are used, PALM offers the possibility of generating time-dependent turbulent inflow data by using a turbulence recycling method. The method follows the one described by [#lund1998 Lund et al. (1998)], with the modifications introduced by [#kataoka2002 Kataoka and Mizuno (2002)]. Figure 2 gives an overview of the recycling method used in PALM. [[Image(03.png,600px,border=1)]] Figure 2: Schematic figure of the turbulence recycling method used for generation of turbulent inflow. The configuration represents exemplary conditions with a built-up analysis area (brown surface) and an open water recycling area (blue surface). The blue arrow indicates the flow direction. The turbulent signal ''φ^'^(y, z, t)'' is taken from a recycling plane which is located at a fixed distance ''x'',,recycle,, from the inlet: {{{ #!Latex \begin{align*} & \varphi^{\prime}(y, z, t) = \varphi(x_{\text{recycle}},y, z, t) - \langle \varphi\rangle_y(z, t), \end{align*} }}} where ''<φ>,,y,,(z, t)'' is the line average of a prognostic variable ''φ ∈ {u, v, w, θ, e}'' along ''y'' at ''x = x,,recycle,,''. ''φ^'^(y, z, t)'' is then added to the mean inflow profile ''<φ,,inflow,,>,,y,,(z)'' at ''x,,inlet,,'' after each time step: {{{ #!Latex \begin{align*} & \varphi_{\text{inlet}}(y, z, t) = \langle \varphi_{\text{inlet}}\rangle_y(z) + \phi(z) \varphi^{\prime}(y, z, t), \end{align*} }}} with the inflow damping function ''Φ(z)'', which has a value of ''1'' below the initial boundary layer height, and which is linearly damped to ''0'' above, in order to inhibit growth of the boundary layer depth. ''<φ,,inlet,,>,,y,,(z)'' is constant in time and either calculated from the results of the precursor run or prescribed by the user. The distance ''x'',,recycle,, has to be chosen much larger than the integral length scale of the respective turbulent flow. Otherwise, the same turbulent structures could be recycled repeatedly, so that the turbulence spectrum is illegally modified. It is thus recommended to use a precursor run for generating the initial turbulence field of the main run. The precursor run can have a comparatively small domain along the horizontal directions. In that case the domain of the main run is filled by cyclic repetition of the precursor run data. Note that the turbulence recycling has not been adapted for humidity and passive scalars so far. Turbulence recycling is frequently used for simulations with urban topography. In such a case, topography elements should be placed sufficiently downstream of ''x'',,recycle,, to prevent effects on the turbulence at the inlet. == Open outflow boundary conditions == At the outflow boundary (outlet), the velocity components ''u,,i,,'' meet radiation boundary conditions, viz. {{{ #!Latex \begin{align*} \frac{\partial u_i}{\partial t} + U_{u_i} \frac{\partial u_i}{\partial n} = 0\,, \end{align*} }}} as proposed by [#orlanski1976 Orlanski (1976)]. Here ''∂/∂n'' is the derivative normal to the outlet (i.e., ''∂/∂x'' in Figure 2 and ''U,,ui,,'' a transport velocity which includes wave propagation and advection. Rewriting the equation above yields the transport velocity {{{ #!Latex \begin{align*} U_{u_i} = -\left(\frac{\partial u_i}{\partial t}\right)\left(\frac{\partial u_i}{\partial n}\right)^{-1} \end{align*} }}} that is calculated at interior grid points next to the outlet at the preceding time step for each velocity component. If the transport velocity, calculated by means of the equation for the transport velocity, is outside the range ''0 ≤ U,,ui,, ≤ '''Δ'''/Δt'', it is set to the respective threshold value that is exceeded. Because this local determination of ''U,,ui,,'' can show high variations in case of complex turbulent flows, it is averaged laterally to the direction of the outflow, so that it varies only in the vertical direction. Alternatively, the transport velocity can be set to the upper threshold value (''U,,ui,, = '''Δ'''/Δt'') for the entire outlet. Both equations mentioned in this section are discretized using an upstream method following [#miller1981 Miller and Thorpe (1981)]. As the radiation boundary condition does not ensure conservation of mass, a mass flux correction can be applied at the outlet. == Topography == = References = * [=#holtslag] '''Holtslag AAM, Bruin HARD.''' 1988. Applied modelling of the night-time surface energy balance over land. J. Appl. Meteorol. 27: 689–704. * [=#panofsky] '''Panofsky HA, Dutton JA.''' 1984. Atmospheric Turbulence, Models and Methods for Engineering Applications. John Wiley & Sons. New York. * [=#paulson] '''Paulson CA''' 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J. Appl. Meteorol. 9: 857–861. * [=#klemp1978] '''Klemp JB, Lilly DK.''' 1978. Numerical simulation of hydrostatic mountain waves. J. Atmos. Sci. 35: 78–107. * [=#davies1976] '''Davies HC.''' 1976. A lateral boundary formulation for multi-level prediction models. Q. J. Roy. Meteor. Soc. 102: 405–418. * [=#lund1998] '''Lund TS, Wu X, Squires KD.''' 1998. Generation of turbulent inflow data for spatially-developing boundary layer simulations. J. Comput. Phys. 140: 233–258. * [=#kataoka2002] '''Kataoka H, Mizuno M.''' 2002. Numerical flow computation around aerolastic 3d square cylinder using inflow turbulence. Wind Struct. 5: 379–392. * [=#orlanski1976] '''Orlanski I.''' 1976. A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys. 21: 251–269. * [=#miller1981] ''' Miller MJ, Thorpe AJ.''' 1981. Radiation conditions for the lateral boundaries of limited-area numerical models. Q. J. Roy. Meteor. Soc. 107: 615–628.