2(a). Momentum balance

The momentum balance includes inertial terms, Coriolis force, wind and water stresses and, most important for this paper, ice interaction. In Cartesian coordinates the momentum balance is

where D/Dt = ∂/∂t + u·Δ the substantial time derivative, is a unit vector normal to the surface, k is the ice velocity, f is the Coriolis parameter, m is the ice mass per unit area, T_a and T_y are forces due to air and water stresses, H is the dynamic height of the sea surface, g is the acceleration due to gravity, and I is the force due to variation in internal ice stress. The air and water stresses are determined by integral boundary layers (see e.g. Reference McPhee and PritchardMcPhee 1980), assuming constant turning angles:

and

where U_g is the geostrophic wind, U_w is the geo-strophic ocean current, ф and θ are air and water turning angles, C_a and C_w are air and water drag coefficients, and p_a and p₅ are air and water densities.

To simplify analysis in these idealized simu-lations, U_w and H have been set equal to zero and the term for momentum advection ignored, (It should be noted, however, that, depending on the location, these current terms may be important in ice-forecasting applications.) As a consequence the momentum balance employed here is

For modeling the ice interaction the ice is con-sidered to ;

obey a nonlinear viscous constitutive law given by

where oij is the two-dimensional stress tensor, ε ij the straln-rate tensor, P/2 a pressure term equal toone half of the ice strength P (taken to be a functionof the ice-thickness characteristics as discussedlater), and ~ and n are nonlinear bulk andshear viscosities. The internal ice stress I isrelated to this constitutive law by

For calculations performed here the dependence of ç and n on εij and P is taken so that the stress state lies on an elliptical yield curve passing through the origin with no-stress conditions applying for pure divergence (Reference HiblerHibler 1979).

2(b). Compactness evolution and ice strength

The main prognostic parameter in the model is the compactness A. Since for forecast purposes we are only concerned with changes of A over, at most, a few days, we neglect thermodynamic terms and consider only the evolution of the ice compactness and thickness due to mechanical effects:

and

where A is the fraction of the area covered by ice, and h is the mean ice thickness averaged over a whole grid cell. Inclusion of thermodynamic terms can make some difference for forecasts longer than a few days. However, within the approximation of the two-level model, for short-term forecasts (<36 h), errors due to the neglect of thermodynamics are felt to be about the same order of magnitude as errors in the initialization of the compactness fields.

We further specify that A<1, which essentially means that once the open water is removed, further convergence of the ice takes place by ice thickening. This is equivalent to adding a term for mechanical redistribution to Equation (4). To relate the ice strength to the compactness (and ice thickness) the following relation is used:

Where P* and C are constants Note that by allowing h and A to go zero, a free ice edge may be treated.

4(a). Analytic analysis of the momentum balance for a one-dimensional system

Consider the one-dimensional case of the momentum balance employing only a linear water drag term Cu, an external constant wind stress F, and a one-dimensional ice interaction stress c:

Consider the region bounded by no slip walls at x = 0 and x = L. For our plastic rheology we will assume that o = -P for ?u/? x<0 and σ = 0 for ?u/?x>0. These assumptions define a rigid plastic rheology with no tensile strength. A solution for this case may be constructed by noticing (i) that for any converging deformation σ = -P, (ii) for no deformation 0> σ >-P, and (iii) the maximum force would be expected to take place at the right-hand boundary since the wind stress has built up at this point. With these considerations in mind we see first that if F L<P no motion of any kind will occur and the system will be rigid. If on the other hand F L>P then a solution of Equation (7) can occur for a σ satisfying the above plastic assumptions and σ(L) = -P. Integrating Equation (7) from x to L we have

However, we also must have ?u/?x = 0 except at the boundary, otherwise σ (X) = -P and Equation (7) cannot be satisfied. Therefore u is a constant, as is F, and we obtain

The value of u0 may be obtained by noting that the total force acting on the moving rigid block is -Cu0L + FL - P, and since there is no acceleration

Another way to see this is to note that σ (x=0) since divergence is occurring there and we nave assumed no tensile strength. The validity of this solution was also verified numerically.

4(b). Analysis of coupled behavior for an idealized one-dimensional system

As in the momentum balance only, consider a one-dimensional system having a linear water drag, rigid plastic interaction and constant wind stress F from left to right. The initial condition consists of a compactness A of 0.9. For the coupling of P to A let us consider P discontinuously changing at A = 1 as follows

