Here, evaluations will be performed on the turbulent fluxes calculated by algorithms from BCC_CSM, ECMWF and CESM. Each model contains a sea-ice module, in which the turbulent fluxes over sea ice are parameterized using a bulk algorithm. These algorithms are widely used by weather prediction and climate models. In the context of the Monin–Obukhov similarity theory, the transfer coefficients for momentum (C _m ), sensible heat (C _h ) and latent heat (C _q ) obey 1 C m = k 2 [ ln ( z z 0 ) - ψ m ] 2 2 C h = k 2 [ ln ( z z 0 ) - ψ m ] [ ln ( z z T ) - ψ h ] 3 C q = k 2 [ ln ( z z 0 ) - ψ m ] [ ln ( z z Q ) - ψ h ]

where k (=0.40) is the von Kármán constant and z is the height of the lowest model layer. The key parameters to estimate C _m , C _h and C _q are the roughness lengths z ₀, z _T and z _Q , as well as the stability functions ψ _m and ψ _h .

All of the models employ the Paulson (1970) functions for ψ _m and ψ _h in unstable stratification and z ₀, z _T and z _Q are set to the same constant in each model. However, they differ in detail from each other. BCC_CSM assumes z ₀=z _T =z _Q =0.5 mm and uses the classic functions from Businger et al. (1971) in stable stratification (see Table 1). ECMWF doubles the roughness lengths (z ₀=z _T =z _Q =1 mm) and uses the functions from Holtslag and de Bruin (1988) for ψ _m and ψ _h in stable stratification (see Table 1). Although CESM uses the same stability functions that ECMWF does, it has the same values of roughness lengths as BCC_CSM. Moreover, CESM includes a windless transfer coefficient for sensible heat in stable conditions (Jordan et al. 1999). In addition, CESM offers an option to calculate the turbulent fluxes based on constant exchange coefficients (hereafter referred to as CESM_const). CESM_const assumes C _m =C _h =1.2×10⁻³ and C _q =1.5×10⁻³.

For a more extensive evaluation, Table 1 presents a collection of stability functions used for stable stratification. The classic log-linear relation is used in BCC_CSM (Businger et al. 1971%3B hereafter referred to as B71). Andreas (2002) assessed several representative functions and recommended the functions developed by Holtslag and de Bruin (1988; hereafter referred to as HDB88), which are used in ECMWF and CESM. Beljaars and Holtslag (1991; hereafter referred to as BH91) altered the HDB88 functions in consideration of different transfer efficiencies between momentum and heat. Zeng et al. (1998; hereafter referred to as Z98) expressed the functions in separate stability regimes. It uses the classic B71 relations along with relations from Holtslag et al. (1990) under very stable conditions. The functions proposed by Cheng & Brutsaert (2005; hereafter referred to as CB05) also take different forms for wind speed and temperature under both weak and very stable conditions.

In terms of roughness lengths, the only theoretically based model that specifically predicts z _T and z _Q over ice- and snow-covered surfaces was proposed by Andreas (1987). This expresses the scalar roughness z _s as a function of the roughness Reynolds number R _* (=u _* z ₀/v), 4 ln ( z s / z 0 ) = b 0 + b 1 ( ln R   ) + b 2 ( ln R   ) 2

where u _* is the friction velocity, v is the kinematic viscosity of air and z _s is either z _T or z _Q . Andreas (1987, 2002) tabulates the polynomial coefficients b ₀, b ₁ and b ₂. This model has been adequately validated and is therefore adopted in the proposed algorithm in this study.

