Introduction
To model pack-ice motions properly requires some constitutive law describing the interaction of an area of pack ice with the surrounding ice. In complete simulation models in which ice thickness characteristics are coupled to ice strengths, plastic or non-linear viscous rheologies appear to hold the best promise for proper modelling of ice stress (e.g. Reference Coon, Coon, Maykut, Pritchard, Rothrock and ThorndikeCoon and others, 1974, Hibler, in press), especially near boundaries. However, as a first step in isolating the effect of ice stress on drift, simple linear viscous calculations with bulk and shear viscosities varying in a seasonal manner can yield considerable insight. Moreover, for periods of several days or more, in a reasonably homogeneous region far from shore, stochastic arguments (Reference HiblerHibler, 1977) suggest that a viscous law with the inclusion of both a bulk and shear viscosity may be a reasonable estimate of a time-averaged plastic rheology. It is also noteworthy that a model with both bulk and shear viscosity can behave quite differently from the Newtonian viscous models (with only a shear viscosity) used by Reference CampbellCampbell (1965) and Reference DoroninDoronin (1970).
In order to assess the effects of ice stress in a relatively simple manner, a number of linear viscous calculations with various types of boundary conditions were carried out. These calculations provide a feeling for the sensitivity of drift to the various forcing fields and boundary effects, and give estimates of approximate strengths needed in both viscous and plastic models. A more detailed description of these calculations is available (Reference HiblerHibler and Tucker, 1977). This paper presents some of the more important results from this study pertaining to ice drift far from shore. In particular we discuss here a series of calculations made with periodic boundary conditions and different viscosities. The calculated drift-rates are then compared with observations over a two-year period. The observed drift-rates were obtained from records of one Russian and two U.S. drifting stations, spanning the time period May 1962 to April 1964. To drive the model, daily atmospheric pressure data compiled by the National Center for Atmospheric Research (NCAR) and steady-current data reported by Reference Coachman, Aagaard and HermanCoachman and Aagaard (1974) were used. In the calculations the viscosity parameter is allowed to vary with time but not spatially. Considerable effort is devoted to examining the empirical seasonally-varying viscosity giving the best fit to observations in order to assess the effects of ice stress on drift-rates. In addition seasonally varying viscosity predictions, both with and without currents, are compared to free-drift results and observations. These comparisons help identify the effect of the increased winter ice interaction on the drift and deformation.
Linear-viscous drift model
For drift calculations, we make use of a linear drift theory similar to that used by Reference HiblerHibler (1974). In this model, steady-state equilibrium drift is considered with the following governing equation:
where u is the ice velocity, f the Coriolis parameter, k a unit vector normal to the x, y plane, m the ice mass per unit area, F the force due to variation in internal ice stress, τ w and τ a the water and air stresses on the ice, and G and T the effect of long term geostrophic currents and ocean tilt on the ice motion. The components of current and ice stresses are given by a modified Ekman layer theory (see, for example, Reference RothrockRothrock, 1975).
where ᶲ and θ are Ekman angles in the air and water respectively, U g and V g are the x and y components of the geostrophic wind, U w and V w are the x and y components of the geostrophic ocean flow beneath the Ekman layer, and B and D are constant wind and water stress parameters. In Equation (2), the G components are those components linear in U w and V w. In the case when G is neglected, the sea ice is effectively considered to be moving across a stagnant ocean. The geostrophic current flow is computed by
The force F due to internal ice stress is obtained from the viscous constitutive law proposed by Reference GlenGlen (1970) relating ice stress and strain-rates in the ice
so that
where η and ζ are the shear and bulk viscosities characterizing the ice cover on the geophysical scale. Substituting the various stresses into Equation (1), a system of two coupled second-order linear partial differential equations are obtained for the drift rates ux (x, y) and uy (x, y). For our purposes here, these equations are solved in the infinite boundary case using a Green function generated by transforms. Following Reference HiblerHibler (1974), the procedure is to Fourier transform the coupled equations and solve for the transformed values of ux and uy . These transformed values are then inverted yielding equations of the form
where S j is the jth component of the external stresses due to wind, tilt, and geostrophic currents. Using an input grid of stresses due to wind, currents, and tilt, the convolutions in Equation (4) are approximated by two-dimensional discrete sums which yield periodic rather than infinite boundaries.
The principal justification for the periodic boundary solution is that the response function Gij (x) falls off to small values for large x at a rate which decreases for increasing viscosities. For example, the response function for drift-rate magnitude due to x stress components for viscosity parameters of η = ζ = 4 × 1011 kg/s, falls off to less than 15% of the maximum at a distance of 875 km. Consequently boundary effects on predicted drift-rates far from shore should be minimal.
Input data
The basic input fields for predicting drift and deformation of sea ice are the atmospheric pressure and the dynamic height of the sea surface. For computational purposes all input fields were reduced to a 16 × 16 square grid of 250 km mesh. This spacing was thought to be reasonable taking into consideration the accuracy of the input pressure data and the areal coverage, and in the interest of computer time.
The grid, shown in Figure 1, is oriented with the x axis parallel to long. 150° W. and the y axis parallel to long. 60°W. In this configuration the geographic North Pole conveniently has the grid coordinates (11, 6). This orientation also has the y axis roughly parallel to the north coast of Alaska.
Fig. 1. Square mesh grid used for numerical calculations.
For determination of geostrophic winds, sea-level barometric pressure data were obtained from the National Center for Atmospheric Research. The data were recorded at 12 h intervals and data points were located every five degrees of latitude and longitude. The pressure data were averaged over 8 d intervals and then sampled at 4 d intervals. The smoothed data were converted to the standard 16 × 16 grid using the cubic spline gridding procedure discussed by Reference Davis and KontisDavis and Kontis (1970) and then differentiated using discrete transforms to obtain geostrophic wind estimates. As a check on the accuracy of the pressure data and on the gridding procedure, comparisons were made between similarly averaged Arlis II observations and a pressure interpolated from the gridded data corresponding to the position of Arlis II. Good agreement was obtained except during times of rapid pressure fluctuations.
For the oceanic dynamic height, average values obtained by Reference Coachman, Aagaard and HermanCoachman and Aagaard (1974) were used. Since in this case only one fixed set of data is used it was convenient to prepare a data grid by band using a visual overlay of the dynamic-height contour map.
For comparisons with model predictions we made use of position measurements of three contemporaneous drifting stations T-3, Arlis II, and NP-10 over a 2 year period from May 1962 to April 1964. The original and final positions of the stations are shown in Figure 8. These data were averaged over 8 d intervals to suppress rather large measurement errors, and then resampled every 4 d. (This is the same type of smoothing as applied to the time series for the gridded geostrophic wind).
Fig. 8. Predicted and observed three-station average drift integrated over 60 d intervals.
Results
Empirical Determination of Seasonal Viscosity Variation
In order to estimate temporal variations for the viscosity parameters in the model, an empirical fit between the observed and predicted drift of the three ice stations over a two-year period was carried out. For this purpose predicted x and y drift-rates for each station were generated every four days over a period of 760 d using the eight different values of viscosity listed in Table I. For this fit as well as for subsequent comparisons, the bulk and shear viscosities were taken to be equal. In general, to predict a drift-rate for a given drifting station, the velocity vector is centered on the station location at the beginning of a 4-day interval by interpolating drift rates from the three nearest grid points. Consequently, the predicted driftrates from currents and tilt will vary in time as the position of the stations changes. However, as this variation was found to be small compared with variations in day-to-day geostrophic winds, the geostrophic current and tilt components were neglected in the calculation of the short-term drift-rates used to estimate the best-fit viscosity.
Table I. Numerical parameters used in drift prediction
The standard values of the numerical parameters used in the calculations are shown in Table I. These values are similar to those used by Reference RothrockRothrock (1975). For each viscosity, we calculated the sum of squared differences between the predicted and observed ice velocity components of the centroid of the three drifting stations. The best-fit viscosity was obtained by minimizing the sum of squares over running 60 d intervals.
The viscosity results are shown in Figure 2. Part b of this figure shows the effect on the best-fit viscosity of modifying the wind and water stress coefficients. The most salient characteristic of this figure is the seasonal variation of the best-fit viscosity. Comparison of the best-fit viscosity (standard case) with 60 d smoothed growth rates of 50 cm ice (taken from Reference Thorndike, Thorndike, Rothrock, Maykut and ColonyThorndike and others, 1975) shows a good correlation. The best correlation (0.80, significant at the 0.01% level) occurred with zero lag time (that is, when the viscosity does not lead or lag behind the growth rate). The regression equation for the best-fit viscosity is given by
where ⟨G⟩ is the 60 d averaged growth rate of 50 cm thick ice in units of cm/d and ζ has units of kg s−1. The dashed line in part a of Figure 2 shows the viscosity generated by the regression equation. Since ⟨G⟩ varies in a seasonal manner, Equation (5) effectively yields a “seasonal” viscosity. The best-fit viscosities for different wind and water stress parameters are similar in form to the standard case with the curves being shifted up for higher wind stress (or lower water stress) and down for lower wind stress.
Fig. 2. a. Viscosity values (solid line) giving the best fit between predicted and observed drift-rates far standard values of the wind and water stress parameters B and D. The dashed line is a regression fit of viscosity upon the 60 d averaged growth rate of 50 cm-thick ice (see Equation (5)).
b. The effect on the best-fit viscosity of increasing (curve 1) and decreasing (curve 3) the wind stress coefficient and increasing the water stress coefficient (curve 2).
c. Mean square error between predicted and observed drift-rates for the best-fit viscosity.
Figure 2C shows the error for the best-fit viscosity. This error tends to be of the order of 1 or 2 cm s−1 (≈1 to 2 nautical miles per day). Interestingly enough, the error tends to get very large around September, when an oscillation in the best-fit viscosity can be seen in Figure 2a. Since September more or less marks the end of the summer melt season, this effect may indicate that certain aspects of the model assumptions are breaking down then. One possibility is that the ocean current structure has become modified due to large amounts of open water.
The dominant effect, which seems to account for much of the seasonal viscosity variation exhibited in Figure 2, is a decrease in drift-rate amplitude as the viscosity is increased. This is illustrated in Figure 3, where we show the observed and predicted three-station average velocity for low and high viscosities of 1.0 × 1010 kg s−1 and 4.0 × 1011 kg s−1. Phsically the increase in viscosity increases the capability of the ice pack to transmit stress. When this occurs the drift tends to be driven by the average wind over a large region rather than the local wind. However, it should be noted here that the linear response function (Equation (4)) for predicting the ice drift is not symmetric. Consequently, increasing the ice interaction is not the same as reducing the wind stress. In particular, increasing the viscosity can cause a counter-clockwise turning of the drift direction of up to 20°. This characteristic will be illustrated later.
Fig. 3. Average drift-rate of the three drifting stations. The curves show predicted drift-rates using constant viscosities of 1 × 1010 kg/s (low viscosity) and 4.0 × 1011 kg/s (high viscosity).
A comparison of seasonal viscosity and free-drift results
In the previous section the best fit viscosity was shown to have a pronounced seasonal variation indicating an increased stress transfer through the pack ice in winter. To help assess the character and magnitude of this effect on the drift-rates it is useful to compare the seasonal viscosity predictions with both the free-drift predictions and observed results. For this comparison, drift predictions were made using Equation (5) to generate seasonal viscosities. The free-drift case was estimated by using a constant low viscosity of 1 × 1010 kg s−1. All predictions and observed data were digitally filtered to assess the predictability of the different time-scale components of the drift-rate for the 2 year period. Figure 4 shows 8 d and 25 d smoothed ice velocity components for the ice station NP-10. While neither the free-drift or seasonal-viscosity prediction is obviously superior, the comparisons show that the free-drift case tends to predict excessive velocity during the winter months. A more quantitative comparison of the free-drift and seasonal-viscosity predictions for NP-10 is given in Table II. This table presents the root-mean-square (r.m.s.) differences between these predictions and the observed drift rates (normalized to the observed standard deviation). For a more direct measure of the amplitude of the variance of the predictions, Table III gives the ratio of the predicted to observed variances.
Fig. 4. Predicted and observed drift-rates of NP-10 for seasonally varying viscosity and free drift for 8d and 25 d smoothing.
Table II. R.M.S. differences between predicted and observed NP-10 drift-rates (Normalized to observed standard deviations)
Table III. Variance of predicted NP-10 drift-rates (Normalized to observed variances)
These tables show that the seasonal viscosity is superior to the free drift. The normalized r.m.s. differences are consistently lower, for both x and y components, for the seasonal viscosity predictions. However, the variance ratios indicate that the viscous predictions appear to be overdamped at high frequencies. The free-drift values on the other hand have a greater amplitude but, as Figure 3 and Table II indicate, it is not always the correct amplitude.
The cumulative drift of the three stations illustrates somewhat more graphically the effect of the ice interaction on the drift. Figure 5 shows cumulative drift plotted every 60 d for the seasonal-viscosity and free-drift cases. The variable-viscosity model gives better magnitudes and directions of displacement than does the free-drift model. This turning-angle effect is particularly important since it is an effect that cannot be obtained by simply reducing the forcing field. Consequently, although the addition of an ad hoc seasonally varying free-drift parameter (which controls the magnitude of the response to the forcing fields) might improve the free-drift magnitudes, it would not predict correct directions.
Fig. 5. Predicted and observed cumulative 60 d drift for each drifting station.
Sensitivity of drift predictions to currents
Since the average geostrophic currents are relatively small compared to the winds, it would be expected that currents would not significantly affect short-term velocity predictions. That this is indeed the case is illustrated in Figure 6.
Fig. 6. Effect of currents and tilt on the 25 d smoothed predicted drift-rate of the Arlis II ice island.
In examining cumulative drift over time periods of several years, however, significant differences occur if currents and tilt are neglected. This is because the current effects, while small, are relatively constant. The wind effects, on the other hand, fluctuate. Over a long time interval, the fluctuations largely average out leaving a relatively small mean value. These long-term current effects are illustrated in Figures 7 and 8. In Figure 7, the predicted x and y drift rates for each station with and without geostrophic currents and tilt have been integrated over a 720 d period and compared with the observed final position. In Figure 8, the predicted drift rate for the center of the triangle formed by the three stations was integrated over 60 d intervals and compared to the observed drift. Part b of this figure shows the low viscosity or free-drift result for the same situation. Overall the geostrophic currents account for about 30% of the drift, with the direction being generally toward the Greenland–Spitsbergen passage, while tilt alone accounts for perhaps 15% of the cumulative effect and has an average direction essentially parallel to the Canadian Archipelago.
Fig. 7. Predicted and observed net drift of the three drifting stations with and without current and tilt effects.
Conclusions
The primary objective of this paper was to simulate the effects of ice stress on sea-ice drift empirically by using a linear-viscous model. Secondary goals were to examine the sensitivity of the ice drift to ocean currents and the magnitude of the wind stress. For this purpose a series of calculations were carried out using idealized periodic boundary conditions.
The results from these calculations indicate that sea-ice movement in the central Arctic Ocean can be simulated by a linear-viscous constitutive law and a seasonally varying viscosity. The low ratio of ice velocity to wind stress seen in winter is effected in the model by a high viscosity. The viscous model also naturally modifies the turning angle between the local wind and ice motion in agreement with observation. In summer, the viscosity is low enough that the model reduces nearly to one of free drift. The viscosity was found to correlate well with the growth rate of 50 cm ice, reflecting the physical notion that rapid growth of thin ice in winter will cause thin ice to be less abundant, and the ice cover to be stronger. Steady currents have a negligible effect on short-term drift (weeks) but are important to the long-term drift (years).
The same success in simulation could probably be achieved with a free-drift model with seasonally varying drag coefficients and turning angles. However, the viscous model has the advantage of treating ice stress directly as a term in the momentum equation. It remains to be seen whether non-linear viscous or plastic models can provide better simulations.
Acknowledgements
We wish to thank Dr W. F. Weeks and Mr S. F. Ackley of U.S. Army CRREL for valuable discussions and comments on this work. Comments by an anonymous reference also helped clarify a number of points and improved the final paper. Drift-station data were kindly provided by Dr Alan Thorndike at the University of Washington. This research was supported by the Office of Naval Research and through Geophysical Fluid Dynamics Laboratory/NOAA Grant Number 04-3-022-33.
