A Package of Momentum and Heat Transfer Coefficients for the Stable Surface Layer Extended by New Coefficients over Sea Ice

The bulk transfer coefficients of momentum, heat, and humidity belong to the main ingredients of numerical weather prediction and climate models. They are needed for the calculation of turbulent fluxes in the surface layer and often rely on the Monin–Obukhov similarity theory requiring universal stability functions. The problem of a derivation of transfer coefficients based on different stability functions has been considered by many researchers over the years but it remains to this day. In this work, dedicated to the memory of S.S. Zilitinkevich, we also address this task, and obtain transfer coefficients from three pairs of theoretically derived stability functions suggested by Zilitinkevich and co-authors for stable conditios. Additionally, we construct non-iterative parametrizations of these transfer coefficients based on earlier work. Results are compared with state-of-the-art coefficients for land, ocean, and sea ice. The combined parametrizations form a package in a universal framework relying on a semi-analytical solution of the Monin-Obukhov similarity theory equations. A comparison with data of the Surface Heat Budget of the Arctic Ocean campaign (SHEBA) over sea ice reveals large differences between the coefficients for land conditions and the measurements over sea ice. However, two schemes of Zilitinkevich and co-authors show, after slight modification, good agreement with SHEBA although they had not been especially developed for sea ice. One pair of the modified transfer coefficients is superior and is compatible to earlier SHEBA-based parametrizations. Finally, an algorithm for practical use of all transfer coefficients in climate models is given.