In this section, the dimensionless gradient function (φ _m , φ _h ) forms of the stability functions listed in Table 1 are compared. Fig. 4 shows plots of these φ _m and φ _h functions. For 0<ζ<1, the five sets of functions show minor differences and approximately increase with ζ in a linear manner. However, as the stratification increases, the φ _m and φ _h functions show significant discrepancies. As for B71, the classic φ _m and φ _h functions clearly make the fastest increase with ζ among any of the sampling functions. HDB88 curves are approximately constant in a flattening range, ca. 3<ζ<7, but continue increasing slowly with ζ. BH91 curves are altered versions of the HDB88 φ _m and φ _h functions. In contrast, the φ _m function of BH91 shows a steeper gradient than HDB88. Compared with φ _m , φ _h curves have larger discrepancies between BH91 and HDB88. Z98 curves are the same as B71 in the range of ζ<1, and show moderate gradients at large ζ. CB05 curves are closer to HDB88 and approach constants at large ζ. Similar plots were considered in the literature (Grachev et al. 2007, their Figs. 4 and 4; Grachev et al. 2008, their Fig. 9; Baas & de Roode 2008, their Fig. 1). The B71 functions are predicted by the Monin–Obukhov similarity theory, which is based on the Kolmogorov turbulence. The applicability of the Monin–Obukhov similarity theory in the stable boundary layer is limited when the stability is below a critical value (Grachev et al. 2013). However, measurements of atmospheric turbulence made during SHEBA contain a large proportion of non-Kolmogorov turbulence data in very stable conditions. Therefore, both φ _m and φ _h functions derived from the SHEBA data increase more slowly with increasing ζ than is predicted by the B71 linear equations.

In Fig. 5, similar behaviour can be seen. For both ψ _m (ζ) and ψ _h (ζ), again the usual B71 functions imply very large values (in magnitude) in very stable conditions but should not be extrapolated for large ζ. Z98 functions, in contrast, were formulated specially to treat very stable stratification.

To judge which of these functions have realistic behaviour, Table 2 lists the values of the four profile metrics, D _m , D _h , R _i and Pr _t , in the limit of large ζ. Andreas (2002) recommended HDB88 functions for representing stratification effects in surface-layer wind speed, temperature and humidity profiles in stable conditions. All metrics reach reasonable limits with HDB88 as ζ increases. In addition to HDB88, Z98 functions also imply reasonable values. In the limit of large ζ, the Deacon numbers D _m and D _h both reaches zero, the limit for laminar flow. The turbulent Prandtl number is always 1, which is the approximate order of the molecular Prandtl number (ca. 0.71). Z98 also predicts a bounded critical Richardson number of 1, in line with Andreas's (2002) discussion that the critical value of R _i could be of order one, despite the existing controversy. B71 functions were never intended to treat very stable stratification. As for BH91, the D _m number, 0, agrees with the laminar limit, but the D _h number, −0.5, does not. BH91 functions also imply unbounded Richardson and Prandtl numbers, both increase as ([2/3]ζ)^1/2. The functions of CB05 predict unusual Deacon numbers of D _m =D _h =1, contrary to the predictions from the laminar flow theory. The CB05 functions also do not produce a critical Richardson number; R _i increases as 0.18ζ. Figs. 6 and 7 suggest large differences among the Deacon number curves as stability increases. The Z98 curves deviate from other curves and seem to monotonously decrease with R _i .

Besides theoretical analysis, an intercomparison was performed based on the observational data. Fig. 8 shows the hourly wind stresses and heat fluxes binned into 0.5K bins of the air-surface potential temperature difference Δθ based on the SHEBA 20-m tower data. The fluxes from the five stability functions with constant roughness lengths (z ₀=z _T =z _Q =0.5 mm) are compared. Positive Δθ values represent stable atmospheric stratification. For wind stresses, observed values decrease as Δθ increases because of turbulence damping by the stable stratification. All stability functions reproduce this stratification effect with little difference.

The magnitude of observed sensible heat flux generally increases with increasing Δθ to a maximum at about 2K and then decreases to a very small value in the very stable regime. The bulk sensible heat fluxes derived from the five stability functions are systematically lower than observed. The relatively large discrepancies between these stability functions are presented in very stable conditions. However, as shown in Fig. 9, the greatest numbers of observations are within Δθ=±3K. Values of Δθ greater than +3K occur infrequently. Differences at very stable Δθ would therefore not contribute much to the mean bulk results of sensible heat flux. Observed latent heat fluxes from the 20-m tower only are positive for weakly stable conditions (Δθ<2K), and near zero for very stable stratification. Unfortunately, all of the five functions produce negative fluxes under stable conditions.

