# Evolution of the Well-Mixed Convective Boundary Layer Described by Eddy-Diffusion Formulations by Eddy-Diffusion Formulatoins

Jun-ichi Yano

► To cite this version:

Jun-ichi Yano. Evolution of the Well-Mixed Convective Boundary Layer Described by Eddy-Diffusion Formulations by Eddy-Diffusion Formulatoins. 2025. hal-05336937

HAL Id: hal-05336937

https://hal.science/hal-05336937v1

Preprint submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Evolution of the Well-Mixed Convective Boundary Layer Described by Eddy-Diffusion Formulations

by

Jun-Ichi Yano

CNRM, UMR 3589 (CNRS), Météo-France, 31057, Toulouse Cedex, France

jun-ichi.yano@cnrs.fr

## Summary

The basic-behaviors of the eddy-diffusion models under various closures in simulating the well-mixed convective boundary layer are examined. Some of the elaborations have inherent difficulties for general applications, such as indeterminate eddy diffusivity. Those are also often implicitly designed to produce an expected behavior of the well-mixed layer, for example, its control by the surface flux or constraining as leading to a nearly neutrally-stratified state as observed. On the other hand, the third-order closure leads to nonlocal transport.

## 1. Purpose

This study examines the capacity of the eddy-diffusion models to simulate the evolution of the well-mixed convective boundary layer predicted from the theory of well-mixed layer.

## 2. Problem Setting

We consider the simple dry atmosphere with the only heat source at the surface. We focus on the evolution of a domain-averaged potential temperature, $\langle\theta\rangle$, in which the domain is taken large enough so that the small-scale eddies controlling the evolution can be considered in quasi-equilibrium, but small enough so that the inhomogeneity in the large scale can be neglected. The evolution of $\langle\theta\rangle$ is dictated by:

$$ \frac{\partial}{\partial t}\langle\theta\rangle = -\frac{\partial}{\partial z}\langle w'\theta'\rangle \qquad (2.1) $$

Here, $w$ is the vertical velocity, and $\langle w'\theta'\rangle$ is the domain averaged vertical eddy heat flux, in which the variables with the prime designate the eddies defined as deviations from the domain averages. Here, $z$ is the vertical coordinate with the surface taken at $z = 0$.

The theory of well-mixed layer predicts that a layer with a constant potential temperature gradually develops with time. The main question to be posed here is how well the eddy-diffusion models can simulate this process.
For simplicity, we consider a system with a prescribed linear growth of the surface potential temperature, i.e.,

$$ \langle \theta \rangle \Big|_{z=0} = \left( \frac{\partial \theta}{\partial t} \right)_s t \qquad (2.2a) $$

with $(\partial\theta/\partial t)_s$ a constant. As it turns out, a more realistic model with a constant surface heat flux leads to a more involved analytical treatment without gaining extra insights as a whole as shown in the Appendix A.2.

We also assume that the atmosphere is under a constant lapse rate initially, i.e.,

$$ \langle \theta \rangle \Big|_{t=0} = \theta_0 + \left( \frac{\partial \bar{\theta}}{\partial z} \right) z \qquad (2.2b) $$

where $\theta_0$ and $(\partial\bar{\theta}/\partial z)$ are constants. Since the absolute value of the potential temperature does not contribute to the problem in the following, we simply set $\theta_0 = 0$ for simplicity in the graphics presentations.

Under those boundary conditions, the theory of well-mixed layer predicts that the system evolves as:

$$ \langle \theta \rangle = \begin{cases} \theta_0 + \Delta\theta & z \le z_i \\ \theta_0 + \left( \frac{\partial \bar{\theta}}{\partial z} \right) z & z > z_i \end{cases} \qquad (2) $$

where

$$ \Delta\theta = \left( \frac{\partial\theta}{\partial t} \right)_s t \qquad (2.3a) $$

is a constant with height, and

$$ z_i = \left( \frac{\partial \theta}{\partial t} \right)_s \left( \frac{\partial \bar{\theta}}{\partial z} \right)^{-1} t \qquad (2.3b) $$

is the height of the well-mixed layer. Furthermore, it may be worthwhile to note the eddy flux itself takes the form:

$$ \langle \theta' w' \rangle = \langle \theta' w' \rangle_0 \left( 1 - \frac{z}{z_i} \right) \qquad (2.3c) $$
<figure><img src="image_1.png" /><figcaption>Fig. 0: Evolution of the well-mixed convective boundary layer initiated by Eq. (2.2a) with the interval of $\Delta t = 0.1$ with the normalization of $(\partial\theta/\partial t)_s = 1$ and $\partial\bar{\theta}/\partial z = 1$.</figcaption></figure>

Any parameterization for the eddy heat flux must lead to the solution (2.3c) under the situation considered. This behavior is well supported both by the observations (Fig. 2 in Betts 1974) and laboratory experiments (Deardorff et al. 1969, 1974). Thus, in the following, the key question to be posed is how well the eddy-diffusion models can reproduce this theoretically predicted behavior.

### 3. Eddy Diffusion Model: the Simplest Case

We attempt to describe the evolution of the well-mixed layer described in the previous section by eddy-diffusion models. Its basic idea is to assume that the eddy heat flux is controlled by the mean potential temperature gradient. Thus,

$$ \langle w'\theta' \rangle = -\kappa \frac{\partial \langle \theta \rangle}{\partial z}, \qquad (3.1a) $$

where $\kappa$ is the eddy-diffusion coefficient, which may be a function of height, more likely through a dependence on physical variables. Substituting the above Eq. (3.1a) into Eq. (2.1), we obtain the eddy-diffusion equation of the form:

$$ \left( \frac{\partial}{\partial t} - \frac{\partial}{\partial z} \kappa \frac{\partial}{\partial z} \right) \langle \theta \rangle = 0 \qquad (3.1b) $$

In the following, we assume that all the equations are nondimensionalized by the characteristic scales of the system. Thus, for example, when the diffusion coefficient is constant, we simply set it as $\kappa = 1$. The characteristic scales of the system is discussed in the Appendix A.1.

Here, it is convenient to separate the total potential temperature, $\langle \theta \rangle$ into the initial mean profile, $\bar{\theta}$, defined by Eq. (2.2b) and the deviation from it, $\theta$:

$$ \langle \theta \rangle = \bar{\theta} + \theta \qquad (3.2a) $$
because the initial state, $\bar{\theta}$, is a trivial solution of Eq. (3.1b) when $\kappa$ is constant, thus we can simply focus on the latter “perturbation” part. Here, the term “perturbation” is adopted for the convenience, but used always with a quotation not to be ocnfused with the eddy components.

As a result, when $\kappa$ is considered constant, the differential equation to work on also reduces from Eq. (3.1b) into:

$$ \left( \frac{\partial}{\partial t} - \frac{\partial}{\partial z} \kappa \frac{\partial}{\partial z} \right) \theta = 0 \qquad (3.2b) $$

However, keep in mind that this simplification is valid only if the eddy diffusivity, $\kappa$, is assumed to be constant. Otherwise, the equation to be solved becomes:

$$ \frac{\partial}{\partial t} \theta - \frac{\partial}{\partial z} \kappa \frac{\partial}{\partial z} (\theta + \bar{\theta}) = 0 \qquad (3.2c) $$

The problem with the boundary condition (2.2a) can be solved analytically by representing the solution by the Fourier series:

$$ \theta = \sum_{n=1}^{\infty} \theta_n \sin k_n z \qquad (3.3a) $$

with the corresponding Fourier transform defined by

$$ \theta_n = \frac{2}{D} \int_0^D \theta \sin k_n z \, dz. \qquad (3.3b) $$

Here, $D$ is the vertical extent of the system, and the discrete wavenumbers are defined by

$$ k_n = \frac{n\pi}{D}. \qquad (3.3c) $$

After the nondimensionalization, the wavenumber is more conveniently given by

$$ k_n = \delta n \pi \qquad (3.4a) $$

where

$$ \delta = d/D \qquad (3.4b) $$

and $d$ is a characteristic vertical scale of the well-mixed layer.

Here, as it turns out, it becomes necessary to assume a system with a finite vertical extent, $D$, because when the Fourier series is replaced by the Fourier integral, the contributions at $k \to 0$ of the Fourier transform induces the singularities. In the following, we set $\delta = 0.1$. We also assume $\theta = 0$ at $z = D$. The sum in the Fourier series (3.3a) will be
truncated at $N = 1000$ in the following. A relatively large number of modes is required to represent the behavior of the solution to the vicinity of $z = 0$ well.

In transforming Eq. (3.2c) by the transform (3.3b), we particularly note that the transform of the second term becomes with partial integrals:

$$ \int_0^D \left[ \frac{\partial}{\partial z} \kappa \frac{\partial}{\partial z} (\theta + \bar{\theta}) \right] \sin k_n z \, dz = \kappa \frac{\partial(\theta + \bar{\theta})}{\partial z} \sin k_n z \Big|_0^D - k\kappa(\theta + \bar{\theta}) \cos k_n z \Big|_0^D - k^2 \int_0^D \kappa(\theta + \bar{\theta}) \sin k_n z \, dz $$

By using the boundary conditions, the above reduces into:

$$ \frac{2}{D} \int_0^D \left[ \frac{\partial}{\partial z} \kappa \frac{\partial}{\partial z} \right] (\theta + \bar{\theta}) \sin k_n z \, dz = \frac{2k}{D} \kappa_s \theta_s - k^2 [\kappa(\theta + \bar{\theta})]_n $$

Thus, Eq. (3.2b) is transformed into:

$$ \frac{d}{dt} \theta_n + k^2 [\kappa(\theta + \bar{\theta})]_n = \frac{2k}{D} \kappa_s \theta_s \quad (3.5a) $$

where the subscript, $s$, to $\kappa$ and $\theta$ suggests the values at the surface ($z = 0$), and $[\kappa(\theta + \bar{\theta})]_n$ is the Fourier transform of $\kappa(\theta + \bar{\theta})$.

The simplest case to consider is when the eddy diffusivity, $\kappa$, is a constant. In this case, $\bar{\theta}$ drops out from the equations, as already remarked. Without loss of generality, setting $\kappa = 1$ along with the other normalizations, the problem (3.5a) to be solved reduces into:

$$ \left( \frac{d}{dt} + k^2 \right) \theta_n = 2\delta kt \quad (3.5b) $$

also introducing a normalization to the boundary condition (2.2a).

This equation can be solved analytically and we obtain:

$$ \theta_n = \frac{2\delta}{k_n^3} (k_n^2 t - 1 + e^{-k_n^2 t}) \quad (3.5c) $$

Note that the above solution can further be simplified by introducing a wavenumber–depending time, $\tau_n = k_n^2 t$: the higher the wavenumbers, the faster the mode evolves with time.

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>
Fig. 1: The solution of the simple eddy diffusion model: evolution of the mean potential temperature profile is shown in the interval of $\delta t = 0.2$ for the “perturbation” (left) and total (right) fields. The initial condition ($t=0$) is in solid, and the varying curves moves upwards in $\theta$ with the increasing time.

Fig. 1 shows the solution (3.5c) for the “perturbation”, $\theta$, (left) and the total, $\langle \theta \rangle = \bar{\theta} + \theta$. A negative “perturbation” slope, $\partial\theta/\partial z < 0$, must be generated to transfer the heat upwards. Furthermore, for increasing the potential temperature, $\theta$, locally with time, this slope must gradually decreases with height, as seen in Fig. 1 (left). As a result, the eddy-diffusion formulation is fundamentally incapable to generating a well-mixed layer as seen in the total field of Fig. 1 (right).

# 4. Eddy Diffusion under the Second-Order Closure

## 4.1 General Solution

The main problem of the formulation with the constant eddy diffusivity is that the heat transport rate only increases proportional to the potential-temperature gradient regardless of the degree of the instability of the system. In order to achieve a well-mixed state, the rate of heat transport must further be enhanced with the increasing instabilities, thus the eddy diffusion itself becomes dependent on the degree of the instability of the system.

For considering those more general cases, thus, the eddy diffusivity is no longer a constant, but depends on both time and space. As a result, the procedure in the last section for obtaining the solution must also be generalized.

To construct such a general solution, we first set in Eq. (3.5a),

$$[\kappa(\theta + \bar{\theta})]_n = \tilde{\kappa}_n \theta_n \quad (4.1a)$$

keeping in mind that $\tilde{\kappa}_n$ evolves with time. As a result, in place of Eq. (3.5b), we obtain:

$$\left(\frac{d}{dt} + k^2 \tilde{\kappa}_n\right) \theta_n = 2\delta k t \quad (4.1b)$$