Preamble and Goals
This paper is written for a special issue of the journal 'Boundary Layer Meteorology' in memory of S. S. Zilitinkevich. Sergey Sergeevich Zilitinkevich was a theoretical meteorologist, oceanographer, and physicist, who always applied his theoretical knowledge to practical problems. At the start of the numerical modeling era in boundary-layer meteorology he became interested in parametrizations for numerical models. In Zilitinkevich et al. (2002) (hereafter ZPK02) he summarized his activity concerning research in numerical modelling as follows: 'As is evident from the foregoing, boundary-layer meteorologists and numerical modellers have rather different perspectives on the subject of flux-profile relationships in the surface layer. The former usually focus on the nature of turbulence and do not bother about such 'technical' aspects as the correction functions; the latter need nothing but these functions. However, modellers do not have time to carefully follow all the latest development in boundary-layer meteorology, not to mention reformulation of recent achievements in terms of the correction functions. As a result, direct empirical information on the correction functions, although highly desirable, is still presented in only a few papers, in particular by King and Connolley (1997), and Derbyshire (1999). The present paper attempts to contribute further to the turbulent-flux parametrization in general circulation models, accounting both for new developments in the theory and for new data representing essentially different geographical sites. ' We remind the reader that the flux-profile relationships are key components of the Monin-Obukhov similarity theory (MOST, Monin and Obukhov (1954), Foken (2006) connecting the mean wind speed and the potential temperature with turbulent fluxes of momentum and heat (defined by Eqs. 9a and 9b in this paper), and here the correction functions ( f m and f h ) are transfer coefficients normalized by their neutral values. We call them normalized transfer coefficients (NTCs) in the following. They also represent important components of bulk parametrizations of turbulent fluxes in numerical models (for the definition of f m and f h see Eq. 5).
In ZPK02 S. S. Zilitinkevich and colleagues pointed to the important role of static stability (Brunt-Väisälä frequency N ) in the free atmosphere adjacent to the boundary layer. They developed new flux-profile relationships and NTCs as a function of the Obukhov length and N . ZPK02 showed (see their Figs. 1 and 2) that results of the new parametrization differ from those currently used, e.g., in the Weather Research and Forecasting (WRF) model (Jiménez et al. 2012), ECHAM6 (Giorgetta et al. 2012), the regional atmospheric climate model HIRHAM5 (Christensen et al. 2007;Dorn et al. 2019), in coupled atmosphere-ocean models like ECHAM6-FESOM (Sidorenko et al. 2015) and in many others (Chang et al. 2020, and references therein). In these models, the parametrizations of Louis et al. (1982) (LTG82) are used, which read as: where c 1 = c 2 = 10, Ri b is the bulk Richardson number (defined by Eq. 21 in this paper). These values of c 1 and c 2 are updated constants (Viterbo et al. 1999), so differ from those originally suggested by LTG82. In this form the approach has been used over many years also in the numerical weather prediction model of the European Centre for Medium Range Weather Forecast (ECMWF) (Integrated Forecasting System IFS, Cy47R3, called ECMWF model in the following). But nowadays, another (iterative) approach is applied in the model, which is given in Sect. 4. ZPK02 indicated that results from models using the turbulent-flux parametrization (1) significantly overestimate the theoretical prediction as well as the results of measurements. Our research, presented here, aims to fill the still existing gap between current achievements in theory and parametrizations of f m and f h in numerical models, indicated by S. S. Zilitinkevich. Our focus is on the parametrizations of f m and f h using the theories of S. S. Zilitinkevich and colleagues suggested after publication of ZPK02 (Zilitinkevich and Esau 2007;Zilitinkevich et al. 2013) (in the following ZE07 and ZEKRE13, respectively). Following the suggestion of S. S. Zilitinkevich that all available data on surface layer fluxes, which 'represent essentially different geographical sites', must be taken into account, we consider not just one parametrization for one particular flux-profile relationship, but a package of relationships (based on six pairs of empirical stability correction functions and three new pairs of theoretical ones of ZE07 and ZEKRE13), which altogether cover the entire range of stability and roughness parameters as observed during the most famous and most comprehensive campaigns for atmospheric surface layer conditions over sea ice, land, and ocean. For all members of the package we establish a hierarchy of new non-iterative parametrizations of turbulent fluxes, which are represented by formulae of the same universal functional form, following the approach proposed by Gryanik and Lüpkes (2018) (GL18) and Gryanik et al. (2020) (GLGS20) on the basis of MOST.
The same strategy has been applied for stability functions by Gryanik et al. (2021) (GLSG21 in the following), but no functions of Zilitinkevich et al. had been included so far in the package. Also, similar to ZPK02, we compare now the NTCs against data for polar regions for all members of the package, both empirically and theoretically derived ones. The main reason for the choice of this region is that stable surface (boundary) layers are very common there. S. S. Zilitinkevich put forward a special term, 'long-lived boundary layer', which then was well accepted by the meteorological community. For the comparison we chose the data of the Surface Heat Energy Budget of the Arctic Ocean campaign (SHEBA) (Andreas et al. 1999;Uttal et al. 2002) (see also GAFGP07) since for conditions over sea ice it is the most comprehensive dataset and of all datasets it has the largest range in stability. Furthermore, the typical conditions over Arctic sea ice represent an ideal measurement environment. We stress that some MOST stability correction functions have been compared with SHEBA data previously (e.g. GAFGP07; Sorbjan and Grachev 2010). However, a comparison of NTCs for sea ice was presented only by GLGS20 for GAFGP07 and GLGS20 stability correction functions. To our knowledge, for the first time, we compare all members of the package (18 NTCs in total) with SHEBA. We underline that this is not equivalent to a validation because all other stability correction functions are originally based on datasets over land. It was already proven that for land conditions and in the stability ranges, for which they were defined, the stability correction functions show good agreement with the data used for their derivation. Thus our comparison aims to point to differences between the NTCs derived from datasets in different regions of the world and to show their differences to the polar data.
Once more repeating ZPK02 'the present paper attempts to contribute further to the turbulent-flux parametrization in GCMs, accounting both for new developments in the theory and for new data representing essentially different geographical sites'. Our research follows his clear methodology and further develops several of his ideas related to the treatment of the stable surface layer.

Introduction
The interaction between atmosphere, land, ocean, and sea ice is governed by radiation fluxes and turbulent fluxes of heat and momentum in the near-surface atmospheric and oceanic layers. In this work the focus is on the turbulent atmospheric fluxes, which are related to the mean wind and surface characteristics as aerodynamic roughness, thermal roughness, and surface temperature, and also to the atmospheric stability. The accuracy of surface layer flux parametrizations describing the impact of unresolved physical processes on resolved ones is still limited, especially when they are applied to weather prediction and climate models, due to their large horizontal grid sizes and low vertical resolution. The key difficulty is the parametrization of the turbulent fluxes in the stable surface layer.
The stable surface layer is often observed in polar regions over sea ice in all seasons as documented, e.g., by the SHEBA measurements, but it occurs also in mid-latitudes over land during night-time, as documented, e.g., by the Cooperative Atmospheric Surface Exchange Study (CASES-99) (Poulos et al. 2002). Furthermore, it is also often observed over ocean during warm air advection across colder water (Tjernström et al. 2015). Stable boundary layer turbulence is weak in comparison with neutral and convective turbulence. Moreover, it is intermittent, sporadic, patchy, and non-stationary (Poulos et al. 2002;Uttal et al. 2002;Grachev et al. 2007). The complicated structure of turbulence during stable stratification is related to anisotropy and large-scale spatial intermittency (e.g., pancake eddies) (Sorbjan and Balsley 2008), and with the non-local interaction with solitary and random internal gravity waves (Sun et al. 2004). Further complexity is related to the weakening of turbulence due to relaminarisation features (Sun et al. 2004), to the formation of nocturnal low-level jets, to effects of density currents (Cheng and Brutsaert 2005), and to isopycnical slope flows (Cheng and Brutsaert 2005). Interaction of the atmospheric flow with surface inhomogeneities affects the stratified turbulence as well (e.g., Lüpkes and Gryanik 2015;Bou-Zeid et al. 2020) and causes additional difficulties.
In numerical atmospheric models the turbulent transport of momentum and heat in the surface layer is usually described by MOST. Important ingredients of MOST are stability functions (SF in the following) depending on the single non-dimensional stability parameter ζ defined as: where z is the distance to the surface and L is the Obukhov length scale, which is positive for stable conditions, negative for unstable conditions and becomes infinite in the neutral case. In Eq. 2, τ and H are the momentum and heat fluxes, u 2 * = τ/ρ and θ * = −H /(ρ c p )u * denote the characteristic friction velocity and characteristic temperature scale, g is acceleration due to gravity; θ v is the virtual potential temperature at some reference level. Using virtual potential temperature assumes that the effect of humidity is considered in the simplest approximation (see, e.g. appendix in GL18). SFs depend on ζ , but also on other non-dimensional parameters (coefficients). In ideal environmental conditions, where the turbulence is statistically quasi-stationary, turbulent fluxes are independent of height above an underlying horizontally homogeneous surface, and tangential stress is aligned with the mean wind. Then, these coefficients become universal numerical constants. Thus, all complexity of stable surface layer turbulence described above is represented by the variety of the MOST stability correction functions ψ m for momentum and ψ h for heat. They are an integral form of the MOST SFs and are abbreviated by SCF in the following, see Eqs. 9a and 9b.
It is a challenge to parametrize the NTCs for all SCFs given in the literature in a universal framework that is well suited for practical use in numerical models. We consider here a package including currently used SCFs, which consists of the empirical SCFs suggested by Businger et al. (1971), Dyer (1974) (BD), Holtslag and De Bruin (1988) (HB88), Beljaars and Holtslag (1991) Cheng and Brutsaert (2005) Grachev et al. (2007) (GAFGP07), Gryanik et al. (2020) (GLGS20), as well as the theoretical SCFs proposed by ZE07 (ZE07-I and ZE07-II in the following) and ZEKRE13. The NTCs differ from each other by the rate of decay with increasing Ri b . Thus, e.g., the BD and HB88 empirical SCFs describe turbulence with finite critical Ri b , while all other SCFs have no critical Ri b . The theoretical SCFs of ZE07 and ZEKRE13 also do not have critical Ri b by construction.
GL18 suggested a semi-analytical method for the non-iterative parametrization of fluxes and derived bulk transfer coeffcients as functions of Ri b , and stability parameters m = z/z m for momentum and t = z/z t for heat, where z 0 is the surface aerodynamic roughness length scale and z t is the corresponding length scale for heat (see also GLGS20 and GLSG21). The most general parametrizations of NTCs f m and f h depend on Ri b and on m and t . Simplified versions of the parametrizations depend on Ri b only, similar to LTG82 (Eq. 1) and are valid for special SCFs only (e.g., Launiainen 1995;Blümel 2000;Li et al. 2014). We focus on the most general parametrizations; the simplifications can be easily derived, if required.
GL18 applied their approach to the SCFs introduced by GAFGP07, and GLGS20 considered the parametrization of their new SCFs. Using the same approach, the parametrizations for BD, HB88, BH91, CB05 were also derived in GLSG21 earlier. Our new goals are: (i) Most importantly, the derivation and analysis of new bulk normalized transfer coefficients for momentum f m (Ri b , m , t ) and heat f h (Ri b , m , t ) based on the theoretical SCFs of ZE07 and ZEKRE13, as well as the comparison with the previously derived ones based on the empirical SCFs of BD, HB88, BH91, CB05, GAFGP07, and GLGS20. (ii) The establishment of universal non-iterative parametrizations for these transfer coefficients and analysis of their advantages and shortcomings. (iii) The comparison of theoretical and all empirical normalized transfer coefficients mentioned above with SHEBA data. Furthermore, the validation of the corresponding non-iterative parametrizations by comparison with the corresponding iterative schemes. (iv) The modification of original SCFs of ZE07 and ZEKRE13 in order to improve the agreement of the related NTCs with SHEBA data. In this study, it is not our goal to tune free constants of other NTCs to improve their agreement with SHEBA measurements. (v) Giving practical recommendations to modellers on the basis of results from (i) to (iv).
This will result in an extension of the package of non-iterative parametrizations of momentum and heat fluxes presented by GLSG21, and will clarify which of the NTC parametrizations are applicable to polar regions and which are rather suggested for an application only to other regions.
The availability of the non-iterative parametrization as a package, in which for each package member transfer coefficients are represented by formulas of the same universal functional form, facilitates the application in NWP and climate models. For example, in the same model, one package member can be applied over sea ice, another one simultaneously over land surfaces, and another one over open ocean depending on the area for which the corresponding SCFs were originally optimized (e.g., GLSG21). Also, the package might be used in a traditional single-function mode, when an existing iterative scheme is replaced by a more efficient non-iterative one. The implementation is especially easy in models, which already use non-iterative parametrizations. When the non-iterative parametrizations are used in combination with iterative parametrizations, non-iterative parametrizations are useful for a proper choice of an initial guess (e.g., Grachev and Fairall 1997). First results applying only two members of the package to the regional Arctic climate model HIRHAM5 are promising (Schneider et al. 2022). They used the GL18 non-iterative parametrization for f m and f h (based on the GAFGP07 stability functions) over sea ice, but the default ones of LTG82 everywhere else. It is obvious that further similar studies, which are based on several members of the package (accounting for the differences in the stability conditions over sea ice, land, and ocean), need to be carried out.
The manuscript is organized as follows: in Sects. 3-5.1 we present the background, and in Sects. 5.2-6.4 the new results, which are summarized in Sect. 7. Finally, for completeness, the technical details and the algorithm for practical use are given in the Appendices.

Bulk Formulation of Turbulent Fluxes Following Monin-Obukhov Similarity Theory
Here, the presentation of the parametrization of turbulent fluxes on the basis of MOST follows GLGS20 (sf. Garratt 1994). The turbulent tangential surface stress τ (absolute value τ , often called vertical flux of momentum) is defined as: with the mean horizontal wind vector U(z) = (U (z), 0), and the mean wind speed U (z) is taken to be zero at the underlying surface. Here, height z is located above the roughness sublayer and below the boundary-layer top. C d (C d ≥ 0) is the bulk transfer coefficient (TC) for momentum, which is often called drag coefficient, and ρ is air density. Similarly, the equation for the heat flux H follows as: with the bulk TC for heat C h , (C h ≥ 0), θ v the virtual potential temperature. Using the virtual potential temperature assumes that the turbulent mixing is the same for moisture as for temperature, thus C h = C t = C q , where C t and C q are the TCs for temperature and for moisture, respectively. c p is specific heat at constant pressure. Index 0 refers to the surface value.
Both C d and C h can be written as: where C dn and C hn are the TCs for neutral stratification: and f m , f h are the normalized transport coefficients for momentum and heat. κ is the von Kármán constant, which is set here to 0.4, Pr 0 is the neutral-limit turbulent Prandtl number defined as the ratio of the momentum diffusion coefficient to the scalar diffusion coefficient in neutral conditions.
In the region of applicability of MOST the normalized transport coefficients f m and f h are obtained as: Here, f m (0) = f h (0) = 1 for neutral stability by definition. Height z has to be located above the roughness sublayer and below the boundary-layer top h (z 0 z h). It is assumed that fluxes are independent of height (with accuracy ∼ 10%), that the surface is homogeneous in horizontal direction, and that the turbulent flow is quasi-stationary. The MOST SCFs (sometimes called integrated SCFs) ψ m (ζ ) for momentum and ψ h (ζ ) for heat depend on the MOST stability parameter ζ = z/L. According to standard conventions, the Obukhov length L is positive in stably stratified layers (L > 0) and negative in convective layers (L < 0). In the stable surface layer L does not differ significantly from the local Obukhov length , which is defined by Eq. 2, where τ and H are the local fluxes. The SCFs ψ m (ζ ) and ψ h (ζ ) define the diabatic corrections to the mean wind speed and the virtual potential temperature θ v (z) profiles as: where the neutral contributions are represented by logarithmic terms depending on the roughness parameters m and t . The SCFs are negative for stable stratification (ψ m < 0, ψ h < 0) and positive for unstable conditions (ψ m > 0, ψ h > 0). The SCFs ψ k (ζ ) are related to SCFs φ k (ζ ) as: where I m = 1 and I h = Pr 0 . Here and in the following the SCFs ψ m (ζ ) and ψ h (ζ ) are normalized in agreement with the normalization φ m (0) = 1 and φ h (0) = Pr 0 of the corresponding SFs. Also, here and in the following a so-called multiplicative formulation of SCFs is used, where the factor Pr 0 is included in the definitions of ψ h , see, e.g., Eq. 11 below. For the relation of this formulation to the other ones, see Appendix A, also Appendix B of GLGS20. The TCs for neutral stratification C dn and C hn depend on the roughness parameters m and t only, as it should be. The SCFs ψ m (ζ ) and ψ h (ζ ) describe the corrections to the mean wind speed and mean temperature profiles due to effect of stability only. However, the normalized transport coefficients f m and f h depend on both the roughness and stability. The interplay between roughness and stability is described by the second terms in brackets in Eqs. 7 and 8, which depend on the parameters ζ , m , and t .
After specifying the similarity functions ψ m (ζ ) and ψ h (ζ ) the system of Eqs. 3-9b with known external variables U, θ v and z 0 , z t can be solved for the unknown variables τ and H .

A Package of SCFs
While many SCFs have been described in the literature, we follow the choice already motivated in Sect. 2. The SCFs included in our package for a detailed study are summarized in the following.
We consider six pairs of empirical SCFs. These are the SCFs of BD, HB88, BH91, CB05, GAFGP07, and GLGS20, and three more theoretically derived pairs of ZE07-I, ZE07-II, and ZEKRE13. The theoretical SCFs are presented here in a functional form, which is obtained by redefining notations for the independent variables and parameters of the original formulas (see Appendix A for details). For all functions, with exception of the empirical functions BD, GLGS20, and the theoretical ones of ZE07-I, ZE07-II, and ZEKRE13, the neutral-limit turbulent Prandtl number Pr 0 = 1.
1. BD SCFs: Businger et al. (1971) and Dyer (1974) SCFs are given by equations: where a m and a h are empirical constants. The most representative values of the empirical constants are a m = 5, a h = 5, and Pr 0 = 1. The functions are based on the landmark 1968 Kansas field experiment (Businger et al. 1971), but are also supported by recent measurements and large-eddy simulation (LES) studies. 2. HB88 SCFs: These SCFs (Holtslag and De Bruin 1988) assume the Reynolds analogy between momentum and heat transport: ψ m (ζ ) = ψ h (ζ ). They are given as: with the four constants a = 0.7, b = 0.75, c = 5, d = 0.35 for both ψ m (ζ ) and ψ h (ζ ). These SCFs were derived using the data from measurements over a grass-covered surface at Cabauw station. 3. BH91 SCFs: Beljaars and Holtslag (1991) suggested the same SCF (12) as HB88 for momentum, but a new function for heat, which is given as: (13) with the same four constants a = 1, b = 0.667, c = 5, and d = 0.35 for both ψ m (ζ ) and ψ h (ζ ). These SCFs were established using the Cabauw station data (as HB88). The approach (13) is used also in the ECMWF model (version Cy47r3). There, ψ m is parametrized by Eq. 12 but with constants as in Eq. 13. We abbreviate the version used by ECMWF as BH91/ECMWF in the following. 4. CB05 SCFs: The SCFs of Cheng and Brutsaert (2005) read as: Here a m = 6.1, b m = 2.5 for momentum (k = m), and a h = 5.3 and b h = 1.1 (k = h) for heat. These functions are based on the measurement campaign CASES-99 (Poulos et al. 2002). 5. GAFGP07 SCFs: The SCFs of Grachev et al. (2007) are given as: , and the optimal constants a m = 5, b m = 0.77, a h = 5, b h = 5, and c h = 3. These functions were derived from SHEBA measurements. 6. ZE07-I SCFs: The function φ m of Zilitinkevich and Esau (2007) depends on one parameter, and φ h on two parameters (see their Eqs. 11a and 11b, which we call the model ZE07-I in the following). Using these functions in Eq. 10 we derived the corresponding SCFs as: where Pr where Pr  SF of BD (Eq. 11), but the SF for heat φ h depends on four constants instead of one for BD (see their Eqs. 70 and 86). Using the SFs of ZEKRE13, we derived the corresponding SCFs ψ m and ψ h as: with Pr 0 = 0.8 and the coefficients a m = 4.0, a h = 4.5, b h = 1.13, c h = −0.0062, and d h = 3.55, which are based on the original constants of ZEKRE13 SFs, as described in Appendix A. 9. GLGS20 SCFs: Gryanik, Lüpkes, Grachev and Sidorenko (2020) extended and improved the functions of GAFGP07 (15a) and (15b) to: with the neutral-limit Prandtl number Pr 0 = 0.98 and coefficients a m = 5.0, a h = 5.0, b m = 0.3, b h = 0.4. These functions were constructed in order to fit the SHEBA data as closely as possible. For this reason their results are very similar to those of the SCFs of GAFGP07.
These functions can be separated into two groups. The BD and HB88 SFs, both describing collapsing turbulence in the range Ri b ≤ Ri b,cr , where Ri b,cr is the critical bulk Richardson number, are included in the first group. The second group contains the functions of BH91, CB05, GAFGP07, and GLGS20, which describe turbulence decaying with Ri b but still persisting at large values of Ri b . The latter group also includes the theoretically derived ZE07-I, ZE07-II, and ZEKRE13 SCFs.
The main idea behind the functions ZE07-I is to avoid a critical Richardson number Ri cr , which-as shown by many datasets-does not exist in reality. ZE07 show that when the BD linear function for φ m (ζ ) is accepted, then the function φ h (ζ ) should grow faster than φ m (ζ ). On the other hand, it should also be linear in ζ for small ζ . Thus ZE07-I postulated a second-order polynomial fit, which guaranteed the required feature avoiding Ri cr .
Another problem was addressed by ZE07 and leads to the SCFs ZE07-II. This is based on the finding that additional length scales need to be accounted for in the calculation of surface fluxes besides ζ . One of these additional scales characterizes the non-local effect of the static stability above the inversion and another one characterizes the effect of the Earth's rotation.
The application of such a concept and LES modelling delivered data for the construction of the ψ-functions of ZE07-II. They represent an excellent fit to their obtained model data.
ZEKRE13 derived SFs on the basis of an energy-and flux-budget (EFB) turbulence closure model. Here, they involved not only the turbulent kinetic energy budget equation as in traditional closures but also an equation for the turbulent potential energy. The latter includes a dependence on the Brunt-Väisälä frequency N , which was also present in the ZPK02 parametrization of the NTCs. The EFB closure equations capture both the range of strong turbulence at small Richardson numbers, Ri, and the range with weak turbulence at large Ri so that no critical Ri exists. The corresponding SCFs ψ m and ψ h given here (Eqs. 18a and 18b) are obtained as described in Appendix A.
For all three pairs of ZE07-I, ZE07-II, and ZEKRE13 SFs no limitation is given explicitly for their validity. The LES data, presented in the corresponding articles, are shown partly until ζ max = 125 and partly until ζ max = 250. However, we keep in mind that φ m of BD, which are to some extent fundamental for the theories, originally was established from measurements only for 0 ≤ ζ < 1, (see Businger et al. 1971;Dyer 1974). Using Eqs. 18a beyond the range of their validity is a hypothesis. This treatment is in line with GAFGP07, Grachev et al. (2013), Kouznetsov and Zilitinkevich (2010), and Casasanta et al. (2021). The region of applicability of φ h cannot be larger than that one for φ m because two functions are closely related by the MOST equations. Thus, we assume the same range of validity 0 ≤ ζ < 1 for φ h . Nevertheless, in Table 1 we include the values for ζ max = 125 and corresponding Ri b,max (in brackets), which rely on the LES data.
We stress that the validity of the stability functions and stability correction functions is not limited to a certain range of roughness because these functions are independent on the values of z 0 and z t . This is in contrast to the transfer coefficients (normalized and non-normalized ones), which both depend on z 0 and z t .
The other physical and mathematical properties of SFs and SCFs presented here are discussed in the original papers, and also by Andreas (2002), Sharan and Kumar (2010), Tastula et al. (2015), Srivastava and Sharan (2019), GLGS20, Casasanta et al. (2021), and references in the above-mentioned papers. Table 1 The package of universal non-iterative parametrizations for stable surface layer transfer coefficients. BD refers to the stability correction functions of Businger et al. (1971), Dyer (1974). HB88, BH91, CB05, GAFGP07 and GLGS20 refer to empirical functions Holtslag and De Bruin (1988), Beljaars and Holtslag (1991), Cheng and Brutsaert (2005), Grachev et al. (2007) and Gryanik et al. (2020), respectively. BH91/ECMWF refers to the parametrization used in the current version (Cy47r3) of the ECMWF model. ZE07-I, ZE07-II and ZEKRE13 refer to theoretical functions of Zilitinkevich and Esau (2007) and Zilitinkevich et al. (2013), respectively. Columns are: (i) the maximal values of the MOST stability parameter ζ max ; (ii) the corresponding Ri b,max , for which the stability functions were originally established; (iii) the constants γ and (iv) ζ a , which are used in Eq. 22a; (v) the equations defining the corresponding parametrization for the transfer coefficients. The last three lines refer to the modified SCFs of ZE07-I, ZE07-II and ZEKRE13 and are abbreviated as ZE07-I/GL, ZE07-II/GL and ZEKRE13/GL for shortening

Universal Approach
As we stated in Sect. 2, we apply the method of GL18 for the derivation of non-iterative parametrizations. In this section we summarize this approach following mainly GLGS20.

Solution of the Governing MOST Equation
The method of GL18 is based on an approximate solution of the governing MOST equation (sometimes known as equation for the Obukhov length), which reads: where the bulk Richardson number combines the external forcing parameters wind and potential temperature according to: As in Eq. 3 wind speed U at z 0 is taken to be zero. An approximate solution ζ = ζ(Ri b , m , t ) of Eq. 20 is obtained as: where andR is an equivalent bulk Richardson number combining the bulk Richardson number Ri b , the neutral-limit Prandtl number Pr 0 , and the roughness parameters m and t . The values of exponent γ and of the parameter ζ a are determined numerically using a least-square fit method. The first term of the solution (22a) approaches asymptotically the exact solution at small Ri b → 0, and the second term approximates the exact solution at large Ri b . Moreover, the solution (22a) coincides with the exact solution at ζ = ζ a . Overall, the parametrization (22a) describes the exact solution in the finite parameter range 0 ≤ Ri b ≤ Ri b,max (practically, Ri b,max ∼ 0.2 − 0.4 or even larger depending on the selected SCFs, see Table 1). It is worth noting that Eq. (22a) with (22c) is similar to Eqs. 42 of GLGS20 where their terms C and A are determined by their Eqs. 43, 46, and 48, where the factor Pr −1 0 is lost in (46) for ψ ha due to a misprint. This time, we do not use the assumptions 1 − Pr 0 1, . These assumptions were motivated earlier (see GL18, GLGS20, GLSG21) by applications to Arctic conditions over sea-ice-covered surfaces. However, the neglected terms ψ m (ζ a / m ) and ψ h (ζ a / t ) can be essential over land surfaces, as well as for proper matching with the roughness sublayer. If so, these terms are not neglected in Eq. 20. Also, in general the neutral-limit Prandtl number Pr 0 is not always close to one (see, e.g., ZE07 and ZEKRE13). We also stress that although Pr 0 is occurring in the coefficients of the right-hand side of Eq. 22a, the coefficients do not really depend on Pr 0 because its cancellation with that one occurring in the definition of SCFs. However, the dependence on Pr 0 via the equivalent bulk Richardson numberRi b forms an essential property of the parametrization. It is obvious that the explicit inclusion of Pr 0 increases the generality of the parametrization (22a) and increases thus the region of applicability of the package considerably.
The flux-profile relationship given by Eq. 22a is a general result, which, in principle, is valid for an arbitrary pair of SCFs ψ m and ψ h . It is also obvious that the parametrizations described by this equation satisfy the quality criteria 1) to 6), see GLSG21, by construction. Note, however, that both values of γ and ζ a are sensitive to the particular functional form of SCFs ψ m (ζ ) and ψ h (ζ ) and to the ranges chosen for the optimisation of the parametrization.
When the solution ζ = ζ(Ri b , m , t ) (Eq. 22a) is known, one can use it in Eqs. 7 and 8 to find the transfer coefficients (5), and so the fluxes of momentum (3) and heat (4). An algorithm for practical use is described in Appendix B.
This parametrization can be used in models of different levels of complexity because for each pair of SCFs Eq. 22a generates a hierarchy of parametrizations. The most general parametrization depends on Ri b , m and t . The parametrizations of intermediate complexity follow from the general formulation when the roughness parameters for momentum and heat are assumed equal (z 0 = z t ), so t = m . This assumption, although not supported by observations, is often used in numerical models. At the lowest level of complexity the parametrizations can be optimized to representative values m and t of roughness parameters. These depend only on Ri b and result in only Ri b -dependent f m (Ri b ) and f h (Ri b ) functions, which are analogous to the parametrization of LTG82 (1) and can be used in numerical models similarly.

