It is desirable that a surface layer scheme in an atmospheric numerical model is consistent with an atmospheric boundary layer scheme incorporated in the same model. In this study, stability functions based on Monin–Obukhov similarity theory for momentum and heat, φ_m and φ_h, in the stable surface layer are derived from the Mellor–Yamada–Nakanishi–Niino (MYNN) scheme. The resulting stability functions are approximated by φ_m = 1 + 4.8z/L and φ_h = 0.74 + 6.0z/L, which can be analytically integrated with respect to height z to obtain momentum and heat fluxes, where L is the Obukhov length. The fluxes thus, obtained are compared with those acquired from stability functions in four previous studies: they are nearly the same as those from two of them, and show better agreement with observational data of the Surface Heat Budget of the Arctic Ocean experiment (SHEBA) over sea ice than those from the other two studies. Detailed comparisons of the results of the MYNN scheme with the SHEBA data suggest that significant variations of the fluxes observed during “winter” when the ice was covered with dry snow may have been caused by those of the surface roughness around the observational site. The stability functions obtained from the MYNN scheme predict that the bulk and flux Richardson numbers approach critical values of 0.26 and 0.21, respectively, in the limit of z/L → ∞. These critical values result from the Kolmogorov hypothesis applied to the turbulent dissipation in the MYNN scheme and are considered to correspond to a transition from Kolmogorov to non-Kolmogorov turbulence.

Weather and climate models often fail to accurately predict near-surface variables over nocturnal land and ice (e.g., Holtslag et al. 2013), implying that they still have challenges in predicting momentum and heat transfers between the surface and the atmosphere, especially in a stable condition. The element that plays the most important role in such predictions is a surface layer scheme. Most of the previously proposed surface layer schemes are based on Monin–Obukhov similarity theory (MOST). MOST defines stability functions (dimensionless gradient functions) for momentum and heat, φ_m and φ_h, which are functions of dimensionless height ζ scaled by the Obukhov length, but it does not provide their functional forms. Thus, various forms of the stability functions have been proposed based on field measurements, laboratory experiments, and theoretical arguments (e.g., Dyer 1974; Beljaars and Holtslag 1991; Cheng and Brutsaert 2005; Gryanik et al. 2020). Their functional forms, however, show considerable differences, especially for stable stratification, reflecting difficulties in measurements in the stable surface layer.

The surface layer scheme gives turbulent fluxes at the surface as the lower boundary conditions to numerical models. If, for example, downward heat fluxes given by a surface layer scheme are larger than those expected from an atmospheric boundary layer (ABL) scheme, turbulent transports could cause artificially excessive cooling of the lowest layer above the surface. Such cooling would then increase the vertical gradient of temperature (the stability) between the first and second layers above the surface and decrease downward heat fluxes obtained from the ABL scheme, resulting in further artificial cooling of the lowest layer above the surface. Thus, the fluxes obtained from the surface layer scheme need to be consistent with those from the ABL scheme.

Based on the governing equations of fluid mechanics and thermodynamics with closure constants determined from laboratory experiments in a neutral stratification, Mellor (1973) proposed a second-order turbulence closure model applicable to a “stratified” ABL. He applied it to a surface boundary layer under MOST and succeeded in reproducing φ_m and φ_h empirically determined from Kansas field experiment data by Businger et al. (1971). Mellor and Yamada (1974) classified this model into several levels according to the degree of anisotropy of turbulence and developed the so-called Mellor–Yamada scheme. To fix several weaknesses of their scheme, such as the underestimation of the depth of the mixed layer and the rapid dissipation of turbulence in the stable stratification, Nakanishi (2001) modified it by introducing more detailed parameterizations for pressure-strain and pressure-temperature-gradient covariances, devising the expression for the master length scale, and determining the closure constants from a large-eddy simulation database. Nakanishi and Niino (2006, 2009) incorporated it into a weather prediction model and termed it the Mellor–Yamada–Nakanishi–Niino (MYNN) scheme. Kitamura (2010) further modified the MYNN scheme by reducing a relaxation time scale of the pressure-temperature-gradient covariance for stable stratification (Canuto et al. 2008).

