1 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 (1972) 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.

2 Formulation of SBL Height by Zilitinkevich (1972)

Using boundary layer scaling arguments, Zilitinkevich (1972) proposed that the height (h) of stationary SBLs can be written as:

$$\begin{aligned} h = \gamma \sqrt{\frac{u_{*0} L}{|f_{cor}|}} = C_h u_{*0}^2 \left| f_{cor} B_s \right| ^{-1/2}. \end{aligned}$$
(1)

Here, the surface friction velocity and surface Obukhov length are denoted by \(u_{*0}\) and L, respectively; \(B_s\) is the surface buoyancy flux. The Coriolis parameter is represented by \(f_{cor} = \text {sgn}(f_{cor}) |f_{cor}|\). In the northern and southern hemispheres the value of \(\text {sgn}(f_{cor})\) is + 1 and − 1, respectively (Garratt 1994). The proportionality constants \(\gamma \) and \(C_h\) are related as follows:

$$\begin{aligned} C_h = \frac{\gamma }{\sqrt{\kappa }} \approx 1.58 \gamma , \end{aligned}$$
(2)

where \(\kappa \) is the von Kármán constant assumed to be equal to 0.4.

Zilitinkevich (1972) assumed \(C_h\) to be order of one. In an analytical study, Businger and Arya (1974) estimated \(\gamma \) to be equal to 0.72. A much lower value of \(\gamma \approx 0.4\) was estimated by Brost and Wyngaard (1978) using a second-order closure model. Garratt (1982) who analyzed observational data from several field campaigns also found \(\gamma \approx 0.4\). In his local-scaling paper, Nieuwstadt (1984) found \(\gamma \approx 0.35\) to be consistent with other equations. Zilitinkevich (1989) 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 \(\gamma \) (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} \rightarrow 0\)); and (ii) for near-neutral conditions (i.e., \(L \rightarrow \infty \)). To circumvent the second issue, a few interpolation approaches have been proposed in the literature (e.g.,Nieuwstadt 1981; Zilitinkevich 1989).

3 Momentum and Sensible Heat Flux Profiles

In SBL flows, the profiles of friction velocity and sensible heat flux are often expressed as follows (Nieuwstadt 1984; Sorbjan 1986):

$$\begin{aligned}{} & {} u_{*L}^2 = \left( \overline{uw}_L^2 + \overline{vw}_L^2\right) ^{1/2} = u_{*0}^2 \left( 1 - \frac{z}{h} \right) ^{\alpha } = u_{*0}^2 f_m, \end{aligned}$$
(3a)
$$\begin{aligned}{} & {} \overline{w\theta }_L = \overline{w\theta }_0 \left( 1- \frac{z}{h} \right) ^{\beta } = \overline{w\theta }_0 f_h. \end{aligned}$$
(3b)

Here, the friction velocity, momentum flux components, and sensible heat flux at height z are denoted by \(u_{*L}\), \(\overline{uw}_L\), \(\overline{vw}_L\), and \(\overline{w\theta }_L\), respectively. The subscript ‘0’ is used to demarcate the corresponding surface values. These equations imply that the magnitudes of turbulent fluxes are maximum near the surface and they monotonically decrease to zero at the top of the boundary layer.

The exponents \(\alpha \) and \(\beta \) in Eqs. 3a and 3b are not universal constants. By utilizing his local-scaling hypothesis, Nieuwstadt (1984) suggested \(\alpha = 1.5\). He also proved that \(\beta \) 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 (Caughey et al. 1979), Sorbjan (1986) estimated \(\alpha = 2\), and \(\beta = 3\) for evolving SBLs. In another observational study, Lenschow et al. (1988) considered the additional effects of radiational cooling and found the optimal \(\alpha \) and \(\beta \) to be equal to 1.75 and 1.5, respectively.

Using the definition of Obukhov length (Stull 1988) in conjunction with Eqs. 3a and 3b, we get:

$$\begin{aligned} \Lambda = L \frac{f_m^{3/2}}{f_h}, \end{aligned}$$
(4)

