# Capillary levelling of thin liquid films of power-law rheology

Michael. C. Dallaston

###### Abstract

We find solutions that describe the levelling of a thin fluid film coating a
substrate, comprising a non-Newtonian power law fluid, under the influence of
surface tension. We consider the evolution from both periodic and localised
initial conditions. These solutions exhibit the generic property that the
profiles are weakly singular (that is, higher order derivatives do not exist),
at points where the pressure gradient vanishes. Numerical solutions of the
thin film equation are shown to approach one of the two cases for appropriate
initial conditions.

###### Contents

  1. 1 Introduction
  2. 2 Weakly singular nature and levelling solutions
    1. 2.1 Weakly singular nature
    2. 2.2 Close-to-uniform approximation
    3. 2.3 Periodic initial condition
      1. 2.3.1 Numerical computation
      2. 2.3.2 Close to Newtonian limit
    4. 2.4 Localised initial condition
  3. 3 Comparison with numerical simulations
  4. 4 Discussion

##  1 Introduction

The importance of free-surface lubrication or coating flows to industry and
science is well documented [Craster2009, Oron1997, Weinstein2004]. Examples of
applications include the construction of smooth, or deliberately patterned,
surfaces, describing the evolution of lava flows and ice sheets, among others.
The mathematical modelling of such phenomena using the lubrication
approximation, which allows the full governing partial differential equations
and boundary conditions to be collapsed to a single equation for the film
thickness, is also ubiquitous [Craster2009, Oron1997].

For small scale industrial applications in which surface tension is a dominant
effect, fluids are frequently better modelled using a non-Newtonian, rather
than Newtonian, rheological model. A very popular model is the power-law
model:

$\tau=\alpha|\dot{\gamma}|^{n-1}\dot{\gamma},$ (1)

where $\tau$ is the stress tensor, $\dot{\gamma}$ is the strain rate tensor,
$n$ is the power law rheology exponent, and $\alpha$ is a constant of
proportionality (see [Myers2005] for a discussion of different rheological
models). The Newtonian case is recovered when $n=1$, in which case $\alpha$ is
the fluid viscosity. When $n<1$, the effective viscosity is decreasing in the
strain rate, while for $n>1$, it is increasing. These two cases are thus known
as shear-thinning and shear-thickening, respectively. Shear thinning fluids
are more common in practice. The importance of rheology in the coating
industry is described in [Eley2005].

From a mathematical perspective, the evolution of a lubricating film,
comprising a power law fluid, under the effect of surface tension, has been
considered with a focus on droplet spreading [King2001], or inclined or
vertical substrates [Allouche2015, Balmforth2003, Sylvester1973]. The
competition between surface tension and other effects such as van der Waals
disjoining pressure, [Garg2017], and horizontally induced thermocapillary
stress [Mantripragada2022].

In this article we focus on the deceptively simple problem that is the
levelling of a film due to surface tension, in which a perturbed uniform film
tends to flatten out over time. In the Newtonian case, this phenomenon has
been studied in the context of viscometry [Benzaquen2013, Benzaquen2014,
McGraw2011, McGraw2012]; by measuring the evolution of a perturbed film over
time, the viscosity of the fluid may be determined. In these works, explicit
similarity solutions are found that model the levelling process. This same
application in the power-law fluid case has been studied to a lesser extent;
[Iyer1996] perform numerical calculations of the levelling problem with a
Carreau fluid and compare to experimental results. [Ahmed2015] perform a
purely numerical study on power law fluids. The similarity solutions
considered in the present study are also a case of nonlinear differential
equations whose existence is established in [Bernis1991]. However, the
explicit computation of similarity solutions that govern levelling behaviour,
and comparison with numerical solutions, has not previously been undertaken in
for power-law rheology.

Suitably nondimensionalised, the equation for a thin film of thickness
$h(x,t)$ and rheology (1) evolving purely due to surface tension in one
dimension is [Garg2017, King2001]:

$\frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=0,\qquad
q=h^{2+1/n}\left|\frac{\partial^{3}h}{\partial
x^{3}}\right|^{1/n}\operatorname{sgn}\left(\frac{\partial^{3}h}{\partial
x^{3}}\right),$ (2)

where we have explicitly defined the flux $q$. Here $n$ is the power law
exponent in the rheology; $n<1$ $n=1$, and $n>1$ correspond to the shear
thinning, Newtonian, and shear thickening cases, respectively. In the
derivation of (2), the third derivative arises as the gradient of the Laplace
pressure (the interface curvature being given by the second derivative in the
lubrication approximation).

The governing equation (2) has as an exact solution a perfectly uniform film
which can be taken to be $h=1$ identically (we are free in the
nondimensionalisation to specify the characteristic thickness as the vertical
length scale). In this work we consider the evolution of perturbations of this
uniform film, distinguishing two cases: the evolution of a perturbation that
is periodic in space, which is equivalent to a perturbation over an entire
finite domain with no-flux boundary conditions; and the evolution of an
initially localised perturbation in an infinite spatial domain (see Fig. 1).
These two cases generalise the results of a Newtonian liquid film, which are
modal linear stability of a flat interface, and the similarity solution
derived in [Benzaquen2013, Benzaquen2014], respectively. We numerically
determine the similarity solutions that describe each of these in the non-
Newtonian case. The computed solutions demonstrate a general property of
solutions to power-law thin film equations, namely, the non-existence of
higher derivatives (depending on the exponent $n$). As well as of mathematical
interest, this property is important to establish from a modelling
perspective, as it affects the accuracy and potentially the correctness of
numerical schemes. We also carry out numerical simulations of the fully time-
dependent thin film equation, and show that solutions with more general
initial conditions are attracted to the periodic or similarity solutions as
appropriate.

Figure 1: A schematic diagram of the two levelling solutions under
consideration: (a) levelling from an initially periodic perturbation of a flat
film, described in Section 2.3, and (b) levelling from an initially localised
perturbation, described by the similarity solutions in Section 2.4. In both
cases $h$ describes the height of the film, while $H$ describes the difference
from the uniform height of unity.

##  2 Weakly singular nature and levelling solutions

###  2.1 Weakly singular nature

In the Newtonian case, (2) is degenerate in that the flux vanishes at points
where the thickness $h$ vanishes. In the non-Newtonian case, (2) is degenerate
even if $h$ is not zero, at a point $x_{0}$ where the pressure gradient (that
is, the third derivative of $h$) is zero. A sensible expansion for $h$ near
such a point is

$h=h_{0}(t)+h_{1}(t)(x-x_{0})+h_{2}(t)(x-x_{0})^{2}+h_{r}(t)|x-x_{0}|^{r},\qquad
r\geq 3,$

where $r$ is an as-yet unknown power. On substitution into (2), the leading
order term in the flux $q$ is

$q=h_{0}^{2+1/n}\left[r(r-1)(r-2)|h_{r}||x-x_{0}|^{r-3}\right]^{1/n}\operatorname{sgn}(x-x_{0})+\ldots$

so that

$\dot{h}_{0}+h_{0}^{2+1/n}\left[r(r-1)(r-2)|h_{r}|\right]^{1/n}\left(\frac{r-3}{n}\right)|x-x_{0}|^{(r-3)/n-1}+\ldots=0.$

Generally $\dot{h}_{0}=h_{t}(x_{0},t)\neq 0$ at such a point, thus, we thus
require $r=3+n$. For non-integer $n$, we expect the film profile will be
weakly singular, such that $h(x,t)$ is $3+\lfloor n\rfloor$ times
differentiable. This is particular important for shear-thinning flows ($n<1$)
as then we only expect weak, rather than classical, solutions of the fourth
order equation (2) to exist. We will observe this singular behaviour in each
of the solutions constructed in this article. We note that as special cases,
if $\dot{h}_{0}=0$, then a singularity at higher order will occur as the
gradient of flux must balance with another term, while if $n$ is an integer
not equal to unity, then there will be higher order noninteger powers in the
flux that need to be balanced, requiring higher order non-integer powers in
the expansion of $h$. We do not consider these special cases further here.

###  2.2 Close-to-uniform approximation

We examine the behaviour of a film that is close to the uniform solution $h=1$
by writing

$h=1+\delta H(x,T),\qquad T=\delta^{1/n-1}t$ (3)

where $\delta\ll 1$. Then to $O(\delta)$,

$\frac{\partial H}{\partial T}+\frac{\partial}{\partial
x}\left[\left|\frac{\partial^{3}H}{\partial
x^{3}}\right|^{1/n}\operatorname{sgn}\left(\frac{\partial^{3}H}{\partial
x^{3}}\right)\right]=0.$ (4)

In the Newtonian case ($n=1$), (4) is linear, and no time rescaling is
necessary, but in the non-Newtonian case ($n\neq 1$), (4) remains nonlinear,
and the amplitude-dependent scaling in time is necessary for the two terms in
the equation to be in balance. The near-level approximation (4) has the same
property of the original equation (2), in that $H$ must be weakly singular at
any point where the third spatial derivative is zero in order for the
evolution to be well defined.

We now construct two special classes of solutions to (4), generalising the
known results for the Newtonian case (see Fig. 1). Firstly, we consider
solutions periodic in $x$, which are relevant for initial conditions that are
periodic or on bounded domains with no-flux conditions. Secondly, we construct
similarity solutions that are appropriate for initial conditions that are
spatially localised.

###  2.3 Periodic initial condition

We start with solutions that are periodic in $x$. In the Newtonian case
($n=1$) this is equivalent to computing the linear stability of the uniform
solution. For a given wavenumber $k$ (so that the period is $2\pi/k$) we have

$H(x,t)=\bar{A}\mathrm{e}^{-k^{4}t}\cos(kx),$

where the amplitude $\bar{A}$ is arbitrary (equivalent to time-translational
invariance). For $n\neq 1$, the near-level behaviour is essentially nonlinear,
and linear stability approaches are not appropriate. Instead, we assume a
solution of (4) more generally in which the time and space dependence may be
separated, or, equivalently, a similarity solution in which the spatial scale
remains fixed. The appropriate ansatz depends on whether the fluid is shear
thinning or thickening:

$H(x,T)=\begin{cases}\displaystyle\left(\frac{n}{1-n}\right)^{n/(1-n)}T^{-n/(1-n)}\bar{H}(x),&n<1\\\
\displaystyle\left(\frac{n}{n-1}\right)^{-n/(n-1)}(-T)^{n/(n-1)}\bar{H}(x),&n>1.\end{cases}$
(5)

where $\bar{H}(x)$ is a periodic, but as yet undetermined, function (the
$n$-dependent constant is included to simplify the problem for $\bar{H}$). The
above ansatz represents infinite-time but algebraic decay for the shear
thinning ($n<1$) case, and finite-time (at $T=0$) levelling for the shear
thickening $n>1$ case; of course, the solutions are invariant to translations
in time, and the finite levelling time will depend on the initial amplitude.

Under the ansatz (5), $\bar{H}$ satisfies the ordinary differential equation

$[|\bar{H}^{\prime\prime\prime}|^{1/n}\operatorname{sgn}(\bar{H}^{\prime\prime\prime})]^{\prime}=\bar{H},$

for both shear-thinning and thickening. Let
$U=|\bar{H}^{\prime\prime\prime}|^{1/n}\operatorname{sgn}(\bar{H}^{\prime\prime\prime})$,
then $U$ satisfies

$U^{\prime\prime\prime\prime}=|U|^{n}\operatorname{sgn}U,$ (6)