In this paper, we will first derive φ_m and φ_h from the MYNN scheme with the modification by Kitamura (2010) and compare them with stability functions proposed by previous studies. We will then approximate them by analytic functions of ζ, which are convenient for practical use such as computing momentum and heat fluxes. The fluxes obtained from the MYNN scheme will be also compared with those obtained from the stability functions in the previous studies, as well as observational data of the Surface Heat Budget of the Arctic Ocean experiment (SHEBA), which were collected on a pack ice that drifted approximately 2700 km in the Beaufort Gyre from October 1997 through September 1998 (Uttal et al. 2002; Andreas et al. 2007). Finally, we will discuss the bulk and flux Richardson numbers in the limit of ζ → ∞ that are predicted by the stability functions in the present and previous studies.

MOST defines φ_m and φ_h as functions of ζ by   

Substitution of Eqs. (1a) and (1b) into the level 2 version of the MYNN scheme gives simultaneous equations for φ_m and φ_h (Nakanishi 2001):   

Kitamura (2010) modified the MYNN scheme for stable stratification after Canuto et al. (2008). The modified simultaneous equations are obtained by replacing the closure constant A₂ by A₂/(1 + Ri) with the gradient Richardson number Ri, which is rewritten as   

We obtain φ_m and φ_h by solving Eqs. (3)–(8) using the Newton–Raphson method.

Several schemes of φ_m and φ_h have been proposed by previous studies and are used in weather and climate models. In the present study, we consider four schemes among them. In what follows, they will be described together with their validity range of ζ.

Högström (1988) reanalyzed the Kansas data and modified φ_m and φ_h of Businger et al. (1971). They increase linearly with increasing ζ:   

Equations (10a) and (10b) give that φ_m and φ_h are proportional to ζ and (2/3)^1/2 ζ^3/2, respectively, for large ζ.

Figure 1 shows dependences of the stability functions on ζ for the four previous schemes described in Subsection 2.2 and MYNN scheme described in Subsection 2.1. As ζ increases, differences among the stability functions become very large. Among the five schemes, CB05 gives the smallest values except for small ζ, while H88 gives the largest values. Compared with φ_h for the original MYNN scheme without the modification of Eq. (8) (Nakanishi 2001; not shown), the slope of the present φ_h is about 28 % larger.

Fig. 1.

Dependences of stability functions (a) φ_m and (b) φ_h on ζ for Högström (1988) (H88), Beljaars and Holtslag (1991) (BH91), Cheng and Brutsaert (2005) (CB05), Gryanik et al. (2020) (G20), and MYNN.

As will be described in the following section, φ_m and φ_h need to be integrated with respect to height. Since the integration of φ_m and φ_h computed from the MYNN scheme is difficult, we attempt to approximate them by analytic functions nearly similar to those of BH91:   

Because the third term on the right-hand side (RHS) of Eq. (13) is expected to be small for a large ζ, a_k and b_k are first determined by applying the first and second terms on the RHS for 1000 < ζ ≤ 2000. c_k, d_k (= 1–2), and e_k are then obtained by using all the terms on the RHS for 0 < ζ ≤ 10. It turned out that a contribution of the third term on the RHS is negligibly small, and φ_m and φ_h are well approximated as   

Statistics such as momentum and heat fluxes that are positive upward, −u²_* and −u_*θ_*, can be computed from Eqs. (1a), (1b), and (2). To avoid losing accuracy in gradient computation, we integrate Eqs. (1a) and (1b) from the roughness height to height z. The resulting relationships are written as   

By giving a value of Θ (z)–Θ (z_0h) = ΔΘ instead of that of U, one can also compute θ_* and u_* from Eqs. (15b) and (2), respectively. In this case, all the resulting θ_* and u_* from the five schemes vary monotonically with ζ.

