Constitutive law

For modeling ice interaction, the ice is considered to be a non-linear viscous compressible fluid obeying the constitutive law

where σ_ij is the two-dimensional stress tensor,

the strain-rate tensor, Ρ is a pressure term (taken as a function of ice-thickness characteristics), and ζ and η are non-linear bulk and shear viscosities. The above equation, with an additional cross tensor term, is known as the Reiner-Rivelin model. The viscosities themselves are functions of the strain-rate invariants. Such a constitutive law enters the momentum balance of sea ice as an internal ice force. For calculations, the dependence of ζ and η on P and Ρ are given by a yield curve in principal stress space. The current study examines some of the effects yield curve shape has on ice flow. Figure 1 shows graphically the four plastic failure criteria used in the present study. Since the system is assumed to be isotropic, the principal stress and strain directions coincide. The normal flow rule states that, when plastic flow occurs (the stress state lies on the yield curve), the strain-rate vector is normal to the yield curve. Of the four different yield curves, three (elliptical, square and cavitating fluid) obey the normal flow rule, whereas the Mohr-Coulomb criterion does not.

Elliptical yield curve

The failure stress for this case lies on an ellipse passing through the origin with a no-stress condition applying for pure divergence. The equations for ζ and η are

with e being the ratio of principal axes of the ellipse. Derivation of the above expression is given by Reference Hibler, III.Hibler (1977). Such a formulation assumes rigid plastic flow together with a normal flow rule. To prevent the values of the viscosities from becoming arbitrarily large (for small strain rates), ζ and η are limited to ceiling values dependent on the ice strength P*. When this occurs, the stress state lies on some smaller ellipse inside the yield surface. The ice cover under such small strain rates is assumed to move as a viscous fluid, or creep — an approximation to the desired rigid behavior. In the present model, the inside ellipse is positioned such that one of the vertices lies on the origin. This is in contrast to a concentric ellipse used in the Reference Hibler,III.Hibler (1979) model. The advantage of this formulation is that, in the absence of external forcing, there will be no motion. The differences in the formulation of the elliptical yield curve are not felt to be critical on the large scale and comparisons show that the two formulations of the elliptical yield curve give essentially the same result.

Fig. 1. Four plastic yield criteria in principal stress space: elliptical (1), cavitation fluid (2), square (3) arid Mohr-Coulomb (4).

The Mohr-Coulomb criterion

The shear strength for this case is taken to be proportional to the compressive stress and the explicit Ρ dependence in the constitutive law is removed. Such a rheology does not have a normal flow rule but, rather, stress states on the two limbs of the Mohr-Coulomb yield curve correspond to pure shear deformation. A cap which limits the maximum allowable compressive stress is introduced to enable ice convergence. The cap strength is parameterized as a function of ice thickness and compactness. The included angle is chosen as 45° for all long-term simulations, which corresponds to an angle of shear resistance, ϕ (in the usual soil mechanics notation), of about 25°. The bulk viscosity, ζ, is zero for pure divergence, and

otherwise. The shear viscosity, η, is given by

where

Cavitating fluid

The cavitating fluid yield curve in principal stress space is that of the Mohr-Coulomb failure criterion with a zero included angle, so no shear stress exists. In this case the normal flow rule is also satisfied. The dependence of ζ is the same as in the Mohr-Coulomb case with α taken to be zero and the shear viscosity, η, also taken to be zero. The cavitating fluid rheology is of interest mainly for its lack of shear stress.

Square rheology

The plastic stress state for this case is confined to four corners of a square. By limiting plastic flow to the vertices, such a yield curve does obey an associative normal flow rule. Figure 2 shows the correspondence between the strain-rate vector and the four failure points. Note that stress-free states correspond to pure divergence. This rheology has the useful physical interpretation that there is a fixed compressive failure stress along any convergent principal axis. For the square yield curve, the value of ζ is not determined explicitly. Instead a bulk viscosity enters Equation (1) as a pressure term Ρ (the first stress invariant) such that

where Ρ is

The dependence of η is given as

where η_c is