Fig. 10 presents monthly averages of the observed turbulent fluxes at the SHEBA 20-m tower, based on the median value observed at the five levels for each hour during the SHEBA year. There is significant seasonality in the wind stress, which is around 0.06 N m⁻². Sensible heat flux is generally above 50 W m⁻² in the terrestrial mid-latitudes and latent heat fluxes could be as large as 200 W m⁻² over tropical oceans (Brunke et al. 2002; Jiménez et al. 2009). In contrast, the magnitudes of both sensible and latent heat fluxes are very small over sea ice in the Arctic. During the non-melting season, the average sensible heat flux is about −5 W m⁻². As the surface temperature rises during the melting season, the sensible heat flux increases and becomes slightly positive in May and August. The data recovery for latent heat flux is low in February, March and September (Persson et al. 2002). In winter, the average latent heat flux is near zero, and it can only reach as high as 2 W m⁻² in summer.

The bulk fluxes obtained from the four model algorithms are shown in Fig. 10. Wind stresses (Fig. 10a) from BCC_CSM and CESM are close to SHEBA observations, while ECMWF's wind stresses are slightly overestimated and CESM_const's wind stresses are badly underestimated. Overestimation by the ECMWF algorithm may be due to enlarged roughness length (1 mm) in the bulk parameterization. Brunke et al. (2006) evaluated the ARCSYM model algorithm using the SHEBA data, in which the roughness length is 60 mm for snow depths less than 5 cm and 40 mm for greater snow depths. They found wind stresses from the ARCSYM model were significantly overestimated compared to observations and other model algorithms. The constant C _m (1.2×10⁻³) in CESM_const corresponds to a roughness length of 0.1 mm, which is at least five times smaller than any other model's roughness length. Monthly mean observed and bulk sensible heat fluxes agree well, even though the bulk fluxes are slightly lower than observed (Fig. 10b). In contrast to the sensible heat flux, the latent heat flux is overestimated by the bulk algorithms in May and June (Fig. 10c). Overestimation of the latent heat flux is also found in similar studies of bulk algorithm intercomparisons based on the SHEBA data (Persson et al. 2002; Brunke et al. 2006). Relatively small values of observed sensible and latent heat fluxes represent a challenge to determine accurate fluxes from bulk algorithms. Even small uncertainties in air–surface temperature and humidity measurements will produce fluxes with uncertainties of ±50%. The inadequacies in measurement technology account for a large part of these discrepancies.

Here, we also used the CHINARE 2010 data to evaluate surface turbulent fluxes from the model algorithms. As shown in Fig. 3, the CHINARE 2010 measurements were obtained mostly in unstable conditions in summertime. The bulk fluxes from BCC_CSM, ECMWF and CESM show little discrepancies, which may be because all the models employ the Paulson (1970) functions for ψ _m and ψ _h in unstable stratification. Fig. 11 shows the simulated fluxes from BCC_CSM along with the corresponding CHINARE 2010 observations. It can be seen that the wind stresses and sensible heat fluxes are generally simulated well. With respect to the latent heat flux, the mean observed value is only about 2.04 W m⁻², but the average algorithm-derived value is up to 12.6 W m⁻².

For more quantitative intercomparisons, Table 3 presents the biases and RMSEs of the four model algorithms to the observed fluxes from SHEBA 20-m tower. The CESM_const algorithm has the largest RMSEs in t (3.94×10⁻² N m⁻²) and sensible heat flux (6.89 W m⁻²), as should be expected since the CESM_const algorithm uses the constant exchange coefficient (1.2×10⁻³ for both momentum flux and latent heat flux) with no stability dependence to obtain the fluxes.

In contrast, the RMSE in latent heat flux from the CESM_const algorithm is comparable to the other three algorithms, though the exchange coefficient is stability independent as well (1.5×10⁻³ for latent heat flux). The bias between computed and observed t is smallest using the CESM algorithm, whereas the bias with respect to sensible heat flux is smallest using the BCC_CSM algorithm.

Of the stability-based algorithms, ECMWF and CESM employ the Holtslag & de Bruin (1988) stability function, while BCC_CSM uses the Businger et al. (1971) function. Furthermore, different roughness lengths are used among these three algorithms (0.5 mm in BCC_CSM and CESM, and 1 mm in ECMWF). Based on the same stability function, the CESM algorithm performs better than the ECMWF algorithm with doubled roughness lengths with respect to t and latent heat flux. However, the BCC_CSM algorithm performs better than CESM in sensible heat flux, both having the same roughness length but differentiating between stability functions. CESM has more biases of negative sensible heat fluxes. Table 3 suggests that the computed fluxes are very sensitive to the representation of the stable conditions and roughness lengths.