To express the above statistics as a function of the bulk Richardson number Rb instead of ζ, we obtain the relationship between Rb and ζ. Substitution of Eqs. (15a), (15b), and (2) into the definition of Rb Becomes   

ψ_m and ψ_h of Eq. (16) for the five schemes (Eqs. 9–12, 14) become as follows:

a. Högström (1988) (H88)

b. Beljaars and Holtslag (1991) (BH91)

c. Cheng and Brutsaert (2005) (CB05)

d. Gryanik et al. (2020) (G20)

e. MYNN scheme (MYNN)

Figure 2a shows dependence of the stability parameter ζ on the bulk Richardson number Rb for the five schemes together with that of SHEBA data (Andreas et al. 2007), where z = 10 m and (z_0m, z_0h) = (0.00077, 0.00054) m. The values of the roughness lengths are adapted from Gryanik et al. (2020). In the range of Rb < ∼0.2, ζ of the SHEBA data tend to increase exponentially with increasing Rb. The dependence of ζ on Rb shows that the schemes are roughly classified into two groups: Group A that consists of H88, G20, and MYNN, and Group B that consists of BH91 and CB05. Group A shows comparatively good agreement with the SHEBA data, although its increase of ζ with Rb is slower than that of the SHEBA data. ζ of Group B, in contrast, increases nearly linearly with Rb. Interestingly, although the stability functions for MYNN and G20 are considerably different (Fig. 1), ζ for both schemes in the range of Rb < 0.2 are very similar.

Fig. 2.

(a) Stability parameter ζ, (b) friction velocity u_*/U, (c) upward heat flux −u_*θ_*/U, and (d) temperature scale θ_* as a function of Rb. Legends are the same as in Fig. 1 except for SHEBA data. Crosses show hourly averaged SHEBA data at 10 m height (Andreas et al. 2007) for the wind speed of 4 ± 1 m s⁻¹, where colors of the crosses are varied from purple to red with increasing wind speed from 3 m s⁻¹ to 5 m s⁻¹. The SHEBA data in (b) and (c) are normalized by the individual wind speed.

Figures 2b–d show dependences of u_*/U, –u_*θ_*/U, and θ_* on Rb for the five schemes together with those of the SHEBA data. In addition to the height and roughness lengths, we prescribe U = 4 m s⁻¹ and Θ₀ = 250 K that nearly correspond to medians of the SHEBA data. Note that the SHEBA data are plotted for U = 4 ± 1 m s⁻¹ and that the plots of −u_*θ_* and θ_* are not normalized by ΔΘ because the normalization makes difficult to distinguish differences among the schemes. Like u_*/U of the SHEBA data, both u_*/U of Groups A and B decrease monotonically with increasing Rb (Fig. 2b), although the decrease of u_*/U of Group B is slow.

−u_*θ_*/U of the SHEBA data tends to have a minimum near Rb = 0.06, but its scatter is extremely large there (Fig. 2c). −u_*θ_*/U of Group A also has a minimum between Rb = 0.06 and 0.08. In contrast, −u_*θ_*/U of Group B has no clear minimum and maintains large negative values in the range of Rb > 0.1. The minimum of −u_*θ_*/U occurs nearly at a maximum of θ_* (Fig. 2d). θ_* < 0.06 K of the SHEBA data between Rb = 0.2 and 0.26 is reasonably predicted only by MYNN. Not only θ_* but also u_*/U and −u_*θ_*/U for MYNN, however, vanish at Rb = 0.26 as will be discussed in Section 4.

Figure 2c demonstrated the very large scatter of −u_*θ_*/U of the SHEBA data especially near Rb = 0.06. We examine whether this large scatter occurred because the SHEBA data have been plotted over a wide range of the mean wind speed between 3 m s⁻¹ and 5 m s⁻¹ (4 ± 1 m s⁻¹).