with $\bar{H}$ readily found from $U$ by $\bar{H}=U^{\prime}$. Since any
solution to (6) may be scaled onto another solution by $x\mapsto\lambda x$,
$U\mapsto\lambda^{4/(1-n)}U$, it suffices to find a solution with period
$2\pi$, which can then be scaled to find the solution for any other period.

####  2.3.1 Numerical computation

Corresponding to the behaviour of the original problem near a singular point
described in Section 2.1 for noninteger $n$, solutions to (6) will be weakly
singular at points $x_{0}$ where $U(x_{0})=0$, with an expansion of the form

$U=U_{1}(x-x_{0})+U_{2}(x-x_{0})^{2}+U_{3}(x-x_{0})^{3}+U_{r}|x-x_{0}|^{4+n}\operatorname{sgn}(x-x_{0}).$

In the shear-thinning case, $U$ will thus be four (but not five)-times
differentiable, corresponding to the profile $\bar{H}$ (and so the time-
dependent profile $h$) being three but not four-times differentiable at points
where the pressure gradient $h_{xxx}=0$. However, since (6) is fourth order,
this singularity is weak enough to be able to proceed (that is, we have
essentially integrated the equation by solving for the antiderivative of
$\bar{H}$).

Assuming symmetry, we find the solution on the half period $0<x<\pi$ with the
conditions

$U(0)=U^{\prime\prime}(0)=U(\pi)=U^{\prime\prime}(\pi)=0.$ (7)

In the half-period, $U$ may be taken to be positive. We calculate solutions to
the boundary value problem (6), (7) numerically using the MATLAB function
bvp4c. An initial guess is provided close to $n=1$ from the asymptotic
solution described below (Section 2.3.2). A continuous branch of solutions is
then constructed by numerical continuation for decreasing and increasing power
law exponent $n$, using the previously computed solution as an initial guess
for the next.

In Fig. 2 we plot the results of this computation, focussing on the values
$n=1/2$ and $n=3/2$. We include both the profiles $\bar{H}=U^{\prime}$, as
well as the fourth derivative
$\bar{H}^{\prime\prime\prime\prime}=n|U|^{n-1}U^{\prime}$. The profiles
themselves are visually difficult to distinguish from the sinusoidal profiles
that are exact for the Newtonian case ($n=1$), but have non-arbitrary,
$n$-dependent amplitude. In plotting the fourth derivative, the difference
between the cases is stark, in particular, that the shear thinning case
($n=1/2$) is singular at the peak and trough of the periodic profile
($x=0,\pi$), while for the shear-thickening case ($n=3/2$), the fourth
derivative is bounded and continuous (although the next derivative will indeed
be singular). We also plot the amplitude $\bar{A}$, calculated as half of the
peak to trough distance $\bar{A}=(\bar{H}(0)-\bar{H}(\pi))/2$; this value
varies only weakly in the range of $n$ values calculated.

Figure 2: Spatial profiles $\bar{H}(x)$ representing levelling solutions for a
periodic initial condition, $n=1/2$ (shear-thinning) and $n=3/2$ (shear-
thickening). (a) The solution profiles (calculated from solutions $U$ to (6)
with $\bar{H}=U^{\prime}$). These profiles are close to sinusoidal. (b) The
fourth derivative $\bar{H}^{\prime\prime\prime\prime}$ show the singular
nature of these profiles, with $\bar{H}$ only three times differentiable for
$n=1/2$, and only four times differentiable for $n=3/2$. (c) The dependence of
amplitude $\bar{A}$ on power law exponent $n$.

####  2.3.2 Close to Newtonian limit

In the above method, the starting point of the numerical solutions was found
by determining the asymptotic solution to (6) in the close-to-Newtonian limit
$n\to 1$. This approximation is also of interest in itself. Defining
$n=1+\epsilon$, where $\epsilon\ll 1$:

$U^{n}=U\mathrm{e}^{\epsilon\log U}=U(1+\epsilon\log
U+\tfrac{1}{2}\epsilon^{2}\log(U)^{2}+\ldots),$

so that if $U=U_{0}+\epsilon U_{1}+\epsilon^{2}U_{2}+\ldots$, we have

$U_{0}^{\prime\prime\prime\prime}=U_{0},\qquad
U_{1}^{\prime\prime\prime\prime}-U_{1}=U_{0}\log(U_{0}).$

Applying the boundary conditions (7) to $U_{0}$, we find

$U_{0}=\bar{A}\sin(x),$

with the amplitude $\bar{A}$ determined by a solvability condition for
$U_{1}$. Applying the same boundary conditions to $U_{1}$, we must have that
(by the Fredholm alternative)

$\int_{0}^{\pi}U_{0}\log(U_{0})\sin(x)\,\mathrm{d}x=\frac{\pi\bar{A}}{2}\left[\log\bar{A}+\frac{1}{2}(1-\log
4)\right]=0.$

Thus either $\bar{A}=0$ (the trivial solution), or

$\bar{A}=\mathrm{e}^{(\log 4-1)/2}=2\mathrm{e}^{-1/2}=1.21\ldots$ (8)

This value agrees with the numerically computed results near $n=1$ (see Fig.
2). We note that in this asymptotic result we have not considered the region
near $U=0$ where $U$ will be weakly singular, as it does not play a role in
selecting the amplitude $\bar{A}$.

###  2.4 Localised initial condition

We now find solutions relevant for an initial condition that is localised in
an infinite spatial domain. In this case, relevant solutions are similarity
solutions that describe how an initial peak spreads over time.

For any value of $n$, similarity solutions to (4), symmetric around the point
$x=0$, take the form

$H(x,t)=H_{0}T^{-\alpha}F(\eta),\qquad\eta=H_{0}^{(1-n)/(n+3)}\alpha^{-n/(n+3)}\frac{x}{T^{\alpha}},\qquad\alpha=\frac{n}{4},$
(9)

where $\eta,F$ are the similarity variables, $H_{0}$ is an arbitrary amplitude
included so that we can specify $F(0)=1$, and the similarity exponent $\alpha$
is determined from the two conditions that the terms in (4) have the same time
dependence, and that mass is conserved. Substituting (9) into (4) results in
the ordinary differential equation for similarity profiles:

$(\eta
F)^{\prime}=[|F^{\prime\prime\prime}|^{1/n}\operatorname{sgn}(F^{\prime\prime\prime})]^{\prime},$

subject to symmetry conditions at zero. This equation may be integrated once
and the boundary conditions used to result in the third-order equation

$F^{\prime\prime\prime}=|\eta F|^{n}\operatorname{sgn}(F),\ \eta>0,\qquad
F(0)=1,\ F^{\prime}(0)=0.$ (10)

A third condition is imposed by requiring solutions to decay as
$\eta\to\infty$. Similarly to the periodic case, since we only deal with third
derivatives of $F$ in (10), the singularities that are present (in particular
in the fourth derivative for $n<1$) will not prevent numerical computation of
the solution.

In the Newtonian case, the similarity solution has been previously identified
in [Benzaquen2013, Benzaquen2014]; In this case the similarity exponent
$\alpha=1/4$, and the equation for the profile (10) is linear:

$F^{\prime\prime\prime}=\eta F$

which, given $F(0)=1$, $F^{\prime}(0)=0$, and decay at infinity, has exact
solution in terms of hypergeometric functions:

$F(\eta)={}_{0}F_{2}\left([],[1/2,3/4],\eta^{4}/64\right)-\frac{\Gamma(3/4)}{4\Gamma(5/4)}\eta^{2}{}_{0}F_{2}\left([],[5/4,3/2],\eta^{4}/64\right).$

For $n\neq 1$, numerical solutions to (10) are obtained by shooting from
$\eta=0$. Since $F(0)=1,F^{\prime}(0)=0$, there is one parameter
$F^{\prime\prime}(0)$ that must be determined numerically by requiring that
$F\to 0$ as $\eta\to\infty$. We use ode45 in MATLAB to solve the equation up
to a sufficiently large value of $\eta$, combined with a root finding method
to determine the appropriate value of $F^{\prime\prime}(0)$. The solutions so
found are plotted in Fig. 3, including both the profiles $F$ and the fourth
derivative, calculated from the numerical solution using

$F^{\prime\prime\prime\prime}=n|\eta F|^{n-1}(F+\eta F^{\prime}).$

As for periodic solutions, plotting this fourth derivative shows the presence
of singularities at higher order. As the similarity solution oscillates an
infinite number of times as $\eta\to\infty$, there are an infinite sequence of
such singularities.

Figure 3: Similarity solutions for thin film levelling of an initially
localised perturbation, for power law exponents $n=1/2$ and $n=3/2$. While the
solution profiles (a) appear qualitatively similar, plotting the fourth
derivative (b) demonstrates the weakly singular nature of the similarity
solutions.

##  3 Comparison with numerical simulations

In order to test the universality of the solutions found in Section 2, we
simulate the time-dependent lubrication equation (2). We use a finite volume
method, which naturally solves the weak formulation of a conservative
transport equation. This property is particularly important in the present
case as for $n<1$, we do not expect (2) to have classical (four-times
differentiable) solutions. We consider a domain $x\in[0,L]$ divided into $N$
cells of width $\Delta x=L/N$, with the $j$th cell centre at $x_{j}=\Delta
x(j-1/2)$, $1<j<N$. The finite difference expression for the third derivative
of $h$ on the face between the $j$th and $(j+1)$th cell is

$h_{xxx}\approx\frac{1}{\Delta
x^{3}}\left(-h_{j-1}+3h_{j}-3h_{j+1}+h_{j+2}\right),$

where $h_{j}$ is taken to be the representative value of $h$ in cell $j$,
while the value of the thickness itself $h$ on the face is given by the
average $(h_{j+1}+h_{j})/2$. These expressions are used to approximate the
fluxes $q_{j}$. The transport equation (2) then gives

$\frac{\mathrm{d}h_{j}}{\mathrm{d}t}=\frac{1}{\Delta
x}\left(q_{j}-q_{j-1}\right).$

No flux conditions are imposed at the two ends $x=0$ and $x=L$, while no-slope
conditions ($h_{x}=0$) are used to define ghost node values that allow the
third derivative to be computed on all internal faces. The numerical scheme is
then advanced in time using MATLAB’s ode15s implicit time-stepping algorithm.

Using this method we compute solutions of (2) for the two values of the power
law exponent we have focussed on thus far ($n=1/2$ and $n=3/2$), and initial
conditions that test the behaviour of periodic and localised solutions,
respectively. To test the behaviour of periodic solutions, we choose a domain
size of $L=2\pi$ and an initial condition

$h(x,0)=1+0.25\cos\left(x\right).$

From the ansatz (5) we predict that the amplitude $A(t)$, which we calculate
numerically as $(\max(h)-\min(h))/2$, will either decay algebraically or
undergo finite-time levelling, depending on $n$:

$A\sim\begin{cases}\displaystyle\bar{A}\left(\frac{n}{1-n}\right)^{n/(1-n)}t^{-n/(1-n)},&n<1\\\
\displaystyle\bar{A}\left(\frac{n}{n-1}\right)^{-n/(n-1)}(t_{0}-t)^{n/(n-1)},&n>1\end{cases},$

where $\bar{A}$ is the ($n$-dependent) amplitude of the solution to the
profile function $\bar{H}(x)$, and $t_{0}$ is the finite levelling time. For
$n=1/2$ and $n=3$, then, these become

$A\sim\begin{cases}1.213t^{-1},&n=1/2\\\
0.044(t_{0}-t)^{3},&n=3/2\end{cases}.$ (11)