Parametrized Normalized Transfer Coefficients
To arrive finally at parametrized NTCs for the whole package, including the theoretical SCFs of ZE07-I, ZE07-II, and ZEKRE13 SCFs, two constants γ and ζ a in Eq. 22a must be calculated for each pair of NTCs. A least-square fit to obtain the best coefficients γ and ζ a has been applied already by GLSG21 to all SCFs of the package except to the SCFs of ZE07-I, ZE07-II and ZEKRE13, which we do here. For optimization we use the same ranges of Ri b , m and t for all six SCFs. The ranges are: 0 ≤ Ri b < 0.4 for stability, and 1.5 × 10 3 ≤ m ≤ 3 × 10 5 and m ≤ t ≤ 10 2 m for the roughness parameters. Based on z = 10 m this corresponds to 3.3 × 10 −5 m ≤ z 0 ≤ 6.7 × 10 −3 m, 0.01z 0 ≤ z t ≤ z 0 . These ranges of Ri b , m , and t are typical for polar sea ice regions, and are documented by SHEBA (see, e.g., GAFGP07) and other campaigns (see, e.g., GL18, GLGS20). It is important to note that these optimal ranges are far away from the roughness values shown in Figs. 1, 2, 3, 4, 5, 6 by the largest and smallest values of m . However, we found that the fit to sea-ice conditions in the given large range results in a sufficient agreement over land as well, while, vice versa, the best fit to land conditions did not result in sufficient accuracy over smooth sea ice.
The solution is unique for each pair of SCFs in the considered stability range. Plugging the constants in Eq. 22a and the result in Eqs. 7 and 8, where ψ m (ζ ) and ψ h (ζ ) are given by Eqs. 16 for ZE07-I, Eq. 17 for ZE07-II and Eqs. 18a, 18b for ZEKRE13 SCFs, we obtain the required parametrizations for the NTCs f m (Ri b , m , t ) and f h (Ri b , m , t ). The parametrizations are summarized in Table 1 for the whole package and visualized in Figs. 1, 2, 3, 4, 5, 6 (see dashed coloured lines) as a function of Ri b for the same values of roughness parameters m and t for all NTCs.
This and the other results are discussed in the next section.