Figure 3 shows dependences of u_*/U and −u_*θ_*/U on Rb for MYNN and G20 for U = 2, 3, and 5 m s⁻¹. The SHEBA data are also plotted for U = 2 ± 0.2, 3 ± 0.2, and 5 ± 0.2 m s⁻¹. u_*/U of the SHEBA data for U = 3 ± 0.2 m s⁻¹ and 5 ± 0.2 m s⁻¹ show large fluctuations (orange crosses and green circles in Figs. 3a, c), but both MYNN and G20 for U = 3 m s⁻¹ and 5 m s⁻¹ reproduce the decrease of u_*/U with increasing Rb reasonably well. Note that the lines for different wind speeds overlap. The scatter of u_*/U of the SHEBA data becomes further large for U = 2 ± 0.2 m s⁻¹ (blue triangles). MYNN cannot reproduce the non-zero values of u_*/U for Rb > 0.26, but G20 reproduces a part of them though the magnitude of u_*/U is small. We consider that u_*/U between 0.01 and 0.03 near Rb = 0.35 might be due to non-stationarity, surface heterogeneity, or both.

Fig. 3.

(a) Friction velocity u_*/U and (b) upward heat flux −u_*θ_*/U as a function of Rb for MYNN. (c) and (d) are the same as (a) and (b), respectively, except for G20. The friction velocities for different U for MYNN and G20 (a and c) are given by a single curve because of the normalization by U. Blue triangles, orange crosses, and green circles show hourly averaged SHEBA data at 10 m (Andreas et al. 2007) for the wind speed of 2 ± 0.2, 3 ± 0.2, and 5 ± 0.2 m s⁻¹, respectively, and are normalized by the individual wind speed.

−u_*θ_*/U of the SHEBA data for U = 5 ± 0.2 m s⁻¹ (green circles in Figs. 3b, d) is distributed between −0.003 K and 0 K in the range of Rb < ∼ 0.06, and it is reproduced by both MYNN and G20 reasonably well (green dashed-dotted lines).

In the case of U = 3 ± 0.2 m s⁻¹, the scatter of −u_*θ_*/U of the SHEBA data (orange crosses in Figs. 3b, d) is still considerably large as in Fig. 2c, and neither MYNN nor G20 (orange solid lines) can explain this large scatter. For U = 2 ± 0.2 m s⁻¹, conversely, −u_*θ_*/U of the SHEBA data (blue triangles) appears to be predicted well by both MYNN and G20 (blue dashed lines); however, the scatter of −u_*θ_*/U of the SHEBA data is actually large, which may reflect the large scatter of u_*/U (blue triangles in Figs. 3a, c). The large scatter of −u_*θ_*/U as in Fig. 2c does also occur for the small range of ± 0.2 m s⁻¹, especially for relatively weak wind speeds.

Gryanik et al. (2020) estimated that the most representative ranges of the roughness lengths for the SHEBA data are given by   

Figures 4a and 4c show dependence of −u_*θ_*/U on Rb for MYNN and G20 together with that of the SHEBA data, where α is varied from 0.01 to 1 for z_0m = 7.7 × 10⁻⁴ m. It turns out that the variations of α little affect a minimum of −u_*θ_*/U both for MYNN and G20, although Rb that gives the minimum slightly becomes smaller with increasing α.

Fig. 4.

Dependence of upward heat flux −u_*θ_*/U on Rb for MYNN (a and b), and for G20 (c and d). The lines in (a) and (c) show dependence of the heat fluxes on α, the ratio of roughness length for heat z_0h to that for momentum z_0m, for fixed z_0m = 0.00077, whereas those in (b) and (d) that on z_0m for fixed α = 0.01. Orange crosses and red circles show hourly averaged SHEBA data at 10 m (Andreas et al. 2007) that have the wind speed of 3 ± 0.2 m s⁻¹ and are normalized by the individual wind speed, where the red circles indicate the data during summer of the SHEBA.