where \(\Lambda \) 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:

$$\begin{aligned}{} & {} f_m' = - \frac{\alpha }{h} \left( 1 - \frac{z}{h}\right) ^{\alpha -1}, \end{aligned}$$
(5a)
$$\begin{aligned}{} & {} f_m'' = \frac{\alpha \left( \alpha - 1\right) }{h^2} \left( 1 - \frac{z}{h}\right) ^{\alpha -2}. \end{aligned}$$
(5b)

We will make use of these derivatives in a later section.

4 Eddy-Viscosity Profile in the Ekman Layer

According to the K-theory, based on the celebrated hypothesis of Boussinesq (1877), turbulent fluxes can be approximated as products of the eddy exchange coefficients and the mean gradients (Lumley and Panofsky 1964):

$$\begin{aligned}{} & {} \overline{uw}_L = - K_M \frac{\partial U}{\partial z}, \end{aligned}$$
(6a)
$$\begin{aligned}{} & {} \overline{vw}_L = - K_M \frac{\partial V}{\partial z}. \end{aligned}$$
(6b)

here \(K_M\) is the so-called eddy-viscosity coefficient. The mean velocity components in x and y directions are denoted as U and V, respectively. Based on these equations, we can write:

$$\begin{aligned} u_{*L}^2 = \left( \overline{uw}_L^2 + \overline{vw}_L^2\right) ^{1/2} = K_M S, \end{aligned}$$
(7)

where S is the magnitude of wind velocity shear.

For boundary layer flows over homogeneous and flat terrain, under steady-state conditions, the averaged equations of motions can be simplified as follows (Zilitinkevich et al. 1967; Brown 1974; Garratt 1994):

$$\begin{aligned}{} & {} \text {sgn}(f_{cor}) |f_{cor}| \left( V - V_g\right) = \frac{\partial \left( \overline{uw}_L\right) }{\partial z}, \end{aligned}$$
(8a)
$$\begin{aligned}{} & {} \text {sgn}(f_{cor})|f_{cor}| \left( U - U_g\right) = - \frac{\partial \left( \overline{vw}_L\right) }{\partial z}. \end{aligned}$$
(8b)

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 (Brown 1974; Nieuwstadt 1983).

In a landmark paper, Ekman (1905) 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.,Wippermann 1973; Brown 1974; Nieuwstadt 1983; Grisogono 1995; Parmhed et al. 2005). 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 conditions, the Ekman equations can be re-written as:

$$\begin{aligned}{} & {} \text {sgn}(f_{cor})|f_{cor}| \frac{\partial V}{\partial z} = \frac{\partial ^2 \left( \overline{uw}_L\right) }{\partial z^2}, \end{aligned}$$
(9a)
$$\begin{aligned}{} & {} \text {sgn}(f_{cor})|f_{cor}| \frac{\partial U}{\partial z} = - \frac{\partial ^2 \left( \overline{vw}_L\right) }{\partial z^2}. \end{aligned}$$
(9b)

Next, by utilizing Eqs. 6a and 6b, these equations are transformed as follows:

$$\begin{aligned}{} & {} - \text {sgn}(f_{cor})|f_{cor}| \frac{\overline{vw}_L}{K_M} = \frac{\partial ^2 \left( \overline{uw}_L\right) }{\partial z^2}, \end{aligned}$$
(10a)
$$\begin{aligned}{} & {} \quad - \text {sgn}(f_{cor}) |f_{cor}| \frac{\overline{uw}_L}{K_M} = - \frac{\partial ^2 \left( \overline{vw}_L\right) }{\partial z^2}. \end{aligned}$$
(10b)

Dividing Eq. 10a by Eq. 10b and rearranging we arrive at:

$$\begin{aligned}{} & {} \overline{uw}_L \frac{\partial ^2 \left( \overline{uw}_L\right) }{\partial z^2} + \overline{vw}_L \frac{\partial ^2 \left( \overline{vw}_L\right) }{\partial z^2} = 0. \end{aligned}$$
(11)

