Stability Functions in the Stable Surface Layer Derived from the Mellor–Yamada–Nakanishi–Niino (MYNN) Scheme

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.8 z / L and φ h = 0.74 + 6.0 z / 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 sta bility 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 corre spond to a transition from Kolmogorov to non-Kolmogorov turbulence.


Introduction
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 underesti-mation 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 andNiino (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.

Stability functions
2.1 Stability functions derived from the MYNN scheme MOST defines φ m and φ h as functions of ζ by where k is the von Kármán constant (= 0.4), U is the mean horizontal wind speed, Θ is the mean temperature, u * is the friction velocity, and θ * is the temperature scale. ζ = z/L is a dimensionless height normalized by the Obukhov length where g is the acceleration of gravity, and Θ 0 is the reference temperature. Substitution of Eqs. (1a) and (1b) into the level 2 version of the MYNN scheme gives simultaneous equations for φ m and φ h (Nakanishi 2001): where q * is a dimensionless turbulent velocity in a local-equilibrium state: l is a turbulent length scale: The values of C 2 and C 3 are based on Olson et al. (2019). Incidentally, Nakanishi (2001) maintained l = kz/3.7 for ζ ³ 1 in the diagnostic equation for the master length scale that consists of three length scales. This procedure was required to fit diagnosed length scales to those of large-eddy simulation data, because his equation gives a master length scale smaller than l in Eq. (5) when another length scale has a comparable magnitude. 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 2 by A 2 /(1 + Ri) with the gradient Richardson number Ri, which is rewritten as by substituting the relationship Ri = ζφ h /φ 2 m . We obtain φ m and φ h by solving Eqs. (3) -(8) using the Newton-Raphson method.

Stability functions proposed by previous studies
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 ζ.
a. Högström (1988) (hereafter H88) Högström (1988) reanalyzed the Kansas data and modified φ m and φ h of Businger et al. (1971). They increase linearly with increasing ζ: where the turbulent Prandtl number for neutral stratification, Pr 0 , is selected to be 0.95. Although similar formulae, φ m = φ h = 1 + 5ζ, proposed by Dyer (1974) are also used widely, they are not plotted in the following to avoid densely packed figures.
d. Gryanik et al. (2020) (hereafter G20) where Pr 0 is selected to be 0.98. Equations (12a) and (12b) give that φ m and φ h approach 11.2ζ 1/3 and 13.5, respectively, in the limit of ζ → ¥. 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.

Comparison of the stability functions
2.4 Approximate forms of φ m and φ h for the MYNN scheme 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: where the subscript k = m or h, φ m (0) = 1, and φ h (0) = Pr 0 (= 0.74; Businger et al. 1971;Mellor 1973).
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 These functional forms are similar to those of H88 and are consistent with MOST in the z-less stratification for extremely large ζ. These forms, however, are considerably different from those of the other schemes, which will be discussed in Section 4. Note that the curves for Eqs. (14a) and (14b) nearly overlap the red lines in Fig. 1.

Computation of statistics
Statistics such as momentum and heat fluxes that are positive upward, -u * 2 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 where z 0k (k = m or h) is the roughness length for momentum or heat, respectively, and ζ 0k = z 0k /L. If U and Θ 0 are prescribed, we can obtain u * from Eq. (15a) and then θ * by substituting it into Eq. (2): where Ψ m is the term inside the square brackets in the RHS of Eq. (15a). The heat flux divided by heat capacity is obtained by multiplying u * by θ * . 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 where Ψ h is the term inside the square brackets in the RHS of Eq. (15b).
c. Cheng and Brutsaert (2005) d. Gryanik et al. (2020) 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.

Comparisons of statistics among the five schemes
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 −1 and Θ 0 = 250 K that nearly correspond to medians of the SHEBA data. Note that the SHEBA data are plotted for U = 4 ± 1 m s −1 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 −1 and 5 m s −1 (4 ± 1 m s −1 ). Figure 3 shows dependences of u * /U and -u * θ * /U on Rb for MYNN and G20 for U = 2, 3, and 5 m s −1 .

Dependence on the mean wind speed
The SHEBA data are also plotted for U = 2 ± 0.2, 3 ± 0.2, and 5 ± 0.2 m s −1 . u * /U of the SHEBA data for U = 3 ± 0.2 m s −1 and 5 ± 0.2 m s −1 show large fluctuations (orange crosses and green circles in Figs. 3a, c), but both MYNN and G20 for U = 3 m s −1 and 5 m s −1 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 −1 (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.
-u * θ * /U of the SHEBA data for U = 5 ± 0.2 m s −1 (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 −1 , 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 −1 , 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 −1 , 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 7.1 ´ 10 −6 £ z 0m £ 5 ´ 10 −2 m, and (24a) , 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 −1 , respectively, and are normalized by the individual wind speed.

Dependence on the roughness length
We examined whether the large scatter of -u * θ * /U of the SHEBA data was caused by the variations of the roughness lengths. Considering the results of Subsection 3.4, we set U to 3 m s −1 and plot the SHEBA data for the mean wind speed of 3 ± 0.2 m s −1 . 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 −4 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 α. Figures 4b and 4d are the same as Figs. 4a and 4c, except that here z 0m is varied from 7 ´ 10 −6 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 −6 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 −6 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 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 −1 and are normalized by the individual wind speed, where the red circles indicate the data during summer of the SHEBA.
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 −1 , z 0m that is averaged in u * bins of 0.02 m s −1 interval is independent of u * and is nearly constant, but (ii) for u * £ 0.15 m s −1 , z 0m fluctuates over a range from 10 −7 m to 0.1 m and the averaged z 0m seems to depend on u * . Since u * £ 0.15 m s −1 for the case of U = 3 m s −1 gives u * /U £ 0.05, almost all the SHEBA data for U = 3 ± 0.2 m s −1 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 −1 . 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 −4 m, illustrating that z 0m for U = 3 ± 0.2 m s −1 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 −1 (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).

Discussion
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 ¥ ).
Pr is expressed using Eqs. (1a) and (1b) as where K m and K h are the turbulent diffusivity coefficients for momentum and heat, respectively. Pr in the limit of ζ → ¥, Pr ¥ , is computed by substituting Eqs.
(9) - (12) and (14) into Eq. (26). Table 1 summarizes Rb ¥ , Rf ¥ , and Pr ¥ obtained from the five schemes. 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 where ε is the dissipation rate of TKE. Equation (27) shows that Rf must be less than 1 in a local-equilibrium state. Rf ¥ £ 1 is satisfied by H88, BH91, and MYNN (Table 1). Rf ¥ < 1 implies that turbulence can persist in the limit of ζ → ¥; for MYNN, however, turbulence does not exist in the limit of ζ → ¥ because substitution of ε = q 3 /B 1 l into Eq. (27) gives Eq. (4), and l vanishes as ζ → ¥ (Eq. 5), which is compatible with the finite value of Rb ¥ beyond which the Kolmogorov turbulence is not sustained. However, Rf ¥ = ¥ for CB05 and G20 may give an unphysical value of ε = -¥ in the limit of ζ → ¥, since Rb for which turbulence vanishes is infinite (Rb ¥ = ¥). By substituting the upper limit of the validity range of ζ instead of ζ → ¥ into the stability functions, we obtain Rf = 0.71 and 1.94 for CB05 and G20, respectively. Rf for G20 is still larger than 1.
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  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 where F Kol is obtained from the stability functions, Eqs. (14a) and (14b), for Kolmogorov turbulence and F nonKol is the contribution of non-Kolmogorov, intermittent, and non-stationary turbulence. A similar idea was introduced by Poulos and Burns (2003), although they considered that F nonKol is a physically based stochastic component to be determined from measurements over a variety of surfaces around the globe.
We plan to obtain a dataset based on direct-numerical and large-eddy simulations, and to design a simplified formulation of F nonKol in the near future.
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.

Conclusions
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.