A novel approach for deriving the stable boundary layer height and eddy viscosity profiles from the Ekman equations

In this study, we utilize a novel approach to solve the Ekman equations for eddy viscosity profiles in the stable boundary layer. By doing so, a well-known expression for the stable boundary layer height by Zilitinkevich (Boundary-Layer Meteorology, 1972, Vol. 3, 141--145) is rediscovered.


I. INTRODUCTION
Sergej Zilitinkevich was one of the giants of atmospheric physics who carried the field of boundary-layer meteorology (BLM) on his shoulder for more than half a century. We, representing the BLM community at large, are indebted forever for his ingenious efforts and lifelong dedication in advancing our field. In the arena of stable boundary layers (SBLs), Zilitinkevich made numerous ground-breaking contributions. As a matter of fact, it would be difficult to find any contemporary article on SBLs which does not make at least one reference to an original publication of Zilitinkevich. The present paper is also following the same tradition.
In this study, we utilize the Ekman equations to analytically derive a stable boundary layer height formulation which was originally proposed by Zilitinkevich [1] based on scaling arguments. During the process, we also derive equations for the vertical profiles of eddy viscosity. To the best of our knowledge, it is the first time that estimates are given for the eddy viscosity profiles directly from the Ekman equations.

II. FORMULATION OF SBL HEIGHT BY
ZILITINKEVICH (1972) Using boundary layer scaling arguments, Zilitinkevich [1] proposed that the height (h) of stationary SBLs can be written as: Here the surface friction velocity and surface Obukhov length are denoted by u * 0 and L, respectively; B s is the buoyancy flux. The Coriolis term is represented by f cor .
The proportionality constants γ and C h are related as follows: where κ is the von Kármán constant assumed to be equal to 0.4. Zilitinkevich [1] assumed C h to be order of one. In an analytical study, Businger and Arya [2] estimated γ to be equal to 0.72. A much lower value of γ ≈ 0.4 was estimated by Brost and Wyngaard [3] using a second-order closure model. Garratt [4] who analyzed observational data from several field campaigns also found γ ≈ 0.4. In his local-scaling paper, Nieuwstadt [5] found γ ≈ 0.35 to be consistent with other equations. Zilitinkevich [6] summarized C h from several studies and found it to vary within the range of 0.55-1.58. A unique aspect of the present study is that we analytically derive γ (and C h ) with limited assumptions. We would like to point out that Eq. 1 predicts physically unrealistic boundary layer heights for two situations: (i) close to the equator (i.e., f cor → 0); and (ii) for near-neutral conditions (i.e., L → ∞). To circumvent the second issue, a few interpolation approaches have been proposed in the literature [e.g., 6,7].

III. MOMENTUM AND SENSIBLE HEAT FLUX PROFILES
In SBL flows, the profiles of friction velocity and sensible heat flux are often expressed as follows: Here the friction velocity, momentum flux components, and sensible heat flux at height z are denoted by u * L , uw L , vw L , and wθ L , respectively. The subscript '0' is used to demarcate the corresponding surface values.
These equations imply that the turbulent fluxes are maximum near the surface and they monotonically decrease to zero at the top of the boundary layer. The exponents α and β in Eqs. 3a and 3b are not universal constants. By utilizing his local-scaling hypothesis, Nieuwstadt [5] suggested α = 1.5. He also proved that β should be equal to 1 under the assumptions of horizontal homogeneity and stationarity. In contrast, based on observational data from the well-known Minnesota field campaign [8], Sorbjan [9] estimated α = 2, and β = 3 for evolving SBLs. In another observational study, Lenschow et al. [10] considered the additional effects of radiational cooling and found the optimal α and β to be equal to 1.75 and 1.5, respectively.
Using the definition of Obukhov length [11] in conjunction with Eqs. 3a and 3b, we get: where Λ is the so-called local Obukhov length at height z.
The first and second derivatives of f m function in Eq. 3a can be written as: We will make use of these derivatives in a later section.