The momentum flux components can be decomposed in terms of local friction velocity as follows (Businger and Arya 1974):

$$\begin{aligned}{} & {} \overline{uw}_L = - u_{*L}^2 \cos \left( \delta \right) , \end{aligned}$$
(12a)
$$\begin{aligned}{} & {} \overline{vw}_L = \text {sgn}(f_{cor}) u_{*L}^2 \sin \left( \delta \right) . \end{aligned}$$
(12b)

Here, \(\delta \) is the angle between the flux vector and the x-axis. By plugging Eq. 3a in Eqs. 12a and 12b, we get:

$$\begin{aligned}{} & {} \overline{uw}_L = - u_{*0}^2 f_m \cos \left( \delta \right) , \end{aligned}$$
(13a)
$$\begin{aligned}{} & {} \overline{vw}_L = \text {sgn}(f_{cor}) u_{*0}^2 f_m \sin \left( \delta \right) . \end{aligned}$$
(13b)

These equations can be differentiated as follows:

$$\begin{aligned}{} & {} \frac{\partial \left( \overline{uw}_L\right) }{\partial z} = - u_{*0}^2 \left[ f'_m \cos (\delta ) - f_m \delta ' \sin (\delta ) \right] , \end{aligned}$$
(14a)
$$\begin{aligned}{} & {} \frac{\partial \left( \overline{vw}_L\right) }{\partial z} = \text {sgn}(f_{cor}) u_{*0}^2 \left[ f'_m \sin (\delta ) + f_m \delta ' \cos (\delta ) \right] . \end{aligned}$$
(14b)

After differentiating one more time we arrive at:

$$\begin{aligned}{} & {} \frac{\partial ^2 \left( \overline{uw}_L\right) }{\partial z^2} = - u_{*0}^2 \left[ f''_m \cos (\delta ) - 2 f'_m \delta ' \sin (\delta ) - f_m \delta '' \sin (\delta ) -f_m \delta '^2 \cos (\delta )\right] , \end{aligned}$$
(15a)
$$\begin{aligned}{} & {} \frac{\partial ^2 \left( \overline{vw}_L\right) }{\partial z^2} = \text {sgn}(f_{cor}) u_{*0}^2 \left[ f''_m \sin (\delta ) + 2 f'_m \delta ' \cos (\delta ) + f_m \delta '' \cos (\delta ) - f_m \delta '^2 \sin (\delta )\right] . \nonumber \\ \end{aligned}$$
(15b)

By combining Eqs. 11, 13a, 13b, 15a, 15b, and simplifying we get:

$$\begin{aligned} f_m f''_m \cos ^2(\delta ) - 2 f_m f'_m \delta ' \sin (\delta )\cos (\delta ) - f_m^2 \delta '' \sin (\delta )\cos (\delta ) - f_m^2 \delta '^2 \cos ^2(\delta ) + \nonumber \\ f_m f''_m \sin ^2(\delta ) + 2 f_m f'_m \delta ' \sin (\delta )\cos (\delta ) + f_m^2 \delta '' \sin (\delta )\cos (\delta ) - f_m^2 \delta '^2 \sin ^2(\delta ) = 0. \end{aligned}$$
(16)

Please note that we have used \(\left[ \text {sgn}(f_{cor})\right] ^2 = 1\).

By invoking the Pythagorean trigonometric identity, we can further simplify this equation to:

$$\begin{aligned} f_m f''_m - f^2_m \delta '^2 = 0. \end{aligned}$$
(17)

Thus,

$$\begin{aligned} \delta ' = \sqrt{\frac{f''_m}{f_m}}. \end{aligned}$$
(18)

Substituting \(f_m\) and \(f_m''\) from Eqs. 3a and 5b in 18, we find:

$$\begin{aligned} \delta ' = \left( \frac{\sqrt{\alpha \left( \alpha -1 \right) }}{h}\right) \left( 1-\frac{z}{h}\right) ^{-1}. \end{aligned}$$
(19)

The second derivative of \(\delta \) is:

$$\begin{aligned} \delta '' = \frac{\sqrt{\alpha \left( \alpha -1 \right) }}{h^2} \left( 1-\frac{z}{h}\right) ^{-2}. \end{aligned}$$
(20)

Please note that these derivatives have real values if \(\alpha \) is greater than one.

Now, we can directly estimate the \(K_M\) profile from Eqs. 10a, 13b, 15a, 19, and 20:

$$\begin{aligned} K_M&= - \frac{\text {sgn}(f_{cor})|f_{cor}| \overline{vw}_L}{\partial ^2 (\overline{uw}_L)/ \partial z^2} \end{aligned}$$
(21a)
$$\begin{aligned}&= + \frac{\left[ \text {sgn}(f_{cor})\right] ^2 |f_{cor}| f_m \sin (\delta )}{f''_m \cos (\delta ) - 2 f'_m \delta ' \sin (\delta ) - f_m \delta '' \sin (\delta ) - f_m \delta '^2 \cos (\delta )} \end{aligned}$$
(21b)
$$\begin{aligned}&= \frac{|f_{cor}| h^2 }{\left( 2 \alpha - 1 \right) \sqrt{\alpha \left( \alpha - 1 \right) } }\left( 1 - \frac{z}{h} \right) ^2. \end{aligned}$$
(21c)

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, 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.,Taylor 1915; Blackadar and Tennekes 1968; Brown 1974; Zilitinkevich 1975) have been proposed to combine outer layer and inner layer (i.e., surface layer) solutions. A nice overview was given by Hess and Garratt (2002). In Sects. 6 and 7, we utilize an unorthodox strategy.

5 Conventional K-profile Approach

One of the most widely used first-order formulation for \(K_M\) is the K-profile approach (Stensrud 2007). O’Brien (1970) 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 (1978) proposed a different K-profile formulation for stably stratified flows:

$$\begin{aligned} K_M = \frac{\left( \kappa z u_{*0} \right) }{\phi _M} \left( 1 - \frac{z}{h} \right) ^{p}. \end{aligned}$$
(22)

Here, \(\phi _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 reduces to \(K_M = \kappa z u_*\); 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.,Businger et al. 1971) from Eq. 22.

The K-profile approach was modified by Troen and Mahrt (1986), Holtslag and Moeng (1991), Holtslag and Boville (1993), Hong and Pan (1996), Noh et al. (2003), Hong et al. (2006), and other researchers 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. Holtslag and Moeng (1991) estimated p to be equal to 2 for unstable conditions.

In the context of numerical stability, the K-profile approach is quite robust (Beljaars 1992). 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 for stable conditions. First, there is uncertainty in the value of the exponent p. Brost and Wyngaard (1978) found p to be equal to 1.5 based on their simulations. On the other hand, based on field campaign data from Minnesota, Sorbjan (1989) found \(p = 1\). In the absence of reliable field observations, Troen and Mahrt (1986) used an integer value of \(p = 2\). The local scaling hypothesis by Nieuwstadt (1984) also leads to \(p = 2\). In the present study, we estimate p analytically.

The other limitation is related to the parameterization of \(\phi _M\). The K-profile formulation uses the following normalized velocity gradient:

$$\begin{aligned}{} & {} \phi _M = \left( \frac{\kappa z S}{u_{*0}} \right) . \end{aligned}$$
(23a)

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. Function \(\phi _M\) is commonly parameterized as (e.g.,Brost and Wyngaard 1978):

$$\begin{aligned}{} & {} \phi _M = 1 + c \left( \frac{z}{L} \right) = 1 + c \left( \frac{z}{h} \right) \left( \frac{h}{L} \right) , \end{aligned}$$
(23b)

where c is a constant 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 conditions (i.e., \(h/L = 0\)).

6 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).

6.1 Option 1

Multiplying both sides of Eq. 7 by \(\left( \kappa z/u_{*0}\right) \), we get:

$$\begin{aligned} \left( \frac{\kappa z}{u_{*0}} \right) u_{*L}^2 = K_M \left( \frac{\kappa z S}{u_{*0}} \right) , \end{aligned}$$
(24a)

By using Eq. 3a, we deduce:

$$\begin{aligned} \left( \kappa z u_{*0} \right) f_m = K_M \phi _M. \end{aligned}$$
(24b)

Thus,

$$\begin{aligned} K_M = \frac{\left( \kappa z u_{*0} \right) f_m}{\phi _M} = \frac{\left( \kappa z u_{*0} \right) }{\phi _M} \left( 1 - \frac{z}{h} \right) ^{\alpha }. \end{aligned}$$
(24c)

This equation is identical to the one proposed by Brost and Wyngaard (1978) based on the second-order turbulence modeling. Except, in case of Eq. 24c, the exponent \(\alpha \) is the same as in Eq. 3a.

6.2 Option 2

An alternate expression for \(K_M\) profile can be found by multiplying Eq. 7 by \((\kappa z/u_{*L})\):

$$\begin{aligned} \left( \frac{\kappa z}{u_{*L}} \right) u_{*L}^2 = K_M \left( \frac{\kappa z S}{u_{*L}} \right) , \end{aligned}$$
(25a)

or

$$\begin{aligned} \left( \kappa z u_{*L} \right) = K_M \phi _{ML}. \end{aligned}$$
(25b)

By using Eq. 3a, we get:

$$\begin{aligned} \left( \kappa z u_{*0} \right) f_m^{1/2} = K_M \phi _{ML}. \end{aligned}$$
(25c)

Hence,

$$\begin{aligned} K_M = \frac{\left( \kappa z u_{*0} \right) f_m^{1/2}}{\phi _{ML}} = \frac{\left( \kappa z u_{*0} \right) }{\phi _{ML}} \left( 1 - \frac{z}{h} \right) ^{\alpha /2}. \end{aligned}$$
(25d)

In these equations, \(\phi _{ML} (= \kappa z S/u_{*L})\) is a local non-dimensional velocity gradient as it utilizes local friction velocity (\(u_{*L}\)) from height z. It is straightforward to show that:

$$\begin{aligned} \phi _M = \phi _{ML} f_m^{1/2}. \end{aligned}$$
(26)

Equation 26 implies that \(\phi _M\) decreases more strongly with height than \(\phi _{ML}\). Several past simulation studies (e.g., Basu and Porté-Agel 2006; Zilitinkevich and Esau 2007; van de Wiel et al. 2008) found the following parameterization for \(\phi _{ML}\):

$$\begin{aligned} \phi _{ML} = 1 + c_L \left( \frac{z}{\Lambda }\right) . \end{aligned}$$
(27)

In those studies, \(c_L\) was found to be between 3 and 5. By regression analysis of field observations, Mahrt and Vickers (2003) found \(c_L\) = 3.7. Here, we assume \(c_L = 4\).

7 Matching of \(K_M\) 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 \gg 1\) and \(z/\Lambda \gg 1\).