Results and Discussion
In this section we derive the NTCs f m and f h depending on Ri b , m , and t for the theoretical stability correction functions (SCFs) of ZE07-I, ZE07-II, and ZEKRE13, and compare them with NTCs based on the empirical SCFs of BD, HB88, BH91, CB05, GAFGP07, and GLGS20, which were derived previously (Sect. 6.1). Also, we compare these exact NTCs with the corresponding non-iterative parametrizations and discuss their accuracy (Sect. 6.2) relative to the iterative schemes. Then we compare all coefficients with the SHEBA measurements (Sect. 6.3). The comparison of NTCs with SHEBA measurements had been shown until now only for the GAFGP07 and GLGS20 functions by GLGS20. Such a comparison is important. Namely, it shows which of the functions and in which stability range they can be applied to sea-ice-covered regions, maybe also beyond the range of ζ , for which they were originally derived (see Table 1). We remind the reader that the extension of the validity to large values of Ri b is a prerequisite for using NTCs in weather prediction and climate models for the determination of near-surface fluxes. At present, there is no practical alternative for such an extension. Finally, in Sect. 6.4 we describe an attempt to modify the stability func- tions (SF) of ZE07-I, ZE07-II, and ZEKRE13 to improve the agreement of the corresponding SCFs and NTCs with SHEBA data. The results are shown for all SCFs in Figs. 1,2,3,4,5,6,7,8,9,10. They contain the exact solutions (solid coloured lines), our non-iterative parametrizations (dashed coloured lines) and parametrizations of LTG82, which are often used in models (solid black lines). All NTCs are presented for three different pairs of roughness parameters and t and as a function of Ri b (0 ≤ Ri b < 0.4). We consider a very smooth ocean surface, a typical sea-ice surface, and a rough land surface. The sea-ice roughness parameters are m = 3 × 10 4 and t = m /0.7 (based on z = 10 m, z 0 = 3.3 × 10 −4 m, z t = 0.7z 0 , see, e.g., GLGS20). The extremely smooth ocean surface is represented by m = 10 7 , z 0 = 10 −6 m (see, e.g., Elvidge et al. 2016) and for the rough land surface we use m = 30, z 0 = 0.33 m (see, e.g., Garratt 1982).