IV. EDDY VISCOSITY PROFILE IN THE EKMAN LAYER
According to the K-theory, based on the celebrated hypothesis of Boussinesq in 1877, turbulent fluxes can be approximated as products of the eddy exchange coefficients and the mean gradients [12].
Here K M is the so-called eddy viscosity coefficient. Based on these equations, we can write: where S is the magnitude of wind speed shear. For boundary layer flows over homogeneous and flat terrain, under steady-state conditions, the averaged equations of motions can be simplified as follows [13,14]: The velocity components in x and y directions are represented as U and V , respectively. Similarly, the geostrophic velocity components are denoted by U g and V g . In the literature, these equations are commonly known as the Ekman layer equations [14,15]. In a landmark paper, Ekman [16] first analytically solved Eqs. 8a and 8b with the assumption of constant eddy viscosity (K M ) in conjunction with appropriate boundary conditions. Over the years, a few more closed form analytical solutions of the Ekman equations have been reported in the literature [e.g., 14,15,[17][18][19]. In all these papers, simplified profiles of K M were always assumed. In this paper, we take a radically different approach. We only assume that Eq. 3a is valid and then deduce K M profile from the Ekman equations as shown below.
For barotropic condition, the Ekman equations can be re-written as: Next, by utilizing Eqs. 6a and 6b, these equations are transformed as follows: Dividing Eq. 10a by Eq. 10b and rearranging we arrive at: The momentum flux components can be decomposed in terms of local friction velocity as follows: Here δ is the angle between the flux vector and the xaxis. By plugging Eq. 3a in Eqs. 12a and 12b, we get: These equations can be differentiated as follows: After differentiating one more time we arrive at: By combining Eqs. 11, 13a, 13b, 15a, 15b, and simplifying we get: By invoking the Pythagorean trigonometric identity, we can further simplify this equation to: Thus, Substituting f m and f ′′ m from Eqs. 3a and 5b in 18, we find: The second derivative of δ is: Please note that these derivatives have real values if α is greater than one. Now, we can directly estimate the K M profile from Eqs. 10a, 13b, 15a, 19, and 20: We would like to emphasize that this formulation of K M is derived directly from the Ekman equations with very limited assumptions. To the best of our knowledge, a similar formulation and derivation have not been reported in the literature.
Similar to other Ekman layer findings, Eq. 21 is only valid in the outer layer. It does not represent surface layer conditions. In the literature, various patching and asymptotic matching approaches [e.g., 14,[20][21][22] have been proposed to combine outer layer and inner layer (i.e., surface layer) solutions. A nice overview was given by [23]. In Sects. VI and VII, we utilize an unorthodox strategy.

V. CONVENTIONAL K-PROFILE APPROACH
One of the most widely used first-order formulation for K M is the K-profile approach [24]. O'Brien [25] was one of the first researchers to propose a K-profile which portrays desirable surface layer behavior, attains a maximum value within the planetary boundary layer (PBL), and decreases to a background diffusion level above the PBL. Based on a second-order closure model, Brost and Wyngaard [3] proposed a different K-profile formulation for stably stratified flows: Here φ M is a type of non-dimensional velocity gradient, defined later. The exponent p is a-priori not known. For neutrally stratified flows in the surface layer, Eq. 22 re-duces to K M = κzu * ; this equation is in complete agreement with the well-established logarithmic law of the wall. Furthermore, for stably stratified surface layers, one can deduce a stability-corrected logarithmic law of the wall [e.g., 26] from Eq. 22. The K-profile approach was modified by several researchers [27][28][29][30][31][32] for its application in the unstable regime. They included a counter-gradient term to include the effects of large-scale eddies on local fluxes. They also considered entrainment fluxes at the top of the PBL.
In the context of numerical stability, the K-profile approach is quite robust [33]. Thus, it is not surprising that it is widely used in numerical weather prediction models. As a matter of fact, the default PBL scheme (called the YSU scheme) of the popular Weather Research and Forecasting (WRF) model uses the K-profile approach for both unstable and stable conditions.
In spite of its wide usage, Eq. 22 suffers from two limitations. First, there is uncertainty in the value of the exponent p. Brost and Wyngaard [3] found p to be equal to 1.5 based on their simulations. On the other hand, based on field campaign data from Minnesota, Sorbjan [34] found p = 1. In the absence of reliable field observations Troen and Mahrt [32] used an integer value of p = 2. The local scaling hypothesis by Nieuwstadt [5] also leads to p = 2. In the present study, we estimate p analytically.
The other limitation is related to the parameterization of φ M . The K-profile formulation uses the following normalized velocity gradient: Please note that in this equation surface friction velocity (u * 0 ) is used. However, in contrast to well-known surface layer formulations, z is not limited to the depth of the surface layer. Instead, z ranges from the surface to the top of the boundary layer. φ M is commonly parameterized as [3]: where c is a constant and often assumed equal to 5. In the surface layer, for z/L < 1, numerous studies documented the validity of this equation. However, its applicability for the outer layer (i.e., above surface layer) is questionable. Furthermore, this equation (incorrectly) implies that the logarithmic law of the wall applies to the entire boundary layer for neutral condition (i.e., h/L = 0).

VI. ALTERNATIVE K-PROFILE APPROACHES
In this section, we derive two competing K-profile formulations. Both the formulations are applicable for the entire SBL (i.e., including the surface layer and the outer layer).

A. Option 1
Multiplying both sides of Eq. 7 by (κz/u * 0 ), we get: By using Eq. 3a, we deduce: Thus, This equation is identical to the one proposed by Brost and Wyngaard [3] based on second-order modeling. Except, in case of Eq. 24c, the exponent α is the same as in Eq. 3a.

B. Option 2
An alternate expression for K M profile can be found by multiplying Eq. 7 by (κz/u * L ): or, (κzu * L ) = K M φ ML .
By using Eq. 3a, we get: Hence, In these equations, φ ML (= κzS/u * L ) is a local nondimensional velocity gradient as it utilizes local friction velocity (u * L ) from height z. It is straightforward to show that: Equation 26 implies that φ M decreases more strongly with height than φ ML . Several past simulation studies [e.g., [35][36][37] found the following parameterization for φ ML : In those studies, c L was found to be between 3 and 5. By regression analysis of field observations, Mahrt and Vickers [38] found c L = 3.7. Here, we assume c L = 4.

VII. MATCHING OF KM PROFILES IN THE OUTER LAYER
Thus far, we have derived 3 different K-profiles. Note that Eq. 21 is only valid in the outer layer. In contrast, Eqs. 24c and 25d are valid in the entire SBL. In the following sub-sections, we match these K M profiles for the outer layer when z/L ≫ 1 and z/Λ ≫ 1.