In this case, the prefactor is determined independently of the initial
condition (although not of the domain size). The levelling time $t_{0}$ is
initial condition-dependent, however. In Fig. 4 we show the results of the
numerical simulation. The observed behaviour of the amplitudes are in
agreement with the prediction (11), indicating that the levelling solutions we
found in Section 2 are accurate, and are attractors of more generic initial
conditions.

Figure 4: Numerical simulation of a periodic solution for (a,b) $n=1/2$, and
(c,d) $n=3/2$, using the method described in Section 3. (a,c) depict the
profile evolution, with initial condition shown as the black dashed curve. In
(b) we observe the predicted power-law decay in amplitude for shear-thinning
flow, with $A\sim 1.246t^{-1}$ (shown as a black dotted line). In (d) we
observe the predicted finite time levelling in the shear-thickening case, with
$A\sim 0.044(t_{0}-t)^{3}$ (shown as a black dotted line); here the finite
levelling time is approximately $t_{0}\approx 1.8$. Blue circles in (b.d)
correspond to the times at which profiles are plotted in (a,c), respectively.

To test the behaviour of localised initial conditions, we start with a
Gaussian initial condition

$h(x,0)=1+0.3\exp[-4(x-L/2)^{2}],$

for both $n=1/2$ and $n=3/2$, with domain size $L$ sufficiently large so as to
avoid the boundary having an effect on the evolution. The similarity ansatz
(9) predicts the amplitude, defined by $A=\max(h)-1$, will go as

$A\sim\text{constant}\cdot t^{-n/4}.$ (12)

As $H_{0}$ is arbitrary in (9), the prefactor in this relation is dependent on
initial condition. We plot the results of these simulations in Fig. 5 for
$n=1/2$, and Fig. 6 for $n=3/2$. In each case, the profiles rapidly evolve
toward the relevant similarity solution (with the appropriate value of $H_{0}$
estimated from the numerical solution at the final time), and the behaviour of
the amplitude $A$ approaches the predicted power law (12). Again, the strong
agreement indicates the similarity solution is an attractor for initially
localised perturbations on an infinite domain.

Figure 5: Numerical simulation of a localised initial perturbation for the
shear-thinning case ($n=1/2$). (a) The solution profiles themselves; (b) the
rescaled profiles (according to (9)) collapse onto a single curve which
matches the similarity profile computed in Section 2.4 (black circles). (c)
The power-law decay of the amplitude $A=\max(h)-1$ also matches the value
predicted by the similarity ansatz. Blue circles in (c) correspond to the
profiles plotted in (a,b). Figure 6: Numerical simulation of a localised
initial perturbation for the shear-thickening case ($n=3/2$), with, similar to
Fig. 5: (a) The profiles, (b) the rescaled profiles which match with the
similarity solution (black circles), and (c) the power-law decay of the
amplitude $A$, which matches the predicted value. Blue circles in (c)
correspond to the profiles plotted in (a,b).

##  4 Discussion

By explicitly calculating similarity solutions that describe levelling, we
have determined the power-law decay of the amplitude of perturbations of a
flat film. As the decay depends on the power in the rheology (1), our results
imply that levelling experiments could be used to determine the rheology of
fluids, in a similar way as they are used to determine the viscosity of
presumed-Newtonian fluids [Benzaquen2013, Benzaquen2014].

We have also demonstrated that solutions to lubrication equations with power-
law rheology will be weakly singular at points where the pressure gradient
vanishes. This will be also true in more complicated models, for instance
those that feature disjoining pressure, Marangoni forces, or gravity. Despite
the weakly singular nature of solutions to (2), previous studies have
successfully computed solutions numerically. Indeed, the distinct infinite-
time and finite-time levelling, for shear-thinning and thickening fluids,
respectively, can be observed in the numerical results of [Ahmed2015],
although they do not interpret them as such. Numerical studies on power-law
thin films are likely to numerically regularise the singularities, or be
formulated in such a way that they correctly solve the weak version of the
equation (as our method does). However, as the singular behaviour becomes more
severe as $n$ is reduced, it is likely that it will be important to take into
account for studies that wish to examine the highly shear-thinning limit
($n\to 0$). On the other hand, studies on inclined or vertical planes, or with
imposed tangent stress, may avoid this issue by never having points at which
the pressure gradient vanishes.

The rheology itself may be explicitly regularised, using, for example, a
Carreau-type model [Myers2005], which acts like a power law fluid except for
small strain rate, where it behaves in a Newtonian manner. For more complex
rheology, the flux becomes more complicated to compute, although it is
possible to make progress using the method described in [Pritchard2015] (see
also [Hinton2022] for the two-dimensional case). Under such a regularisation,
a thin film may be expected to level in the same way as described in this
paper, until the film becomes very close to level, at which point the
Newtonian regime will be reached, and perturbations will switch to exponential
decay according to the linear stability analysis of the Newtonian case.

Extending to two spatial dimensions is also of interest. Presumably, there
exists a radially symmetric similarity solution that would describe the decay
of localised perturbations. For periodic perturbations, one could imagine a
radial solution for a thin film in a finite circular domain, or, for a more
general shape, one would have to solve the nonlinear elliptic equation that is
the counterpart to (6) in higher dimensions. The power law decay of the
amplitude in time, however, would be universal.

### Author address


Published by Journal of the Atmospheric Sciences gray6.5,-200.7501.2

# The mean wind and potential temperature flux profiles in convective boundary
layers

###### Abstract

We develop innovative analytical expressions for the mean wind and potential
temperature flux profiles in convective boundary layers (CBLs). CBLs are
frequently observed during daytime as the Earth’s surface is warmed by solar
radiation. Therefore, their modeling is relevant for weather forecasting,
climate modeling, and wind energy applications. For CBLs in the convective-
roll dominated regime, the mean velocity and potential temperature in the bulk
region of the mixed layer are approximately uniform. We propose an analytical
expression for the normalized potential temperature flux profile as a function
of height, using a perturbation method approach in which we employ the
horizontally homogeneous and quasi-stationary characteristics of the surface
and inversion layers. The velocity profile in the mixed layer and the
entrainment zone is constructed based on insights obtained from the proposed
potential temperature flux profile and the convective logarithmic friction
law. Combining this with the well-known Monin-Obukhov similarity theory allows
us to capture the velocity profile over the entire boundary layer height. The
proposed profiles agree excellently with large-eddy simulation results over
the range of $-L/z_{0}\in[3.6\times 10^{2},0.7\times 10^{5}]$, where $L$ is
the Obukhov length and $z_{0}$ is the roughness length.

##  1 Introduction

Convective boundary layers (CBLs) are frequently observed during daytime as
the Earth’s surface is warmed by solar radiation (Stull, 1988). Due to their
frequent occurrence, the fundamental understanding of CBLs is highly relevant
to agriculture, architectural design, aviation, climate modeling, weather
prediction, and wind energy applications, to name a few. The modern scientific
literature on CBLs goes back over 100 years. Initially, the focus was on low-
altitude measurements, and with the introduction of more advanced measurement
techniques, the focus gradually shifted upwards. However, only after the
introduction of large eddy simulations (LES) in the early seventies, it has
become widely accepted that thermodynamic indicators are most suitable to
identify the different CBL regions (LeMone et al., 2019). However, obtaining
analytical profiles that describe the wind and potential temperature flux in
the entire CBL has remained challenging due to the different flow physics in
the various CBL regions.

The CBL can be subdivided into three layers (excluding the roughness
sublayer), i.e. the surface layer, the mixed layer, and the entrainment zone
(see figure 1). The surface layer is characterized by a superadiabatic
potential temperature gradient and a strong wind shear, which is usually
described by the Monin-Obukhov similarity theory (MOST, Monin and Obukhov,
1954). According to the MOST the non-dimensional wind speed and potential
temperature gradient profiles are universal functions of the dimensionless
height $z/L$, where $z$ is the height above the surface and $L$ is the surface
Obukhov length (Obukhov, 1946; Monin and Obukhov, 1954). Many studies have
pointed out that the MOST does not explain all important surface-layer
statistics under convective conditions (Panofsky et al., 1977; Khanna and
Brasseur, 1997; Johansson et al., 2001; McNaughton et al., 2007; Salesky and
Anderson, 2020; Cheng et al., 2021) or very stable conditions (Mahrt et al.,
1998; Cheng et al., 2005). In particular, the normalized wind gradient
$\phi_{m}=(\kappa z/u_{*})(\partial U/\partial z)$ depends both on $z/L$ and
$z_{i}/L$ (Khanna and Brasseur, 1997; Johansson et al., 2001), where $z_{i}$
is the height of the inversion layer (see figure 1). Nevertheless, MOST is
still widely used in numerical weather prediction and climate models (Salesky
and Anderson, 2020), and thus will be used in the theoretical analysis and
numerical simulations of this study. MOST applies only to the surface layer,
and for it to be applicable, the absolute value of the Obukhov length $L$ must
be smaller than the height of the surface layer. Therefore we only consider
the CBL with $-z_{i}/L\gg 1$. In particular, we focus on the convective-roll
dominant regime with $-z_{i}/L\gtrsim 10$ (Salesky et al., 2017). Furthermore,
we focus on dry and cloud-free CBLs to avoid complications due to physical
processes like evaporation, precipitation, and cloud formation.

\begin{overpic}[width=346.89731pt]{fig1} \end{overpic} Figure 1: Profiles of
the potential temperature $\Theta$, wind speed $U_{\rm mag}$, and potential
temperature flux $q$ in the CBL. The vertical lines with double arrows
indicate different length scales in the CBL, namely, from left to right, the
Obukhov length $L$, the lowest height where the potential temperature flux
first reaches zero, $h_{1}$, the inversion layer height at which the potential
temperature flux reaches its minimum value, $z_{i}$, and the height where the
potential temperature flux first recovers zero, $h_{2}$. The background color
indicates the magnitude of the potential temperature flux $q(z)$ for case 2,
see Table 1.

The mixed layer is characterized by intense vertical mixing caused by warm air
thermals rising from the ground. Within the mixed layer, the magnitude of the
mean velocity is much larger than the variations in the mean velocity. Thus,
for applications where the mean velocity gradient is unimportant, the wind
speed and potential temperature can be regarded as uniform (Kaimal et al.,
1976; Salesky et al., 2017). This insight is incorporated in various CBL
models (Lilly, 1968; Deardorff, 1973; Stull, 1976; Deardorff, 1979; Tennekes
and Driedonks, 1981; Garratt et al., 1982). The entrainment zone is
characterized by entrainment of air from the free atmosphere. Deardorff et al.
(1980) found in laboratory experiments that the ratio of the entrainment zone
thickness to the depth of the mixing layer decreases asymptotically with
increasing Richardson number $Ri$ as follows
$(h_{2}-h_{1})/h_{1}=0.21+1.31/Ri$ (see figure 1 for the definitions of
$h_{1}$ and $h_{2}$). This ratio is essential for developing entrainment
models and has been studied extensively (Lilly, 1968; Sullivan et al., 1998;
Zilitinkevich et al., 2012; Haghshenas and Mellado, 2019). The potential
temperature flux profile decreases linearly with height and becomes negative
in the entrainment zone. The entrainment flux ratio $\Pi_{m}$, which is
defined as the ratio of the potential temperature flux at the inversion layer
height to its value at the ground, turns out to be nearly constant, i.e.
$\Pi_{m}\approx-0.2$ (Stull, 1976; Sorbjan, 1996; Conzemius and Fedorovich,
2006; Sun and Wang, 2008; LeMone et al., 2019). Note that the inversion layer
is the upper region of the entrainment zone in which the potential temperature
flux increases steeply from its minimum value at $z=z_{i}$ to zero at
$z=h_{2}$ (see figure 1).

