One-Dimensional Sea Ice-Ocean Model Applied to SHEBA Experiment in 1997-1998 Winter

A one-dimensional sea ice-ocean model with its application in the Arctic Ocean is presented. The model includes a mixed-layer ocean model, a multi-layer snow/ice model, and the interfaces among atmosphere, snow/sea ice, and sea water. The observational data from the measurements at the ice station of the Surface Heat Budget of the Arctic Ocean (SHEBA) field experiment between November 1997 and January 1998 were used to drive and validate the model. The energy budget of the stand-alone simulations shows that the longwave radiative cooling is mainly balanced by the heat released of freezing at the bottom of the sea-ice. The results also show that the effect of ventilation and blowing snow are required to reproduce the detailed observed surface temperature, thickness of the sea ice, sensible heat flux and upward longwave radiation.


IntrODuctIOn
In the Polar Regions the energy budgets of atmosphere and ocean are strongly affected by the presence of snow/ sea ice. As a result of high albedo, snow and sea ice dramatically reduce the amount of shortwave radiative energy available at the surface. Because of its low thermal conductivity, snow/sea ice also significantly restricts heat exchange between the ocean and the atmosphere. The variation of sea ice distribution, concentration, and thickness has been widely recognized as one of the strongest signals in climate changes (Houghton et al. 1990). In the past, several sea-ice models and observational studies Washington et al. 1976;Mellor 1986;Mellor et al. 1986;Price et al. 1986;Josberger 1987;Maykut and Perovich 1987;Morison et al. 1987;Omstedt 1990;McPhee 1992;Jin et al. 1994;Kiehl et al. 1996;McPhee et al. 1999; and others) have been presented to study the thermodynamics and dynamics of sea ice and interactions among the atmosphere, sea ice and sea water. Most numerical models were applied to simulate the seasonal or annual variation of sea ice instead of the daily evolution due to the limitation of observational data.
Recently, SHEBA (Moritz and Perovich 1996) provided daily observations, including upward and downward radiative fluxes, precipitation, wind, humidity, air and snow temperatures, sensible heat flux, as well as temperature, salinity, and density of sea water in 1997 -1998. The schematic diagram of snow and ice sites on SHEBA is shown in Fig. 1. SHEBA also provided data of snow and sea ice depths. Here, we present a one-dimensional snow/sea ice-ocean model and comparisons between the model simulations and observations during November 1997 and January 1998, when the observed temperature, salinity, and density of sea water are available at UCAR Joint Office for Science Support (JOSS). The simulations show that the radiative cooling at surface is mainly balanced by the heat released from freezing at the bottom of sea ice. Sensitivity tests show that simulated sea ice thickness and other fields are strongly affected by snow thickness. Simulations also reveal that ventilation can affect the heat transfer and temperature in snow.
The multi-layer snow/sea ice-ocean model is based on conservation of energy, mass, and momentum (Sun and Chern 2005). The mixed-layer ocean model is used to predict water temperature, salinity, density, and turbulent kinetic energy, while the snow/sea ice model is used to predict the thickness, temperature, and heat transfer of snow/ sea ice. Mathematical formulation and physical processes of the mixed-layer ocean model are presented in section 2. The formulas of the one-dimensional snow/sea ice model are in section 3. A simple parameterization for snow ventilation is presented in section 4. The comparisons between the modeled simulations and observations are discussed in section 5, which is followed by Summary.

