RESEARCH/REVIEW ARTICLE

Validation of parameterizations for the surface turbulent fluxes over sea ice with CHINARE 2010 and SHEBA data

Yixiong Lu,¹ Mingyu Zhou² & Tongwen Wu¹

¹ Beijing Climate Centre, China Meteorological Administration, Beijing, China
² National Marine Environmental Forecasting Centre, Beijing, China

Abstract

This study examines the modelled surface turbulent fluxes over sea ice from the bulk algorithms of the Beijing Climate Centre Climate System Model (BCC_CSM), the European Centre for Medium-Range Weather Forecasts (ECMWF) model and the Community Earth System Model (CESM) with data from the fourth Chinese National Arctic Research Expedition (CHINARE 2010) and the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment. Of all the model algorithms, wind stresses are replicated well and have small annual biases (−0.6% in BCC_CSM, 0.2% in CESM and 17% in ECMWF) with observations, annual sensible heat fluxes are consistently underestimated by 83–141%, and annual latent heat fluxes are generally overestimated by 49–73%. Five sets of stability functions for stable stratification are evaluated based on theoretical and observational analyses, and the superior stability functions are employed in a new bulk algorithm proposal, which also features varying roughness lengths. Compared to BCC_CSM, the new algorithm can estimate the friction velocity with significantly reduced bias, 84% smaller in winter and 56% smaller in summer, respectively. For the sensible heat flux, the bias of the new algorithm is 30% smaller in winter and 19% smaller in summer than that of BCC_CSM. Finally, the bias of modelled latent heat fluxes is 27% smaller in summer.

Keywords
Ice–atmosphere coupling; turbulent flux parameterization; Monin–Obukhov similarity theory; SHEBA; stable boundary layer.

Correspondence
Yixiong Lu, Beijing Climate Centre, 46 Zhongguancun Nandajie, Beijing 100081, China. E-mail: luyx@cma.gov.cn

(Published: 23 September 2013)

