16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 Influence of the ice - ocean heat flux on the ice thickness evolution in Saroma-ko lagoon, Hokkaido, Japan Matti Leppäranta1*, Kunio Shirasawa2 Division of Geophysics, Department of Physical Sciences, University of Helsinki, P.O. Box 64 (Gustaf Hällströminkatu 2) FI-00014 Helsinki, Finland; *Corresponding author, e-mail: matti.lepparanta@helsinki.fi 2 Pan-Okhotsk Research Center, Institute of Low Temperature Science, Hokkaido University, Kita-19, Nishi-8, Kita-Ku, Sapporo, Hokkaido 060-0819 Japan; e-mail: kunio@lowtem.hokudai.ac.jp ABSTRACT Sea ice grows and decays as forced by the fluxes through the boundaries. In particular, the flux at the lower boundary ­ the heat flux from the water body into the bottom of the ice sheet ­ is not very well known quantity. A thermodynamic sea ice model is employed to examine the influence of the oceanic heat flux on the thickness of the ice. Saroma-ko ice station data is used to analyse the physics and calibrate the model. The oceanic heat flux is normally 5­10 W/m2, and the ice thickness ranges in 30­50 cm; being that thin, the ice has a very active role in the thermodynamics. KEYWORDS: Ice - ocean heat flux; Ice thickness; Thermodynamic sea ice model; Measurements; Saroma-ko lagoon 1. INTRODUCTION Thermodynamic evolution of sea ice thickness is forced by the heat fluxes at its upper and lower boundaries and penetration of solar radiation into the ice sheet. Much research work has been done on the heat fluxes at the upper surface by measurement campaigns and analytical and numerical modelling. The lower boundary, i.e. the heat flux from the liquid water body into the bottom of the ice sheet is less understood. According to measurements this heat flux is of the order of 1­100 W/m2, which makes a significant contribution to the ice thickness evolution (e.g., Shirasawa et al., 2006). It is also as such important for water body itself being the principal mechanism of heat loss during the ice-covered season. In numerical models of sea ice thermodynamics the oceanic heat flux has often be used as a tuning factor. Maykut and Untersteiner (1971) employed a constant oceanic heat flux of 2 W/m2 to obtain the best-fit equilibrium thickness cycle for multiyear ice in the Arctic Ocean in their classical model. Modelling in the Antarctic seas has indicated that there the heat flux can be one order of magnitude larger. In Saroma-ko lagoon, on the northern coast of Hokkaido, Shirasawa et al. (2005) obtained the best fit 10-year statistics for the ice and snow thickness using a fixed heat flux of 5 W/m2. Ice models usually take a fixed oceanic heat flux and use that as a tuning parameter. It is, however, the next question that what is the influence of the variability of oceanic heat flux on the ice thickness and on the structure of ice as well. This problem is examined here by mathematical modelling using analytical and numerical models (Leppäranta, 1993; Shirasawa et al., 2005). This work belongs to the long-term joint research of the authors on ice­ocean thermodynamics (Shirasawa et al., 2006). 2. MATHEMATICAL MODELLING Simple analytic models provide rather good results for the climatology of sea ice thickness (Leppäranta, 1993). A generalized form including the oceanic heat flux and air ­ ice interaction reads in differential form 83 1 16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 dh T f - T a Qw =a (T a T f ) - dt h+d L (1) where h is ice thickness, t is time, a = /(L) 5.5 cm2/(ºC·day) is freezing-degree-day coefficient, is thermal conductivity of ice, is ice density, L is latent heat of freezing, Tf is the freezing point temperature, Ta is air temperature, d 10 cm is effective insulating thickness of the near surface air­ snow buffer, and Qw is oceanic heat flux. The Stefan solution h2 = 2aS, where S = ot max(0, Tf ­Ta)ds, comes from d = Qw = 0 and can be taken as an upper bound for favourable ice growth conditions. In numerical modelling the temperature profile and resulting heat flow is solved in a dense grid across the ice sheet (Maykut and Untersteiner, 1971; Shirasawa et al., 2005). The salinity of the ice is prescribed and together with the temperature determines the brine volume. Comparing with the analytical models, more realistic boundary conditions can be applied and the thermal inertia of the heat flow through the ice can be included. The heat conduction equation is c T 2T = 2 - q( z ) t z (2) where c is the specific heat of ice, z is vertical coordinate, and q is the absorption of solar radiation by the ice. The boundaries of the system move according to the heat fluxes resulting in changing of the thickness of the ice, at the bottom /z = Ldh/dt + Qw and at the top similarly except that for growth liquid water or slush must be available. 3. RESULTS Model investigations are made for idealized cases and comparing with observations in the Saroma-ko lagoon located on the northern coast of Hokkaido (Fig. 1). Its surface area is 150 km2 and its mean depth is 8.7 m. The lagoon is connected to the Sea of Okhotsk via two inlets, and the salinity is about 32 psu. At the beginning of the freeze-up, atmospheric forcing keeps the lagoon open and frazil ice is formed. Ice reaches a maximum thickness of 35­62 cm annually in the southeastern part of the lagoon. There are large spatial variations in ice thickness due to the oceanic heat flux from the inlets. An example of the data set is shown in Fig. 2. The period of observations of water current, temperature and salinity covered the ice season with freeze-up and breakup dates included. The deployment of the current meter was made on 6 December 1999 when the lagoon was open and the water temperature was about 3oC. Ice freeze-up started at late January, and the almost entire surface of the lagoon was covered with sea ice. The onset of ice breakup was at early April. The current speed was as high as 10 cm/s before ice freeze-up and reduced to 2-3 cm/s during the ice-covered period. The solid ice cover causes reduction of vertical mixing of the water under the ice, and thus there is not much momentum transfer from the wind to the water body. The main forcing of the circulation comes from the two inlets and some major rivers, and heating of the water body by the heat flux from the bottom sediments and the solar heating becomes active as soon as the snow has melted. Fig. 3 shows results from measurements of the oceanic heat flux in the Saroma-ko lagoon. The level is at about 10 W/m2 with fluctuations of 5 W/m2 around the mean. High frequency variations seem to be quite small. If the oceanic heat flux is a non-zero constant and air temperature is Ta = constant, the analytic model gives the steady-state solution h + d = (Tf ­Ta)/Qw (3) 84 16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 Figure 1. A location map of Saroma-ko lagoon. The picture shows sea ice areas in Saroma-ko lagoon and in the coastal region of the Sea of Okhotsk in white colour. Since h 0 and d 0, it is seen that the heat flux of Qw = (Tf ­Ta)/d is sufficient to prevent ice formation; this works as a simple polynomial equation. If the heat flux is less than that, the steady state ice thickness becomes h* = (Tf ­Ta)/Qw ­ d. With Ta = constant, the time evolution is obtained from Eq. (1) as (T f - T a ) h h t + log1 - =- 2 h* h * h* (4) In the beginning (h << h*) the ice grows as proportional to the square root of time. To understand the role of the variable oceanic heat flux, perturbation techniques can be utilized. Writing the ice thickness as h = h*+ , and assuming << h*, the ice growth equation can be approximated as d Qw T f -Ta =a 1 - - dt h * h * L = q1 (5) For cyclic oceanic forcing Qw/L = qo + q1sin(t), the thickness disturbance has a persistent cyclic part + 2 2 sin(t + ) (6) where = ­a(Tf ­ Ta)/h*2 is a relaxation constant for ice growth and is the phase shift. When both relaxation constant and frequency of oceanic heat flux are small, considerable thickness variations result. The amplitude q1 could be of the order of 1 cm/day, and therefore with , <<1/day the thickness variations would be 10 cm or more. Consequently, tidal heat transport may show up in the fortnightly cycle but not significantly in diurnal or semidiurnal cycles. For a numerical modelling effort, a thermodynamic sea ice model (Shirasawa et al., 2005) was used to examine the evolution of ice thickness. As an example, simulations for a 100-day case are shown in Fig. 4. Simulations were performed with zero, nonzero constant, and cyclic oceanic heat flux. The constant was taken as 10 W/m2 corresponding to a typical level, while in the cyclic case the mean was 85 16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 the same and the amplitude also the same, i.e. the heat flux varied between zero and 20 W/m2. It is clear that much oceanic heat influences the ice thickness cycle much. In the growth season the differences are small, but when the growth weakens the oceanic heat starts to melt the ice, and the ice break-up comes three weeks earlier. But what is remarkable that the cyclic oceanic heat flux did not much differ from the constant heat flux case as long as the averages are equal. With amplitude of 5 W/m2 the differences were smaller, and also the higher the frequency was the smaller the difference. The cyclic case of Fig. 4 is therefore exaggerated for illustrative purposes: the amplitude is 10 W/m2 and the cycle length so 30 days. Current Direction ( ° ) Current Speed (cm/s) Water Temp. (°C) Salinity (psu) Off Horoiwa , Saroma-ko 1999/12/6`2000/04/25 360 270 180 90 0 20 15 10 5 0 4 2 0 -2 33 32 31 30 29 28 8.0 Complete Ice Coverage Depth (m) 7.5 7.0 6.5 6.0 Dec.1 1999 Jan.1 2000 Feb.1 Mar.1 Apr.1 May1 Figure 2. Current speed and direction, temperature, salinity and depth obtained from the mooring station at the central area of Saroma-ko lagoon during the period from 6 December 1999 through 25 April 2000. The complete ice coverage indicates that the almost entire surface of the lagoon is covered with sea ice, except the areas near the inlets. [after Shirasawa and Leppäranta, 2003]. 86 16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 Saroma-ko 1999 18 16 14 12 10 8 6 4 2 0 1051 1226 1401 1576 1751 1926 2101 2276 2451 2626 2801 2976 3151 3326 3501 3676 3851 4026 176 351 526 701 876 1 W/m2 Feb 10 - Mar 11, 1999 Figure 3. Measurements of the oceanic heat flux in Saroma-ko lagoon during one month based on ice and near surface layer ocean measurements. Influence of oceanic heat flux on ice thickness 0,25 0,2 0,15 0,1 0,05 0 1 30 59 88 117 146 175 204 233 262 291 320 349 378 407 436 465 494 -0,05 523 Time (unit = 0.25 day) Figure 4. Model simulations for ice thickness (m) in a selected example case: zero heat flux (top thin curve, constant heat flux of 10 W/m2 (thick curve), and cyclic oceanic heat flux with 10 W/m2 mean, 10 W/m2 amplitude and 30-day period (lower thin curve). 87 16th International Northern Research Basins Symposium and Workshop Petrozavodsk, Russia, 27 Aug. ­ 2 Sept. 2007 4. CONCLUSIONS The oceanic heat flux from the water body to sea ice has been examined by mathematical modelling. The motivation is that ice grows and decays as forced by the fluxes through the upper and lower boundaries, and in particular the flux at the lower boundary ­ i.e. the oceanic heat flux ­ is not very well known quantity. Measurement data of Saroma-ko Lagoon on the northern coast of Hokkaido are also utilised in the analysis. Analytical models can be used to examine the equilibrium and spin-up of sea ice thickness under fixed external conditions and with perturbation techniques for cyclic oceanic forcing, too. In the cyclic case low frequencies (0.1 cycles per day or less) and low relaxation constants (i.e., fast response times for ice) can lead to remarkable ice thickness cycles. An important feature of the oceanic heat flux is that it influences not only the total ice thickness but also the stratification of the ice sheet, since it melts congelation ice at the bottom and provides thus potential for more snow-ice formation. Modelling work is ongoing in collaboration between the present authors to come up with more extensive conclusions and consider requirements for advanced 1-dimensional ice­ocean thermodynamic modelling. ACKNOWLEGMENTS This work belongs to the project "Ice Climatology of the Okhotsk and Baltic Seas" supported by the Japan Society for the Promotion of Science, the Japanese Ministry of Education, Culture, Sports, Science and Technology, and the Academy of Finland. REFERENCES Leppäranta, M. (1993) A review of analytical sea ice growth models. Atmosphere­Ocean 31(1), 123­138. Maykut, G.A. and Untersteiner, N. (1971) Some results from a time-dependent, thermodynamic model of sea ice, Journal of Geophysical Research, 76, 1550­1575. Shirasawa, K. and Leppäranta, M. (2003) Hydrometeorological and sea ice conditions at Saroma-ko lagoon, Hokkaido, Japan. Proceedings of the Seminar "Sea Ice Climate and Marine Environments in the Okhotsk and Baltic Seas ­ The Present Status and Prospects", Seili, Finland, 10-13 September 2001, Report Series in Geophysics, Division of Geophysics, University of Helsinki, No. 46, pp. 161-168. Shirasawa, K., Leppäranta, M., Saloranta, T., Kawamura, T., Polomoshnov, A. and Surkov, G. (2005) The thickness of coastal fast ice in the Sea of Okhotsk. Cold Regions Science and Technology, 42, 25­40. Shirasawa, K., Leppäranta, M., Kawamura, T., Ishikawa, M., and Takatsuka, T. (2006) Measurements and modelling of the water ­ ice heat flux in natural waters. Proceedings of The 18th IAHR International Symposium on Ice, Vol. 1, Sapporo, Japan, 28 Aug. ­ 1 Sept. 2006, pp. 85­91. Stefan, J. (1891) Über die Theorie der Eisbildung, insbesondere über Eisbildung im Polarmeere. Ann. Physik, 3rd Ser. 42, 269­286. 88