OnE-DIMEnSIOnAl MIxED-lAyEr OcEAn MODEl
The models are applied to simulate the response of the change of snow, sea ice, and sea water due to the atmospheric forcing. It is noted that the current velocity is excluded due to the lack of observational data and the horizontal pressure gradient. It is also noted that shortwave radiation is not calculated due to the high latitude of SHEBA floe in the winter, although equations are included here for completeness. The equations for a one-dimensional ocean model, which includes 27 layers, are where T is water temperature, c w is the specific heat capacity of sea water, Δz is the layer-thickness, S is salinity, w t is the density of seawater, H is the heating, w S l l is salinity flux, respectively. δ(φ) is defined as Here subscripts k -1/2 and k + 1/2 stand for the value at the top and bottom of the k-layer. Because the advection is not included, we may include the Newtonian forcing in equations for temperature ( ) = 1000 Kg m -3 , and coefficients of a t and b t are as follows: where T c = T -273.15 K, H f -m is the net latent heat released from freezing or melting. The heat flux inside the water is where w T l l is the eddy heat flux, α w , albedo at the open water, is a function of solar zenith angle, etc. (Kraus and Businger 1994) (hereafter KB), and Rs w . is downward shortwave radiation inside water. The absorption of shortwave radiation inside water may be approximated in the form where (a 1 , a 2 ) = (0.4, 0.6), and (b 1 , b 2 ) = (15, 0.5) m -1 (Price et al. 1986). The downward and upward longwave radiations at the surface are .
where R air , . is downward longwave radiance reaching the sea surface, emittance ε sw is 0.98, and Stefan-Boltzmann constant σ is 5.67 × 10 -8 W m -2 K -4 . Hs, H, , and Hp are the surface sensible heat flux, latent heat flux, and heat transfer associated with precipitation at the surface, respectively (Lynch-Stieglitz 1994, hereafter L-S). They follow the similarity equations: where V a is the wind speed near the surface, and in which Pr and Ps are precipitation rate for rain and snow (in m s -1 ), c w and c sn are specific heat capacity of water and snow, respectively. T a is the air temperature, q sfc is 0.98 of the saturated mixing ratio at ocean surface, and kg m 170 is the density of the new fallen snow. From Eq. (3), we obtain the buoyancy The eddy coefficients ĸ m and ĸ T (m 2 s -1 ), dissipation and transport term, are parameterized according to turbulence closure of Sun (1989, 1993 a, and. The mass of the top-layer water will change due to evaporation, precipitation, melt or freeze: the snow/sea ice model, whenever the sea water is frozen, it was automatically taken as sea ice and resulting in the decrease of the top-layer water. On the other hand, melt of sea ice increases the mass of water beneath sea ice according to Ep. (17).
The freezing point of the seawater T sw (fre) is a function of salinity: T sw (fre) = 273.15 -55S, S in kg/kg (Maykut 1985) ( Similarity equations are applied at the atmospheric surface layer to calculate u * , θ * , and q * . The heat flux between sea ice and water was proposed by Josberger (1987) and McPhee (1992): where c hsw = 0.006, u *sw = 0.5 cm s -1 is friction velocity under sea ice, T ml is the mixed layer temperature. u *sw varies from 0 to 2 cm s -1 and T ml -T sw (

BASIc EquAtIOn fOr OnE-DIMEnSIOnAl SnOW/SEA IcE MODEl
The basic equations for snow pack are similar to those presented in Sun and Chern (2005), except ventilation and sea ice being included here. The mass, thickness, and density of a snow layer are defined as where subscripts "sn", "s", "w" are related to entire snow (= dry snow + liquid water), dry snow, and liquid water, respectively. The temperature (T sn ) can be calculated from equation of enthalpy within snow where enthalpy is: The heat fluxes at the top and bottom of layer k (except at the surface) is defined as where ventilation flux ( ) c w T a a z t d l l will be discussed in the next section, C sw ≈ C w = 4.218 × 10 6 J(m 3 K) -1 according to Verseghy (1991) and L-S. The heat capacity of solid phase of snow/ice C i is where kg m 920 The heat flux at the snow surface is where longwave and shortwave radiations are .
The snow emissivity ε sn = 0.98. Albedo of snow/ice is a function of wavelength and cloudiness (Grenfell and Perovich 1984) and varies from 0.2 to 0.9. Following Anderson (1976) and Liston and Hall (1995), we assume ] ] The decrease of shortwave radiance within the k th layer of snow/ice is calculated as The attenuation α sn of solar energy is a function of depth, property of snow/ice, cloudiness, solar zenith angle, etc. (Grenfell and Maykut 1977;Warren 1982;Jin et al. 1994;Ebert et al. 1995). For simplicity, α sn is assigned a representative value of 20 m -1 for snow and 7 m -1 for ice (Gu and Stefan 1990;Verseghy 1991). We finally get the change of mass within each layer M f -m is the change of phase caused by freezing-melting, the second term is applied at snow/ice surface, and (net water flow) is the net water accumulation from water flowing through that layer. Following Kojima (1967), Pitman et al. (1991), and L-S, the snow compactness between time interval Δt (= t n+1t n ) due to vertical stress and viscosity is parameterized: where gN k is the weight of the snow pack above the midpoint of layer k. Snow and ice start melting, when the temperature is above 0°C. A portion or entire water trapped inside snow/ sea-ice can freeze, if it is below 0°C. At the sea ice-water interface, seawater freezes and leaves salinity behind, when the temperature is below freezing point according to Eq. (18).

EffEct Of SnOW VEntIlAtIOn
Ventilation can be important to the exchange of heat, vapor, and chemical constituents between the atmosphere and snow (Cunningham and Waddington 1993;Waddington et al. 1996), but it is difficult to formulate. Gjessing (1977) and Colbeck (1989) suggested that pressure patterns associated with wind flow over topographic obstacles such as sastrugi and drifts could cause steady airflow within the snow. Colbeck (1989) concluded that mechanism of this flow could generate large air flow in the top few meters of a snow pack. He also estimated that the advection of sensible heat was sufficient to alter the near surface thermal region of snow pack. For simplicity, we assume that snow ventilation is simply caused by a wavy surface pressure perturba-tion, which is generated by topographic obstacles or snow drifts (Cunningham and Waddington 1993). The pressure perturbation for a steady flow within the snow is governed by a Laplace equation, i.e., If the surface perturbation is given by a sinusoidal wave in which λ is the horizontal wavelength , h sn is snow depth (not including sea ice), z is the location of the snow (0 ≥ z ≥ -h sn ). Once the pressure distribution is known throughout the snow pack, then the volumetric flux of air can be found by Darcy's law with viscosity of air μ =1.6 × 10 -5 Pa at -15°C and snow permeability K = 10 -9 m 2 (Waddington et al. 1996). The vertical velocity of air parcel is where the porosity (φ ≈ 0.62) in the polar region (Cunningham and Waddington 1993). The temperature perturbation inside snow can be estimated by where τ is the time scale. We then obtain the average heat transfer due to ventilation Ventilation increases heat transfer within and at the top of snow. Although the observed thickness of snow/sea ice varies considerably at different gauges during SHEBA experiment, it is difficult to determine the values of p o , λ, and τ in Eq. (37). Here we set τ = 5 min. which is the time scale of gravity waves in nocturnal stable layer, λ = 2 m, and p o = 5 Pa (with ventilation) or 0 (without ventilation). Albert (1996) used λ = 1.7 and 3.3 m and p o = 5 Pa to simulate the transfer of heat and species in polar firn.
The snow/sea ice model that consists of up to 6 layers will be activated when combination of solid phase of snow and ice reaches a thickness of 1.0 -5 m and over. We also separate the snow layer (snow and liquid water) from the sea ice layer (sea ice and liquid water) in the model. The thickness of the top and bottom layers is limited to 6 cm in order to have a better response to the external forcings from the atmosphere and the ocean.

nuMErIcAl SIMulAtIOnS BEtWEEn nO-VEMBEr 1997 AnD JAnuAry 1998
The data at Data Management Center include the vertical profiles of water temperature, salinity and density between 287.14 and 393.23 on SHEBA Days (January 1, 1997 is Day 1). The towers provided air temperature, dew point, wind direction and speed, pressure, long-and shortwave radiations, precipitation, and sensible heat flux started from SHEBA Day 304. The coupled snow/sea ice and ocean models were integrated from October 31 (SHEBA Day 304) to January 27, 1998 (SHEBA Day 392) for 89 days continuously, with the observed air temperature, wind, precipitation (or the observed snow depth, Fig. 2), and downward longwave radiation at 2 m height as the external forcing from the atmosphere. The measured sea water temperature and salinity on Day 304 were used as the initial condition for the oceanic mixed-layer model to calculate the sea water temperature, salinity, density, and turbulent kinetic energy. But we did not calculate current velocity due to the lack of the pressure gradient force and the initial current data. The snow and sea ice measured at Baltimore and Seattle sites will be used as the initial conditions to simulate the evolution of snow/sea ice and interactions between atmosphere and snow/sea ice.

Simulations for young Ice
The thickness of snow and sea ice measured at Gauge #28 on Day 304 is used as the initial condition. The average thickness of snow was 19.30 cm at Baltimore sites, while the thickness of sea ice was 53 cm with 920 kg m I 3 t = measured at Gauge # 28, which also represented four other gauges (#4, 52, 148, and 37). The initial snow density is assumed equal to 240 kg m -3 at the top 25% of snow pack and 360 kg m -3 in the lower layer.

case A: case without Blowing Snow
The surface temperature, heat flux, and freezing rate of sea water are strongly influenced by the thickness of snow and sea ice. Unfortunately, snow depth can be quite different from the accumulation of snow precipitation due to snow drifting. Case A is designed to show the results when the snow thickness is mainly calculated from the observed precipitation, which is shown in Fig. 3. Other control parameters used in this simulation include no ventilation (p o = 0) and γ = 1/(6 hr) in Eqs. (1) and (2). It is noted that the property of snow/ice is not sensitive to the values of γ (γ = 1/6 hr, 1/60 hr, and 0 being tested in the simulations), because freezing temperature T sw (fre) remains about the same. Figure 4a shows the number and thickness of layers inside snow (thickness > 0) and sea ice (thickness < 0). Layers of snow increase from two initially to three around December 13, while layers of ice decrease from four to three. It is also noted that the thickness of those layers changes with time, except the top and bottom layers remaining constant (= 6 cm), so the model can respond to the external forces efficiently. The enthalpy and mass are conserved while the number and thickness of layers change. Figure 4a also shows that the calculated snow thickness (indicated by the solid line above thickness zero) increases with time and reaches 46.7 cm on Day 392, which is much deeper than the measurement, 26.6 cm. A thicker snow pack effectively reduces the heat transfer across it, hence less sea water freezes compared with the measurement. The simulated ice thickness (indicated by solid line below thickness zero), 93 cm, on Day 392, is much less than 117 cm observed (indicated by crosses). The simulated average sea ice thickness is 74.5 cm versus the observed one, 81.6 cm (Table 1a). Meanwhile, the calculated surface temperature is colder than measurement, especially during    the cold episodes after December (Fig. 4b). Consequently, it also reduces the upward longwave radiation, as shown in Table 1a. The snow is quite light with a mean density of 279.6 kg m -3 , because a lot of fresh snow is present.

case B: control case for young Ice
Here we have included the effects of ventilation, Newtonian forcing of sea water temperature and salinity (γ = 1/6 hr), and the effective precipitation rate (Fig. 5), which is derived from the observed snow depth form Baltimore site shown in Fig. 2. In Fig. 6a, obviously, the modeled depth of sea ice is comparable with the observed depth, although the simulated ice thickness shows fewer fluctuations due to missing the detailed variation of snow thickness in Fig. 2. The observed ice thickness is 117 cm versus the simulated one, 113 cm, on Day 392. The mean sea ice thickness is 81.6 cm and the simulated one is 81.5 cm after 89 days integration. If we define error (Err) and root mean square (rms) of a variable f as: The calculated mean surface temperature is 244.10 K versus the observed temperature, 244.12 K. The error of the surface temperature is -0.02 K and the root mean square is 1.5 K. The modeled mean sensible heat flux is -3.35 W m -2 versus the measurement of -5.2 W m -2 ( Fig. 6c) with an error (Hs) of 1.85 and rms (Hs) of 5.9 W m -2 . Be noted that ventilation contributes -1.29 W m -2 to the surface sensible heat flux. Figure 6d and Table 1a show that the observed upward longwave radiation is well reproduced by the simulation. The upward longwave radiation is much larger than the downward longwave radiation with an average net radiative cooling of -31 W m -2 ( Fig. 6a and Table 1a), which is the driving mechanism for the development of sea ice shown in Fig. 6a. The radiative cooling at the snow surface penetrates through snow and sea ice, and causes water beneath sea ice to freeze. The mean latent heat released from water freezing is 24.1 W m -2 and the sensible heat flux transported from water to sea ice is about 1.7 W m -2 . Since the snow surface temperature is slightly lower than air temperature, a small amount of sensible heat flux (-3.35 W m -2 ) is transported into snow but the surface latent heat flux is negligible (< 0.1 W m -2 ).
The observed and simulated sea water temperature, salinity, and density are shown in Figs. 7 and 8. It is not surprising that the simulated profiles follow the observed closely, because the Newtonian forcing is applied in both temperature and salinity equations. Note that variation of water temperature in the mixed layer is less than 0.1°C because the freezing point remains about the same.

case c: Same as control case but without newtonian forcing in Sea Water temperature and Salinity Equations
Without Newtonian forcing, the one-dimensional ocean model produces a mixed layer with uniform temperature and salinity (Fig. 9). We can also see that the modeled salinity in the mixed layer increases with time as water freezes and leaves salt behind. The temperature and the freezing point drop slightly according to Eq. (18). However, the observed salinity and water temperature change with time irregularly. Water temperature at the top mixed layer increased from -1.5 to -1.45°C in December, then decreased to -1.53°C in late January. Furthermore, the difference between freezing point and the temperature in the upper mixed layer, T ml -T sw (fre) , is 2.4 × 10 -3 K with γ = 0 here, and T ml -T sw (fre) = 3.3 × 10 -3 K with γ = 1/(6 hr) in Case B after 89 days integration. Because the simulated temperature is very close to the observed one in the oceanic mixed layer, the simulated ice thickness, snow temperature, and other properties are almost identical to those in Case B. More detailed comparison between Cases B and C is presented in Table 1a.

case D: Same as case B but without Ventilation
The simulated ice thickness on Day 392 is 109.3 cm   and the mean thickness after 89 days integration is 79.7 cm (Fig. 10a). They are thinner than the observed values, 117 and 81.6 cm, respectively, or the simulated results presented in Cases B and C (Table 1a). On the other hand, the simulated mean surface temperature, 243.45 K, is colder than the mean measured temperature of 244.12 K (Fig. 10b). Although the simulated mean surface temperature is less than that of Case B, the mean downward sensible heat is only 3.01 W m -2 compared with 3.35 W m -2 in Case B. Hence, we may well conclude that ventilation can increase the heat transfer in snow. Also, the model with ventilation predicts the growth of sea ice better than the model without it.

Simulations for Multiyear Ice
The sea ice thickness of 145.5 cm measured at Gauge #34 and snow depth of 14.2 cm measured at Seattle Site on Day 304 were used as the initial condition for this simulation. The density of snow was 240 kg m -3 at the top 25% and 360 kg m -3 in the lower 75% of snow pack initially.

case E: control case for Multiyear Sea Ice
Compared with Case B, here we have used a thicker initial sea ice and a different effective precipitation rate, which is derived from the snow depth at Seattle site in Fig. 2. The simulated sea ice thickness grows from 145.5 cm on Day 304 to 185.0 cm on Day 392, which is in good agreement with 185.8 cm measured on Day 392 ( Fig. 11a and Table 1b). The modeled and the observed mean ice depths are 163.2 and 164.8 cm, respectively. The simulated mean surface temperature is 243.4 K and the observed one is 244.1 K with an error of -0.68 K and root mean square of 1.57 K (Fig. 11b). The calculated mean surface sensible heat flux is -7.18 W m -2 versus the measurement of -5.2 W m -2 . Ventilation contributes -2.47 W m -2 to the sensible heat flux. The mean latent heat flux calculated at surface is -0.28 W m -2 . A thicker sea ice can also hinder the heat exchange between the atmosphere and the sea water beneath sea ice, hence, the sea ice grows slower than that of Case B. The average heat transfer from water to sea ice is 1.69 W m -2 and the latent heat released from water freezing is 15.96 W m -2 (Table 1b). It is noted that the errors of surface temperature, sensible heat flux, and sea ice thickness simulated from this case are larger than those from the Case B over the young ice. However, the error of surface temperature remains less than 1 K, the error of the surface sensible heat flux is about 2 W m -2 , and the deviation of sea ice thickness is less than 2 cm after 89 days integration, which might be acceptable for a regional climate model. Figures 12a -b show that the modeled surface tem-   perature is colder than the measurement, especially during the cold periods. The mean surface temperature, 242.6 K, is 1.5 K colder than observed value of 244.1 K. Since the temperature difference between atmosphere and snow is larger, the calculated sensible heat flux, -5.53 W m -2 , which is close to observation, -5.20 W m -2 , but is larger than that of Case E. Although the snow surface is colder than Case E, less water freezes at the sea ice base, as is seen that a sea ice thickness of 183.6 cm on Day 392 versus 185.2 cm for Case E, and the mean sea ice thickness of 162.5 cm compared with 163.2 cm for Case E. The latent heat released by water freezing, 15.35 W m -2 , is also slightly less than 15.93 W m -2 in Case E. Therefore, ventilation, which contributes about one third of the simulated surface heat fluxes, enhances the heat transfer inside the snow. It also reduces the cold bias at snow surface where a nocturnal inversion can easily form by radiative cooling in winter. Table 1b shows that numerical results of Cases E and G are almost identical, which may confirm that the develop-ment of snow/sea ice depends upon the radiative cooling at the surface, freezing point temperature and freezing rate of sea water instead of the detailed structure inside the ocean.

SuMMAry
A one-dimensional snow/sea ice-ocean model with its application to the Arctic Ocean is presented. The model includes a mixed-layer ocean model, a multi-layer snow/ice model, and the interfaces.
Numerical simulations of the SHEBA experiment between November 1997 and January 1998 show that net longwave radiation deficit caused cooling of snow/sea ice and freezing of the sea water beneath the sea ice. The latent heat released from freezing at ice-water interface was transferred upward through snow/sea ice with a strong temperature gradient. The atmosphere and sea water also provided a small amount of sensible heat fluxes to snow/sea ice. The simulations also show that the accurate snow thickness is crucial to predict the freezing rate of the sea water beneath. Meanwhile, including ventilation can increase the heat transfer inside the snow. The results also indicate that the horizontal advection and drifting of sea ice should be included in the ocean model to reproduce the observed property of sea wa- ter and current velocity, even though the snow and sea ice are not sensitive to the detailed profiles of salinity and water temperature. Finally, more observations and advanced models are needed to simulate the drifting snow and ventilation in the Arctic realistically.