The geostrophic wind $(U_{g},V_{g})$ and the friction velocity $u_{*}$ are
usually connected through the well-known geostrophic drag law, which was
initially derived for neutral boundary layers (Rossby and Montgomery, 1935;
Blackadar and Tennekes, 1968; Tennekes and Lumley, 1972) and later extended to
include buoyancy effects (Zilitinkevich, 1969). To include the effect of
unsteadiness, Zilitinkevich and Deardorff (1974) and Arya (1975) proposed to
replace the Ekman depth $u_{*}/|f|$ in the geostrophic draw law, where $f$ is
the Coriolis parameter, with the time-dependent inversion layer height
$z_{i}$. However, significant disparities were observed between the
geostrophic drag law for CBLs and measurement data (Zilitinkevich, 1975).
Garratt et al. (1982) derived a relationship for the velocity defects in the
mixed layer using a three-layer CBL model, which accounts for the effects of
entrainment, baroclinity, advection, and local acceleration. In their
formulation, the velocity defects are defined as the differences between the
mixed-layer averaged winds and the geostrophic winds. They proposed a
geostrophic drag law to relate the geostrophic winds and the friction velocity
based on the assumption that the mean velocity at the top of the surface layer
is equal to that in the mixed layer. In addition, an empirical stability
function $\psi_{m}$, which may be inaccurate for large values of $-z/L$, is
employed.

Recently, Tong and Ding (2020) analytically derived the convective logarithmic
friction law from first principles. They identified three scaling layers for
the CBL with $-z_{i}/L\gg 1$: the outer layer, the inner-outer layer, and the
inner-inner layer. The characteristic length scales for these three layers are
the inversion layer height $z_{i}$, the Obukhov length $L$, and the roughness
length $z_{0}$, respectively. The mixed-layer mean velocity scale $U_{m}$ and
the geostrophic wind component $V_{g}$ are the characteristic streamwise and
spanwise velocity scales in the outer layer. The difference between the
horizontally and temporally averaged velocity $U(z)$ and $U_{m}$ is the mixed-
layer velocity-defect law, which has a velocity scale of $u_{*}^{2}/w_{*}\ll
U_{m}$ with $w_{*}$ the convective velocity. This indicates that $U_{m}$ is
very close to the mean velocity $U(z)$ in the mixed layer. For the inner-outer
layer they derived the surface-layer velocity-defect law, which states that
the velocity defect $U-U_{m}$ scales with $u_{*}$. The convective logarithmic
friction law is derived from matching the law of the wall in the inner-inner
layer with the velocity-defect law in the surface-layer. This exact leading-
order result relates the friction velocity $(u_{*})$ to the mixed-layer
velocity scale $(U_{m})$. The difference between $U_{g}$ and $U_{m}$ scales as
$(u_{*}^{2}w_{e})/(fz_{i})^{2}$, where $w_{e}$ is the entrainment velocity and
$V_{g}$ scales as $-u_{*}^{2}/(fz_{i})$. Thus, up to non-dimensional
coefficients, one can relate the geostrophic velocities $(U_{g},V_{g})$ to
$u_{*}$. As Tong and Ding (2020) do not consider the effects of the
entrainment zone the velocity profiles are only valid for $z/z_{i}<0.4$.

Various time-dependent models have been developed to explicitly account for
entrainment processes at the top of CBLs (Troen and Mahrt, 1986; Noh et al.,
2003; Hong et al., 2006). For example, the counter-gradient transport method
(Holtslag and Moeng, 1991) and the eddy-diffusivity mass-flux approach
(Siebesma et al., 2007; Li et al., 2021) are widely used in coarse-resolution
climate models. In general, the potential temperature is time-dependent
(Lilly, 1968) and the entrainment velocity can affect the mean wind speed in
the mixed layer (Tong and Ding, 2020). However, the velocity and potential
temperature flux profiles are quasi-stationary, and therefore similarity
theory can be employed to obtain analytical expressions for these profiles
shapes (Zilitinkevich and Deardorff, 1974; Arya, 1975; Zilitinkevich et al.,
1992).

In this study, we focus on the the derivation of analytical expressions for
the mean velocity and potential temperature flux profiles in cloud-free CBLs.
We use a perturbation method approach to construct an analytical expression
for the normalized potential temperature flux profile as a function of height,
taking into account the characteristics of both the surface layer and the
capping inversion layer. The depth of the entrainment zone is connected to the
convective logarithmic friction law to obtain analytical expressions for the
velocity profile in the mixed layer and the entrainment zone. As remarked
previously, the surface layer is still described by the MOST.

The organization of the paper is as follows. In Section 2 we obtain analytical
expression for the potential temperature flux and wind profiles. In Section 3
we validate the proposed profiles against LES. The conclusions are given in
Section 4.

##  2 Theory

###  2.1 Potential temperature flux profile

The potential temperature flux profile provides a precise and convenient
demarcation between the mixed layer and the entrainment zone of the CBL
(Deardorff, 1979; Deardorff et al., 1980). Figure 1 shows a definition of the
various length scales in the CBL. Previous studies (Kaimal et al., 1976;
Deardorff et al., 1980; Moeng and Sullivan, 1994; Noh et al., 2003; Garcia and
Mellado, 2014; Haghshenas and Mellado, 2019) showed that the potential
temperature flux (including both the turbulent part and the diffusive part) in
CBLs decreases linearly from its maximum value at the surface to a minimum
value at $z=z_{i}$, and then increases steeply to zero in a narrow region
$z_{i}\leq z\leq h_{2}$ at the top of the boundary layer (figure 1). For
typical CBLs the condition $|\textrm{d}z_{i}/\textrm{d}t|\ll w_{*}$ holds,
which implies that the boundary layer is quasi-stationary (Nieuwstadt et al.,
2016, Section 7.6). Besides, the potential temperature flux $q$ is fixed at
the surface, and its value at the inversion layer height is nearly a constant
fraction of the value at the ground $q_{w}$. Therefore, the normalized
potential temperature flux $q(z,t)/q_{w}$ only depends on the similarity
variable $\xi=z/h_{2}(t)$, i.e.

$q(z,t)/q_{w}=\Pi(\xi),$ (1)

where the form of $\Pi$ remains to be determined. Using the potential
temperature equation, we derive below an ordinary differential equation (ODE)
for the determination of $\Pi$. However, it is important to emphasize that
this does not mean $\Pi$ is independent of time as the similarity variable
$\xi=z/h_{2}(t)$ is still time-dependent.

Under the assumption of horizontal homogeneity, the potential temperature
equation reduces to

$\frac{\partial\Theta}{\partial t}=-\frac{\partial q}{\partial
z}=-\frac{q_{w}}{h_{2}}\Pi^{\prime},$ (2)

where $\Pi^{\prime}=\partial\Pi/\partial\xi$. In the mixed layer the potential
temperature $\Theta$ is almost spatially uniform and hence the left-hand side
of Eq. (2) is independent of $z$. Therefore, the governing equation of the
potential temperature flux in the mixed layer can be approximated as

$-\Pi^{\prime}=c_{\Pi},$ (3)

where $c_{\Pi}$ is the gradient of the normalized potential temperature flux
in the mixed layer (see figure 1). In the mixed layer, it is well-known that
the eddy-diffusivity approach cannot adequately describe the potential
temperature flux as the gradient of the potential temperature nearly vanishes
(Wyngaard, 2010). In contrast, in the inversion layer the potential
temperature gradient is dominant such that the potential temperature flux can
be approximated by $q=-\nu_{\theta}\partial\Theta/\partial z$. Here
$\nu_{\theta}\propto|w_{e}|(h_{2}-z_{i})$ is the eddy diffusivity and
$w_{e}=\textrm{d}z_{i}/\textrm{d}t$ is the entrainment velocity. That the
ratio $h_{2}/z_{i}$ is approximately constant implies that
$\partial\xi/\partial
t=-(z/h_{2}^{2})(h_{2}/z_{i})(\textrm{d}z_{i}/\textrm{d}t)\approx-
w_{e}/z_{i}$. From the zero-order jump model $\textrm{d}z_{i}^{2}/\textrm{d}t$
is independent of time (Nieuwstadt et al., 2016, Section 7.6), such that
$\nu_{\theta}$ is approximately constant. Then, by taking the vertical
derivative of Eq. (2), we obtain that

$\epsilon\Pi^{\prime\prime}-\Pi^{\prime}=0,$ (4)

where $\epsilon$ is a small dimensionless parameter,

$\epsilon\equiv c\frac{h_{2}-z_{i}}{h_{2}}\ll 1,$ (5)

with $c=(\nu_{\theta}z_{i})/[-w_{e}h_{2}(h_{2}-z_{i})]=1/2$ being an empirical
constant that is determined by comparing the model profiles to the simulation
results. This indicates that the parameter $\epsilon$ represents the half
thickness of the inversion layer normalized by the boundary layer depth (see
figure 1). To get the potential temperature flux profile in the entire
boundary layer, we combine Eqs. (3) and (4), which leads to the following
second-order ODE for the potential temperature flux,

$\epsilon\Pi^{\prime\prime}-\Pi^{\prime}=c_{\Pi},\quad\Pi(0)=1,\quad\Pi(1)=0.$
(6)

The solution of Eq. (6) reads

$\displaystyle\Pi=1-c_{\Pi}\xi+(c_{\Pi}-1)\frac{e^{\xi/\epsilon}-1}{e^{1/\epsilon}-1},\quad
0\leq\xi\leq 1.$ (7)

Since $\epsilon\ll 1$, the value of $\Pi$ in the bulk of the mixed layer can
be approximated as

$\displaystyle\Pi\approx 1-c_{\Pi}\xi.$ (8)

Similarly, since $\Pi=0$ at $\xi=h_{1}/h_{2}$, the slope $c_{\Pi}$ reduces to

$\displaystyle c_{\Pi}=h_{2}/h_{1}>1.$ (9)

Therefore, the ratio of the entrainment zone thickness to the mixing layer
depth is

$R\equiv({h_{2}-h_{1}})/{h_{1}}=c_{\Pi}-1>0.$ (10)

Deardorff et al. (1980) found in laboratory experiments that the value of $R$
is between 0.2 and 0.4. Furthermore, the entrainment flux ratio $\Pi_{m}$,
i.e. the minimum value of $\Pi$, can be approximated as

$\Pi_{m}\approx 1-c_{\Pi}(1-2\epsilon)\approx-(R-2\epsilon).\ $ (11)

Stull (1976) and Sorbjan (1996) found that $-0.3\leq\Pi_{m}\leq-0.1$. This
results is consistent with the LES results of Sullivan and Patton (2011) with
$\Pi_{m}\approx-0.2$, the empirical results of Lenschow (1974) with
$\Pi_{m}=-0.1$, and the direct numerical simulation results of Garcia and
Mellado (2014) with $\Pi_{m}\approx-0.12$.