New Normalized Transfer Coefficients and Comparison with Earlier Ones
The normalized transfer coefficients f m and f h based on ZE07-I, ZE07-II, and ZEKRE13 SCFs are obtained numerically by an iterative solution ζ = ζ(Ri b , m , t ) (which we call iterative or exact solution interchangeably) of the governing MOST Eq. 20 and substituting Several important properties of these NTCs can be understood by analysing the asymptotes of Eqs. 7, 8 with Eqs. 20 and 16-18a. The analysis reveals that at Ri b → 0 we have the common linear decrease of both f m and f h for ZE07-I and ZEKRE13 as: This linear dependence on Ri b can be explained by the fact that the corresponding SCFs follow the Businger-Dyer limit (11) at ζ → 0, and that ζ ∼ Ri b (see Eq. 22a) in this limit. The other NTCs based on empirical SCFs (11)-(15b) and (19a), (19b) have similar asymptotes, but the constants a m and a h are different for different functions. In contrast, ZE07-II SCFs (16) do not have this linear limit, and the corresponding NTCs decrease as: Also, we found that in the opposite limit of Ri b → ∞ all six functions f m and f h decrease following a power-law. The analytical asymptotes read as: with the exponents Δ k as: ZE07 − II : According to Eqs. 28 with (29) and (30) the NTCs f m for ZE07-I and ZEKRE13 SCFs represent long-tail functions, while the f h are short-tail functions. We remind the reader that if a function is decaying as 1/Ri 3 b at Ri b → ∞, or faster, it is called a short-tail function, but it is called a long-tail function if the decay is slower. Actually, the choice of the value Δ = 3 as transitional value is rather arbitrary, but we use it in our study as a reasonable working hypothesis, sf. LTG82, GLSG21. Thus, both f m and f h based on the ZE07-II SCFs belong to the same group of short-tail functions. It is obvious that the turbulent fluxes at Ri b → ∞ are larger for long-tail NTCs then for short-tail ones.
By comparing the theoretical NTCs of ZE07-I and ZEKRE13 (Figs. 5, 6 and in most concise form for sea-ice conditions only in Figs. 7 and 8) with the empirical ones of GLGS20 (Figs. 3, 4), it becomes obvious that the qualitative behaviour of the ZE07-I function f m is very similar to the corresponding GLGS20 function. This concerns the dependence on Ri b and on m . The reason for this similarity can be explained by the fact that the asymptote of f m for ZE07-I coincides with that one of GLGS20 in both limits of small and large Ri b . In the limits of small Ri b the coefficients a m in Eq. 26a are equal (a m = 5 for both). Also, the Δ m values are equal as well, because Eq. 28 with (29) and (30) is identical to Eqs. (60) in GLGS20, which reads f m ∼ 1/Ri 2 b . As a result, the NTCs are similar in both limits of small and large Ri b . This is sufficient for the similarity of f m in the whole range The quantitative agreement of ZEKRE13 and GLGS20 is also good, but slightly less pronounced due to the difference in the values of the coefficients: a m = 4 and Pr 0 = 0.80 for ZEKRE12, but a m = 5 and Pr 0 = 0.98 for GLGS20. The agreement can be improved by equalizing the values of these constants. We come back to this idea of a modification of the ZE07-I and ZEKRE13 NTCs in Sect. 6.4.
Finally, we stress that the described similarity of ZE07-I and ZEKRE13 f m with the GLGS20 f m is an unexpected result keeping in mind the difference of the asymptotes of the corresponding SFs at large ζ : ζ vs ζ 1/3 for ψ m (ζ ) and ζ 2 vs const for ψ h (ζ ). This difference in the SFs is reflected in a larger difference of NTCs for heat f h . Nevertheless, the agreement for f h is also reasonably good, as shown in the figures.
NTCs of ZE07-I and ZE07-II are closer to the BD curves than the other NTCs, except GLGS20. The figures also contain results of the LTG82 NTCs. As already discussed by GLSG21, the BH91 curves are those with the smallest differences to the LTG82 curves. In contrast, the difference between all theoretical NTCs and LTG82 is the largest.
Concerning the dependence of NTCs on roughness, we highlight the following. In contrast to the NTCs of LTG82, all coefficients presented in Figs. 1, 2, 3, 4, 5, 6 show a strong dependency on the surface roughness parameters m and t , except the ones based on BD ( Figs. 1 and 2, upper rows). The dependence becomes stronger with increasing Ri b , which is visible by the increasing spread of curves towards large Ri b . However, we see also that this growing dependency starts at different values for each function. For example, for BH91 the curves for f m start to diverge at Ri b ≈ 0.08 while for GAFGP the corresponding value is Ri b ≈ 0.1 and for GLSG20 it is Ri b ≈ 0.12. For the theoretical NTCs the dependence on m is the largest for ZEKRE13 and the smallest for ZE07-II. Figures 7 and 8 also contain the solutions of the BH91/ECMWF scheme. The corresponding NTCs differ only slightly from those of the HB88 and BH91 schemes. This can be expected because the same equations are used as in BH91 but just with other constants for ψ m . Thus, conclusions for the schemes HB88 and BH91 are valid also for the BH91/ECMWF scheme. The corresponding results of the non-iterative scheme for BH91/ECMWF are not shown but required values of constants are included in Table 1. The quality of agreement between the iterative and non-iterative scheme is the same as for BH91.