Andreas and co-workers developed a new bulk turbulent flux algorithm over sea ice for winter based on the SHEBA data set (Andreas, Persson et al. 2010), and they proposed a companion algorithm for summer (Andreas, Horst et al. 2010). The distinctive feature of these algorithms is incorporation of both new stability functions and varying roughness lengths. Here, we also propose an alternative algorithm by introducing the Z98 stability functions and incorporating Andreas's (1987) theoretical scalar roughness model. Values of z ₀ tend to be less than 0.5 mm over sea ice in winter and are prone to be greater than 0.5 mm in summer (Andreas, Horst et al. 2010; Andreas, Persson et al. 2010). For simplicity, z ₀ is treated in the new algorithm as follows: 12 z 0 ( mm ) = { 0.1 for T s ≤ - 2 ° C 0.8 for T s > - 2 ° C

To examine the new algorithm's performance, Figs. (12–14) show scatterplots of observed u _*, H _s and H _l, along with modelled values from the new, BCC_CSM, ECMWF and CESM algorithms, respectively. The u _* and H _s values are based on the CHINARE 2010 data and the Atlanta data and the H _l values are based on the CHINARE 2010 data and the 20-m tower data. In addition, Table 4 lists the metrics from the Atlanta data and the 20-m tower data for a quantitative comparison with the results of Andreas and co-workers (Andreas, Horst et al. 2010; Andreas, Persson et al. 2010) in winter (October 1997–14 May 1998 and 15 September 1998 to the end of the SHEBA deployment in early October 1998) and summer (15 May–14 September 1998), respectively. We chose to use the Atlanta data in this section for a comparison with the results of Andreas and co-workers, which are based on the same data set.

In Fig. 12 for the friction velocity u _*, the fitting line for the new algorithm is visually closer to the 1:1 line than those fits based on the other three models. In winter (Table 4), the bias error for the new algorithm is 84% smaller than that for the BCC_CSM algorithm, and the RMSE is 20% smaller. Errors associated with the Andreas, Persson et al. (2010) algorithm are even smaller but it is based on the same Atlanta data used to develop the algorithm, while our new algorithm uses independent stability functions and simpler z ₀ treatments.

Table 4 shows that in summer, the bias errors are also smaller with the new algorithm than that with the BCC_CSM model, by 56%. The RMSE of the new algorithm is, however, slightly larger than BCC_CSM. This is mainly because, in summer, the factors affecting z ₀ are more complex than in winter, and simply parameterizing z ₀ to an increased constant (0.8) does not adequately represent the drag effects. Andreas, Horst et al. (2010) developed physically based parameterizations to describe the behaviour of the drag coefficient with the fractional ice concentration. The algorithm with this parameterization for z ₀ provides better results presented in Table 4.

For the sensible heat flux H _s, modelled values from BCC_CSM, ECMWF and CESM have large negative biases shown in Fig. 13. The agreement between modelled and observed sensible heat flux for the new algorithm is noticeably better than for other models. H _s values tend to be closer in magnitude to observations in stable stratification. In other words, negative fluxes are modelled more accurately with the new algorithm. The fitting line is thus rotated clockwise from the 1:1 line, as we see in Fig. 13a. The values for bias and RMSE in Table 4 confirm the better performance of the new algorithm. The bias error for the new algorithm is 30% smaller than for the BCC_CSM algorithm, and the RMSE is 17% smaller in winter. In summer, the bias error for the new algorithm is 19% smaller than for the BCC_CSM algorithm and the RMSE is 18% smaller. Our new algorithm presents good performances comparable to the results of Andreas and co-workers (Andreas, Horst et al. 2010; Andreas, Persson et al. 2010) results.

