The Surface Heat Budget of the Arctic Ocean (SHEBA) experiment produced 18 000 h of turbulence data from the atmospheric surface layer over sea ice while the ice camp drifted for a year in the Beaufort Gyre. Multiple sites instrumented during SHEBA suggest only two aerodynamic seasons over sea ice. In “winter” (October 1997 through 14 May 1998 and 15 September 1998 through the end of the SHEBA deployment in early October 1998), the ice was compact and snow covered, and the snow was dry enough to drift and blow. In “summer” (15 May through 14 September 1998 in this dataset), the snow melted, and melt ponds and leads appeared and covered as much as 40% of the surface with open water. This paper develops a bulk turbulent flux algorithm to explain the winter data. This algorithm predicts the surface fluxes of momentum, and sensible and latent heat from more readily measured or modeled quantities. A main result of the analysis is that the roughness length for wind speed z0 does not depend on the friction velocity u* in the drifting snow regime (u* ≥ 0.30 m s−1) but, rather, is constant in the SHEBA dataset at about 2.3 × 10−4 m. Previous analyses that found z0 to increase with u* during drifting snow may have suffered from fictitious correlation because u* also appears in z0. The present analysis mitigates this fictitious correlation by plotting measured z0 against the corresponding u* computed from the bulk flux algorithm. Such plots, created with data from six different SHEBA sites, show z0 to be independent of the bulk u* for 0.15 < u* ≤ 0.65 m s−1. This study also evaluates the roughness lengths for temperature zT and humidity zQ, incorporates new profile stratification corrections for stable stratification, addresses the singularities that often occur in iterative flux algorithms in very light winds, and includes an extensive analysis of whether atmospheric stratification affects z0, zT, and zQ.
The turbulent momentum flux from the atmosphere to compact sea ice forces the ice to move and, in turn, drives ocean currents. It also creates pressure ridges where the ice converges, opens leads where the ice diverges, and redistributes deposited snow through blowing and drifting. The turbulent surface fluxes of sensible and latent heat, in contrast, are typically secondary terms to the radiative components in the surface energy budget of sea ice (e.g., Jordan et al. 1999; Persson et al. 2002; Andreas et al. 2004a; Huwald et al. 2005). Nevertheless, because these turbulent heat fluxes couple the sea ice and the near-surface atmosphere, they dictate the thermal structure of the lower atmosphere and thereby crucially influence the momentum flux.
One of the goals of our participation in the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment was to develop a bulk flux algorithm—comparable in style and simplicity to the Coupled Ocean–Atmosphere Response Experiment (COARE) algorithm (Fairall et al. 1996, 2003)—to accurately predict the surface fluxes of momentum and sensible and latent heat (Andreas et al. 1999). Here, we describe the bulk flux algorithm for compact, snow-covered winter sea ice that we have developed. Andreas et al. (2003, 2004b) and Jordan et al. (2003) have reported our preliminary work on this algorithm. Brunke et al. (2006) offer an alternative bulk flux algorithm based on the SHEBA dataset, but we believe some of their results reflect fictitious correlation rather than real physics—an analysis problem that we address here.
Our SHEBA observations ran for almost a year—from late October 1997 through the end of September 1998. In terms of the turbulent exchange, the SHEBA year separated into two seasons based strictly on aerodynamic considerations: winter and summer (Andreas et al. 2003; cf. Brunke et al. 2006). “Winter” was the period during SHEBA when the snow was dry enough and deep enough to drift and blow. During “summer,” the snow at SHEBA was too wet to move under wind forcing and eventually melted and disappeared altogether. Winter ran from the start of the experiment, October 1997, through 14 May 1998, resumed on 15 September 1998, and continued through the end of the measurement period in late September. Summer ran from 15 May through 14 September 1998. Here, we focus strictly on the winter period, when the sea ice was compact and covered with dry snow.
We base our bulk flux algorithm on Monin–Obukhov similarity theory (e.g., Garratt 1992, 49–56). From the SHEBA dataset, Grachev et al. (2007a,b) developed new Monin–Obukhov similarity functions that specifically treat the very stable stratification that we encountered during the winter. We incorporate those into our flux algorithm.
Other essential features of a bulk flux algorithm are modules that predict the roughness lengths for the wind speed (z0), temperature (zT), and humidity (zQ) profiles. In our preliminary work on this algorithm, we had speculated that drifting and blowing snow influenced z0; it therefore increased with the friction velocity u*. Here, however, we explain that this effect was likely fictitious correlation caused by the need to use u* to compute z0. For our current analysis, we plot z0 versus a u* value derived from our bulk flux algorithm, an approach that mitigates fictitious correlation. In these plots, z0 is constant for u* values above 0.15 m s−1.
We also present plots of zT/z0 and zQ/z0 versus the roughness Reynolds number R* = u*z0/ν, where ν is the kinematic viscosity of air. These plots generally follow Andreas’s (1987) theoretical model; but, again, such plots are prone to fictitious correlation. We circumvent that problem by using the flux data to compute zT/z0 and zQ/z0 but using our bulk flux algorithm to compute R*. We also test z0, zT/z0, and zQ/z0 for the influence of stratification. Both z0 and zQ/z0 seem independent of stratification; the behavior of zT/z0 is less clear.
2. Bulk flux algorithm
Energy budget studies or atmospheric models with snow as the lower boundary almost always estimate the surface fluxes of momentum τ, sensible heat Hs, and latent heat HL from a bulk flux algorithm (e.g., Brun et al. 1989; Jordan et al. 1999; Bintanja 2000; Lehning et al. 2002; Briegleb et al. 2004). In our algorithm, the relevant flux equations take the form
In these equations, u, w, θ, and q are turbulent fluctuations in longitudinal wind speed, vertical wind speed, temperature, and specific humidity, respectively; the overbar indicates a time average. Also, ρ is the air density; cp is the specific heat of air at constant pressure; Lυ is the latent heat of sublimation; Sr is the effective wind speed at reference height r; Θr and Qr are the potential temperature and specific humidity at r, respectively; and Θs and Qs are the temperature and specific humidity at the snow surface. We evaluate Qs as the saturation value at Θs. Equation (2.1a) also defines the friction velocity u*, which we use henceforth as a surrogate for surface stress.
The crux of any bulk flux algorithm is evaluating the transfer coefficients for momentum, sensible heat, and latent heat appropriate for height r—that is, CDr, CHr, and CEr, respectively, in (2.1). These generally derive from Monin–Obukhov similarity theory and formally are (e.g., Garratt 1992, 52–55; Andreas 1998)
In these equations, k (=0.40) is the von Kármán constant, and ψm and ψh are empirical functions of the Obukhov length
Here, g is the acceleration of gravity; Θ, Θυ, and Q are surface-layer averages of the air temperature, virtual temperature, and specific humidity, respectively; and is the flux of virtual temperature.
For the stratification corrections ψm and ψh in (2.2), we use the functions from Paulson (1970) in unstable (Un) stratification and the functions from Grachev et al. (2007a) in stable (St) stratification. These latter functions are based on our SHEBA dataset and include proper treatment of a heretofore unrecognized scaling regime in very stable stratification.
The z0, zT, and zQ in (2.2) are the roughness lengths for wind speed, temperature, and humidity, respectively. Developing a new parameterization for z0 and testing Andreas’s (1987) theoretical model for zT and zQ are the main subjects of this paper.
Last, Sr in (2.1) is an effective mean wind speed. For compatibility with the COARE algorithm (Fairall et al. 1996, 2003) and other recent flux algorithms (e.g., Andreas et al. 2008), we acknowledge that in unstable stratification gustiness enhances turbulent exchange and therefore model
Here, Ur is interpreted as the measured or modeled magnitude of the mean wind vector at reference height r, βg = 1.25 (Fairall et al. 1996), and w* is Deardorff’s (1970) convective velocity scale (Godfrey and Beljaars 1991):
where zi is the depth of the convective boundary layer. We take zi as a constant, 600 m (Kahl 1990; Serreze et al. 1992; Bradley et al. 1993; Tjernström and Graversen 2009), because variability in it does not have much effect on our calculations.
We adopt the suggestion by Jordan et al. (1999) that a similar “windless” coefficient is necessary for stable stratification but express it as
Here, both Ur and Sr are in meters per second. Equation (2.6) has some similarity with Mahrt’s (2008) recent parameterization that includes a term to quantify mesoscale meandering flow where we have the sech term.
In effect, both (2.4) and (2.6) prevent a singularity in the bulk flux algorithm by making the transfer coefficients well behaved in light winds (cf. Godfrey and Beljaars 1991; Fairall et al. 1996; Zeng et al. 1998; Andreas et al. 2008). If we use Ur instead of Sr, then the transfer coefficients must approach infinity to maintain finite fluxes as Ur approaches zero. But as Ur increases from zero, the Sr values calculated with either (2.4) or (2.6) quickly approach Ur for Arctic conditions. For instance, when Ur reaches 2.2 m s−1, Sr is already within 5% of Ur.
Because (2.1) and (2.2) are coupled through the Obukhov length (2.3), they must be solved iteratively using the mean measured or modeled conditions, namely, Ur, Θr, Qr and Θs. Moreover, (2.4) and our algorithms for z0, zT, and zQ also include u*, as we will show. These equations must also be part of the iteration. That iteration usually converges in 3–5 steps.
3. The SHEBA data
The SHEBA ice camp drifted approximately 2700 km in the Beaufort Gyre between 2 October 1997 and 11 October 1998 (Uttal et al. 2002).
During our SHEBA deployment, we had one central site in the SHEBA ice camp and, usually, four remote sites that ranged in distance from 0.4 to 20 km from the main camp. We serviced these remote sites about once a week. Andreas et al. (1999, 2002, 2006), Persson et al. (2002), Grachev et al. (2005, 2007a), and Brunke et al. (2006) describe the instruments that our Atmospheric Surface Flux Group (ASFG) deployed during SHEBA and review our data processing. Persson et al. show pictures of the instruments at our main site. Most of these papers, however, discuss data only from our main site. So far, the paper by Brunke et al. (2006) is the only one to make extensive use of the turbulence data from our remote sites. Here, we concentrate equally on data from our main and remote sites.
The centerpiece of our site in the main SHEBA camp was a 20-m tower instrumented at five levels with identical sonic anemometer/thermometers [K-type sonics from Applied Technologies, Inc. (ATI)] and Vaisala HMP235 temperature and humidity sensors. The tower also held one Ophir hygrometer that was mounted at 8 m, near the sonic at that level.
Through eddy-covariance measurements using standard turbulence processing, as described in Persson et al. (2002), Grachev et al. (2005, 2007a), and Andreas et al. (2006), we measured τ and Hs at each of the five tower levels and HL at one level [Eq. (2.1)]. This latter was the only direct, long-term measurement of latent heat flux from SHEBA.
About 30 m from the ASFG tower was a full suite of radiometers for measuring incoming and outgoing longwave and shortwave radiation, along with several additional sensors for measuring Θs (e.g., Claffey et al. 1999). Generally, for Θs we used the value implied by the emitted () and incoming () longwave radiation measured by our Eppley pyrgeometers (cf. Fairall et al. 1998):
Here, ɛ (=0.99; Warren 1982; Dozier and Warren 1982; Jordan et al. 1999) is the surface emissivity and σ (=5.670 51 × 10−8 W m−2 K−4) is the Stefan–Boltzmann constant. See the SHEBA data archive (http://www.eol.ucar.edu/projects/sheba) for tabulations and descriptions of these surface temperature data and the other datasets that we use in this study.
Our remote sites were instrumented with Flux-PAM (portable automated mesonet) stations from the National Center for Atmospheric Research instrument pool (Militzer et al. 1995; Horst et al. 1997). Our first four sites, which we deployed in October 1997, were named Atlanta, Baltimore, Cleveland, and Florida. The Cleveland station was damaged by a pressure ridge in early February 1998, removed from service, and redeployed at a new site called Seattle in mid-April 1998. Seattle, however, became an untenable site because of ice motions, and this PAM station was again redeployed to a site named Maui in mid-June 1998. That Maui site lasted until late September 1998, as did the original Atlanta, Baltimore, and Florida sites.
We have found the data from Seattle to be disturbed by a pressure ridge just upwind of the station and, thus, do not include data from that site in our analysis. We do use data from the other five PAM sites, however.
Each Flux-PAM station measured at one height the same quantities that we measured at the ASFG tower—that is, wind speed and direction, temperature, relative humidity, and the turbulent fluxes of momentum and sensible heat. The Web site (http://www.eol.ucar.edu/rtf/projects/SHEBA) contains instrument details, a history of each of the PAM sites deployed during SHEBA, and information on data processing.
Briefly, the PAM stations provided hourly averaged data, as did the main tower site. Each PAM station used a sonic anemometer/thermometer mounted at a height between 2.5 and 3.5 m, depending on instrument type, to measure the turbulent fluxes of momentum and sensible heat by eddy-covariance techniques [i.e., τ and Hs in (2.1)]. We used sonics from both Gill (Solent R2) and Applied Technologies, Inc. (K-type), for the turbulence measurements at these sites.
Sonic anemometer/thermometers do not measure the turbulent air temperature θ exactly. Rather, the measured turbulent “sonic” temperature θs is a combination of θ and q, the turbulent fluctuation in specific humidity (Schotanus et al. 1983; Kaimal and Gaynor 1991; Larsen et al. 1993):
This θs is very close to the fluctuation in virtual temperature
Consequently, the covariance between θs and the vertical velocity (w) is
This covariance is almost identical to the flux of virtual temperature required in the Obukhov length, (2.3), and it can be used there directly with negligible error (e.g., Kaimal and Finnigan 1994, p. 224). Furthermore, because of the low humidity and small Bowen ratio over winter sea ice, wθ and are typically within 5% of each other and are essentially interchangeable in our dataset (cf. Andreas et al. 2005; Grachev et al. 2005).
The Flux-PAM sites used Vaisala HMD50Y sensors in mechanically aspirated radiation shields to measure mean temperature and relative humidity (Andreas et al. 2002). Each PAM station also measured barometric pressure with a Vaisala PTB 220B digital barometer. Because our main ASFG site had no measurement of barometric pressure, we also used the pressure data from the Florida PAM site at the tower site because only about 400 m separated the two sites.
Each PAM station included sets of up-looking and down-looking radiometers to measure shortwave (Kipp and Zonen model CM21) and longwave (Eppley model PIR) radiation. We obtained surface temperature at the PAM sites exclusively from these longwave radiometers through (3.1).
We saw early in our SHEBA deployment that sensor riming was a problem at the PAM sites. Rime ice occasionally coated the domes of the radiometers (especially the up-looking ones) and thereby ruined the radiation measurements and our estimates of surface temperature. Rime also collected on the sonics; when it coated the transducers, we lost wind information. In early March 1998, we installed effective heaters and blowers on all four radiometers at each site. This heating minimized the effects of dome icing for the rest of the experiment.
To mitigate riming on the PAM sonics, in mid-January 1998 we fixed heating tape around each transducer and at several locations along the support arms. The computer-controlled data acquisition system that monitored the PAM stations turned this heating on, however, only when it noticed that data returns from a sonic were deteriorating; this heating remained on only until the data returns again reached 100%. We flagged these heating periods and did not use any flux data that were collected when the heaters were on.
The sonics and radiometers at our main site did not suffer as much from riming as the PAM instruments. First, someone was always attending the instruments at our main site and cleaning them as necessary. Second, the radiometers at our main site had effective blowers and heaters from the beginning of the experiment. Third, our 20-m tower was fitted exclusively with ATI sonics, which seemed to resist riming much better than the Gill sonics that were originally mounted on each PAM station. We later replaced the Gill sonics with ATI sonics at the three most remote PAM sites.
From the turbulent fluxes and mean meteorological quantities measured at multiple levels on our main tower and at the Flux-PAM sites, we could compute the turbulent transfer coefficients CDr, CHr, and CEr from (2.1). These, in turn, give us estimates of the roughness lengths from (2.2):
Here, the roughness lengths are in meters if r is in meters. Moreover, solving these involves no iteration because our data also provided L.
4. Uncertainty analysis and quality controls
Turbulence data generally have a lot of random scatter. Roughness lengths evaluated from these data, in turn, are commonly quite scattered because they derive from (3.5) and (2.1), and thus rely on several mean and turbulence variables. The only way to overcome this scatter is to collect a lot of high-quality data.
Table 1 summarizes typical uncertainties in the quantities that we measured at the ASFG tower and PAM sites during SHEBA and in the variables that we calculate from these data. We base these estimates on previous similar summaries by Fairall et al. (1996) and Persson et al. (2002), on data analyses by Andreas et al. (2002), and on similar uncertainty analyses by Andreas et al. (2005, 2006). The essential messages of Table 1 are that evaluations of z0 are uncertain by a factor of 3 and evaluations of zT and zQ are uncertain by a factor of 200.
In light of these uncertainties, we screened the data that we used in our analyses to ensure a reasonable signal-to-noise ratio. If any hour of data from a PAM station or from any level on the ASFG tower met the following criteria, we excluded the data as inadequate:
Only data that fail (4.1a) or (4.1b) prevent our computing z0. Remember, we need Hs—actually ρcp—to compute the Obukhov length. Data that fail (4.1a), (4.1b), or (4.1d) prevent our computing zT. A failure to pass any of (4.1a), (4.1b), (4.1c), or (4.1e) prevent our computing zQ. These data requirements explain the number of available z0, zT, and zQ values that we summarize in Table 2.
As additional screening, if any calculated roughness lengths met the following criteria, we assumed the result was unrealistic and ignored it:
We instituted (4.2a) because roughness lengths simply cannot be this large over snow-covered, compact sea ice (e.g., Banke et al. 1980; Overland 1985; Guest and Davidson 1991; Andreas 1995). The second limit, (4.2b), is the approximate mean free path of air molecules at sea level. We presume that the surface exchange of heat and moisture cannot occur at scales smaller than this (cf. Andreas and Emanuel 2001; Andreas et al. 2008).
In effect, the screening criteria (4.1) and (4.2) excluded from our analysis situations in which Monin–Obukhov similarity theory breaks down. In very stable boundary layers, a host of phenomena occur that violate similarity theory; for example, the boundary layer may be so thin that a constant-flux layer does not exist, the turbulence may be only intermittent, gravity waves can confound the turbulence series, or a low-level jet rather than the surface may be the source of the turbulence (e.g., Mahrt 1998, 1999; Grachev et al. 2005). Tests (4.1) and (4.2) tended to keep us out of these regimes.
Because the five ASFG tower levels each had a sonic anemometer/thermometer and a temperature and humidity sensor, any hour of data could yield from none to five independent estimates of z0 and zT from this site. We did use the same Θs for all estimates, though. We do not report all of these values. Rather, we took the median value for all results that passed our screening for that hour. The median is the “most common robust and resistant measure of central tendency” (Wilks 2006, p. 26). For example, later we will show plots of measured z0 versus a bulk estimate of u*. For the tower data, the plotted z0 value will be the median from all tower levels reporting z0 for that hour. We also ran our bulk flux algorithm for all tower levels with sufficient data for it; the bulk u* is, again, the median from all tower levels that produced a bulk u*. Hence, some of our hourly tower estimates are based on data from only one level, but some are based on data from all five levels.
Because of this “averaging,” we tend to have more confidence in the results from the ASFG tower than from the Flux-PAM stations, which did not have the luxury of redundant measurements. Taking the median value of the tower data also tends to mitigate the effects of fictitious correlation because, for example, the tower levels that yielded median measured values of zT and z0 did not always yield the median bulk estimates of u* and z0 that we used in comparing measured z0 against bulk u* and measured zT/z0 against bulk R* = u*z0/ν in plots to follow.
5. Roughness length for wind, z0
Large and Pond (1982) reported that, over the open ocean, zT depends on the atmospheric stratification. In the next section, we will test whether our zT and zQ data depend on stratification. As a prelude to those tests, we look here at the related question of whether z0 depends on atmospheric stratification.
Figure 1 therefore shows the hourly z0 data collected at all six SHEBA winter sites and sorted according to stratification—determined from the measured Obukhov length. The horizontal axis is the u* value computed with our bulk flux algorithm for the hour plotted.
In preliminary analyses (Andreas et al. 2003, 2004b), we had plotted measured z0 against measured u*. If u* is measured to be erroneously large, then propagating this error through (3.5a) convinced us that z0 would be evaluated to be erroneously large. Likewise, if u* is measured to be erroneously small, then z0 will be evaluated to be erroneously small. As a result, plots of z0 versus measured u* tend to have a positive slope because of shared errors. We call this fictitious correlation because it does not result from real physics. This fictitious correlation probably explains why Andreas et al. (2003, 2004b) and Brunke et al. (2006) found z0 to generally increase with u* in earlier analyses of the SHEBA dataset.
By plotting measured z0 versus an estimate of u* from a bulk flux algorithm in Fig. 1, we eliminate the shared errors and anticipate minimal effects from fictitious correlation. Moreover, the objective of a bulk flux algorithm is to estimate an accurate roughness length from a bulk flux estimate of u*—that relationship is what Fig. 1 shows.
The black circles in Fig. 1 are geometric means of the hourly z0 values in u* bins that are mostly 2 cm s−1 wide. Measured roughness lengths have distributions that are approximately lognormal (cf. Vickers and Mahrt 2006). The geometric mean, which is the exponential of the average of the ln(z0) values within a bin, is a better indicator of the central tendency of this distribution than is the arithmetic mean, which preferentially weights the largest z0 values. The error bars in Fig. 1 are likewise computed as two standard deviations in the means of these ln(z0) values.
For reference in Fig. 1, the curve is
Here, ν is again the kinematic viscosity of air [see Andreas (2005) for functions to compute this and other constants that we use here] and u*,B is the friction velocity from our bulk flux algorithm. In (5.1), z0 is in meters when the other variables have mks units.
The two panels in Fig. 1 display some obvious differences but also many similarities. To evaluate whether the z0 values depicted in the two panels differ and thus imply a stratification dependence, we make statistical tests bin by bin. Table 3 lists the geometric mean z0 values from Fig. 1 for the stable and unstable bins that the two panels in the figure have in common.
We test the difference between these bin averages using Student’s t statistic, which we calculate as
Here, xs and xu are the bin averages of the ln(z0) values collected in stable and unstable stratification, respectively; σs and σu are the corresponding standard deviations; and Ns and Nu are the numbers of samples in the stable and unstable bins. The t statistic has Ns + Nu − 2 degrees of freedom (DOF).
The t values in Table 3 are almost evenly spread between positive and negative—that is, for some bins, the mean z0 in the stable bin is larger than the mean z0 in the unstable bin; however, in other bins, the opposite is true. When |t| is greater than 1.64, we can say with 95% statistical confidence that the stable and unstable bins have different means. When |t| is smaller, statistical differences are less certain. Of the 22 bins in Table 3, 12 meet the criterion of |t| > 1.64. But for 7 of these 12 bins, t is positive (stable mean larger than unstable mean); whereas for 5 bins, t is negative (unstable mean larger than stable mean). Consequently, these tests yield no compelling evidence that z0 depends on stratification. In our subsequent analysis and in our bulk flux algorithm, we presume that z0 does not depend on stratification.
For Fig. 2, we therefore bin-averaged the z0 data collected in both stable and unstable stratification at the five SHEBA sites with the longest records. As in Fig. 1, the bin-averaged z0 values here are geometric means. Table 2 summarizes the amount of data that went into Fig. 2.
We fitted (5.1) by eye to the data in Fig. 2. The first term on the right of (5.1) represents aerodynamically smooth scaling. Here, z0 decreases with increasing u*,B. The data in Fig. 1 collected in stable stratification and three of the five sites depicted in Fig. 2 follow aerodynamically smooth scaling for u*,B between 0 and 0.15 m s−1. The scanty data from Cleveland (i.e., large error bars) and the Florida data do not show the minimum in z0 associated with aerodynamically smooth scaling, however. The Cleveland site was chosen because it was in very rough ice; Florida was also visually rougher than our other sites. Perhaps, these two sites rarely experienced aerodynamically smooth flow.
Beyond about u*,B = 0.04 m s−1, z0 increases with u*,B; however, that increase is effectively over when u*,B reaches approximately 0.15 m s−1, and z0 becomes constant at about 2.3 × 10−4 m. All five sites represented in Fig. 2 support the conclusion that z0 is independent of u*,B for u*,B between 0.15 and 0.65 m s−1, the upper limit of our data. Our preliminary analyses had shown z0 to continue increasing when we plotted it against the measured value of u*. We henceforth adopt (5.1) as the z0 parameterization in our bulk flux algorithm.
Although some theoretical arguments predict that z0 should increase with u* when snow begins drifting (e.g., Owen 1964), other theories imply that z0 should reach a plateau and even decrease with increasing u* because more and more particles are suspended in a turbulent flow. This latter argument speculates that the suspended particles increase the density of the near-surface fluid enough to decouple the surface from higher-level winds, thereby reducing the vertical momentum exchange. Wamser and Lykossov (1995) applied these ideas to drifting snow, and Makin (2005) and Barenblatt et al. (2005) considered these stabilizing effects as a means by which sea spray could cause the drag coefficient over the ocean to level off in high winds.
By obviating the misleading effects of fictitious correlation with our analysis, we argue against the explanation that drifting snow explains the increase in z0 with u* that is commonly reported (e.g., Chamberlain 1983; Pomeroy and Gray 1990; Pomeroy et al. 1993; Bintanja and Van den Broeke 1995). But we cannot say from our data whether the combined stabilizing effects of airborne snow and the mechanical increase in z0 associated with saltating snow explain why z0 is apparently constant in the drifting snow regime in our dataset.
6. Scalar roughness lengths for zT and zQ
Our candidate expression for the roughness lengths zT and zQ is Andreas’s (1987) model. This expresses the ratio zs/z0, where zs is either scalar roughness (zT or zQ) as a function of the roughness Reynolds number R* = u*z0/ν:
Andreas (1987, 2002) tabulates the polynomial coefficients b0, b1, and b2. Many sources corroborate that the scalar roughness should depend on the roughness Reynolds number (e.g., Garratt and Hicks 1973; Liu et al. 1979; Brutsaert 1982, 89–97) because R* quantifies how far the physical roughness elements protrude above the molecular sublayer (Andreas 1998).
Andreas (2002) tested (6.1) with data from the literature. None of these data were collected over sea ice but, rather, came from snow-covered ground or glaciers. Denby and Snellen (2002) likewise tested the zT prediction from (6.1) over an Icelandic glacier and found good agreement. Reijmer et al. (2004) tested (6.1) in a regional climate model to predict the heat fluxes over the Antarctic continent and concluded it was a “good option.”
Evidently, only Andreas et al. (2005) and Brunke et al. (2006), however, have tested (6.1) with flux data collected over sea ice. Andreas et al. used the data from Ice Station Weddell (ISW), and Brunke et al. used these same SHEBA data. Both teams found that (6.1) fit their data reasonably well. Here, we redo some of the analyses in Brunke et al. for several reasons: They did not study the behavior of zQ; we introduce new techniques to mitigate the effects of fictitious correlation; and we consider possible stratification effects on zT and zQ.
Andreas (2002) demonstrated that, when all quantities are measured, plots of zT/z0 and zQ/z0 versus R*, necessary to test (6.1), suffer from fictitious correlation such that the plotted data naturally tend to follow the slope that (6.1) predicts. If we use a bulk flux algorithm to compute R*, however, the same measured u* and z0 do not appear in both dependent and independent variables, and the fictitious correlation is mitigated. Andreas et al. (2006) first suggested using a bulk estimate of R*, denoted R*,B = u*,Bz0,B/ν, where u*,B and z0,B come from our bulk flux algorithm, to parameterize another set of turbulence measurements. Here, we exclusively use R*,B to study the behaviors of zT/z0 and zQ/z0. In fact, this is the proper way to develop a bulk flux algorithm because, in practice, modelers would have only the bulk estimate of R* but need to predict accurate values of zT/z0 and zQ/z0.
Large and Pond (1982; also Large et al. 1994) suggested that zT over the open ocean depends on the atmospheric stratification; they found zT to be four orders of magnitude larger in unstable stratification than in stable stratification. To our knowledge, no other studies of oceanic zT have corroborated this result (e.g., Fairall et al. 1996, 2003; Zeng et al. 1998; Brunke et al. 2003). Moreover, we know of no studies that have even considered the possible stratification dependence of z0, zT, and zQ over snow or sea ice. We have enough data from SHEBA to consider this question.
Figures 3 and 4 introduce our analysis of stratification effects on zT. Figure 3 shows all the hourly zT/z0 data from the ASFG tower, separated into stable and unstable cases. Figure 4 is a similar plot of data collected at Florida. We present these two plots to assess whether the one-level data from the Flux-PAM stations are significantly different from the median of the multilevel data from the ASFG tower. In these plots, zT and z0 were measured, while the independent variable is R*,B, an estimate of the roughness Reynolds number from our bulk flux algorithm. In both plots, the data collected in stable stratification tend to straddle the curve representing Andreas’s (1987) model, (6.1). Both plots, however, suggest that zT/z0 values from unstable stratification tend to be above that curve.
As Table 2 shows, we collected far fewer zT data during unstable stratification than during stable stratification. Because of the uncertainties associated with measuring zT and z0, as tabulated in Table 1, we need to bring as much data as possible to bear on this question before reaching a conclusion. Figure 5 is therefore like Fig. 1: It shows all our winter zT/z0 data in separate plots for stable and unstable stratification. The figure also shows bin averages and error bars. The bin averages and error bars are, again, based on the geometric mean.
The bin averages in the stable panel in Fig. 5 agree very well with Andreas’s (1987) theoretical model, (6.1), except for small R*,B, where the two deviant bins contain few data points. Some of the points in the unstable panel also agree with Andreas’s model, but other points are significantly above it. In Table 4, we use Student’s t statistic to test for statistical differences between the bin-averaged values of zT/z0 in stable and unstable stratification in Fig. 5, as we did with z0 in Table 3. As in Table 3, the means and standard deviations used in (5.2) are based on natural logarithms—of zT/z0 in this case.
All t values in Table 4 are negative: The bin-averaged zT/z0 values tend to be smaller in stable stratification than in unstable stratification. For the bins in Table 4 for which t < −2.6, we can reject the hypothesis that the bin-averaged zT/z0 values in stable and unstable stratification are the same with 99% statistical confidence. For the t value of −2.10, we can reject this null hypothesis with 97.5% confidence. But for the three bins with smaller |t|, the difference between stable and unstable stratification is less statistically clear.
Figure 6 shows the winter zQ data from SHEBA. As in Figs. 3 and 4, we distinguish data collected in stable and unstable stratification. Contrary to these figures, though, most of the zQ data came from unstable stratification. As Fig. 6 and Table 2 show, few latent heat flux values that were measured in stable stratification survived our screening. When the sensible heat flux is downward (i.e., in stable stratification), the absolute value of the latent heat flux is often so small that screening condition (4.1c) or (4.1e) told us to ignore that data.
The stable and unstable cases in Fig. 6 are fairly evenly distributed on either side of Andreas’s (1987) model. Hence, we discount the idea that stratification influences zQ and average both stratification classes together to get the bin-averaged data in the figure. These bin averages agree adequately with Andreas’s model.
Because Fig. 6 suggests that stratification does not significantly affect zQ and because of the inconclusive results concerning stratification effects on zT, we retain (6.1) for predicting both zT/z0 and zQ/z0 in our bulk flux algorithm. But because of the possible stratification effects hinted by Fig. 5 and Table 4, we encourage others who measure z0, zT, and zQ to analyze these for stratification effects. We can conclude, nevertheless, that our data do not support Large and Pond’s (1982) result that zT is four orders of magnitude larger in unstable stratification than in stable stratification.
Although Andreas (1987) shows predictions for zT/z0 and zQ/z0 for R* values up to 1000, Figs. 3 –6 actually span the normal R*,B range that exists over snow-covered sea ice. The maximum hourly averaged wind speed that we measured at the 3-m level during winter on the ASFG tower was about 15 m s−1 (Andreas et al. 2002). Wind speeds at higher levels would have been higher, but we can presume that near-surface winds over Arctic sea ice during winter will almost always be less than 20 m s−1. We therefore estimate that u*,B will rarely be greater than 0.85 m s−1 and that R*,B will rarely be as large as 20.
7. Testing the algorithm
Equations (2.1)–(2.6), (5.1), and (6.1) are the main equations in our bulk flux algorithm. Andreas (2005) gives the additional equations that we use to compute other dynamic and thermodynamic variables in the algorithm, such as kinematic viscosity, specific humidity, and latent heat of sublimation.
To test the effectiveness of this algorithm, we compare it against the Community Ice Code (CICE; Hunke and Lipscomb 2008). The CICE module is an integral part of the Community Climate System Model (Briegleb et al. 2004) and, as such, is a crucial tool for climate studies (e.g., Kiehl and Gent 2004; Collins et al. 2006).
Although CICE is much more than a bulk turbulent flux algorithm, surface flux calculations are essential in it; it parameterizes these fluxes using the same formalism that our algorithm does. CICE differs in detail from our algorithm, however. It assumes z0 = zT = zQ = 5.0 × 10−4 m. Although it uses the Paulson (1970) functions for ψm and ψh in unstable stratification, as does our algorithm, it uses the functions from Holtslag and De Bruin (1988) in stable stratification. Remember, the SHEBA algorithm uses the new functions from Grachev et al. (2007a). Finally, CICE does not have a gustiness parameterization, like our (2.4), in unstable stratification but does parameterize windless transfer in stable stratification—for sensible heat only—following Jordan et al. (1999) [compare our (2.6)]. In effect, for stable stratification, CICE adds to all computed sensible heat fluxes 1 W m−2 K−1 (Θs − Θr). Consequently, all negative fluxes become more negative.
Figure 7 compares u* and Hs values measured at the Flux-PAM site called Atlanta with values modeled with both our new SHEBA algorithm and the surface flux algorithm in CICE. We chose Atlanta for this comparison because it has a long record but, unlike the ASFG tower, had only one measurement height and, therefore, provides a data source similar to what a model would be replicating. In our plots of Hs and HL, a positive flux is upward—from surface to air [Eqs. (2.1b) and (2.1c)].
To evaluate the performance of the two flux algorithms, we computed two metrics, following Willmott (1982). Let Mi be a measured flux (i.e., u*, Hs, or HL) and let Bi be the corresponding estimate from the SHEBA or CICE bulk flux algorithm. The mean bias error in the modeled (bulk) values is
where N is the number of observations (see Table 2). The root-mean-square error (RMSE) in the modeled values is
In Fig. 7, the fitting lines for the SHEBA algorithm are visually closer to the 1:1 line than are the fits based on the CICE algorithm. Table 5 corroborates these visual results: both the bias errors and the root-mean-square errors are smaller with the SHEBA algorithm than with the CICE algorithm.
You may complain that Fig. 7 is not a meaningful comparison because it is based on the same data that we used to develop the SHEBA algorithm. For the u* comparison, this is a valid concern. Table 2 shows that we used 2263 Atlanta measurements to develop our z0 parameterization and 2340 u* measurements to test it in Fig. 7. However, we used only 401 Atlanta zT/z0 values to test our scalar roughness model in Fig. 5 but 2553 measurements of Hs from Atlanta to test our algorithm in Fig. 7. Furthermore, we did no tuning of the zT/z0 algorithm based on any SHEBA data but, rather, used the algorithm as originally published in Andreas (1987).
Nevertheless, to make independent tests of the SHEBA and CICE algorithms, we repeat the comparisons in Fig. 7 using data from ISW. Ice Station Weddell was a 4-month deployment over sea ice in the western Weddell Sea in 1992 in winter conditions (e.g., Andreas and Claffey 1995; Andreas et al. 2004a, 2005). The dataset includes mean meteorological quantities and hourly eddy-covariance measurements of u*, Hs, and HL with an ATI sonic anemometer/thermometer and a Lyman-alpha hygrometer (from Atmospheric Instrumentation Research) at one height, 4.65 m. The ISW site was on a large second-year ice floe that had a snow cover of about 50 cm. In essence, ISW provided more than three months of data in conditions like those representing Atlanta in Fig. 7.
Figures 8 –10 show scatterplots of measured and modeled u*, Hs, and HL values, respectively, based on the ISW data. Again, the modeled values come from our SHEBA bulk flux algorithm and from the CICE algorithm. Table 5 lists the metrics for the comparisons.
The CICE algorithm does better in representing u* for the ISW dataset than does the SHEBA algorithm. Andreas et al. (2005) explained that ISW seemed to be an aerodynamically rougher site than SHEBA. The larger, constant value of z0 = 5.0 × 10−4 m in the CICE algorithm does better, on average, in representing this surface than does our SHEBA algorithm, which predicts a maximum z0 of 2.3 × 10−4 m.
Because of its constant z0, the CICE algorithm is poorer than the SHEBA algorithm at predicting u* for small u*. For small u*, the SHEBA algorithm allows z0 to follow aerodynamically smooth scaling. The small z0 values that result produce better agreement between measured and modeled u* with the SHEBA algorithm than with the CICE algorithm for u* < 0.1 m s−1.
Because the CICE algorithm assumes zT = zQ = z0 = 5.0 × 10−4 m, however, the SHEBA algorithm does better than it in representing the sensible and latent heat flux data in Figs. 9 and 10 (also see Table 5). Values of zT and zQ simply tend to be smaller than 5.0 × 10−4 m over sea ice (Andreas et al. 2004a,b). Because of these large scalar roughness lengths in the CICE algorithm, both Hs and HL tend to be larger in magnitude than they should be in both stable and unstable stratification. In other words, negative fluxes are more negative than they should be and positive fluxes are more positive. The fitting line is thus rotated counterclockwise from the 1:1 line, as we see in the CICE panels in Figs. 9 and 10.
The metrics for the ISW heat flux comparisons in Table 5 confirm the better results with the SHEBA algorithm. For sensible heat flux, the CICE algorithm has a bias error that is 10 W m−2 more negative than with the SHEBA algorithm and a root-mean-square error that is 50% larger. Similarly, for the latent heat flux, the bias error for the CICE algorithm is 50% larger than for the SHEBA algorithm, and the root-mean-square error is 20% larger.
Figures 7, 9, and 10 emphasize how small the sensible and latent heat fluxes are over winter sea ice and, thus, how difficult they are to both measure and model. The sensible heat flux is generally measured to be between −20 and +20 W m−2, and the latent heat flux is typically between −5 and +5 W m−2.
The scatter in Figs. 7 –10 is probably typical of the uncertainty in surface fluxes when only one set of sensors is available for modeling or measuring these fluxes. Therefore, the values in these panels might be examples of comparisons that would result from large-scale model predictions of fluxes that are based on only the surface temperature and the wind speed, air temperature, and humidity from the lowest model level. And we cannot tell whether the scatter results from random measurement errors in the fluxes or uncertainties in the variables used in the bulk flux algorithm. Most likely, uncertainties in both measured and modeled quantities contribute approximately equally to the scatter.
Under this assumption that half of the scatter results from the measurements and half from the model, we can use the root-mean-square errors in Table 5 to estimate the typical accuracy in flux estimates derived from our flux algorithm. For u*, the entries in that table imply an accuracy of 0.02–0.03 m s−1; for Hs, 6–8 W m−2; and for HL, 2 W m−2.
The u* comparisons in Figs. 7 and 8 have further implications for creating a bulk flux algorithm for sea ice. The CICE algorithm, which uses a constant and relatively large value for z0, overestimates u* for the Arctic data (Fig. 7) but does well for the Antarctic data (Fig. 8). The SHEBA algorithm, for which z0 is relatively smaller, does very well in representing the Arctic u* data, but it underestimates u* for the Antarctic data. It may ultimately be necessary to introduce a parameterization for z0 that is site specific or depends on the topography of the sea ice and the characteristics of the snow cover.
Besides developing (5.1) and confirming (6.1), we calculated average values for z0, zT, and zQ based on all the winter SHEBA data. These values—z0 = 2.1 × 10−4 m, zT = 2.0 × 10−4 m, and zQ = 3 × 10−4 m—can be substituted for (5.1) and (6.1) in our algorithm and therefore would constitute an alternative flux algorithm that requires fewer computations than our full algorithm.
Although this algorithm might be appropriate for some applications, we do not advocate it in general. First, these average values for the roughness lengths may be appropriate only for the SHEBA dataset. Meanwhile, our full formulations for z0, zT, and zQ are grounded in theory and, therefore, can be extrapolated to conditions not present in the SHEBA dataset. Second, because a constant z0 of 2.1 × 10−4 is too large when u*,B is small but is too small when u*,B is large (see Fig. 2), it leads to u*,B estimates that are biased high compared with data when u*,B is small and to estimates that are biased low when u*,B is large.
Lastly, for selected SHEBA sites, we calculated the mean bias errors and the root-mean-square errors for u*, Hs, and HL measurements and corresponding values computed with this version of the algorithm that uses constant roughness lengths (i.e., as in Table 5). We do not tabulate these error metrics but simply report that they generally indicate slightly poorer agreement between the data and this version of the algorithm than for our full algorithm.
Using data from the year-long experiment that studied the Surface Heat Budget of the Arctic Ocean, we have developed and tested parameterizations for the turbulent fluxes of momentum and sensible and latent heat for winter conditions over sea ice. Ours is an aerodynamic definition of winter as the period during which the sea ice is compact and snow covered, and the snow is dry enough to drift and blow in response to the wind.
The bulk flux algorithm that we developed from the SHEBA turbulence data consists of the usual flux equations, (2.1) and (2.2); an equation for the Obukhov length, (2.3); and expressions for the effective wind speed Sr, (2.4)–(2.6). We also include the new corrections for stable stratification, ψm and ψh, that Grachev et al. (2007a) deduced from the SHEBA data. The other major components of our algorithm are a new parameterization for z0, (5.1), and Andreas’s (1987) model for zs/z0, (6.1), where zs is the scalar roughness, either zT or zQ. All these equations are coupled and must be solved iteratively for the three fluxes; the solution typically converges within 3–5 iterations.
Fictitious correlation is a common but often unrecognized problem in attempts to validate parameterizations for z0 as a function of u* and zs/z0 as a function of the roughness Reynolds number R* = u*z0/ν. As a way to mitigate the fictitious correlation, Andreas et al. (2006) suggested using an R* value based on results from a bulk flux algorithm rather than on forming R* from measurements of u* and z0. We implement that practice here in our plots of zT/z0 and zQ/z0 versus R*.
We also realized that plots of measured ln(z0) versus measured u* almost always produce a positive slope because of the fictitious correlation. That is, z0 appears to increase with u*—long believed to be a signature for the influence of drifting snow. Hence, here we plot measured ln(z0) versus an estimate of u* from our bulk flux algorithm to mitigate this fictitious correlation.
Such plots, based on more than 9000 h of data from six SHEBA sites, show no z0 dependence on the bulk u* in the presumed blowing and drifting snow regime—u* values of 0.30 m s−1 and above (Andreas et al. 2004b). This analysis thus yielded our new parameterization for z0 in terms of u*,B, (5.1). In this relation, z0 does still depend on u*,B but only for u*,B below 0.15 m s−1, where three of five sites suggest that z0 obeys aerodynamically smooth scaling.
We investigated the possibility that z0, zT/z0, and zQ/z0 depend on atmospheric stratification. Although we see no theoretical reason why such stratification dependence should exist (i.e., Andreas 1987), one dataset collected by Large and Pond (1982) over the open ocean suggests that zT depends strongly on atmospheric stratification. To our knowledge, no one has considered this question for data collected over snow or ice. Neither the z0 nor zQ/z0 data, however, provide any compelling evidence that these quantities depend on atmospheric stratification. The results for zT/z0 are not as conclusive; therefore, this ratio requires more research. On average, both zT/z0 and zQ/z0 corroborate Andreas’s (1987) theoretical model, (6.1), which includes no stratification effects. We therefore retain this model as our parameterization for zs/z0 in our bulk flux algorithm (cf. Jordan et al. 1999; Andreas et al. 2004a, 2005).
To evaluate the benefits of our new algorithm, we compared its predictions of the turbulent surface fluxes with predictions from the surface flux algorithm in the Community Ice Code (CICE). CICE is a module in the Community Climate System Model and therefore plays a crucial role in climate studies, including those done for the Intergovernmental Panel on Climate Change (IPCC) climate assessments. For these tests, we supplemented our SHEBA data with the comparable dataset from Ice Station Weddell.
On visually inspecting scatterplots (i.e., Figs. 7 –10) and calculating mean bias errors and root-mean-square errors (Table 5), we find that the new SHEBA algorithm is better than the CICE algorithm for representing both the Arctic and Antarctic turbulent fluxes—with the exception of u* for the ISW dataset. The CICE algorithm assumes a constant value for z0 of 5.0 × 10−4 m, whereas the maximum z0 value predicted by the SHEBA algorithm is about 2.3 × 10−4 m. Because the surface at ISW was aerodynamically rougher than at SHEBA (Andreas et al. 2005), the CICE algorithm does better than the SHEBA algorithm in representing the ISW u* measurements. In all other comparisons of measured and modeled u*, Hs, and HL values, our new SHEBA algorithm performed better.
Those comparisons of u* predictions suggest that the ultimate parameterization for z0 may need to be site specific. For instance, Banke et al. (1980) presented a rudimentary parameterization for z0 that depends on the measured physical roughness of the sea ice. Unfortunately, such measurements of physical roughness are still not routine and likely would not be available in a modeling environment. We are thus still looking for an effective way to include local sea ice topography in a z0 parameterization.
In closing, we have FORTRAN code for the bulk flux algorithm that we have described here and are willing to share it.
The U.S. National Science Foundation (NSF) supported our initial participation in SHEBA with awards to the U.S. Army Cold Regions Research and Engineering Laboratory, NOAA’s Environmental Technology Laboratory (now the Earth System Research Laboratory), the Naval Postgraduate School, and the Cooperative Institute for Research in Environmental Sciences. NSF also supported our use of the Flux-PAM stations from the facilities pool at the National Center for Atmospheric Research. Both NSF (Award 06-11942) and the National Aeronautics and Space Administration (Award NNX07AL77G) supported ELA at NorthWest Research Associates during the preparation of this manuscript. Lastly, we thank William Gutowski and three anonymous reviewers for their thoughtful comments on this manuscript.
Corresponding author address: Dr. Edgar L Andreas, NorthWest Research Associates, Inc. (Seattle Division), 25 Eagle Ridge, Lebanon, NH 03766-1900. Email: firstname.lastname@example.org