Figures 4b and 4d are the same as Figs. 4a and 4c, except that here z_0m is varied from 7 × 10⁻⁶ to 0.03 m for α = 0.01. Both minima of −u_*θ_*/U for MYNN and G20 decrease as z_0m increases. If z_0m around the site would have varied between 7 × 10⁻⁶ and 0.03 m for some reason, −u_*θ_*/U for MYNN in Fig. 4b nearly covers the variations of the SHEBA data except around Rb = 0.26, where its magnitude is smaller than the observed values of the SHEBA data. −u_*θ_*/U for G20 in Fig. 4d is nearly similar to that of MYNN. Except for z_0m = 7 × 10⁻⁶ m, however, its magnitude is larger than the observed values of the SHEBA data around Rb = 0.26.

Andreas et al. (2010) classified the period of the SHEBA into two seasons: “winter” and “summer”. Winter was the period when the ice was covered with dry snow that was deep enough to be blown off and drifted, while summer was the period when the snow melted, and ponds and leads appeared. For the period of winter, they showed variations of measured z_0m with u_* estimated from a bulk flux algorithm (Fig. 1 of their paper): (i) for u_* > 0.15 m s⁻¹, z_0m that is averaged in u_* bins of 0.02 m s⁻¹ interval is independent of u_* and is nearly constant, but (ii) for u_* ≤ 0.15 m s⁻¹, z_0m fluctuates over a range from 10⁻⁷ m to 0.1 m and the averaged z_0m seems to depend on u_*.

Since u_* ≤ 0.15 m s⁻¹ for the case of U = 3 m s⁻¹ gives u_*/U ≤ 0.05, almost all the SHEBA data for U = 3 ± 0.2 m s⁻¹ correspond to (ii) (orange crosses in Figs. 3a, c). The range of the fluctuations of z_0m in (ii) is nearly the same as that in Figs. 4b and 4d. Both results of MYNN and G20 (all lines in Figs. 4b, d) suggest that the scatter of −u_*θ_*/U of the SHEBA data may be due to that of z_0m for u_* ≤ 0.15 m s⁻¹. In contrast, −u_*θ_*/U of the SHEBA data for the period of summer (red circles in Fig. 4) is distributed comparatively near the orange line with z_0m = 7.7 × 10⁻⁴ m, illustrating that z_0m for U = 3 ± 0.2 m s⁻¹ is considered to be comparatively constant during summer. Thus, it appears that the scatter of z_0m may have originated from dry snow during winter. When we interpret the results from the SHEBA data, it may be important to pay more attention to the estimate of z_0m especially during winter.

In addition, the majority of the SHEBA data for U = 5 ± 0.2 m s⁻¹ (green circles in Figs. 3a, c) are in the range of u_*/U > 0.03 and correspond to (i). The relatively small scatter of the SHEBA data (green circles in Fig. 3) is indeed consistent with the characteristics of z_0m in (i).

The critical values of the bulk and flux Richardson numbers, Rb and Rf, and the turbulent Prandtl number Pr are often estimated by assuming the applicability of the stability functions in the limit of ζ → ∞ (e.g., Wyngaard 2010). Because these critical values do not necessarily coincide with the true ones, the subscript ∞ will be used instead of c to indicate a critical value (e.g., Rb_∞).

Rb in the limit of ζ → ∞, Rb_∞, can be computed from Eqs. (18)–(23). As is inferred from Figs. 2–4, Rb_∞ for H88 and MYNN are 0.22 and 0.26, respectively. However, Rb_∞ for the other schemes are infinite.

Rf is expressed using Eqs. (1a), (1b), and (2) as   

Pr is expressed using Eqs. (1a) and (1b) as   