Latent heat flux H _l is compared with the data from the CHINARE 2010 data and the 20-m tower. Fig. 14 shows that the overestimation of H _l with BCC_CSM, ECMWF and CESM is reduced by the new algorithm, though the scatter in Fig. 14a is not very different from the scatter in Fig. 14b–d. Again, the values for bias and RMSE in Table 4 reiterate the better performance of the new algorithm; the bias is 27% smaller and the RMSE value is 21% smaller than BCC_CSM. The results of Andreas, Horst et al. (2010) have even smaller errors, which may be attributed to the more complex z ₀ parameterization and the stability functions from Grachev et al. (2007) derived from the same data set.

This work aims to evaluate the parameterizations for the turbulent surface fluxes over sea ice, which features distinctive treatments of roughness lengths and correction functions for stable conditions. Five sets of stability functions are compared with the theoretically expected values. Four metrics, D _m , D _h , R _i and Pr _t , in the limit of laminar flow are used to judge which of these functions have realistic behaviour. Andreas (2002) made a comprehensive assessment of different stability functions using this technique and recommended the expressions of HDB88 due to their best performance in the four metrics. All of these metrics of HDB88 reach reasonable limits as ζ increases: D _m and D _h both reaches zero; the limit for laminar flow (Pr _t ) is always 1, which is the approximate order of the molecular Prandtl number; R _i is 1.43, a critical Richardson number of order 1. This study suggests that the Z98 functions could provide theoretically expected values also. Both the Deacon numbers D _m and D _h are zero. The turbulent Prandtl number of Z98 is 1, which is the same as that of HDB88. Finally, Z98's functions predict that the critical Richardson number is 1.

Compared with observations in different stability regimes, bulk fluxes derived from the functions of Z98 are comparable to that from HDB88's functions. For very stable conditions (Δζ>2K), the bulk fluxes from Z98's functions are even closer to observations than those from HDB88's functions. These results imply that the stability functions of Z98 can be the best expressions for representing stratification effects in stable conditions over sea ice.

Based on the observational data from the CHINARE 2010 experiment and the 20-m tower data from the SHEBA experiment, the bulk algorithms from BCC_CSM, ECMWF and CESM are compared. Results show that these algorithms generally perform better in calculating the momentum flux than the sensible and latent heat fluxes. Based on the SHEBA 20-m tower data, the biases of calculated wind stresses are very small in BCC_CSM (ca. 0.63%) and CESM (ca. 0.16%). The doubled roughness lengths in the ECMWF algorithm increase the relative bias to 16.8%; the bulk wind stresses are very sensitive to the roughness length specified.

With respect to the heat fluxes, these algorithms generally produce underestimated sensible heat fluxes and overestimated latent heat fluxes. For sensible heat flux, the models are biased low by 1.6–2.7 W m⁻². Because the observed mean is small (−1.92 W m⁻²), these biases seem to be relatively too large.

Modelled latent heat fluxes have large discrepancies in May and June. These errors may be attributed to the different assumptions of constant roughness lengths in models. Measurements are rare over sea ice and the correction functions for stable stratification used in most current models are based on measurements in the traditional nocturnal boundary layer in mid-latitudes. Furthermore, some unknown physics may account for these biases.

We propose a new algorithm focusing on improving flux parameterizations over sea ice with Eqn. 4 for varying scalar roughness lengths and the stability functions of Z98. The roughness lengths for wind speed in the new algorithm differentiate between surface temperatures instead of being always constant in BCC_CSM. For surface temperature less than −2°C, z ₀ is 0.1 mm, but is increased to 0.8 mm for warmer sea-ice surfaces. When compared with BCC_CSM, ECMWF and CESM, the results from the new algorithm confirm that these are improvements. For friction velocity u _*, the bias of the new algorithm is 0.0077 m s⁻¹ in winter and 0.0014 m s⁻¹ in summer, which is about 84% and 56% smaller than that of BCC_CSM, respectively. Downward sensible heat fluxes are modelled more accurately with the new algorithm, with significantly reduced biases 30% in winter and 19% in summer smaller than BCC_CSM, respectively. Finally, for the latent heat flux, the new algorithm has a bias error of only 1.71 W m⁻² in summer, 27% smaller than BCC_CSM. These results are comparable to the results of Andreas and co-workers (Andreas, Horst et al. 2010; Andreas, Persson et al. 2010), which are based on representations of roughness lengths and stability functions derived from the same SHEBA data.