The above solution can readily be written in a closed form, and

$$\theta_n = 2\delta k \left[ e^{-k^2 K_n(t)} \int_0^t e^{k^2 K_n(t)} dt \right] \quad (4.2a)$$

where

$$K_n(t) = \int_0^t \tilde{\kappa}_n dt \quad (4.2b)$$

Here, we need to explicitly evaluate the Fourier transform, $(\kappa\theta)_n$, which is performed in a discrete form presented in the Appendix B. The time integrations in Eqs. (4.2a, b) is
based on the Adams-Bashforth third-order time stepping (cf., Durran 1991) with the step of $\Delta t = 10^{-5}$.

However, it turns out that the transformation between the Fourier and physical spaces that is required for evaluating $(\kappa\theta)_n$ is rather time consuming. Thus, a direct use of the formulation outlined above is rapidly discarded, and replaced by a first-order differentiation in the physical space. As for an alternative more efficient scheme, see Teixeria (1999).

## 4.2 Second-Order Closure

A standard procedure to alleviate the problem with the simple eddy diffusivity considered in the last section is to enhance the eddy diffusivity with the increasing instability. The so-called second order turbulence closure (cf., Sommeria 1975, Redelsberger and Sommeria 1981) attempts to achieve this goal by presenting the vertical heat flux in a more elaborated form. Outlines of the turbulence closure formulation are presented in the Appendix C, because this approach is not accessible to wider audience.

A standard form of the second-order closure assuming horizontal homogeneity is:

$$ \langle w' \theta' \rangle = -\frac{2}{3} l \frac{E^{1/2}}{C_s} \frac{\partial \langle \theta \rangle}{\partial z} + \frac{2}{3} l \frac{\alpha g}{C_s E^{1/2}} \langle \theta'^2 \rangle \quad (4.3) $$

(cf., Eq. 26 of Sommeria 1976). Here, $\alpha = 1/\theta_0$ is the coefficient of the thermal expansion with $\theta_0 \simeq 300$ K a reference potential temperature, $g$ the acceleration of gravity, $E$ the turbulent kinetic energy, $l$ is the mixing length, and $C_s$ is a standard constant introduced. Here, the term with $E^{1/2}$ can also be replaced by various other measures of the turbulent velocity, for example, its standard deviation (e.g., Cuiipers and Holtslag 1998).

The eddy flux formulation given by Eq. (4.3) is designed to alleviate the problem with a constant eddy diffusivity. First, by introducing the dependence of the eddy diffusivity on the turbulent kinetic energy, $E$, in the first term, the eddy diffusivity itself increases with the increasing convective instability, because it generates a more turbulent kinetic energy. Second, by adding an extra term (the second) to the original eddy-diffusion formulation, the eddy heat flux is directly enhanced by the increasing variability of the eddy potential temperature, $\theta'$.

Thus, how much those modifications can alleviate the problem with the constant eddy-diffusivity formulation? To answer this question, we need to derive a closed expression for the eddy heat flux from the above formulation (4.3) by introducing closed expressions for $E$ and $\langle \theta'^2 \rangle$. They are, respectively, presented under the standard second-order turbulence closure by

$$ E = \frac{l}{C_E} \left[ \frac{4}{15} \frac{l}{C_m} \left( \frac{\partial \langle u_i \rangle}{\partial x_j} + \frac{\partial \langle u_j \rangle}{\partial x_i} \right) \frac{\partial \langle u_i \rangle}{\partial x_j} + \frac{\alpha g}{E^{1/2}} \langle w' \theta' \rangle \right] \quad (4.4a) $$
$$ \langle \theta'^2 \rangle = - \frac{2l}{C_{\theta} E^{1/2}} \langle w' \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial z} \qquad (4.4b) $$

(cf., Eqs. 29 and 30 of Sommeria 1976). Here, $u_i$ is the i-th velocity component, and the sum over the indices, i and j, is kept implicit by following the literature. Further constants, $C_E$ and $C_{\theta}$, have also been introduced for the closures. In the present case, as being implicitly assumed, there is no background flow, thus Eq. (4.4a) is further simplified into:

$$ E = \frac{\alpha g l}{C_E E^{1/2}} \langle w' \theta' \rangle \qquad (4.4c) $$

By further substituting the eddy diffusion formula (3.1a) into the above, it reduces to:

$$ E = \kappa \frac{\alpha g l}{C_E E^{1/2}} \left( - \frac{\partial \langle \theta \rangle}{\partial z} \right) \qquad (4.4d) $$

First, substituting Eq. (4.4b) into Eq. (4.3), after some re-arrangements, it reduces to the form comparable to Eq. (3.1a) with the eddy diffusion defined by

$$ \kappa = \left[ 1 + \frac{4}{3} l^2 \frac{\alpha g}{C_s C_{\theta} E} \left( - \frac{\partial \langle \theta \rangle}{\partial z} \right) \right]^{-1} \frac{2}{3} l \frac{E^{1/2}}{C_s} \qquad (4.5a) $$

Note further from Eq. (4.4d),

$$ \frac{1}{E} \left( - \frac{\partial \langle \theta \rangle}{\partial z} \right) = \kappa \frac{\alpha g l}{C_E E^{1/2}} $$

thus the factor on the left-hand side of the above found in the denominator can be replaced by a term proportional to $\kappa$ on the right-hand side. Re-arrangement of this result leads to:

$$ \kappa = \frac{2}{3} l \left( \frac{2C_E}{C_{\theta}} + 1 \right) \frac{E^{1/2}}{C_s} \qquad (4.5b) $$

The remaining factor, $E^{1/2}$, can be eliminated by noting from Eq. (4.4d),

$$ E = \left[ \kappa \frac{\alpha g l}{C_E} \left( - \frac{\partial \langle \theta \rangle}{\partial z} \right) \right]^{2/3} . $$

Thus, we finally obtain

$$ \kappa = \left( \frac{2}{3C_s} \right)^{3/2} \frac{l^2}{C_E^{1/2}} \left( \frac{2C_E}{C_{\theta}} + 1 \right)^{3/2} \left( -\alpha g \frac{\partial \langle \theta \rangle}{\partial z} \right)^{1/2} . \qquad (4.5c) $$

It also follows that vertical heat flux, $\langle w' \theta' \rangle$, now becomes proportional to $(\partial \langle \theta \rangle / \partial z)^{3/2}$ rather than to $\partial \langle \theta \rangle / \partial z$, a slightly stronger sensitivity to the potential-temperature gradient, but only by a factor. Thus, one may foresee no dramatic effect by this modification even without evaluating the evolution under this formulation.
A more serious practical problem with the result (4.5c) is that it is valid only for negative gradients (i.e., unstable stratifications), $\partial\langle\theta\rangle/\partial z < 0$, and with positive gradients, $\partial\langle\theta\rangle/\partial z > 0$, it becomes singular, i.e., the eddy diffusive becomes an imaginary number. This singularity stems from the fact that the turbulent kinetic energy, $E$, is expected to be maintained only under unstable stratifications through Eq. (4.4d). This equation itself becomes invalid under stable stratifications.

In reality, the eddy flow still survives under stable stratifications. A simple remedy to taken into account of this issue would be to set a minimum, $E_{min}$, to the turbulent kinetic energy, thus Eq. (4.4d) is replaced by

$$E = \begin{cases} \kappa \frac{\alpha g l}{C_E E^{1/2}} \left(-\frac{\partial \langle \theta \rangle}{\partial z}\right), & \kappa \frac{\alpha g l}{C_E E^{1/2}} \left(-\frac{\partial \langle \theta \rangle}{\partial z}\right) > E_{min} \\ E_{min}, & \kappa \frac{\alpha g l}{C_E E^{1/2}} \left(-\frac{\partial \langle \theta \rangle}{\partial z}\right) \le E_{min} \end{cases} \quad (4.6a)$$

It follows that the eddy diffusivity formula (4.5c) is replaced by with the corresponding minimum eddy diffusivity,

$$\kappa = \kappa_{min} \quad (4.6b)$$

when

$$\left(\frac{2}{3C_s}\right)^3 \frac{l^4}{C_E} \left(\frac{2C_E}{C_\theta} + 1\right)^3 \left(-\alpha g \frac{\partial \langle \theta \rangle}{\partial z}\right) \le \kappa_{min}^2,$$

and otherwise, the definition (4.5c) is maintained.

The second-order closure formulation just introduced is examined under the normalizations, as would be accomplished by proper nondimensionalizations, in a similar manner as in the last section. After the normalization, the eddy-diffusion formulation presented by Eqs. (4.5c) and (4.6b) reduces to:

$$\kappa = \begin{cases} \left(-\frac{\partial \langle \theta \rangle}{\partial z}\right)^{1/2}, & \text{when } \frac{\partial \langle \theta \rangle}{\partial z} < -\kappa_{min}^2 \\ \kappa_{min}, & \text{otherwise.} \end{cases} \quad (4.7)$$

Here, $\kappa_{min}$ remains an arbitrary parameter to choose, but we set in the following to be $\kappa_{min} = 0$ so that the eddy diffusion is completely suppressed under stable stratifications.

Here, in normalizing Eq. (4.5c) into the form of Eq. (4.7), the mixing length, $l$, is assumed to be constant with height. Note that the standard formulation following Blackadar (1962) introduces a height-dependence of the mixing length over the surface layer. This dependency is neglected here: this simplification is consistent with the fact that the surface layer is effectively eliminated in the present study by directly prescribing the evolution of the surface potential temperature.
<figure><img src="image_4.png" /><figcaption>Fig. 2: Evolution of the potential temperature, $\langle\theta\rangle = \bar{\theta} + \theta$, with the eddy diffusion model under the second-order closure and with $\kappa_{min} = 0$, and the time steps, $\Delta t = 10^{-5}$ (left), $\Delta t = 10^{-6}$ (middle), and $\Delta t = 10^{-7}$ (right). The plotted interval is $\delta t = 0.2$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3: Evolution of the eddy diffusion, $\kappa$, under the second-order closure and with $\kappa_{min} = 0$, and the time steps, $\Delta t = 10^{-5}$ (left), $\Delta t = 10^{-6}$ (middle), and $\Delta t = 10^{-7}$ (right).</figcaption></figure>

Fig. 2 shows the evolution of the potential-temperature profile with time under this eddy-diffusion formulation with varying time steps, $\Delta t = 10^{-5}-10^{-7}$. In spite of the stability dependence of the eddy diffusivity, which enhances the diffusivity with the increasing instabilities, the gradient of the potential temperature in the unstable layer is more enhanced than with the case of the constant diffusivity. This ironical behavior is partially explained the fact that the height-dependence of the eddy diffusivity no longer requires the potential-temperature gradient to be curved with height. As a result, the system chooses (but why?) to achieve the required heat transfer rate by enhancing the potential-temperature gradient rather than its curvature as in the previous case. [On the contrary to the earlier results, no sensitivity on the time step, $\Delta t$, is seen.]
Fig. 3 shows the time evolution of the eddy-diffusivity profile: it is the largest at the surface, and rapidly decreases to zero towards the top of the convectively unstable layer. Here, we also note a weak sensitivity on the time step, $\Delta t$, which tends to decrease the surface eddy diffusivity slightly with the decreasing time steps. As a major improvement from the constant eddy-diffusivity case is that due to the decreasing eddy diffusivity with height, as a result, a well-defined top of the well-mixed layer is now identified (Fig. 2), a feature that was missing with the constant eddy diffusivity (Fig. 1).

The resulting large potential-temperature gradient may partially be attributed to a relatively weak dependence of the eddy diffusivity on the latter. Partially motivated from this shortcoming, Teixeria and Cheinet (2004) propose an additional dependence of the mixing length on the turbulent kinetic energy, which is expected to further enhance the dependence of the eddy diffusivity on the potential-temperature gradient through the relation (4.4d). However, as going to be shown in Sec. 5, this formulation is ill-posed under the free-convection regime considered herein.

## 4.3 Generalized Eddy-Diffusivity Dependence on $\partial\langle\theta\rangle/\partial z$

<figure><img src="image_6.png" /><figcaption>Fig. 4: The same as Fig. 2, but with $\nu = 1$ (left), $\nu = 2$ (middle), and $\nu = 3$ (right). The adopted time step are $\Delta t = 10^{-5}$ on the left, otherwise $\Delta t = 10^{-6}$.</figcaption></figure>

The failure of the eddy diffusivity under the second-order closure to produce a well-mixed boundary layer, even to the extent achieved with the constant eddy diffusivity, may be attributed to its rather weak dependence on the potential-temperature gradient. It may be speculated that a stronger dependence may alleviate this problem. To examine this possibility, the eddy-diffusion formulation (4.7) of the last subsection is generalized
into:

$$ \kappa = \begin{cases} \left(-\frac{\partial \langle \theta \rangle}{\partial z}\right)^{\nu}, & \text{when } \frac{\partial \langle \theta \rangle}{\partial z} < -\kappa_{min}^2 \\ \kappa_{min}, & \text{otherwise.} \end{cases} \qquad (4.8) $$

with an exponent, $\nu$, to be specified. The formulation of the last section is recovered with $\nu = 1/2$.

Fig. 4 show the same as Fig. 2, but with the choice of $\nu = 1$ (left), 2 (middle), and 3 (right). Against the anticipation, the strong potential-temperature gradient does not diminish with the increasing exponent, $\nu$. Nevertheless, we may note that the curvature slope diminishes with the increasing exponent, $\nu$. The same analysis is also repeated with the decreasing $\nu$, but no visible departure from the result of Fig. 2 even with $\nu = 0.1$ (not shown).

<figure><img src="image_7.png" /><figcaption>Fig. 5: The parameter choices as for Fig. 4, but the vertical profiles of $\kappa$. Note that the time step for $\nu = 1$ is also set $\Delta t = 10^{-6}$ here.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6: The same as for Fig. 5, but with $\nu = 0.1$.</figcaption></figure>
For aiding the interpretations in the next subsection, we further plot in Figs. 5 and 6, the vertical profiles of $\kappa$ with the varying exponent, $\nu$. Fig. 5 suggests that with the increasing exponent, $\nu$, the $\kappa$-profile asymptotically approaches a linear profile, whereas the comparisons of Figs. 3 and 6 suggest that as the exponential approaches to $\nu \to 0$, the diffusivity, $\kappa$, asymptotically approaches to constant with height, with a sudden decrease to zero at the top of the well-mixed layer.

## 4.4 Interpretation

To interpret the results so far in a concise manner, let us consider a solution of the form:

$$ \langle \theta \rangle = t - t_0 + \Delta\theta(z_0 - z)^q \quad (4.9a) $$

for the convective layer below $z = z_i$. Here, $t_0$, $\Delta\theta$, $z_0$, and $q$ are assumed to be constants. To satisfy the boundary condition of $\theta$ at $z = 0$ as given by Eq. (2.2a) with the normalization of $(\partial\theta/\partial t)_s = 1$, we need to set:

$$ t_0 = \Delta\theta z_0^q \quad (4.9b) $$

Obviously,

$$ \frac{\partial \langle \theta \rangle}{\partial z} = -q \Delta \theta (z_0 - z)^{(q-1)} $$

By substituting the above into Eq. (4.8), we find

$$ \kappa = (q \Delta \theta)^{\nu} (z_0 - z)^{\nu(q-1)} \quad (4.10a) $$

and furthermore,

$$ \frac{\partial}{\partial z} \kappa z = (\nu + 1)(q - 1)(q \Delta \theta)^{\nu+1} (z_0 - z)^{(\nu+1)(q-1)-1} \quad (4.10b) $$

By substituting the result (4.10b) into the diffusion equation (3.2b), also noting $\partial\langle\theta\rangle/\partial t = 1$, we obtain the relation:

$$ (\nu + 1)(q - 1)(q \Delta \theta)^{\nu+1} (z_0 - z)^{(\nu+1)(q-1)-1} = 1 \quad (4.11) $$

For the above condition (4.11) to be satisfied at every height, $z$, it requires to satisfy the following two subconditions:

$$ (\nu + 1)(q - 1) - 1 = 0, \quad (4.12a) $$

$$ (\nu + 1)(q - 1)(q \Delta \theta)^{\nu+1} = 1. \quad (4.12b) $$
From the first subcondition (4.12a), we can define the exponent, $q$, in the assumed solution (4.9a) as:

$$q = \frac{\nu + 2}{\nu + 1} \qquad (4.13a)$$

With the help of the first subcondition (4.12a), the second subcondition (4.12b) further reduces into:

$$(q\Delta\theta)^{\nu+1} = 1$$

that can be immediately re-interpreted as the condition for the constant, $\Delta\theta$ in the solution (4.9a), which is given by:

$$\Delta\theta = 1/q \qquad (4.13b)$$

The first result (4.13a) shows how the curvature of the potential-temperature profile changes with the changing exponent, $\nu$, in the definition of the diffusion (Eq. 4.8). As the exponent decreases towards $\nu \to 0$, the potential-temperature profile asymptotically approaches to a parabolic form as $q \to 2$, as suggested by Fig. 2. On the other hand, as the exponent, $\nu$, increases, asymptotically towards $\nu \to +\infty$, the potential-temperature profile asymptotically approaches to a linear with height with $q \to 1$ as suggested by Fig. 4.

On the other hand, $\Delta\theta$ measures a relative vertical extent of the well-mixed layer: it decreases with the increasing, $\nu$, counting on $q$-dependence on $\nu$, as just discussed. This tendency is also consistent with that of Fig. 4.

Furthermore, by substituting Eq. (4.13a,b) into the solution (4.10a) for the diffusion, we find the dependence of $\kappa$ with height is given by

$$\kappa = (z_0 - z)^{\nu/(\nu+1)} \qquad (4.13c)$$

The result predicts that the eddy diffusivity profile approaches to constant as $\nu \to 0$ and to a linear dependence as $\nu \to +\infty$. These predictions are also consistent with the tendencies unidentified in Figs. 3, 5, and 6.

## 5. Countergradient Term Formulations

Only a particular form of the second-order closure has been considered in Sec. 4.2. Yet, more or less all the variants take a generalization of the form (4.3), which can be written as:

$$\langle w'\theta' \rangle = -\kappa_0 \left( \frac{\partial \langle \theta \rangle}{\partial z} - \gamma \right) \qquad (5.1)$$

Here, the first term on the right-hand side is the eddy diffusion, and the second term is called the countergradient, because by setting $\gamma > 0$, this term can transfer heat in a counter gradient manner. Here, we have added the subscript 0 to the eddy diffusivity
above, because its definition differs from the original given by Eq. (3.1a): in the latter definition, the countergradient term is also considered a part of eddy diffusivity. The literature emphasizes that this term is required to represent the inversion layer at the top of the well-mixed convective boundary layer, an effect that has been neglected so far. However, it should also be emphasized that the countergradient term should also help to achieve a well-mixed state itself.

Here, in the case considered in Sec. 4.2, the eddy diffusion and the countergradient are, respectively, given by

$$ \kappa_0 = \frac{2}{3} l \frac{E^{1/2}}{C_s} \qquad (5.2a) $$

$$ \gamma = \frac{\alpha g}{C_s E} \langle \theta'^2 \rangle \qquad (5.2b) $$

Examination of this particular case already suggests that the above-stated aimed goal is not that easy to accomplish.

In this case, the countergradient term, $\gamma$, is positive definite, as should be. Thus, one may expect that this term may counteract well against the tendency of the eddy diffusion term to produce a finite potential-temperature slope. However, this anticipation is simply betrayed by substituting the additional standard closures in the countergradient term, (5.2b). As it turns out, this term merely contribute to introduce a weak dependence on the potential-temperature gradient in the modified “effective” eddy diffusivity, as seen in Eq. (4.5a).

The purpose of this section is to review some alternative countergradient formulations, and to examine how they have tried to overcome aforementioned difficulty.

## 5.1 Formulation by Cuiipers and Holtslag (1998)

Alternatively, Cuijpers and Holtslag (1998) propose the countergradient term given by

$$ \gamma = \frac{\beta_0 \tau w_*^2 \theta_*}{\kappa_0 z_i} \qquad (5.3) $$

Here, $\beta_0$ is a constant, $\tau$ is the relaxation time scale introduced in the pressure covariance closure, and $w_*$ and $\theta_*$ are the scales for the velocity and potential temperature, respectively. More precisely, $w_*\theta_*$ is a vertically averaged surface flux defined by

$$ w_* \theta_* = \frac{1}{z_i} \int_0^{z_i} \langle w' \theta' \rangle dz \qquad (5.4a) $$

whereas the velocity scale is defined by

$$ w_* = \left[ \frac{2}{1 - \mu} \int_0^{z_i} \alpha g \langle w' \theta' \rangle dz \right]^{1/3} \qquad (5.4b) $$
Here, $\mu$ is the faction of the surface flux that observed at the top of the well-mixed layer, i.e.,

$$ \langle w'\theta' \rangle \Big|_{z=z_i} = -\mu \langle w'\theta' \rangle_0 \qquad (5.4c) $$

and the subscript 0 on the right-hand side suggests a surface value.

Note that with a finite inversion and $\mu > 0$, the heat flux solution (2.3c) is modified into:

$$ \langle \theta' w' \rangle = \langle \theta' w' \rangle_0 \left[ 1 - (1 + \mu) \frac{z}{z_i} \right] \qquad (5.5) $$

Substituting Eq. (5.5) into Eqs. (5.4a, b), we find:

$$ w_* \theta_* = \frac{1 - \mu}{2} \langle w' \theta' \rangle_0 \qquad (5.6a) $$

$$ w_* = [\alpha g z_i \langle w' \theta' \rangle_0]^{1/3} \qquad (5.6b) $$

Thus, $w_*$ reduces to the Deadorff's (1970) velocity scale.

We may introduce an approximation, $w_* = E^{1/2}$, and substituting the definition of the eddy diffusivity (5.2a) and Eq. (5.6a) into the countergradient term (5.3), we find:

$$ \gamma = \frac{\gamma_0}{l} \langle w' \theta' \rangle_0 \qquad (5.7) $$

where $\gamma_0 = 3C_s\beta_0\tau/2$. Thus, the countergradient term is solely controlled by the surface flux. Note that qualitatively the same conclusion follows even when the flux profile does not follow Eq. (5.5).

It transpires that Cuiipers and Holtslag's countergradient-term formulation reduces to a kind of redundancy. Especially, when the eddy diffusion term can be neglected, substitution of Eq. (5.7) into Eq. (5.1) leads to:

$$ \langle w' \theta' \rangle = \kappa_0 \gamma = \kappa_0 \frac{\gamma_0}{l} \langle w' \theta' \rangle_0 \qquad (5.8) $$

Thus, by tuning the prefactor, $\kappa_0\gamma_0/l$, we can obtain the desired result (5.5) at will.

## 5.2 Formulation by Lumley (1975)

Yet, another alternative formulation is proposed by Lumley (1975):

$$ \gamma = - \frac{S_w \sigma_w \tau_L}{2 \kappa_0} \frac{\partial \langle w' \theta' \rangle}{\partial z} \qquad (5.8) $$

Here, $S_w = \langle w'^3 \rangle / \sigma_w^3$ is the vertical velocity skewness. Its substitution into Eq. (5.1) leads to:

$$ \langle w' \theta' \rangle = - \kappa_0 \frac{\partial \langle \theta \rangle}{\partial z} - \frac{S_w \sigma_w \tau_L}{2 \kappa_0} \frac{\partial \langle w' \theta' \rangle}{\partial z} \qquad (5.9) $$
Thus, the problem of defining the heat flux reduces to that of solving the differential equation (5.9). Although its precise solution depends on the closure for the skewness, $S_w$, by solving it assuming the coefficients to be prescribed, we find that the heat flux solution takes an integral form, including the nonlocal effects.

## 6. Modification to the Mixing length

Alternative means of enhancing vertical mixing with the increasing unstable stratification is to incorporate a dependency on a stability measure into the mixing length. Such a formulation is proposed by Teixeria and Cheinet (2004), and further tested by Teixeria et al. (2004). In order to introduce such a formulation, a particular attention is also paid so that the given formulation represents a certain nonlocal nature of heat transport. Under this consideration, they introduce a dependence of the mixing length on the turbulent kinetic energy, $E$.

From the dimensional consistency, the most reasonable formulation would be to set:

$$l = \tau E^{1/2} \qquad (6.1)$$

where $\tau$ is a characteristic time scale.

The effect on the solution by introducing the dependence of the form (6.1) into the mixing length can be considered for the two cases: when the eddy-diffusion formulation is considered in the standalone manner as in Sec. 3 and when the countergradient term is introduced. For the latter case, we specifically consider the second-order closure of Sec. 4.

### 6.1 Standalone Eddy-Diffusion Formulation

We first consider the eddy diffusivity in the standalone manner. In this case, the eddy diffusion is defined by Eq. (5.2a). Here, from the definition (6.1) of the mixing length, we obtain

$$lE^{1/2} = \tau E$$

Substituting this expression into Eq. (5.2a),

$$\kappa_0 = \frac{2\tau}{3C_s} E \qquad (6.2a)$$