Many field observations, laboratory experiments, and numerical simulations show that turbulence exists even when Ri (≈ Rb) and Rf exceed the critical value of 0.20–0.25 (e.g., Grachev et al. 2013). In this regard, based on the SHEBA data, Grachev et al. (2013) demonstrated that, as Ri and Rf exceed the critical values, turbulence decays rapidly and the slope of energy spectra of velocity components becomes larger than the value of −5/3 for the Kolmogorov power law, arguing that the critical values of Ri and Rf do not correspond to a transition from turbulent to laminar flows, but to that from Kolmogorov to non-Kolmogorov turbulence, and coincide with upper limits of the applicability of MOST. Note that Kolmogorov hypothesis is applied to the turbulent dissipation in the MYNN scheme (Mellor 1973).

Using Eqs. (1a), (1b), and (25), the equation for turbulent kinetic energy (TKE) in a local-equilibrium state is written as   

Based on the SHEBA data, Grachev et al. (2007) showed that Pr decreases with increasing Rb and is smaller than unity in the very stable case. Based on field experiments conducted in central United States, Howell and Sun (1999) showed that Pr is scattered around unity as a function of ζ in the range of 0.01 < ζ < 10, but mentioned that its behavior for ζ > 1 remains uncertain. Pr_∞ for MYNN is 1.25 (Table 1). Kitamura (2010) showed that Pr = K_m/K_h for the MYNN scheme with his modification becomes infinite in the limit of “Ri → ∞”. It is noted, however, that Ri does not become infinite but approaches a finite critical value in the limit of ζ → ∞.

We pointed out that MYNN had a critical value Rb_∞ of 0.26 in the limit of ζ → ∞. Many of the SHEBA data are in the range of Rb £ 0.26 (Figs. 2–4) and are reproduced well by MYNN, but a few of the SHEBA data for Rb > 0.26 are not. In contrast, G20 has no critical value of Rb in the limit of ζ → ∞. Although it predicts unphysically large Rf_∞, Gryanik et al. (2020) explained that, for practical purposes, they were forced to use the stability functions in the supercritical (non-Kolmogorov) regime.

The stability functions for G20 are based on the SHEBA data in both Kolmogorov and non-Kolmogorov regimes, to which Eqs. (1a) and (1b) are equally applied. The distributions of φ_m and φ_h estimated from the SHEBA data, however, seem to have considerable fluctuations (Fig. 1a of Gryanik et al. 2020). These fluctuations are likely to result from the difficulty of evaluating dU/dz and dθ/dz in extremely stable conditions and the applicability limit of MOST as was stated by Grachev et al. (2013). In any case, the stability functions applicable to all ζ have complicated nonlinear forms. We therefore, consider that the stability functions or turbulent fluxes for non-Kolmogorov, intermittent, and non-stationary flows should be predicted separately. For example, a total turbulent flux at the surface, F, may be written as   

Recently, Łobocki and Porretta-Tomaszewska (2021) showed that gradient-based similarity functions obtained from the MYNN scheme are in good agreement with empirical functions proposed by Sorbjan (2017) for Ri < 0.2 in which the Kolmogorov regime prevails. They also stated that further progress requires a reliable parameterization of the turbulent dissipation in the non-Kolmogorov regime.

The stability functions for momentum and heat in the stable surface layer were derived from the MYNN scheme with the modification by Kitamura (2010) and were shown to be approximated by φ_m = 1 + 4.8ζ and φ_h = 0.74 + 6.0ζ, respectively, for practical applications. The fluxes computed from the functions for the MYNN scheme were nearly the same as those from the functions by Högström (1988) and Gryanik et al. (2020), and they show better agreement with the SHEBA data than those from the functions by Beljaars and Holtslag (1991) and Cheng and Brutsaert (2005).

Like the stability functions of Högström (1988), the approximated stability functions for the MYNN scheme have an advantage that the momentum and heat fluxes are calculated without iterations. They predict that the bulk and flux Richardson numbers, Rb and Rf, in the limit of ζ → ∞ are 0.26 and 0.21, respectively, which are about 20 % and 25 % larger than those of Högström (1988), respectively. Turbulence has been reported to exist for Ri (≈ Rb) and Rf exceeding the critical value of 0.20–0.25, although it becomes non-Kolmogorov, intermittent, and non-stationary (Grachev et al. 2013). In fact, a large part of the SHEBA data are in the range of Rb ≤ 0.26, but some of them are outside the range, i.e., Rb > 0.26. To realistically reproduce the data for Rb > 0.26, we plan to formulate the term expressing the effect of non-Kolmogorov, intermittent, and non-stationary turbulence.

