1. Introduction
Pack-ice growth and decay are controlled by both the dynamic and thermodynamic characteristics of sea ice. Owing to the unconstrained nature of the Antarctic sea-ice cover, it is likely that dynamic effects play a dominant role in its behavior. The Weddell Sea pack ice composes an important portion of the Antarctic ice cover and has both representative and unique characteristics. The eastern portion tends to be largely unconstrained, with a cycle of advance and retreat that is representative of much of the remainder of the Antarctic pack. The western portion, on the other hand, contains a region of perennial sea ice which forms a base for advective advance of the pack in the fall and causes a residual tongue of ice to remain during the spring and summer. High productivity is often associated with such a residual ice edge, making sea-ice fluctuations particularly relevant to biological activity (Reference Ackley, Buck and TaguchiAckley and others 1979, Reference AlexanderAlexander 1980).
In addition to these features, sea-ice drift and air-temperature data from drifting buoys(Reference AckleyAckley 1981[b]), make the Weddell Sea a particularly useful region for a first test of an Antarctic ice model. In this paper we examine some initial results from a seasonal dynamic-thermodynamic simulation of the Weddell Sea pack ice. In these simulations we make use of an extension of the dynamic-thermodynamic sea-ice model developed by Reference HiblerHibler (1979). This study differs from previous seasonal Antarctic simulations (Reference PeasePease 1975, Reference Washington, Semtner, Parkinson and MorrisonWashington and others 1976, Reference Parkinson and WashingtonParkinson and Washington 1979) in two important respects:(1) real-time varying forcing is used to drive the model instead of climatological data, and (2) a realistic rheology and thickness coupling is used to describe the ice dynamics instead of eitlier no ice motion or ad hoc ice-velocity corrections without reference to an ice rheology.
2. Description of the model
Basically the model consists of a momentum balance coupled to equations describing the evolution of ice thickness characteristics. Except for the surface heat budget (Reference HiblerHibler 1980[b]) a complete description of these equations is given in Reference HiblerHibler (1979). We summarize briefly the main characteristics here. In Cartesian coordinates the governing equations are:
where D/Dt = ∂/∂t + u·∇ is the substantial time derivative, k a unit vector normal to the surface, u the ice velocity, f the Coriolis parameter, m the Tee mass per unit area, τa and τW forces due to air and water stress, H the sea-surface dynamic height, g the acceleration due to gravity, A the fraction of area covered by ice, h the average areal ice thickness (equal to the ice-mass per unit area divided by the ice density), and F the force due to variations in internal ice stress. SA and Sh are thermodynamic growth and decay terms and are described below. For modeling the ice interaction the ice is considered to have a nonlinear viscous-plastic constitutive law σjj = σjj(⋵ij, P*) where σjj is the two-dimensional stress tensor, ⋵ijthe strain-rate tensor, and P* the ice strength. The constitutive law employs an elliptical yield curve passing through the origin with a no-stress condition applying for pure divergence (see Reference HiblerHibler (1977) for explicit equations). Rigid plastic behavior is approximated in the rheology by allowing the ice to flow plastically for normal strain-rates and to creep in a linear viscous manner for very small strain-rates. The viscosity used for small strain-rates is, however, very large. As a consequence, in the linear regime, simulated deformation rates are small enough to approximate rigid behavior well.
The coupling of the ice strength to the ice thickness characteristics is given by
where P0 and C are strength constants
This formulation makes the ice strength P* strongly dependent on the amount of thin ice while also allowing the ice to strengthen as i t becomes thicker. For the calculations here we take Pρ = 2.75x104 N m−2 and C = 20. This P0 value is about five times larger than the P0 value used by Reference HiblerHibler (1979). This increase was necessitated by the fact that daily wind fields were used here as opposed to the eight-day average wind fields used by Reference HiblerHibler (1979). On the average, these daily wind fields yielded higher wind stresses than the eight-day averaged data. As a consequence a higher ice strength was needed to obtain realistic drift rates. The value of P* used here is the same as the best fit value obtained by Hibler and Walsh (to be published) for the Arctic.
The thermodynamic terms Sh and SA are given by
with f(h) the growth rate of ice of thickness h, and h0 a fixed demarcation thickness between thin and thick ice. The Sh term gives the net growth or melt of ice while SA characterizes the way in which growth and decay change the relative areal extents of open water and ice. The h0 parameter dictates how rapidly open water is removed by growth. A more complete discussion of these equations is given by Reference HiblerHibler (1979). For these Antarctic simulations h0 is set equal to 1.0, as opposed to 0.5 used in Arctic simulations (Reference HiblerHibler 1979, Hibler and Walsh to be published). This larger value is based on field observations reported in Reference Ackley, Gow, Buck and GoldenAckley and others (1980) and Reference Gow, Ackley, Weeks and GovoniGow and others (1982) which suggest that frazil ice formation may well prolong ice production in open-water regions. Other than h0 and P0, the numerical parameters used in the simulations are identical to those used by Reference HiblerHibler (1979).
The vertical growth rates f(h) are estimated by employing Reference SemtnerSemtner's (1976) sea-ice thermodynamic model in conjunction with a surface heat-budget calculation similar to that of Reference Parkinson and WashingtonParkinson and Washington (1979) and Reference Manabe, Bryan and SpelmanManabe and others (1979). The basic surface energy-budget calculation (for details see Reference HiblerHibler 1980[b]: appendix B) includes incoming shortwave and long-wave radiation and sensible and latent heat fluxes, and provides for the determination of the surface temperature of the ice by numerical iteration.
Following Reference Manabe, Bryan and SpelmanManabe and others (1979), the effects of snow cover are approximated by allowing the icesurface albedo to be that of snow when the calculated surface temperature is below freezing, and that of snow-free ice when the surface temperature of the ice is at the melting point. The model also includes a mixed layer of fixed depth (30 m) and assumes a constant oceanic heat flux of 2 W m−2 (the same as that used in Arctic simulations (Reference HiblerHibler 1980[b]) into this mixed layer from below. Following Equation 5, when open water is absorbing heat, all of the absorbed heat is used to melt the ice further. If there is no ice the absorbed heat raises the temperature of the mixed layer. Under growth conditions, no ice is allowed to form until the mixed layer reaches the freezing temperature of sea-water.
These coupled equations are solved numerically using a time-marching procedure with one-day time steps. In this procedure the momentum equations are solved implicitly using relaxation methods, and the thickness equations are solved explicitly. More details on the numerical scheme are presented by Reference HiblerHibler (1979). In addition documentation for all portions of the code except the heat budget is available in Reference HiblerHibler (1980[a]).
3. Simulation results
To investigate the applicability of this model to the Weddell Sea pack ice, a series of two-year seasonal simulations employing different input values of the atmospheric temperature and humidity fields were carried out. Daily atmospheric data from 1 January 1979 to December 1979 were used to drive the model with the same atmospheric forcing being repeated the second year. In all cases the initial conditions were chosen to be 2 m of ice everywhere. After two years the model was found to be close enough to equilibrium to yield useful sensitivity analyses. More specifically, after this time the model results were found to be relatively independent of initial conditions and to vary little between sequential seasonal cycles. The particular time interval used in these experiments coincided with the deployment of several drifting buoys in the Weddell Sea (Reference AckleyAckley 1981[b]) which supplied valuable drift, temperature, and atmosphericpressure data. Agrid with 2° latitudinal (222 km) resolution based on an azimuthal equidistant projection (Fig.l) was used for the simulations. Zero ice velocities were assumed at all boundaries. However, at the shaded grid cells ttie ice strength was taken to be zero. In addition, ice is only allowed to be transferred into or out of these cells by advection and once there i s considered to have flowed out of the basin. In practice these cells provide a natural outflow and/or inflow boundary condition, and are discussed in more detail by Reference HiblerHibler (1979, Reference Hibler1980[a]).
Fig.1. Grid used for numerical simulations. Dashed lines enclose area where temperature and humidity fields were re-analyzed. Hatched grid cells represent open boundaries and positions of two drifting buoys on 1 January 1979 are shown (buoy 527 and buoy 1433).
(a) Input fields
Input fields to drive the model consisted of average daily surface atmospheric air temperatures, mixing ratios, and surface pressure taken from the 1979 Australian analysis. Preliminary simulations with these analyzed atmospheric temperature fields yielded unrealistically small ice concentrations in summer. Comparison of analyzed surface temperatures with buoy observations (Fig.2) indicated that the main problem was the excessively warm temperatures in the late summer and fall(February to April) in the Australian surface analysis, which did not include these buoy observations. Consequently, a more realistic surface analysis was constructed by not allowing (except for January) air temperatures and mixing ratios to rise above 271.2 K and 0.25 g kg−1in the region designated by the dashed lines in Figure 1.
Fig.2. Comparison of observed surface temperatures with those derived from Australian analysis temperatures at location of buoy 1433.
Experiments with slightly different regions had little effect on the results.
For the calculation of geostrophic ocean-current fields, mean dynamic topography values reported by Reference Gordon, Molinelli and BakerGordon and others (1978) were used. In the radiation calculations, parameterizations similar to those employed by Reference Parkinson and WashingtonParkinson and Washington (1979) were used. Specifically, daily global solar radiation under cloudless skies was obtained by integrating an empirical equation by Reference Zillman and HayesZillman (1972) over solar zenith angles for any particular day. Zenith angles at half-hour intervals for this purpose were obtained from a numerical solution of Kepler's equation (Holloway personal communication). Incoming long-wave radiation was obtained using Reference Idso and JacksonIdso and Jackson's (1969) formula for clear skies. For cloud cover, a constant value of 0.85 was assumed at all points. This value is commensurate with maximum estimates by Reference Van Loon and NewtonVan Loon (1972). Experiments with other values were found to modify the simulation results minimally.
(b) Sensitivity of ice-edge fluctuations to ice dynamics
Using the input fields discussed above, two-year simulations were carried out using either the full coupled model or, for comparison, only the thermodynamic portion of the model. Comparison of these results gives some insight into the role of ice dynamics on the ice-edge fluctuations.
Some ice-extent results from these simulations are illustrated in Figure 3, which shows plots on two specific days of the compactness for the full model and thin-ice contours for the thermodynamics-only case. The observations are based on charts of northern ice limits (prepared by the US Navy-NOAA Joint Ice Center, Suitland, Maryland, USA) which report ice having a concentration of at least 50%. In most cases this 50% limit is effectively at the ice edge estimated in the charts. In this figure two different days were chosen to represent salient features of both the simulated and observed ice extents. Day 18 (January 18) represents a period close to the end of the rapid decay and shows a characteristic geometric tongue of ice often seen at this time of year. Day 298 (October 25), on the other hand, represents a typical maximum extent at a time just prior to the beginning of the rapid decay phase.
Fig.3. Simulated and observed ice concentration and extent for full coupled model ((a),(b)) and for thermodynamics-only case ((c),(d)). Dashed line denotes the observed limit of ice having a concentration of at least 50%. Dotted line shows interior polynyas having concentrations less than 87.5%.
As can be seen from Figure 3 on day 18 the full model reproduces the tongue effect and ice extent much better than the thermodynamics-only case. Analysis of the ice-velocity field (see section (d)) shows that this tongue is largely due to eastward and northward advectlon of ice by the ice-velocity field. (The tongue is, however, slightly far south, partly due to insufficient transport of ice northward as discussed in the next section.) In the thermodynamics-only case the structure of the ice edge is totally incorrect and there is too much ice.
On day 298 the results are, however, quite similar, with both the full model and the thermodynamic case showing extents slightly less than observed. This in turn suggests that the expansion to maximum is largely thermodynamic in nature and controlled (in both simulations) by the temperature field.
The effect of dynamics on the seasonal variation is illustrated more graphically in Figure 4, which shows the observed time series of the total grid area covered by ice with concentrations greater than 15%. Also shown in this figure are the time series of the thermodynamics-only case where the area given is for ice of non-zero thickness. As can be seen, the coupled model with the dynamics included gives a bigger seasonal swing of extent and a more rapid decay of the ice cover. Both these effects are primarily due to the large amount of deformation and advection caused by the dynamics, which in turn creates open water. As the periods of short-wave radiation increase during spring and summer this open water absorbs large amounts of heat and causes rapid decay of the ice cover. However, one difficulty with the full model (and the thermodynamics-only case) is a lack of persistence of the ice cover in the spring (day 250 on). Specifically, the simulated ice cover starts decaying sooner and in a smoother manner than observed. A feedback parameterization is proposed in the next section which modifies this feature considerably.
Fig.4. Time series of total area covered by ice within simulation region. Data were plotted every two days. In both simulated and observed cases, the resolution element was one 220 × 220 km grid cell (see Fig.l).
An additional feature of some interest is the effective area covered by ice (area x compactness) versus total area. Time series of these two quantities are compared in Figure 5 for the full model.
Fig.5. Total area and effective area covered by ice for standard simulation. Effective area is total area x average compaction
The difference emphasizes the large amounts of open water in the Antarctic pack ice in all seasons. This feature is in contrast to the ice charts of northern ice limits prepared by the US Navy-NOAA Joint Ice Center, but is in agreement with concentration estimates made using electrically scanning microwave radiometer Imagery (Reference Zwally, Parkinson, Carsey, Gloersen, Campbell and RamseierZwally and others (1979) as referenced in Reference AckleyAckley (1981[a]).
(C) An ice-atmosphere feedback mechanism
To improve the persistence of the simulated ice cover in spring, a simple feedback mechanism between ice extent and atmospheric temperature was constructed. The dynamic-thermodynamlc coupled simulations including this mechanism yielded quite realistic results, and suggest that such a feedback may play much more of a role in sea-Ice fluctuations than had previously been considered. However, the feedback with the thermodynamics only model produced highly unrealistic results, a fact which further indicates that realistic dynamics must be included 1n sea-ice models in order to obtain useful results.
The feedback mechanism consisted of specifying that the air temperatures and humidity are not allowed to rise above 271.2 K and 0.25 g kg−1respectively if the ice-cover compactness is greater than 0.50. This correction was applied at all grid cells. The input fields used in this case were then taken to be the standard Australian analysis data without the simple re-analysis described In the previous section. It is notable that such a feedback mechanism is consistent with field observation of an abrupt change of air temperature as the ice edge is crossed. Also, measurements of atmospheric boundary layer modification (Reference Andreas, Paulson, Williams, Lindsay and BusingerAndreas and others 1979) support the assertion of a feedback between atmospheric temperature and ice extent.
Two-year simulations employing this feedback with both the full dynamlc-thermodynamic model and a thermodynamics-only case were carried out. In the thermodynamics-only case the feedback was turned off when the 1ce thickness was less than 10 mm. The seasonal variation of ice extent from these simulations is shown in Figure 6. As can be seen, the full model reproduces a realistic seasonal variation and geographical extent. In addition it produces the persistence in the ice extent in the spring (day 250 onward) that is in better agreement with observation than the standard case discussed above without significantly modifying the rapid decay seen in the summer period (December to January). Also, although not shown, the feedback produces a realistic ice tongue similar to the standard case discussed previously. The thermodynamics-only case, on the other hand, effectively produces a constant ice edge with very little seasonal cycle.
Fig.6. Total area covered by ice for simulation including feedback between atmospheric temperatures and presence of ice. The dashed “feedback” curve here refers to full dynamic-thermodynamlc simulation.
(d) Simulated ice drift
Some of the ice-velocity characteristics are illustrated in Figure 7.This figure compares the 30-day simulated and observed ice velocities of two drifting buoys initially deployed in the western Weddell Sea for the standard case. Similar results were obtained for the feedback cases. Also shown in this figure is the final location of a simulated buoy trajectory. For buoy 527 the trajectory and velocity vectors are similar, indicating a relatively smoothlyvarying velocity field. For buoy 1433 slower velocities near the coast substantially reduce the trajectory simulation to about day 120; these velocities are commensurate with other data from drifting buoys obtained during this time (Reference AckleyAckley 1981[b]).
Fig.7. Simulated and observed 30 d average velocity vectors for two drifting buoys (a),(b). Simulated velocity vectors were obtained by finding, through interpolation, simulated velocity vectors at actual buoy locations. In addition to velocity vectors, final position of simulated trajectory is shown as open circle. (Trajectory uses simulated positions of buoys rather than actual positions as used in velocity vector plot.) Velocity vectors have been converted to displacements by a multiplication by time.
While there are some differences, Figure 7 shows the simulated and observed drift to have similar magnitude and direction. There is, however, a larger observed northward displacement near the coast than simulated. The most likely reason for this discrepancy is a strong thermal wind due to damming of cold air against the Antarctic Peninsula as suggested by Reference SchwerdtfegerSchwerdtfeger (1979). Such thermal wind effects are not considered in the simple geostrophic wind analysis performed here.
4. Discussion
The initial simulations reported here have yielded several insights Into modeling Antarctic seaice fluctuations. The main effect is the role of dynamics in creating a much larger seasonal swing of sea-ice extent in agreement with observation. A second effect of the dynamics is to decrease the sensitivity of the model to feedback between ice extent and the atmosphere. This result has important implications for climate models.
For the forcing employed here it is not necessary to increase the oceanic heat flux above typical Arctic values to obtain a realistic seasonal cycle when realistic dynamics are included. This is in contrast to calculations based primarily on thermodynamic considerations by Reference Parkinson and WashingtonParkinson and Washington (1979) and Reference GordonGordon (1981), which have suggested the need for a larger oceanic heat flux to explain the ice decay. However, further sensitivity studies and analyses are needed to examine the behavior of this dynamicthermodynamic model of the Weddell Sea pack ice in more detail. Such work is in progress.
Acknowledgements
The authors thank Betsy Holt for her assistance in reducing and analyzing the buoy data used in the paper. Comments on the manuscript of this paper by I Allison were most helpful. The work was supported by US National Aeronautics and Space Administration grant S-79806-B and US National Science Foundation grants DPP-77-24528 A02 and DPP-8006922, and by the US Army Cold Regions Research and Engineering Laboratory. This support is gratefully acknowledged.