Note that Eq. (6.2a) is equivalent to Eq. (2) of Teixeria and Cheinet (2004).

Alternatively, from Eq. (6.1),

$$lE^{-1/2} = \tau$$

Substituting the above on the right-hand side of Eq. (4.4d),

$$E = \kappa \frac{\alpha g \tau}{C_E} \left( -\frac{\partial \langle \theta \rangle}{\partial z} \right) \qquad (6.2b)$$
Substituting Eq. (6.2b), also setting $\kappa = \kappa_0$, into Eq. (6.2a), we find:

$$ \kappa_0 = \frac{2\tau}{3C_s}\kappa_0\frac{\alpha g\tau}{C_E}\left(-\frac{\partial\langle\theta\rangle}{\partial z}\right) $$

Consequently, the eddy diffusivity, $\kappa_0$, found on both sides of the equation cancels out. This means that the eddy diffusivity remains indeterminate under the closures introduced, i.e., Eqs. (6.1), (5.2a), and (4.4d). After this cancel out, we obtain the constraint:

$$ 1 = \frac{2 \alpha g \tau^2}{3 C_s C_E} \left( -\frac{\partial \langle \theta \rangle}{\partial z} \right) $$

which can be re-written as the constraint on the potential-temperature gradient:

$$ \frac{\partial \langle \theta \rangle}{\partial z} = -\frac{3 C_s C_E}{2 \alpha g \tau^2} \qquad (6.3) $$

Most importantly, the constrained potential-temperature gradient given above a negative slope, an unstable state. Furthermore, by substituting the standard values from Sommeria (1976), $C_s \simeq 4$ and $C_E \simeq 0.7$, as well as $\alpha = 1/300 \text{ 1/K}$, we obtain an estimate of

$$ \frac{\partial \langle \theta \rangle}{\partial z} = -10^{-4} \text{ K/m} \qquad (6.4) $$

which is a very weakly unstable state, close the neutral stratification expected from the well-mixed layer theory.

## 6.2 Second-Order Closure

When the second-order closure as in Sec. 4.1 is assumed, the eddy diffusion is now given by Eq. (4.5a), into that we substitute the relations, $l^2/E = \tau^2$ and $lE^{1/2} = \tau E$ so that the mixing-length dependence is eliminated, and:

$$ \kappa = \left[ 1 + \frac{4 \alpha g \tau^2}{3 C_s C_\theta} \left( -\frac{\partial \langle \theta \rangle}{\partial z} \right) \right]^{-1} \frac{2 \tau}{3 C_s} E \qquad (6.5) $$

Further substituting Eq. (6.2b), and eliminating $E$, we obtain:

$$ \kappa = \left[ 1 + \frac{4 \alpha g \tau^2}{3 C_s C_\theta} \left( -\frac{\partial \langle \theta \rangle}{\partial z} \right) \right]^{-1} \frac{2 \tau}{3 C_s} \kappa \frac{\alpha g \tau}{C_E} \left( -\frac{\partial \langle \theta \rangle}{\partial z} \right) $$

Again the eddy diffusion cancels out from both sides, thus $\kappa$ remains indeterminate. Re-arranging the above, we obtain an alternative constraint:

$$ \frac{\partial \langle \theta \rangle}{\partial z} = \frac{3 C_s C_E}{2 \alpha g \tau^2} \qquad (6.6) $$
The constraint is identical to the original constraint (5.3), but flipping the sign: it leads to a very weakly stably stratified state. Such a state may appear to be unintuitive, yet it is constant with the observation from a laboratory experiment (Deardorff et al. 1969, 1974) as well as observations (Fig. 2 in Betts 1974): the well-mixed convective boundary layer leads to a weakly stable stratified state. Thus, the second-order closure leads to a more realistic result that with the standalone eddy diffusion formulation.

## 6.3 With the Prognostic Turbulent Kinetic Energy

Considerations of the last two subsections with the two different closures show that the modified mixing-length formulation given by Eq. (6.1) leads to an indeterminate eddy diffusion. However, in defense of the original authors, they have adopted a different closure the turbulent kinetic energy. We now consider their closure: substituting the closures for the various associated variables, it reduces to:

$$ \frac{\partial E}{\partial t} = \frac{\partial}{\partial z} K_e \frac{\partial E}{\partial z} - \kappa \alpha g \frac{\partial \langle \theta \rangle}{\partial z} - \frac{C_E}{l} E^{3/2}, \quad (6.7a) $$

where $K_e$ is the eddy diffusivity for the turbulent kinetic energy.

To get a sense of the general behavior of the system (6.7a), we drop the eddy-diffusion term on the right-hand side, and consider a simplified system:

$$ \left( \frac{\partial}{\partial t} + \frac{C_E}{\tau} \right) E = \kappa \alpha g \frac{\partial \langle \theta \rangle}{\partial z} \quad (6.7b) $$

also substituting the definition (6.1) of the eddy diffusion on the last term in Eq. (6.7a).

This equation (6.7b) can easily be solved analytically, when we assume that right-hand side is constant with time, and this solution can also be considered an asymptotic tendency of the system when the evolution of the right-hand side is slow enough. Under these assumptions, we see that the system exponentially approaches to an equilibrium state with the adjustment time scale of $\tau/C_E$, and the final equilibrium state becomes identical to Eq. (6.2b). Further substituting the definition of the eddy diffusivity given by Eq. (6.2a), we recover the same constraint (6.3) as before. The same conclusion also follows when the form (6.5) of the eddy diffusivity given under the second-order closure, and in this case, the final equilibrium state leads to the constraint (6.6).

These considerations suggests that the actual value of the eddy diffusion is inconsequential so long as the transient behavior is short enough. Under this condition, the effective consequence of the modified mixing-length formulation (6.1) is nothing other than imposing a particular potential-temperature gradient given either by Eq. (6.3) or Eq. (6.6).
The prediction with Eq. (6.3) is consistent with the results of Teixeria and Cheinet (2004: see their Fig. 2) and Teixeria et al. (2004: see their Fig. 6) achieving a well mixed state of the convective boundary layer. At the same time, the analysis here suggests that the actual vertical profile of the eddy diffusivity is inconsequential, because the eddy diffusivity becomes indeterminate at a final equilibrium, although the original studies tend to emphasize its significance (cf., Fig. 4 of Teixeria and Cheinet 2004, Fig. 7 of Teixeria et al. 2004).

## 6.4 Alternative Mixing–Length Formulation

A more elaborated formulation for the mixing length is proposed by Bougeault and Lacarrére (1989). We briefly review their formulation in this subsection. As it turns out, their formulation leads to a similar difficulty as the case with that by Teixeria and Cheinet (2004).

Bougeault and Lacarrére (1989) take a parcel–based argument to define the mixing length. More precisely, the two parcel–based characteristic lengths are considered. The first is the distance, $l_u$, that the parcel must be lifted from a reference to height, $z$, upwards to attain the turbulent kinetic energy, $E$, at the given height by a work by the parcel buoyancy, $b$. Thus,

$$ \int_{z}^{z+l_u} b \, dz' = E, \qquad (6.8a) $$

where