Polar Research 2013. © 2013 Y. Lu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial 3.0 Unported License (http://creativecommons.org/licenses/by-nc/3.0/), permitting all non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Citation: Polar Research 2013, 32, 20818, http://dx.doi.org/10.3402/polar.v32i0.20818

Parameterizing surface turbulent fluxes is particularly important in coupled models that are used for seasonal forecasting and climate simulation (Ahlgrimm & Randall 2006). In the polar regions, the atmosphere and sea-ice interplay through the exchange of surface turbulent fluxes. Momentum flux forces sea ice to move and in turn, drives ocean currents and surface mixing. It also creates pressure ridges and leads of open water, and moves the snow around on top of sea ice. Turbulent heat fluxes dictate the thermal structure of the lower atmosphere and thereby crucially influence the momentum flux and general atmospheric circulation (Andreas, Persson et al. 2010).

For unstable conditions, there is some controversy about parameterizing surface turbulent fluxes. Over sea ice, however, the atmospheric surface layer is mostly stably stratified. Although much effort has been made to research the stable boundary layer, a unified picture or theory does not exist. The weakly stable boundary layer has been extensively studied based on measurements in the typical nocturnal boundary layer in mid-latitudes. For example, the log-linear relation (Businger et al. 1971; Dyer 1974) for stratification effects is commonly used in many global climate models, even for very stable conditions. However, the processes of the very stable boundary layer are very complex. It is characterized by intermittent turbulence, meandering motions, gravity waves and drainage flows (Grachev et al. 2005). Many nonlinear models have been developed for the very stable stratification (Holtslag & de Bruin 1988; Cheng & Brutsaert 2005). Grachev et al. (2007) developed new mathematical forms for stable conditions based on the field experiment known as Surface Heat Budget of the Arctic Ocean (SHEBA), which made comprehensive measurements of the stable surface layer and collected ample data for the very stable case. Up to now, the issue of how the stability functions perform for bulk parameterization over sea ice has still received little attention (Andreas 2002; Brunke et al. 2006).

Another crucial issue involving parameterizing surface turbulent fluxes over sea ice is how to treat the roughness length for wind speed (z₀), and the so-called scalar roughness length for temperature (z_T ) and humidity (z_Q ). As is well known, there is a strong seasonality to the roughness lengths, resulting from the surface characteristics of sea ice. In winter, the ice is compact and snow covered, and the snow is dry enough to blow and drift. In summer, the deposited snow melts, and the resulting melt ponds and leads create vertical surfaces that increase momentum roughness (Andreas, Horst et al. 2010; Andreas, Persson et al. 2010). Pressure ridges also change the roughness as compared to undeformed ice (Brunke et al. 2006).

Parameterizing the roughness lengths of sea ice is complicated, and efforts have been made to propose mathematical formulations. The roughness lengths of sea ice are generally considered to be constant and equal to each other in climate models. Birnbaum & Lüpkes (2002) and Lüpkes & Birnbaum (2005) developed physically based parameterizations for the surface drag of summer sea ice, but the models are too complex to be used in global climate models. Andreas et al. (2005) derived a parameterization for snow-covered sea ice using data from Ice Station Weddell that combine three regimes. Brunke et al. (2006) offered a fit of an exponential form based on the SHEBA data set that is simpler than the Andreas et al. (2005) equation, but the results were considered to be fictitious correlation caused by the analysis method. To minimize the effects of fictitious correlation, Andreas and co-workers (Andreas, Persson et al. 2010) gave an alternative formulation for winter based on the same SHEBA data set. A companion formulation was proposed for summer (Andreas, Horst et al. 2010). However, these parameterizations do not have a unified form applicable throughout the year, and may need to be tuned for a particular location accounting for local sea-ice topography. Simple and practical z₀ parameterizations for global climate models still await an intensive study. In contrast, finding z_T and z_Q over sea ice is not as complicated and there is a theoretically based model developed by Andreas (1987) for which predictions agree well with observations (Denby & Snellen 2002; Andreas et al. 2005; Andreas, Horst et al. 2010; Andreas, Persson et al. 2010).

The objective of this study is to investigate the role of stable stratification and roughness lengths in parameterizing surface turbulent fluxes over sea ice based on the data from the fourth Chinese National Arctic Research Expedition (CHINARE 2010) and the SHEBA field campaign. First, the stability functions and roughness length formulations are described and the observational data and verification methods used in this study are then described. An intercomparison is made based on five sets of stability functions and the bulk algorithms widely used by global climate models are examined. These models are the Beijing Climate Centre Climate System Model (BCC_CSM), the European Centre for Medium-Range Weather Forecasts model (ECMWF) and the Community Earth System Model (CESM). A new algorithm characterized by the superior stability functions and varying roughness lengths is tested. Finally, a summary and discussion are presented.

Bulk aerodynamic algorithms

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

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⁻³.

Table 1 A summary of five stability functions used for stable conditions.

Algorithm	Reference	Stability functions	Comments
B71	Businger et al. 1971		Used in the Beijing Climate Centre Climate System Model
HDB88	Holtslag & de Bruin 1988	where a=0.7, b=0.75, c=5 and d=0.35	Used in the European Centre for Medium-Range Weather Forecasts model and the Community Earth System Model
BH91	Beljaars & Holtslag 1991	where a=1, b=0.667, c=5 and d=0.35	A revision of HDB88
Z98	Zeng et al. 1998		Used in the new algorithm in this article
CB05	Cheng & Brutsaert 2005	where a=6.1, b=2.5, c=5.3 and d=1.1

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; 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),

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.

Data and methods

Observational data

As part of CHINARE 2010, an ice station was established close to the North Pole and the atmospheric surface layer was monitored for about one week. Serving as the platform for the expedition, the RV Xuelong sailed up to 87°N and an ice station was deployed on a large ice floe, north of the SHEBA region (Fig. 1). The hourly turbulent fluxes used in this study were measured on a 3-m tower for about one week (12–19 August). Fast measurements of wind velocities and temperature were made using a CSAT3 three-dimensional sonic anemometer/thermometer at a height of 2.65 m (Campbell Scientific, Logan, UT). The data were rotated in three dimensions to stream-wise flow. At the same level, the tower also held a LI-7500 fast-response humidity sensor (Li-Cor Biosciences, Lincoln, NE). The sensible heat fluxes were obtained applying small moisture corrections. A minor correction for the Webb effect (Webb et al. 1980) was made to the flux data. In addition, the data were quality controlled on the basis of objective and subjective criteria including validity limits. Fluxes were removed if the data were flagged as bad. The air temperature and relative humidity were measured with HMP45C sensors (Vaisala, Vantaa, Finland) at a height of 2 m. The upward and downward components of the shortwave and longwave radiation were measured by the CNR1 net radiometer (Kipp & Zonen, Delft, The Netherlands). Surface temperatures were derived from the radiometric measurements of the emitted (R_lu ) and incoming (R_ld ) longwave radiation:

Fig. 1  The fourth Chinese National Arctic Research Expedition (CHINARE 2010) ice station drift indicated by the line (A→B), along with the Surface Heat Budget of the Arctic Ocean (SHEBA) ice camp drift (X→Y). The CHINARE 2010 drift lasted from 12 to 19 August 2010, and the SHEBA drift started in October 1997 and ended in October 1998.

Here, ɛ (=0.99) is the surface emissivity and σ (=5.67051×10⁻⁸ W m⁻² K⁻⁴) is the Stefan–Boltzmann constant.

Figure 2 shows the near-surface observations during the CHINARE 2010 period. The mean air temperature was −1.5°C with a range from −4.8 to 1.1°C. After 17 August, the wind speed increased caused by a low in the New Siberian Islands. There was a snowfall event that began in the afternoon of 16 August, and continued until the evening of 17 August. This was associated with unusually high air temperature. The diurnal cycle of both upward and downward shortwave radiation was pronounced, while the diurnal cycle of longwave radiation was absent. In this article, upward (or downward) fluxes are defined to be positive (or negative).

Fig. 2  Hourly mean (a) air pressure, (b) surface and air temperature, (c) relative humidity (with respect to saturation over water), (d) wind speed and (e) radiation measurements made during the fourth Chinese National Arctic Research Expedition. The downward and upward shortwave radiations are indicated by Sd and Su, respectively, and the downward and upward longwave radiations are indicated by Ld and Lu, respectively. Rn represents the net radiation.

