Surface Layer Boundaries
The surface layer provides an abstraction layer for users to directly specify diffusive fluxes at a boundary via a multitude of methods. More specifically, the surface layer condition applies an impenetrable condition for the boundary normal velocity but higher order extrapolation for all the other variables and then allows a user to specify a method for calculating the diffusive fluxes
erf.surface_layer.flux_type = STRING #flux types (moeng, custom, bulk_coeff)
The bulk_coeff flux type employs bulk drag coefficients to compute the diffusive stresses while the moeng type
employs Moeng’s formulation for Monin-Obukhov similarity theory (MOST) and custom allows the user to directly
specify the fluxes through ustar; tstar; qstar. Currently, the MOST pathway is the primary flux type employed
in ERF simulations and will be the focus in subsequent sections.
MOST Theory
Monin-Obukhov similarity theory (MOST) is used to describe the atmospheric surface layer (ASL), the lowest part of the atmospheric boundary layer. The implementation of MOST in ERF follows that in AMR-Wind, which is based on the surface layer profiles presented in P. van der Laan, et al., Wind Energy, 2017 and D. Etling, “Modeling the vertical ABL structure”, 1999. MOST theory assumes that the ASL is in a steady state and horizontally homogeneous, and kinematic fluxes due to turbulent transport (\(\overline{u^{'}w^{'}}\), \(\overline{v^{'}w^{'}}\), and \(\overline{\theta^{'}w^{'}}\)) are constant with height. \(\Phi_m\) and \(\Phi_h\) are the nondimensional wind shear and temperature gradient, respectively, which are assumed to follow universal similarity laws based on dimensional arguments. With these assumptions, the MOST theory can be written as:
where the nondimensional gradients are expressed in terms of the MOST stability parameter, \(\zeta = \frac{z}{L} = -\frac{\kappa z}{u_{\star}^{3}} \frac{g}{\overline{\theta}} \overline{w^{'}\theta^{'}}\), which serves as a surface layer scaling parameter. Here, \(L\) is the Monin-Obukhov length, \(u_{\star}\) is the friction velocity (defined for \(u\) aligned with the wind direction), \(\theta_{\star}\) is the surface layer temperature scale, \(\overline{\theta}\) is the reference virtual potential temperature for the ASL, and \(\kappa\) is the von Karman constant (taken to be \(0.41\)).
Integration of the MOST assumption equations give the classical MOST profiles of mean velocity, potential temperature, and water vapor
where \(\theta_0\) and \(q_{v,0}\) are the surface values and \(z_0\) is a characteristic roughness height. The integrated similarity functions,
are calculated analytically from empirical gradient functions \(\Phi_m\) and \(\Phi_h\), which are defined piecewise for stable and unstable values of the stability parameter.
Unstable: \((\zeta < 0)\)
Stable: \((\zeta > 0)\)
where the constants take the values proposed in Dyer, Boundary Layer Meteorology, 1974:
Inverting the equations above, the MOST stability parameter,
is determined by the friction velocity
and the characteristic surface layer temperature
MOST Implementation
As noted in Surface Layer Boundaries, the boundary conditions for velocity, temperature, and water vapor do not change with the flux type. Therefore, the MOST implementation in ERF is a specific method for computing the diffusive fluxes, which are directly written into the stress tensor/vector. The MOST pathway is structured to allow either the surface temperature (\(\theta_0\)) or surface temperature flux (\(\overline{w^{'}\theta^{'}}\)) to be enforced. To compute the MOST flux, the following algorithm is applied:
Horizontal (planar) averages \(\bar{u}\), \(\bar{v}\) and \(\overline{\theta}\) are computed at a reference height \(z_{ref}\) assumed to be within the surface layer.
Initially, neutral conditions (\(L=\infty, \zeta=0\)) are assumed and used to compute a provisional \(u_{\star}\) using the equation given above. If \(\theta_0\) is specified, the above equation for \(\theta_{\star}\) is applied and then the surface flux is computed \(\overline{w^{'}\theta^{'}} = -u_{\star} \theta_{\star}\). If \(\overline{w^{'}\theta^{'}}\) is specified, \(\theta_{\star}\) is computed as \(-\overline{w^{'}\theta^{'}}/u_{\star}\) and the previous equation is inverted to compute \(\theta_0\).
The stability parameter \(\zeta\) is recomputed using the equation given above based on the provisional values of \(u_{\star}\) and \(\theta_{\star}\).
The previous two steps are repeated iteratively, sequentially updating the values of \(u_{\star}\) and \(\zeta\), until the change in the value of \(u_{\star}\) on each iteration falls below a specified tolerance.
Once the MOST iterations have converged, and the planar average surface flux values are known, the approach from Moeng, Journal of the Atmospheric Sciences, 1984 is applied to consistently compute local surface-normal stress/flux values (e.g., \(\tau_{xz} = - \rho \overline{u^{'}w^{'}}\)):
\[ \begin{align}\begin{aligned}\left. \frac{\tau_{xz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(u - \bar{u})\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} + \bar{u}\sqrt{u^2 + v^2} }{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}}^2},\\\left. \frac{\tau_{yz}}{\rho} \right|_0 &= u_{\star}^{2} \frac{(v - \bar{v})\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} + \bar{v}\sqrt{u^2 + v^2} }{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}}^2},\\\left. \frac{\tau_{\theta z}}{\rho} \right|_0 &= \theta_\star u_{\star} \frac{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} ({\theta} - \overline{\theta}) + \sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ \overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} (\overline{\theta} -\theta_0) } = u_{\star} \kappa \frac{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} ({\theta} - \overline{\theta}) + \sqrt{u^2+v^2} (\overline{\theta} - \theta_0) }{ \overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} [ \mathrm{ln}(z_{ref} / z_0)-\Psi_{h}(z_{ref}/L)] }\\\left. \frac{\tau_{q_{v} z}}{\rho} \right|_0 &= q_{v,\star} u_{\star} \frac{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} ({q_{v}} - \overline{q_{v}}) + \sqrt{u^2+v^2} (\overline{q_{v}} - q_{v,0}) }{ \overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} (\overline{q_{v}} - q_{v,0}) } = u_{\star} \kappa \frac{\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} ({q_{v}} - \overline{q_{v}}) + \sqrt{u^2+v^2} (\overline{q_{v}} - q_{v,0}) }{ \overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}} [ \mathrm{ln}(z_{ref} / z_0)-\Psi_{h}(z_{ref}/L)] }\end{aligned}\end{align} \]where \(\bar{\psi}\) are plane averaged values (at \(z_{ref}\)) and \(\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}}\) is the plane averaged magnitude of horizontal velocity (plane averaged wind speed). We note a slight variation in the denominator of the velocity terms from the form of the equations presented in Moeng. This difference is due to how the stress components are computed — i.e., the stress componen in Moeng’s work, \(\tau_{xz,0}\), is given by \(\tau_{xz,0} = u_{\star}^{2} \bar{u}/\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}}\), where \(\bar{u}/\overline{{\scriptstyle|}\mathbf{u}{\scriptstyle|}}\) is the unit vector component applied to the total stress \(u_{\star}^{2}\).
Finally, it must be noted that using terrain-fitted coorindates will modify the surface normal and tangent vectors. Consequently, the MOST implementation with terrain-fitted coorindates will formally require local vector rotations. Stress rotations with MOST are a work in progress but may be activated with
erf.use_rotate_surface_flux = true. Therefore, running with terrain (erf.terrain_type = StaticFittedMesh) and with MOST (surface_flux.flux_type = "moeng") is not recommended.
MOST Inputs
To evaluate the fluxes with MOST, the surface rougness parameter \(z_{0}\) must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity \(u_{\star}\). ERF supports five methods for parameterizing the surface roughness: constant, charnock, modified_charnock, donelan, and wave_coupled. The latter four methods parameterize \(z_{0} = f(u_{\star})\) and are described in Jimenez & Dudhia, American Meteorological Society, 2018 and Warner et. al, Ocean Modelling, 2010. The rougness calculation method may be specified with
erf.most.roughness_type_land = STRING # Z_0 type over land (constant)
erf.most.roughness_type_sea = STRING # Z_0 type over water (charnock, coare3.0,
# donelan, modified_charnock, wave_coupled,
# constant)
Land and water are set separately. The older single erf.most.roughness_type
input was replaced by these two and now aborts the run if it is present. The
complete list of surface-layer inputs, with acceptable values and defaults, is in
Surface Layer (MOST).
If the charnock method is employed, the \(a\) constant may be specified with erf.most.charnock_constant (defaults to 0.0185). If the modified_charnock method is employed, the depth \(d\) may be specified with erf.most.modified_charnock_depth (defaults to 30 m). If the wave_coupled method is employed, the user must provide wave height and mean wavelength data.
While the MOST methods relevant to air-sea interfaces (charnock, modified_charnock, and wave_coupled) dynamically compute \(z_{0}\), the constant case can also allow for spatially varying \(z_{0}\) if an inputs file is employed. More specifically, one may specify
erf.most.roughness_file_name = STRING #Name of file that contains (x,y,z_0)
in the inputs file and ERF will populate the 2D \(z_{0}\) array with values contained in the text file.
When computing an average \(\overline{\phi}\) for the MOST boundary, where \(\phi\) denotes a generic variable, ERF supports a variety of approaches. Specifically, planar averages and local region averages may be computed with or without time averaging. With each averaging methodology, the query point \(z\) may be determined from the following procedures: specified distance \(z_{ref}\) inward from the selected boundary, specified \(k_{index}\), or (when employing terrain-fitted coordinates) specified normal vector length \(z_{ref}\). The available inputs to the MOST boundary and their associated data types are
erf.most.average_policy = INT #POLICY FOR AVERAGING
erf.most.use_normal_vector = BOOL #USE NORMAL VECTOR W/ TERRAIN?
erf.most.use_interpolation = BOOL #INTERPOLATE QUERY POINT W/ TERRAIN?
erf.most.time_average = BOOL #USE TIME AVERAGING?
erf.most.z0 = FLOAT #SURFACE ROUGHNESS [m]
erf.most.zref = FLOAT #QUERY DISTANCE (HEIGHT OR NORM LENGTH) [m]
erf.most.pot_surf_temp = FLOAT #SPECIFIED SURFACE POTENTIAL TEMP [K]
erf.most.abs_surf_temp = FLOAT #SPECIFIED SURFACE ABSOLUTE TEMP [K]
erf.most.surf_temp_flux = FLOAT #SPECIFIED SURFACE TEMP FLUX [K-m/s]
erf.most.surf_heating_rate = FLOAT #SPECIFIED RATE OF SURFACE TEMP CHANGE [K/h]
erf.most.surf_moist = FLOAT #SPECIFIED SURFACE MOISTURE [-]
erf.most.surf_moist_flux = FLOAT #SPECIFIED SURFACE MOISTURE FLUX [m/s]
erf.most.k_arr_in = INT #SPECIFIED K INDEX ARRAY (MAXLEV)
erf.most.radius = INT #SPECIFIED REGION RADIUS [grid cells]
erf.most.time_window = FLOAT #WINDOW FOR TIME AVG [s]
On the surface-temperature moisture path, when no positive
most.surf_moist is available, the surface moisture used in the
thermodynamic MOST state falls back to the atmospheric value.
When neither erf.most.zref nor erf.most.k_arr_in is given, a mesh with terrain-fitted coordinates queries \(z_{ref} = 10\) m above the surface, while a mesh without them uses the center of the first cell. On a vertically stretched mesh (erf.grid_stretching_ratio or erf.terrain_z_levels) that center, and the cell holding a specified \(z_{ref}\), come from the actual z levels rather than from a uniform spacing. A query height that lies exactly on a cell face belongs to the cell above the face, or is interpolated across it when erf.most.use_interpolation is true. The height in use is printed at start-up (MOST reference height at level ...). It is always the height above the surface; on the low z face of a constant-dz mesh, however, a value given in erf.most.zref is read as an absolute z, the height plus geometry.prob_lo[2], and erf.most.k_arr_in as an absolute cell index.
We now consider two concrete examples. To employ an instantaneous planar average at a specified vertical height above the bottom surface, one would specify:
erf.most.average_policy = 0
erf.most.use_normal_vector = false
erf.most.time_average = false
erf.most.z0 = 0.1
erf.most.zref = 1.0
erf.most.surf_temp_flux = 0.0
erf.most.surf_moist_flux = 0.0
By contrast, local region averaging would be employed in conjunction with time averaging for the following inputs:
erf.most.average_policy = 1
erf.most.use_normal_vector = true
erf.most.use_interpolation = true
erf.most.time_average = true
erf.most.z0 = 0.1
erf.most.zref = 1.0
erf.most.surf_temp_flux = 0.0
erf.most.surf_moist_flux = 0.0
erf.most.radius = 1
erf.most.time_window = 10.0
In the above case, use_normal_vector utilizes the a local surface-normal vector with length \(z_{ref}\) to construct the positions of the query points. Each query point, and surrounding points that are within erf.most.radius from the query point, are interpolated to and averaged; for a radius of 1, 27 points are averaged. The time average is completed by way of an exponential filter function whose peak coincides with the current time step and tail extends backwards in time
Due to the form of the above integral, it is advantageous to consider \(\tau\) as a multiple of the simulation time step \(\Delta t\), which is specified by erf.most.time_window. As erf.most.time_window is reduced to 0, the exponential filter function tends to a Dirac delta function (prior averages are irrelevant). Increasing erf.most.time_window extends the tail of the exponential and more heavily weights prior averages.
The state of this filter is written to and read from checkpoint files, so a restarted run continues the average rather than beginning it again from the instantaneous value; see Checkpoint / Restart. Note that regridding rebuilds the averaging containers and therefore restarts the filter on the affected level.
Low-speed corrections
The following options are available:
erf.most.include_wstar = true
erf.most.include_subgrid_vel = true
erf.most.pblh_calc = MYNN25
These correspond to a mean surface velocity of
where \(w_*\) is the (Deardorff) convective velocity scale and
\(\beta=1.2\) (Beljaars 1995, QJRMS). Because \(w_*\) is computed from the
boundary layer depth, erf.most.include_wstar requires erf.most.pblh_calc to
name a PBL height scheme; leaving it at its default of none is an error, since no
boundary layer depth would ever be diagnosed.
On a refined level, the MYNN25 diagnostic runs on the level’s own columns only when its
grids run from the ground to the top of the domain in every column. A refined region that
ends lower (a patch near the ground, say) cannot see the top of a boundary layer deeper than
itself, and one that does not reach the ground (a region aloft) has no column to scan, so
either takes the PBL height of the next coarser level at each of its columns instead: each
fine column gets the height of the coarse column that holds it. With subcycling that is the
coarse height of the start of the coarse step. With EB terrain every level still diagnoses
its own.
The subgrid velocity scale \(V_{sg}\) handles weak large-scale flow that is underresolved (Mahrt & Sun 1995, MWR). This is parameterized as
which vanishes for grid spacings of \(\Delta x < 5\) km.
The MYNN25 PBL height is found by scanning each column upwards from the ground, so
it does not depend on how the grids are decomposed: where boxes are stacked in z, the
scan runs on the whole column they make up. The scan ends at the top of the grids, so
on a refined level that stops below the top of the domain a PBL top above the refined
region is not found. The PBL height is zero wherever it is not found, and wherever
the grids of a level do not reach the ground. With erf.terrain_type = EB the PBL
height is stored on the three-dimensional grids, and every grid must start at the
bottom of the domain (no decomposition in z); the run aborts otherwise.
Surface Layer boundary on all faces
While the surface layer boundary condition is primarily used on the bottom z face (zlo), it can also be applied to all other domain sides. This can allow MOST theory to be used on faces for simulating closed environments (e.g, cloud chamber).
To use the surface layer boundary on different faces, all existing MOST options are prefixed with the corresponding face name (xlo, ylo, zlo, xhi, yhi zhi). For example,
xlo.type = "surface_layer"
ylo.type = "surface_layer"
erf.xlo.surface_layer.flux_type = MOENG
erf.xlo.most.pot_surf_temp = 289.5
erf.xlo.most.surf_moist = 0.00765
erf.xlo.most.z0 = 0.000035
erf.xlo.most.zref = 0.03125
erf.xlo.most.average_policy = 1
erf.xlo.most.radius = 0
erf.ylo.surface_layer.flux_type = MOENG
erf.ylo.most.pot_surf_temp = 289.5
erf.ylo.most.surf_moist = 0.00765
erf.ylo.most.z0 = 0.000035
erf.ylo.most.zref = 0.03125
erf.ylo.most.average_policy = 1
erf.ylo.most.radius = 0
The face name is omitted when zlo is the only surface-layer face, so an input
file written before this capability existed keeps working unchanged: with
zlo.type = "surface_layer" and no other surface-layer face, erf.most.* and
erf.surface_layer.* are read as before.
As soon as a second surface-layer face is declared, every face is expected to
qualify its inputs, zlo included. Unqualified inputs are still honored in
that case – they are applied to zlo, with a warning asking that they be
migrated to erf.zlo.most.* and erf.zlo.surface_layer.* – so that adding a
wall to a working deck cannot silently fall back to the built-in defaults. Giving
both spellings at once (erf.most.z0 together with erf.zlo.most.z0, say) is
an error and aborts, since there is no way to tell which one was meant.
Note that not all existing options are supported when using MOST on other faces (such as interpolation and time averaging).
Currently the MOENG flux type is supported on all faces; bulk_coeff,
custom and rico abort anywhere else. Of the surface thermodynamic
pathways, only the prescribed surface potential temperature (most.pot_surf_temp) carries the
orientation through to the flux, so a prescribed surface heat flux
(most.surf_temp_flux), the fully adiabatic case, and any variable sea
roughness (most.roughness_type_sea other than constant over water cells)
abort on a face other than zlo rather than run with the zlo sign convention. The
MOST PBL-height diagnostic (most.pblh_calc) and the free-convection correction
it feeds (most.include_wstar) are likewise zlo-only.
Lateral (x/y) surface-layer boundaries currently require every grid on that
level to span its full vertical extent. Grids decomposed in z and
partial-height refined grids are not supported; configure
erf.max_grid_size_z accordingly.
Deardorff LES, RANS, and PBL models, including SHOC, support SurfaceLayer
only at zlo; they cannot be combined with lateral or upper surface-layer
boundaries. Smagorinsky LES, including its 2-D variant, uses the generic
all-face diffusion path and is supported.
Prescribed time-varying surface forcing
Instead of computing fluxes directly through MOST, the surface land/ocean fluxes can be prescribed over time from a text file.
Additionally, the sea-surface temperature can be prescribed while allowing MOST to compute the surface fluxes.
For prescribing surface fluxes from H, LE, and TAU:
erf.most.use_sfc_fluxes = true # flag for prescribing surface fluxes
erf.most.sfc_file = sfc.txt # path to file containing surface data
where erf.most.sfc_file is the path to a text file with the following example format:
day sst(K) H(W/m2) LE(W/m2) TAU(m2/s2)
205.500 293.65 4.43 22.73 0.00
205.542 295.63 43.55 88.67 0.00
205.583 298.40 104.08 154.39 0.00
H and LE are dynamic fluxes and are applied as \(H/c_{p}\) and \(LE/L_{v}\);
TAU is the kinematic stress \(u_{\star}^{2}\) and is multiplied by the local density, so the
applied surface stress is \(\rho \, \mathrm{TAU}\).
This mode implies erf.most.flux_type = custom and is not compatible with an active land surface model, EB,
or erf.most.rhosurf.
The existing text-file forcing restriction with EB remains in effect.
For prescribing sea-surface temperature:
erf.most.use_sfc_sst = true # flag for prescribing sea-surface temperature
erf.most.sfc_file = sfc.txt # path to file containing surface data
where erf.most.sfc_file is the path to a text file with the same format as above,
or purely time and SST information, such as:
day sst(K)
203.000000 292.800
204.000000 294.580
The SST column in this file is absolute temperature in kelvin. ERF converts
it to the SurfaceLayer’s canonical potential-temperature field using the
physical surface pressure diagnosed from the lowest atmospheric cell. In
contrast, erf.most.pot_surf_temp and erf.most.surf_heating_rate are already
specified in the SurfaceLayer potential-temperature convention, while
erf.most.abs_surf_temp (with erf.most.surf_heating_rate then a rate of change of
the absolute temperature) is converted every step with the same diagnosed surface
pressure. erf.most.surf_temp is the deprecated spelling of erf.most.pot_surf_temp. SST and
skin-temperature data from WRFInput, wrflowinp, and Metgrid are normalized to
that same convention at ingestion; coupled SST is converted when it enters
the SurfaceLayer, while Noah-MP’s absolute t_sfc remains owned by the
radiation path and is not adopted as the SurfaceLayer field.
Notes
With erf.terrain_type = EB, SurfaceLayer does not support planar SST/TSK
from real-data input, coupled-ocean SST, or LSM surface-temperature or
surface-flux providers because those LSM flux arrays are planar and the helper
only writes their k=0 slab, while EB MOST consumes surface parameters on
arbitrary cut-cell k; no 2-D-to-cut-cell mapping is defined. Custom/file-
driven roughness (most.roughness_file_name) is also unsupported under EB
for the same planar-provider compatibility policy. Ordinary EB SurfaceLayer
operation using SurfaceLayer-native prescribed or default state remains
supported.
For both flux and sst modes, the time column is reset relative to the simulation elapsed time (t=0s). This might cause issues with simulations initialized with real data, such as WRFInput.
The sst mode only applies the surface temperature where there are ocean cells (landmask=0). This can be forced with
erf.is_land = 0.The surface fluxes are currently applied regardless of land mask (both ocean and land).