$$ b = \alpha g [\langle \theta \rangle(z) - \langle \theta \rangle(z')]. \qquad (6.8b) $$

The second is the distance, $l_d$, that the parcel must be lowered to lose the given turbulent kinetic energy due to a work against the negative parcel buoyancy:

$$ - \int_{z-l_d}^{z} b \, dz' = E. \qquad (6.8c) $$

Under this formulation, the mixing length is defined as $l = (l_u l_d)^{1/2}$.

To obtain an insight on this formulation, let's assume a simpler situation that the potential temperature profile can locally be approximated to linear with height. Thus,

$$ \langle \theta \rangle = \theta_0 + \frac{\partial \langle \theta \rangle}{\partial z} (z' - z). \qquad (6.9) $$

In this case, the two characteristics lengths become identical, and $l_u = l_d = l$. The integrals on the left–hand side of Eqs. (6.8a, b) can be performed analytically, and we obtain:

$$ -\alpha g \frac{\partial \langle \theta \rangle}{\partial z} \frac{l^2}{2} = E. \qquad (6.10a) $$
Thus,

$$l = \left(\frac{2E}{\alpha g}\right)^{1/2} \left(-\frac{\partial\langle\theta\rangle}{\partial z}\right)^{1/2}. \qquad (6.10b)$$

Furthermore, by comparing the two expressions for E given by Eqs. (6.10a) and (4.4d), the mixing length is further reduced into

$$l = \frac{2\kappa}{C_E E^{1/2}} \qquad (6.10c)$$

However, the mixing length defined by Eq. (6.10b) or (6.10c) leads to an indefinite eddy diffusivity with the two closures considered in Secs. 6.1 and 6.2.

First, when the standalone eddy diffusivity is considered as in Sec. 6.1, by substituting the result (6.10c) into the definition (5.2a) of the eddy diffusivity, also setting $\kappa = \kappa_0$ in Eq. (6.10c), we obtain:

$$\kappa_0 = \frac{4\kappa_0}{3} \frac{1}{C_s C_E}$$

Again, the eddy diffusion cancels out on both sides, and we obtain the constraint:

$$C_s C_E = \frac{4}{3}. \qquad (6.11a)$$

Unfortunately, the constraint is not consistent with the standard values of $C_s \simeq 4$ and $C_E \simeq 0.7$.

When the second-order closure is adopted as in Sec. 6.2, we substitute the relation (6.10c) in the form (4.5a) of the eddy diffusivity in this case. Again, the factor, $\kappa$, cancels out on both sides, we obtain the constraint

$$1 + \frac{4}{C_s C_e} = 0. \qquad (6.11b)$$

In this case, the given constraint is impossible satisfy as both $C_s$ and $C_e$ are positive definite.

A more serious problem with the mixing-length formulation by Bougeault and Lacarrére (1989) is that the equilibrium potential-temperature profile, $\langle\theta\rangle$, itself remains indeterminant. This happens, because the given consistency condition no longer constraints this profile.

## 7. Third-Order Closure

The well-mixed convective layer is expected to be maintained at the state of neutrality, $\partial\langle\theta\rangle/\partial z = 0$, by maintaining a nonvanshing vertical heat flux, i.e., $\langle w'\theta' \rangle \neq 0$. First, the simple eddy-diffusion formulation given by Eq. (2.1) never satisfies this condition.
A similar difficulty is also faced with the second-order closure due to the only limited turbulent variables, namely, the turbulent kinetic energy, $E$, and the potential-temperature variance, $\langle \theta'^2 \rangle$. In the steady limit, as considered so far, those two variables only depends on the mean potential temperature gradient, $\partial \langle \theta \rangle / \partial z$ (cf., Eqs. 4.4b, d), as a result the heat flux, $\langle w' \theta' \rangle$, is also controlled solely by the mean potential temperature gradient, $\partial \langle \theta \rangle / \partial z$. Thus, the formula for the heat flux reduces to a type of eddy-diffusion formulation in one way or another (e.g., Eq. 4.5c) in one way or another. We have also found that not only that those formulations fail to satisfy the expected asymptotic tendency of $\langle w' \theta' \rangle \neq 0$ as $\partial \langle \theta \rangle / \partial z \to 0$, but also face a quantitative difficulty of maintaining the mean potential temperature gradient, $\partial \langle \theta \rangle / \partial z$, weak enough.

We have furthermore, reviewed some of the existing efforts overcoming this difficulty by introducing an additional countergradient term, $\gamma$, or a state-dependence in the mixing length, $l$. However, those efforts turn out to be rather unsatisfactory: especially the proposed elaborations to the mixing length considered here lead to the eddy diffusivity indefinite in the steady limit, although it may approach to an almost neutral state by an adjustment process.

For those reasons, we now turn our attention to the third-order closure. Here, we follow the third-order closure formulation proposed by André et al. (1978), but keeping the consistency with the second-order closure considered in Sec. 4. More specifically, by following Sommeria (1976), we set $\langle w'^2 \rangle = 2E/3$, and translate their equation for $\langle w'^2 \rangle$ for that of $E$. Furthermore, we focus on the steady state as before, although the original formulation by André et al. (1978) is time dependent.

Under the third-order closure, additional terms arising from the third-order turbulent variables are added to Eqs. (4.3), (4.4c), and (4.4b):

$$ \langle w' \theta' \rangle = -\frac{2l}{3} \frac{E^{1/2}}{C_s} \frac{\partial \langle \theta \rangle}{\partial z} + \frac{2l}{3} \frac{\alpha g}{C_s E^{1/2}} \langle \theta'^2 \rangle - \frac{l}{C_s E^{1/2}} \frac{\partial}{\partial z} \langle w'^2 \theta' \rangle \quad (7.1a) $$

$$ E = \frac{\alpha g l}{C_E E^{1/2}} \langle w' \theta' \rangle - \frac{l}{C'_E E^{1/2}} \frac{\partial}{\partial z} \langle w'^3 \rangle \quad (7.1b) $$

$$ \langle \theta'^2 \rangle = -\frac{l}{C_{\theta} E^{1/2}} \left[ 2 \langle w' \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial z} + \frac{\partial}{\partial z} \langle w' \theta'^2 \rangle \right] \quad (7.1c) $$

Here, $C'_E$ is an additional constant to be introduced. The derivation of the above starting from Eqs. (C.1.1), (C.2.1a, b), and (C.3.1a, b) is analogous to that of the second-order closure considered in the Appendix C, but retaining the third-order terms.

The higher-order turbulent variables introduced in Eqs. (7.1a, b, c) are defined by

$$ \langle w'^3 \rangle = \frac{3\alpha g l}{C_{w\theta\theta} E^{1/2}} \langle w'^2 \theta' \rangle - \frac{4l}{3C_{w\theta\theta}} E^{1/2} \frac{\partial E}{\partial z} \quad (7.2a) $$
$$ \langle w'^2 \theta' \rangle = \frac{2\alpha g l}{C_{ww\theta} E^{1/2}} \langle w' \theta'^2 \rangle - \frac{l}{C_{ww\theta} E^{1/2}} \left[ \langle w'^3 \rangle \frac{\partial \langle \theta \rangle}{\partial z} + \frac{4}{3} \langle w' \theta' \rangle \frac{\partial E}{\partial z} + \frac{2}{3} E \frac{\partial \langle w' \theta' \rangle}{\partial z} \right] \quad (7.2b) $$

$$ \langle w' \theta'^2 \rangle = \frac{\alpha g l}{C_{w\theta\theta} E^{1/2}} \langle \theta'^3 \rangle - \frac{2l}{C_{w\theta\theta} E^{1/2}} \left[ \langle w'^2 \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial z} + \frac{E}{3} \frac{\partial \langle \theta'^2 \rangle}{\partial z} + \langle w' \theta' \rangle \frac{\partial \langle w' \theta' \rangle}{\partial z} \right] \quad (7.2c) $$

$$ \langle \theta'^3 \rangle = -\frac{3l}{C_{\theta\theta\theta} E^{1/2}} \left[ \langle w' \theta'^2 \rangle \frac{\partial \langle \theta \rangle}{\partial z} + \langle w' \theta' \rangle \frac{\partial \langle \theta'^2 \rangle}{\partial z} \right] \quad (7.2d) $$

by introducing the closures in analogous manner as for the second-order variables. Here, $C_{w\theta\theta}$, $C_{ww\theta}$, and $C_{\theta\theta\theta}$ are the constants arising from the third-order closures.

## 7.1 Heat–Flux Solution

Unfortunately, this seven-variable system above is rather involved, and it is no longer easy to derive an explicit closed solution for the heat flux, as in the case with the second-order closure. In the following, we infer the behavior of this system from the two perspectives. First, in this subsection, we derive a closed form of the solution for the heat flux from Eqs. (7.1a, b, c), assuming that the higher-order variables can be considered given. In the next subsection, in turn, we derive the equations for solving those higher-order variables separately, assuming that the second-order variables are given. In principle, iterations for those two separate procedures lead to the full solution of this system.

In deriving a closed solution for the heat flux in terms of given higher-order variables, we first note that the turbulent kinetic energy is, from Eq. (7.1b), given in a closed form in terms of the other second-order variables as:

$$ E = \left[ \frac{\alpha g l}{C_E} \langle w' \theta' \rangle - \frac{l}{C'_E} \frac{\partial}{\partial z} \langle w'^3 \rangle \right]^{2/3} \quad (7.3) $$

In remainder of this subsection, the turbulent kinetic energy, $E$, is considered given, and skip the final step of substituting Eq. (7.3) into the solution for the heat flux to obtain a final closed form.

By substituting Eq. (7.1c) into Eq. (7.1a), we obtain a standalone differential equation to be solved for $\langle w' \theta' \rangle$ as:

$$ \langle w' \theta' \rangle^2 + \lambda_1 \langle w' \theta' \rangle + \lambda_2 = 0 \quad (7.4) $$

where

$$ \lambda_1 = \frac{2 l C_E E^{1/2}}{3 C_s} \left( \frac{2}{C_\theta} + \frac{1}{C_E} \right) \frac{\partial \langle \theta \rangle}{\partial z} - \frac{C_E}{\alpha g C'_E} \frac{\partial \langle w'^3 \rangle}{\partial z} + \frac{l}{C_s E^{1/2}} \mathcal{F} \quad (7.5a) $$

$$ \lambda_2 = -\frac{l C_E}{\alpha g C_s C'_E} \left[ \frac{2}{3} E^{1/2} \frac{\partial \langle \theta \rangle}{\partial z} + \frac{1}{E^{1/2}} \mathcal{F} \right] \frac{\partial \langle w'^3 \rangle}{\partial z} \quad (7.5b) $$
We have furthermore introduced a term, $\mathcal{F}$, representing a type of flux divergence:

$$ \mathcal{F} = \frac{2\alpha g l}{3C_{\theta}E^{1/2}} \frac{\partial \langle w' \theta'^2 \rangle}{\partial z} + \frac{\partial \langle w'^2 \theta' \rangle}{\partial z} \qquad (7.5c) $$

Eq. (7.4) is, being quadratic, readily solved in a closed form. Yet, the obtained solution is rather involved to be interpreted easily. Thus, instead, we consider some imitating cases as follows.

i) Assume that all the higher-order variables are small, and re-set as $\langle w'^3 \rangle = \epsilon \langle w'^3 \rangle$ and $\mathcal{F} = \epsilon \mathcal{F}$ with a small parameter, $\epsilon$. Expanding the solution as we as the parameters as:

$$ \begin{aligned} \langle w' \theta' \rangle &= \langle w' \theta' \rangle^{(0)} + \epsilon \langle w' \theta' \rangle^{(1)} + \dots \\ \lambda_1 &= \lambda_1^{(0)} + \epsilon \lambda_1^{(1)} + \dots \\ \lambda_2 &= \lambda_2^{(0)} + \epsilon \lambda_2^{(1)} + \dots, \end{aligned} $$

we find that $\lambda_2^{(0)} = 0$. Thus, to the leading order, Eq. (7.4) reduces to

$$ \langle w' \theta' \rangle^{(0)2} + \lambda_1^{(0)} \langle w' \theta' \rangle^{(0)} = 0 $$

It is immediately solved into:

$$ \langle w' \theta' \rangle^{(0)} = -\lambda_1^{(0)} \qquad (7.6a) $$

Substituting the leading-order parameter,

$$ \lambda_1^{(0)} = \frac{2 l C_E E^{1/2}}{3 C_s} \left( \frac{2}{C_\theta} + \frac{1}{C_E} \right) \frac{\partial \langle \theta \rangle}{\partial z}, $$

into the above, we recover the standard eddy-diffusion formula (3.1) with the eddy diffusivity defined by Eq. (4.5b).

To $O(\epsilon)$, the equation to be solved becomes

$$ [2\langle w' \theta' \rangle^{(0)} + \lambda_1^{(0)}] \langle w' \theta' \rangle^{(1)} + \lambda_1^{(1)} \langle w' \theta' \rangle^{(0)} + \lambda_2^{(1)} = 0. $$

After eliminating the leading-order solution, $\langle w' \theta' \rangle^{(0)}$, by Eq. (7.6a), we obtain

$$ \langle w' \theta' \rangle^{(1)} = -\lambda_1^{(1)} + \lambda_2^{(1)} / \lambda_2^{(0)}. \qquad (7.6b) $$

Here, the first-order parameters are

$$ \begin{aligned} \lambda_1 &= -\frac{C_E}{\alpha g C'_E} \frac{\partial \langle w'^3 \rangle}{\partial z} + \frac{l}{C_s E^{1/2}} \mathcal{F}, \\ \lambda_2 &= -\frac{2l C_E E^{1/2}}{3\alpha g C_s C'_E} \frac{\partial \langle \theta \rangle}{\partial z} \frac{\partial \langle w'^3 \rangle}{\partial z}. \end{aligned} $$
Substitution of them into Eq. (7.6b) leads to:

$$ \langle w' \theta' \rangle^{(1)} = \frac{C_E}{\alpha g C'_E} \left(1 + \frac{C_\theta}{2C_E}\right)^{-1} \frac{\partial \langle w'^3 \rangle}{\partial z} - \frac{l}{C_E E^{1/2}} \mathcal{F}. \quad (7.7) $$

ii) When the leading-order evolution is dictated by the nonlinear flux divergence, $\mathcal{F}$, with an almost neutral state. Thus, we re-set as: $\partial\langle\theta\rangle/\partial z = \epsilon\partial\langle\theta\rangle/\partial z$ and $\langle w'^3 \rangle = \epsilon\langle w'^3 \rangle$.

As the last case, $\lambda_2^{(0)} = 0$, thus the leading-order solution is given by Eq. (7.6a), and substitution of an explicit form for $\lambda_1^{(0)}$ finds:

$$ \langle w' \theta' \rangle^{(0)} = - \frac{l}{C_s E^{1/2}} \mathcal{F}. \quad (7.8a) $$

Note that this is the identical to the second term of the solution (7.7).

To $O(\epsilon)$, the parameters are:

$$ \begin{aligned} \lambda_1^{(1)} &= \frac{2}{3} l \frac{E^{1/2}}{C_s} \left( \frac{2C_E}{C_\theta} + 1 \right) \frac{\partial \langle \theta \rangle}{\partial z} - \frac{C_E}{\alpha g C'_E} \frac{\partial \langle w'^3 \rangle}{\partial z}, \\ \lambda_2^{(1)} &= - \frac{l C_E}{\alpha g C_s C'_E E^{1/2}} \mathcal{F} \frac{\partial \langle w'^3 \rangle}{\partial z}. \end{aligned} $$

Substitution of them into Eq. (7.6b) leads to:

$$ \langle w' \theta' \rangle^{(1)} = - \frac{2 l E^{1/2}}{3 C_s} \left( \frac{2 C_E}{C_s} + 1 \right) \frac{\partial \langle \theta \rangle}{\partial z} \quad (7.8b) $$

This has the form identical to the leading-order solution for the case i), and also consistent with Eqs. (3.1a) and (4.5b).

iii) Assume an almost neutral state as for the case ii), but now assume that the vertical-velocity skewness is dominant to the leading order. Thus, we re-set as: $\partial\langle\theta\rangle/\partial z = \epsilon\partial\langle\theta\rangle/\partial z$ and $\mathcal{F} = \epsilon\mathcal{F}$.

Reduction proceeds in the analogous manner as in the last two cases, we obtain:

$$ \langle w' \theta' \rangle^{(0)} = \frac{C_E}{\alpha g C'_E} \frac{\partial \langle w'^3 \rangle}{\partial z}, \quad (7.9a) $$

$$ \langle w' \theta' \rangle^{(1)} = - \frac{4}{3} l \frac{C_E}{C_s C_\theta} E^{1/2} \frac{\partial \langle \theta \rangle}{\partial z}. \quad (7.9a) $$

The leading-order solution, $\langle w' \theta' \rangle^{(0)}$, can be compared with the first term in the solution (7.7), whereas the first-order solution, $\langle w' \theta' \rangle^{(1)}$, can be compared with the solution (7.8b).

In the above analysis, the last two cases are more significant, because it shows that under the third-order closure, a finite nonvanishing heat flux is found to the leading order
even in the asymptotic limit of neutral state, i.e., $\partial\langle\theta\rangle/\partial z \to 0$ (cf., Eqs. 7.8a and 7.9a). As a whole, it is found that the positive vertical-velocity skewness gradient, $\partial\langle w'^3\rangle/\partial z$, and the negative nonlinear flux divergence, $\mathcal{F}$, contribute to the positive heat flux, the sign expected for the well-mixed convective layer.

Note further that the vertical-velocity skewness gradient, $\partial\langle w'^3\rangle/\partial z$, is defined by Eq. (7.2a) through the flux, $\langle w'^2\theta'\rangle$, and the turbulence kinetic energy, $E$:

