1. Introduction
Sea ice is a geophysical component playing an important role in the weather and climate system. Sea-ice formation is due to the thermal interaction between the atmosphere and the underlying ocean. The ice acts as a barrier, reducing the transfer of moisture, heat and momentum between the atmosphere and the oceans. Because of its high surface albedo, ice reflects a large part of the incoming solar radiation. The thermal regime within the ice is dominated by the variation in its thermal properties and by the external forcing. Freezing at the ice bottom tends to reject salt, whereas the melting of sea ice will decrease the ocean surface salinity. The ice–ocean interaction can affect the stability of the upper oceans and alter the water circulation throughout the world oceans.
In recent decades, numerical sea-ice thermodynamic models (e.g. Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference SemtnerSemtner, 1976; Reference GabisonGabison, 1987; Reference Ebert and CurryEbert and Curry, 1993; Reference Flato and BrownFlato and Brown, 1996; Reference Launiainen and ChengLauniainen and Cheng, 1998; Reference Bitz and LipscombBitz and Lipscomb, 1999; Reference WintonWinton, 2000) have become established as a useful tool for understanding the mechanism of the heat exchange between the air and the ice, the thermal variations in ice growth and melting, as well as the ice–ocean interactions. The numerical integration of the heat-conduction equation is one of the main components of ice thermodynamic models. Before applying any numerical scheme, one has to define the model’s spatial (grid sizes) and temporal (time-step) resolution. In general mathematical terms, increasing the numerical resolution of a scheme should yield a better accuracy in the result and vice versa. In the case of process studies, fine-resolution models are usually applied, in order to obtain details. A prototype of such a model was developed by Reference Maykut and UntersteinerMaykut and Untersteiner (1971). For climate studies, coarse-resolution models are often used, in order to enable faster long-time simulations. A typical example is the model of Reference SemtnerSemtner (1976).
Physically, the thermal variation of sea ice is a complex process. The interfacial surface fluxes between the ice and the atmosphere and the in-ice heat conduction are strongly affected by the ice thickness (Reference MaykutMaykut, 1978). The in-ice thermal regime responds to the heat-flux balance at its surface and bottom. Its multi-phase constituency (ice crystal, solid and liquid brine, inclusion of air and other impurities) varies concurrently with the changes in heat flow and energy storage inside the ice. One has to define a thin interfacial layer (surface layer) with a finite mass and heat capacity to study the surface heat balance, especially for a fine-resolution model. The ice mass and heat storage clearly depend on the thickness of the surface layer, and this dependence, in general, is non-linear because the brine pockets and vapour bubbles affect the physical and thermal properties of the ice. The solar radiation flux, which penetrates into the interior of the ice or snow, is strongly attenuated in the surface layer. The surface heat balance therefore tends to be established differently depending on the various assumptions made for the model spatial resolution. Below the surface, the absorbed solar radiation is an internal heat source and plays an important role in the heat conduction of the upper ice layer. The estimation of this term using a finite-difference scheme is directly affected by the spatial resolution used. To systematically demonstrate how spatial resolution affects model results is one of our motivations for this study.
The effect of spatial resolution on ice growth, or in general on the parabolic diffusion equation, is not a new problem, but such studies often appear in technical reports not easily available to the general public. Reference SemtnerSemtner (1976) tested the effect of spatial resolution on the annual cycle of ice thickness. He concluded that the asymptotic limit of the modelled annual cycle of the equilibrium ice thickness is already approached with two model layers. Increasing the number of model layers from three to nine, results in a deviation of only 2% from the result of a three-layer model. These numerical tests were performed with simplified forcing conditions without taking into account the absorbed solar radiation in the ice layer. The effects of the number of model layers on the surface temperature, surface heat balance and ice-temperature profiles were not mentioned. Reference SavijärviSavijärvi (1992) evaluated several numerical schemes for the soil-diffusion equation in atmospheric models. Special attention was paid to the number of soil levels needed. The analytical solution of a parabolic soil-diffusion equation was compared with results obtained using various spatial resolutions. The simulated temperatures at different depths produced by a Crank–Nicholson (CN) numerical scheme with 3–10 layers were quite close to the analytical solutions.
The effects of the time-step (temporal resolution) on model results are often discussed mathematically or associated with the model spatial resolution. For example, Reference PatankarPatankar (1980) pointed out that with a small time-step the six-point symmetrical CN scheme has a better accuracy than the fully implicit scheme. Reference SmithSmith (1985) indicated that oscillations in the result for a CN scheme may occur if long timesteps with small gridcells are used. Physically, Reference Hanesiak, Barber and FlatoHanesiak and others (1999) have pointed out that an ice-growth model yields considerably different results using external forcing with a variable length scale that is linked with the model time-step. They compared simulations of the surface heat balance and the annual cycle of sea-ice thermodynamics by running Reference Flato and BrownFlato and Brown’s (1996) ice model with hourly and daily average forcing data. They found significantly different results for break-up dates, open-water duration and snow ablations. These differences are probably due to the physical effect of the forcing data, i.e. the model response to the synoptic external forcing is different from that to daily or monthly external forcing.
In section 2, the ice model applied in this paper is introduced; a two-layer parameterization scheme of solar radiation penetrating through the sea ice is described, and the numerical resolutions of various ice models with their grid systems are summarized. In section 3, the model results are first compared with the possible analytical solutions. The purpose of such tests is (1) to verify the numerical solution and (2) to systematically demonstrate the accuracy of the results when the resolution is increased. The full numerical simulations are then presented and analyzed. Model runs are performed with constant weather-forcing data so that results depend only on changes in resolution. The conclusions are presented in section 4.
2. Theoretical Background
2.1. The sea-ice model
The ice model presented by Reference Launiainen and ChengLauniainen and Cheng (1998) is used in this study. The basic model equations read:
The ice temperature is governed by the heat-conduction Equation (1), while Equations (2) and (3) represent the boundary conditions for the temperature, and the heat and mass balance at the surface and bottom.
In Equation (1), T is the temperature, ρ is density, c is specific heat, k is thermal conductivity, q(z, t) is the amount of solar radiation penetrating below the surface, t is time, and z is the vertical coordinate below the surface (positive downwards). Subscripts i and s denote ice and snow, respectively. The thermal conductivity and heat capacity of sea ice are given as: k i = k i0 + βs i/(T i − 273.15) and (ρc)i = ρ 0 c 0 + γs i/(T i − 273.15)2, following Reference Maykut and UntersteinerMaykut and Untersteiner (1971), where k i0, ρ 0 and c 0 are the thermal conductivity, the density and the specific heat of pure ice, respectively, si is the ice salinity, and β and γ are constants.
In Equation (2), α is the surface albedo. The downward shortwave radiation (Q s) is calculated by an empirical formula (Reference ShineShine, 1984) with the cloudiness factor (C) of Bennett (1982): Q s = S 0 cos2 Z/[(cos Z + 1.0)10−3 e + 1.2cos Z + 0.0455] (1 − 0.52C), where S 0 is the solar constant, Z is the local solar zenith angle, e is the vapour pressure and C is 0–1.
The solar radiation penetrating below the surface (I 0) is that part of the energy that contributes to internal heating of the ice/snow. It is the portion of solar radiation that does not contribute directly or immediately to changes in mass or temperature at the surface.
The incoming atmospheric longwave radiation (Q
d) is calculated by the formula of Reference EfimovaEfimova (1961) with the cloud effect according to Reference Jacobs, Barry and JacobsJacobs (1978): 
, where σ is the Stefan–Boltzmann constant and T
a is the air temperature. The outgoing longwave radiation from the surface is estimated by the Stefan–Boltzmann law 
with a constant surface emissivity (ε = 0.97), where T
sfc is the surface temperature.
The turbulent sensible-heat (Q
h) and latent-heat (Q
le) fluxes are calculated by the bulk formulae: Q
h = −ρ
a
c
p
C
H (T
sfc − Tz
a)Vz
a and 
, where ρ
a is the air density, c
p is the specific heat of air, L
v is the enthalpy of vaporization, CH
and C
E are the turbulent transfer coefficients, V is wind speed, 
is the specific humidity, and the subscripts sfc and z
a refer to the surface and a height of z
a in the air, respectively. The coefficients C
H and C
E are estimated using the Monin–Obukhov similarity theory with stability effects based on Reference HögströmHögström (1988) in unstable cases and Reference Holtslag and de BruinHoltslag and de Bruin (1988) in stable cases. An aerodynamic roughness value of 0.001 m is used and the thermal roughness length is calculated according to Reference AndreasAndreas (1987).
The surface conductive heat flux is Fc , and the heat flux due to surface melting is F m. When T sfc tends to be larger than the freezing temperature (T f), T sfc remains as T f, and the heat used for melting is F m = ρ i,s L f dH i,s/dt, where L f is the latent heat of fusion assumed to be constant, and H i,s is the thickness of ice or snow. The heat-balance Equation (2) can also be seen as a polynomial of the surface temperature Γ(T sfc). The temperature T sfc is calculated iteratively from Equation (2) and used as the upper boundary condition for the numerical scheme of Equation (1).
In Equation (3a) and (3b) the ice bottom temperature (T bot) is constrained to be the freezing temperature, and F w is the oceanic heat flux, assumed constant.
An implicit six-point symmetrical CN scheme, which was derived on the basis of the integral interpolation method (Reference ChengCheng, 1996), is used to solve the heat-conduction equation. The scheme derived by this method has been mathematically proven to be conservative (Reference Li and FengLi and Feng, 1980). In addition, such a numerical technique ensures that the scheme can be easily extended to an uneven spatial grid size (Reference ChengCheng, 2000).
The model parameters are listed in Table 1. The sea-ice properties are based on average values from the literature and field measurements in the Baltic Sea. For the seasonal ice evolution, this model has been verified using historical climatic weather-forcing data, while for a short-term ice simulation, field measurements have been used (Reference Cheng, Launiainen, Vihma and UotilaCheng and others, 2001). In these studies, the snow and ice were divided into 5 and 10 layers, respectively, and the time-step was 600 s. The model yielded results in good agreement with the observations.
Table 1. Model parameters based on field measurements in the Baltic Sea and values found in the literature
2.2. Penetrating solar radiation in ice and snow
The solar radiation penetrating into sea ice depends strongly on the wavelength of the irradiance, on its angle of incidence, on the structure of the sea ice and on sky conditions. For ice-modelling purposes, however, simple parameterizations based on the Bouguer–Lambert law are often used (e.g. Reference UntersteinerUntersteiner, 1964; Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference Grenfell and MaykutGrenfell and Maykut, 1977). In these studies 
, z > z
i, is described as an exponential decay through the ice depth, where i
0 = F(C, C
i) is defined as the fraction of the wavelength-integrated incident irradiance transmitted through the top z
i = 0.1 m of the ice, and parameterized according to the sky conditions (C) and sea-ice colour (C
i) (e.g. i
0 = 0.18 (1 − C) + 0.35C for white ice, and i
0 = 0.43(1 − C) + 0.63C for blue ice (Reference Grenfell and MaykutGrenfell and Maykut, 1977; Reference PerovichPerovich, 1996)), and κ
i is the bulk extinction coefficient below z
i varying from 1.1 to 1.5 m−1 (Reference UntersteinerUntersteiner, 1961). Near the surface, κ
i can be one or two orders of magnitude larger than 1.5 m−1 (Reference Grenfell and MaykutGrenfell and Maykut, 1977). Accordingly, a two-layer parameterization scheme for q
i(z, t) is assumed in our ice model. In the top 0.1 m of the ice, we applied 
, 0 < z < z
i. The extinction coefficient κ
1 is calculated as κ
1 = −10 × ln(i
0), i.e. κ
1 is valid for the very uppermost layer by fitting the values of i
0 of Reference Grenfell and MaykutGrenfell and Maykut (1977) observed for the 0.1 m level in the ice. For example, κ
1 = 17 m−1 for clear-sky (C = 0) and white-ice conditions. Below 0.1 m in the ice, we used 
, z ≥ z
i according to Reference Maykut and UntersteinerMaykut and Untersteiner (1971), where κ
2 = 1.5 m−1. A two-layer scheme for q
i (z, t) has been used earlier by Reference SahlbergSahlberg (1988) with a linear profile fitting the i
0 of Reference Grenfell and MaykutGrenfell and Maykut (1977) at 0.1 m in the ice.
More sophisticated expressions for q i(z, t) derived from radiative transfer models have also been suggested and incorporated into ice models. For example, Reference GrenfellGrenfell (1979) applied a radiative transfer model to derive a parameterization for q i (z, t). Reference Brandt and WarrenBrandt and Warren (1993) introduced a downward bulk extinction coefficient κ i (z), which describes the local attenuation rate of absorbed solar radiation, and applied a radiative-transfer model to estimate absorbed solar radiation. The latter has been used by Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999) to calculate the radiation penetration in the Antarctic ice.
In snow, the variation of solar radiation absorbed with depth in snow follows simply the Bouguer–Lambert law, i.e. 
, where the extinction coefficient of snow (κ
s) varies from 5 m−1 for dense snow up to almost 50 m−1 for newly fallen snow, depending on the snow density and grain-size (Reference PerovichPerovich, 1996).
Figure 1 shows the variation of solar radiation in ice and snow, non-dimensionalized by the net solar radiation at the surface, as obtained from various schemes for q(z, t). Examples of this ratio for blue ice estimated by our two-layer scheme and that of Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999) are also included in the figure. One can realize that estimation of ∂q(z, t)/∂z using a finite-difference approach would yield different results using different spatial resolutions.
Fig. 1. Variation of absorbed solar radiation in ice (a) and snow (b) non-dimensionalized by the total net solar radiation at the surface. In (a) the thick line below 0.1 m represents the q(z, t) of Reference Grenfell and MaykutGrenfell and Maykut (1977), and is connected to those of Reference SahlbergSahlberg (1988) (thin line) and Reference Cheng and LauniainenLauniainen and Cheng (1998) (dotted line) within the top 0.1 m ice. The dashed line denotes q(z, t) estimated by Reference GrenfellGrenfell (1979). The lowermost solid line and circles represent q(z, t) for blue ice estimated by the scheme of Reference Launiainen and ChengLauniainen and Cheng (1998) and that given by Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999), respectively. In (b) the snow extinction coefficients used for the various curves are 40 m−1 (dotted line), 25 m−1 (dot-dashed line), 15 m−1 (dashed line) and 5 m−1 (solid line). The circles are model results from Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999).
2.3. Numerical resolution of sea-ice models
The numerical resolution of ice models varies depending on the numerical scheme used and the time-scales of the questions modelled. For example, Reference Maykut and UntersteinerMaykut and Untersteiner (1971) solved the full heat-conduction equation using a primitive form of the Alternating Direction Explicit (ADE) method, i.e. an explicit approximation of the CN scheme. The scheme comprised some 40 gridpoints (10 cm apart) for the vertical resolution together with a time-step of 12 hours. This model was applied to study the equilibrium ice thickness in the Arctic. Reference SemtnerSemtner (1976) considered only three layers in his model with a time-step of 8 hours or so. Due to the coarse resolution, the absorbed solar radiation contributed mainly to the ice-surface heat balance. The source term was thus left out of Equation (1). A simple forward-difference scheme, which significantly reduces the computer resources required, was used. Reference GabisonGabison (1987) constructed a three-layer ice model. The numerical method applied was an implicit CN scheme which satisfies the absolute stable condition and also has only a small (second-order) truncation error. The model time-step was 12 hours. A climatological sea-ice thermodynamic model introduced by Reference Ebert and CurryEbert and Curry (1993) had 10 layers with a time-step of 8 hours. The model followed Reference SemtnerSemtner’s (1976) assumption when calculating in-ice and snow temperatures. The model presented by Reference Flato and BrownFlato and Brown (1996) had a similar overall structure to that of Reference SemtnerSemtner (1976) and was even closer to that of Reference Ebert and CurryEbert and Curry (1993). An implicit CN scheme with a spatial resolution of 10 layers and a time-step of 1 day was used.
It is possible to use a model grid with a fixed number of gridpoints, i.e. a Lagrangian grid (Fig. 2a), as in Reference SemtnerSemtner (1976), Reference Ebert and CurryEbert and Curry (1993), Reference Schramm, Holland, Curry and EbertSchramm and others (1997) and Reference Launiainen and ChengLauniainen and Cheng (1998). Because of the freezing and melting of the ice, the gridpoints in the Lagrangian grid are not at fixed depths. In order to adopt moving gridpoints, Reference SemtnerSemtner (1976) used a linear interpolation of the in-ice temperature from the previous coordinate to the current one at each time-step, subject to the conservation of the total heat content of the ice in his model. However, the temperature dependence of the specific heat of sea ice was ignored. To accommodate the change in grid size in their model, Reference Schramm, Holland, Curry and EbertSchramm and others (1997) employed a procedure assuming the conservation of enthalpy to calculate the interior ice temperature. The volumetric heat capacity was therefore given as a function of ice salinity and temperature. This has also been considered by Reference WintonWinton (2000) in his reformulated version of Semtner’s three-layer model. Reference Cheng and LauniainenCheng and Launiainen (1998) applied an iterative procedure to calculate the in-ice temperature at the model gridpoints for each time-step with a piecewise interpolation of the values from those points given by the previous step. Numerical tests indicated convergence of the result within three steps. Another alternative for the model grid is to maintain a constant inner grid size, i.e. an Eulerian grid (Fig. 2b), as in Reference Maykut and UntersteinerMaykut and Untersteiner (1971). Each inner gridpoint has a fixed coordinate, and only the boundary grid size will vary according to the variation in ice thickness. The total number of gridpoints increases or decreases depending on ice growth or melting.
Fig. 2. Definitions of Lagrangian (a) and Eulerian (b) grid systems with spatial (j) and temporal (k) steps. The black dots are the gridpoints defined by the current step. The short line segment in (a) marks the gridpoints defined by the previous step; the number of gridpoints remains constant. The grid size at each time-step is slightly different (e.g. Δh k−1 , Δh k ). In (b) circles indicate the gridpoints of the previous step, and the black squares are the new gridpoints of the current step. The interior grid size (Δh) remains constant at each time-step. The boundary grid sizes are time-dependent, i.e. Δh sfck−1 , Δh botk−1 .
In this paper an Eulerian grid is employed. The geometric surface variation, due, for example, to snowfall, snowdrift and ice ridging, is not considered for the model run, i.e. the change in the upper boundary grid size is minimized. Thus, for a model run with a specified spatial resolution, the increase in the number of gridpoints and change in the boundary grid size mainly occur at the ice bottom. The model spatial resolution refers to the grid size, not to the number of layers. It is defined as Δh i = H i0/N i, where H i0 is the initial ice thickness and N i is the number of gridpoints. In the simulations we let N i take every integer from 3 to 10, and then even integers from 12 to 32, i.e. a total of 19 spatial resolutions is used. The time-step (τ) we used here is from 600 s up to 6 hours.
3. Simulations and Results
3.1. Model runs compared with analytical solutions
Here the numerical results for simplified ice-model cases are compared with their analytical solutions. The effect of numerical resolution on ice growth, as well as ice and snow temperature profiles, is studied for these simple idealized cases.
During the ice-growth season, Stefan’s law (Reference StefanStefan, 1891) yields the first-order result of ice growth, i.e. 
, where α
2 = 2k/ρL, and 
. Gorres-ponding to an assumption of Stefan’s law (e.g. Reference LeppärantaLeppäranta, 1993), a simplified numerical ice model would be Equation (1) without the source term; a fixed surface temperature T
sfc replaces the surface boundary Equation (2), and F
w = 0 in Equation (3), together with constant ice thermal properties (ρ, c, k and L). Assuming T
sfc = −10°C, T
bot = 0°C, H
i0 = 0.05 m, and using constant ρ, c, k, L (Table 1), we then run both Stefan’s law and the ice model for a 15 day period. The numerical calculations are made using all 19 spatial resolutions, with the initial ice temperature obtained from a linear interpolation between T
sfc and T
bot.
Figure 3 shows the calculated ice-growth rate. The agreement between the analytic solution and the numerical results is very good. The total ice thickness reaches 0.42 m after a 15 day period. With τ = 600 s, the mean error in ice thickness between all 19 averaged model runs and Stefan’s law is only 0. 0034 m. This value increases only slightly to 0.0042 when τ = 6 hours. This procedure has often been used to verify the consistency of a numerical model (e.g. Reference SalorantaSaloranta, 1998). The numerical results showed some oscillations on increasing the time-step (Fig. 3b). This can be improved by using a fully implicit numerical scheme. Since the initial ice formation is very thin, the thermal inertia is small, and the linear ice-temperature profile yielded by Stefan’s law is very close (not shown here) to those given by the numerical simulations, indicating the good accuracy of the numerical results. This suggests that the simple classic linear temperature profile model can be used during the initial ice-growth phase.
Fig. 3. Ice-growth rate calculated by Stefan’s law (dotted line) and the ice model (solid line). There are 19 curves corresponding to the various spatial resolutions within the thickness of the solid line in each sub-plot. The model time-step is 600 s in (a) and 6 hours in (b).
During the ice thermal equilibrium stage, the ice thickness (H i) may be regarded as a constant. If we ignore the source item q(z, t) in Equation (1), the ice model simplifies to ρc∂T/∂t = k∂2 T/∂z, where the thermal properties remain constant. For a given surface temperature (Dirichlet boundary) of a single sinusoidal wave with frequency ω and amplitude T A, i.e. T sfc(0, t) = T bot + T A sin(2πωt), T → T bot for z → ∞, the analytic solution is T(z, t) = T bot + T Ae−(z/d) sin(2πωt − z/d), where d = [(k/(ρcπω)]1/2 (Reference LeppärantaLeppäranta, and others, 1995). For a given surface net flux (Q n) gradient (Neumann boundary) such as k∂T/∂z = Q n(t), and when Q n is further assumed to be of the form Q n = Q A sin(2πωt), T → T bot for z → ∞, the analytic solution becomes T(z, t) = T bot + T A e(−bz) sin(2πωt − bz − π/4), in which the amplitude T A = Q A (2πωρck)(−1/2) and b = [2πωρc/(2k)](1/2) (Reference SavijärviSavijärvi, 1992). Let us assume that H i = 1.2 m, T bot = −5°C, T A = 2°C, ω = 1/86 400 s−1 (diurnal cycle), and constant ice properties; for the Neumann boundary, we take Q A = 33 W m−2 in order to produce a diurnal temperature wave of 2°C amplitude. All 19 spatial resolution runs were made with τ = 600 s. In order to minimize the effect of the initial condition, the ice temperatures at various depths corresponding to different spatial resolutions are read directly from the analytic solution T(z, 0). The model runs are made for a 5 day period. Figure 4 shows comparisons for the vertical ice temperatures between the analytical solutions and the numerical results. It shows that, for both the Dirichlet and Neumann boundaries, the accuracy of the model result is systematically increased if the spatial resolution is improved.
Fig. 4. The vertical ice-temperature profile obtained by an analytic solution (solid line) and the numerical model using spatial resolutions of N i = 3 (dashed line) 10 (dot-dashed line) and 30 (dotted line). The boundary conditions are of the Dirichlet (a) and Neumann (b) type. Each ice-temperature profile corresponds to a time near midday.
During the warm-up season, the term q(z, t) is important and its effect may dominate the temperature profile in the upper ice layer. Assuming that a stationary ice-temperature profile is reached in a short time (∂T(z, t)/∂t = 0), the ice model is simplified to d/dz[kdT(z)/dz − q(z)] = 0, where q(z) = q(0)e−κz, and q(0) =(1− α)Q s, if the lower boundary is given as T = T bot for z = H i (ice thickness).
The analytic solutions for Dirichlet and Neumann boundaries are, respectively:
and
Such solutions illustrate the asymptotic state where the solar radiation tends to drive the temperature profile (Reference LeppärantaLeppäranta, 1995).
For the Dirichlet boundary, the numerical results are quite close to the analytic solution for Ni > 10, similarly to Figure 4a. For the Neumann boundary, the accuracy of the numerical results depends on the magnitude of Q n and q(0) for a given Ni . For example, if Q n is negative and its magnitude ≫ q(0), which means a cold surface (note that Q n represents the sum of surface radiative and turbulent heat fluxes), a large vertical gradient in the ice-temperature profile can be generated. In such a case, the numerical results converge to the analytic solution very well, even for N i = 3 (not shown here). If Q n remains negative, but its magnitude is comparable to that of q(0), this indicates a small vertical gradient of in-ice temperature. Calculations in such a situation for ice and snow temperatures are given in Figure 5 assuming: Q n = [−18, 8] W m−2, q(0) = [30, 10] W m−2, H = [1.2, 0.4] m, T bot = [−2.5, −2]°C and κ = [1.5, 15] m−1, the values in parentheses corresponding to those for ice and snow, respectively. For N i,s = 3, the numerical results were quite inaccurate. This is because the upper boundary condition approximated by Q n = k∂T/∂z ≈ k[T(1) − T sfc]/Δh gives only a crude estimation of the derivative of temperature, where T(1) is the inner ice or snow temperature. Large values of N i,s are necessary for such a weak temperature-gradient solution in order to obtain accurate numerical results. If Q n is positive, it will enhance the effect of q(0) since the latter is always ≥ 0. As a consequence, the calculated ice temperature may tend to rise above freezing point. In such a case, the ice temperature has to remain at T f and the upper ice layer becomes isothermal (dT/dz = 0). The warm-up season turns into a melting season, and the energy absorbed by the ice is used entirely for producing the phase change.
Fig. 5. Vertical temperature profile in ice (a) and snow (b) obtained by the analytic solution (solid line) and a numerical model using as its spatial resolution N i,s = 3 (dashed line), 10 (dot-dashed line) and 30 (dotted line).The boundary condition is of the Neumann type.
We should emphasize that for numerical sea-ice modelling, Q n is actually a complex polynomial of the surface temperature (cf. Equation (2)). Therefore, instead of directly applying a Neumann flux boundary, T sfc is usually solved iteratively from Equation (2) and used as the upper boundary for the numerical scheme of Equation (1) (e.g. Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference GabisonGabison, 1987; Reference Ebert and CurryEbert and Curry, 1993; Reference Launiainen and ChengLauniainen and Cheng, 1998).
3.2. Numerical experiments
The full heat-conduction equation can not be solved analytically. In the following, we present and analyze groups of numerical simulations in terms of various resolutions.
3.2.1. Forcing data
The forcing data are taken from the Baltic Sea ice climate database IDA (Reference HaapalaHaapala and others, 1996). Three winters representing a normal (1983/84), a severe (1986/87) and a mild (1991/92) ice season were selected as typical Baltic winter climate scenarios.
The meteorological data (air temperature T a, wind speed V a, relative humidity Rh, and cloudiness C) observed at Kemi, a meteorological station in the northern Baltic Sea (latitude 65.6° N), were used. Each parameter is simply averaged over the three winters and for each month. If we look at the observed average ice-growth rate from IDA, December, March and April may be referred to as the ice-freezing, ice thermal equilibrium, and ice warm-up seasons, respectively (Reference LauniainenCheng and Launiainen, 1999). The model runs are accordingly made for these three months with fixed meteorological data (Table 2). The solar radiation is incorporated as a monthly average diurnal cycle. The oceanic heat flux is assumed to be 1.0 W m−2, based on field measurements (Reference Shirasawa, Kobinata, Kawamura, Launiainen and VihmaShirasawa and others, 2001).
Table 2. Weather-forcing data for the numerical model runs in section 3.2
For the ice-growth season, the initial ice thickness is taken to be 0.05 m without snowfall. For the ice equilibrium stage, the initial snow and ice thicknesses are 0.3 and 0.6 m, respectively, in accordance with the monthly average from IDA. The spatial resolution of the snow layer is defined as Δh s = H s0/N s, where H s0 is the initial snow thickness; we let N s vary from 3 to 15, while the snow thickness is assumed to remain constant. For the warm-up season, only bare ice with an initial thickness H i0 = 0.7 m is considered.
3.2.2. Ice -growth season
We ignored the diurnal solar radiation in this example. The initial ice temperature is set as a linear interpolation between the air temperature and the freezing point. Figure 6 shows the mean surface temperature and the total ice formation calculated after a 5 day period. With a small Δh i, the model results are sensitive to the value of the time-step τ, this sensitivity rapidly decreasing with increasing Δh i. The larger τ is, the stronger is the impact of Δh i on the results. The sensitivity of the model predictions to Δh i and τ may be explained by the ice mass-balance (ice-growth) adjustment process in response to the external forcing and by the non-linearity effect of ice thermal properties. From the beginning of the model run, the ice layer starts to adapt to its external forcing via phase change. This mass balance is expected to result in a steady evolution process when the external forcing remains constant. The time taken to establish such a steady evolution process depends on the initial ice thickness, the ice thermal properties and the model resolution.
Fig. 6. Modelled average surface temperature (a) and total ice formation (b) during the ice-growth season for a 5 day period vs spatial resolution. Each symbol indicates a given time-step, i.e. 600 s (o) 1 hour (×) 3 hours (+) and 6 hours (*).
We calculated ice-growth rate over a 15 day period. Different Δh i yield results that converge towards each other (Fig. 7a and b). The influence of the time-step τ is studied by inspecting the average ice-growth rate with respect to spatial resolution (Fig. 7c). It appears that the average growth rate converges with increasing spatial resolution and that after decoupling from the initial conditions the grid-size sensitivity decreases.
Fig. 7. Modelled ice-growth rate for all 19 spatial resolution runs with a time-step of (a) 600 s and (b) 6 hours, (c) Average ice-growth rate of all 19 runs for various time-steps: 600 s (solid line), 1 hour (dashed line), 3 hours (dotted line) and 6 hours (dot-dashed line).
The differences between the predictions in Figure 6 are due to the accumulating effects of the ice-mass phase change. Further sensitivity studies indicated that the effect of resolution on model predictions is much less when the ice thickness remains constant. From Figure 6, one may also confirm that the results involving a coarse resolution (both Δh i and τ) are close to those with a fine resolution. This is because, during the freezing season, the solar radiation is small, the heat conduction in the ice is strong and the in-ice temperature appears to be a linear profile. Therefore, the classic linear-temperature ice model may still be valid with reasonable accuracy. Figure 7 indicates that the effect of spatial resolution on the ice-growth rate is significant for a short-time-scale prediction, and in such a case one should pay attention to the model initialization. A large time-step should be avoided under conditions of rapid ice growth. The oscillations of the model results obtained with a large ratio of τ/Δh i can be reduced by adapting the numerical scheme to a fully implicit form.
3.2.3. Thermal equilibrium stage
In this example the model is run for a 5 day period. The initial temperature profiles in the snow and ice are obtained by pre-running the model with the forcing data of March with fixed snow and ice initial thicknesses until a stable temperature profile is achieved. The time-step is 600 s. For comparison, we varied the total extinction coefficient from 5 m−1 for compact snow to 25 m−1 for new snow. The term q i (z, t) in ice is ignored.
The daily surface temperature vs Δh s after 5 days is shown in Figure 8. The daily maximum surface temperature obtained near midday is sensitive to the grid size Δh s. This sensitivity is due to the strong gradient in the absorbed solar radiation near the surface in conjunction with the variation in Δh s. With a large value of Δh s, the absorbed solar radiation tends to contribute more to the surface heat balance and thus leads to a higher surface temperature. The daily minimum surface temperature reached at night is not affected by the daytime solar heating effect, and shows a similar dependency to that in Figure 6a. The diurnal mean surface temperature, however, remains quite stable, indicating its lower sensitivity to spatial resolution. In general, the low-resolution model yields larger diurnal amplitude of surface temperature.
Fig. 8. Modelled daily variation of surface temperature in March (on the 5th day modelled) vs spatial resolution for compact snow (a) and new snow (b). The stars indicate the daily maximum (upper) and minimum (lower) surface temperatures, while the circles give the average daily surface temperature.
Increasing the total extinction coefficient of snow yields similar results (Fig. 8b); however, the maximum surface temperature is less sensitive to spatial resolution. This is because the solar radiation penetrating into the new snow is mostly absorbed in the uppermost layer due to the large extinction coefficient. Under such conditions, the model yields higher maximum surface temperatures compared with Figure 8a, and accordingly produces a larger daily variation.
We analyzed the various surface heat fluxes of Equation (2) modelled as daily means (Fig. 9). The fluxes are quite independent of the grid size, except for the surface conductive heat flux and the surface absorbed shortwave radiation flux. For example, the surface net longwave radiation and turbulent flux (sensible and latent heat) are approximately −25 and 9 W m−2, respectively. The net downward solar radiation at the surface is about 13 W m−2. The magnitudes of the absorbed solar radiation in the surface layer (shortwave radiation contributing to the surface heat balance) and the upward surface conductive heat flux tend to complement each other so that their sum remains quite stable whatever the value of Δh s is. The total net surface flux is not sensitive to changes in spatial resolution.
Fig. 9. Modelled daily average surface heat fluxes vs model spatial resolution using a snow-extinction coefficient of 5 m−1 (a) and 25 m−1 (b). Each symbol refers to a term in the surface heat flux: net longwave radiation flux (*), sensible-heat flux (×), latent-heat flux (+) absorbed solar radiation in the surface layer (o) conductive-heat flux (•) and net surface heat flux (⊕).
For bare ice, the surface temperature and average surface heat fluxes vs spatial resolution using the same weather and initial ice conditions and applying the two-layer scheme of q i(z, t) are shown in Figure 10. The calculated net longwave radiation is about −38 W m−2. The net downward solar radiation at the surface is about 20 W m−2, and the part absorbed by the surface layer together with an upward conductive heat flux in the surface layer approximately balances the outgoing longwave radiation, while the turbulent fluxes are small. The effect of grid size on the absorbed solar radiation in the ice and on the surface conductive heat flux for bare ice is qualitatively similar to but stronger than that on those in snow. Low spatial resolution leads to a large variability in the in-ice temperature in the upper part of the ice due to the greater amount of solar radiation available there for absorption. This suggests that the impact of the spatial resolution on model results is likely to be greatest near the surface because of the strong gradient in the absorbed solar radiation. The sensitivity study indicated that simulations using the q i(z, t) of Reference GrenfellGrenfell (1979) and Reference SahlbergSahlberg (1988), discussed in section 2.2, yield similar results to those obtained using our two-layer scheme for q i(z, t).
3.2.4. Warm-up season
In this case, the air temperature is close to the ice melting point. The initial in-ice temperature is assumed to be a linear interpolation between −8° and −0.5°C. This corresponds to the average night-time steady minimum profile obtained from runs with the March weather data for bare ice. The simulations are made for a 5 day period with a time-step of 600 s. In April, the monthly average diurnal solar radiation is already so large that on the third day of the simulations, surface melting has occurred when using a coarse spatial resolution. The surface temperature and surface heat fluxes vs spatial resolution are shown in Figure 11. The net longwave radiation is −35 W m−2 and the turbulent fluxes of sensible and latent heat are some 5 and −12 W m−2, respectively. The net downward solar radiation at the surface is 67 W m−2. The daily minimum surface temperature is sensitive to the spatial resolution in this case. Because of more solar radiation available to the surface heat balance with a large Δh i, the surface temperature rises and approaches the melting point very quickly so that a downward conductive heat flux is produced. During the night, the heat loss from the surface leads to a low surface temperature. On the other hand, for a small Δh i, the large amount of solar radiation absorbed is considered an internal heat source. This energy is used to heat the subsurface ice rather than to contribute directly to the surface heat balance. Accordingly, an upward below-surface conductive heat flux prevents low minimum surface temperatures during the night by the thermal inertia effect. The subsurface temperature may even increase up to the melting point. The magnitude of solar radiation absorbed in the surface layer is relatively small and may not complement the surface conductive heat flux (Fig. 11b).
Fig. 11. Modelled surface temperature and surface heat fluxes as in Figure 10, but using weather data for April.
Figure 12 shows the in-ice temperature variation on the second day. The subsurface heating can be clearly seen at higher spatial resolutions. Under suitable weather-forcing conditions, a model using a high spatial resolution can yield an in-ice maximum temperature reaching the melting point, thus inducing subsurface melting. This has been found both in field measurements and in numerical sea-ice modelling (Reference Winther, Elvehøy, Bøggild, Sand and ListonWinther and others, 1996; Reference Liston, Winther, Bruland, Elvehøy and SandListon and others, 1999). Our model runs shown in Figure 12 indicate that a high resolution gives a better estimate of the in-ice temperature profile. This suggests that the full numerical model results are consistent with the solutions obtained from the analytical studies. The spatial resolution may affect the location where the ice first starts to melt.
Fig. 12. Modelled in-ice temperature profile using spatial resolutions of (a) 23 cm, (b) 8.75 cm, (c) 3 cm and (d) 2.5 cm. The corresponding values of Ni are: 3, 8, 18 and 28, respectively. In each panel, the three lines indicate the daily minimum (left), daily average (middle) and daily maximum (right) in-ice temperature.
During the melting season, the average air temperature is normally well above the freezing temperature. For example, the average weather data from IDA for May are T a = 4.6°C, V a = 4.1 m s−1, Rh = 70%, C = 0.6, and the daily mean of Q 0 = 310 W m−2. For an initial H i0 = 0.4 m, a melting surface temperature is obtained that then remains unchanged for each spatial resolution applied. An isothermal layer developed from the surface and the subsurface temperature maximum may no longer exist. The in-ice temperature eventually reaches an isothermal stage, i.e. a uniform melting temperature throughout the ice layer. It should be emphasized that the thermal regime of sea ice is affected by various positive feedbacks in response to the warm weather during early spring. For example, surface melting leads to a decrease in surface albedo, and thus more energy is absorbed, resulting in more surface melting. Actual ice melting may therefore be much more complex in response to daily variations in weather forcing.
4. Conclusions
The impact of numerical resolution on the results of a one-dimensional thermodynamic sea-ice model has been studied using average climatic weather-forcing data corresponding to ice-freezing, ice thermal equilibrium and ice warm-up conditions. The numerical results were also compared with traceable analytical solutions for simple model cases. The study focused on the effect of spatial resolution.
During the freezing season, the model yields quite accurate results compared with analytical solutions, indicating that the classical linear ice-temperature model is still valid with reasonable accuracy. The sensitivity of the results to numerical resolution occurs for large values of the ratio of the time-step (s) to the grid size (m). This sensitivity can be decreased by using a fully implicit scheme or by reducing the value of this ratio. Resolution tends more to affect model results for short-term predictions. After adjustment to the initial conditions, the ice growth rate converges with decreasing grid size, and resolution has no significant influence on model predictions.
For the ice-equilibrium stage, comparison with analytical solutions shows that the ice model gives good results for both Dirichlet and Neumann boundary conditions if N i is >10. When the warm-up starts, the downward net solar radiation at the surface becomes large and even comparable with other upward surface fluxes. In this case, the model should have a high spatial resolution in order to obtain good results, especially for a Neumann boundary. For these simplified ideal cases, the accuracy of the result increases as Δh i decreases.
When the strong attenuation of transmitted solar radiation in the surface layer causes the absorbed solar radiation to act effectively as an internal heating source, the calculation of surface temperature and surface fluxes by the full numerical ice model is affected by the grid resolution. With low spatial resolution, the absorbed solar radiation mostly contributes to the surface heat balance. The heat-conduction Equation (1) can be simplified to a parabolic diffusion equation without an internal heating source. For such a condition, the accuracy of the result can be increased by having more layers in the ice or snow. A model with a large Δh i yields a large daily variation in the surface temperature. Results remain sensitive to spatial resolution, especially for fine-resolution cases.
For most thermodynamic ice models (e.g. Reference Maykut and UntersteinerMaykut and Untersteiner, 1971; Reference GabisonGabison, 1987; Reference Ebert and CurryEbert and Curry, 1993), the contribution of the absorbed solar radiation is divided explicitly by assuming a certain portion (0.17−0.35) of the total amount used in Equation (1) and the rest used for the surface heat balance (Equation (2)). This corresponds to a relatively weak effect of penetrating solar radiation q(z, t) inside the ice. On the other hand, however, the grid spacing in these models was >0.1 m, satisfying the physical distribution of q(z, t), i.e. a strong exponential decay near the surface. However, with a grid spacing <0.1 m, the contribution of absorbed solar radiation should be adapted from the partition assumed above. This is particularly important for process studies of subsurface melting of ice or snow. For example, Reference Bøggild, Winther, Sand and ElvehøyBøggild and others (1995) modelled the subsurface melting by implementing an ice model using a resolution of 0.03 m near the surface. The results of Figure 12 implicitly indicate that subsurface melting cannot be modelled with a coarse-spatial-resolution model.
Finally, we should emphasize that our studies correspond to somewhat simplified constant external forcing conditions so that the results on the effect of resolution depend on the validity of such an assumption. Other simplifications are made in connection with the optical properties of snow or ice, such as considering separately the surface albedo and the extinction coefficient. The results presented were obtained using an Eulerian grid, but additional sensitivity studies indicated that similar characteristics can be obtained using a Lagrangian grid. We suggest that for process studies, an ice model should apply a time-step of about 600 s and a spatial resolution of 2−5 cm in ice or snow, if possible. For climatological studies, a relaxation of resolution is certainly needed, but the ratio of time-step to grid size should remain small.
Acknowledgements
J. Launiainen is warmly thanked for his contribution and guidance in developing the ice model. The author is grateful to T.Vihma for fruitful discussions and helpful comments on the manuscript. Suggestions and constructive criticisms by M. Leppäranta, S. Joffre and an anonymous referee produced significant improvements in the manuscript. M. A. Lange was the Scientific Editor. This work was supported by the Baltic Air–Sea–Ice Study of the European Commission funded by contract MAS3-CT97-0117.
