A Revised Scheme for the WRF Surface Layer Formulation
PEDRO A. JIME
´
NEZ,* JIMY DUDHIA,
1
J. FIDEL GONZA
´
LEZ-ROUCO,
#
JORGE NAVARRO,
@
JUAN P. MONTA
´
VEZ,
&
AND ELENA GARCI
´
A-BUSTAMANTE**
* Divisio
´
n de Energı´as Renovables, CIEMAT, Madrid, Spain, and Mesoscale and Microscale Meteorology Division,
NCAR, Boulder, Colorado
1
Mesoscale and Microscale Meteorology Division, NCAR, Boulder, Colorado
#
Departamento de Astrofı´sica y Ciencias de la Atmo
´
sfera, Universidad Complutense de Madrid, Madrid, Spain
@
Divisio
´
n de Energı´as Renovables, CIEMAT, Madrid, Spain
&
Departamento de Fı´sica, Universidad de Murcia, Murcia, Spain
** Departamento de Astrofı´sica y Ciencias de la Atmo
´
sfera, Universidad Complutense de Madrid, and Divisio
´
n
de Energı´as Renovables, CIEMAT, Madrid, Spain
(Manuscript received 24 February 2011, in final form 17 August 2011)
ABSTRACT
This study summarizes the revision performed on the surface layer formulation of the Weather Research
and Forecasting (WRF) model. A first set of modifications are introduced to provide more suitable similarity
functions to simulate the surface layer evolution under strong stable/unstable conditions. A second set of
changes are incorporated to reduce or suppress the limits that are imposed on certain variables in order to
avoid undesired effects (e.g., a lower limit in u
*
). The changes introduced lead to a more consistent surface
layer formulation that covers the full range of atmospheric stabilities. The turbulent fluxes are more (less)
efficient during the day (night) in the revised scheme and produce a sharper afternoon transition that shows
the largest impacts in the planetary boundary layer meteorological variables. The most important impacts in
the near-surface diagnostic variables are analyzed and compared with observations from a mesoscale network.
1. Introduction
The lowest part of the planetary boundary layer (PBL)
wherein the turbulent fluxes vary less than 10% of their
magnitude is known as the atmospheric surface layer
(Stull 1988). Meteorological variables experience a sharp
variation with height within this layer that exhibits the
most significant exchanges of momentum, heat, and
moisture (Arya 1988). The surface layer state determines
the land–atmosphere interaction and, thus, its accurate
formulation is crucial to provide an adequate atmo-
spheric evolution by numerical models.
The Monin–Obukhov similarity theory (Obukhov 1946;
Monin and Obukhov 1954) is a widely used framework
to compute the surface turbulent fluxes (Beljaars and
Holtslag 1991). The theory also provides information
of the profiles within the surface layer that are used
to diagnose meteorological variables at their typical
observational height such as the wind at 10 m or the
temperature and moisture at 2 m. A limitation, however,
is that the predicted similarity functions (f
h,m
)necessary
to compute both the fluxes and the profiles need to be
determined empirically.
The Kansas field program (Izumi 1971) provided esti-
mations of the similarity functions for a limited range of
atmospheric stabilities (Businger et al. 1971; Dyer 1974;
Hicks 1976). For this reason, extensions to highly stable
situations (e.g., Webb 1970; van Ulden and Holtslag 1985;
Holtslag and de Bruin 1988; Beljaars and Holtslag
1991; Cheng and Brutsaert 2005) as well as highly un-
stable conditions (e.g., Brutsaert 1992; Fairall et al. 1996;
Grachev et al. 2000; Wilson 2001; Fairall et al. 2003) have
been proposed. For instance, Fairall et al. (1996, hereafter
F96) used the asymptotic behavior predicted by the
theory to extend the Kansas type of similarity functions
to higher instabilities. The proposed similarity functions
are therefore valid from neutral to free convective situ-
ations. Similarly, Cheng and Brutsaert (2005, hereafter
CB05) found an asymptotic behavior of the similarity
functions for the stable part and derived functions
valid from neutral situations to very stable conditions.
Corresponding author address: Pedro A. Jimenez, Mesoscale and
Microscale Meteorology Division, National Center for Atmo-
spheric Research, 3450 Mitchell Ln., Boulder, CO 80301.
E-mail: jimenez@ucar.edu
898 MONTHLY WEATHER REVIEW VOLUME 140
DOI: 10.1175/MWR-D-11-00056.1
Ó 2012 American Meteorological Society
Hence, combining similarity functions such as the ones
proposed by F96 and CB05 allows one to cover more
accurately the full range of atmospheric stabilities by the
Monin–Obukhov similarity theory.
The purpose of this investigation is to improve the
surface layer formulation of the Weather Research
and Forecast model (WRF; Skamarock et al. 2008), in
particular, the surface layer scheme based on the fifth-
generation Pennsylvania State University–National Cen-
ter for Atmospheric Research Mesoscale Model (MM5)
parameterization (Grell et al. 1994). Although the scheme
is widely used for quite different atmospheric investi-
gations (e.g., Weisman et al. 2008; Jingyong et al. 2008;
Jime
´
nez et al. 2010a), it uses Kansas-type similarity
functions with their limited coverage of atmospheric
stabilities. Here, the similarity functions are replaced
by those proposed by F96 and CB05 in order to provide
the scheme with a more appropriate framework for
strongly stable/unstable conditions.
An additional target of this investigation is to review the
limits that are imposed to certain variables (e.g., u
*
,the
friction velocity) in order to prevent undesirable effects
from the formulation. It has been found herein that these
limits can be reduced or removed in order to provide a less
restrictive, and more consistent, surface layer formulation
that covers the full range of atmospheric stabilities.
The impact that these changes produce in the surface
fluxes, the diagnostic surface meteorological variables,
and the PBL dynamics is analyzed. The most important
impacts in the near-surface variables have been tested
against observations from a mesoscale network located
in the northeast of the Iberian Peninsula (Jime
´
nez et al.
2010b). The observations have been used in previous
studies with quite different orientations (e.g., Jime
´
nez
et al. 2009a; Garcı
´
a-Bustamante et al. 2011). A complete
summer season has been simulated herein (at a high
horizontal resolution of 2 km) in order to obtain a statis-
tically robust characterization of the changes introduced
by the new formulation. The standard WRF model
output has been complemented by recording the surface
layer variables every time step at the observational sites
in order to provide a detailed evolution of the atmo-
sphere within this layer.
2. Surface layer parameterization
This section describes the current WRF surface layer
formulation (surface layer physics option 1 in WRF,
section 2a), its limitations (section 2b), and the revised
formulation here proposed in order to overcome the
current problems (section 2c).
The computation of the fluxes in WRF not only de-
pends on the surface layer physics but also in the land
surface model. The present description of fluxes follows
the definition of the soil scheme used in this investigation
(surface-physics option 1 in WRF; Dudhia 1996; Dudhia
et al. 2004), but the changes described apply equally to
other land surface options. More details of the WRF
configuration used in this investigation will be provided
in section 3.
a. WRF surface layer
The surface layer is assumed to be the first vertical
layer and the surface fluxes are parameterized as follows:
t 5 ru
2
*
5 rC
d
U
2
, (1)
H 52rc
p
u
*
u
*
52rc
p
C
h
U(u
a
2 u
g
), (2)
LH 5 L
e
ru
*
q
*
5 L
e
rMC
q
U(q
g
2 q
a
), (3)
where t, H, and LH are the fluxes of momentum, sen-
sible heat, and latent heat, respectively; u
*
and q
*
are
the temperature and moisture scales, respectively; r is
the air density in the surface layer; c
p
is the specific heat
capacity at constant pressure; and U is the wind speed in
the lower layer enhanced by a convective velocity fol-
lowing Beljaars (1995) and a subgrid velocity following
Mahrt and Sun (1995). This last correction only applies
for horizontal grid resolutions higher than 5 km. Here L
e
is the latent heat of vaporization; M is the soil moisture
availability; u
a
and u
g
are the air and ground surface
potential temperature, respectively; q
g
is the saturated
specific humidity at the ground; q
a
is the specific hu-
midity in the surface layer; and C
d
, C
h
, and C
q
are the
dimensionless bulk transfer coefficients (Stull 1988) for
momentum, heat, and moisture, respectively.
The Monin–Obukhov similarity theory is used to cal-
culate the transfer coefficients. The dimensionless wind
shear and potential temperature gradient are usually
expressed as (e.g., Arya 1988)
kz
u
*
›u
a
›z
5 f
m
z
L
;
kz
u
*
›u
›z
5 f
h
z
L
,
where k 5 0.4 is the von Ka
´
rma
´
n constant, u
a
is the wind
speed at level z, and L is the Obukhov length (Obukhov
1946). Integrating the equations with respect to height z,
leads to
u
a
5
u
*
k
ln
z
z
0
2 c
m
z
L
1 c
m
z
0
L
, (4)
(u
a
2 u
g
) 5
u
*
k
ln
z
z
0
2 c
h
z
L
1 c
h
z
0
L
, (5)
MARCH 2012 J I M E
´
NEZ ET AL. 899
where z
0
is the roughness length and c
m,h
are the in-
tegrated similarity functions for momentum and heat
that are defined as follows (e.g., Panofsky 1963):
c
m,h
z
L
[
ð
z/L
0
[1 2 f
m,h
(z)]
dz
z
.
Combining Eq. (1) and Eq. (4) and neglecting the
contribution of c
m
(z
0
/L) allows one to obtain the bulk
transfer coefficient for momentum:
C
d
5
k
2
ln
z
z
0
2 c
m
z
L
2
. (6)
Analogously, combining Eq. (2) with Eqs. (4) and (5),
and neglecting the contribution of c
h
(z
0
/L) allows one
to obtain the bulk transfer coefficient for heat:
C
h
5
k
2
ln
z
z
0
2 c
m
z
L
ln
z
z
0
2 c
h
z
L
, (7)
where it has been assumed that u
a
5 U.
For the case of moisture the surface layer formulation
follows Carlson and Boland (1978). The existence is as-
sumed of a viscous sublayer from the ground to a height
z
l
(z
l
5 0.01 m over land and z
0
over water), and a tur-
bulent layer wherein Monin–Obukhov theory is appli-
cable from z
l
to z. A similar derivation to the one used to
obtain the transfer coefficients for momentum and heat
leads to obtain the value of the bulk transfer coefficient
for moisture (Carlson and Boland 1978; Grell et al. 1994):
C
q
5
k
2
ln
z
z
0
2 c
m
z
L
"
ln
rc
p
ku
*
z
c
s
1
z
z
l
!
2 c
h
z
L
#
,
(8)
where c
s
is the effective heat transfer coefficient for
nonturbulent processes. Note that it has been assumed
that the dimensionless similarity function for moisture is
the same as heat. This hypothesis is based on experi-
mental evidence (e.g., Dyer 1967; Dyer and Bradley
1982) but it has been recently questioned (e.g., Park
et al. 2009).
The integrated similarity functions are calculated
according to four stability regimes (Zhang and Anthes
1982) defined in terms of the bulk Richardson number:
Ri
b
5
g
u
a
z
u
va
2 u
vg
U
2
, (9)
where g is the gravitational acceleration, u
va
is the virtual
potential temperature of the air in the surface layer, and
u
vg
is the virtual potential temperature of the ground. To
prevent Ri
b
from being inordinately high, a lower limit
of 0.1 m s
21
is applied to U.
The first regime, Ri
b
$ 0.2, is associated with stable
(nighttime) conditions and
c
m
5 c
h
5210 ln
z
z
0
. (10)
The second one, 0 , Ri
b
, 0.2, corresponds with a dam-
ped mechanical turbulence regime wherein
c
m
5 c
h
525Ri
b
ln
z
z
0
1:1 2 5Ri
b
. (11)
The third regime, Ri
b
5 0, is associated with forced
convection:
c
m
5 c
h
5 0, (12)
and the fourth one, Ri
b
, 0, with free convection:
c
m
5 2ln
1 1 x
2
1 ln
1 1 x
2
2
2 2 tan
21
x 1
p
2
,
(13)
c
h
5 2ln
1 1 x
2
2
, (14)
wherein x 5 [1 2 16(z/L)]
1/4
and the Monin–Obukhov
stability parameter:
z
L
5 k
g
u
a
z
u
*
u
2
*
(15)
is calculated using the friction velocity, u
*
5 kU/[ln(z/
z
0
) 2 c
m
(z/L)], and the temperature scale, u
*
5 k(u
a
2
u
g
)/[ln(z/z
0
) 2 c
h
(z/L)], from the previous numerical
time step. Note u
*
is negative for unstable conditions.
The functions for the stable regime 2 come from a
slight modification of the linear relationship [f
m,h
5
25(z/L)] found in the Kansas program (Arya 1988) in
order to ensure continuity with the functions of the more
stable regime 1. A limit of 210 is used for both c
h
and c
m
in order to avoid the use of the Kansas-type functions for
very stable conditions. The unstable functions of regime
4 are also from the Kansas field experiment (Paulson
1970). A lower limit of 210 is imposed to z/L to prevent
the use of these functions for very unstable conditions.
900 MONTHLY WEATHER REVIEW VOLUME 140
The wind, temperature, and moisture are diagnosed at
their typical observational height using the integrated
dimensionless equations and assuming that u
*
, u
*
, and
q
*
are constant with height. For instance, using Eq. (4)
to obtain an expression for the wind at z 5 10 m and
dividing by the general form of the same Eq. (4) leads to
u
10m
5 u
a
ln
10
z
0
2 c
m
10
L
ln
z
z
0
2 c
m
z
L
,
where, as for the case of the transfer coefficients, the
contribution of c
m
(z
0
/L) has been neglected. An anal-
ogous derivation is used to diagnose the temperature
and moisture at 2 m:
u
2m
5 u
g
1 (u
a
2 u
g
)
ln
2
z
0
2 c
h
2
L
ln
z
z
0
2 c
h
z
L
q
2m
5 q
g
1 (q
a
2 q
g
)
ln
rc
p
ku
*
2
c
s
1
2
z
l
!
2 c
h
2
L
ln
rc
p
ku
*
z
c
s
1
z
z
l
!
2 c
h
z
L
.
Additional restrictions to the allowable values of
certain variables used to compute the fluxes and the
near-surface variables are introduced in order to avoid
undesired effects. For instance, ln(z/z
0
) 2 c
h
(z/L) is not
allowed to be lower than 2 in order to avoid a high heat
exchange coefficient, defined as the right-hand side of
Eq. (2) except for the temperature differences, during
unstable conditions in very thin surface layers with high
roughness length. For similar reasons, c
h,m
is forced to
be lower or equal to 0.9 ln(z/z
0
). In addition, the friction
velocity is arithmetically averaged with its previous value
in order to prevent large oscillations, and a lower limit of
u
*
5 0.1 m s
21
is imposed in order to prevent the heat
flux from being zero under very stable conditions. It was
considered that a smaller u
*
could potentially decouple
the temperature of the atmosphere from the ground that
starts cooling by radiation faster than observed in what is
known as the runaway cooling effect (e.g., Louis 1979).
b. Limitations of the present formulation
Aside from the problems associated with the use of the
Kansas similarity functions already mentioned, some other
limitations in the above surface layer formulation can be
pointed out. For instance, the lower limit u
*
5 0.1 m s
21
used to avoid a potential runaway cooling effect, or the
limits ln(z/z
0
) 2 c
h
(z/L) . 2andc
m,h
(z/L) # 0.9 ln(z/z
0
)
to avoid undesired effects in unstable conditions and
very thin surface layers affect the self-consistency be-
tween the surface layer variables. For the case of u
*
the
influence of the limit produces another negative im-
pact, it prevents u
*
from reproducing the observed be-
havior since u
*
values below 0.1 m s
21
are common
during the night (e.g., Shin and Hong 2011).
There are some other limitations that are not obvious
at first glance, but are responsible for inconsistencies in
the formulation. These limitations become evident in
the dispersion diagram of Ri
b
versus z/L at one location
wherein both variables were recorded every numerical
time step (Fig. 1a). Only information regarding the un-
stable part is displayed since the formulation does not
require computing z/L in the two stable regimes (see
previous section).
Theoretically, Ri
b
and z/L share the following re-
lationship (e.g., Arya 1988):
Ri
b
5
z
L
ln
z
z
0
2 c
h
z
L
ln
z
z
0
2 c
m
z
L
2
, (16)
FIG. 1. (a) The Ri
b
vs z/L from the WRF output and (b) Ri
b
vs
z/L diagnosed [Eq. (15)] with data from the WRF output. The gray
line in (b) shows the theoretical relationship.
M
ARCH 2012 J I M E
´
NEZ ET AL. 901
which indicates that no scatter should be expected in the
dispersion diagram of Fig. 1a. This is obviously not the
case since noticeable dispersion can be appreciated (inset
in Fig. 1a). The restrictions associated with the limit of
z/L 5210 are clearly evident. A tendency to report
zero z/L values can also be appreciated. More information
regarding the dispersion diagram between Ri
b
and z/L,
including the part associated with the stable regime, can
be obtained using Eq. (15) to diagnose z/L from the model
output. The dispersion diagram between Ri
b
and the di-
agnosed z/L is shown in Fig. 1b. The tendency to report
zero values of z/L has disappeared but the scatter is even
larger, especially for stable situations. There are also some
situations wherein the surface layer is unstable accord-
ing to Ri
b
(negative values), but stable according to z/L
(positive values). In addition, there is a large discrepancy
with the theoretical relationship between both variables
(see inset of Fig. 1b).
Figure 2 provides a further understanding of the
sources of the scatter. The dispersion diagram of Ri
b
and
the diagnosed z/L after removing those instances wherein
u
*
reaches its limit of 0.1 m s
21
isshowninFig.2a.A
large part of the scatter has been eliminated. This in-
dicates that the limit in u
*
was altering the values of z/L
through Eq. (15). However, inconsistencies in the stabil-
ity definition are still evident since there are instances
wherein both variables present different signs (Fig. 2a).
The origin of this inconsistency is the use of the virtual
potential temperature in the Ri
b
calculation [Eq. (9)]
and the potential temperature for u
*
in the z/L calcu-
lation [Eq. (15)]. This is evident in Fig. 2b, which re-
moves from the dispersion diagram those instances
with u
*
5 0.1 m s
21
and diagnoses z/L with the virtual
potential temperature. Inconsistencies in sign, and
therefore stability, no longer appear. This change also
reduces the scatter significantly. The inconsistencies in
the stability definition are the reason for the large number
of cases with z/L 5 0 during unstable conditions (Fig. 1a),
since the surface layer scheme reports z/L 5 0when
Ri
b
indicates a unstable surface layer (positive) and
z/L reports a stable surface layer (negative). Although
a large part of the scatter has disappeared after the
previous corrections, a noticeable dispersion in the
unstable part is still evident (Fig. 2b). The reason for
the largest dispersion that still remains is the limit ln(z/
z
0
) 2 c
h
(z/L) . 2 since the scatter is reduced after the
situations exceeding the limit are removed (Fig. 2c). A
better relationship between Ri
b
and z/L is obtained in
stable conditions if those instances with c 5210 are also
removed (Fig. 2d). This suppresses a change in the
slope in the stable part and leads to a relationship that
is in agreement with the theoretical behavior (gray line
in Fig. 2d). Some scatter around the theoretical line is
still evident, which is ultimately related to the use of
information from the previous time step to calculate
z/L [Eq. (15)].
FIG. 2. The Ri
b
vs z/L after removing the effects of (a) u
*
5 0.1 m s
21
, (b) including the
effects of moisture in the diagnosis of z/L, (c) removing the cases wherein ln(z/z
0
) 2 c
h
(z/L) 5
2, and (d) removing the instances with c
h,m
5210. The gray line in (d) shows the theoretical
relationship.
902 MONTHLY WEATHER REVIEW VOLUME 140
