of 29
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
A Generalized Mixing Length Closure for
1
Eddy-Diffusivity Mass-Flux Schemes of Turbulence and
2
Convection
3
Ignacio Lopez-Gomez
1
, Yair Cohen
1
, Jia He
1
, Anna Jaruga
1
,
2
, Tapio
4
Schneider
1
,
2
5
1
California Institute of Technology, Pasadena, California, USA.
6
2
Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, US
7
Key Points:
8
EDMF schemes represent boundary-layer (BL) turbulence and convection sepa-
9
rately, yet consistently.
10
A mixing length model based on kinetic energy constraints represents turbulent
11
fluxes in EDMF schemes well.
12
The resulting EDMF scheme captures dynamic regimes ranging from stable BLs
13
to stratocumulus-topped BLs.
14
Corresponding author: Tapio Schneider,
tapio@caltech.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Abstract
15
Because of their limited spatial resolution, numerical weather prediction and climate mod-
16
els have to rely on parameterizations to represent atmospheric turbulence and convec-
17
tion. Historically, largely independent approaches have been used to represent bound-
18
ary layer turbulence and convection, neglecting important interactions at the subgrid scale.
19
Here we build on an eddy-diffusivity mass-flux (EDMF) scheme that represents all subgrid-
20
scale mixing in a unified manner, partitioning subgrid-scale fluctuations into contribu-
21
tions from local diffusive mixing and coherent advective structures and allowing them
22
to interact within a single framework. The EDMF scheme requires closures for the in-
23
teraction between the turbulent environment and the plumes and for local mixing. A second-
24
order equation for turbulence kinetic energy (TKE) provides one ingredient for the dif-
25
fusive local mixing closure, leaving a mixing length to be parameterized. A new mixing
26
length formulation is proposed, based on constraints derived from the TKE balance. It
27
expresses local mixing in terms of the same physical processes in all regimes of bound-
28
ary layer flow. The formulation is tested at a range of resolutions and across a wide range
29
of boundary layer regimes, including a stably stratified boundary layer, a stratocumulus-
30
topped marine boundary layer, and dry convection. Comparison with large eddy sim-
31
ulations (LES) shows that the EDMF scheme with this diffusive mixing parameteriza-
32
tion accurately captures the structure of the boundary layer and clouds in all cases con-
33
sidered.
34
Plain Language Summary
35
Turbulence and convection transport heat and moisture in the atmosphere and are
36
ultimately responsible for the formation of clouds. However, they act on scales far too
37
small to be resolved in current global atmosphere models. Instead, parameterizations have
38
to be used to approximate their average effect on the finite volumes that are resolved in
39
a global model. These parameterizations are often tailored to specific atmospheric con-
40
ditions and fail when those conditions are not met. Here we propose a parameterization
41
that aims to reproduce the average effect of turbulent heat and moisture transport un-
42
der arbitrary atmospheric conditions. Numerical simulations demonstrate the accuracy
43
of the parameterization in simulating turbulence in atmospheric boundary layers under
44
stable and convective conditions, including the simulation of stratocumulus clouds.
45
1 Introduction
46
Turbulence is ubiquitous in the planetary boundary layer. Small-scale chaotic air
47
motions enhance mixing, homogenizing temperature and water content in the lower tro-
48
posphere. Under statically unstable conditions, convective updrafts and downdrafts fur-
49
ther increase the vertical transport of heat and moisture between the surface and the air
50
aloft. Together, turbulence and convection shape the vertical distribution of tempera-
51
ture and water vapor that sustains clouds. However, these processes act on scales far too
52
small to be resolved in global climate models (GCMs), with resolutions constrained by
53
current computational power (Schneider et al., 2017). Although the unabated increase
54
in processing power will make resolving deep convective processes possible in the com-
55
ing years (Kajikawa et al., 2016), resolving turbulent mixing and shallow convection will
56
remain an intractable problem for decades. Instead, parameterizations have to be used
57
to approximate the average effect of these subgrid-scale processes on the grid scale.
58
Conventional parameterizations consider atmospheric turbulence and convection
59
as independent processes, neglecting interactions that alter their combined effect on the
60
large scale. These parameterizations are often regime-dependent, leading to models that
61
artificially split the spectrum of atmospheric conditions into a discrete number of cases.
62
Examples of such case-dependent approaches include parameterizations of cumulus (Arakawa,
63
2004) and stratocumulus clouds (Lilly, 1968; Schubert, 1976). However accurate, the use
64
–2–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
of disparate schemes for different conditions complicates a seamless representation of subgrid-
65
scale processes in the lower troposphere.
66
Several approaches to obtain a unified model of turbulence and convection have been
67
proposed (Lappen & Randall, 2001; Park, 2014; Thuburn et al., 2018). Here we focus
68
on the extended formulation of an eddy-diffusivity mass-flux (EDMF) scheme developed
69
in Tan et al. (2018), which in turn built on work by Siebesma and Teixeira (2000); Soares
70
et al. (2004); Siebesma et al. (2007) and Angevine et al. (2010), among others. In the
71
EDMF framework, the flow within each grid cell is decomposed into several distinct sub-
72
domains, representing coherent convective structures and their relatively isotropic tur-
73
bulent environment. Convective transport is captured by mass flux terms that depend
74
on differences between subdomain-mean properties; more isotropic turbulent transport,
75
associated with small-scale fluctuations within each subdomain, is captured by eddy dif-
76
fusion closures.
77
The extended EDMF framework uses additional prognostic equations for subdo-
78
main variables, such as the environmental turbulence kinetic energy, and it requires clo-
79
sures for local turbulent fluxes and for the mass exchange between subdomains (Tan et
80
al., 2018). Even though the EDMF framework arises from the need for a unified model
81
of turbulence and convection, the parameterizations used for entrainment and turbulent
82
mixing are usually defined differently for each regime (Suselj et al., 2013; Witek et al.,
83
2011). The development of regime-independent parameterizations for the required clo-
84
sures is the last step in the construction of a unified model of atmospheric turbulence
85
and convection.
86
Here, a regime-independent closure for turbulent mixing within the EDMF frame-
87
work is proposed. Section 2 reviews the decomposition of subgrid-scale fluxes in the ex-
88
tended EDMF scheme. Section 3 introduces the formulation of the closure. Section 4 il-
89
lustrates the performance of the EDMF scheme with the turbulent mixing closure in bound-
90
ary layer regimes where vertical transport is strongly dependent on the turbulence clo-
91
sure used: the stable boundary layer (SBL), the stratocumulus-topped boundary layer
92
(STBL), and dry convection. The performance of the extended EDMF scheme with this
93
closure in moist-convective cases is demonstrated in a companion paper (Cohen et al.,
94
2020). Finally, Section 5 summarizes the results and conclusions.
95
2 EDMF Framework
96
In the EDMF framework, each grid-cell volume is decomposed into
n
updrafts or
downdrafts (labeled by index
i
= 1
,...,n
) and an environment (labeled by index
i
=
0) in which they are embedded. Following this decomposition, the grid-mean value of
variable
ψ
may be written as
ψ
=
n
i
0
a
i
̄
ψ
i
.
(1)
Here, angle brackets
〈·〉
denote the grid mean,
̄
ψ
i
denotes the Favre average of
ψ
over
subdomain
i
, and
a
i
is the mean horizontal cross-sectional area covered by subdomain
i
within the grid cell. This partition is motivated by the anisotropy of turbulent convec-
tive flows, in which isotropic turbulent eddies coexist with coherent columnar structures
that induce a strong vertical transport (Bjerknes, 1938). The subdomain decomposition
is simplified for the horizontal velocity
u
h
, which is taken to have the same mean value
for all subdomains,
̄
u
h,i
=
u
h
. Applying the subdomain decomposition to higher-order
moments introduces additional terms associated with the difference between grid and
subdomain means. For the vertical subgrid-scale flux of
ψ
, this leads to
w
ψ
=
n
i
0
a
i
(
w
i
ψ
i
+ ̄
w
i
̄
ψ
i
)
.
(2)
–3–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Here,
w
is the vertical velocity,
ψ
=
ψ
−〈
ψ
,
ψ
i
=
ψ
̄
ψ
i
, and
̄
ψ
i
=
̄
ψ
i
−〈
ψ
. The de-
97
composition (2) partitions the subgrid-scale flux into contributions from small-scale fluc-
98
tuations, associated with turbulence, and subdomain-mean terms, representative of con-
99
vection. In the following, we will refer to these contributions as turbulent and convec-
100
tive fluxes, respectively.
101
The subdomain-mean terms can be explicitly solved for by introducing
n
prognos-
tic subdomain equations for each variable and an additional equation for each plume area
fraction
a
i
, which may be diagnostic or prognostic (Tan et al., 2018). Turbulent fluxes
within each subdomain are modeled as downgradient and proportional to an eddy dif-
fusivity
K
ψ,i
, where
ψ
is the property being transported. For the vertical turbulent flux
in (2), this gives
w
i
ψ
i
=
K
ψ,i
̄
ψ
i
∂z
.
(3)
The eddy diffusivity
K
ψ,i
is proportional to a characteristic velocity scale and the length
102
scale of the eddies driving the transport, both of which must be parameterized.
103
Proposed closures for the eddy diffusivity vary from simple diagnostic expressions
to second-order models that introduce prognostic equations for both scales (Umlauf &
Burchard, 2003). The 1.5-order turbulence kinetic energy (TKE) model
1
is a particu-
larly popular choice due to its balance between accuracy and computational efficiency
(Mellor & Yamada, 1982). The 1.5-order model makes use of a prognostic equation for
TKE and a diagnostic expression for the mixing length. In the EDMF framework, the
grid-mean TKE
e
can be decomposed following expression (2) for second-order moments
as
e
=
n
i
0
a
i
(
̄
e
i
+
̄
w
i
̄
w
i
2
)
,
(4)
where ̄
e
i
is the TKE of subdomain
i
. This expression can be simplified by assuming that
for the updrafts and downdrafts (
i >
0), the contribution to the grid-mean TKE from
small-scale turbulence is negligible compared to the convective term, an assumption com-
monly made in EDMF schemes:
e
=
a
0
̄
e
0
+
n
i
0
a
i
̄
w
i
̄
w
i
2
.
(5)
The TKE decomposition (5) can also be obtained by assuming a small updraft and down-
104
draft area fraction and similar turbulence intensity in all subdomains (Siebesma et al.,
105
2007). However, the equations derived for the subdomain second-order moments with
106
these two approaches differ in the source terms that appear due to entrainment processes
107
between subdomains. The former approximation is favored here to allow for the use of
108
this framework in high-resolution models, where the assumption of slender updrafts may
109
become inadequate (Randall, 2013).
110
Given an updraft area fraction
a
i
, the grid-mean TKE is determined by the envi-
ronmental TKE ̄
e
0
and the subdomain-mean vertical velocities ̄
w
i
. The subdomain-mean
vertical velocity equation for subdomain
i
is
(
ρa
i
̄
w
i
)
∂t
+
(
ρa
i
̄
w
2
i
)
∂z
+
h
·
(
ρa
i
̄
u
h,i
̄
w
i
) =
(
ρa
i
w
i
w
i
)
∂z
−∇
h
·
(
ρa
i
u
h,i
w
i
)
+
j
6
=
i
[
E
ij
̄
w
j
ij
̄
w
i
+
ˆ
E
ij
( ̄
w
j
̄
w
i
)
]
+
ρa
i
̄
b
i
ρa
i
̄
Ψ
i
∂z
,
(6)
1
This model is also referred to as the Level 2.5 model in the Mellor-Yamada hierarchy (Mellor & Ya-
mada, 1982).
–4–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
where
h
is the horizontal gradient operator, Ψ =
p/ρ
is the pressure potential and
the turbulent transport terms on the right-hand side are negligible for all subdomains
except the environment (
i
= 0). Subgrid density changes are only considered in the buoy-
ancy term, such that
ρ
=
ρ
in the previous equation, in order to avoid creation of spu-
rious acoustic modes through the subdomain decomposition (Cohen et al., 2020). The
buoyancy
̄
b
i
and the pressure potential anomaly
̄
Ψ
i
are defined with respect to a refer-
ence hydrostatic pressure profile
p
h
(
z
) and density
ρ
h
(
z
), related by
z
p
h
=
ρ
h
g
:
̄
b
i
=
g
̄
ρ
i
ρ
h
ρ
,
̄
Ψ
i
∂z
=
∂z
(
̄
p
i
ρ
)
+
g
ρ
h
ρ
.
(7)
Here, ̄
p
i
is the subdomain-mean pressure. Density appears inside the pressure gradients
111
in (6) and (7) to ensure thermodynamic consistency of the subgrid-scale anelastic ap-
112
proximation (Cohen et al., 2020). Interactions between subdomains are captured by en-
113
trainment and detrainment fluxes. In the vertical velocity equation (6), ∆
ij
is the dy-
114
namical detrainment of air mass from subdomain
i
into subdomain
j
, and
E
ij
and
ˆ
E
ij
115
are the dynamical and turbulent entrainment from subdomain
j
into subdomain
i
, re-
116
spectively. It is assumed that entrainment events occur over timescales much shorter than
117
the eddy turnover rate
K
ψ,i
/
̄
e
i
, so that entrained air carries the properties of the sub-
118
domain it detrains from. In addition, for now we assume entrainment occurs only be-
119
tween convective plumes and the environment, not among plumes.
120
The prognostic equation for environmental TKE can be written in non-conservative
form as
̄
e
0
∂t
+ ̄
w
0
̄
e
0
∂z
+
u
h
〉·∇
h
̄
e
0
=
w
0
u
0
u
∂z
w
0
v
0
v
∂z
w
2
0
̄
w
0
∂z
+
w
0
b
0
−P
1
ρa
0
∂z
(
ρa
0
w
0
e
0
)
+
i>
0
[
i
0
ρa
0
(
( ̄
w
i
̄
w
0
)
2
2
̄
e
0
)
ˆ
E
i
0
ρa
0
( ̄
w
0
( ̄
w
i
̄
w
0
) + ̄
e
0
)
]
−D
1
ρa
0
h
·
(
ρa
0
u
h,
0
e
0
)
u
h,
0
u
0
·∇
h
u
〉−
u
h,
0
v
0
·∇
h
v
〉−
u
h,
0
w
0
·∇
h
̄
w
0
.
(8)
Here,
u
and
v
are the components of
u
h
,
P
is the velocity pressure-gradient cor-
relation, and
D
is the turbulent dissipation. All sources and sinks of ̄
e
0
account for un-
resolved processes on the grid scale, so they must be parameterized. Subdomain covari-
ances in (8) are modeled diffusively, with the environmental eddy diffusivity
K
ψ
defined
as
K
ψ
=
c
ψ
l
̄
e
1
/
2
0
,
(9)
where
l
is the mixing length, and
c
ψ
is a fitting parameter. The subscript 0 in the eddy
121
diffusivity is dropped to simplify notation. The coefficient
c
ψ
is taken to be equal to
c
h
122
for the diffusion of all fields except for momentum, for which
c
ψ
=
c
m
. The eddy vis-
123
cosity
K
m
is related to
K
h
through the turbulent Prandtl number Pr
t
, such that
K
m
=
124
Pr
t
K
h
.
125
Under the assumption that subgrid-scale pressure work on the grid-mean is neg-
126
ligible,
P
is taken as opposite to the pressure work on the plumes (Tan et al., 2018). Hence,
127
P
acts as a return-to-isotropy term on the full grid, transferring momentum from the
128
strongly anisotropic coherent structures into the relatively isotropic eddies in the envi-
129
ronment:
130
P
=
[
w
0
(
Ψ
∂z
)
0
+
u
0
(
Ψ
∂x
)
0
+
v
0
(
Ψ
∂y
)
0
]
=
i>
0
a
i
a
0
( ̄
w
i
̄
w
0
)
̄
Ψ
i
∂z
,
(10)
The pressure work on the plumes is formulated in terms of contributions from a virtual
131
mass term (Gregory, 2001), an advective term (Jia He, personal communication) and a
132
–5–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
drag term (Romps & Charn, 2015), yielding the following expression for the velocity pressure-
133
gradient correlation:
134
P
=
i>
0
a
i
a
0
( ̄
w
i
̄
w
0
)
(
α
b
̄
b
i
α
a
̄
w
i
̄
w
i
∂z
+
α
d
( ̄
w
i
̄
w
0
)
|
̄
w
i
̄
w
0
|
H
i
)
,
(11)
where
α
a
and
α
d
are constant parameters,
H
i
is the plume height and
α
b
is a function
of the aspect ratio of the plume (Jia He, personal communication). Finally, assuming
statistical equilibrium at scales
l
(Vassilicos, 2015), turbulent dissipation can be estimated
from the spectral transport relation that follows from Kolmogorov’s theory of inertial
turbulence, giving Taylor’s dissipation surrogate
D
=
c
d
̄
e
3
/
2
0
l
.
(12)
Here,
c
d
is an empirical coefficient and
l
is the dissipation length, taken to be equal to
135
the mixing length in our model. Expressions (3) and (5)–(12) provide closure to a 1.5-
136
order model of turbulence within the EDMF framework, given diagnostic expressions for
137
the mixing length
l
and for entrainment and detrainment.
138
3 Mixing Length Formulation
139
We seek to obtain a regime-independent eddy diffusivity closure that provides an
140
accurate representation of turbulent subgrid-scale fluxes, over a wide range of host model
141
resolutions. Thus, the eddy diffusivity should reduce to an LES-type closure at high res-
142
olution, while being able to account for the processes that modify turbulent fluxes at larger
143
scales. The formulation of the closure is organized following this logic.
144
3.1 Minimum Dissipation of Environmental TKE
145
As in Verstappen (2011) and Abkar and Moin (2017), we assume that at the small
scales represented by the environment in the EDMF scheme, TKE is dissipated at least
at the rate at which it is produced. This condition translates into an inequality for the
production and dissipation terms in the environmental TKE budget:
w
0
b
0
w
0
u
0
u
∂z
w
0
v
0
v
∂z
w
2
0
̄
w
0
∂z
u
h,
0
u
0
·∇
h
u
〉−
u
h,
0
v
0
·∇
h
v
u
h,
0
w
0
·∇
h
̄
w
0
+
i>
0
[
i
0
ρa
0
(
( ̄
w
i
̄
w
0
)
2
2
̄
e
0
)
ˆ
E
i
0
ρa
0
( ̄
w
0
( ̄
w
i
̄
w
0
) + ̄
e
0
)
]
≤D
.
(13)
Here, the terms involving TKE injection from entrained air are also taken to be locally
balanced by dissipation, consistent with the assumption that entrainment events occur
over timescales much smaller than the eddy turnover time
K
ψ,i
/
̄
e
i
. The inequality (13)
is a local condition for the environment, so it does not preclude net subgrid-scale energy
production due to processes such as convection represented by plumes. The evolution
of the grid-mean TKE that follows from (5) and (13) is
e
∂t
+
u
h
〉·∇
h
e
+
w
e
∂z
+
w
e
∂z
=
i
a
i
(
̄
w
i
̄
b
i
̄
w
2
i
w
∂z
)
a
0
γ
0
+
a
0
(
w
0
w
0
̄
w
0
∂z
+
u
h,
0
w
0
·∇
h
̄
w
0
)
i>
0
[
i
0
ρ
(
( ̄
w
i
̄
w
0
)
2
2
)
ˆ
E
i
0
ρ
̄
w
0
( ̄
w
i
̄
w
0
)
]
,
(14)
where
γ
0
is the net environmental TKE dissipation with which the TKE production–dissipation
146
inequality (13) becomes an equality. Under the net dissipation closure, grid-mean TKE
147
production occurs through the first two terms on the right-hand side of (14): the con-
148
vective buoyancy flux and the subdomain-scale shear production. The subgrid-scale ki-
149
netic energy pathways in the extended EDMF scheme are described in Appendix B.
150
–6–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Using Taylor’s dissipation surrogate (12) and downgradient closures for the shear
and buoyancy terms, the net dissipation condition (13) leads to a condition for the max-
imum value of the mixing length
l
at which the net dissipation
γ
0
is still positive semidef-
inite:
{
3
k
=1
[
(
u
∂x
k
)
2
+
(
v
∂x
k
)
2
+
(
̄
w
0
∂x
k
)
2
]
1
Pr
t
̄
b
0
∂z
}
l
2
+
i>
0
[
i
0
ρa
0
(
( ̄
w
i
̄
w
0
)
2
2
̄
e
0
)
ˆ
E
i
0
ρa
0
( ̄
w
0
( ̄
w
i
̄
w
0
) + ̄
e
0
)
]
l
c
d
c
m
̄
e
0
.
(15)
Here, the environmental buoyancy gradient is computed following Tan et al. (2018). From
this inequality, an expression for the mixing length that minimizes turbulent dissipation
can be obtained by solving for
l
. For the resulting value of the mixing length, produc-
tion and dissipation of TKE are locally balanced:
l
tke
=
−I
2(
S
l
+
B
l
)
=
I
2(
S
l
+
B
l
)
+
I
2
+ 4(
S
+
B
)
D
2(
S
l
+
B
l
)
.
(16)
Here, ∆ is the discriminant and the different terms are given by
S
l
+
B
l
=
c
m
̄
e
1
/
2
0
{
3
k
=1
[
(
u
∂x
k
)
2
+
(
v
∂x
k
)
2
+
(
̄
w
0
∂x
k
)
2
]
1
Pr
t
̄
b
0
∂z
}
,
I
=
i>
0
[
i
0
ρa
0
(
( ̄
w
i
̄
w
0
)
2
2
̄
e
0
)
ˆ
E
i
0
ρa
0
( ̄
w
0
( ̄
w
i
̄
w
0
) + ̄
e
0
)
]
,
S
+
B
= (
S
l
+
B
l
)
l.
(17)
In (16), the product (
S
+
B
)
D
is independent of the mixing length, so
l
tke
can be read-
ily evaluated. Although the term (
S
+
B
) is sign-indefinite, the discriminant
∆ =
I
2
+ 4(
S
+
B
)
D
in (16) can be shown to remain positive semidefinite even when the shear and buoyancy
terms result in TKE destruction, provided that the inequality (13) holds. This is because
the minimum dissipation balance requires
I
=
D−
(
S
+
B
)
,
(18)
so that the expression for the discriminant ∆ is of the form
∆ = [
D−
(
S
+
B
)]
2
+ 4(
S
+
B
)
D
= [
D
+ (
S
+
B
)]
2
0
.
(19)
The mixing length
l
tke
depends on local characteristics of the environment and on
151
the vertical velocity difference between subdomains, which enter the entrainment and
152
detrainment terms. Thus, convection modifies the environmental diffusive transport through
153
entrainment processes. This approach can also be applied to turbulence models that re-
154
tain covariance terms
w
i
ψ
i
for other subdomains, and not only for the environment. In
155
this case, the minimum dissipation condition may be used to obtain a characteristic mix-
156
ing length
l
tke
,i
for each subdomain. However, variance within plumes can also be ac-
157
counted for by variance among plumes when the number of subdomains is increased.
158
In stably stratified boundary layers, where convection is inhibited, pressure work
and entrainment fluxes in (6) act to homogenize the different subdomains, such that
̄
ψ
i
0 for any variable
ψ
and
a
0
1 (i.e., there are no convective plumes). Under these con-
ditions, the minimum dissipation mixing length (16) reduces to the expression proposed
–7–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
by Grisogono (2010) for steady-state stable boundary layer (SBL) flow:
l
tke
=
(
S
+
B
)
D
(
S
l
+
B
l
)
=
c
d
c
m
e
{
3
k
=1
[
(
u
∂x
k
)
2
+
(
v
∂x
k
)
2
+
(
w
∂x
k
)
2
]
1
Pr
t
b
∂z
}
1
/
2
.
(20)
The balance between shear production, destruction due to stratification, and dissipation,
159
which arises when using this mixing length, is a well-known leading-order state in neu-
160
tral (Spalart, 1988) and moderately stable boundary layer flows (Li et al., 2016).
161
3.2 Limitations of the Minimum-Dissipation Closure
162
Expression (16) for the mixing length
l
tke
captures the leading-order balance in the
environmental TKE budget at small scales. However, a model with a diffusive closure
based on
l
tke
cannot fully describe the dynamics of the boundary layer at the coarse res-
olutions typical of GCMs, on the order of 10
4
m in the horizontal and 10
100 m in the
vertical. At these scales, the resolved horizontal gradients are weak, and the environmen-
tal TKE equation (8) can be simplified using the boundary layer approximation (neglect-
ing horizontal relative to vertical derivatives):
̄
e
0
∂t
+ ̄
w
0
̄
e
0
∂z
=
1
ρa
0
∂z
(
ρa
0
w
0
e
0
)
(
w
0
u
0
u
∂z
+
w
0
v
0
v
∂z
+
w
2
0
̄
w
0
∂z
)
+
w
0
b
0
+
ρ
i>
0
[
i
0
ρa
0
(
( ̄
w
i
̄
w
0
)
2
2
̄
e
0
)
ˆ
E
i
0
ρa
0
( ̄
w
0
( ̄
w
i
̄
w
0
) + ̄
e
0
)
]
−P −D
.
(21)
In stable conditions, using
l
tke
for the mixing length and integrating the conservative form
of (21) from the surface layer (
z
s
) to the free troposphere above (
z
i
), the evolution equa-
tion for the vertically integrated environmental TKE reduces to
z
i
z
s
(
ρa
0
̄
e
0
)
∂t
dz
=
[
ρa
0
w
0
e
0
]
z
i
z
s
≈−
ρa
0
K
m
∂e
0
∂z
z
s
.
(22)
Note that in stable conditions,
a
0
1 and
̄
ψ
i
0 for any variable
ψ
. From (22), it
163
follows that the evolution of the vertically integrated TKE under the minimum dissipa-
164
tion closure only depends on the flux from the unresolved surface layer. But unbalanced
165
TKE dissipation has been observed to become increasingly important as stratification
166
develops in field studies of the atmospheric boundary layer (Li et al., 2016), and it can
167
be expected to play a role in conditions of strong surface cooling. The budget (22) can-
168
not capture unbalanced TKE destruction within the boundary layer due to stratifica-
169
tion. Furthermore, the minimum dissipation mixing length
l
tke
leads to enhanced eddy
170
diffusion with increasing stratification, contrary to the evidence of turbulent mixing be-
171
ing inhibited in strong stratification, such as near strong inversions.
172
The limitations of a minimum dissipation model also become apparent in convec-
tively unstable boundary layers. The use of expression (16) for the mixing length results
in a simplified form of the TKE balance (21) because of the strict balance of all produc-
tion and dissipation terms in (13). Integrating the TKE equation in the vertical, the evo-
lution of the vertically integrated environmental TKE in convective conditions reads
z
i
z
s
(
ρa
0
̄
e
0
)
∂t
dz
=
ρa
0
w
0
e
0
|
z
s
+
z
i
z
s
ρ
i>
0
a
i
( ̄
w
i
̄
w
0
)
̄
Ψ
i
∂z
dz
+
z
i
z
s
i>
0
(∆
i
0
E
i
0
) ̄
e
0
dz.
(23)
–8–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Here, the last term only accounts for changes in environmental area fraction and does
173
not result in a source or sink of ̄
e
0
. A major difference between the SBL budget (22) and
174
the convective budget (23) is the contribution of the velocity pressure-gradient correla-
175
tion. From the velocity pressure-gradient relation (11), pressure work captures the im-
176
portant energization of turbulence in the environment owing to ascending or descend-
177
ing plumes (Schumann & Moeng, 1991). At the grid-scale, the source of this subgrid-
178
scale energy term is the convective buoyancy flux in (14), which accelerates the plumes
179
in convective conditions.
180
The TKE balance (23) shows that, in convective conditions, the source of environ-
181
mental TKE from updrafts or downdrafts can only be compensated by the flux from the
182
unresolved surface layer. This is often a source term rather than a sink term, because
183
shear production is surface intensified. Thus, the TKE balance (23) suggests an unbal-
184
anced growth of TKE in convective boundary layers. This continuous TKE increase in
185
convective conditions is inconsistent with LES results showing quasi-stationary TKE lev-
186
els in convective boundary layers (Nieuwstadt et al., 1993).
187
The TKE balances (22) and (23) highlight the shortcomings of the minimum dis-
188
sipation balance (16) as a general closure for diffusive mixing in the boundary layer. The
189
lack of net dissipation mechanisms in the vertically integrated TKE balance hinders the
190
correct representation of important processes, such as the shallowing of the boundary
191
layer in the late afternoon or the sharp mixing inhibition near inversions. Moreover, it
192
precludes reaching a quasi-stationary state in statically unstable boundary layers. Nev-
193
ertheless, the limitations of the minimum dissipation model can be used to inform the
194
construction of a generalized master length scale based on the TKE production–dissipation
195
inequality (13).
196
The limitations of the minimum dissipation balance showcased in this section are
197
not necessarily applicable to other turbulence models. For example, He et al. (2019) use
198
the production-dissipation condition to diagnose TKE and eddy diffusivity from a mix-
199
ing length
l
. This allows an instantaneous adjustment of TKE to a new balanced state,
200
at the cost of representing convection with an empirical parameterization that has no
201
subgrid interaction with turbulent diffusion.
202
3.3 Constrained Minimization of TKE Dissipation
203
A master length scale that corrects the limitations of the minimum-dissipation model
can be constructed by taking dissipation to be higher than production under certain cir-
cumstances. Using closures of the form (12) for the dissipation and (9) for turbulent dif-
fusion, it follows from the production–dissipation inequality (13) that excess dissipation
occurs for
l < l
tke
. Hence, unbalanced TKE dissipation arises naturally in regions of
the boundary layer where the characteristic size of environmental eddies is constrained
to be smaller than
l
tke
. A general mixing length capturing this condition can be writ-
ten as
l
=
s
min
(
l
tke
,l
1
,l
2
,...
)
,
(24)
where
l
j
(
j
= 1
,
2
,...,N
) are candidate mixing lengths based on flow constraints, and
s
min
(
x
) is a smooth minimum function defined in Appendix A. The TKE production–
dissipation inequality (15) with the closures substituted implies that the minimum length
scale provides maximum TKE dissipation. Thus, the use of the minimum length scale
(24) is equivalent to the minimization of TKE dissipation in (13) subject to the constraint
that dissipation exceeds the candidate dissipation rates,
D ≥ D|
l
=
l
j
j,
(25)
where
D|
l
=
l
j
is the candidate dissipation rate evaluated at length scale
l
j
.
204
Our suggestion for choosing a general mixing length as a smooth minimum of var-
ious candidates contrasts with the common practice (e.g., He et al., 2019; Han & Brether-
–9–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
ton, 2019) to use the expression suggested by Blackadar (1962),
l
h
=
(
1
l
1
+
1
l
2
)
1
,
(26)
for a master length scale
l
h
. This length scale
l
h
, proportional to the harmonic mean of
205
the candidates
l
1
and
l
2
, is smaller than both
l
1
and
l
2
. If closures similar to (9) and (12)
206
are used in a prognostic equation for TKE, the mixing length (26) results in an unre-
207
alistic intensification of TKE dissipation in regions where the candidate length scales
l
1
208
and
l
2
are similar. This undesirable characteristic is avoided by using the smooth min-
209
imum (24).
210
We consider two limiting factors for the characteristic length scale of turbulent mo-
211
tion in boundary layer flows: stable stratification and the distance to solid boundaries.
212
3.3.1 Stratification Constraints
213
Environmental stratification constrains the size of turbulent eddies by inhibiting
the vertical displacement of air masses. Stably stratified turbulence is known to show
high vertical variability and reorganization into layered structures, with most mixing oc-
curring within the layers (Waite, 2011). The thickness of these layers is determined by
the vertical scale at which the governing dynamic equations become self-similar (Billant
& Chomaz, 2001; Augier et al., 2012), known as the buoyancy scale
l
b
. For a flow with
an imposed stratification given by the Brunt-V ̈ais ̈al ̈a frequency
N
e
, this length scale is
l
b
=
c
b
( ̄
e
0
)
1
/
2
N
e
,
(27)
where
c
b
is an empirical coefficient. It is important to note that imposing
l
b
as an up-
per bound for the size of eddies is similar to doing so by the Ozmidov scale
l
o
D
/N
3
e
only if turbulent motions at the scale in question are assumed to be in the inertial sub-
range, such that (12) holds. In this case, an expression equivalent to (27) for the Ozmi-
dov scale is
l
o
=
(
c
3
b
c
d
D
N
3
e
)
1
/
2
.
(28)
However, recent experimental studies suggest that under strong stratification, turbulence
214
may not display an inertial subrange (Grachev et al., 2013). In that case, expression (12)
215
and the Ozmidov scale (28) may not be applicable (Li et al., 2016), whereas the buoy-
216
ancy scale (27) still holds.
217
The buoyancy frequency of moist air depends on the latent heat release and evap-
orative cooling associated with the vertical displacement of air parcels. In general, the
effective static stability
N
e
lies between the dry and the moist adiabatic limits. Follow-
ing O’Gorman (2010), we use an effective static stability of the form
N
2
e
=
g
̄
θ
v,
0
(
̄
θ
v,
0
∂z
λ
̄
θ
v,
0
∂z
θ
vl,
0
)
=
g
̄
θ
v,
0
[
(1
λ
)
̄
θ
v,
0
∂z
+
λ
̄
θ
v,
0
̄
θ
vl,
0
̄
θ
vl,
0
∂z
]
,
(29)
where
θ
v
is the virtual potential temperature and
λ
represents the area fraction of en-
vironmental air undergoing phase change. In the non-precipitating cases considered here,
λ
is given by the environmental cloud fraction
f
c,
0
. Cloud fraction diagnosis is cloud-
type dependent in many current GCMs (Collins et al., 2004). In our EDMF scheme, we
use a regime-independent probabilistic cloud scheme (see Appendix C). The liquid-water
virtual potential temperature
θ
vl
appearing in the effective static stability measures the
buoyancy of cloudy air parcels when moist-adiabatically returned to clear conditions,
̄
θ
vl
(1 +
η
̄
q
t
)
̄
θ
l
̄
θ
v
exp
(
L
v
̄
q
l
c
p
̄
T
)
.
(30)
–10–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Here,
η
=
R
v
/R
d
1,
L
v
is the latent heat of vaporization,
c
p
is the specific heat of
218
air,
q
t
and
q
l
are the total and liquid water specific humidities,
θ
l
is the liquid water po-
219
tential temperature,
T
is the temperature and
R
v
,R
d
are the gas constants for water va-
220
por and dry air, respectively. Note that the effective static stability (29) converges to
221
the dry limit when
q
l
0 for all values of
λ
; it reduces to
N
2
e
= (1
λ
)
N
2
, with dry
222
buoyancy frequency
N
, in conditions that are well mixed in
θ
l
and
q
t
.
223
3.3.2 Wall Constraints
224
The presence of boundaries also imposes an upper limit on the size of eddies near
them. Following Monin and Obukhov (1954), the eddy diffusivity in the surface layer
has the form
K
ψ,w
=
u
κz
φ
ψ
(
ξ
)
(31)
where
ξ
=
z/L
,
φ
ψ
(
ξ
) is an empirical stability function,
κ
is the von K ́arm ́an constant,
L
is the Obukhov length, and
u
is the friction velocity. The upper bound for the mix-
ing length near the surface is obtained by matching this eddy diffusivity with the expres-
sion (9) for the eddy diffusivity:
l
w
=
κ
c
ψ
κ
φ
ψ
(
ξ
)
z.
(32)
Here,
κ
= ̄
e
1
/
2
0
/u
is the ratio of rms turbulent velocity to the friction velocity in the
225
surface layer. The friction velocity in our model is diagnosed using the flux-profile re-
226
lationships of Byun (1990), except in free convective conditions. When the conditions
227
for free convection are satisfied, the diagnostic of
u
, which is a function of the horizon-
228
tal wind speed at the lowest model level, is modified following Beljaars (1995).
229
The choice of a common master length for momentum and tracer diffusion implies
c
h
φ
h
=
c
m
φ
m
, such that
φ
h
= Pr
t
φ
m
. In our formulation, the turbulent Prandtl num-
ber is taken to be a function of the gradient Richardson number Ri, based on a simpli-
fied cospectral budget of momentum and heat transport (Katul et al., 2013; Li, 2019):
Pr
t
=
2Ri
1 +
ω
2
Ri
4Ri + (1 +
ω
2
Ri)
2
Pr
t,
0
.
(33)
Here,
ω
2
= 40
/
13 is a phenomenological constant, and Pr
t,
0
is the turbulent Prandtl
number in neutral conditions. The stability function
φ
m
is often written in the form (Businger
et al., 1971; Nakanishi, 2001)
φ
m
= [1 +
a
1
(
ξ
)
ξ
]
a
2
(
ξ
)
, a
i
=
a
i
+ (
a
+
i
a
i
)
H
(
ξ
)
,
(34)
where
H
(
·
) is the Heaviside function and
a
i
,a
+
i
are empirical functions. The values of
230
a
i
are taken as negative definite to reflect the convective elongation of eddies in unsta-
231
ble conditions. In stable conditions, self-similarity of the flow requires that
a
+
2
= 1 and
232
a
+
1
>
0, such that under strong stratification, the mixing length (32) becomes indepen-
233
dent of
ξ
. As shown by Monin and Obukhov (1954), the asymptotic limit of
φ
m
under
234
strong stratification also requires that
a
+
1
= Pr
t
(Ri
st
)
/
Ri
st
. Here, Ri
st
is the asymp-
235
totic Richardson number at
ξ