ζ_max and η_max are specified values of viscosity. When either of the maximum values are used, the ice is considered to be in a creep state. When this occurs, the stress state lies either on the side of the square or inside. For comparison purposes, the state of pure convergence is located at the maximum extent of the elliptical yield case.

Fig. 2. Strain-rate, vectors and plastic stress state correspondences for square rheology.

Fig. 3. Comparison between computation B-Grid and C-Grid.

Computational grid

The present numerical scheme uses the Arakawa C-Grid geometry instead of the B-Grid of the Reference Hibler,III.Hibler (1979) model. The difference between these grids is illustrated in Figure 3. Here, scalar quantities for both cases lie at the center of the cell. For the B-Grid, the ice velocities are placed at the corners, but for the C-Grid they are located on the sides. The reason for using the C-Grid configuration is two-fold. Firstly, the C-Grid allows for a simpler finite differencing formulation for most of the terms in the internal force calculation. Secondly, by using the C-Grid, an alternating grid point instability in the pressure field is avoided and the total kinetic energy fluctuation with time is also eliminated.

Ice dynamics calculation

The time-stepping numerical scheme replaces the relaxation routine for the momentum time step in the Reference Hibler,III.Hibler (1979) model with an explicit forward Euler time-stepping scheme. A similar technique, using a rheology exhibiting bulk viscosity only, was used by Scmtner (1987). Such a technique allows direct calculation of the stress states and the ice forces once a specific yield curve is given. It is the use of such a method which allows the treatment of more generalized yield curves with case. However, while more general, the procedure is less computationally efficient than any semi-implicit procedure. For the present scheme to be numerically stable, the maximum time step is of the order of seconds its opposed to the one-day time step for theReference Hibler,III.Hibler (1979) model. Figure 4 contrasts the total kinetic energy evolution for both a stable and an unstable case, demonstrating the importance of time-step selection.

Fig. 4. Time-series of total kinetic energy for stable and unstable time step. The two curves are results for 1.5 s and 15.0s time steps using the Mohr-Coulomb rheology. The inset shows the same time-series with an expanded scale so that the oscillatory nature of the unstable solution can be. seen.

* The time-step limit here is shown for an included angle, α, of 45°.

From the time-stepping point of view, Table 1 provides a listing of the time-step stability limit for the four rheologies using the present numerical scheme. Included in the table is the CPU time required for a one-day simulation on an Alliant FX-40 mini-supercomputer. For the present computer code, the computational time required on the Alliant supercomputer is comparable to the time required using a Cyber 205. Note that, whereas the square rheology is least constrained in terms of the time step, it is still computationally inefficient compared to the relaxation method of Reference Hibler,III.Hibler (1979). In the Mohr-Coulomb yield curve case, the matter is further complicated by the time-step limit dependence on the included angle, α. The amount of shear stress allowable in the Mohr-Coulomb criterion is parameterized by the included angle: when α equals zero, the caviating fluid is obtained; when α equals 90°, a portion of the square yield curve is obtained. Figure 5 shows graphically the dependence of time-step limit on the value of α. From the figure, it appears that the time-step limit is reduced dramatically with an increase in the value of the included angle (or amount of shear stress allowed at failure).

Fig. 5. Stability requirement (included angle dependence) for Mohr-Coulomb yield curve.

As indicated, the explicit numerical scheme, although useful, is extremely computationally intensive due to small time-step requirement. The actual CPU time required for a one-year integration varies from 28 h for the square rheology to about 56 h for the Mohr-Coulomb case. For a study of the rheological effects on sea-ice simulations, it is desirable to use other less computationally intensive methods if available. Other more efficient relaxation methods do exist for the elliptical, cavitating fluid and Mohr-Coulomb rheologies. The time required for a one-year integration for these models ranges from 1 to 3h. The three faster models are: model for elliptical case; Flato and Reference Flato and HiblerHibler (1990) cavitating fluid model for cavitating fluid case; and Flato and Reference Flato and HiblerHibler (1990) model for Mohr-Coulomb case (Reference Flato and Hibler, III.Flato and Hibler, 1991). The last method is an extension of the cavitating fluid model with added shear stress in which the finite differencing scheme for handling the shear components is identical to that of the explicit time-stepping method described here. Note that the square rheology cannot be readily implemented with a relaxation scheme. The comparison of direct time-stepping implementation of the various rheologies to the more efficient relaxation schemes is briefly discussed below.