where P>P1 Because of this discontinuous change the local momentum balance will not be modified by any changes in compactness until the compactness goes to 1. Based on the analysis of the momentum balance only, presented earlier, all deformation takes place at a right-hand boundary; consequently, at some arbitrary time At after the initiation of motion we would expect the geometry in Figure 1 to apply. By the earlier analysis of momentum balance, for (P-P ₁ < FL₃ all velocities in the L3 region (compactness of 1) will be zero. Also by this analysis, the velocities in the L2 region will be given by

so that as L₂ becomes smaller U_gwill slow down unless P₁ = 0. Since, with the L₃ region motionless, all the deformation is occurring at the L₂, L₃ boundary it is clear that the L_g, L₃ boundary will propagate outward (to the left) at a speed of U₂2/0.1. However, note that for L₃ large enough (i.e. 1₃(P-P₁)/F) the L3 region will start moving again with a speed

Fig. 1.

Geometry of idealized compactness conditions.

When this occurs, the propagation of the ridging front will slow down to [U_g-U₃)/0,l. Also, once all the ice is converted to 100% compactness, the speed will then slowly decrease as the length L₃ decreases due to removal of ice by ridging at the right-hand boundary. Based on these considerations one can construct a qualitative speed versus time of two arbitrary points, say p and p', where p' is located further left than p (Fig.2.)

Fig. 2.

Qualitative plot of velocity versus time for two points in example of idealized build-up.

Obviously by varying strengths and/or the wind forcing the various slopes, magnitudes and spacing in time of these profiles can be modified. One of the more interesting features of the system is that due to the complexities of the thickness coupling, for appropriate strength differences, points near the right-hand coast will actually stop for a while and then speed up, even though external forcing is constant. The essential physics here is that, because of the discontinuous nature of the ridging front and the nonlinear rheology, points in the L₃ region will not move at all until 1I3 becomes large enough to accumulate an adequate wind fetch to overcome the differential of plastic stresses at the boundaries.

4(c). Numerical investigation of response characteristics of a two-dimensional system

A series of two-dimensional numerical simulations yielded results that generally agree with the one-dimensional analysis. However, there are signify-cant differences associated with the fact that P does not change discontinuously but is instead a smooth function of A. For the standard case study, the coupled equations were integrated for 18 h at time steps of 15 min on a 9 × 9 grid. The initial conditions consisted of 10% open water together with a mean ice thickness of 0.5 m. To simplify analysis the ice strength is taken to depend only on the compactness, not on the ice thickness. As a consequence the ice strength in these particular idealized studies is given by

To make the test comparable to the one-dimensional analysis, the turning angles and the Coriolis parameter were set equal to zero. A constant wind speed of 9.23 m s”1 in the positive x direction was used. Other parameters are the same as in Table I. Tests using different values of C and P* (as well as a linear dependence of P upon A) are discussed in Hibler and others (in press).

The basic characteristics of this idealized numerical simulation are shown in Figures 3 and 4, which show x velocity components versus time for several grid points, together with velocity and compactness profiles for selected times. Both these plots were taken from grid points and grid cells centered in the y direction. The general behavior of the velocity time series is commensurate with the idealized one-dimensional analysis, with the points nearest the coast rapidly decreasing in speed with a later increase as the region of ridging moves further out. However, there are significant differences from the analytic case associated with the fact that the motion near the coast is balanced by a gradient in compactness (see Fig.4.(a)) rather than a discontinuous change. This gradient means that (i) the points do not ever totally stop (since some convergence must occur everywhere to allow the gradient to move outward) and (ii) kinematic waves (see Fig.3.(a)) are allowed to propagate upstream (or more precisely down the stress gradient).

Fig. 3.

Time series of x component of ice velocity at grid points progressively further from right-hand boundary of build-up experiment. The grid points used were centered in the y direction and the distances from the right-hand boundary are labeled.

Fig. 4.

(a) The profit Velocity and (b) compactness profiles at different times in the numerical experiments on build-up. Fîles were taken from grid cells centered in the y direction.

These waves arise from the fact that the velocity of ice is balanced in part by the stress gradient while the stress gradient is maintained by a convergence of the ice velocity field. As a consequence, a perturbation of either one of these quantities can cause a wave to propagate down the stress gradient slope. This wave can naturally develop in the evolution of the system since the initial build-up can be out of equilibrium with regard to kinematic waves (when the dependence of strength on compactness is nonlinear). Such a kinematic wave is apparent in Figure 4 mostly notably at hour 6 in the velocity time series of the point nearest the coast. Tests at smaller time steps verified that this wave was not a numerical artifact. However, while this wave is noticeable in the velocity time series, its effect on the compactness profiles is not major (see Fig.4.(b)). Closer examinations of this wave show it to slow down and diffuse as it progresses outward. While beyond the scope of this paper it should be noted that both these features can be explained using a one-dimensional analytical model.

An important ramification of compactness gradient is that there will also be a spatial velocity gradient (see Fig.4.(b)) with smaller velocities closer to the coast. This velocity gradient will, however, gradually decrease as the ice becomes more compact due to convergence. Observations in the Bay of Bothnia (see e.g. Omstedt and Sahlberg 1977, Leppäranta 1980) indicate that such velocity gradients are observed, and hence are a useful feature of this forecasting model. Specifically, observations from three different ships during a joint Swedish-Finnish field experiment in 1977 (Omstedt and Sahlberg 1977) showed that under onshore winds the drift of ships nearest the coast decreased more rapidly than did the drift of ships further out. Simulations with a modification of this model representing the Bay of Bothnia reproduced the main features of this observation with the ice build-up playing a significant role in the drift reduction (Hibler and others in press).