Detailed descriptions of the SHEBA programme can be found in Andreas et al. (1999), Persson et al. (2002), Uttal et al. (2002) and Persson (2012). Only the relevant information about the SHEBA sites and flux data is given here. The SHEBA field experiment began in October 1997 around the Canadian Coast Guard ice-breaker Des Groseilliers, which was frozen into the Arctic ice pack north of Alaska at 75°N, 142°W, and ended in October 1998 at about 80°N, 160°W, as shown in Fig. 1. Hourly turbulent flux measurements and mean meteorological observations from the 20-m tower (http://data.eol.ucar.edu/codiac/projs?SHEBA) and from the portable, automated mesonet station named Atlanta (http://www.eol.ucar.edu/isf/projects/sheba) are used in this study. The 20-m tower collected turbulent and mean meteorological data at five levels, nominally 2.2, 3.2, 5.1, 8.9 and 18.2 m (or 14 m during most of winter). The latent heat flux was measured at only one level (8.1 m) at the tower. At the Atlanta station, where there were no measurements of latent heat flux, measurements of temperature and relative humidity were made at a height of 1.5–2 m, and a sonic anemometer/thermometer sensor was mounted at a height between 2.5 and 3.5 m.

The year-long SHEBA data offer several advantages for studying different climate regimes in detail. Figure 3 shows surface temperature versus the air–surface temperature difference (ΔT=T_a−T_s) for hourly SHEBA observations on the 20-m tower, along with the CHINARE 2010 data. In the combined data set, the temperatures are lowest in winter (December, January and February) and rise to summer levels of about 0°C from June to August. Also, the stability is mainly stable (ΔT>0) or weakly unstable while some are quite stable.

Fig. 3  The surface temperature (T_s) in spring (March—May 1998, green dots), summer (June—August 1998, red dots), autumn (October 1997 and September—October 1998, orange dots) and winter (December 1997—February 1998, blue dots) plotted as a function of the air–surface temperature difference (ΔT=T_a−T_s) using data from the Surface Heat Budget of the Arctic Ocean (SHEBA) 20-m tower. Data from the fourth Chinese National Arctic Research Expedition (CHINARE 2010) ice station are represented by black dots. T_a is the air temperature measured at a height of 2.2 m for SHEBA and 2 m for CHINARE 2010.

Verification methodology

Two metrics were computed to evaluate the performance of bulk flux algorithms. They are the bias and root-mean square error (RMSE) between the modelled and observed time series. The bias is defined as

where m_i and o_i are the modelled and observed values, respectively, for measurement i. N represents the number of observations. The RMSE in the modelled fluxes is defined as follows:

Furthermore, the method by Andreas (2002) is adopted to examine the behaviour of stability functions in Table 1 with four metrics in the limit of laminar flow. These metrics include the gradient Richardson number R_i , the Deacon numbers for wind speed D_m and potential temperature D_h , and the turbulent Prandtl number Pr_t , which are defined as follows:

where φ_m and φ_h are stability functions expressed in the form of the dimensionless gradient functions. As ζ approaches infinity, the metrics can be deduced with Eqns. 8–11 for a stability function and compared with the theoretically expected values.

Intercomparison of stability functions

In this section, the dimensionless gradient function (φ_m , φ_h ) forms of the stability functions listed in Table 1 are compared. Figure 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 figures 4, 5; Grachev et al. 2008, their figure 9; Baas & de Roode 2008, their figure 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.

Fig. 4  Intercomparison of the dimensionless gradient functions for (a) wind speed φ_m (ζ) and (b) potential temperature φ_h (ζ) for stable stratification. B71, HDB88, BH91, Z98 and CB05 denote the functions proposed by Businger et al. (1971), Holtslag & de Bruin (1988), Beljaars & Holtslag (1991), Zeng et al. (1998) and Cheng & Brutsaert (2005), respectively.

Fig. 5  Intercomparison of the stability functions for (a) wind speed ψ_m (ζ) and (b) potential temperature ψ_h (ζ) for stable stratification. B71, HDB88, BH91, Z98, and CB05 denote the functions proposed by Businger et al. (1971), Holtslag & de Bruin (1988), Beljaars & Holtslag (1991), Zeng et al. (1998) and Cheng & Brutsaert (2005), respectively.

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ζ. Figures 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 .

Fig. 6  Theoretical relationships between the Deacon number for wind speed profile D_m and gradient Richardson number R_i for (a) R_i less than 0.5 and (b) R_i greater than 0.1.

Fig. 7  Theoretical relationships between the Deacon number for potential temperature profile D_h and gradient Richardson number R_i for (a) R_i less than 0.5 and (b) R_i greater than 0.1.

Table 2 A comparison of five sets of φ_m and φ_h functions in the limit of very stable stratification. The laminar case is a theoretically expected limit.

	Lim ζ →∞
φm , φh	Dm	Dh	Ri	Prt
Businger et al. 1971	0	0	0.20	1
Holtslag & de Bruin 1988	0	0	1.43	1
Beljaars & Holtslag 1991	0	−0.5
Zeng et al. 1998	0	0	1	1
Cheng & Brutsaert 2005	1	1	0.18κ	0.89
Laminar	0	0		∼0.71

Besides theoretical analysis, an intercomparison was performed based on the observational data. Figure 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.

Fig. 8  Hourly (a) wind stresses, (b) sensible heat (SH) fluxes and (c) latent heat (LH) fluxes from Surface Heat Budget of the Arctic Ocean (SHEBA) 20-m tower observations and bulk algorithms plotted as a function of the air–surface potential temperature difference Δθ (Δθ=θ_a−θ_s) in 0.5K bins on the SHEBA 20-m tower. θ_a is the potential temperature calculated at the lowest tower level (2.2 m).

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. 9  The number of hourly observations for all the 5 levels (2.2 m, 3.2 m, 5.1 m, 8.9 m, and 18.2 m or 14 m) at the Surface Heat Budget of the Arctic Ocean 20-m tower in each of the 0.5K Δθ bins between −4 and+12K.

Model validation

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

Fig. 10  Monthly means of (a) wind stress, (b) sensible heat (SH) flux and (c) latent heat (LH) flux measured at the 20-m Surface Heat Budget of the Arctic Ocean (SHEBA) tower along with the bulk fluxes from BCC_CSM, ECMWF, CESM and CESM_const during the SHEBA year, with positive fluxes being upward.

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. Figure 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⁻².

Fig. 11  Hourly mean (a) wind stress, (b) sensible heat flux and (c) latent heat flux during the period of the fourth Chinese National Arctic Research Expedition. The in situ measurements are indicated by blue lines, and the algorithm fluxes derived from BCC_CSM are indicated by red lines.

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 τ (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 τ 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 τ 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.

Test of an alternative algorithm

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:

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.

Fig. 12  Scatterplots of friction velocity u_* measured at the fourth Chinese National Arctic Research Expedition ice station (green circles) and the portable automated mesonet station named Atlanta (red dots) and modelled with the bulk algorithms of (a) the new algorithm proposed in this article (NEW), (b) the Beijing Climate Centre Climate System Model (BCC_CSM), (c) the European Centre for Medium-Range Weather Forecasts model (ECMWF) and (d) the Community Earth System Model (CESM). In each panel, the black dashed line is 1:1. The blue line is the best fit through the data, taken as the bisector of y-vs.-x and x-vs.-y least-square fits (e.g., Andreas 2002).

Fig. 13  Scatterplots of the sensible heat flux H_s at the fourth Chinese National Arctic Research Expedition ice station (green circles) and the portable automated mesonet station named Atlanta (red dots) and modelled with the bulk algorithms of (a) the new algorithm proposed in this article (NEW), (b) the Beijing Climate Centre Climate System Model (BCC_CSM), (c) the European Centre for Medium-Range Weather Forecasts model (ECMWF) and (d) the Community Earth System Model (CESM). In each panel, the black dashed line is 1:1. The blue line is the best fit through the data taken as the bisector of y-vs.-x and x-vs.-y least-square fits (e.g., Andreas 2002).

Fig. 14  Scatterplots of the latent heat flux H_l at the fourth Chinese National Arctic Research Expedition ice station (green circles) and on the 20-m tower (red dots) and modelled with the bulk algorithms of (a) the new algorithm proposed in this article (NEW), (b) the Beijing Climate Centre Climate System Model (BCC_CSM), (c) the European Centre for Medium-Range Weather Forecasts model (ECMWF) and (d) the Community Earth System Model (CESM). In each panel, the black dashed line is 1:1. The blue line is the best fit through the data taken as the bisector of y-vs.-x and x-vs.-y least squares fits (e.g., Andreas 2002).

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

Summary and discussion

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.

Acknowledgements

This study was supported by the Ministry of Science and Technology of the People's Republic of China (grant no. 2010CB951902) and by the National Natural Science Foundation of China (grant no. 41175010). The authors are very grateful to three reviewers for their careful review and valuable comments that improved this article.

References

Ahlgrimm M. & Randall D.A. 2006. Diagnosing monthly mean boundary layer properties from reanalysis data using a bulk boundary layer model. Journal of Atmospheric Sciences 63, 998–1012. Publisher Full Text

Andreas E.L. 1987. A theory for the scalar roughness and the scalar transfer coefficients over snow and sea ice. Boundary-Layer Meteorology 38, 159–184. Publisher Full Text

Andreas E.L. 2002. Parameterizing scalar transfer over snow and ice: a review. Journal of Hydrometeorology 3, 417–432.

Andreas E.L., Fairall C.W., Guest P.S. & Persson P.O.G. 1999. An overview of the SHEBA atmospheric surface flux program. In: 5th Conference on Polar Meteorology and Oceanography, 10–15 January 1999, Dallas, Texas. Pp. 411–416. Boston, MA: American Meteorological Society.

Andreas E.L., Horst T.W., Grachev A.A., Persson P.O.G., Fairall C.W., Guest P.S. & Jordan R.E. 2010. Parametrizing turbulent exchange over summer sea ice and the marginal ice zone. Quarterly Journal of the Royal Meteorological Society 136, 927–943. Publisher Full Text

Andreas E.L., Jordan R.E. & Makshtas A.P. 2005. Parameterizing turbulent exchange over sea ice: the Ice Station Weddell results. Boundary-Layer Meteorology 114, 439–460. Publisher Full Text

Andreas E.L., Persson P.O.G., Jordan R.E., Horst T.W., Guest P.S., Grachev A.A. & Fairall C.W. 2010. Parameterizing turbulent exchange over sea ice in winter. Journal of Hydrometeorology 11, 87–104. Publisher Full Text

Baas P. & de Roode S.R. 2008. The scaling behavior of a turbulent kinetic energy closure model for stably stratified conditions. Boundary-Layer Meteorology 127, 17–36. Publisher Full Text

Beljaars A.C.M. & Holtslag A.A.M. 1991. Flux parameterization over land surfaces for atmospheric models. Journal of Applied Meteorology 30, 327–341.

Birnbaum G. & Lüpkes C. 2002. A new parameterization of surface drag in the marginal sea ice zone. Tellus A 54, 107–123. Publisher Full Text

Brunke M.A., Zeng X. & Anderson S. 2002. Uncertainties in sea surface turbulent flux algorithms and data sets. Journal of Geophysical Research—Oceans 107, article no. 3141, doi: 10.1029/2001JC000992. Publisher Full Text

Brunke M.A., Zhou M., Zeng X. & Andreas E.L. 2006. An intercomparison of bulk aerodynamic algorithms used over sea ice with data from the Surface Heat Budget for the Arctic Ocean (SHEBA) experiment. Journal of Geophysical Research—Oceans 111, C09001, doi: 10.1029/2005JC002907. Publisher Full Text

Businger J.A., Wyngaard J.C., Izumi Y. & Bradley E.F. 1971. Flux–profile relationships in the atmospheric surface layer. Journal of Atmospheric Sciences 28, 181–189.

Cheng Y.G. & Brutsaert W. 2005. Flux–profile relationships for wind speed and temperature in the stable atmospheric boundary layer. Boundary-Layer Meteorology 114, 519–538. Publisher Full Text

Denby B. & Snellen H. 2002. A comparison of surface renewal theory with the observed roughness length for temperature on a melting glacier surface. Boundary-Layer Meteorology 103, 459–468. Publisher Full Text

Dyer A.J. 1974. A review of flux–profile relationships. Boundary-Layer Meteorology 7, 363–372. Publisher Full Text

Grachev A.A., Andreas E.L., Fairall C.W., Guest P.S. & Persson P.O.G. 2007. SHEBA flux–profile relationships in the stable atmospheric boundary layer. Boundary-Layer Meteorology 124, 315–333. Publisher Full Text

Grachev A.A., Andreas E.L., Fairall C.W., Guest P.S. & Persson P.O.G. 2008. Turbulent measurements in the stable atmospheric boundary layer during SHEBA: ten years after. Acta Geophysica 56, 142–166. Publisher Full Text

Grachev A.A., Andreas E.L., Fairall C.W., Guest P.S. & Persson P.O.G. 2013. The critical Richardson number and limits of applicability of local similarity theory in the stable boundary layer. Boundary-Layer Meteorology 147, 51–82. Publisher Full Text

Grachev A.A., Fairall C.W., Persson P.O.G., Andreas E.L. & Guest P.S. 2005. Stable boundary-layer scaling regimes: the SHEBA data. Boundary-Layer Meteorology 116, 201–235. Publisher Full Text

Holtslag A.A.M., de Bruijn E.I.F. & Pan H.-L. 1990. A high-resolution air mass transformation model for short-range weather forecasting. Monthly Weather Review 118, 1561–1575.

Holtslag A.A.M. & de Bruin H.A.R. 1988. Applied modeling of the nighttime surface energy balance over land. Journal of Applied Meteorology 27, 689–704.

Jiménez C., Prigent C. & Aires F. 2009. Toward an estimation of global land surface heat fluxes from multisatellite observations. Journal of Geophysical Research—Atmospheres 114, D06305, doi: 10.1029/2008JD011392. Publisher Full Text

Jordan R.E., Andreas E.L. & Makshtas A.P. 1999. Heat budget of snow-covered sea ice at North Pole 4. Journal of Geophysical Research—Oceans 104, 7785–7806. Publisher Full Text

Lüpkes C. & Birnbaum G. 2005. Surface drag in the Arctic marginal sea-ice zone: a comparison of different parameterization concepts. Boundary-Layer Meteorology 117, 179–211. Publisher Full Text

Paulson C.A. 1970. The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. Journal of Applied Meteorology 9, 857–861.

Persson P.O.G. 2012. Onset and end of the summer melt season over sea ice: thermal structure and surface energy perspective from SHEBA. Climate Dynamics 39, 1349–1371. Publisher Full Text

Persson P.O.G., Fairall C.W., Andreas E.L., Guest P.S. & Perovich D.K. 2002. Measurements near the Atmospheric Surface Flux Group tower at SHEBA: near-surface conditions and surface energy budget. Journal of Geophysical Research—Oceans 107, article no. 8045, doi: 10.1029/2000JC000705. Publisher Full Text

Uttal T., Curry J.A., McPhee M.G., Perovich D.K., Moritz R.E., Maslanik J.A., Guest P.S., Stern H.L., Moore J.A., Turenne R., Heiberg A., Serreze M.C., Wylie D.P., Persson O.G., Paulson C.A., Halle C., Morison J.H., Wheeler P.A., Makshtas A., Welch H., Shupe M.D., Intrieri J.M., Stamnes K., Lindsey R.W., Pinkel R., Pegau W.S., Stanton T.P. & Grenfeld T.C. 2002. Surface Heat Budget of the Arctic Ocean. Bulletin of the American Meteorological Society 83, 255–275.

Webb E.K., Pearman G.J. & Leuning R. 1980. Correction of flux measurements for density effects due to heat and water vapor transport. Quarterly Journal of the Royal Meteorological Society 106, 85–100. Publisher Full Text

Zeng X., Zhao M. & Dickinson R.E. 1998. Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using TOGA COARE and TAO data. Journal of Climate 11, 2628–2644.