Comparison of relaxation method with explicit time-stepping scheme

Before proceeding to the long-term simulations, some short-term simulations will be discussed to compare the two numerical schemes. One-day simulations with a constant initial thickness field were carried out for the elliptical, cavitating fluid and Mohr-Coulomb rheologies for both the explicit time-stepping procedure and the relaxation methods. The resulting ice-velocity pattern and pressure field were compared and found to be very similar. Furthermore, particle trajectories were also calculated and compared. Here, several initial starting points on the grid were chosen and the particle paths were traced for 300 d using the one-day velocity field result. In general, the magnitude of the vector difference between the daily drift paths calculated using the two different numerical schemes did not exceed 1%. Thus, the use of the more efficient relaxation method instead of current

explicit time-stepping procedure for long-term simulation of elliptical, Mohr-Coulomb and cavitating fluid for the present study is justified.

In addition, for the explicit time-stepping scheme, the one-day simulation results were also used to determine the locations of the stress states in principal space. It was found that about 70–80% of the states for the entire grid lie on the corresponding yield curves. Thus, the explicit time-stepping procedure docs indeed seem to simulate plastic flow.

Boundary and initial conditions

For the long-term simulations discussed below, the models were forced with daily wind, surface air temperature and relative humidity provided by Walsh (personal communication) on a 5° × 5° grid for 1981–82, and mean annual ocean currents and heat flux from the diagnostic ice-ocean calculation of Hibler and Bryan (1987). For the square rheology, the time step used in the simulation is 4.5 s with a single day distorted so that 5 h is taken to be one day. Numerical experiments have been carried out to show that after 5 h the ice velocity has nearly reached the end of day value with the total kinetic energy reaching 99.5% of the end of day value. The time distortion reduced computer run time for one-year integration to about 30 h. The maximum bulk viscosity is reduced by a factor of five from the Hibler (1979) model value to improve the stability of the explicit time-stepping method. One-day time steps are used for all the other models. It is also important to note that all models used in the study have the same thermodynamic treatment as that of Hibler and Walsh (1980), which is essentially that of

Fig. 6. Computation mesh used in simulations.

Parkinson and Washington (Reference Parkinson and Washington1979) with minor modifications.

Velocity and ice-thickness build-up comparison

In order to investigate the effect of different rheologies, the models were applied to the Arctic Basin and integrated for two years. Four separate simulations were carried out, each with a different ice rheology: (a) square, (b) elliptical, (c) Mohr-Coulomb with a cap and (d) cavitating fluid.

Monthly velocity fields were compared for all four cases and found to follow a similar pattern. Figure 7 shows the average monthly velocity pattern for March 1982 for all the different rheologies. It is noted that, for the time period studied, the square and elliptical rheologies produce a velocity field similar to one another as do the Mohr-Coulomb and cavitating fluid rheologies. However, there are differences when the results are compared between the two sets (elliptical and square to Mohr-Coulomb and cavitating fluid). Examining closely Figure 7a–d, it is noted that the velocity patterns for the Mohr-Coulomb and cavitating fluid cases reveal a small gyre near the North Pole which is absent for both the square rheology and the elliptical cases. In addition, the magnitudes of ice velocities for both the elliptical and square cases are found to be lower in the coastal regions especially near the Canadian archipelago. To show the importance of ice interaction, Figures 7e and f show the case of an incompressible Mohr-Coulomb (Mohr-Coulomb failure with no cap) and for free drift (no ice interaction). Notice that the incompressible Mohr-Coulomb rheology produces an entirely different velocity field. On the other hand, free drift produces a velocity pattern similar to the other four rheologies studied, but it allows considerable convergence near the coast in contrast to actual observation. Note also that the ice-velocity magnitude is greater for the free drift case, as expected.

