Comparisons of Sea-Ice Velocity Fields from ERS-1 SAR and a dynamic model

Abstract

Numerical models for sea-ice thickness distribution and velocity are used for ice-dynamics research and ice forecasting. In the modeling work, ERS-1 SAR is an excellent tool, in particular by providing spatial ice-velocity fields as described in the present Baltic Sea study. Ice velocities were extracted from SAR data with 3 and 6 day time intervals using the optical-flow method. A considerable stiffening of the ice pack was observed due to the change in the character of ice déformation under compression from rafting to ridging as the minimum ice thickness increased from 10 to 3O cm. The coastal alignment was strong in the ice motion and the coastal boundary layer width was 20-30 km. An analysis of the SAR data with an ice-dynamics model showed that the observed overall ice-velocity field could be produred using the Hibler viscous-plastic ice rheology. The compressive strength of the ice (over 10 km sc ales) was 2.5x 10⁴ N m⁻² ±50% for ridging and negligible for rafting of very thin ice. The shear strength was significant and the normal yield ellipse aspect ratio of 2 was valid. The 3 day time interval is val id for updating an ice model but for detailed ice-dynamics investigations a data frequency of 1 d⁻¹ or higher would be preferable.

1 Introduction

Models for sea-ice dynamics describe the motion of ice and the evolution of ice-thickness distribution driven by wind and ocean current. The key model element is a formulation how the ice itself reacts to the external forcing fields by redistribution of ice thickness depending on its internal friction. The physics is understood in general but there are still open questions about the details and the parameterization schemes for modeling. Satellite-borne synthetic aperture radars (SAR) are a powerful tool for sea-ice-modeling investigations providing all-weather, high-resolution information on the ice type, concentration, roughness, floe size and velocity. In particular, the ice velocities are valuable data both for qualitative validation of model physics and for quantitative model tuning.

Ice-kinematics algorithms for SAR imagery were widely examined in the 1980s (e.g. Fily and Rothrork, 1987; Vesecky and others, 1988; Kwok and others, 1990). The results were quite good and consequently ice-kinematics products began to serve ice-model development (e.g. Leppäranta and Sun, 1991). ERS-1 SAR provided the first possibility of obtaining satellite SAR time series and, in ice applications, indeed the kinematics data from ERS-1 have been very useful (e.g. Drinkwater and others, 1995; Stern and others, 1995).

The BEERS (Baltic Experiment for ERS-1) program has examined the use of ERS-1 SAR for sea-ice mapping in the Baltic Sea (Leppäranta and others, 1991). The present work is based on data collected during the winter of 1994 when a 3 day repeat cycle was employed for ERS-1 (i.e. images were collected over exactly the same areas at 3 day intervals). This frequency is suitable for the ice-kinematics and modeling investigations. A similar study was made earlier from 1992 ERS-1 SAR data but the results were not as good, because the winter was very mild (Leppäranta and Zhang, 1992); there was only thin ice in die northern Baltic Sea which was frequently destroyed by storms leaving few features lo track and therefore the SAR data could only be utilized for overall changes in ice conditions.

In the winter of 1994, the ice conditions were normal and the SAR data series contains highly interesting ice-dynamics phenomena. Also, in March 1994, a field campaign based from an ice camp was carried out in the Bay of Bothnia, the northernmost basin of the Baltic Sea (Carlslröm, 1994). This paper gives the results of three specific cases of ice cover including an estimation of the ice-displacement fields from the SAR imagery, an analysis of the ice kinematics and an examination of the observed kinematics with a numerical ice-dynamics model of the Hibler type. All the original ice and forcing data are available on the World Wide Web (Haapala and others, 1997; http: //geophysics. helsinki.ri/publ ications/wngl/saLcal. htmf).

2 Sea-Ice Dynamics

2.1. Modeling

Sea-ice mass and motion are described by the thickness distribution Π and velocity u, all are functions of time t and horizontal coordinates x and y. The evolution of these quantities is determined by the conservation laws of momentum and mass, given in general form by, for example, Hibler (1986)

(1)

(2)

where ρ is the ice density, h is the mean ice thickness, d/dt is the material time derivative, f is the Coriolis parameter, k is the unit vector vertically upward, σ is the ice-stress tensor, T _a and T _w are the tangential stresses of air and water on ice, u_g is the surface geostrophic current, Ψ is the mechanical ice-thickness redistributor and Φ is the thermodynamic growth decay term. Bold face type indicates a vector or tensor variable. Equation (1) is obtained by integrating over a gridcell with open water taken as ice with zero thickness. Another formulation is to integrate over ice only and then the wind and water stresses become multiplied by A (e.g. Gray and Morland, 1994). In the present work, the advective terms and the sea-surface slope are neglected, which is a good approximation for the Baltic Sea ire dynamics (Leppäranta, 1981). Because the time-scales considered here are short, the thermodynamic term is neglected.

The mechanical properties of the pack-ire medium show up in the ice-stress divergence or the internal friction and in the ice-thickness redistributor. In general, it is known that

(3)

where Є = 1/2{△u + △u ^T} is the strain-rate tensor, the superscript T standing for transpose. This symmetric tensor has two invariants, normally chosen as Є₁ = Є₁₁ + Є₂₂ and Є₁₁ = {(Є₁₁ — Є₂₂)₂ + 4 Є122}_1/2 equal to the divergence and to ice the maximum rate of shear, respectively. The surface stresses are given by

(4a)

(4b)

where ρ_a, and ρ_w are the air and water densities, C_a and C_w are the air-ice and ice water drag coefficients, u_a and u_w are the air and water velocities, and C_a and C_w are the boundary-layer turning angles in air and water. Equation (4a) assumes that the ice velocity is negligible compared to the wind velocity. In the Baltic Sea, representative values of the drag-law parameters are (Leppäranta and Omstedt, 1990): for the surface wind C_a = 1.8 x 10 ⁻³ and θ_a = θ°, and for the geostrophic current C_w = 3.5 x 10 ⁻³ and θ_w =20°.

In general, given the wind and current fields, the dynamic response of the ice pack lies between two extreme cases: (1) the ice strength is high compared with the forcing and there is no motion, or (2) there is no ire strength and consequently no internal friction which results in so-called “free drift”. Free drift can be expressed as u = u’_a + u_wg where u’ is the wind-driven part of velocity. In the Baltic Sea, as the ire thickness is small, the Coriolis effect is small and the wind-driven part of velocity is about 2.5% ≈ (ρ_aC_a/ ρ_wC_w)^1/2 of the wind speed and directed 20° to the right from the wind I Leppäranta and Omstedt, 1990).