1
/a
+
1
in the surface layer.
236
The empirical function (34) has been shown to become increasingly inaccurate with
stability for
ξ >
0
.
5 (Sorbjan & Grachev, 2010; Optis et al., 2016). Moreover, extend-
ing the use of the limiting scale
l
w
above the surface layer precludes the use of
a
+
1
6
= 0
in stable conditions, since the Obukhov length characterizes stratification only in the con-
stant flux layer near the surface. Although the use of
l
w
in expression (24) mandates
a
+
1
=
0, the effect of stability in eddy diffusion is still captured by
l
b
. In the constant flux layer,
the limiting length
l
b
is equivalent to the use of the empirical function (34) in the strongly
stable limit, with
a
+
1
=
1
(
κ
2
c
m
c
b
)
2
Pr
t
, ξ

1
a
+
1
.
(35)
–11–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Under weaker stratification, turbulence in the surface layer can reach a quasi-steady state
(Spalart, 1988). In this case, the limiting scale
l
w
should converge to
l
tke
. Assuming that
entrainment processes are limited to dynamical entrainment by the plumes in the sur-
face layer, the ratio of the two length scales can be written as
l
w
l
tke
̄
e
0
=
(1
Ri
/
Pr
t
)
1
/
2
(
c
d
c
m
)
1
/
2
κ
2
,
(36)
which is constant under neutral stratification and is slowly varying with Ri due to the
237
opposing effect of the Prandtl number (33). From (36), the convergence of
l
tke
to
l
w
in
238
the surface layer is satisfied for (
c
d
c
m
)
1
/
2
κ
2
1.
239
The use of a soft minimum function for the mixing length (24) allows for a smooth
240
transition from Monin-Obukhov similarity theory near the surface to a local turbulent
241
closure farther away from it, where the use of Monin-Obukhov scaling may be inaccu-
242
rate (Optis et al., 2016). In addition, the expressions (35) and (36) show that this tran-
243
sition is asymptotically consistent.
244
3.3.3 Master Mixing Length
245
Finally, the smooth minimum of the three candidate length scales (16), (27), and
(32) determines the mixing length,
l
=
s
min
(
l
tke
,l
w
,l
b
)
.
(37)
The mixing length (37) depends on a group of nondimensional parameters
C
that
must be obtained empirically:
C
=
{
c
m
,c
d
,c
b
,κ,κ
,a
1
,a
2
,
Pr
t,
0
}
.
(38)
Values for these parameters are reported in studies of boundary layer turbulence, obtained
246
from field observations (Businger et al., 1971) or LES (Nakanishi, 2001). However, the
247
direct use of some of these values in the EDMF scheme is not justified due to the decom-
248
position of the subgrid-scale flow into different subdomains. Because of the large size of
249
the parameter space
C
and the presence of other parameters in the EDMF scheme, we
250
limit the parameter optimization process to
C
=
{
c
m
,c
d
,c
b
}
in this study.
C
contains
251
the parameters that appear in the closures that are most strongly affected by the do-
252
main decomposition. All other parameters in
C
, except Pr
t,
0
, arise from similarity the-
253
ory arguments for the unresolved surface layer. Here, it is assumed that similarity ar-
254
guments are valid outside convective updrafts, and all values are taken from Nakanishi
255
(2001). For the simulations reported in the next section, the parameter space used is shown
256
in Table 1. The rest of parameters used in the EDMF scheme, which do not appear ex-
257
plicitly in the formulation of the mixing length closure, are reported in Cohen et al. (2020).
258
4 Results for Single-Column Simulations
259
Here we focus on case studies targeting the simulation of the Arctic stable bound-
260
ary layer (SBL), stratocumulus clouds, and dry convection. The performance of the ex-
261
tended EDMF scheme in moist-convective cases is explored in Cohen et al. (2020), us-
262
ing the same set of parameters. The extended EDMF scheme is tested for horizontal res-
263
olutions typical of GCMs. Invoking the boundary layer approximation (neglecting hor-
264
izontal derivatives), we perform simulations in a single-column model (SCM). The SCM
265
is a one-dimensional vertical model that aims to represent a single atmospheric column
266
within a GCM. Results from single-column simulations using the extended EDMF scheme
267
are then compared to horizontal averages obtained from LES over the same domain. LES
268
are set up by further discretizing the atmospheric column horizontally and using hor-
269
izontal doubly-periodic boundary conditions.
270
–12–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.
manuscript submitted to
Journal of Advances in Modeling Earth Systems (JAMES)
Table 1.
Parameters in the mixing length closure and values used in this study.
Symbol
Description
Value
c
m
Eddy viscosity coefficient
0
.
14
c
d
Turbulent dissipation coefficient
0
.
22
c
b
Static stability coefficient
0
.
63
κ
von K ́arm ́an constant
0
.
4
κ
Ratio of rms turbulent velocity to friction velocity
1
.
94
a
1
Empirical stability function coefficient
100
a
2
Empirical stability function coefficient
0
.
2
Pr
t,
0
Turbulent Prandtl number in neutral conditions
0
.
74
The EDMF scheme used here differs from the one described in Tan et al. (2018)
271
in the parameterizations of the eddy diffusivity
K
ψ
, the vertical pressure anomaly gra-
272
dients in (6) and (10), entrainment and detrainment, and the addition of turbulent en-
273
trainment
ˆ
E
ij
. The parameterization of the eddy diffusivity follows (9) and (37). The
274
entrainment parameterization is described in Cohen et al. (2020), and the treatment of
275
the pressure anomaly term is shown in (11). In addition, although the theoretical frame-
276
work presented here allows for the use of downdrafts, the implementation used in this
277
section decomposes the domain solely into one updraft and its turbulent environment.
278
LES are performed using PyCLES, an anelastic fluid solver in which the subgrid-
279
scale fluxes are treated implicitly by the WENO scheme used to discretize the prognos-
280
tic equations (Pressel et al., 2015). Implicit LES using WENO numerics have been shown
281
to result in higher effective resolution than other combinations of numerics and explicit
282
SGS closures (Pressel et al., 2017). Finally, LES results from previous model intercom-
283
parison projects are also reported where available.
284
4.1 Stable Boundary Layer
285
Statically stable conditions in the boundary layer inhibit convection, reducing the
286
EDMF scheme to a diffusive closure. In the implementation of the scheme, this trans-
287
lates to conditioning the surface updraft area fraction on the sign of the surface buoy-
288
ancy flux, such that it becomes zero in conditions of surface cooling. With no updrafts
289
or downdrafts, the only contribution to the subgrid-scale flux (2) comes from the envi-
290
ronmental downgradient turbulent flux (3). This leads to a high sensitivity of SCM sim-
291
ulations to changes in the mixing length formulation. Here we focus on the GEWEX At-
292
mospheric Boundary Layer Study (GABLS), discussed in Beare et al. (2006).
293
4.1.1 Simulation Setup
294
The initial and boundary conditions of the simulation are adapted from observa-
295
tions during the Beaufort and Arctic Seas Experiment (Curry et al., 1997) and follow
296
Beare et al. (2006). The velocity field is initialized as (
u
,
v
) = (
u
g
,
0), where the geostrophic
297
velocity is
u
g
= 8 ms
1
. The initial temperature sounding is given by a mixed layer
298
with potential temperature
θ
= 265 K up to 100 m, overlain by an inversion with a po-
299
tential temperature gradient of 10 K km
1
. The surface boundary condition is given by
300
constant cooling,
̇
θ
z
=0
=
0
.
25 K h
1
.
301
For both the SCM and LES, the domain height is 400 m. In the LES configura-
302
tion, the domain spans 400 m in both horizontal directions as well. The LES data is gen-
303
erated using an isotropic mesh with ∆
x
i
= 3
.
125 m resolution, which translates into
304
2
×
10
6
degrees of freedom. The full range of LES results from Beare et al. (2006), us-
305
–13–
ESSOAr | https://doi.org/10.1002/essoar.10502906.1 | Non-exclusive | First posted online: Sat, 2 May 2020 02:38:16 | This content has not been peer reviewed.