The 4-month drift of Ice Station Weddell (ISW) produced over 2000 h of nearly continuous measurements in the atmospheric surface layer and in the snow and sea ice in the western Weddell Sea. This paper reports simulations, based on these data, of processes in the air, snow, and sea ice at ISW using SNTHERM, a one-dimensional mass and energy balance model. An earlier version of SNTHERM had to be adapted, however, to treat the flooding that often occurs on sea ice in the western Weddell Sea. To treat this layer of slush and brine, SNTHERM holds the brine salinity constant at its initial value of 31.5 psu until 80% of this slush layer freezes. The current version of SNTHERM also incorporates a new parameterization for the roughness length for wind speed, z0, derived from analyses of ISW eddy-covariance data. SNTHERM's simulations are validated with temperature measurements within the ice and snow and with eddy-covariance measurements of the surface momentum and sensible and latent heat fluxes. The simulated turbulent fluxes agree fairly well with the measured fluxes, except the simulated sensible heat flux is biased low by 4–5 W m−2 for both stable and unstable stratification. The simulated temperature profiles in the snow and ice also agree well with the measured temperatures. In particular, allowing seawater to flush the slush layer until it is 80% frozen delays the freezing of this layer such that its behavior mirrors the data.
Late in the drift of Ice Station Weddell, seawater began seeping into our instrument hut. This flooding was an unpleasant reminder of the main difference between Arctic sea ice and Antarctic sea ice: Even the perennial sea ice in the western Weddell Sea, where we had deployed Ice Station Weddell (ISW), is much thinner than perennial Arctic sea ice. The weight of our meteorological hut and the drifting snow that had accumulated around it eventually depressed the sea ice surface below sea level, and the ice flooded.
Even without the weight of our meteorological hut, though, ice in the Weddell Sea often floods (e.g., Eicken 1992; Ackley and Sullivan 1994; Lytle and Ackley 1996). The combination of thin ice and a few tens of centimeters of snow is often enough to push the snow/ ice interface below sea level. Because the thermodynamic effects associated with such flooding are not treated in many sea ice models, we incorporate them here in SNTHERM, a one-dimensional mass and energy balance model (Jordan 1991), and use it to produce a unified picture of snow, ice, and surface-level meteorological conditions on Ice Station Weddell.
Several other detailed models for temperate snow covers have come into use over the past decade, for example, CROCUS (Brun et al. 1989; Martin and Lejeune 1998), DAISY (Bader and Weilenmann 1992; Morris et al. 1997), and SNOWPACK (Bartelt and Lehning 2002). These models share similar parameterizations of in-snow processes with SNTHERM, and all simulate snow-cover properties such as density, grain size, and liquid water as these respond over time to meteorological forcing. The consensus after a decade of validation (mostly with observations of snow depth, density, temperature, and snow water equivalent mass) is that these models are consistent and perform reasonably well.
At higher latitudes, though, extreme cold and high winds lead to wind-packed snow covers with highly developed depth hoar. Thus, questioning the legitimacy of applying models developed for warmer climates to polar snow is reasonable. For example, Martin and Lejeune (1998) generally underestimate the snow surface temperature and show biased sensible heat fluxes in their 1-month CROCUS simulations of snow-covered Arctic sea ice. Nevertheless, successful simulations of a continental Antarctic snow cover with DAISY (Morris et al. 1997) and of Arctic and Antarctic sea ice with SNTHERM (Jordan et al. 1999, 2001, 2003) demonstrate the universal applicability (with minor modifications) of these physically based snow models. Their more sophisticated handling of the snow could supplement or replace the simple one-layer snow parameterizations now used in most sea ice models (e.g., Maykut and Untersteiner 1971; Ebert and Curry 1993; Flato and Brown 1996; Vihma et al. 2002). To this end, Maksym and Jeffries (2000) recently added a more detailed snow module to their sea ice model and, with it, modeled flooded sea ice in the Ross Sea.
Jordan et al. (2001) report on our preliminary efforts to simulate the ice and snow processes and the surface energy budget on Ice Station Weddell using SNTHERM. Here we expand that analysis by incorporating a new parameterization for the roughness length for wind speed, z0, that Andreas et al. (2004b, manuscript submitted to Bound.-Layer Meteor., hereafter AJM) develop on the basis of the ISW dataset. We also adapt SNTHERM to treat the flooded layer at the snow/ice interface that is commonly found in the Weddell Sea.
Finally, using measured incoming longwave radiation, incoming and reflected shortwave radiation, and mean meteorological variables as forcing, we use SNTHERM to simulate the turbulent surface fluxes and the temperatures in the ice and snow at ISW. We compare both the simulated turbulent fluxes and the in-snow temperature profiles with comparable ISW measurements. Evidently, only Martin and Lejeune (1998) and Jordan et al. (2003) before us have compared simulated turbulent surface fluxes from detailed snow models with actual eddy-covariance measurements of these fluxes. We close by showing time series of the components of the surface energy budget during our nearly 4-month deployment on ISW.
2. Measurements on Ice Station Weddell
Ice Station Weddell was a joint Russian–American, multidisciplinary study of the climatologically important western Weddell Sea (Gordon and Lukin 1992; ISW Group 1993). Many papers based on ISW have already appeared, for example, in oceanography, Gordon et al. (1993), McPhee and Martinson (1994), Muench and Gordon (1995), and Levine et al. (1997); in ocean chemistry, Weppernig et al. (1996); in biology, Fritsen et al. (1994); in ice physics, Lytle and Ackley (1996), Drinkwater and Lytle (1997), and Geiger et al. (1998); and in meteorology, Andreas and Claffey (1995), Makshtas et al. (1998, 1999), and Andreas et al. (2000). Here we report more results from the ISW meteorology program.
We maintained two meteorological sites on ISW. The main site included our instrument hut, our 5-m meteorological tower, and the Russian radiometers. Site Andy was about 100 m from this main site (Ackley and Lytle 1992; Lytle and Ackley 1996, their site A). Here we had the American radiometers and a thermistor string from the air, through the snow, and down into the sea ice.
Both the Russian and American teams on ISW installed broadband radiometers to measure longwave and shortwave radiation. Claffey et al. (1995) describe these measurements and, as part of our quality control, show scatterplots that compare results from similar instruments. We briefly summarize that material here.
We measured incoming and reflected shortwave radiation at the main meteorological site with the M-115M, a thermocouple pyranometer constructed by the Russian Main Geophysical Observatory. This instrument has a spectral band of 0.34–2.2 μm and a 180° field of view. This pyranometer is the standard shortwave instrument that has been used for many years at meteorological stations in the former Soviet Union and at Russian field sites and drifting stations in the polar regions.
We measured incoming and emitted longwave radiation with an experimental narrow-angle radiation balance gauge that uses thermocouple sensors and was built by the Leningrad Electrotechnical Institute. The instrument has a 45° field of view in both up and down directions and is sensitive to light in the spectral band 2– 40 μm. At the time of the experiment, this instrument was a relatively new design and had been used only experimentally in the polar regions.
Still at the main meteorological site, we had a Barnes Precision Radiation Thermometer (PRT-5) looking down at the snow surface from a height of about 1 m. This instrument has a 2° field of view and is sensitive to radiation in the range 9.5–11.5 μm. We used this instrument primarily for measuring the snow-surface temperature, but it also provides an independent measurement of the emitted longwave radiation.
At site Andy, we measured incoming shortwave radiation with an Eppley Precision Spectral Pyranometer (model PSP) and incoming longwave radiation with an Eppley Precision Infrared Radiometer (model PIR). Both instruments have a 180° field of view. The PSP is sensitive to light in the band 0.285–2.80 μm, while the PIR is sensitive in the 4–50-μm range. Both radiometers were calibrated after their deployment on ISW and were still within 1% of their original calibrations.
As forcing for SNTHERM, we use the American measurements of incoming longwave radiation (QL↓) from the Eppley PIR. We use the Russian measurements, however, to fill gaps in the series but corrected these using the correlation relation that Claffey et al. (1995) report.
Likewise, we use the PRT-5 measurements—extrapolated over the full infrared band—as our primary emitted longwave series (QL↑). We filled gaps in this series with the Russian data, corrected again using correlation results from Claffey et al. (1995). Note, however, that we do not use the emitted longwave radiation in our simulations; SNTHERM actually predicts this. We use the emitted longwave only in evaluating the surface temperature (as described in the next section).
We also use the American measurement of incoming shortwave radiation (QS↓) and filled gaps in this series with the Russian measurements, corrected again using the analysis by Claffey et al. (1995). Our only measurement of reflected shortwave radiation (QS↑) was with the Russian instrument.
b. Surface temperature
We measured the snow surface temperature Ts in several ways. Our primary measurement was with the Barnes PRT-5. The blackbody temperature that the PRT-5 measures, TBB, is related to the components of the net longwave balance of the snow surface according to
Here σ (=5.670 51 × 10−8 W m−2 K−4) is the Stefan– Boltzmann constant, and both TBB and Ts must be in kelvins. For the surface emissivity ɛ, we use 0.99 (Warren 1982; Dozier and Warren 1982). On rearranging (2.1), we could compute Ts. The emitted longwave radiation measured by the Russian radiometer provided an independent, radiative measurement of Ts through an equation like (2.1).
For surrogate measurements of surface temperature, we also deployed at the main site a General Eastern 1200MPS cooled-mirror, dewpoint hygrometer at a height of only 8–10 cm above the snow surface. This instrument measures air temperature with a platinum resistance thermometer and dewpoint (actually frost point) temperature with a second platinum resistance thermometer placed just under a cooled mirror. Both sensors are in an aspirated radiation shield. The near-surface air temperature approximates the surface temperature. The near-surface dewpoint temperature should be an even better estimate of the surface temperature because the near-surface humidity over snow is at the saturation value at Ts and the vertical humidity gradient over snow is typically very small (Andreas 1986).
c. Meteorological variables
The focal point of our main meteorological site was a 5-m tower. Here, at a height of 4.65 m above the snow surface, we had several instruments for measuring the mean meteorological quantities needed to drive SNTHERM. A second General Eastern 1200MPS measured air temperature and dewpoint temperature. We also had an R. M. Young 35003 propeller–vane on the tower. This yielded the mean wind speed and direction. Our 5-m tower rotated a full 360°. We used this vane to periodically align our instruments with the wind for best exposure.
AJM describe the turbulence instruments we also had on this tower at 4.65 m and how we processed their data. Briefly, we used a K-type sonic anemometer–thermometer made by Applied Technologies, Inc. (ATI), to measure the turbulent fluctuations in the velocity vector and in the air temperature. Within 25 cm of this sonic was a Lyman-α hygrometer made by Atmospheric Instrumentation Research (AIR) for measuring the turbulent humidity fluctuations. By doing covariance calculations using these turbulence quantities, we obtained direct measurements of the turbulent surface fluxes of momentum and sensible and latent heat that we will use to validate SNTHERM's simulations.
d. Subsurface temperature measurements
We use the in-snow and in-ice temperature profile measurements that Lytle and Ackley (1996) made at site Andy to initialize and validate SNTHERM. These temperatures came from a thermistor string constructed with Yellow Springs, Inc., 44033 thermistors embedded with 5–20-cm vertical spacing in a solid plastic pipe. These thermistors have a nominal resistance of 10 kΩ at −5°C and were read through a dc half bridge using a Campbell CR-10 datalogger that supplied a 100-mV excitation voltage. The reference resistor in this bridge was nominally 5 kΩ, and the circuitry was designed for this bridge and reference to be common to all the thermistors along the string. Hence, we estimate the absolute accuracy of these temperatures to be ±0.25°C. Lytle and Ackley (1996) provide other details on these temperature measurements, especially on the care they took in initially placing these thermistor strings in the ice.
3. Turbulence algorithm
SNTHERM predicts the turbulent surface fluxes of momentum (τ) and sensible (Hs) and latent (HL) heat through a bulk flux algorithm such that
Here, Ur, Θr, and Qr are average values of the wind speed, potential temperature, and specific humidity at reference height r; Ts and Qs are average surface temperature and specific humidity, the latter evaluated as the saturation value at Ts; and ρa, cp, and Lυ are, respectively, the air density, specific heat of air at constant pressure, and latent heat of sublimation. Equation (3.1a) also defines the friction velocity u∗.
The key to the bulk flux algorithm is evaluating the drag coefficient (CDr) and the transfer coefficients for sensible (CHr) and latent (CEr) heat at reference height r. These are related to the roughness lengths for wind speed (z0), temperature (zT), and humidity (zQ) according to (e.g., Garratt 1992, p. 52–56; Andreas 1998)
Here, k (=0.40) is the von Kármán constant, and ψm and ψh are known corrections to the wind speed and scalar profiles that account for stratification effects through the stability parameter r/L, where L is the Obukhov length. For unstable stratification (i.e., L < 0), we use Paulson's (1970) functions. For stable stratification (i.e., L > 0), we use Holtslag and De Bruin's (1988) formulation because it has the best properties in very stable stratification (Jordan et al. 1999; Andreas 2002). We make the usual assumption that ψh is the same for both the temperature and humidity profiles.
Once z0, zT, and zQ are known, computing the turbulent fluxes using (3.1) and (3.2) is straightforward. For these roughness lengths, we use the parameterizations that AJM develop or verify from analyzing the ISW turbulence data.
Specifically, for z0 we use
This equation gives z0 in m when u∗ is in m s−1, ν, the kinematic viscosity of air, is in m2 s−1, and g, the acceleration of gravity, is in m s−2.
This parameterization acknowledges that three processes influence the roughness of snow-covered sea ice and, therefore, builds on similar algorithms for the open ocean that Smith (1988) and Fairall et al. (1996) evaluate. The first term on the right-hand side of (3.3) shows that molecular viscosity sets the roughness length when u∗ is small and the flow is aerodynamically smooth. At the opposite extreme, in high winds with large u∗, snow drifts and blows. Through an energy conservation argument, the height of the resulting saltation layer scales with u2∗/g and sets z0 (e.g., Owen 1964); the 1 in curly brackets in (3.3) corresponds with this regime. Midway between these two extremes, the “fundamental” or permanent roughness of the surface sets z0 (Andreas and Claffey 1995; Andreas 1995; Andreas et al. 2003; AJM). The exponential (Gaussian) term in (3.3) operates over only a limited u∗ range around 0.18 m s−1 and models this permanent roughness. Coincidently, Andreas et al. (2004) report that (3.3) also fits a very large dataset collected during the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment. The only difference in the SHEBA data is that the ice seemed fundamentally smoother than at ISW; the SHEBA parameterization replaces the 5 in front of the exponential term in (3.3) with a 1.
SNTHERM has long used Andreas's (1987) theoretical model for the scalar roughness lengths zT and zQ. This predicts
where R∗ = u∗z0/ν is the roughness Reynolds number. Table I in Andreas (1987; or Table 2 in Andreas 2002) gives the polynomial coefficients in (3.4) when zs is either zT or zQ. The ISW analysis by AJM confirms that (3.4) is a good model for zT. Their associated derivation of zQ values provides tentative—though not conclusive—evidence that (3.4) is also useful for predicting zQ. We therefore continue to use (3.4) to model both zT and zQ in our SNTHERM simulations here.
4. SNTHERM for Antarctic sea ice
Jordan et al. (1999) describe modifications we made to SNTHERM to model the snow cover and sea ice at Russian drifting station North Pole 4. Strong winds, the cold, and the sea ice substrate necessitated these changes. Our simulations for Ice Station Weddell also required that, in addition, we treat the flooding seawater that infiltrates the snowpack. Jordan et al. (2001) report our preliminary adaptations to SNTHERM. Here we summarize them more fully. All units in what follows are SI.
Following Maykut and Untersteiner (1971), we parameterize the thermal properties of sea ice in terms of its temperature depression TD below 0°C and its bulk salinity S (in psu). The thermal conductivity ksi and apparent heat capacity csi of sea ice are thus,
In these, ki is the thermal conductivity and ci is the specific heat of pure ice, cw is the specific heat of pure water, Lf is the latent heat of fusion of water, and fb is the mass fraction of brine.
The final term in (4.2) treats the change in latent heat with temperature through a phase relationship between TD and fb. The partial derivative ∂fb/∂TD usually depends only on salinity and therefore becomes zero for pure ice. This derivative also becomes zero when the brine salinity is constant. Using the Clausius–Clapeyron equation, we derive the more general function
Besides the standard salinity term, (4.3) now contains an interfacial curvature term that extends the function to pure ice. The curvature term also provides a means for modeling the heat capacity of brine-flooded snow, where circulating seawater keeps the brine salinity close to 35 psu.
Because (4.3) assumes that the brine is trapped and that the salinity thus increases in the brine as it freezes, the flooded slush layer would cool too rapidly. To correct this deficiency, we assume that seawater continually flushes the slush while the brine volume is above a cutoff fraction fbc. The bulk salinity in the slush is then fbSsw, where Ssw is the salinity of seawater, and the temperature depression becomes
The thermal conductivity of snow ks varies widely for the snow types found at ISW: from as little as 0.06 W m−1 K−1 for new, fluffy snow to as much as 1.0 W m−1 K−1 for flooded snow. As Sturm et al. (1997) explain, most parameterizations of ks are empirical functions of porosity and, to a lesser extent, temperature. SNTHERM, for instance, routinely uses Yen's (1965) function to compute ks from the square of snow density. But because Yen's function overpredicts ks for the dense snow typical of ISW and does not treat wet snow, we switched to a modified version of the analytical model of Pitman and Zuckerman (1968) that Crocker (1984) developed for brine snow. Here we follow J. M. Hanesiak (2001, personal communication), who used this method in a model that combined SNTHERM with Flato and Brown's (1996) sea ice model.
Because the Pitman–Zuckerman (1968) model underpredicts heat flow at ISW, we set the dimensionless bond size (γ in their nomenclature) to 0.009, which is somewhat above their range of 0.005– 0.0071. Their model is also limited to snow densities less than 480 kg m−3. For denser snow, we extrapolate ks to ki with a quadratic function of snow density.
To handle drifting snow and ablation at the very windy North Pole 4 site, Jordan et al. (1999) compute a wind advection flux (J) as the gradient in wind transport near the surface,
Here, τt is the threshold surface stress (or momentum flux) at which saltation begins, and x is in the along-wind direction. The fractional stress gradient τ−1∂τ/∂x is “tuned” to adjust the advection flux. A positive gradient removes snow; a negative gradient deposits it. Because SNTHERM computes the threshold friction velocity [u∗t = (τt/ρa)1/2] from updated snow characteristics, it realistically predicts more advection for newer snow.
Although (3.1b) is the standard bulk parameterization for sensible heat flux, in SNTHERM we compute this as
Using SNTHERM, we simulated snow and sea ice thermal and mass variables at ISW between 25 February (Julian day 56) and 29 May 1992 (Julian day 150). We drove SNTHERM with hourly measurements of air temperature, relative humidity, wind speed (all at 4.65 m), incoming longwave radiation, and incoming and reflected shortwave radiation. Figure 1 shows time series of the meteorological conditions on ISW.
Precipitation data are crucial to the simulation but were not available from ISW. Jordan et al. (2001) therefore use Koschmieder's equation and a linear relationship between the extinction coefficient and the precipitation rate (Koh 1989) to estimate snowfall from the routine, 6-hourly visibility observations on ISW. The results of this estimate match observed snow depths (as judged from the thermistor profiles) but somewhat underpredict depth during blowing snow events toward the end of the experiment. To compare modeled in-snow temperatures with observations, however, we need more accurate estimates of snow depth. We therefore refined the estimates of Jordan et al. by increasing all precipitation by 30% in snow water equivalent and by ablating snow with the wind advection routine after 20 April (day 111).
Model results are very sensitive to the density of newly fallen snow. We estimate new-snow density using a function based on wind speed and air temperature (Jordan et al. 1999) but modified here such that the density increases as air temperature decreases below −10°C (Jordan et al. 2003). We set a lower limit of 150 kg m−3 in this parameterization since snow pit observations at ISW do not support values below this limit.
We simulated the snow cover at site Andy on the ISW floe. Snow pit data, ice cores, and thermistor profiles provided initial values of snow density, salinity, and temperature, the latter at 19 levels within the snow and sea ice. A negative freeboard at the start of our simulation caused flooding of the lower 0.18 m of snow (referred to as the slush layer) and produced isothermal ice of temperature −1.8°C. Flooding most likely occurred through numerous 2–3-cm drainage channels that formed in the ice during the previous summer (Lytle and Ackley 1996). The slush was overlain by 0.39 m of snow and underlain by 0.92 m of lightly deformed old ice. Initially, this slush was about 50% ice crystals and 50% brine and had a salinity of 31.5 psu. Upward-wicking brine gave the lower 10 cm of snow a bulk salinity of about 0.5 psu. The underlying ice had a bulk salinity of 2–3 psu (V. N. Churun 2000, personal communication). Because the slush layer later froze to become sea ice, we will refer to the top of this layer as the snow/ice interface.
a. Temperature traces in the snow
Figure 2 compares measured and simulated temperature traces at selected thermistor levels within the snow and at a 3-cm depth within the sea ice. Depths refer to the snow/ice interface, which was 2 cm above sea level when our measurements began. At the start of the simulation, the top three thermistors were above the snow and therefore show the air temperature. Between 28 April (day 119) and 24 May (day 145), the simulated snowpack covers the top thermistor. Most snow fell in two series of storms, between 3 and 16 April (days 94– 107) and between 23 and 29 April (days 114–120). The damping in the lower snow traces after these dates reflect the deeper snow.
SNTHERM closely replicates measured temperatures for most of the simulation. This close agreement between measured and modeled profile temperatures validates SNTHERM's ability to correctly simulate the thermal processes. Predictions, however, were slightly too warm during three periods—for 4–6 March (days 64–66), 26 March–3 April (days 86–94), and 23–28 May (days 144–149)—and too cold from late on 18 April (day 109) through around 27 April (day 118). Interestingly, the periods that we predict to be too warm coincide with (or slightly follow) three events of high winds and low air temperatures.
The most obvious cause for misestimated heat fluxes in the snow is an error in thermal conductivity. Such an error results from problems either with our predictive function for ks or from errors in snow density, on which ks depends strongly. Because of uncertainty in the predictions of new-snow density, the latter seems to be more likely. Unpublished density measurements from the upper 5–14 cm of the snowpack between 26 February and 10 April ranged between 197 and 340 kg m−3 (V. N. Churun 2000, personal communication; V. I. Lytle 2000, personal communication). SNTHERM predicts an average of 238 kg m−3 for the top 10 cm of snow during this period, a density that is within the bounds but slightly low. Increasing the density, however, led to poorer root-mean-square (rms) errors in snow temperature.
Averaged thermal conductivities for the upper and lower snowpack for the experiment (excluding the frozen brine snow) are, respectively, 0.21 and 0.37 W m−1 K−1. Compared with the constant ks of 0.33 W m−1 K−1 suggested by Sturm et al. (2002), our simulation, thus, predicts lesser and greater heat flow for the two respective levels. While an alternative simulation that used their ks and a constant snow density of 330 kg m−3 everywhere predicted snow temperatures more accurately during windy periods, temperatures overall were too cold. This was particularly the case after 9 April (day 100), when simulated temperatures were consistently 2°–3°C below the observations.
In a similar analysis of heat flow during the SHEBA experiment, SNTHERM also failed to predict enough cooling during windy periods (Jordan et al. 2003). Because the SHEBA site was windier (with speeds reaching 14–15 m s−1) and had shallower snow, wind effects were more pronounced. Jordan et al. suggest that wind pumping, a mechanism not treated in SNTHERM, caused the colder air to penetrate to depths in the snowpack of around 20 cm. Using Cunningham and Waddington's (1993) formula, they estimate in-snow velocities near the surface to be on the order of millimeters per second.
While the role of wind pumping in polar firn is well established (Albert 2002; Albert and Hawley 2002), there is some question as to whether winds can penetrate the dense, wind-packed snow on high-latitude sea ice (Sturm et al. 2002; Waddington et al. 1996). The anomalous temperature traces at ISW between days 144 and 149, when wind speeds were the highest, however, suggest wind pumping. When the winds decreased after day 147, modeled and observed temperature traces come back together. Further, the rms error for the upper snowpack increased to 2.17°C for periods when wind speeds exceeded 8 m s−1, while that error was 1.51°C for the entire simulation.
b. Temperature profiles in the snow
While Fig. 2 presents our results as temperature time series for the duration of ISW, we also show our results as temperature profiles in Fig. 3. These highlight the vertical stratification of thermal properties, for example, the thermal conductivity and heat capacity.
Figure 3 compares simulated and measured profiles from the ice through the slush and snow. Following Lytle and Ackley (1996, their Fig. 4), we plot the profiles at midnight at around 20-day intervals. The negative slopes of the profiles and their convex shapes indicate that the ocean is warming the ice and the atmosphere is cooling the snow surface. We also note an increasingly negative slope up the temperature profiles from the high-conductivity ice to the low-conductivity snow, consistent with Fourier's heat flow equation.
The profiles in Fig. 3 for days 60, 80, and 101 (29 February, 20 March, and 10 April) show the slush and brine freezing and the −1.8°C isotherm moving downward between about 29 February and 22 March (days 60 to 82). An important feature of Fig. 3 is the 17-day delay in cooling of the brine slush, which the circulating seawater kept near −1.8°C. The high heat capacity simulated by using (4.4) in (4.2) correctly captures this delay.
We model the temperature depression in the brine-slush layer with (4.4) until the ice fraction (1 − fb) reaches about 0.80. We established this ice cutoff for brine circulation (=1 − fbc) by comparing the initial brine salinity of 31.5 psu observed by Lytle and Ackley (1996) to their final bulk salinity of 6 psu. By contrast, the −1.8°C isotherm moved downward through the slush layer in about 3 days when we used the standard temperature depression, (4.3), in the heat capacity model (4.2).
In simulating flooding and snow/ice formation on Antarctic sea ice, Maksym and Jeffries (2000) model brine flow with two different approaches. The first, or standard model, assumes Darcian flow of brine through the ice. The second, or simple model, assumes that the brine percolates through cracks and that there is no interaction with the ice–brine network. This simple model produced the more reasonable results and is similar to ours in that they held the brine salinity of the slush constant until the slush froze. They, however, assumed the slush to be frozen when the ice fraction exceeds 0.60.
Continued air temperatures between −20° and −30°C cooled the ice surface to around −7°C by day 101. This is about 1.5°C warmer than the minimum temperature reached on day 94 (3 April), after which a series of snowfalls rapidly increased the snow depth and further insulated the ice. Minimal cloud cover (10%) and low wind speeds (3.5 m s−1) contributed to the large simulated temperature gradient on day 101 of 91°C m−1 at the surface—somewhat larger than the thermistors recorded.
The profile for day 118 (27 April) reflects a warming period that began around 23 April (day 114). Moderate snowfall at this time and drifting snow during the previous 6 h could account for the poorer fit in the snow on day 118. A comparison of modeled and measured temperature profiles for days 117 and 118 also suggests that predicted snow depths are several centimeters too low. We also note that the overly cool predicted temperatures continue a trend that began on day 109 (18 April), as discussed in the previous section.
The profiles for days 140 and 150 (19 and 29 May) reflect a return to colder temperatures with the approach of austral winter. The fit for day 140 is fairly close, except in the lower snowpack, where the slope of the simulated trace is too steep—suggesting that the thermal conductivity is too high. At this point in the season, older snow at the bottom of the pack has compacted to around 400 kg m−3. The Pitman–Zuckerman (1968) function that we use predicts, for snow densities above about 300 kg m−3, a steep rise in ks that may not be totally realistic.
Notice the decrease in snow depth between days 140 and 150. During this period, we eroded snow within SNTHERM using a stress gradient of 0.004 m−1. As we suggested in the previous section, the colder observed temperatures in the snow on day 150 reflect possible wind pumping on days 144–146, when wind speeds exceeded 9 m s−1. Strong negative curvature of the temperature profiles in the upper snowpack on day 144 (not shown) indicates convective heat transfer and, therefore, suggests possible wind pumping. Flatter curvature in the observed trace on day 150, compared with the simulated trace, possibly reflects this earlier cooling.
To match the snow depth to the thermistor profiles, we plotted diurnal in-snow temperatures for about half the days in the experiment. Figure 4 compares our simulations of snow temperatures with measured temperatures on 29 March (day 89). This day illustrates a possible near-surface snow process not included in SNTHERM: thermal disequilibrium between ice and air.
The time in each panel in Fig. 4 is with respect to solar noon. This day had wind speeds typical of ISW, 4–5 m s−1, but relatively strong incoming shortwave radiation around noon. Early in the day and late in the day, the thermistors in the air match well with the aspirated platinum resistance thermometer 10 cm above the snow surface. The simulated in-snow temperature profiles also agree well with the surface temperature measured by the Barnes PRT-5 for the times depicted in Fig. 4. The measured air temperatures, however, are 2°–4°C warmer than the measured and simulated surface temperatures. Because the deeper snow is both measured and simulated to be warmer than the surface too, the simulated temperature profile from the air, through the air/snow interface and down into the snow, looks like a lazy V (i.e., <), with the apex of the V (the coldest point) at the air/snow interface. We saw similar lazy V profiles on many other days, too.
Such a profile is, thermodynamically, very unstable. Heat will diffuse downward from the warm air to bring the surface temperature closer to the air temperature. And the very unstable profile in the snow, which shows a gradient of −169°C m−1 in Fig. 4, will foster convective mixing within the near-surface pore space that will likewise try to warm the surface. Such heat conduction toward the surface in both the snow and air could explain the effectiveness of the windless coefficient, E0, in (4.6).
In other words, Fig. 4 hints that the air and snow are not in equilibrium across the air/snow interface, contrary to the common belief for such entirely water-phase surfaces (e.g., Andreas and Cash 1996). SNTHERM, like most other one-dimensional snow models, assumes that the ice matrix and the interstitial air are at the same temperature and that the air is at rest. But the in-snow temperature profiles in Fig. 4 suggest that the simulated near-surface temperatures and the measured snow-surface temperature are often colder than the bracketing thermistor measurements. Both the simulation and the Barnes PRT-5 radiometer capture the ice temperature of the snow. The thermistors, on the other hand, more likely record the temperature of the air, especially near the surface where the porosity is large and contact with the snow grains may not be perfect.
Realize in this context, too, that the PRT-5 measurements of surface temperature and the simulations in Fig. 4 are independent. SNTHERM predicts the snow-surface temperature from the driving variables; it is not forced to fit our measured snow-surface temperature.
c. Turbulent surface fluxes
SNTHERM uses its prediction of the snow-surface temperature Ts to compute the sensible and latent heat fluxes via (4.6) and (3.1c), respectively. It also computes the surface stress with (3.1a). Beside using the algorithms for z0, zT, and zQ, (3.3) and (3.4), in our simulation, as a sensitivity study, we also carried out another simulation of the ISW dataset with the version of SNTHERM we had used for our North Pole 4 simulations (Jordan et al. 1999). In this, we also used (3.4) but set z0 to a fixed value of 1 mm.
Although we tuned (3.3) with u∗ data from Ice Station Weddell (AJM), Figs. 5 and 6 are not very different. The measured and simulated u∗ values in the two plots exhibit essentially the same correlation coefficient, 0.93, and similar fitting lines. In Fig. 5, the simulated values represent the measured values well at small u∗ but overpredict them at large u∗. In Fig. 6, the simulated u∗ values tend to underpredict the measured values at small u∗ but to overpredict them less than in Fig. 5 at large u∗. As a result, the fitting lines in the two figures are similar.
Equation (3.3) shows that z0 is smaller than 1 mm when u∗ is smaller than about 0.2 m s−1 and is larger than 1 mm when u∗ is larger than about 0.55 m s−1. This realization explains some of the differences between Figs. 5 and 6. What it does not explain, however, is why z0 = 1 mm (i.e., Fig. 5) does a reasonable job at predicting u∗ when u∗ is less than 0.2 m s−1 while the z0 parameterization (Fig. 6) is not as good.
AJM base (3.3) on a subset of the u∗ data depicted in Figs. 5 and 6, and these data clearly show that z0 is, on average, less than 1 mm (often by an order of magnitude) when u∗ is less than 0.2 m s−1. AJM therefore develop (3.3) to capture this behavior. The fact, evidenced by comparing Figs. 5 and 6, that our SNTHERM simulations for small u∗ are better with z0 = 1 mm than with z0 accurately modeled is problematic. One explanation may be that we did not measure u∗ very well when u∗ was small. Another may be that, in low winds (which generally mean very stable stratification at ISW), our stability corrections ψm and ψh in (3.2) are too large. These would thus make CDr, CHr, and CEr too small; a positive feedback loop then results in which SNTHERM predicts stronger stable stratification and even smaller CDr, CHr, and CEr values.
Figure 7 compares sensible heat fluxes measured by eddy covariance and simulated by SNTHERM. Figure 8 is a similar comparison for the latent heat fluxes. Figure 7 contains 1113 data pairs; Fig. 8, 1043. For the simulations presented in these two figures, we modeled z0 with (3.3). Similar plots (not shown) based on z0 = 1 mm were not appreciably different.
The relatively poor correlation coefficients for Figs. 7 and 8, 0.735 and 0.644, respectively, emphasize how difficult it is to both measure and model the turbulent heat fluxes over snow and ice surfaces at low temperatures. Because the vertical temperature and water vapor gradients in the air above the surface are often small and the fluxes are therefore also often small, the measurement and the simulation sometimes do not even agree on the sign of the flux. The small gradients and small fluxes, quite simply, are near the limits of our current measurement technology.
The simulated sensible heat fluxes in Fig. 7 are low by 4–5 W m−2, on average. The measured and simulated latent heat fluxes (Fig. 8), though more scattered, agree better, on average. The low bias in the simulated sensible heat flux could result from SNTHERM's underpredicting surface temperature. Recall, the temperature profiles in Fig. 4 suggest this could be a possibility.
Equation (4.6) shows that Hs depends linearly on Ts − Θr. If Ts − Θr is negative (stable stratification) and SNTHERM simulates Ts to be too low, Hs will be too negative. If Ts − Θr is positive (unstable stratification) and SNTHERM simulates Ts to be too low, Hs will be positive but too small. That is, regardless of the stratification, if Ts is biased low, the fitting line will be shifted downward, as in Fig. 7. A low bias in the sensible heat flux of 4–5 W m−2 corresponds to a low bias in surface temperature of about 0.6°C for a mean wind speed of 5 m s−1.
Martin and Lejeune's (1998) simulations with the snow model CROCUS and the SNTHERM SHEBA simulations of Jordan et al. (2003) show biases in sensible heat flux, too. As in Fig. 7, Martin and Lejeune and Jordan et al. predict the sensible heat flux to be biased low when Hs is negative. Curiously, though, Martin and Lejeune predict Hs to be biased high when Hs is positive, in contrast with Fig. 7, while Jordan et al. predict more scattered values here than in either Fig. 7 or Martin and Lejeune.
6. Components of the surface energy budget
The Weddell Sea is climatologically important because it is the source of Antarctic Bottom Water (e.g., Gordon et al. 1993; Meredith et al. 2000; Fahrbach et al. 2001), which is the lowest water mass in much of the World Ocean. Antarctic Bottom Water forms as a consequence of heat and salt exchange at the ocean surface. Modeling its rate of formation thus requires knowing the components of the surface energy budget. Through our simulations, we have evaluated these components over compact sea ice and present time series of them in Fig. 9 for the duration of ISW. Although heat exchange through leads and polynyas is another important element of bottom water formation, our work has not treated such surfaces.
Figure 9 shows daily averages of the net radiative fluxes, the turbulent surface fluxes, and the conductive flux at the snow/ice interface. Recall that we measured the incoming longwave radiation and the incoming and reflected components of shortwave radiation. SNTHERM computed the remaining fluxes depicted in Fig. 9.
The air temperature averaged −9.6°C for the 5 days of measurements we had in February 1992 but fell precipitously to an average of −20.1°C in March, −17.4°C in April, and −24.0°C in May. Low temperatures and relatively few clouds combined to produce large longwave losses and rapid cooling of the snow cover during March. More persistent cloud cover after 3 April mitigated these losses and, together with a deeper snow cover, led to actual warming of the lower snowpack and the sea ice (Fig. 2).
The sixth trace in Fig. 9 is the surface balance, which is
where the terms in parentheses here are the net shortwave and net longwave fluxes in Fig. 9. In our convention, positive B values mean the snow surface is loosing heat to the atmosphere. Figure 9, therefore, demonstrates that the snowpack was almost always cooling during our deployment. The flux from the ocean to the atmosphere, through the ice and snow, was also always upward during this period, as evidenced by the snow– ice interface trace in Fig. 9. Though they did not compute this exact same flux, Lytle and Ackley's (1996) calculations of the ocean heat flux and the conductive flux through the snow at several ISW sites reiterate that the conductive flux through the ice and snow was always upward.
On comparing the bottom two traces in Fig. 9, we see that this conductive flux at the snow–ice interface is almost large enough to explain the surface balance. The daily averaged imbalance between the total surface losses and the conductive flux is only 0.6 W m−2 during our deployment. Cooling in the snow accounts for this imbalance.
We see from Fig. 9 that the sensible heat flux predominantly warmed the surface to make up for the longwave losses. Modeled daily averages of sensible heat flux ranged from −37 W m−2 on nearly clear days to 7 W m−2 on overcast days. During our deployment, the daily averaged latent heat flux was small—1–3 W m−2 in magnitude—and changed from predominantly cooling the surface early in the experiment to warming it as the relative humidity rose to near 100% for the last 30 days of the experiment.
In SNTHERM, a one-dimensional mass and energy balance model, we have implemented parameterizations for the roughness lengths z0, zT, and zQ based on the eddy-covariance measurements on Ice Station Weddell reported by AJM. We then used SNTHERM to simulate snow, ice, and atmospheric processes for each hour for the duration of ISW. As a sensitivity study, we also did complete ISW simulations with a previous version of SNTHERM in which z0 is fixed at 1 mm.
Although the analysis of ISW roughness lengths by AJM clearly demonstrates that z0 is significantly less than 1 mm when u∗ is less than 0.2 m s−1 and greater than 1 mm when u∗ is above about 0.55 m s−1, our comparisons of measured and simulated u∗, Hs, and HL values with the variable and constant z0 formulations did not show much difference. Although computations with constant z0 are faster, we recommend the variable z0 formulation of (3.3) because it represents the relevant physical processes better and should therefore be more accurate for simulating conditions outside the domain we have tested.
Ice in the Weddell Sea commonly floods. We thus had to modify SNTHERM to treat this flooding since the intruding seawater is an advective, subsurface heat source. To treat the freezing of this slush and brine layer, we held the brine salinity constant at its initial value, 31.5 psu, until about 80% of the slush freezes. Physically, this procedure represents the continual flushing of the slush with fresh seawater until enough slush has frozen to close this connection with the ocean.
We replaced SNTHERM's standard thermal conductivity algorithm with Pitman and Zuckerman's (1968) routine, modified according to Crocker (1984) to handle wet and brine snow. This new conductivity algorithm let us treat brine snow at the base of the snowpack and reduced SNTHERM's tendency to overpredict the conductivity for dense snow. We also modified the new-snow density function used by Jordan et al. (1999) such that the density now increases with decreasing temperature below −10°C (Jordan et al. 2003); we set a minimum of 150 kg m−3 for new-snow density. These changes increased the thermal conductivity at the top of the snowpack and, thereby, improved SNTHERM's prediction for the depth of diurnal cooling.
With these modifications, SNTHERM matched recorded snow and ice temperatures exceptionally well. SNTHERM, however, failed to predict enough cooling in the snowpack during (or slightly after) three periods of high winds. We speculate that wind pumping caused this cooling.
Finally, we used SNTHERM to produce an ISW-long series of the components of the surface heat budget for compact sea ice in the western Weddell Sea. From these series (Fig. 9), we see that the radiative terms dominate in this area from late February through May. The net shortwave flux, of course, always warms the surface. The net longwave flux, on average, cools the surface during this period. The sensible heat flux tends to mirror the net longwave flux, offsetting it and thereby forcing the surface energy balance toward zero. The latent heat flux is small and can be either positive or negative.
We thank Boris V. Ivanov and Kerry J. Claffey for assisting with the data collection and processing, Victoria I. Lytle for supplying the ISW thermistor data, Vladimir N. Churun for providing his ISW measurements of ice and snow properties, Stephen F. Ackley for helpful discussions about ice and snow processes at ISW, Emily B. Andreas for helping with the data processing, and three anonymous referees for fair and thorough reviews. The National Science Foundation supported this research with Awards OPP-98-14014 and OPP-00-84190, and the U.S. Department of the Army supported it through Project 4A1611AT24.
Corresponding author address: Dr. Edgar L Andreas, U.S. Army Cold Regions Research and Engineering Laboratory, 72 Lyme Road, Hanover, NH 03755-1290. Email: email@example.com