A new Baltic Sea ice model is used to examine the ice velocities derived from ERS-1 SAR (Haapala and Leppäranta, 1996). This model is a seasonal coupled ice- ocean model with a three-level ice-thickness distribution including ice compactness A, undeformed ice thickness h₁, and deformed ice thickness h₁,. The model ice rheology is the Hibler 1979: viscous-plastic law

(5a)

(5b)

where ζ and η are the bulk and shear viscosities, d_ij = 1 (i = j) or 0≠(i j) is the Kronecker delta, P is the ice-strength function, P is the compressive strength of compact ice of unit thickness. C is the compaction constant for strength, ζsub>max is the maximum bulk viscosity, Δ = [Є₄² + Є_Π²/ e² and e is the aspect ratio of the yield ellipse. The parameter ζ_max is introduced for the transition between viscous and plastic stress slates. At ζ_mmax ≤ ?/2Δ, a linear viscous law results, otherwise the flow is plastic. The “standard” rheology parameters are P_* = 10¹² N m ⁻², C = 20, e = 2 and C_max = 10¹²kgs⁻¹ (Hibler, 1979). For these parameters, several sensitivity studies have been made by modeling (e.g. Ip and others, 1991; Holland and others, 1993; Wang and others, 1994). There are few detailed data available of real ice-velocity fields but validation has been based on ice compactness and thickness.

The grid size normally used in the Baltic Sea has been 10 km (5° in latitude and 10° in longitude) with a time step of 30 min. Reduction of the grid size to this level is possible without severely violating the continuum approximation in the Hibler model. There should be many ice floes in a single gridcell; in the Baltic Sea the larges floes may be 10 20 km in diameter which is too great but a more common size of large floes is 1-3 km. These sizes are smaller by a factor of about 5 compared with the Arctic Ocean (e.g. Rothrock and Thoradike, 1984). In general, the park-ice thickness and geometry also scale by this factor and therefore the Baltic Sea 10 km grid would correspond to. 50 km in the Arctic Ocean, which is acceptable there. The original Hibler (1979) model for the Arctic Ocean employed a 100 km grid size.

The model-tuning variables are the drag law parameters for the air and water stresses, which are the dominant external forces, and the ice-rheology parameters. The drag-law parameters are fairly well known but large uncertainties exit in the rheology. The constant P_* sets the ice-stress scale. It is the maximum stress level of an ice sheet which unit thickness may reach; as P_* changes from low to high values the ice drift changes from the free-drill state to a very slow viscous creep. The parameter r is highly important for understanding the stress modes and will be examined below. Results from pack-ice rheology studies have in general shown that the shear strength is less than the compressive strength and thus e > 1; as e becomes very large, a cavitating fluid is approached (Flato and Hibler, 1992). An interesting recent question is whether the Baltic Sea icebreakers affect the geophysical ice strength due to the orientation of ship channels with the coastline, the icebreaking activity might show up in the ratio e.

To study C would require mon¹ accurate ice-compactness data than available here. For example, the ice strength drops by one order of magnitude as ice compactness drops from unity to 1 – 2.3i/C(= 0.88 for C = 20). The linear viscous regime, described by the parameter ζ_max is mainly an aid for numerical solution and it is questionable whether any real physics are involved. The authors are not aware of results based on good observational data about this problem; however, the SAR interference technique looks promising in this respect, providing ice-velocity information on a very fine space scale (Dammert arid others, 1997}. No information can be obtained from (be present experiment, since for ζ _max = 10¹²kgs⁻¹, a 100 km line could creep about 100 m in 3 days which is less than the measurement accuracy. The form of the yield curve would be a very interesting question but it is omitted from the present work. A modeling study around that question has been made by Ip and others (1991).

2.2. Ice velocity from SAR

Sea-ice dynamics can be observed from sequential satellite images in terms of spatial displacement or velocity fields. For the retrieval of ice motion, satellite-borne Synthetic Aperture Radar (SAR) images are of special interest, because of their weather-independence and high resolution. There are several techniques for obtaining sea-ice velocity from consecutive SAR imagery. There are a number of geometric features: ice-floe edges, individual floes, ice ridges, etc., which SARs see and which are not all destroyed as the ice drifts. The swath width is a limitation, as during a 3 day period ice may drift 10-50 km which is up to one-half of the ERS-1 SAR image width (100 km). The location of common image features can be quite different between successive images which causes difficulties in feature tracking.

Manual ice-drift analysis from SAR data is a laborious task. This method best produces the ice-motion field and was also used by Lcppäranla and Sun (1991). However, for statistical work, large datasets have to be processed and for operational work a last method is needed. Consequently, much research has been done to develop automatic methods for ice-drift detection (e.g. Fily and Rothrock, 1987; Veseeky and others, 1988; Daida and others, 1990; Kwok and others, 1990; McConnell and others, 1991).

A very recently developed automatic algorithm is used in this paper (Sun, 1996). The method is based on optical-flow calculation and has advantages in being able to handle simultaneous rotation and deformation, and in reduction of computation time. Moreover, the accuracy of the resulting ice-velocity field can be verified by a statistical method instead of a manual one. This method is a two-step algorithm. The first step is to derive the first-order motion vectors which are related to rigid motion (translation and rotation) and which describe the large-scale displacement of the ice pack. The second step is to obtain the higher-order motion or the deformation which is due to non-rigid motion. Normally, the higher-order motion is one order of magnitude smaller than the rigid displacement; only in the coastal boundary layer does the deformation become large so as to meet the no-slip boundary condition.

The first-order motion information is obtained by a cross-correlation technique. A notable characteristic of the technique is the detection of rotation where the power spectra of two images are transformed from rectangular Cartesian coordinates (w,z) to polar coordinates (r,S), so that the rotation parameter in the (w, z) plane can be converted to the translation parameter in the (r, θ) plane. The axes are cast and north for the (w, z) coordinates and radius and counterclockwise angle from cast for the (r, θ) coordinates. The basic principle is described here; for the details see Sun (1994).

Let P_r and P_T present the power spectra of the reference image and target image, respectively. The two images differ by a rotation angle θ₀. The relation between the power spectra is, respectively, in Cartesian and polar coordinates

(6a)

(6b)

where r = [w² + z²]^1/2- and θ = aretan(z/w). The input is a pair of images with coarse resolution, and the output is the translation and rotation parameters, x₀, y₀ and θ₀ which correspond to the mean drift and rotation of the ice cover. With these parameters,x₀, y₀ and θ₀, the two images can be overlaid so that the common image features are located in similar coordinates. The two colocated images are then used as the input to the second step of the algorithm.

The changes in brightness values over small neighborhoods in both spatial and time domains are used to derive the motion vectors. The raw vectors are then refined step-wise by least-squares estimation with a bilinear transformation model. If the deformation scale of the ice cover is large, a pyramid grid structure is required. Consequently, the common features need to be further colocated at a higher-resolution level according to the refined vectors obtained in the preceding level (i.e. the coarser-resolution level) and the optical-flow calculation and the refinement procedures have to be repeated at each pyramid level. Finally all the motion vectors, from the first, part and all levels of the pyramid, are summed to give the final result of the displacement field.

In the second step, the higher-order motion vectors are derived in a dense grid based on the optical-flow calculation (Fig. 1). The objects imaged at different times consecutively move in space-time. Their trajectories form a continuous flow of brightness which is the so-called optical flow. The orientation of the flow surface is an indicator of object motion. For instance, after time i_p the object α in Figure 1 is steady, so that the surface of the flow is parallel with the t axis. Pure translation corresponds to a surface oblique to the t axis (Fig. 1; object α before time t_p) and rotation is related by a surface whose orientation varies during the motion (Fig. 1, object b.)

Fig. 1. Illustration of the optical flow concept.

Thus, the velocity of the objects can be calculated, based on some assumptions of the optical flow. In this paper, the motion velocity u is derived assuming that the contrast (defined as the spatial gradient of the brightness I = I{x,y,t)) is conserved during the object motion. With reference to Figure 1, this assumption means that total time derivative of the contrast is zero:

(7)

By solving this equation for each point of an image pair, a densely gridded motion field is obtained. The details of the optical-flow calculation can be found in Sun (1996).

With the ice velocity from SAR data, the following key questions can lie examined for, sea-ice models: (1) the physical problem of ice stress, i.e. the rheology problem; (2) the drag-law parameters; and (3) the thickness redistributor, through solving the mass-conservation equation, The thickness and compactness field must also be available from another source.

3 Description of the Selected Cases

3.1. Weather and ice conditions

The winter of 1993-94 was slightly more severe than normal in the Baltic Sea area. The maximum ice extent was about 60% of the whole area of the Baltic Sea, while the long-term average is 45%. Three cases were selected for the present study from the early and middle winter, and they represent 3 and 6 day changes in the ice pack. Figure 2 shows the ice conditions based on the ice reports and charts of the Finnish Institute of Marine Research. In all cases, the basin was almost completely ice-covered but the Ice thickness, increasing during the winter, was different. The accuracy of the ice boundaries depends on how dynamic the situation is, since each ire chart describes the ice information collected during the previous 24 hours or so and die accuracy of the ice thicknesses is estimated as 50%.

Figure 3 presents the wind histories from the regional marine wind estimates of the Finnish Meteorological Institute based on observations and model calculations. There are three wind areas for the Bay of Bothnia and the time interval of the data is 6 hours. The differences in the winds between the three regions were quite small and therefore most of the differential motion must be due to the geometry of the fast-ice boundary and differences in ice characteristics. The accuracy of the wind velocity is at each recording time 2—3 m s⁻¹ and, when averaged over 3 or 6 days, the mean wind is estimated at about 0.5 m s⁻¹ accuracy.

Fig. 2. Ice situations in the Bay of Bothnia in winter 1994: (a) 24 January, (b) 3 February, and (c) 7 March. The rectangles show the location of the SAR images.

Fig. 3. The wind history in the study cases; wind direction is clockwise from north.

Case 1 (25-28 January). On 24 January, the whole Bay of Bothnia was covered by ice (Fig. 2a). The fast-ice thickness was 10-60 cm and the thickness of undeformed pack ice ranged from less than 10 to 50 cm. The western side of the basin was covered by thin (5 cm) new ice. In the central part there was a patch 100 km across of 30-50 cm thick ridged ice. During 25-28 January, the weather was cold: the daily-mean air temperature was 17° to 15° G. Until 27 January, moderate easterly winds prevailed and there was just thermal ice growth; the thin ice grew to 10 cm in thickness. Then northerly winds started and the magnitude increased up to 12 m s⁻¹, forcing the ice to drift. A lead opened up at the northern fast-ice boundary (it froze over rapidly) and the western thin ice was compressed with new ridges forming. The southward ice motion continued until the ERS-1 overflight time on 28 January and further thereafter.

Case 2 (28 January-3 February). This case is a direct continuation of case 1. The daily mean air temperature was 11 to -22 °C. On 28-29 January, the wind still blew from the north with a magnitude of 10-16ms⁻¹. It calmed down on 30 January but then became strong again blowing from the southeast. The ice maps show that in this case the ice thickness increased by about 10 cm and experienced first an overall southwest displacement of 20-30 km and then another smaller displacement westwards. A wide lead opened at the eastern and northern fast-ice boundary and rapidly became tilled with new ice (Fig. 2b). Thus, in this 6 day period there were two separate storm events which makes it more difficult to interpret the modeling comparisons.

Case 3 (5—8 March). The whole Bay of Bothnia was ice-covered and the ice chart of 7 March is well representative for the whole period (Fig. 2c). The ice thickness was 45-80 cm in the fast-ice /one and 30 70 cm in the pack-ice field. The ice pack was considerably stronger and less mobile than in the previous cases, as the minimum thickness was higher. The mean daily air temperature was -3° to -5°C and, since the initial thickness was 30 cm or more, ice growth was insignificant. From 3 March onwards, an intense southern wind developed with a speed more than 12 m s (maximum 18 m s⁻¹) for 5-6 March. Then, the wind ceased but arose again from the southwest. During these storms, the ice situation changed even though the ice pack was initially compact and somewhat thick. Later ice charts show small leads opening in the southern basin and some cracking in the north. The BEERS-94 ice camp in the northern Bay of Bothnia was active during this period and. according to their observations, the ire displacement there was about 8 km north during 2-9 March (Carlström, 1994). No more detailed displacement time series are available for the ice camp but the wind history suggests that the timing of the displacements was during 5-7 March.

3.2. SAR products

In January-March 1991, the 3 day repeat-cycle orbits were employed for ERS-1. For the present study, Fast Delivery (FD) images from descending orbit No. 31 have been used from 25 and 28 January, 3 February, 5 and 8 March. The satellite passage time was 0948 GMT (1148h Finnish time). Each image covers a 250 x 100 km area centered at 64.5° N, 22.5° E, and consists of 2.5 scenes of size 100 x 100 km² (see Fig. 2). Hereafter, the first and second images will be called image 1 and image 2 for each pair of images. All the images have been averaged over 20 x 25 pixels which results in a pixel size of 400 x 400 m². Since the grid used by the ice-dynamics model is about 10 km, it is not necessary to use the full-resolution images.

The images for the first case are displayed in Figure 4a and b (see also the ice map in Figure 2a). On the left side, there is land (the Swedish coast) which is excluded by using a land mask. This mask is made by first taking the difference between the two images, then calculating the local variance of the difference image and finally using a threshold to separate the land area from the moving-ice area. The fast-ice zone also becomes a part of the land mask. To verify the results, the displacement vectors were used to reconstruct image 1 backwards from image 2 in each case (Figure 4c for case 1). The correlation coefficient of the reconstructed image and image 1 (the land area excluded) was 0.57 for case I, 0,53 for case 2 and 0.71 for case 3. The lower values for cases 1 and 2 are due to the fact that the ice was much thinner then, resulting in more destruction of the ice landscape.

Fig. 4. The SAR data for case 1: (a) Image 1 (25 Jan). (b) Image 2 (28 Jan), and (c) reconstruction of image from image 2 using the displacement vectors.

Figure 5 shows the correlogram for test sub-scenes (15 x 15 pixels). Each pair of sub-scenes was matched using the cross-correlation method. Then, the correlation coefficient for each pair of matched sub-scenes was calculated for 0 5 pixel shifts. The correlation nearly vanishes over 3pixel shifts, and from that an error estimate of 2 pixels or O.8 km is obtained for the displacement vectors. The géographie accuracy in the overlaying of consecutive images is estimated as 0.5 km, and the resulting total error becomes 0.9 km by the rms addition. The actual displacements were an order of magnitude larger. In terms of the mean velocity, the accuracy is 0.4 and 0.2 cm s⁻¹ over 3 and 6 days, respectively. Finally, the confidence distributions were produced by calculating the correlation coefficient within a moving window (Fig. 6). A brighter pixel indicates a higher accuracy of the displacement vector at the pixel. Low confidences were found over areas where much ice was destroyed by mechanical deformation. In cases 1 and 2, there are large areas but in case 3, when the ice was thicker and more rigid, the less confident areas are narrow zones at fast-ice boundaries.

4 Comparison with the Sea-Ice Model

4.1. Numerical experiments

The model domain was the whole Gulf of Bothnia (60-66° N); the grid size was 10 km and the time step was 30 min. The boundary conditions were defined with the Fast-ice zone as the area where the sea depth is less than 10 m and setting u there equal to zero. The fast-ice zone was found all around the Gulf of Bothnia except in the south where there was a narrow drill-ice channel out from the Gulf of Bothnia, but in the model domain this was also closed. This channel is located at 60° N, which is significantly far away from the present study area (63 30°—66° N) not to contaminate the model calculations.

The initial ice compactness was set equal to 0.99 in all cases and the ice-thickness fields were initialized on the basis of the ice-chart information. Where a range was given for the ice thickness, the median was taken for the model initial field. The ocean-current velocity was set equal to zero. The regional wind data were linearly interpolated into the ice-model grid.

Fig. 5. Correlogram between the original and reconstructed images.

Based on the initial conditions and wind information, evolution of ice velocity, ice compactness and ice thickness were produced by the model. For the comparisons with SAR products, the modeled velocity time series were averaged over the satellite-pass intervals. Cases 1 and 2 were simulated by a continuous model run beginning on 24January (34 hours before the first SAR image) and continuing until 3 February. Case 3 was simulated from 28 February to 10 March. Three different types of numerical experiments were made: (1) A free-drift model run (P* =0); (2) analysis of the strength constant P*, P* = 10⁴ Nm⁻² (standard), 2.5x10⁴Nm⁻² (high) or P* =5 x l0⁴N m⁻² (very high); and (3) analysis of the aspect ratio e, e = 2 (standard), e = 1.5 (low) or e = 20 (high).The third type was examined only for the March case. The standard values were used for the other model constants, C=20 and ζ _max =10^l2kgs⁻¹. The runs are described in Table 1.

Fig. 6. The confidence images. Left: case 1 (23-28 Jan); middle: case 2 (28 Jan-3 Feb); and right: case 3 (5—8 Mar). The correlation increases with increasing brightness.

The ice strength is proportional to P* and the ice thickness according to Equation (5b). Since the thickness is known to 50%, the runs with the different P* values here should be taken more as sensitivity studies; the chosen P* value is good within the 50%. Below, it will be seen that, within the range 1 x 10⁴ -5 x 10⁴ N m⁻² very large differences appear in the ice-velocity fields; the solution is a highly non-linear function of the ice strength.

Table 1. The varied ice-rheology constants in the numerical experiments

The observed data illustrate the evolution of the icepack mobility during the winter. In cases 1 and 2, the ice is still rather thin and consequently the resulting displacement fields have a large spatial variability. Then, in case 3, the ice pack moves as a rigid body over most of the basin and deformation occurs in narrow zones at the fast-ice boundary. The free-drift velocities reflect the wind-velocity variations.

4.2. Case 1

In case 1, the observed ice velocities were quite high (Fig. 7). The whole basin was ice-covered but the ice at the western side was thin. In the central basin, the magnitude of the ice velocities was close to free-drift but the differential motion was remarkable. The northern side was moving more towards the coast but elsewhere the geometry of the coastline seems to have been strongly forcing the ice flow to the southwest. At the northern fast-ice boundary, the offshore motion increased with distance from the fast-ice boundary across a 20-30 km zone. The largest velocities were about 20 cm s⁻¹.

In the western basin, there was a 20-40 km wide ice-boundary zone which showed strong deformation, and at the boundary between thin and thick ice there was a sharp velocity change. The minimum ice thickness was less than 10 cm and, when compressed, such ice undergoes rafting rather than ridging. The P*≠ 0 cases give too low velocities in the middle of the basin and the whole ice pack is far too stiff in terms of the differential ice motion (Figs 7 and 8). For P* = 10⁴ N m⁻², the modeled and observed velocity fields are close in the coastal zone. Simulations with large P* produce velocities that are 50% lower than observed. The conclusion is that frictional losses due to ice deformation were small, i.e. internal friction of the ice was small and therefore the resulting velocities in the main ice pack were practically in the free-drift mode. The coastline effect is, however, notable as the free-drift direction is biased.

Fig. 7. Mean velocity fields for case 1 (25-28 Jan): (a) from ERS-1 SAR, (b) free drift, (c) standard strength (P* = 2.5 x 10⁴ N m⁻², (d) high strength (P* —2.5 x 10⁴N m⁻²).The points where time series wen taken are shown in (a).

The overall ice compactness and thickness evolution is similar in all simulations. But there are large differences in the extent of open water and ridging. With no resistance to deformation, the free-drift simulation over estimates the deformation near the coast. The P* ≠ 0 runs also produced large ice-thickness changes in the western basin (Fig. 9). In the standard strength run, the thicknesses increased up to a maximum of 100 cm from an initial 10 cm, which is abnormally high. However, the ice-chart data on ice thickness are not good for deformed ice and they cannot be compared with the modeled ice-thickness evolution.

Fig. 8. The longitudinal and transverse velocity components in case 1 along 65° Nas observed by ERS-1 SAR and as produced by the model runs 1—4 (see Table 1).

Fig. 9. The thickness fields in case 1 ; (a) initial, (b) at the end of run 2 (standard), and (c) at the end of run 3 (high strength). Units are in cm.

For the ice-velocity field observed by the SAR, the displacement divergence shows values down to -1. If we assume a constant rate, then the mean ice thickness must have increased by a factor of e or almost three-fold (the compactness was initially 0.99). As the ice was originally thin, the result must have been 2.5 times thicker rafted ice or rubble fields rather than ridges, which is thinner than predicted by the model (Fig. 9). The eastern and northern side of the SAR swath experienced opening with the displacement divergence up to one which means that the ice compactness decreased to 0.99 — 1/e ≈ 0.62. There was intensive cracking and lead opening and, furthermore, due to the low air temperature, rapid freezing and new ice growth. The ice charts suggest that openings formed here froze up to 10 cm in thickness during this period; but, this thickness seems to be too small to stiffen significantly the ice pack. The 3 day rotation of the ice field was al largest 0.5 rad, clockwise (negative) in the north and anticlockwise (positive) in the south. The gyres were 30 - 40 km m size. The main cause of the northern gyre was the land-boundary condition but, for the southern opposite gyre, there was an ice-drift speed maximum in the coastal boundary layer and consequently an anticlockwise gyre. These gyres do not appear well in the model output. Probably, a more detailed fast-ice boundary configuration and initial ice-thickness field would have led to better agreement. Also, the passive ocean treatment here may be a reason for lack of the gyres in the model.

The time series of the model ice velocities at 64°N, 22°E and 65° N, 23° E are shown in Figure 10 (see also Figure 7). Clearly, the ice followed the wind all the time. In weak winds, the high ice-strength simulation drops to a creep stale but in strong winds the difference between the runs becomes smaller — a feature of the non-linear mechanical behavior of pack ice.

4.3. Case 2

Case 2 was 6 days long with variable winds and the observed net displacement was southward, at a maximum speed of about 5 cm s⁻¹ (Kg. 11). This was much less than for case 1. The small total displacement was dur to the variability of winds which caused motion back and forth. The real vector field shows a much different motion as compared m the model. Again, the sharpest velocity gradients are at the boundary of old and new ice and thus connected with ice compactness and thickness fields. The offshore motion in the north has a weak spatial gradient.

Fig. 10. The ice-velocity time series for selected points from cases 1 and 2 (25 Jan-3Feb) in the model runs 1-4 (see Table 1). The point are 64°N, 22° E (top) and65°N,23°E (below).

The free-drift solution (Fig. l1b) is non-uniform within the basin and suggests ice drift mainly towards the northwest, very different from the observed drift. In variable wind conditions, a large discrepancy may result due to averaging effects, because ice drift is a highly non-linear function of wind while the free-drift model is linear (except at the unrealistically high ice-accumulation zones). In this case, there was much more resistance to die ice motion to the west than to the south. The simulations with P* ≠ 0 give better agreement but still the general pattern looks different (Fig. la and d). All plastic solutions show in the north an anticlockwise gyre and strong coastal drift in the west which are absent from the SAR data. Because of the large variability of the wind history, it is difficult to say what is the best choice for die strength constant but P* = 1.0-2.5 x 10⁴ N m⁻² is reasonable. Increasing P* would direct the motion more southward but then also decrease the magnitude of velocity ton much. In fact, the theology becomes crude when examined in detail as here, and different P* values might even be needed for the westward and southward displacements.

Differences in die simulations are well illustrated by the time series I Fig. 10). In the north, the results were similar because of the low ice compactness. In the south, the compactness was close to 1 and the ice becomes more immobile with increasing P*. The opening in the north and ridging in the south continued. The maximum thicknesses are: 80-100 cm for P* = 10⁴ N m ⁻², 60-80 cm for P* = 25 x l0⁴ N m ⁻² and 45-60 cm for P* = 5 x 10⁴ N m ⁻². The ice charts do not show ice-thickness build-up, because of the dynamics and therefore a comparison does not tell much (Fig. 2b).

The strain magnitudes were smaller than in case 1. In this period, there were temporal variations in the ice drift as discussed above and therefore the total strains were small. Because ice ridging and rubble formation are irreversible, the production of deformed ice may be seriously underestimated if it is based on averaging of temporally highly variable ice dynamics.

4.4. Case 3

Case 3 was 1 month later than case 2, when the ice had grown much thicker. The observations show a quite stiff northward displacement (Figs 12 and 13). The minimum pack-tee thickness was 30cm except for a small spot in the south. The displacements were mostly around 10km; in the north, the SAR product agrees well with the observed ice-camp displacement of 8 kits (Carlström, 1994). The drift was uniform over most of the ice Held and the deformation /one was only about one gridcell (10km) wide in the north. The behavior of this thick ice pack is characteristic of a plastic medium under constant forcing with deformation focused in narrow slip lines.

The modeled ice-drift directions are consistent in all simulations with the observations, since the period was governed by a rather uniform southerly wind. However, the ice-velocity magnitudes are very different. The free drift (Fig. 12b and the standard ice strength (P* = 10⁴N m⁻²) result in serious overestimation of the ice velocity; they give magnitudes of about 20 and 10 cm s ⁻¹, respectively, while the observed values are only around 3 cm s⁻¹ (Fig. 12a). Note that the free-drill solution shows odd features at the fast-ice boundary in the north. Because of no resistance to deformation, much ice accumulates al the boundary from the free-drift, which consequently results in an unrealistically high Coriolis force.

Even the high ice-strength [P* = 2.5 x 10’ X m ²) simulations result in too large motions (Fig. 12c) and a very high value (P* = 5x10⁴ N m⁻²) is needed to explain the case (Fig. 12d). Time series of the model ice velocities at 65 °N, 23° E are shown in Figure 14. Comparing with information from the field base (Carlström, 1994), the very high ice-strength solution seems to be the most realistic, although it slightly underestimates the velocities. In this model run, the ice was mobile for l day with maximum drift speed about 10 cm s⁻¹and stationary for the rest of the time. The ice-thickness changes (Fig. 15) show clearly how in the high P^A run the thicknesses increase up to three-fold, while the very high P* run gives small changes -25% at most. This thick ice-field transforms to thick ridges in compression, leaving spots between them unaffected, the three-folding concerning just spatial averages over gridcells. The field group reported ridging but not to as large an extent as modeled (Carlström, 1994).

Fig. 11. Mean velocity fields for case 2 (28 Jan-3Feb): (a) from ERS-1 SÂR. (b) free drift, (c) standard strength (P* = 10⁴ Nm⁻²). (d) high strength (P* =2.5x10⁴ Nm⁻²).The points where time series were taken are shown in (a).

The shear-strength parameter e was examined by using the high overall strength level (P* = 2.5 x 10⁴ N m⁻²) and performing the simulation with low (e = 1.5) and high (e = 20) values (Fig. 12e and f) the choice of e affects the differential motion. For large e, as the shear strength disappears, more spatial-drift variations appear and the flow can follow the coastline geometry more smoothly. The model outcome becomes qualitatively different from the observed, almost rigid displacement, and the existence of a significant shear strength is concluded. Thus, the cavitating fluid approach (Flato and Hibler, 1992) is not good in the Baltic-Sea when the ice has grown to 30 cm or more in thickness. The low-e case results are close to the standard but in a run with e = 1 the motion was in very slow viscous creep during the whole period. The data are too few to establish a single value for e but it is in the range 1.5-4; the normally used value e = 2 is therefore good.

The model produces significant ice motion only when the yield strength is achieved. A one-dimensional plastic model has a simple yield criterion: T_aL > P*h for compact ice cover. In the present ease, L ≈250 km and the product T_aL is about 5xl0⁴N m⁻² for U_a and so 0.5 m thick ice would be stationary if P* were 10⁵ N m⁻² or more. There is observational evidence that such strengths may occur (Zhang and Leppäranta, 1995) but in the present cases the strengths were lower.

Fig. 12. Mean velocity fields for case 3 (5-8 Mar): (a) from ERS-1 SAR, (b) free drift, (c) high strength (P* = 2.5x10⁴N m⁻²), (d) very high strength (P* =5 x 10⁴ N m⁻², (e) high shear strength (e = 1.5). and (f) low shear strength (e = 20). The points where time series were token are shown in (a).

Fig. 13. The longitudinal and transverse velocity components in case 3 along 65°N as observed by ERS-1 SAR and as produced by the model runs 1-6 (see Table 1)

Fig. 14. The ice-velocity time series for one selected point (65°N, 23°E) from case 3 (3-5 Mar) in the model runs 1-6 (see Table 1).

Fig. 15. The thickness fields in case 3: (a) initial, (b) at the end of run 3 (high strength hand (c) at the end of run 4 (very high strength). Units are in cm.

The mechanical deformation of the ice pack, as observed by the SAR, shows notable values only at boundaries, since the ice pack moved almost as one block. With this stiff north-ward ice flow, the absolute value of the divergence remained small. The southern ice held experienced 10%, openings which were also seen in the ice charts. The west and north sides converged but the convergence was less than 5%. This corresponds to the formation of new ridges and again the mechanical ice-thickness production is much less than in the model output in Figure 15. The 5% convergence corresponds here to an increase of about 1.5 cm in the mean ice thickness. Converting to normal Baltic Sea ridge statistics (Leppäranta, 1981), this means that the number of ridges increased by about 0.5- 1 km⁻¹. The rotation field was also simple, showing 0.1 rad anticlockwise turning al 64.5°N.

5 Conclusions

An analysis of ice displacements in the Bay of Bothnia, Baltic Sea, has been made based on ERS-I SAR imagery and a sea-ice dynamics model. The basin was almost completely covered by ice with the ice thickness increasing as time advanced. During the study period, the 3 day repeat-cycle orbits were employed for the ERS-I SAR. An excellent ice-data- time series was produced by the satellites and the ice velocities were extracted from the SAR data using the optical-flow method. Three cases, two for 3 day and one for 6 day displacements, were selected for the present study.

A new automatic algorithm has recently been developed to obtain the ice velocity from SAR data (Sun, 1996) and is used in this work. This algorithm is based on the optical-flow method and has advantages in being able to handle simultaneous rotation and deformation, and in the reduction of computation time. This method is a two-step algorithm: the first step is to derive the first-order motion vectors (rigid translation and rotation) and the second step is to obtain the deformation. In the present dataset, the geographic accuracy in the overlaying of consecutive images was estimated as 0.5 km and the total error in displacements was 0.9 km. The actual displacements were an order of magnitude larger.

The observed ice velocities showed a considerable stiffening of the ice pack as the minimum ice thickness increased from 10 to 30cm.This is due to the change in the character of ice deformation under compression from rafting to ridging. Thin-ice compressions (rafting) were one order of magnitude larger than thick-ice compressions (ridging). The study eases showed openings and closings of tens of per cents during 3 6day time periods. Such large deformations highly affect the ice-mass budget of the basin, creating thicker ice areas due to rafting and ridging and opening leads which are potential areas for rapid new-ice production. The observed data snow that, for highly compact ice fields, the coastal alignment of ice drift is strong and dominates over the Ekman mining angle.

An analysis was made- of differences between the ice velocities from the SAR data and from an ice-dynamics model for the Baltic Sea. This model uses the viscous-plastic ice rheology of Hibler (1979). Comparisons showed that the observed ice-velocity field could be produced with the model. In all cases, the ice field experienced heavy compression. The results supported the assumption of a plastic rheology for thick (more than 30em and compact ice. The rheology parameters were examined through several model experiments and the result was as follows: the strength constant P* (equal to the compressive strength of compact ice of unit thickness) is best represented by the value 2.5 xl0⁴ N m⁻² ±50%. The resulting velocity Held was sensitive to P’ and reduction to 1 xl0⁴ N m ⁻² was close to free drift but, with doubling of P* the motion dropped remarkably and would for even larger P* soon become a slow creep state. Also, the model experiments suggested that the best P* varied in different ice conditions. A preponderance of thin ice or leads produces almost free-drift motion (P* =0), while for thick compact mid-winter ice the best value for P* was 5 x 10⁴ N m⁻². Therefore, the direct proportionality of the strength of compact ice to mean ice thickness needs modification; the present cases are, however, not a large enough dataset to examine this question further.

The parameter e (the ratio of compressive strength to shear strength) was found to be probably within the range 1.5-4. The parameter C (1/C is the e-folding of ice-strength dependency of ice compactness) could not be studied in detail, because of the lack of accurate ice-compactness data, hut the model results proved that the sensitivity of the strength to compactness is high, which means that C »l and C = 20 is reasonable. The fourth parameter of the Hibler (1979) rheology is the maximum viscosity, which is mainly an aid for numerical solution: to examine whether there are any physics involved, would need much better displacement accuracy than here. An elliptic yield curve was employed in our model, as it was in the original Hibler model, but no comparisons were made with other yield curves. Such a study would, however, be very useful.

ERS-1 SAR is an excellent tool in sea-ice modeling work, in particular by providing spatial ice-velocity fields. The repeat cycle of 3 days is good for updating an ice model but for detailed ice-dynamics investigations a data frequency of once per day would be preferable. Over very short time-scales, inaccuracy is important and also instantaneous velocity fields are less informative, and over very long time-scales feature identification and result interpretation become difficult. The cases presented here serve as a good validation test for sea-ice model development in the future.