We note that the perturbation method approach to model the potential
temperature flux profile was recently introduced by Liu et al. (2021b) for
conventionally neutral atmospheric boundary layers where the surface potential
temperature flux is always zero. However, it should be noted that in the CBLs
under consideration, the surface is heated and thermal plumes are generated at
the ground, resulting in significantly different turbulence generation
mechanisms. The applicability of the perturbation method approach to model the
strong inversion layer relies on its ability to capture the strong gradients
in the inversion layer. Here we used the second-order ODE defined by Eq. (6)
to model the potential temperature flux profile as our a posteriori tests
confirm that this is sufficient to capture the inversion layer accurately.
Higher-order terms could be incorporated, but this is not considered here to
keep the obtained profiles relatively simple. An important observation is that
the perturbation method approach is consistent with the finding of Garcia and
Mellado (2014). They showed that the vertical structure of the entrainment
zone is best described by two overlapping sublayers characterized by different
length scales, namely the mean penetration depth of an overshooting thermal
for the upper sublayer and the thickness of the CBL for the lower sublayer.
Similarly, the second-order ODE, i.e. Eq. (6), indicates that there are two
distinct length scales for the description of the entrainment zone (figure 1):
one is the upper sublayer with $z_{i}\leq z\leq h_{2}$, where the gradient of
potential temperature flux is proportional to $-(q_{w}\Pi_{m})/(2\epsilon
h_{2})$, and the other is the lower sublayer with $h_{1}\leq z\leq z_{i}$,
where the gradient of potential temperature flux is proportional to
$(q_{w}\Pi_{m})/z_{i}$. Since $\epsilon\ll 1$ the potential temperature flux
varies much steeper in the upper sublayer than in the lower sublayer.

###  2.2 Wind profile

We consider MOST to describe the wind speed profile in the surface layer
(Monin and Obukhov, 1954). In surface-layer coordinates, it states that the
non-dimensional streamwise velocity gradient can be written as

$\frac{\kappa
z}{u_{*}}\frac{{\textrm{d}}U}{{\textrm{d}}z}=\phi_{m}\left(\frac{z}{L}\right),$
(12)

where $\phi_{m}$ is the dimensionless stability correction function and
$L=-u_{*}^{3}/(\kappa\beta q_{w})$ is the Obukhov length with $\beta$ the
buoyancy parameter. By integrating Eq. (12), one can obtain the explicit
formula for the streamwise velocity $U$,

$\frac{\kappa
U}{u_{*}}=\ln\left(\frac{z}{z_{0}}\right)-\psi_{m}\left(\frac{z}{L}\right).$
(13)

Here $\kappa=0.4$ is the von Kármán constant, $z_{0}$ is the roughness length,
and

$\psi_{m}=\int_{z_{0}/L}^{z/L}\frac{1-\phi_{m}(\zeta)}{\zeta}{\textrm{d}}\zeta$
(14)

is the stability correction function for the momentum. We use the well-known
Businger-Dyer expression (Paulson, 1970; Businger et al., 1971; Dyer, 1974;
Brutsaert, 1982)

$\psi_{m}=\ln\frac{(1+x^{2})(1+x)^{2}}{8}-2\arctan x+\frac{\pi}{2},$ (15)

with $x=(1-16z/L)^{1/4}$ to model the stability correction function, but we
note that other parameterizations exist (Katul et al., 2011). Note that
$\psi_{m}=0$ for $x=1$ (or $L=\infty$), which reduces Eq. (13) to the
classical logarithmic law for neutral boundary layers.

To model the wind profile $U$ in the mixed layer and entrainment zone, we
again employ a second-order ODE. In general, the detailed wind profile evolves
when stability changes. However, the variations of the mean velocity are small
compared to the magnitude of the mean velocity in the mixed layer when the
stability parameter $-z_{i}/L\gtrsim 10$ (e.g. Salesky et al., 2017), which
covers the range of stability conditions considered in this study. Thus, we
can assume that the ODE is dominated by the $U^{\prime}=0$ term in most of the
domain. Recently, Tong and Ding (2020) derived the convective logarithmic
friction law from first principles, which connects the mixed-layer mean
velocity scale $U_{m}$ and the friction velocity $u_{*}$ in the convective-
roll dominant regime ($-z_{i}/L\gg 1$) as follows,

$\frac{U_{m}}{u_{*}}=\frac{1}{\kappa}\ln\left(-\frac{L}{z_{0}}\right)-C,$ (16)

where $C=1$ is an empirical constant determined from our LES database (see
below).

From the potential temperature flux profile modeling we learned that the ODE
should have a second-order derivative term $\epsilon U^{\prime\prime}$ to
model the entrainment zone near the top of the boundary layer. The top
boundary condition is given by the geostrophic wind component $U_{g}$. The
lower boundary condition is given by equaling Eqs. (13) and (16), namely
$U(\xi_{0})=U_{m}$, since Tong and Ding (2020) showed that $U_{m}$ is very
close to $U(z)$ in the mixed layer. Here $\xi_{0}$ represents the height of
the top of surface layer, which can be determined using Eqs. (13) and (16),

$\ln\left(-\frac{h_{2}}{L}\xi_{0}\right)-\psi_{m}\left(\frac{h_{2}}{L}\xi_{0}\right)=-\kappa
C\quad\Rightarrow\quad\xi_{0}=\xi_{0}\left(\frac{h_{2}}{L}\right).$ (17)

Because $\epsilon\ll 1$, the solution obtained from $U(\xi_{0})=U_{m}$ is
almost the same as from $U(0)=U_{m}$, while the expression of the latter is
much simpler. Therefore, we model the profile of the streamwise velocity $U$
in the mixed layer and the entrainment zone as:

$\epsilon U^{\prime\prime}-U^{\prime}=0,\quad U(0)=U_{m},\quad U(1)=U_{g}.$
(18)

The solution of Eq. (18) is

$U=U_{m}+(U_{g}-U_{m})\frac{e^{\xi/\epsilon}-1}{e^{1/\epsilon}-1}.$ (19)

We note that Eq. (19) is only valid in the mixed layer and entrainment zone as
the wind speed in the surface layer is still modeled using the MOST. By
combining Eqs. (13) and (19) and recalling that Eq. (13) increases
monotonically as $z$ increases, we obtain the following analytic description
of the streamwise velocity profile $U(z)$ for the entire CBL,