$$ \frac{\partial\langle w'^3\rangle}{\partial z} = \frac{3\alpha g l}{C_{w\theta\theta}} \frac{\partial}{\partial z} \left( \frac{1}{E^{1/2}} \langle w'^2\theta' \rangle \right) - \frac{4l}{3C_{w\theta\theta}} \frac{\partial}{\partial z} \left( E^{1/2} \frac{\partial E}{\partial z} \right) $$

Thus, the negative nonlinear-flux divergence, $\partial\langle w'^2\theta'\rangle/\partial z$, tends to also lead to negative $\partial\langle w'^3\rangle/\partial z$ unless $E$ strongly decreases with height, thus the sign can be inverted.

## 7.2. Third-Order Moments

We now turn attention to the third-order moments defined by Eqs. (7.2a, b, c, d). We adopt the identical strategy as the last subsection, but now the second-order moments are considered given, and analyze the third moments under this assumption. More specifically, we focus on deriving a closed equation for $\langle w'\theta'^2\rangle$. Paying the attention to the fact that $\langle\theta'^2\rangle$ explicitly depends on $\langle w'\theta'^2\rangle$, we eliminate $\langle\theta'^2\rangle$, and derive a closed equation with the two other second-order moments, $\langle w'\theta'\rangle$ and $E$.

We first note that $\langle\theta'^2\rangle$ appears in Eqs. (7.2a, b, c) with the vertical derivative. It is given by:

$$ \frac{\partial\langle\theta'^2\rangle}{\partial z} = -\frac{l}{C_{\theta}} \left[ 2\frac{\partial}{\partial z}\left(\frac{\langle w'\theta'\rangle}{E^{1/2}}\frac{\partial\langle\theta\rangle}{\partial z}\right) + \frac{\partial}{\partial z}\left(\frac{1}{E^{1/2}}\frac{\partial}{\partial z}\langle w'\theta'^2\rangle\right) \right]. \quad (7.10) $$

The basic procedure to follow is to first eliminate $\langle w'^3\rangle$ and $\langle\theta'^3\rangle$, respectively, from $\langle w'^2\theta'\rangle$ and $\langle w'\theta'^2\rangle$, deriving equations for the latter given in terms of $\langle w'\theta'\rangle$ and $E$ only. The final step is to derive a closed equation for $\langle w'\theta'^2\rangle$.

First, substituting Eqs. (7.2a, d), respectively, into Eqs. (7.2b, c), we obtain:

$$ \langle w'^2 \theta' \rangle = \left[ 1 + \frac{3\alpha g l^2}{C_{ww\theta} C_{w\theta\theta} E} \frac{\partial \langle \theta \rangle}{\partial z} \right]^{-1} \left\{ \frac{2\alpha g l}{C_{ww\theta} E^{1/2}} \langle w' \theta'^2 \rangle + \frac{4l^2}{3C_{ww\theta} C_{w\theta\theta}} \frac{\partial E}{\partial z} \frac{\partial \langle \theta \rangle}{\partial z} - \frac{l}{C_{ww\theta} E^{1/2}} \left[ \frac{4}{3} \langle w' \theta' \rangle \frac{\partial E}{\partial z} + \frac{2}{3} E \frac{\partial \langle w' \theta' \rangle}{\partial z} \right] \right\} \quad (7.11a) $$

$$ \langle w' \theta'^2 \rangle = - \left[ 1 + \frac{3\alpha g l^2}{C_{w\theta\theta} C_{\theta\theta\theta} E} \frac{\partial \langle \theta \rangle}{\partial z} \right]^{-1} \left\{ \frac{l}{C_{w\theta\theta}} \left[ \frac{3\alpha g l}{C_{\theta\theta\theta} E} \langle w' \theta' \rangle + \frac{2}{3} E^{1/2} \right] \frac{\partial \langle \theta'^2 \rangle}{\partial z} \right. \\ \left. + \frac{2l}{C_{w\theta\theta} E^{1/2}} \left[ \langle w'^2 \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial z} + \langle w' \theta' \rangle \frac{\partial \langle w' \theta' \rangle}{\partial z} \right] \right\} \quad (7.11b) $$
We next substitute Eq. (7.10) into Eq. (7.11b) to eliminate $\partial\langle\theta'^3\rangle/\partial z$ from the latter to obtain:

$$ \left\{ \frac{l^2}{C_\theta C_{w\theta\theta}} \left[ \frac{3\alpha gl}{C_{\theta\theta\theta}E} \langle w'\theta' \rangle + \frac{2}{3}E^{1/2} \right] \frac{\partial}{\partial z} \frac{1}{E^{1/2}} \frac{\partial}{\partial z} - \left[ 1 + \frac{3\alpha gl^2}{C_{w\theta\theta}C_{\theta\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial z} \right] \right\} \langle w'\theta'^2 \rangle \\ = - \frac{l^2}{C_\theta C_{w\theta\theta}} \left[ \frac{3\alpha gl}{C_{\theta\theta\theta}E} \langle w'\theta' \rangle + \frac{2}{3}E^{1/2} \right] \frac{\partial}{\partial z} \left( \frac{\langle w'\theta' \rangle}{E^{1/2}} \frac{\partial\langle\theta\rangle}{\partial z} \right) + \frac{2l}{C_{w\theta\theta}E^{1/2}} \left[ \langle w'^2\theta' \rangle \frac{\partial\langle\theta\rangle}{\partial z} + \langle w'\theta' \rangle \frac{\partial\langle w'\theta' \rangle}{\partial z} \right] \quad (7.11c) $$

Further substituting Eq. (7.11a) into the above, $\langle w'^2\theta' \rangle$ can be eliminated, and we obtain a closed equation for $\langle w'\theta'^2 \rangle$ as

$$ \left\{ \frac{l^2}{C_\theta C_{w\theta\theta}} \left[ 1 + \frac{3\alpha gl^2}{C_{ww\theta}C_{w\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial z} \right] \left[ \frac{3\alpha gl}{C_{\theta\theta\theta}E} \langle w'\theta' \rangle + \frac{2}{3}E^{1/2} \right] \frac{\partial}{\partial z} \frac{1}{E^{1/2}} \frac{\partial}{\partial z} - \mathcal{D} \right\} \langle w'\theta'^2 \rangle \\ = \frac{4l^2}{3C_{ww\theta}C_{w\theta\theta}E^{1/2}} \frac{\partial\langle\theta\rangle}{\partial z} \left\{ 2 \left[ \frac{l}{C_{w\theta\theta}} \frac{\partial\langle\theta\rangle}{\partial z} \right] \frac{\partial E}{\partial z} - \frac{E^{1/2}}{3} \frac{\partial\langle w'\theta' \rangle}{\partial z} \right\} \\ + \left[ 1 + \frac{3\alpha gl^2}{C_{ww\theta}C_{w\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial z} \right] \left\{ \frac{2l}{C_{ww\theta}E^{1/2}} \langle w'\theta' \rangle \frac{\partial\langle w'\theta' \rangle}{\partial z} \right. \\ \left. - \frac{l^2}{C_\theta C_{w\theta\theta}} \left[ \frac{3\alpha gl}{C_{\theta\theta\theta}E} \langle w'\theta' \rangle + \frac{2}{3}E^{1/2} \right] \frac{\partial}{\partial z} \left( \frac{\langle w'\theta' \rangle}{E^{1/2}} \frac{\partial\langle\theta\rangle}{\partial z} \right) \right\} (7.12a) $$

where

$$ \mathcal{D} = \left[ 1 + \frac{3\alpha gl^2}{C_{ww\theta}C_{w\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial t} \right] \left[ 1 + \frac{3\alpha gl^2}{C_{w\theta\theta}C_{\theta\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial t} \right] + \frac{4\alpha gl^2}{C_{ww\theta}C_{w\theta\theta}E} \frac{\partial\langle\theta\rangle}{\partial t} \quad (7.12b) $$

When the right-hand side is considered given, the above equation takes a form of the Poisson equation, involving the second-order vertical derivative on the left-hand side, which takes a form of the Helmholtz differential operator as a whole. Since the inverse of the Helmholtz operator is nonlocal, the response of the third-order flux, $\langle w'\theta'^2 \rangle$, to the given second-order variables on the right-hand side, also becomes nonlocal. This nonlocal response of $\langle w'\theta'^2 \rangle$ further makes the behavior of the heat flux, $\langle w'\theta' \rangle$, also nonlocal through Eq. (7.4). In this manner, on the contrary to the common understanding, the moment-based turbulence description with the third-order closure describes the nonlocal transport process.

Eq. (7.12a) itself is too involved to obtain any further insights, but the solution to Eq. (7.12a) can further be investigated by taking various limiting cases, as considered in the last subsection.

## 8 Conclusions
• A simple problem of the well-mixed convective layer is considered with a formulation of the constant eddy diffusivity. This model, indeed presents some tendency for convective neutrality in lower layers, but it is fundamentally incapable of simulating the same in the upper layers. As a result, no well-defined top of the mixed layer is identified in simulation. Instead, the stratification only asymptotically approaches to the background state with height.

• The standard second-order closure, as proposed by Sommeria (1976), merely modifies the eddy diffusivity to be weakly depending on the stratification, which tends to enhance the heat transport with the more convective instability. Furthermore, this modification merely deteriorates the incapacity of the eddy diffusivity to well mix the convectively-unstable layer.

• Attempts have been made in the literature to overcome the limitations of the second-order closure approaches by introducing a countergradient term or an explicit state-dependence to the mixing length. However, our analysis shows that those existing proposals contain various unsatisfactory features. More specifically:

i) The countergradient formulation by Cuiipers and Holtslag (1998) embeds the surface-flux control of the potential-temperature, more or less, explicitly into the formulation. Thus, the problem reduces to that of properly tuning the height-dependence of the remaining undefined parameters to fit the height-dependence of the surface flux consistent with that of the well-mixed convective layer. This approach works only when the boundary layer is strongly controlled by the surface fluxes, but that is not necessarily the case under strong stratifications.

ii) Proposed elaborations of the mixing length proposed by Teixeria and Cheinet (2004) and Bougeault and Lacarrére (1989) tend to make the behavior of the system singular in the sense that in the equilibrium limit, the eddy diffusivity becomes indeterminate. In the case of Teixeria and Cheinet (2004), the consistency condition predicts an almost neutrally-stratified state, and we have also shown that their formulation can adjust the system towards an almost neutral state by explicitly introducing a time-dependent adjustment process. On the other hand, the formulation by Bougeault and Lacarrére (1989) even fails to constrain the mean potential-temperature profile by consistency at the equilibrium.

• The third-order closure leads to a nonlocal description of the vertical eddy transport, as demonstrated for the heat flux herein.

The last conclusion is particularly important, because it is often argued that the standard moment-based description of the turbulence is fundamentally local, thus it is incapable of representing long-range transport as realized by plumes in convective layers.
Yet, the work is still be performed to more convincingly demonstrate that a higher-order closure model can indeed represent the plume-type long-range transport.

# 9 Discussions

Eddy-flux formulations and their higher-order closure extensions for global atmospheric models have already been extensively tested in the literature, but mostly in the context of operational implementations. Thus, the vertical resolutions tend to be the crude (e.g., Holtslag and Boville 2993), as often the case in operational settings. Consequently, the reported results do not necessarily represent theoretically-predicted performance of those schemes. Rather, the results tend to be about a performance of a model under a given particular vertical resolution as well as numerical implementations.

On the other hand, by taking an idealized setting with more careful considerations of numerical algorithms, the present study intends to elucidate the theoretically-predicted performance of those eddy-flux schemes. Such an exercise is crucial, because without distinguishing it from numerical issues, it is not possible to well identify the limitations of those existing models. However, somehow such fundamental investigations are missing from the literature. The present work have attempted to establish such baseline behaviors of the eddy-diffusion models and their extensions. A major advantage of the present study is that by adopting a simple well-mixed convective layer system without any large-scale background flow, it becomes much easier to elucidate the behavior of the system analytically.

The performed analysis also points to some general issues in examining the existing proposed formulations systematically. First is an importance distinction between the steady and time-dependent formulations. Especially, some of the remedies to the second-order closures proposed under time-dependent formulations lead to singular behaviors under the steady formulation.

Second, although it is commonly argued that those extensions of the eddy-diffusion based formulations are fundamentally incapable of representing the nonlocal transport processes, the present analysis has demonstrated otherwise. This argument provides a justification for combining the eddy-diffusion (ED) and mass-flux (MF) formulations together, leading to the notion of EDMF. Yet, the present analysis suggests that it may not be crucial to introduce MF to ED, but the limitations of ED can be overcome by proceeding to higher-order closures.

It may also be remarked that rich behaviors of the turbulence-closure models are missing out from the present analysis. For example, it is known that when horizontal homogeneity is introduced parameterically with $\lambda^2 > 0$, eddy diffusivity can be augmented singularly
with Ri (Sommeria 1976). Unfortunately, we miss to address this aspect. Yet, we also emphasize that this great deal of simplification leads to greater elucidations of the model behaviors

**acknowledgments** I acknowledge helpful discussions with Vince Larson.

# Appendix A Basic Considerations of the Eddy-Diffusion Model

## A.1 Scales of the System

Eq. (3.2b) can be nondimensionalized by introducing the characteristic scales for the time and in vertical, $\Delta t$ and $d$, respectively. As a result, the nondimensionalized eddy-diffusion coefficient, $\tilde{\kappa}$, is given by

$$ \tilde{\kappa} = \frac{\kappa \Delta t}{d^2} \qquad (A.1.1a) $$

By setting $\tilde{\kappa} = 1$, the characteristic time scale, $\Delta t$, is defined by

$$ \Delta t = \frac{d^2}{\kappa} \qquad (A.1.1b) $$

A typical eddy diffusivity in the well-mixed convective boundary layer is estimated as $\kappa \sim 1 \text{ m}^2/\text{s}$ based on a typical turbulent kinetic energy of $10^{-1} \text{ m}^2/\text{s}^2$ along with standard parameters given in e.g., Redelsperger and Sommeria (1981) as well as the mixing length, $l = 10$ m. [NB; Holtslag and Boville (1993) assume $l = 30$ m as an asymptotic value for the free troposphere.] Also assuming a typical vertical scale to be $d \sim 10^3$ m, we find

$$ \Delta t \sim (10^3)^2/1 \text{ m}^2/\text{s} \sim 10^6 \text{ sec} \qquad (A.1.2) $$

which is about ten days, too slow to be realistic. To obtain a more realistic time scale of $\Delta t \sim 10^5$ day $\sim 1$ day, we need to assume rather a large eddy diffusivity, $\kappa \sim 10 \text{ m}^2/\text{s}$. [e.g., Holtslag and Boville (1993) show in their tests of the eddy-diffusion mode, the eddy diffusivity can go up to $30 \text{ m}^2/\text{s}$.]

The prescribed heating rate, $(\partial\theta/\partial t)_s$, must generate a well mixed layer of the depth about $d \sim 10$ m over $\delta t \sim 10^5$ sec $\sim 1$ day. Assuming that the well-mixed layer deepens as predicted from the well-mixed layer theory,

$$ \left(\frac{\partial\theta}{\partial t}\right)_s \sim \frac{d}{\Delta t} \frac{\partial\bar{\theta}}{\partial z} \sim \frac{10^3 \text{ m}}{10^5 \text{ s}} \times 3 \times 10^{-3} \text{ K/m} \sim 3 \times 10^{-5} \text{ K/s} \qquad (A.1.3a) $$

Alternatively, noting that $\partial\bar{\theta}/\partial z = \theta_*/d$ in terms of the characteristic scale, $\theta_*$, of the potential temperature, we find:

$$ \left(\frac{\partial\theta}{\partial t}\right)_s \sim \frac{\theta_*}{\Delta t} \qquad (A.1.3b) $$
where $\theta_* \sim 3$ K to be consistent with a typical stratification, $\partial\bar{\theta}/\partial z$, of the atmosphere.

On the other hand, to normalize the right-hand side of Eq. (3.5b), we need to assume:

$$ \frac{\kappa}{d^2\theta_*} \left(\frac{\partial\theta}{\partial t}\right)_s (\Delta t)^2 = 1 $$

By further substituting the relation (A.1.3b), we obtain:

$$ \frac{\kappa}{d^2\theta_*} \frac{\theta_*}{\Delta t} (\Delta t)^2 = 1 $$

which is equivalent to Eq. (A.1.1b). Thus, the normalization of the inhomogeneous term in Eq. (3.5b) is automatically accomplished by an already-introduced normalization.

## A.2 System Driven by Constant Heat Flux

### A.2.1 Formulation

A more realistic case to be considered is a system driven by a constant heat flux, $H$. However, in this case, the treatment of the system becomes slightly more involved. First, with the initial condition (2.2b) with a stable stratification, $\partial\bar{\theta}/\partial z > 0$, the flux at the surface, $H$, never transported upwards. To initiate this process, we need to assume an existence of an unstable surface layer of the depth, say, $\Delta z$. Over this layer, $\partial\bar{\theta}/\partial z < 0$ so that the heat can be transported upwards by the eddy diffusion. We can assume a constant stratification, say, $(\partial\bar{\theta}/\partial z)_b(< 0)$, and a stable stratification with $\partial\bar{\theta}/\partial z(> 0)$ is found from $z = \Delta z$ to the above.

As a result, we find a discontinuity of the potential temperature gradient at the height, $z = \Delta z$. This discontinuity is alleviated by integrating Eq. (3.1b) over a thin layer of $\varepsilon(\ll \Delta z)$ extending over $\Delta z - \varepsilon/2$ to $\Delta z + \varepsilon/2$ to obtain:

$$ 2\varepsilon \frac{\partial}{\partial t} \theta_b = \kappa \frac{\partial \langle\theta\rangle}{\partial z} \Big|_{\Delta z - \varepsilon/2}^{\Delta z + \varepsilon/2} $$

The subscript, $b$, suggests the value at $z = \Delta z$. The above can be re-arranged into:

$$ \frac{\partial}{\partial t} \theta_b = \frac{1}{\varepsilon} \left[ \kappa \left( \frac{\partial \langle \theta \rangle}{\partial z} \right)_{b+} - \kappa \left( \frac{\partial \bar{\theta}}{\partial z} \right)_b \right] $$

Here, the subscript $b+$ suggests a value just above the surface layer top. Further noting that $-\kappa(\partial\bar{\theta}/\partial z)_b = H$, the temporal tendency at $z = \Delta z$ is defined by:

$$ \frac{\partial \theta_b}{\partial t} = \frac{1}{2\varepsilon} \left[ \kappa \left( \frac{\partial \langle \theta \rangle}{\partial z} \right)_{b+} + H \right] \qquad (A.2.1a) $$
Thus, the problem to be solved is the evolution of the system (3.2b) with the boundary condition (A.2.1) given at $z = \Delta z$. Here, note further that

$$ \left( \frac{\partial \langle \theta \rangle}{\partial z} \right)_{b+} = \left( \frac{\partial \bar{\theta}}{\partial z} \right)_{b+} + \left( \frac{\partial \theta}{\partial z} \right)_{b+} \qquad (A.2.1b) $$

in which the first term on the right-hand side remains constant of time, but the second term evolves with time.

It transpires that the solution for Eq. (3.2b) can be sought in an analogous manner as in the main text, but by shifting the height coordinate by $\Delta z$ so that $z = \Delta z$ becomes $z = 0$. As the result, the procedure of solving the problem becomes identical to that of the main text, but changing the boundary condition from Eq. (2.2a) to Eq. (A.2.1a, b). Note importantly, the depth, $\Delta z$, of the surface layer no longer explicitly enters the problem as a result.

A first technical difficulty is to evaluate $\partial\theta/\partial z$ at $z = 0$, because the Fourier series defined by Eq. (3.3a) is singular at $z = 0$. The problem is mitigated by taking the position close enough to $z = 0$, but with a finite departure of, say, $\delta z$. As the result, we obtain

$$ \left( \frac{\partial \theta}{\partial z} \right)_{b+} = \sum_{n=1}^{N} k_n \theta_n \cos(k_n \delta z) \qquad (A.2.2) $$

In the following numerical presentations, we rather arbitrary sets, $\delta z = 1/N$.

### A.2.2 Typical Scales

By substituting the boundary condition (A.2.1a, b) into Eq. (3.2b), the problem to be solved becomes:

$$ \left( \frac{d}{dt} + k^2 \right) \theta_n = 2\delta k \left[ \mu \int_0^t \tilde{\kappa} \left( \frac{\partial \theta}{\partial z} \right)_{b+} dt + \lambda \left( \frac{1}{H_*} \frac{\partial \bar{\theta}}{\partial z_*} \int_0^t \tilde{\kappa} dt + t \right) \right] \qquad (A.2.3) $$

in place of Eq. (3.5b), after the nondimensionalizations. Here, the subscript, *, has been added to $H$ and $\partial\bar{\theta}/\partial z$ to suggest that they are dimensional, whereas ~ has been added to $\kappa$ for nondimensionality. The two nondimensional parameters have been introduced, which are defined by

$$ \mu = \frac{\kappa_*^2 (\Delta t)^2}{\varepsilon d^3} \qquad (A.2.4a) $$

$$ \lambda = \kappa_* \frac{H_* (\Delta t)^2}{\varepsilon \theta_* d^2} \qquad (A.2.4b). $$

where $\kappa = \kappa_* \tilde{\kappa}$. By further combining the definitions of those two parameters, we can estimate the vertical scale, $d$, to be

$$ d = \left( \frac{\lambda}{\mu} \right) \frac{\kappa_* \theta_*}{H} \qquad (A.2.4c) $$
Here, we first note that comparing the two boundary conditions in the main text and herein, we identify that the last term in Eq. (A.2.1a) becomes equivalent to Eq. (2.2a) by setting $H/\epsilon = (\partial\theta/\partial t)_s$. However, these two parameters take substantially different values. Taking a typical heat flux value, 30 W/m², we obtain $H \sim 3 \times 10^{-2}$ mK/s. Furthermore assuming $\epsilon \sim 0.1$ m, we find $H/\epsilon \sim 3 \times 10^{-1}$ K/s, which is substantially larger than the estimate surface temperature tendency by Eq. (A.1.2). It means that the surface layer is heated extremely intensively, and the excess of heat must be transported upwards efficiently through the first term on the right-hand side of Eq. (A.2.1a). Furthermore, to maintain the normalized state with $\lambda = 1$ in Eq. (A.2.4b), the compensating factor, $(\Delta t)^2/\theta_*d^2$, must be proportionally smaller, either by reducing $\Delta t$ or by increasing $d$.

We next note that the potential-temperature scale, $\theta_*$, and the vertical scale, $d$, are constrained by the lapse rate:

$$ \frac{\theta_*}{d} = \frac{\partial \bar{\theta}}{\partial z} \sim 3 \times 10^{-3} \text{ K/m} \qquad (A.2.5) $$

By substituting Eq. (A.2.5) into Eq. (A.2.4c), we find:

$$ \kappa_* = H \left( \frac{\partial \bar{\theta}}{\partial z} \right)^{-1} \sim 3 \times 10^{-2} \text{ mK/s} / 3 \times 10^{-3} \text{ K/s} \sim 10 \text{ m}^2/\text{s.} $$

This estimate of the eddy diffusivity is consistent with that of the main text.

More generally, by setting $\mu = 1$ and $\lambda = 1$ in Eqs. (A.2.4a, b), we find that the time and space scales are constrained by

$$ (\Delta t)^2 = \frac{\epsilon}{\kappa_*^2} d^3 \qquad (A.2.6a) $$

$$ (\Delta t)^2 = \frac{\epsilon \theta_*}{\kappa_* H} d^2 \qquad (A.2.6b). $$

where the involved factors are estimated as

$$ \frac{\epsilon}{\kappa_*^2} \sim \frac{0.1 \text{ m}}{(10 \text{ m}^2/\text{s})^2} \sim 10^{-3} \text{ s}^2/\text{m}^3 \qquad (A.2.7a) $$

$$ \frac{\theta_*}{\kappa_* H / \epsilon} \sim \frac{\theta_*}{10 \text{ m}^2/\text{s} \times 3 \times 10^{-1} \text{ K/s}} \sim \frac{\theta_*}{3 \text{ K}} \text{ s}^2/\text{m}^2 \qquad (A.2.7b). $$

by taking the standard values already introduced. Here, $\theta_*$ remains arbitrary.

When we set $\theta_* = 3$ K, from Eq. (A.2.6b), $\Delta t$ [s] = $d$ [m]. Substituting this result into Eq. (A.2.6a), we obtain $d = 10^3$ m and $\Delta t = 10^3$ sec. Thus, the same vertical scale as before, but much faster time scale. Note that with the decreasing $\theta_*$, the vertical scale,
$d$, also decreases proportionally by Eq. (A.2.5), and that also decrease the time scale, $\Delta t$, further.

More generally, from Eq. (A.2.7b),

$$ \theta_* = \kappa_* \frac{H}{\varepsilon} $$

thus the potential-temperature scale decreases with the increasing $\varepsilon$, keeping $\Delta t$ [s] = $d$ [m]. Furthermore, the increasing $\varepsilon$ decreases the characteristic scales, $d$ and $\Delta t$ from Eq. (A.2.7a). Alternatively, by holding $\theta_*$ constant and increasing $\varepsilon$, the time scale, $\Delta t$ can be increases relative to the length scale. More specifically, by setting $\varepsilon = 1$ m, we obtain a more realistic time scale of $\Delta t = 10^4$ sec.

### A.2.3 Solution

We now solve Eq. (A.2.3) by by the normalizations of $\mu = \lambda = (\partial\bar{\theta}/\partial z)_*/H_* = 1$, also setting $\tilde{\kappa} = 1$. Recall that the evolution of $(\partial\theta/\partial z)_{b+}$ is predicted by Eq. (A.2.2). The solution for Eq. (A.2.3) can be presented in an integral form, but no longer in a closed analytical form. The necessary time integral has been performed numerically with the Euler time stepping with the time step of $10^{-4}$.

<figure><img src="image_9.png" /><figcaption>Fig. A.1: The same as Fig. 1, but when a more sophisticated surface-layer description is included, but also with a short time interval of $\delta t = 0.1$.</figcaption></figure>

The resulting solutions are presented in Fig. A.1 in the same format as in Fig. 1. It is seen that under this formulation, the convective layer develops much faster than in the formulation presented in the main text, even after nondimensionalization with a shorter time scale.

## Appendix B: Discrete Sine Transform
When a function, $f(z)$, is defined at $N$ points defined at

$$z_j = \frac{jD}{N+1} \qquad (B.1)$$

for $j = 1, \dots, N$ over the domain, $0 \le z \le D$, by taking the analogy with the fast Fourier transform, the discrete sine transform is defined by

$$\hat{f}_n = \frac{2}{N} \sum_{j=1}^{N} f(z_j) \sin k_n z_j \qquad (B.2a)$$

where

$$k_n = \frac{n\pi}{D} \qquad (B.2b)$$

for $n = 1, \dots, N$. The Fourier series (3.3a) is also truncated into:

$$f(z) = \sum_{n=1}^{N} \hat{f}_n \sin k_n z_j \qquad (B.3)$$

## Appendix C: Second-Order Closure

The purpose of the Appendix C is to derive the closures, Eqs. (4.3) and (4.4a, b), for $\langle w'\theta' \rangle$, $E$, and $\langle \theta'^2 \rangle$. As seen in the following stepwise derivations, all those closures are derived from the corresponding budget equations of the given variables, but setting the temporal tendency terms to vanish. Thus, without further closure assumptions, there is no longer a term explicitly containing the budget variables. A key step to recover the dependence of the budgets on the budget variables is to assume that the dissipation rates of the given variables take forms of relaxation processes, thus they become $\langle w'\theta' \rangle / \tau_1$, $E/\tau_2$, and $\langle \theta'^2 \rangle / \tau_3$ where $\tau_1$, $\tau_2$, and $\tau_3$ are unspecified damping time scales. Since those closures are heuristic at the best, we need to see those closures with some cautions.

### C.1 Heat Flux

A complete heat-flux equation is given by

$$\frac{\partial}{\partial t} \langle w' \theta' \rangle + \frac{\partial}{\partial x_i} [\langle u_i \rangle \langle w' \theta' \rangle + \langle \theta' w' u'_i \rangle] + \langle \theta' u'_i \rangle \frac{\partial \langle w \rangle}{\partial x_i} + \langle w' u'_i \rangle \frac{\partial \langle \theta \rangle}{\partial x_i} = \\ \alpha g \langle \theta'^2 \rangle + \frac{1}{\rho} \langle p' \frac{\partial \theta'}{\partial z} \rangle + \frac{1}{\rho} \frac{\partial \langle p' \theta' \rangle}{\partial z} + \frac{\partial}{\partial x_i} [\nu \langle \theta' \frac{\partial w'}{\partial x_i} \rangle + \mu \langle w' \frac{\partial \theta'}{\partial x_i} \rangle] - (\nu + \mu) \langle \frac{\partial w'}{\partial x_i} \frac{\partial \theta'}{\partial x_i} \rangle \qquad (C.1.1)$$

(cf., Eq. 8 of Mellor 1973) in which Einstein's sum rule is applied.

We simplify this equation by assuming that: i) the higher-order terms are negligible, i.e., $\langle u_i \rangle \langle w' \theta' \rangle = 0$ and $\langle u_i \rangle \langle w' \theta' \rangle = 0$; ii) the pressure perturbation is not correlated with
the potential-temperature perturbation, i.e., $\langle p'\theta' \rangle = 0$. We further assume that non-isotropic eddy stress, $\tau_{ij}$, does not contribute, thus $\langle w'u'_i \rangle = \tau_{3i} + (2/3)\delta_{3i}E \simeq (2/3)\delta_{3i}E$, and the system is quasi-steady, thus the time derivative can be dropped, and $\partial\langle w'\theta' \rangle/\partial t = 0$. Here, the nonisotropic stress is defined by $\tau_{ij} = \langle u'_i u'_j \rangle - (2/3)\delta_{ij}E$, and $\delta_{ij}$ is the Kronecker's delta. Thus, the full budget (C.1.1) reduces to:

$$ -\frac{2}{3}E \frac{\partial \langle \theta \rangle}{\partial z} + \frac{1}{\rho} \langle p' \frac{\partial \theta'}{\partial z} \rangle + \alpha g \langle \theta'^2 \rangle = 0 \qquad (C.1.2) $$

Furthermore by introducing a generalization of Rotta's (1951) "return-to-isotropy" hypothesis to the pressure term, we set

$$ \frac{1}{\rho} \langle p' \frac{\partial \theta'}{\partial z} \rangle = -C_s \frac{E^{1/2}}{l} \langle w' \theta' \rangle - \frac{\alpha g}{3} \langle \theta'^2 \rangle \qquad (C.1.3) $$

Substitution of Eq. (C.1.3) into Eq. (C.1.2) leads to Eq. (4.3).

## C.2 Turbulent Kinetic Energy

A full equation for the turbulent kinetic energy, $E$, is:

$$ \frac{\partial E}{\partial t} + \frac{\partial}{\partial x_j} [\langle u_j \rangle E + \frac{1}{2} \langle u_i'^2 u_j' \rangle + \frac{1}{\rho} \langle p' u_j' \rangle] + \tau_{ij} \frac{\partial \langle u_i \rangle}{\partial x_j} = \alpha g \langle w' \theta' \rangle - \varepsilon_K \qquad (C.2.1a) $$