Comparison of Iterative and Non-iterative Solutions
The comparison of the iterative and non-iterative solutions for BD, HB88, BH91, CB05, GAFGP07, and GLGS20 has been shown already in earlier work (GL18, GLGS20, and GLSG21, see Figs. 1, 2, 3, 4). The main findings are that the non-iterative schemes are in most cases very close to the iterative solutions. The non-iterative solutions reproduce both long-tail functions and short-tail functions with the same quality of agreement, but there are some exceptions. Considering the panels with linear vertical axes, it seems that there is not much difference between the iterative and non-iterative solutions in the range 0 ≤ Ri b < 0.1, but the panels with logarithmic axes, focusing on the very stable stratification (large Ri b ), illustrate some differences (Figs. 12, 3, 4, 5, 6). The largest difference occurs for the very rough case (green lines). For the GLGS20, ZE07-I, and ZEKRE13 functions the non-iterative schemes underestimate both f m and f h , especially in the range 0.05 < Ri b < 0.15. But one must keep in mind the large scatter of measurements (see Sect. 6.3) and that the TCs C d = C dn f m and C h = C hn f h are goverened not only by f m and f h but also by C dn and C hn (thus by z 0 , z t ) with usually high uncertainty. The latter has a larger effect on TCs than the uncertainty of the non-iterative scheme for the determination of the NTCs.
Overall, the strong dependency of NTCs on the surface roughness parameters mentioned above is reproduced by the non-iterative solutions. And similar to the iterative solutions, in the non-iterative solution the dependency is also stronger at larger Ri b . Additionally, as for the iterative schemes, this growing dependency on Ri b starts at different values for each function.
Finally, all parametrizations approximate the exact NTCs most accurately for values of the roughness parameters for momentum m (blue curves in the figures) and for heat t , which are typical for sea-ice conditions. One can expect this result because the ranges, which were used for an optimization of the parameters ζ a and γ (Eq. 22a), were based on these typical values.
Overall, the non-iterative parametrizations reproduce the iteratively determined NTCs reasonably well. The largest biases exist for large surface roughness. Differences between the exact NTCs and the parametrizations are comparable with the scatter among the results for different functions. Thus, we consider the differences between results of the iterative and non-iterative schemes as being small (see the discussions of this issue in GL18, GLGS20, and GLSG21). They are also small as compared with the scatter in the measurements (see Sect. 6.3).
An analysis of the figures also reveals an interesting fact that several parametrizations, such as HB88, BH91, CB05, GAFGP07, GLGS20, ZE07-I (see Figs. 1, 2, 3, 4, 5, 6), have systematic biases (especially for smooth surfaces, see red curves) resulting in an underestimation of the exact NTCs at large Ri b . This finding has important consequences, which we discuss in the next section.