The velocity averaged over each month and over the entire Arctic Basin, and the monthly average outflow through the Fram Strait, are compared next for the four rheologies. The basin velocity and the outflow through the Fram Strait averaged over a two-year period (1981–82) for the different rheologies is listed in Table 2. Again similarities exist for the elliptical and square cases and for the cavitating fluid and Mohr-Coulomb model. It is found that both the cavitating fluid and Mohr-Coulomb cases result in higher velocities than the other two cases. The same holds true for the outflow through the Fram Strait. The smaller velocity and outflow for the elliptical and the square yield curves are believed to be a result of the greater shear stress that these two cases allow.

Spatial patterns of ice-thickness build-up for each month of the two-year integration were also examined for the four cases and were found to be very similar. Consequently no significant effect of rheology could be ascertained. However, total Arctic Basin ice mass (volume) did differ. Figure 8 gives the time series of total ice volume retained in the Arctic Basin for 1982. Here, it is seen that a greater thickness build-up occurs for the square and elliptical rheologies. This result appears to be due to the greater ice outflow through Fram Strait for the Mohr-Coulomb and cavitating fluid rheologies (Table 2).

Fig. 7. Average monthly velocity fields for March 1982: (a) elliptical; (b) square; (c) Mohr-Coulomb; (d) cavitating fluid; (e) incompressible Mohr-Coulomb; (f) free drift. A velocity vector one grid cell long is approximately 0.12 m s⁻¹.

* All results are calculated based on monthly values.

Buoy-drift comparison

Comparison of the overall drift patterns was treated in a previous section. This section will deal with some comparison of observed and simulated buoy drift. Several observed buoy-drift tracks in the Arctic Basin are shown in Figure 9. These buoy-drift data were obtained from World Data Center A for Glaciology and are based on buoy data collected by the Arctic Buoy Program (Colony, personal communication). The lines here connect the monthly observed locations of six individual buoys beginning at June 1981. Figure 10 shows a close-up of two particular buoy tracks along with predicted tracks for the different rheologies involved. Here, the location of the observed buoy was first identified and the corresponding velocity predicted from the model was calculated for a month. The resulting monthly drift vectors were then added vectorally. The model drift shown is thus a cumulative drift track. There are two points to notice in Figure 10. First, the buoy tracks calculated for the elliptical and square rheologies tend to cluster together as do the Mohr-Coulomb and cavitating fluid. Second, the square and elliptical cases exhibit slower drift compared to the other two cases. Such a slow-down may be caused by the larger allowable shear stress present in the former rheologies. On the same figure, predicted buoy drifts for the incompressible Mohr-Coulomb and free drift are shown on the right. Here, the free drift case exhibits the largest drift distance whereas drift in the incompressible case is considerably reduced. The clustering of the elliptical and square cases and of the cavitating Huid and Mohr-Coulomb cases supports the earlier comparison of velocity pattern and ice-thickness build-up.

Fig. 8. Total ice volume retained in the Arctic Basin for 1982.

Fig. 9. Observed buoy-drift tracks in the Arctic Basin. The buoys all begin to drift at June 1981.

Fig. 10. Comparison of observed and simulated buoy drift. The axes are in terms of grid cell coordinates. (a) Buoy no. 3805. (b) Buoy no. 3808.

In terms of statistical analysis, correlation coefficients between monthly observed and computed drift were calculated based on the six buoy tracks studied. These correlation coefficients are shown in Table 3. We see here that the square rheology is the best correlated with the observed buoy drift whereas the free drift case is the worst. However, considering the limited data, the coefficients are too close for all the rheologies to distinguish between them. Table 4 presents some further buoy-drift statistics. In this case the average error radius is the average magnitude of the vector difference between the observed and computed monthly drift. The model-drift ratio to buoy-drift ratio is the ratio of model predicted drift to average monthly buoy-drift distance. Again, the values reported are not significantly different among the rheologies although the square rheology is the best case.

It is believed that a better comparison of rheologies lies in examination of velocity pattern, thickness build-up, and individual model predicted drift tracks. In addition, statistical results and drift comparison based on daily values may be more useful than the monthly average values used in this study.