7.1 Option 1

In the outer layer, for \(z/L \gg 1\), Eq. 24c simplifies to:

$$\begin{aligned} K_M&= \frac{\left( \kappa z u_{*0} \right) f_m}{c \left( \frac{z}{L}\right) }, \end{aligned}$$
(28a)
$$\begin{aligned}&= \left( \frac{\kappa }{c}\right) \left( u_{*0} L \right) f_m, \end{aligned}$$
(28b)
$$\begin{aligned}&= \left( \frac{\kappa }{c}\right) \left( u_{*0} L\right) \left( 1 - \frac{z}{h} \right) ^{\alpha }. \end{aligned}$$
(28c)

If we assume \(\alpha = 2\), then we get:

$$\begin{aligned} K_M = \left( \frac{\kappa }{c}\right) \left( u_{*0} L \right) \left( 1 - \frac{z}{h} \right) ^2. \end{aligned}$$
(29)

Similarly, by plugging in \(\alpha = 2\) in Eq. 21, we get:

$$\begin{aligned} K_M = \frac{|f_{cor}| h^2}{3\sqrt{2}} \left( 1 - \frac{z}{h} \right) ^2. \end{aligned}$$
(30)

By directly matching Eq. 29 with Eq. 30 we arrive at:

$$\begin{aligned} \frac{|f_{cor}| h^2}{3\sqrt{2}} = \left( \frac{\kappa }{c}\right) \left( u_{*0} L \right) . \end{aligned}$$
(31)

