Introduction
A typical sea-ice cover in the polar regions contains a variety of ice thicknesses which evolve in response to both dynamic and thermodynamic forcing. The variable thickness of the ice cover is created by deformation which simultaneously causes formation of thick ice through ridge building and thin ice through lead creation. Since the energy expended in deformation is largely determined by the ridging process, investigating the energetics of pressure ridging is critical for the determination of ice strength on a geophysical scale.
A framework for describing the variable thickness character of sea ice was developed by Thorndike and others (Reference Thorndike, Rothrock, Maykut and Colony1975). Within this framework, pressure ridging is treated statistically by a redistribution process whereby thin ice is transferred to thick ice categories. The amount of energy used for this process can be related to the large-scale strength of the pack ice (Rothrock, Reference Rothrock1975). A key issue here is the ratio of the total energy consumed to the potential energy created in the ridging process. Most estimates have typically placed the total energy losses at about twice the potential energy (Rothrock, Reference Rothrock1975), a factor that was based largely on a kinematic ridge model developed by Parmerter and Coon (Reference Parmerter and Coon1972). However, while modeling the kinematics of the ridging process, this ridge model did not explicitly calculate frictional and other non-potential energy losses. The factor of two was also inconsistent with large-scale simulation results. In a seasonal simulation of the Arctic basin using a variable thickness sea-ice model, Hibler (Reference Hibler, III.1980) estimated the ratio of total work to the change in potential energy to be between two and ten.
In their classical paper, Parmerter and Coon (Reference Parmerter and Coon1972) constructed a two-dimensional kinematic model of the ridging process. They began with two ice floes of equal thickness separated by a lead uniformly filled with rubble.
As the floes converged, an arbitrary fraction of the displaced rubble was thrust above the sheet and the remainder below. The rubble thrust above and below the floe was shaped by two mechanisms. A lateral transfer coefficient was used to define the fraction of rubble carried with the floe toward the center of the lead and the fraction which remained in place. An assumed angle of repose shaped the rubble piles forming above and below the water level. The floes themselves were treated structurally as plates on an elastic foundation.
Hopkins and others (Reference Hopkins, Hibler, III. and Flato.1991) adopted a similar picture of ridge growth. However, in place of the arbitrary ridge building mechanisms used by Parmerter and Coon, a particle simulation model was used in which the dynamics of the blocks of rubble were calculated explicitly. Since frictional and inelastic contacts were included, this model furnished a more complete picture of the energetics of the ridging process than the lower bound provided by Parmerter and Coon, which was based on potential energy alone. This type of ridge growth characteristically shows a small sail to keel ratio, because the layer of rubble remains in hydrostatic equilibrium until the multi-year floes approach each other.
Whereas instances of ridges grown from rubble-filled leads probably do occur in the central Arctic pack, this model seems better suited to areas of the pack in the vicinity of land. A more likely model for ridge growth in the central pack is for a plate of thin ice covering a newly frozen lead to be driven against a multi-year floe by the converging pack. With enough lead ice and a sustained period of convergence, blocks broken from the leading edge of the lead ice will collect above and beneath the multi-year floe to form the characteristic ridge structure seen in the central Arctic.
In the present work, the latter type of ridge growth is modeled. The simulation begins with a thin plate of ice moving at a constant speed toward a thick multi-year floe (see the first snapshot in Fig. 2). The lead ice is modeled as a plate on an elastic foundation using standard solutions (Hentenyi, Reference Hetenyi1946) for vertical and angular deflection and moment at regular points in the plate. As the lead ice collides with the floe, it bends. When the moment at a point in the lead ice exceeds its flexural strength, it breaks at that point. As the lead ice continues to move toward the floe, additional blocks break off and contribute to the ridge structure. At the same time, the motions of the blocks which make up the growing ridge are calculated based on contact forces between blocks and body forces. The individual blocks may also fail in flexure. The energy dissipated at every contact is calculated explicitly, as well as the changes in potential energy of the rubble blocks and the thin ice plate. The overall energy consumed is found by calculating the work performed by the horizontal force pushing the plate of lead ice into the ridge. The multi-year floe is treated as a rigid plate.
The Particle Simulation
In general terms, a particle simulation is a computer program which models the Newtonian dynamics of a large system of discrete particles. The simulation stores the shape and the instantaneous position and velocity of each particle. The contact and body forces on each particle are calculated at each time step and the particles are moved to new locations with new velocities which depend on the resultant of the forces. In the following paragraphs, the general details of the particle simulation used in this paper are described. A detailed description of simulation methods using polygonal block-shaped particles may be found in Hopkins and Hibler (in press) and Walton (Reference Walton1980).
The domain of the simulation is defined sufficiently far above, below and to each side of the lead that the domain will contain the entire ridge at all times. The rubble is composed of block-shaped polygons of uniform thickness which were broken from the lead ice. In order to model the contact forces between blocks, a list of neighboring pairs of blocks is compiled periodically. The method used to compile the list begins by overlaying the domain with a square grid. Each cell of the grid is small enough that at most one block center may lie in the cell at any time. A two-dimensional array is dimensioned to correspond with this grid and the index of each block is assigned to the array element corresponding to the cell in which its center lies. In addition, two inverse address arrays are defined. One holds the i location of a block and the other the j location, where i denotes a column and j a row in the grid. These arrays enable the program to find the block located at (i,j) in the domain and the (i,j) address of a given block. With this information, it is a straightforward matter to find the close neighbors of each block. The pairs of blocks which are discovered less than some arbitrarily small distance apart are placed in an array of near neighbors, which is refreshed periodically.
At each time step, Δt, each pair of blocks in the array of near neighbors is checked for contact. Contact exists when the blocks overlap. If a pair of blocks is found to overlap (Fig. 1a), a local orthogonal coordinate frame (n,t) is established with its origin at the centroid of the area of
(a) A pair of overlapping ice blocks, (b) The contact force model.
overlap. The normal axis n is perpendicular and the tangential axis t is tangent to the surface of (at least one of) the blocks.
A representation of the interparticle force model is shown in Figure 1b. The normal contact force F n has an elastic component proportional to the area of contact A c and a viscous or dissipative component proportional to the normal component of the relative velocity at the point of contact V rel
where k ns is an elastic spring constant, k nd is a viscous damping constant, and the superscript n denotes the current time step. Tensile forces are not allowed.
In the tangential direction, the blocks are coupled by a sliding surface to a linear spring. The sliding surfaces are characterized by a constant coefficient of friction μ. The initial relative tangential motion at the point of contact compresses the spring creating a force. The magnitude of the force increases incrementally from time step to time step until it reaches the limiting value of μ∣F n∣. The magnitude of the frictional force may not exceed that value; instead, sliding occurs at the contact surface. The equation defining the tangential force F t (Walton, Reference Walton1980) is
with ∣F t∣ ≤ μ∣F n∣ In Equation (2) k ts is an elastic spring constant. The moments of F n and F t, about the center of each block, follow from the geometry of the contact.
After the body forces and the forces and moments exerted on each block by surrounding blocks have been found, the equations of motion are derived from a Taylor series expansion about the current time. The equations are second order accurate with velocity dependent forces. The velocity of a block at time n + 1/2 is
In Equation (3) the subscript l denotes direction in Cartesian space, F is the resultant force on the block, and M is the mass of the block. The velocity of a block at time n + 1 is
The position of a block at time n + 1 is
Similar equations determine the angular orientation θ and velocity ω of the block.
Because of the tendency for polygonal blocks to form stable load-bearing structures by interlocking or stacking, some form of breakage mechanism is necessary. A simple breakage model suitable for slender blocks, in which the blocks fail in flexure, was developed by Hopkins and Hibler (in press). In this model, breakage is caused by flexure in the x, z plane. The d'Alembert force and torque due to the contact forces and unbalanced body forces on the block are calculated. Given the d'Alembert force and torque the calculation of the moment at a point in the block is reduced to a statics problem. If the moment at a point exceeds the flexural strength of the block, then the block is broken at that point. The flexural strength of the block is (h 2/6) σ cr (Timoshenko, Reference Timoshenko1956), where h is the block thickness and σ cr is the tensile strength of sea ice.
The vertical and angular deflections of the lead ice are calculated using standard solutions (Hetenyi, Reference Hetenyi1946) in which plate rigidity replaces beam stiffness. These are equilibrium solutions which are assumed to apply to the slowly moving plate. At regular time intervals during the simulation, the deflections are calculated at regular points in the lead ice. The calculations use forces which are the averages of the forces exerted on the lead ice by the rubble during the preceding interval. Vertical and angular velocities are calculated to bring each point in the plate to its new position over a second time interval. In the ridge simulation discussed in this paper, deflections were calculated at 0.1 s intervals and velocities were calculated to bring each point in the plate to the new position in 2 s. This lag was necessary in order to couple the lead ice to the rubble, which has a slower response time. At the same time, moments were calculated at each point. If the moment at a point exceeded the flexural strength of the lead ice, it was broken at that point and the broken piece was added to the rubble.
The Energetics
The work performed in building the ridge is transformed into the potential energy of the ridge structure and lost through dissipative mechanisms. The dissipative mechanisms consist of frictional contacts between sliding blocks and inelasticity, which implies crushing at the point of contact, implemented by the dashpot in the normal force (Equation 1). The stored elastic energy which is lost when the lead ice breaks is implicitly accounted for as unrestored potential energy. The overall ridge building work is calculated directly from the horizontal forces Fx exerted by the rubble blocks in contact with the lead ice moving with speed u lead
The summations in Equation (6) are over the forces exerted by the rubble on the lead ice (Σi) at each time step (Σn) in the experiment. The change in the potential energy of the ridge structure is
The summations in Equation (7) are over each block in the rubble at each time step. A is the area and w the vertical component of the velocity of block i; As is the submerged area and w s the vertical component of the velocity of the center of mass of the submerged area of block i; ρ i is the sea-ice density and ρ w is the sea-water density. The frictional dissipation, Φ f, is the work performed by the tangential component of the contact forces
Similarly, the dissipation Φ d by the dashpot in the normal direction, which is equivalent to particle inelasticity, is the net work performed by the normal component of the contact forces
In both Equations (8) and (9), the summations are over the pairs of blocks in contact at each time step. The changes in kinetic energy of the rubble are also calculated and in general remain negligible compared to the above processes.
Discussion of the Results
The results given here are from a single preliminary experiment with the ridge simulation described above.
Fifteen snapshots of simulated ridge growth at 50 s intervals.
Snapshots of the simulated ridge growth are shown in Figure 2. The following parameters were used in the simulation:
Separate values of the friction coefficient were used for dry contacts (above the water) and wet contacts (underwater). The value of Young's modulus used in the simulation is much less than values typically cited. Young's modulus occurs in the elastic beam equations where it controls the length of the blocks broken from the lead ice plate. The value used results in blocks with a typical length of about 1.5 m. Young's modulus also affects the rigidity of the plate, which may have a secondary effect on the energetics.
The rates of work, dissipation, and increase of potential energy are shown in Figure 3 for the simulated ridge growth in Figure 2. The overall rate of work increases in an approximately linear fashion throughout the experiment, whereas the rate of increase of potential energy remains relatively constant. This behavior may be explained in the following way. Consider the section of rubble directly beneath the plate of lead ice. The overall work is largely a function of the frictional force exerted on the underside of the lead ice. The frictional force, which depends on the buoyancy of the rubble, is proportional to the area of the rubble. Since the lead ice is fed into the ridge at a constant rate, the area of rubble beneath the
The rates of work, dissipation and potential energy change for the ridge shown in Figure 2.
plate, the frictional force, and the overall work, will increase linearly with time. The potential energy of the rubble beneath the plate is also proportional to the area of the rubble. Therefore, the increase of potential energy will be linear, and its derivative (plotted in Fig. 3) will be constant. The ratio of the overall rate of work to the rate of increase of potential energy in Figure 3 varied between 15:1 and 20:1. In the ridges grown from a rubble lead (Hopkins and others, Reference Hopkins, Hibler, III. and Flato.1991) the ratio was between 3:1 and 5:1 depending on the coefficient of friction.
Conclusions
The intent of this work has been to give a brief description of a model for the simulation of the formation of pressure ridges in situations where relatively thin, intact lead ice is driven against a thick floe. The ridging model uses a dynamic particle simulation in which the motions of ice blocks are a consequence of explicitly modeled contact and body forces on each block. A preliminary numerical experiment with the ridging model furnished pictures of simulated ridge growth (Fig. 2) and an estimate of the energetics of the ridging process (Fig. 3). In this experiment, the total energy consumed in ridging was much higher than the value of two times the increase of potential energy estimated by Rothrock (Reference Rothrock1975). It was also significantly higher than the estimates of Hopkins and others (Reference Hopkins, Hibler, III. and Flato.1991) of three to five times the increase of potential energy (depending on friction), for ridges grown from a rubble-filled lead. The result implies that ridges grown from an intact frozen lead require much more energy per unit increase of potential energy in the ridge structure than ridges grown from a rubble-filled lead.
Since these are preliminary results, they are offered with some trepidation and should be interpreted with caution. The results may depend strongly on the speed of the lead ice, which was very high in the simulated ridge described above or on other factors yet to be considered. Future experiments will examine the sensitivity of the ridging energy budget to variations in speed, friction, flexural strength, Young's modulus, and plate thickness.
Acknowledgements
The authors would like to acknowledge Greg Flato's help in clarifying questions about the behavior of the energetics of this ridge model. This work was supported by the Office of Naval Research through the University Research Initiative in Arctic Ice-Ocean Dynamics and Mechanics, contract number N00014-86-K-069.