$U=\left\\{\begin{split}&\frac{u_{*}}{\kappa}\left[\ln\left(\frac{z}{z_{0}}\right)-\psi_{m}\left(\frac{z}{L}\right)\right],\quad&&\xi\leq\xi_{0},\\\
&U_{m}+(U_{g}-U_{m})\frac{e^{\xi/\epsilon}-1}{e^{1/\epsilon}-1},\quad&&\xi_{0}<\xi\leq
1,\end{split}\right.$ (20)

where $\xi_{0}$ is given by Eq. (17). As remarked in Section 1 the surface
layer profile contains two length scales, i.e. $z_{0}$ for the inner-inner
layer and the Obukhov length $L$ for the inner-outer layer. Similarly, the
velocity profile for the mixed layer and entrainment zone contains two length
scales, i.e. $z_{i}$ to describe the mixed layer and $h_{2}-z_{i}=2\epsilon
h_{2}$ to describe the upper sublayer of the entrainment zone. This confirms
the view presented by Tong and Ding (2020) that the entrainment zone has a
different scaling than the surface and mixed layers, and can therefore be
considered as another inner layer in the overall CBL problem. We note that the
proposed analytical profile is empirical, similar to the MOST, and that the
parameter $\epsilon$ parameterizes the effect of various physical processes.
We further note that $U_{m}$ and $u_{*}$ are related as given by Eq. (16), and
that the difference $U_{g}-U_{m}$ scales as $(u_{*}^{2}w_{e})/(fz_{i})^{2}$
(Tong and Ding, 2020). Thus, Eq. (20) is predictive if the entrainment
velocity $w_{e}$ is given as an input parameter. To determine the value of
$w_{e}$, one may need to revisit the entrainment processes at the top of CBLs
(e.g. Garcia and Mellado, 2014). In addition, the velocity predicted by Eq.
(20) is continuous throughout the boundary layer and applicable for the
considered ranges (see figure 5 below). However, its first derivative is
discontinuous at the patching location $\xi=\xi_{0}$. This is a typical
character of low-order models (Garratt et al., 1982). To capture the smooth
transition, a high-order model is needed (Tong and Ding, 2020). We leave these
for future work.

To model the wind profile $V$ in the the mixed layer and entrainment zone, we
use a similar ODE as Eq. (18). The top boundary condition is given by the
geostrophic wind component $V_{g}$. As the spanwise velocity $V$ is small
compared to the streamwise velocity $U$ in the mixed layer (Tong and Ding,
2020), the lower boundary condition is given by $V(\xi_{0})=V(0)=0$.
Therefore, we model the profile of the spanwise velocity $V$ in the entire
boundary layer using

$\epsilon V^{\prime\prime}-V^{\prime}=0,\quad V(0)=0,\quad V(1)=V_{g}.$ (21)

The solution of Eq. (21) is

$V=V_{g}\frac{e^{\xi/\epsilon}-1}{e^{1/\epsilon}-1}.$ (22)

Since $V_{g}$ scales as $-u_{*}^{2}/(fz_{i})$ (e.g. Wyngaard, 2010; Tong and
Ding, 2020), the geostrophic wind component $V_{g}$ can be connected to
$u_{*}$, up to a non-dimensional coefficient $-V_{g}fz_{i}/u_{*}^{2}=0.66$,
which is determined from our LES database (see Table 1).

##  3 Numerical validation

###  3.1 Numerical method and computational setup

We use LES to simulate the CBL flow over an infinite flat surface with
homogeneous roughness. We integrate the spatially-filtered Navier-Stokes
equations and the filtered transport equation for the potential temperature
(Albertson, 1996; Albertson and Parlange, 1999; Gadde et al., 2021; Liu et
al., 2021a, b; Liu and Stevens, 2021). Molecular viscosity is neglected as the
Reynolds number in the atmospheric boundary layer flow is very high, and we
use the advanced Lagrangian-averaging scale-dependent model to parameterize
the sub-grid scale shear stress and potential temperature flux (Bou-Zeid et
al., 2005; Stoll and Porté-Agel, 2008). We note that the Lagrangian-averaging
scale-dependent model has been extensively validated and widely used in the
literature (Bou-Zeid et al., 2005; Stoll and Porté-Agel, 2008; Calaf et al.,
2010; Wu and Porté-Agel, 2011; Zhang et al., 2019; Gadde et al., 2021).

Our code is an updated version of the one used by Albertson and Parlange
(1999). The grid points are uniformly distributed, and the computational
planes for horizontal and vertical velocities are staggered in the vertical
direction. The first vertical velocity grid plane is located at the ground.
The first gridpoint for the horizontal velocity components and the potential
temperature is located at half a grid distance above the ground. We use a
second-order finite difference method in the vertical direction and a pseudo-
spectral discretization in the horizontal directions. Time integration is
performed using the second-order Adams-Bashforth method. The projection method
is used to enforce the divergence-free condition. At the top boundary, we
impose a constant potential temperature lapse rate, zero vertical velocity,
and zero shear stress boundary condition. At the bottom boundary, we employ
the classical wall stress and potential temperature flux formulations based on
the MOST (Moeng, 1984; Bou-Zeid et al., 2005; Stoll and Porté-Agel, 2008;
Gadde et al., 2021).

We perform eleven LES to verify the validity of the derived wind speed and
potential temperature flux profiles for CBLs. The computational domain is
$5\,\textrm{km}\times 5\,\textrm{km}\times 2\,\textrm{km}$ and the grid
resolution is $256\times 256\times 256$. Due to large computational expense,
only several external parameters are varied in the simulations. The flow is
driven by the geostrophic wind of $G=\sqrt{U_{g}^{2}+V_{g}^{2}}=10$ m/s, the
buoyancy parameter is
$\beta=0.0325\,\textrm{m}/(\textrm{s}^{2}\cdot\textrm{K})$, and the Coriolis
parameter is $f=1\times 10^{-4}\,\textrm{rad/s}$ (Moeng and Sullivan, 1994;
Abkar and Moin, 2017; Gadde et al., 2021). To ensure the CBLs are in the
convective-roll dominant regime with $-z_{i}/L\gtrsim 10$, the surface
potential temperature flux is set to $q_{w}=0.12\sim
0.24\,\textrm{K}\cdot\textrm{m/s}$. Note that the convective logarithmic
friction law (Eq. (16)) is derived very recently by Tong and Ding (2020) and
tested only in a relatively narrow range of $-L/z_{0}$, namely
$-L/z_{0}\in[2.5\times 10^{2},1.5\times 10^{3}]$. To evaluate the performance
of this law in much wider range, i.e. $-L/z_{0}\in[3.6\times 10^{2},0.7\times
10^{5}]$, the roughness length is varied between $z_{0}=0.0002~{}\textrm{m}$
and $z_{0}=0.16~{}\textrm{m}$, where the lower bound of $z_{0}$ is set to a
representative value of the sea surface (Wieringa et al., 2001). The vertical
potential temperature gradient is varied between $\Gamma=1$ K/km and
$\Gamma=9$ K/km to capture the relevant range observed in atmospheric
measurements (Sorbjan, 1996). The velocity field is initialized with the
geostrophic wind $G=10$ m/s. The initial potential temperature is 300 K up to
937 m and increases with 8 K in the next 126 m above. Above $1063$ m the
constant vertical derivative of the potential temperature $\Gamma$ is
specified. The simulations are run for about 25 large-eddy turnover times
$T=z_{i}/w_{*}$, where $w_{*}=(\beta q_{w}z_{i})^{1/3}$ is the convective
velocity scale, and the statistics are computed from the time interval of
$12T$ to $25T$ when the boundary layer is quasi-stationary (Ding and Tong,
2021). We note that inertial oscillation develops as the flow is initialized
with a profile in geostrophic equilibrium (Schröter et al., 2013). However, it
should have only negligible effect on the simulated statistics since the
typical large-eddy turnover time $T\approx 10$ mins is two order smaller than
the natural inertial periodicity $2\pi/f\approx 17.5$ h.

A summary of all simulated cases is presented in Table 1. Note that the cases
in Table 1 are arranged such that the value of $-L/z_{0}$ increases
monotonically. Furthermore, we note that case 2 has been validated against
atmospheric observations, and the simulation results obtained using different
sub-grid scale models and grid resolutions is very similar (Gadde et al.,
2021). To show the simulated results are independent of the computational
domain size, we have performed an additional simulation for case 2 in a larger
computational domain ($12\,\textrm{km}\times 6\,\textrm{km}\times
2\,\textrm{km}$) on a mesh with $600\times 300\times 240$ nodes such that the
grid spacings are nearly identical. Figure 2 shows the simulated wind speed
and potential temperature flux profiles for case 2. The good agreement between
the results obtained with different computational domain size confirms that
the simulated results are independent of the computational domain size.

Table 1: Summary of all simulated cases. Here $\Gamma=\partial\Theta/\partial z$ is the vertical derivative of the mean potential temperature in the free atmosphere, $q_{w}$ is the surface potential temperature flux, $z_{0}$ is the roughness length, $u_{*}$ is the friction velocity, $U_{m}$ is the wind speed in the mixed layer, $L$ is the Obukhov length, $z_{i}$ is the inversion layer height, $|V_{g}|$ is the magnitude of the spanwise geostrophic wind, $\epsilon$ and $c_{\Pi}$ are dimensionless parameters calculated by Eqs. (5) and (9), respectively, and $Ri=\beta\Delta\Theta z_{i}/w_{*}^{2}$ is the Richardson number, where $\Delta\Theta=\Theta(h_{2})-\Theta(h_{1})$ is the potential temperature difference across the entrainment zone and $w_{*}=(\beta q_{w}z_{i})^{1/3}$ is the convective velocity. Case |  $\Gamma$ (K/km) |  $q_{w}$ (K$\cdot$m/s) |  $z_{0}$ (m) |  $u_{*}$ (m/s) |  $U_{m}$ (m/s) |  $|V_{g}|$ (m/s) | $\epsilon$ | $c_{\Pi}$ | $Ri$ | $-z_{i}/L$ | $-L/z_{0}$
---|---|---|---|---|---|---|---|---|---|---|---
1 | 9 | 0.24 | 0.16 | 0.562 | 7.60 | 2.00 | 0.044 | 1.32 | 56.1 | 19.2 | $3.6\times 10^{2}$
2 | 3 | 0.24 | 0.16 | 0.563 | 7.59 | 1.87 | 0.052 | 1.34 | 51.0 | 19.1 | $3.6\times 10^{2}$
3 | 1 | 0.24 | 0.16 | 0.562 | 7.59 | 1.84 | 0.055 | 1.34 | 47.9 | 19.2 | $3.6\times 10^{2}$
4 | 3 | 0.12 | 0.16 | 0.533 | 7.70 | 2.20 | 0.046 | 1.34 | 94.2 | 11.0 | $0.6\times 10^{3}$
5 | 3 | 0.24 | 0.016 | 0.463 | 8.44 | 1.17 | 0.050 | 1.33 | 51.0 | 34.5 | $2.0\times 10^{3}$
6 | 3 | 0.20 | 0.02 | 0.468 | 8.36 | 1.32 | 0.046 | 1.31 | 59.6 | 27.5 | $2.0\times 10^{3}$
7 | 3 | 0.12 | 0.016 | 0.444 | 8.45 | 1.43 | 0.038 | 1.31 | 91.5 | 19.0 | $3.5\times 10^{3}$
8 | 3 | 0.20 | 0.002 | 0.392 | 8.93 | 0.86 | 0.033 | 1.28 | 55.5 | 47.3 | $1.2\times 10^{4}$
9 | 3 | 0.24 | 0.0016 | 0.389 | 8.94 | 0.75 | 0.044 | 1.30 | 50.4 | 58.3 | $1.2\times 10^{4}$
10 | 3 | 0.12 | 0.0016 | 0.375 | 8.94 | 0.88 | 0.036 | 1.30 | 93.9 | 31.6 | $2.1\times 10^{4}$
11 | 3 | 0.20 | 0.0002 | 0.334 | 9.24 | 0.57 | 0.041 | 1.30 | 58.6 | 75.8 | $0.7\times 10^{5}$
|  |  |  |  |  |  |  |  |  |  |
\begin{overpic}[width=325.215pt]{fig2} \end{overpic} Figure 2: The comparison
of simulated (a) normalized wind speed $\sqrt{U^{2}+V^{2}}/G$ and (b)
normalized potential temperature flux $q/q_{w}$ profiles for case 2 (see Table
1) with different computational domain sizes.

###  3.2 Validation of analytical profiles

Figure 3 shows that with increasing height the normalized potential
temperature flux profile $q/q_{w}$ first decreases linearly from unity at the
surface to a minimum at $z_{i}/h_{2}=1-2\epsilon$, before it rapidly increases
to zero in the inversion layer ($1-2\epsilon\leq z/h_{2}\leq 1$). The
normalized thickness of the inversion layer, which is parameterized by
$\epsilon$, is expected to depend on the Richardson number (Deardorff et al.,
1980), potential temperature gradient (Sorbjan, 1996), and wind shear
(Conzemius and Fedorovich, 2006). However, we find that for the parameter
range under consideration, the variation in the normalized thickness of the
inversion layer is limited (see Table 1). Therefore, we use a fixed
representative value $\epsilon=0.044$ to model the potential temperature flux
profile, and the figure confirms that this ensures that the potential
temperature flux profile obtained from the model agrees excellently with all
available simulation data, which validates the chosen approach. To further
confirm the validity of the potential temperature flux profile, we also
compare our results in figure 3 with previous LES from Mason (1989), Sorbjan
(1996), and Abkar and Moin (2017), the direct numerical simulations data by
Garcia and Mellado (2014), and the empirical models by Lenschow (1974) and Noh
et al. (2003). Clearly, the model predictions agree well with these previous
studies.

\begin{overpic}[width=216.81pt]{fig3} \end{overpic} Figure 3: Vertical profile
of normalized potential temperature flux $q/q_{w}$. Circles: LES data (Table
1); up-triangle: LES data of Mason (1989); down-triangle: LES data of Sorbjan
(1996); square: direct numerical simulations data of Garcia and Mellado
(2014); diamond: LES data of Abkar and Moin (2017); red line: prediction given
by Lenschow (1974); blue line: prediction given by Noh et al. (2003); black
line: prediction given by Eq. (7) with $\epsilon=0.044$ and $c_{\Pi}=1.32$.
The figure shows that the proposed model captures the simulation trends very
well. \begin{overpic}[width=325.215pt]{fig4} \end{overpic} Figure 4: The
dependence of (a) the normalized mixed-layer mean velocity scale $U_{m}/u_{*}$
and (b) the ratio of the predicted and simulated mixed-layer mean velocity
scale $U_{m}^{\rm pred}/U_{m}^{\rm LES}$ against the normalized Obukhov length
$-L/z_{0}$. Circles: LES data (Table 1); squares: LES data of Tong and Ding
(2020); solid line: prediction given by Eq. (16) with $C=1$ determined using a
fit to the present simulation data.  \begin{overpic}[width=325.215pt]{fig5}
\end{overpic} Figure 5: Vertical profile of the mean streamwise velocity $U$.
Filled circles: LES data (Table 1); dashed line: prediction given by the MOST;
solid line: prediction given by Eq. (20) with $\epsilon=0.044$. The prediction
by the MOST is plotted outside the surface layer region to demonstrate the
difference with the new profile.

In figure 4(a) we compare the wind speed in the bulk of the mixed layer
($0.4\leq z/h_{2}\leq 0.6$) against the normalized Obukhov length $-L/z_{0}$
with results from Tong and Ding (2020). The figure shows that our simulations
convincingly confirm the validity of the convective logarithmic friction law
for the wind speed (Eq. (16) with $C=1$) over a much wider range of
$-L/z_{0}\in[3.6\times 10^{2},0.7\times 10^{5}]$ than previously considered
($-L/z_{0}\in[2.5\times 10^{2},1.5\times 10^{3}]$). To further confirm the
convective logarithmic friction law, figure 4(b) shows the ratio of the
predicted and simulated mixed-layer mean velocity scale $U_{m}^{\rm
pred}/U_{m}^{\rm LES}$ against the normalized Obukhov length $-L/z_{0}$. The
figure shows that the predicted value is within 5% of the simulation result
for all considered cases.

\begin{overpic}[width=325.215pt]{fig6} \end{overpic} Figure 6: Vertical
profile of the mean spanwise velocity $V$. Filled circles: LES data (Table 1);
solid line: prediction given by Eq. (22) with $\epsilon=0.044$.

In figure 5 we compare the vertical profile of the mean streamwise velocity
$U$ for four typical cases with different surface potential temperature flux
and roughness length with the simulation results. The filled symbols are the
present LES data, the dashed line is the theoretical prediction given by the
MOST, and the solid line is the prediction given by Eq. (20) with
$\epsilon=0.044$. The figure shows that the MOST accurately captures the
surface layer’s wind profile (lowest $20\%$ of the boundary layer). However,
in the mixed layer, the prediction of the MOST deviates significantly from the
LES data. In particular, the discrepancy from the MOST increases as the
surface potential temperature flux $q_{w}$ decreases (figure 5b,c) or the
roughness length $z_{0}$ increases (figure 5a,d). Therefore, MOST is seldom
used to specify wind profiles outside of the surface layer. Figure 5 shows
that the proposed wind profile given by Eq. (20) accurately captures the
velocity profile throughout the entire boundary layer. This excellent
agreement confirms the validity of our proposed wind profile of Eq. (20) for
atmospheric boundary layers in the range of studied parameters.

Figure 6 shows the corresponding profiles of the mean spanwise velocity $V$.
The filled symbols are the present LES data and the solid line is the
prediction given by Eq. (22) with $\epsilon=0.044$. Overall, the agreement
between the proposed wind profile given by Eq. (22) and the LES data is
reasonably good in the entire boundary layer. This agreement confirms the
validity of our proposed wind profile of Eq. (22) for CBLs in the range of
studied parameters (i.e. $-L/z_{0}\in[3.6\times 10^{2},0.7\times 10^{5}]$). We
note that the figure confirms that the spanwise velocity $V$ is much smaller
than the streamwise velocity $U$. The figure also indicates that the magnitude
of the geostrophic wind component $|V_{g}|$ increases as the surface potential
temperature flux $q_{w}$ decreases (figure 6b,c) or the roughness length
$z_{0}$ increases (figure 6a,d).

##  4 Conclusions

This work uses a perturbation method approach in conjuncture with the
convective logarithmic friction law and the Monin-Obukhov similarity theory to
develop analytical expressions of the wind and potential temperature flux
profiles in convective atmospheric boundary layers. The validity of the
proposed wind (given by Eqs. (20) and (22)) and potential temperature flux
profiles (given by Eq. (7)) has been confirmed by their excellent agreement
with large-eddy simulations results for atmospheric boundary layers in the
convective-roll dominant regime with $-z_{i}/L\gtrsim 10$, where $L$ is the
Obukhov length and $z_{i}$ the inversion layer height. Furthermore, our
simulations confirm that the convective logarithmic friction law of Eq. (16),
which was originally proposed by Tong and Ding (2020) for the mixed-layer mean
velocity scale, is valid for an extensive range of $-L/z_{0}$, namely
$-L/z_{0}\in[3.6\times 10^{2},0.7\times 10^{5}]$, where $z_{0}$ is the surface
roughness length. Since accurate capturing the coupling between meso and
microscale processes is a long-standing challenge in numerical weather
predictions (Wyngaard, 2004; Larsén et al., 2018; Veers et al., 2019), the
proposed analytical profiles may be relevant for climate modeling and weather
forecasting to better understand the effect of convective atmospheric boundary
layers on, for example, wind farms. Possible future work will involve
investigating models to predict the entrainment velocity at the top of CBLs
and developing a high-order model that can capture the transition between the
entrainment zone and free atmosphere. The latter may require a formal
asymptotic series expansion of the governing equations, allowing for the
separation into a time-dependent and steady-state problem at different orders.

###### Acknowledgements.

This work was supported by the Hundred Talents Program of the Chinese Academy
of Sciences, the National Natural Science Fund for Excellent Young Scientists
Fund Program (Overseas), the National Natural Science Foundation of China
Grant (No. 11621202), the Shell-NWO/FOM-initiative Computational sciences for
energy research of Shell and Chemical Sciences, Earth and Live Sciences,
Physical Sciences, Stichting voor Fundamenteel Onderzoek der Materie (FOM) and
STW, and an STW VIDI Grant (No. 14868). This work was sponsored by NWO Domain
Science for the use of the national computer facilities. We acknowledge PRACE
for awarding us access to Irene at Très Grand Centre de Calcul du CEA (TGCC)
under PRACE project 2019215098, and the advanced computing resources provided
by the Supercomputing Center of the USTC. The data that support the findings
of this study are available from the corresponding author upon reasonable
request.

## References

  * Abkar and Moin (2017) Abkar, M., and P. Moin, 2017: Large eddy simulation of thermally stratified atmospheric boundary layer flow using a minimum dissipation model. Boundary-Layer Meteorol., 165 (3), 405–419, 10.1007/s10546-017-0288-4.
  * Albertson (1996) Albertson, J. D., 1996: Large eddy simulation of land-atmosphere interaction. Ph.D. thesis, University of California.
  * Albertson and Parlange (1999) Albertson, J. D., and M. B. Parlange, 1999: Surface length-scales and shear stress: implications for land-atmosphere interaction over complex terrain. Water Resour. Res., 35, 2121–2132, 10.1029/1999WR900094.
  * Arya (1975) Arya, S. P. S., 1975: Comments on “similarity theory for the planetary boundary layer of time-dependent height”. J. Atmos. Sci., 32 (4), 839–840, 10.1175/1520-0469(1975)032¡0839:COTFTP¿2.0.CO;2.
  * Blackadar and Tennekes (1968) Blackadar, A. K., and H. Tennekes, 1968: Asymptotic similarity in neutral barotropic planetary boundary layers. J. Atmos. Sci., 25 (6), 1015–1020, 10.1175/1520-0469(1968)025¡1015:ASINBP¿2.0.CO;2.
  * Bou-Zeid et al. (2005) Bou-Zeid, E., C. Meneveau, and M. B. Parlange, 2005: A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows. Phys. Fluids, 17, 025 105, 10.1063/1.1839152.
  * Brutsaert (1982) Brutsaert, W., 1982: Evaporation into the Atmosphere: Theory, History and Applications. Springer, Netherlands.
  * Businger et al. (1971) Businger, J. A., J. C. Wyngaard, Y. Izumi, and E. F. Bradley, 1971: Flux-profile relationships in the atmospheric surface layer. J. Atmos. Sci., 28 (2), 181–189, 10.1175/1520-0469(1971)028¡0181:FPRITA¿2.0.CO;2.
  * Calaf et al. (2010) Calaf, M., C. Meneveau, and J. Meyers, 2010: Large eddy simulations of fully developed wind-turbine array boundary layers. Phys. Fluids, 22, 015 110, 10.1063/1.3291077.
  * Cheng et al. (2021) Cheng, Y., Q. Li, D. Li, and P. Gentine, 2021: Logarithmic profile of temperature in sheared and unstably stratified atmospheric boundary layers. Phys. Rev. Fluids, 6, 034 606, 10.1103/PhysRevFluids.6.034606.
  * Cheng et al. (2005) Cheng, Y., M. B. Parlange, and W. Brutsaert, 2005: Pathology of Monin-Obukhov similarity in the stable boundary layer. J. Geophys. Res., 110, D06 101, 10.1029/2004JD004923.
  * Conzemius and Fedorovich (2006) Conzemius, R. J., and E. Fedorovich, 2006: Dynamics of sheared convective boundary layer entrainment. Part I: Methodological background and large-eddy simulations. J. Atmos. Sci., 63 (4), 1151–1178, 10.1175/jas3691.1.
  * Deardorff (1973) Deardorff, J. W., 1973: An explanation of anomalously large Reynolds stresses within the convective planetary boundary layer. J. Atmos. Sci., 30 (6), 1070–1076, 10.1175/1520-0469(1973)030¡1070:Aeoalr¿2.0.Co;2.
  * Deardorff (1979) Deardorff, J. W., 1979: Prediction of convective mixed-layer entrainment for realistic capping inversion structure. J. Atmos. Sci., 36 (3), 424–436, 10.1175/1520-0469(1979)036¡0424:Pocmle¿2.0.Co;2.
  * Deardorff et al. (1980) Deardorff, J. W., G. E. Willis, and B. H. Stockton, 1980: Laboratory studies of the entrainment zone of a convectively mixed layer. J. Fluid Mech., 100, 41–64, 10.1017/S0022112080001000.
  * Ding and Tong (2021) Ding, M., and C. Tong, 2021: Multi-point Monin-Obukhov similarity of turbulence cospectra in the convective atmospheric boundary layer. Boundary-Layer Meteorol., 178, 185–199, 10.1007/s10546-020-00571-4.
  * Dyer (1974) Dyer, A. J., 1974: A review of flux-profile relationships. Boundary-Layer Meteorol., 7 (3), 363–372, 10.1007/BF00240838.
  * Gadde et al. (2021) Gadde, S. N., A. Stieren, and R. J. A. M. Stevens, 2021: Large-eddy simulations of stratified atmospheric boundary layers: Comparison of different subgrid models. Boundary-Layer Meteorol., 178 (3), 363–382, 10.1007/s10546-020-00570-5.
  * Garcia and Mellado (2014) Garcia, J. R., and J. P. Mellado, 2014: The two-layer structure of the entrainment zone in the convective boundary layer. J. Atmos. Sci., 71 (6), 1935–1955, 10.1175/jas-d-13-0148.1.
  * Garratt et al. (1982) Garratt, J. R., J. C. Wyngaard, and R. J. Francey, 1982: Winds in the atmospheric boundary layer-prediction and observation. J. Atmos. Sci., 39, 1307–1316, 10.1175/1520-0469(1982)039¡1307:Witabl¿2.0.Co;2.
  * Haghshenas and Mellado (2019) Haghshenas, A., and J. P. Mellado, 2019: Characterization of wind-shear effects on entrainment in a convective boundary layer. J. Fluid Mech., 858, 145–183, 10.1017/jfm.2018.761.
  * Holtslag and Moeng (1991) Holtslag, A. A. M., and C.-H. Moeng, 1991: Eddy diffusivity and countergradient transport in the convective atmospheric boundary layer. J. Atmos. Sci., 48, 1690–1698, 10.1175/1520-0469(1991)048¡1690:Edacti¿2.0.Co;2.
  * Hong et al. (2006) Hong, S.-Y., Y. Noh, and J. Dudhia, 2006: A new vertical diffusion package with an explicit treatment of entrainment processes. Mon. Weather Rev., 134 (9), 2318–2341, 10.1175/mwr3199.1.
  * Johansson et al. (2001) Johansson, C., A.-S. Smedman, U. Högström, J. G. Brasseur, and S. Khanna, 2001: Critical test of the validity of Monin-Obukhov similarity during convective conditions. J. Atmos. Sci., 58, 1549–1566, 10.1175/1520-0469(2001)058¡1549:Ctotvo¿2.0.Co;2.
  * Kaimal et al. (1976) Kaimal, J. C., J. C. Wyngaard, D. A. Haugen, O. R. Coté, Y. Izumi, S. J. Caughey, and C. J. Readings, 1976: Turbulence structure in the convective boundary layer. J. Atmos. Sci., 33 (11), 2152–2169, 10.1175/1520-0469(1976)033¡2152:TSITCB¿2.0.CO;2.
  * Katul et al. (2011) Katul, G. G., A. G. Konings, and A. Porporato, 2011: Mean velocity profile in a sheared and thermally stratified atmospheric boundary layer. Phys. Rev. Lett., 107 (26), 268 502, 10.1103/PhysRevLett.107.268502.
  * Khanna and Brasseur (1997) Khanna, S., and J. G. Brasseur, 1997: Analysis of Monin-Obukhov similarity from large-eddy simulation. J. Fluid Mech., 345, 251–286, 10.1017/S0022112097006277.
  * Larsén et al. (2018) Larsén, X. G., E. L. Petersen, and S. E. Larsen, 2018: Variation of boundary‐layer wind spectra with height. Q. J. R. Meteorol. Soc., 144, 2054–2066, 10.1002/qj.3301.
  * LeMone et al. (2019) LeMone, M. A., and Coauthors, 2019: 100 years of progress in boundary layer meteorology. Meteorological Monographs, 59, 9.1–9.85, 10.1175/amsmonographs-d-18-0013.1.
  * Lenschow (1974) Lenschow, D. H., 1974: Model of the height variation of the turbulence kinetic energy budget in the unstable planetary boundary layer. J. Atmos. Sci., 31 (2), 465–474, 10.1175/1520-0469(1974)031¡0465:Mothvo¿2.0.Co;2.
  * Li et al. (2021) Li, Q., Y. Cheng, and P. Gentine, 2021: Connection between mass flux transport and eddy diffusivity in convective atmospheric boundary layers. Geophys. Res. Lett., 48, e2020GL092 073, 10.1029/2020gl092073.
  * Lilly (1968) Lilly, D. K., 1968: Models of cloud-topped mixed layers under a strong inversion. Q. J. R. Meteorol. Soc., 94 (401), 292–309, 10.1002/qj.49709440106.
  * Liu et al. (2021a) Liu, L., S. N. Gadde, and R. J. A. M. Stevens, 2021a: Geostrophic drag law for conventionally neutral atmospheric boundary layers revisited. Q. J. R. Meteorol. Soc., 147, 847–857, 10.1002/qj.3949.
  * Liu et al. (2021b) Liu, L., S. N. Gadde, and R. J. A. M. Stevens, 2021b: Universal wind profile for conventionally neutral atmospheric boundary layers. Phys. Rev. Lett., 126, 104 502, 10.1103/PhysRevLett.126.104502.
  * Liu and Stevens (2021) Liu, L., and R. J. A. M. Stevens, 2021: Effects of atmospheric stability on the performance of a wind turbine located behind a three-dimensional hill. Renewable Energy, 175, 926–935, 10.1016/j.renene.2021.05.035.
  * Mahrt et al. (1998) Mahrt, L., J. Sun, W. Blumen, T. Delany, and S. Oncley, 1998: Nocturnal boundary-layer regimes. Boundary-Layer Meteorol., 88 (2), 255–278, 10.1023/A:1001171313493.
  * Mason (1989) Mason, P. J., 1989: Large-eddy simulation of the convective atmospheric boundary layer. J. Atmos. Sci., 46, 1492–1516, 10.1175/1520-0469(1989)046¡1492:LESOTC¿2.0.CO;2.
  * McNaughton et al. (2007) McNaughton, K. G., R. J. Clement, and J. B. Moncrieff, 2007: Scaling properties of velocity and temperature spectra above the surface friction layer in a convective atmospheric boundary layer. Nonlin. Processes Geophys., 14, 257–271, 10.5194/npg-14-257-2007.
  * Moeng (1984) Moeng, C.-H., 1984: A large-eddy simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052–2062, 10.1175/1520-0469(1984)041¡2052:ALESMF¿2.0.CO;2.
  * Moeng and Sullivan (1994) Moeng, C.-H., and P. P. Sullivan, 1994: A comparison of shear- and buoyancy-driven planetary boundary layer flows. J. Atmos. Sci., 51, 999–1022, 10.1175/1520-0469(1994)051¡0999:ACOSAB¿2.0.CO;2.
  * Monin and Obukhov (1954) Monin, A. S., and A. M. Obukhov, 1954: Basic laws of turbulent mixing in the surface layer of the atmosphere. Tr. Akad. Nauk SSSR Geophiz. Inst., 24 (151), 163–187.
  * Nieuwstadt et al. (2016) Nieuwstadt, F. T. M., J. Westerweel, and B. J. Boersma, 2016: Turbulence: Introduction to Theory and Applications of Turbulent Flows. Springer, Switzerland.
  * Noh et al. (2003) Noh, Y., W. G. Cheon, S. Y. Hong, and S. Raasch, 2003: Improvement of the K-profile model for the planetary boundary layer based on large eddy simulation data. Boundary-Layer Meteorol., 107 (2), 401–427, 10.1023/A:1022146015946.
  * Obukhov (1946) Obukhov, A. M., 1946: Turbulence in an atmosphere with inhomogeneous temperature. Trans. Inst. Teoret. Geoz. Akad. Nauk SSSR, 1, 95–115.
  * Panofsky et al. (1977) Panofsky, H. A., H. Tennekes, D. H. Lenschow, and J. C. Wyngaard, 1977: The characteristics of turbulent velocity components in the surface layer under convective conditions. Boundary-Layer Meteorol., 11, 355–361, 10.1007/BF02186086.
  * Paulson (1970) Paulson, C. A., 1970: The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. J. Appl. Meteorol. Climatol., 9, 857–861, 10.1175/1520-0450(1970)009¡0857:TMROWS¿2.0.CO;2.
  * Rossby and Montgomery (1935) Rossby, C. G., and R. B. Montgomery, 1935: The layers of frictional influence in wind and ocean currents. Pap. Phys. Oceanogr. Meteor., 3 (3), 1–101, 10.1575/1912/1157.
  * Salesky and Anderson (2020) Salesky, S. T., and W. Anderson, 2020: Coherent structures modulate atmospheric surface layer flux-gradient relationships. Phys. Rev. Lett., 125, 124 501, 10.1103/PhysRevLett.125.124501.
  * Salesky et al. (2017) Salesky, S. T., M. Chamecki, and E. Bou-Zeid, 2017: On the nature of the transition between roll and cellular organization in the convective boundary layer. Boundary-Layer Meteorol., 163, 41–68, 10.1007/s10546-016-0220-3.
  * Schröter et al. (2013) Schröter, J. S., A. F. Moene, and A. A. M. Holtslag, 2013: Convective boundary layer wind dynamics and inertial oscillations: the influence of surface stress. Q. J. R. Meteorol. Soc., 139, 1694–1711, 10.1002/qj.2069.
  * Siebesma et al. (2007) Siebesma, A. P., P. M. M. Soares, and J. Teixeira, 2007: A combined eddy-diffusivity mass-flux approach for the convective boundary layer. J. Atmos. Sci., 64, 1230–1248, 10.1175/jas3888.1.
  * Sorbjan (1996) Sorbjan, Z., 1996: Effects caused by varying the strength of the capping inversion based on a large eddy simulation model of the shear-free convective boundary layer. J. Atmos. Sci., 53 (14), 2015–2024, 10.1175/1520-0469(1996)053¡2015:ECBVTS¿2.0.CO;2.
  * Stoll and Porté-Agel (2008) Stoll, R., and F. Porté-Agel, 2008: Large-eddy simulation of the stable atmospheric boundary layer using dynamic models with different averaging schemes. Boundary-Layer Meteorol., 126, 1–28, 10.1007/s10546-007-9207-4.
  * Stull (1976) Stull, R. B., 1976: The energetics of entrainment across a density interface. J. Atmos. Sci., 33 (7), 1260–1267, 10.1175/1520-0469(1976)033¡1260:Teoead¿2.0.Co;2.
  * Stull (1988) Stull, R. B., 1988: An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Boston.
  * Sullivan et al. (1998) Sullivan, P. P., C.-H. Moeng, B. Stevens, D. H. Lenschow, and S. D. Mayor, 1998: Structure of the entrainment zone capping the convective atmospheric boundary layer. J. Atmos. Sci., 55, 3042–3064, 10.1175/1520-0469(1998)055¡3042:SOTEZC¿2.0.CO;2.
  * Sullivan and Patton (2011) Sullivan, P. P., and E. G. Patton, 2011: The effect of mesh resolution on convective boundary layer statistics and structures generated by large-eddy simulation. J. Atmos. Sci., 68, 2395–2415, 10.1175/JAS-D-10-05010.1.
  * Sun and Wang (2008) Sun, J., and Y. Wang, 2008: Effect of the entrainment flux ratio on the relationship between entrainment rate and convective Richardson number. Boundary-Layer Meteorol., 126, 237–247, 10.1007/s10546-007-9231-4.
  * Tennekes and Driedonks (1981) Tennekes, H., and A. G. M. Driedonks, 1981: Basic entrainment equations for the atmospheric boundary layer. Boundary-Layer Meteorol., 20 (4), 515–531, 10.1007/BF00122299.
  * Tennekes and Lumley (1972) Tennekes, H., and J. L. Lumley, 1972: A First Course in Turbulence. The MIT Press, Cambridge, Massachusetts.
  * Tong and Ding (2020) Tong, C., and M. Ding, 2020: Velocity-defect laws, log law and logarithmic friction law in the convective atmospheric boundary layer. J. Fluid Mech., 883, A36, 10.1017/jfm.2019.898.
  * Troen and Mahrt (1986) Troen, I. B., and L. Mahrt, 1986: A simple model of the atmospheric boundary layer; sensitivity to surface evaporation. Boundary-Layer Meteorol., 37 (1), 129–148, 10.1007/BF00122760.
  * Veers et al. (2019) Veers, P., and Coauthors, 2019: Grand challenges in the science of wind energy. Science, 366 (6464), eaau2027, 10.1126/science.aau2027.
  * Wieringa et al. (2001) Wieringa, J., A. G. Davenport, C. S. B. Grimmond, and T. R. Oke, 2001: New revision of Davenport roughness classification. 3EACWE (3rd European $\&$ African Conference on Wind Engineering), Eindhoven, The Netherlands.
  * Wu and Porté-Agel (2011) Wu, Y.-T., and F. Porté-Agel, 2011: Large-eddy simulation of wind-turbine wakes: Evaluation of turbine parametrisations. Boundary-Layer Meteorol., 138, 345–366, 10.1007/s10546-010-9569-x.
  * Wyngaard (2004) Wyngaard, J. C., 2004: Toward numerical modeling in the “Terra Incognita”. J. Atmos. Sci., 61, 1816–1826, 10.1175/1520-0469(2004)061¡1816:TNMITT¿2.0.CO;2.
  * Wyngaard (2010) Wyngaard, J. C., 2010: Turbulence in the Atmosphere. Cambridge University Press, Cambridge.
  * Zhang et al. (2019) Zhang, M., M. G. Arendshorst, and R. J. A. M. Stevens, 2019: Large eddy simulations of the effect of vertical staggering in extended wind farms. Wind Energy, 22 (2), 189–204, 10.1002/we.2278.
  * Zilitinkevich (1969) Zilitinkevich, S. S., 1969: On the computation of the basic parameters of the interaction between the atmosphere and the ocean. Tellus, 21 (1), 17–24, 10.1111/j.2153-3490.1969.tb00414.x.
  * Zilitinkevich (1975) Zilitinkevich, S. S., 1975: Resistance laws and prediction equations for the depth of the planetary boundary layer. J. Atmos. Sci., 32 (4), 741–752, 10.1175/1520-0469(1975)032¡0741:Rlapef¿2.0.Co;2.
  * Zilitinkevich and Deardorff (1974) Zilitinkevich, S. S., and J. W. Deardorff, 1974: Similarity theory for the planetary boundary layer of time-dependent height. J. Atmos. Sci., 31 (5), 1449–1452, 10.1175/1520-0469(1974)031¡1449:Stftpb¿2.0.Co;2.
  * Zilitinkevich et al. (1992) Zilitinkevich, S. S., E. E. Fedorovich, and M. V. Shabalova, 1992: Numerical model of a non-steady atmospheric planetary boundary layer, based on similarity theory. Boundary-Layer Meteorol., 59 (4), 387–411, 10.1007/BF02215460.
  * Zilitinkevich et al. (2012) Zilitinkevich, S. S., S. A. Tyuryakov, Y. I. Troitskaya, and E. A. Mareev, 2012: Theoretical models of the height of the atmospheric boundary layer and turbulent entrainment at its upper boundary. Izv. Atmos. Ocean. Phys., 48 (1), 133–142, 10.1134/S0001433812010148.