Further simplification leads to the SBL height parameterization of Zilitinkevich (1972):

$$\begin{aligned} h = \gamma _1 \sqrt{\frac{u_{*0} L}{|f_{cor}|}}, \end{aligned}$$
(32a)

where:

$$\begin{aligned} \gamma _1 = \sqrt{3 \sqrt{2} \left( \frac{\kappa }{c}\right) }. \end{aligned}$$
(32b)

For \(c = 5\), we have: \(\gamma _1 = 0.583\).

7.2 Option 2

Similar to option 1, in the outer layer, for \(z/\Lambda \gg 1\), Eq. 25d simplifies to:

$$\begin{aligned} K_M&= \frac{\left( \kappa z u_{*0} \right) f_m^{1/2}}{c_L \left( \frac{z}{\Lambda }\right) }, \end{aligned}$$
(33a)
$$\begin{aligned}&= \left( \frac{\kappa }{c_L}\right) \left( u_{*0} \Lambda \right) f_m^{1/2}, \end{aligned}$$
(33b)
$$\begin{aligned}&= \left( \frac{\kappa }{c_L}\right) \left( \frac{u_{*0} L f_m^{2}}{f_h} \right) . \end{aligned}$$
(33c)

In this derivation, we used Eq. 4 in the conversion of \(\Lambda \) to L. If we assume \(\alpha = 3/2\) and \(\beta = 1\), we get:

$$\begin{aligned} K_M = \left( \frac{\kappa }{c_L}\right) \left( u_{*0} L \right) \left( 1 - \frac{z}{h} \right) ^2. \end{aligned}$$
(34)

Comparing this equation with Eq. 21, we arrive at:

$$\begin{aligned} h = \gamma _2 \sqrt{\frac{u_{*0} L}{|f_{cor}|}}, \end{aligned}$$
(35a)

where:

$$\begin{aligned} \gamma _2 = \sqrt{\sqrt{3} \left( \frac{\kappa }{c_L}\right) }. \end{aligned}$$
(35b)

If \(c_L = 4\), we get \(\gamma _2 = 0.416\).

Interestingly, both the options 1 and 2 lead to the SBL height parameterization of Zilitinkevich (1972). The estimated proportionality constants, \(\gamma _1\) and \(\gamma _2\), fall within the range of previous observation-based and simulation-based empirical values. Furthermore, based on the literature, both \(\alpha = 3/2\) and \(\alpha = 2\) are plausible in SBLs. However, using \(\phi _M = 1 + 5 z/L\) for the entire SBL, as commonly done in practice, does not seem physically meaningful and is in contrast with observations (e.g, Holtslag 1984, among many others). From this perspective, option 2 seems to be a better option.

8 Conclusion

Fifty years ago, Zilitinkevich (1972) proposed a formulation for the SBL height by using boundary-layer scaling arguments. In this study, we derive the same formulation from an analytical approach involving the Ekman layer equations. In addition, we provide novel derivations for eddy-viscosity profiles in the SBL. Our approach makes use of the following assumptions: (i) the magnitude of turbulent fluxes decrease monotonically with height and go to zero at the top of the boundary layer; (ii) the K theory is applicable for stably stratified conditions; and (iii) under steady-state condition, for flows over homogeneous and flat terrain, geostrophic balance holds in the boundary layer. None of these assumptions are unorthodox. In our future work, we hope to extend our analytical approach to derive geostrophic drag laws for SBLs.