where

$$ \varepsilon_K = -\nu \left[ \nabla^2 E - \left\langle \left( \frac{\partial u_i}{\partial x_j} \right)^2 \right\rangle \right] \qquad (C.2.1b) $$

is the energy dissipation rate (cf., Eq. 71, Mellor 1973).

We simplify Eq. (C.2.1a) in an equivalent manner as for the heat equation: i) neglect the higher-order terms, i.e., $\langle u_j \rangle E = \langle u_i'^2 u_j' \rangle = 0$, ii) the pressure perturbation is not correlated with the eddy velocity, i.e., $\langle p' u'_i \rangle = 0$, and the quasi-steady state i.e., $\partial E/\partial t = 0$. Furthermore, assuming a closed form for the energy dissipation,

$$ \varepsilon_K = C_E \frac{E^{3/2}}{l} \qquad (C.2.2) $$

we arrive at:

$$ -\tau_{ij} \frac{\partial \langle u_i \rangle}{\partial x_j} + \alpha g \langle w' \theta' \rangle - C_E \frac{E^{3/2}}{l} = 0 \qquad (C.2.3) $$

We assume that the wind stress, $\tau_{ij}$, takes the form

$$ \tau_{ij} = -K_m \left( \frac{\partial \langle u_i \rangle}{\partial x_j} + \frac{\partial \langle u_j \rangle}{\partial x_i} \right), \qquad (C.2.4a) $$
where

$$K_m = \frac{4}{15} l \frac{E^{1/2}}{C_m}. \qquad (C.2.4b)$$

Substituting Eqs. (C.2.4a,b) into Eq. (C.2.3), and after re-arrangements, we obtain Eq. (4.4a).

### C.3 Potential–Temperature Variance

A full equation for the potential–temperature variance (cf., Eq. 9, Mellor 1973) is:

$$\frac{\partial}{\partial t} \langle \theta'^2 \rangle + \frac{\partial}{\partial x_i} [\langle u_i' \rangle \langle \theta'^2 \rangle + \langle u_i' \theta'^2 \rangle] = -2 \langle u_i' \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial x_i} - \chi \qquad (C.3.1a)$$

where the dissipation term, $\chi$, is defined by

$$\chi = -\mu[\nabla^2 \langle \theta'^2 \rangle - 2 \langle \left( \frac{\partial \theta'}{\partial x_i} \right)^2 \rangle] \qquad (C.3.1b)$$

As before, neglecting the higher–order terms and also assuming the quasi–steadiness, the above is simplified into:

$$2 \langle u_i' \theta' \rangle \frac{\partial \langle \theta \rangle}{\partial x_i} - \chi = 0 \qquad (C.3.2)$$

to which we introduce a closure

$$\chi = C_{\theta} \frac{E^{1/2}}{l} \langle \theta'^2 \rangle \qquad (C.3.3)$$

That leads to Eq. (4.4b).

### References

* André, J. C., G. De Moor, P. Lacarrére, G. Therry, and R. Du Vachat, 1978: Modeling the 24-hour evolution of the mean and turbulent structures of the planetary boundary layer. *J. Atmos. Sci.*, **35**, 1861–1883.

* Betts, A. K., 1974: Reply to comment on the paper ‘Non-precipitating cumulus convection and its parameterization’/ Quator. J. Roy. Meteor. Soc., **100**, 469–471.

* Blackadar, A., 1962: The vertical distribution of wind and turbulent exchange in a neutral atmosphere. *J. Geophys. Res.*, **67**, 3095–3102.

* Bougeault, P., and P. Lacarrére, 1989: Parameterization of orography-induced turbulence in a mesobeta-scale model. *Mon. Wea. Rev.*, **117**, 1872–1890, https://doi.org/10.1175/1520-0493(1989)117:1872:POOITI¿2.0.CO;2.

* Cuiipers, J. W. M., and A. A. M. Holtslag, 1998: Impact of skewness and nonlocal effects on scalar and buoyancy fluxes in convective boundary layer. *J. Atmos. Sci.*, **55**, 151–162.
Deardorff, J. W., 1970: Preliminary results from numerical integrations of the unstable planetary boundary layer. *J. Atmos. Sci.*, **27**, 1209-1211.

Deardorff, J. W., G. E. Willis, and D. K. Lilly, 1969: Laboratory investigation of non-steady penetrative convection. *J. Fluid Mech.*, **35**, 7-31.

Deardorff, J. W., G. E. Willis, and D. K. Lilly, 1974: Comment on the paper by A. K. Betts 'Non-precipitating cumulus convection and its parameterization' / Quator. J. Roy. Meteor. Soc., **100**, 122-123.

Durran, D. R., 1991: The third-order Adams-Bashforth–LMethod: An attractive alternative to leapfrog time differencing. *Mon. Wea. Rev.*, **119**, 702–720, https://doi.org/10.1175/1520-0493(1991)119;0702:TTOABM;2.0.CO;2.

Holtslag, A. A. M., and B. A. Boville, 1993: Local versus nonlocal boundary-layer diffusion in a global climate model. *J. Climate*, **6**, 1825–1842.

Lumely, J. L., 1975: Modeling turbulent flux of passive scalar quantities in inhomogeneous flows. *Phys. Fluids*, **18**, 619–621.

Mellor, G. L., 1973: Analytic prediction of the properties of stratified planetary surface layers. *J. Atmos. Sci.*, **30**, 1061–1069.

Redelsperger, J. L., and G. Sommeria, 1981: Methode de representation de la turbulence d'echelle inferiure q lq mqille pour un nmodelel tri-dimensionel de convection nuageuse. *Boundary-Layer Meteorol.*, **21**, 509-530.

Rotta, J. C., 1951: Statistische Theorie nichthomogenener Turbulenz. *Z. Phys.*, **129**, 547-572, **131**, 51-77.

Sommeria, G., 1976: Three-dimensional simulation of turbulent processes in an undisturbed trade wind boundary layer. *J. Atmos. Sci.*, **33**, 216–241.

Teixeria, J., 1999: Stable schemes for partial differential equations: The one-dimensional diffusion equation. *J. Comp. Phys*, **153**, 403-417.

Teixeria, J., J. P. Perreira, P. M. A. Miranda, T. Haack, J. Doyle, A. P. Siebesma, and R. Salgado, 2004: A new mixing-length formulation for the parameterization of dry convection: Implementation and evaluation in a mesoscale model. *Mon. Wea. Rev.*, **132**, 2698-2707.