Comparison of NTCs with SHEBA Data
As we already mentioned in Sects. 1 and 2 SHEBA provides an extraordinary dataset because it covers a large range of stability (0 ≤ ζ < 100). This might be related to the special conditions over sea ice with often surface-based inversions during winter. The data represent measurements over about one year with frequent stable conditions. They were obtained from sonic anemometers installed at a 20 m tower in five levels (Andreas et al. 1999;Uttal et al. 2002) (GAFGP07). Due to the high quality of the data we use them to investigate the applicability of NTCs of all members of the package of stability correction functions over sea ice (see Table 1). As mentioned already, a disagreement with SHEBA data of the NTCs, derived from data in mid latitudes, means only that they have drawbacks for sea-ice-covered regions, but they should be part of the parametrization package covering more than just polar regions.
NTCs f m for momentum and f h for heat based on bin-averaged data are shown in Figs. 1, 2, 3, 4, 5, 6 (squares) as a function of the bulk Richardson number Ri b using linear-linear scaling (left column) and log-linear one (right column). Note that only the blue lines can be compared with the SHEBA data since the prescribed roughnesses for the green and red curves are outside of the SHEBA roughness range. The colours of the symbols represent SHEBA measurements in different heights but not for different roughnesses.
For small Ri b (see the range 0 ≤ Ri b < 0.08 in the left columns of the figures) all NTCs and corresponding parametrizations represent the data well. However, the analysis of the figures reveals that for large Ri b (practically Ri b > 0.08) not all NTCs result in a good agreement with measurements.
The best agreement with the SHEBA measurements among the empirical NTCs is obtained for the GLGS20 ones, which is slightly better than that obtained by the GAFGP07 NTCs. This improvement is not surprising because the functional form and parameters of the GLGS20 and GAFGP07 SCFs had been optimized using SHEBA measurements. As concerns the results of the theoretical NTCs, f m based on ZE07-I shows even a slightly better agreement with the SHEBA data in the range Ri b > 0.25 than the GLGS20 results. ZEKRE13 results for f m overestimate the measured data for Ri b > 0.1. For f h , GLGS20 shows the best agreement with SHEBA but also ZEKRE13 results agree well with the observations while ZE07-I results underestimate f h for Ri b > 0.2. ZE07-II results for both f m and f h differ much more from the observations for both weak and strong stability, which becomes especially obvious from the figures with logarithmic axes (middle rows and right columns in Figs. 5 and 6). This can be expected because the corresponding SCFs do not fulfill the BD limit at Ri b → 0. The functions with the largest discrepancy are the ones by LTG82, BD, CB05, and ZE07-II. We stress, however, that the discrepancy is relatively large only for the range Ri b > 0.08. But it shows that over sea ice the application of NTCs to the full range 0 ≤ Ri b < 0.5 is critical for the functions BD, HB88, BH91, and ZE07-II because apparently they do not agree well with the SHEBA data.
Caution is necessary when the green curves are considered. They look for some NTCs ( f m and f h of CB05 in Figs. 3 and 4, and f h of GAFGP07 in Fig. 4) as if they were good representations of the SHEBA data. But this is not the case because the roughness for the green curves ( m = 30, z 0 = 0.33 m) is much higher than any roughness observed during SHEBA.
Finally, the comparison of the non-iterative parametrizations with measurements points to an unexpected benefit, at least when the parametrizations are applied over sea ice. As we mentioned in the previous section, several non-iterative parametrizations (,e.g., HB88, BH91, CB05, GLGS20, and ZE07-I) underestimate the iteratively determined NTCs (see, e.g., Figs. 1,2,3,4,5,6), but this underestimation compensates an overestimation of the SHEBA data by actual NTCs for rough surfaces, and, thus, results in better agreement of the non-iterative parametrizations with the data. In this respect, for those cases the non-iterative parametrizations of NTCs are even better suited for their practical use in models over sea ice than their iterative counterpart. Therefore, one can consider the parametrizations as a modification of the iterative NTCs representing solutions of the MOST equations, which are better adjusted to SHEBA data.
It is interesting that the NTCs of ZE07-I and ZEKRE13 show a relatively good agreement with the SHEBA data, although the derivation of the corresponding stability functions SF and stability correction functions SCF is based on theory and LES data, forming a completely independent data source. The question arises again, why other functions, and thus the data they are based on, differ from the SHEBA data. It has been shown already by GLSG20 that the so-called Ranchi data (Srivastava et al. 2020) agree well with the SHEBA data, at least in the stability range where they show the highest accuracy. These data are valid for land conditions far from polar regions. It underlines that the question, why datasets show differences to each other, is not yet solved. The reason cannot be the surface roughness because the possible range of modification by roughness is much smaller than the difference between the NTCs of different authors (see Figs. 1,2,3,4,5,6).
In the next section we show that by slight modifications of the ZE07-I, ZE07-II, and ZEKRE13 SCFs the agreement of the related NTCs with the SHEBA measurements can be improved.

Modifications of ZE07-I and ZEKRE13 by Adjustment to SHEBA Data
A modification of the MOST SFs is common practice. It can be motivated by measurements, which become newly available, or by new theoretical physical ideas. The modification can concern the functional form of the stability functions (e.g., GLGS20) and the values of the empirical constants (e.g., Dyer 1974;Srivastava et al. 2020). Dyer (1974) modified the constants of the SFs proposed by BD, and Srivastava et al. (2020) changed the constants of GAFGP07 using the Ranchi data for land surfaces. GLGS20 modified the functional form of the GAFGP07 SFs using the same SHEBA data, but improving the description of the ζ − Ri b relationship. In this study, in line with Dyer (1974) and with Srivastava et al. (2020), but in contrast to GLGS20, we attempt to improve the agreement of the ZE07-I, ZE07-II, and ZEKRE13 NTCs with the SHEBA data keeping the functional form of SCFs unchanged and optimizing their constants.
By the trial and error method we found that an improved agreement to SHEBA data can be achieved by a Here and in the following we use the abbreviation ZE07-I/GL, ZE07-II/GL, and ZEKRE13/ GL for modified SCFs and NTCs for shortness. Results are presented in Fig. 9 for f m and Fig. 10 for f h . Apparently, these modifications improve the agreement of the NTCs obtained by the schemes ZE07-I/GL, ZE07-II/GL, and ZEKRE13/GL with the observations. For ZE07-II/GL, there is still a large difference for Ri b > 0.25 but an improvement is visible for 0.1 < Ri b < 0.2. We obtain now for ZE07-I/GL and ZEKRE13/GL a similar quality of agreement as for GLGS20 for both f m and f h . The best agreement for f m with respect to sea-ice conditions is found for ZE07-I/GL showing almost an optimal agreement with the measurements. It is even slightly better than for GLGS20 in case of f m while it is slightly worse than GLSG20 with respect to f h .
Actually, the close agreement with the results of measurements can be expected because in Sect. 6.1 the similarity of the asymptotes belonging to the ZE07-I and ZEKRE13 NTCs with those belonging to the GLGS20 NTCs was shown in both limits of small Ri b → 0 and of large Ri b → ∞ (see Eqs. 26a, 26b and Eq. 28 with 29). Since the SFs of GLGS20 were derived to approximate the SHEBA data as closely as possible, it is not surprising that ZE07-I/GL and ZEKRE13/GL with constants (31) and (33) also approximate the SHEBA data well.
Another reason for the good agreement between the parametrizations ZE07-I/GL and ZEKRE13/GL is also that the logarithmic term in the ψ h function of ZEKRE13 (Eq. 18b) is small against both other terms. Thus we have almost the same functional dependency on ζ or Ri b respectively for ZE07-I/GL and ZEKRE13/GL. However, it is also important to note that the neutral-limit Prandtl numbers Pr 0 of both schemes differ strongly after the modification. For ZE07-I/GL it coincides with that of the empirical GLGS20 SCF ψ h (Pr 0 = 0.98), but for ZEKRE13/GL with that of the theoretical Sukoriansky (2008) SCF (Pr 0 = 0.7). The latter was derived rigorously under several well established assumptions on the basis of spectral quasi-normal-scale-elimination (QNSE) theory for stably stratified turbulence (Sukoriansky and Galperin 2013, and references therein).
The modified values of the constants (31)-(33) lead to the modification of the fitting parameters γ and ζ a , which specify the parametrization. The optimal parameters are as follows: ZE07 − II/GL : γ = 4.30, ζ a = 14.0, The dashed coloured curves in Figs. 9 and 10 show the parametrized NTCs f m for momentum and f h for heat. The accuracy of the parametrization is similar to that for the NTCs based on the original SCFs.
Although the obtained results are impressive, the trial and error method for finding optimal coefficients cannot really replace a systematic sensitivity study leading to the 'best possible' modification of the SCFs. Keeping in mind that four or five parameters must be fitted (see, e.g., eqs. 18a and 18b), and that the coefficients are mutually dependent ones (see Appendix A), a more rigorous method is required for such a study. Moreover, the modification of the functional form might be necessary also, e.g., for achieving a consistency of ZE07-II SCFs with the Businger-Dyer limit (11) at ζ → 0. We leave this finding of the 'best' modification to future research.
Summarizing, at least for applications to sea-ice-covered regions we can recommend the choice of parameters (31), (32), and (33) more than the original ones. For this reason we include these new modified SCFs in Table 1. According to the table, the region of applicability is enlarged as compared to the original theoretical one 0 ≤ ζ < 1 because the new NTCs are derived on the basis of SHEBA data, which are available in the larger range 0 ≤ ζ < 100. Future studies will show to which extent our modification of the constants in the SCFs is compatible with the allowed possible range of constants following from the EFB turbulence closure model. If an appropriate choice is impossible, the modified functions must be considered as new SCFs, which are based on SHEBA data. These functions are complementary to the GAFGP07 and GLGS20 SFs.