It is recommended that the present stability functions, Eqs. (14a) and (14b), derived from the MYNN scheme are used to give the surface boundary conditions in weather and climate models in which the MYNN scheme is incorporated to represent the ABL, because they not only provide fluxes consistent with those predicted by the ABL parameterization, but also attain the good performance against the SHEBA data and the numerical efficiency without iterative computation.

Finally, the present study also suggests that the large scatter of heat fluxes for relatively weak wind speeds in the SHEBA data is likely to have been caused by significant variations of the surface roughness during winter. This needs to be confirmed by future field observations over snow-covered ice.

The authors thank two anonymous reviewers for their valuable suggestions to improve the manuscript. Data were provided by NCAR/EOL under the sponsorship of the National Science Foundation. https://data.eol.ucar.edu/. This work was partly supported by JSPS KAKENHI Grant Number 19K03967 and by MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Large Ensemble Atmospheric and Environmental Prediction for Disaster Prevention and Mitigation; Project ID:hp210166).

•  Andreas,  E. L.,  C. W.  Fairall,  P. S.  Guest, and  P. O. G.  Persson, 2007: Tower, interpolated hourly measurements to 2.5 and 10m at Met City (ASFG). EOL Dataset archive, UCAR/NCAR—Earth Observing Laboratory, doi:10.5065/D6862DVD.
•  Andreas,  E. L.,  P. O. G.  Persson,  R. E.  Jordan,  T. W.  Horst,  P. S.  Guest,  A. A.  Grachev, and  C. W.  Fairall, 2010: Parameterizing turbulent exchange over sea ice in winter. J. Hydrometeor., 11, 87-104.
•  Beljaars,  A. C. M., and  A. A. M.  Holtslag, 1991: Flux parameterization over land surfaces for atmospheric models. J. Appl. Meteor., 30, 327-341.
•  Businger,  J. A.,  J. C.  Wyngaard,  Y.  Izumi, and  E. F.  Bradley, 1971: Flux-profile relationships in the atmospheric surface layer. J. Atmos. Sci., 28, 181-189.
•  Canuto,  V. M.,  Y.  Cheng,  A. M.  Howard, and  I. N.  Esau, 2008: Stably stratified flows: A model with no Ri(cr). J. Atmos. Sci., 65, 2437-2447.
•  Cheng,  Y., and  W.  Brutsaert, 2005: Flux-profile relationships for wind speed and temperature in the stable atmospheric boundary layer. Bound.-Layer Meteor., 114, 519-538.
•  Dyer,  A. J., 1974: A review of flux-profile relationships. Bound.-Layer Meteor., 7, 363-372.
•  Grachev,  A. A.,  E. L.  Andreas,  C. W.  Fairall,  P. S.  Guest, and  P. O. G.  Persson, 2007: On the turbulent Prandtl number in the stable atmospheric boundary layer. Bound.-Layer Meteor., 125, 329-341.
•  Grachev,  A. A.,  E. L.  Andreas,  C. W.  Fairall,  P. S.  Guest, and  P. O. G.  Persson, 2013: The critical Richardson number and limits of applicability of local similarity theory in the stable boundary layer. Bound.-Layer Meteor., 147, 51-82.
•  Gryanik,  V. M.,  C.  Lüpkes,  A.  Grachev, and  D.  Sidorenko, 2020: New modified and extended stability functions for the stable boundary layer based on SHEBA and parametrizations of bulk transfer coefficients for climate models. J. Atmos. Sci., 77, 2687-2716.
•  Högström,  U., 1988: Non-dimensional wind and temperature profiles in the atmospheric surface layer: A reevaluation. Bound.-Layer Meteor., 42, 55-78.
•  Holtslag,  A. A. M.,  G.  Svensson,  P.  Baas,  S.  Basu,  B.  Beare,  A. C. M.  Beljaars,  F. C.  Bosveld,  J.  Cuxart,  J.  Lindvall,  G. J.  Steeneveld,  M.  Tjernström, and  B. J. H.  Van De Wiel, 2013: Stable atmospheric boundary layers and diurnal cycles: Challenges for weather and climate models. Bull. Amer. Meteor. Soc., 94, 1691-1706.
•  Howell,  J. F., and  J.  Sun, 1999: Surface-layer fluxes in stable conditions. Bound.-Layer Meteor., 90, 495-520.
•  Kitamura,  Y., 2010: Modifications to the Mellor–Yamada–Nakanishi–Niino (MYNN) model for the stable stratification case. J. Meteor. Soc. Japan, 88, 857-864.
•  Łobocki,  L., and  P.  Porretta-Tomaszewska, 2021: Prediction of gradient-based similarity functions from the Mellor–Yamada model. Quart. J. Roy. Meteor. Soc., early view, doi:10.1002/qj.4161.
•  Mellor,  G. L., 1973: Analytic prediction of the properties of stratified planetary surface layers. J. Atmos. Sci., 30, 1061-1069.
•  Mellor,  G. L., and  T.  Yamada, 1974: A hierarchy of turbulence closure models for planetary boundary layers. J. Atmos. Sci., 31, 1791-1806.
•  Nakanishi,  M., 2001: Improvement of the Mellor–Yamada turbulence closure model based on large-eddy simulation data. Bound.-Layer Meteor., 99, 349-378.
•  Nakanishi,  M., and  H.  Niino, 2006: An improved Mellor–Yamada Level 3 model: Its numerical stability and application to a regional prediction of advection fog. Bound.-Layer Meteor., 119, 397-407.
•  Nakanishi,  M., and  H.  Niino, 2009: Development of an improved turbulence closure model for the atmospheric boundary layer. J. Meteor. Soc. Japan, 87, 895-912.
•  Olson,  J. B.,  J. S.  Kenyon,  W. A.  Angevine,  J. M.  Brown,  M.  Pagowski, and  K.  Sušelj, 2019: A description of the MYNN-EDMF scheme and the coupling to other components in WRF-ARW. NOAA Technical Memorandum OAR GSD, 61, 42 pp.
•  Poulos,  G. S., and  S. P.  Burns, 2003: An evaluation of bulk Ri-based surface layer flux formulas for stable and very stable conditions with intermittent turbulence. J. Atmos. Sci., 60, 2523-2537.
•  Sorbjan,  Z., 2017: Assessment of gradient-based similarity functions in the stable boundary layer derived from a large-eddy simulation. Bound.-Layer Meteor., 163, 375-392.
•  Uttal,  T.,  J. A.  Curry,  M. G.  McPhee,  D. K.  Perovich,  R. E.  Moritz,  J. A.  Maslanik,  P. S.  Guest,  H. L.  Stern,  J. A.  Moore,  R.  Turenne,  A.  Heiberg,  M. C.  Serreze,  D. P.  Wylie,  O. G.  Persson,  C. A.  Paulson,  C.  Halle,  J. H.  Morison,  P. A.  Wheeler,  A.  Makshtas,  H.  Welch,  M. D.  Shupe,  J. M.  Intrieri,  K.  Stamnes,  R. W.  Lindsey,  R.  Pinkel,  W. S.  Pegau,  T. P.  Stanton, and  T. C.  Grenfeld, 2002: Surface heat budget of the Arctic ocean. Bull. Amer. Meteor. Soc., 83, 255-276.
•  Wyngaard,  J. C., 2010: Turbulence in the atmosphere. Cambridge University Press, New York, 393 pp.