Summary and Future Work
In accordance with the five main goals of our study (i) to (v) stated in the Introduction, the main results are summarized as follows: (i) Momentum and heat transfer coefficients based on the three pairs of theoretical SCFs of ZE07-I, ZE07-II and ZEKRE13 were derived, analyzed and compared with six pairs of the normalized transfer coefficients derived earlier.
1. The normalized transfer coefficients for momentum f m based on ZE07-I SCFs are long-tail functions, while all other normalized transfer coefficients are short-tail functions. 2. Relative to observations over sea ice (SHEBA data) the parametrizations of the normalized transfer coefficients based on the theoretical SCF of ZE07-I and ZEKRE13 are superior to many other empirical ones derived from other datasets. They are very similar to the coefficients of GLGS20, both qualitatively and quantitatively.
(ii) Non-iterative parametrizations for these new transfer coefficients are established.
1. The parametrizations are derived in a universal framework relying on an approximate semi-analytical solution of the MOST equations. These have a simple analytical functional form, which is well suited for practical use in numerical weather prediction and climate models. 2. Results from the non-iterative parametrizations have sufficient accuracy relative to results from iterative schemes. For all normalized transfer coefficients the difference between the exact iterative solution and the parametrization is smaller than the scatter in measurements as is obvious from the SHEBA data. 3. The new transfer coefficients are included as members of an extended package of non-iterative parametrizations suggested earlier by GLSG21 (see Table 1).
(iii) The verification of 18 members of the package versus SHEBA data is performed.
1. The results confirm previous findings that there are large differences in the performance of the transfer coefficients based on different pairs of stability functions relative to the SHEBA observations. This is not a failure of the parametrizations adjusted to other datasets, but it demonstrates the large variability of datasets, which until today cannot be fully explained. 2. For sea ice covered regions, the transfer coefficients based on the theoretical SCFs of ZE07-I and ZEKRE13 are superior to ones based on empirical SCFs of BD, HB88, BH91, and CB05, and only slightly worse in comparison to those based on empirical SCFs of GLGS20, which were derived on the basis of SHEBA data.
(iv) Modified stability functions based on stability functions of ZE07-I, ZE07-II, and ZEKRE13 are introduced. The modified functions improve the agreement of the corresponding normalized transfer coefficients with SHEBA data. (v) Based on our results we can recommend the NTCs by GLSG20 over sea ice, but also the modified NTCs by ZE07-I/GL for practical use in NWP and climate models, when traditional single mode parametrization schemes are used. But the entire package is recommended, if an aggregation scheme is available in the model.
Overall, the new package of non-iterative parametrizations is well suited for a systematic study of the impact of a different treatment of the turbulent surface layer (represented by different SCFs) in NWP and climate models. Also, the comparison versus SHEBA data of all members of the package opens the door for a possible tuning of NTCs in future studies, to improve turbulent fluxes over sea ice when a model uses just one member of the scheme for the flux determination everywhere.
An implementation of the package of universal non-iterative parametrizations in current Earth system models, especially in those relying on aggregated schemes of transfer coefficients, is possible (GLSG21, Schneider et al. 2022). Work in this direction is in progress (Khosravi et al. 2020).

Epilogue
We both had the privilege to work in close contact with S. S. Zilitinkevich during his long stay at Alfred-Wegener-Institute Bremerhaven from 1994-1997. Sergej always generated a warm, friendly and productive working atmosphere. He was not only an outstanding researcher, but also an excellent teacher. This work, following his clear methodology and developing some of his brilliant ideas, is dedicated to the bright memory of Sergej.

ZE07-I SFs:
The ZE07-I SFs m for momentum and h for heat in the traditional functional form read as: with the coefficients a m = 5, Pr 0 = 0.85, a h = 4, b h = 1.25.
ZEKRE13 SFs: The SFs h (ζ ) and h (ζ ) of ZEKRE13 in the original notations read where the coefficients a 1 = 0.18, a 2 = 0.16, a 3 = 1.42, Ri ∞ = 0.25 (41) are based on the coefficients of the EFB second-order closure model (see their Eqs. 83-85 with the values of the empirical constant C 0 = 0.125, C r = 1.5, Ri ∞ = 0.25). The von Kármán constant has the classical value κ = 0.4, but the constant for temperature is κ t = 0.5. The correction of the original SF φ m function (39) is the same as for ZE07-I and leads to the similar equation for SCF ψ m (18a). However, a m = 1/Ri ∞ = 4.0 ( 4 2 ) in that case. The traditional SF φ h (ζ ) follows after correcting Eq. 39 according to the definition of the Obukhov length (2) and introducing the neutral-limit Prandtl number Pr 0 = κ/κ t instead of the von Kármán constant for temperature κ t , all similar as for ZE07-I and ZE07-II. The result reads as: where the new coefficients are: Pr 0 = 0.8, c 1 = a 1 /κ = 0.45, c 2 = a 2 /κ 2 = 1.0, c 3 = a 3 /κ = 3.55.
Plugging SF (43) into Eq. 10 and calculating the integral, we derive the SCF ψ h (ζ ) of the functional form given by Eq. 18b in the main text, where the coefficients a h , b h , c h , d h are related to the coefficients c 1 , c 2 , c 3 , and c 4 as: The numerical values of these constants follow after using the constants (41) in Eq. 44 and substituting the result in Eqs. 45a-45d.
In the limit ζ → 0 (Businger-Dyer limit) the equation for ψ h (18b) is reduced to the equation of the BD SCF (11), but with a h = 4.45. Also, if c h = 0 in Eq. 18b, the equations for the SCFs of ZEKRE13 are reduced to the equations of ZE07-I, but with different values of constants a m = 4.0, a h = 4.45 and b h = 1.13. These two facts are widely used in the main text.
Finally, it is worth to note that the SFs of Zilitinkevich and co-autors are used very often without proper correction to the Obukhov length and neutral-limit Prandtl number. For example, Casasanta et al. (2021) overlooked the factor Pr 0 (see their Eq. 7b). Tastula et al. (2015) overlooked the factors Pr 0 and also the necessity to renormalize the stability parameter ζ (see their Eq. 7). Although this does not lead to qualitative changes of the results, this causes quantitative errors in the estimation of turbulent fluxes.
4. The NTCs f m and f h can be calculated now as functions of Ri b using ζ from step 3 in Eqs. 7 and 8.
This results finally in a package of different non-iterative parametrizations, where the package members differ by the functional form of the SCFs and by values of two constants γ and ζ a given in Table 1.
The algorithm is well suited for an implementation in numerical models, in particular because the explicit formulas have simple analytical forms, see, e.g. Eqs. 22, 23 by Gryanik and Lüpkes (2018), and Eqs. 55, 56 by Gryanik et al. (2020). Furthermore, the algorithm is flexible in the sense that all parametrizations are derived by the same method and have the same universal functional form for all SCFs, so that new functions can be easily added to the package. Finally, it is worth noting that the new algorithm is more general in comparison with GLSG21 because the assumptions Pr 0 ≈ 1 and ma (ζ / m ) ψ ma (ζ ) and ha (ζ / t ) ψ ma (ζ ) are not used.
For the treatment of the turbulent flux of specific humidity we refer to GL18 (see their Appendix).