Numerical Modelling of Realistic Ice Floes in Ocean Waves

Abstract

Results from a numerical model for the computation of ice-floe motions in ocean waves are presented and discussed for floes of various sectional shapes. It is shown that the beam-to-draft ratio is a crucial factor in determining the behaviour of each floe, and that ridges and keels can substantially affect the roll and sway characteristics, particularly for thick floes. Undercuts beneath floes appear to have little effect, but a protruding sill can lead to decreased motion at certain frequencies and a reduced capability for making waves. As the underwater character of an ice floe cannot easily be measured, it is necessary to infer draft from freeboard estimates using the method proposed by Ackley and others (1976). Although most of this paper is concerned with rigid body motions, a brief account is included of a modification to the model which allows each floe to bend on the passing wave. By this means, wave-induced flexure is discussed for floes of non-simple underwater shape.

1. Introduction

The Scott Polar Research Institute has studied the interaction of ocean waves with pack ice since 1959 when Robin (1963) first used a shipborne wave recorder to measure waves in Antarctica. Robin’s approach and that of subsequent early studies was to consider the marginal ice zone (MIZ) as a whole and to evaluate its effect on the incident ocean-wave spectrum. The presence of the ice cover was found to affect the passage of ocean waves substantially. Waves penetrating from the open ocean were preferentially attenuated by the sea ice, and their observed spectrum was seen to change in shape and content with distance from the ice edge. No detailed measurements of the motion of the individual ice floes making up the MIZ were collected, however. More recently, the importance of determining these motions has been recognized and accelerometers and strainmeters have been used on single ice floes to collect data on heave, sway, surge, roll and bending (Wadhams 1979, wadhams and Squire 1980). The results from these experiments will be presented soon. The present work introduces a theoretical model which has been developed to interpret and to simulate the results of these studies. More details may be found in Squire (1981).

2. The Numerical Model

We expect two results from an ice floe in waves. Firstly, the waves will diffract upon the floe and will scatter in all directions and, secondly, they will set the floe into motion and this motion will itself generate waves moving out in all directions. Although the first breakthrough in this work came in the last century (Froude 1861), it is only recently that a full solution including both contributions has been possible. Even the formulation of the problem is difficult, since the coordinate frame which is suited best to describe the motion of a floe moves with the floe, while the coordinate frame which describes the fluid most easily is an inertial system which moves with the fluid. Another difficulty is associated with added mass. Added mass is a contribution to the mass or inertia of the body due to the fluid entrained by the motion. In waves it is a frequency-dependent quantity which is not easy to calculate but cannot be neglected. In a linear theory the body’s ability to create and sustain outgoing surface waves is equivalent to a frequency-dependent damping coefficient in the equations of motion. Damping is therefore also crucial to the full solution. When added mass and damping have been determined for all required frequencies. the full solution may proceed.

There are two different approaches to solving the hydrodynamic problem of a body {floating or otherwise) interacting with ocean waves: Morison’s equation and potential theory. Each method has been discussed in detail in the literature (see, for example, Sarpkaya and Isaacson 1981) and has been shown to be valid over a different range of applicability, Morison’s equation is used consistently by offshore engineers whenever inertia and viscous forces dominate, whereas potential theory becomes important when diffraction influences the fluid motions (in the small body limit, the potential theory tends rigorously to Morison’s equation when viscous effects are neglected). Since potential theory is used throughout the present study, there is some question regarding its validity when viscous drag becomes significant. This occurs either when the incoming wave height 2a is larger than the body’s diameter 2b or when b is much smaller than the wavelength λ. This is shown diagrammatically in Figure 1 reproduced with modification from Garrison (1979), where only those values beneath the dashed curve are realistic due to restrictions of wave steepness. In the overlap region both theories are valid but become less exact in some undefined way which depends on the shape of the body, so that potential theory appears at least to be reasonable over a wide range of floe diameters and incident wavelengths. If the floe to be modelled has sharp corners, viscous drag effects will no doubt be more serious since the flow will separate and development of a wake will occur. This will be particularly significant in roll. Unfortunately, no reliable data are available to establish the magnitude of the wake-induced drag and hence we are unable to estimate the necessary non-linear drag coefficient.

Fig. 1. Regions of validity of Morison’s equation and potential theory. Wave amplitude: a, floe beam (diameter): 2b, wavelength: λ (after Garrison 1979).

There are two principal methods by which added mass, damping and the wave-induced motions may be evaluated using potential theory. The first, developed by Ursell (1949) for semi-immersed circular cylinders and later extended to more general sectional shapes by conformai mapping, uses a multipole series expansion. The second is a numerical technique which models the cross-sectional contour of the body beneath the waterline by a finite number of nodal points (Fig.2.). The unknown velocity potential is then represented by a distribution of wave sources whose varying strengths are to be determined fron an integral equation found by applying the usual kinematic boundary condition. This equation is solved algebraically by assuming the source strengths to be constant along each nodal segment but to vary from segment to segment. Frank (1967[b]) and Lee (1976) developed this technique to interpret the motions of a ship by replacing the ship with a series of two-dimensional sections across the beam, and then integrating lengthwise to gain a pseudo three-dimensional solution. They, and other authors, have found good agreement between theory and model experiments (Vugts 1968, personal communication from C M Lee 1980). In contrast, Faltinsen and Michel sen (1974) report disagreements between the two procedures for bulbous-shaped bodies where the length, beam and draft are comparable in size. There is no doubt that a three-dimensional solution is desirable. However, the simpler approach discussed in this paper can be defended on three counts: first, the author knows of no similar comparison for floating bodies of small draft where perhaps a more serious criticism is the tacit assumption of rigidity; second, the full solution is impractical for the large number of calculations required for MIZ modelling due to prohibitive computing times; and, third, typical ocean-wave records are usually smoothed by between 5 and 20 contiguous spectral values implying errors between 45 and 22% on the smoothed estimate, so there is therefore little point in using overly complex solutions when the two-dimensional estimate is usually much better than this.

Fig. 2. Schematic diagram of floe showing coordinate axes and discretization scheme used in model, and defining the in-plane bodily motions.

3. Geometrically Simple Floes

In this paper we assume that long-crested waves of typical ocean-wave periods are available. We shall adopt the nomenclature of naval engineering throughout in order to emphasize the similarity between this work and studies to evaluate motions of ship in a seaway. Only beam-on waves are considered so that pitch, surge and yaw will not be induced, and the floe will respond in heave, roll and sway alone. We begin by discussing these motions for an ice floe whose sectional shape is a simple rectangle of beam 2b and draft d. Since the response of geometrically similar floating bodies may be scaled while the Froude number remains constant, it is natural to present the motions in the form of a response amplitude operator (RAO), [n such a plot, the gain factor, defined as motion amplitude/wave amplitude, is plotted against the wavelength λ of the incoming ocean wave normalized by the half-beam b of the body. Hence all ice floes of geometrically similar shape will appear as a single curve, and deviations from that curve will only be Introduced by changes in geometric similitude. For the simple rectangular section, we therefore consider the effect of various beam-to-draft ratios it on the induced motions. Extreme values of П=1 and П=100 with four intermediate values of П=2, 3, 10, and 30 are used in the simulations. The case of П=l is presented separately, as in this case the stability of the floe is marginal. Also, calculations of response for long, thin floes at small values of λ/b are likely to be approximate as the varying wave-pressure field beneath the floe will induce a significant flexural response. This bending may be calculated and will be discussed briefly later in the paper.

The heave RAO for various values of n is shown in Figure 3. For very short waves, most of the energy is reflected and the floe experiences little movement. At very long periods, on the other hand, where the length of the floe is small compared to the wavelength, the floe response is perfect (1:1). Between these two extremes lies a region where n is critical in determining the induced motion. For long thin floes, heave increases monotonically with period, but as the thickness becomes significant so a mid-range resonant peak appears and begins to dominate the motion. The peak becomes more severe as Π decreases. A spar-loaded wave buoy, for example, would possess a very distorted RAO and would tend to give a biased estimate of wave energy at mid-range periods.

Fig. 3. RAO in heave at beam-to-draft ratios of П=2, 3, 10, 30, 100. Non-dimensional wavelength: λ/b, gain factor is the amplitude of heaving motion of the floe/amplitude of forcing wave.

The equivalent roll RAO is shown in Figure 4. Once again short waves have little effect on the floe, but in this case there is also little motion when the wavelength is much greater than the beam of the floe. This is to be expected since the relative slope of the sea surface due to ocean waves becomes less as the non-dimensionalized wavelength increases. Within the short and long wave limits, the rolling of the floe increases rapidly to a resonant maximum and then decreases. The point at which roll resonance occurs is determined by the beam-to-draft ratio Π long, thin floes resonate when the incoming wave is a few times longer than the length of the half-beam; short, thick floes require longer waves to resonate. Furthermore as the value of II becomes larger so the damping becomes less and the peak becomes narrower, or alternatively, in the jargon of electrical engineers, the “Q” becomes larger. Unlike heave, the resonant character of the rolling of the floe is always present and becomes more dominant as the beam-to-draft ratio is increased.

Fig. 4. RAO in roll at beam-to-draft ratios of Π=2, 3, 10, 30, 100. Non-dimensional wavelength: λ/b, gain factor is the amplitude of rolling motion of the floe/amplitude of forcing wave.

In sway, the predicted motion is extremely complex (Fig.5.). Once again the motion is small for short waves and perfect for waves which are long compared to the beam of the floe, but between these limits we encounter resonance and a sway minimum. The initial peak is critically determined by n and is nearly absent when H assumes low values. The sway minimum also depends on n, but to a lesser extent, and it shifts little for the two orders of magnitude change in II. Complicated as these curves are, experimental evidence (Vugts 1968) suggests that they are correct, and are due in part to coupling between the roll and sway motions. It is significant to note that the predicted minimum corresponds to a 180° phase shift between roll and sway which has also been verified experimentally.

Fig. 5. RAO in sway at beam-to-draft ratios of Π=2, 3, 10, 30, 100. Non-dimensional wavelength: λ/b, gain factor is the amplitude of swaying motion of the floe/amplitude of forcing wave.

The final figure of this section (Fig.6.) shows the three motions when the beam-to-draft ratio It is unity. For an approximately homogeneous body with a density roughly that of ice, a value of Π=l is a lower limit since smaller values will render the floe unstable. Moreover, ice floes with this beam-to-draft ratio will tend to be small in breadth since sea ice will not grow naturally more than a few metres in thickness. We therefore suppose Figure 6 to represent the motion response of ice cakes and rubble to incident waves, and include it for completeness. All three motions exhibit severe resonance at a period which is very long compared to their half-beam. It is significant that the peak of the heave occurs before those of roll and sway, and though no data are available to substantiate this phenomenon on ice floes, it has been seen in the motion of icebergs (Kristensen and others 1982). It is an indication that a floating body is close to being unstable, an observation that is borne out by the size of the gain factor in the figure. No attempt has been made to include viscous roll damping in the simulation. Although this can be done, it is not sensible to do so until reliable estimates of the drag factor for roll damping are measured experimentally. Its effect would be to slightly decrease and broaden the roll and sway peaks.

Fig. 6. RAO for heave, roll and sway motions of an ice floe with Π=1. Non-dimensional wavelength: λ/b, gain factor is the amplitude of motion of the floe/amplitude of forcing wave.

4. Floes with Ridges, Keels and Sills

A major problem in any simulation study involving real ice floes is that each floe is unlikely to be simple in shape. This is particularly true of floes which have been subjected to intense wave activity (Squire and Moore 1980), or floes which have spent some of their life in a region where compressional forces can buckle and deform the ice cover heavily. Most of the floes in east Greenland waters are known to originate within the Arctic Ocean, and a significant proportion possess pressure ridges and corresponding keels (Wadhams 1981), Undercut floes, and floes with sills beneath the waterline, are also commonplace.

The effect of various ridge-keel configurations on the wave-induced behaviour of ice floes has been studied in detail. In all simulations a mean ridge/ keel slope angle of 32° (Wadhams 1978) is assumed, and the ratio of ridge height to keel depth is kept constant. The presence of these features can substantially affect the rotatory motions of the floe and it is essential to include an appropriately modified moment of inertia in the calculation. This is achieved using the numerical method proposed by Hall (1979) to compute the moment of inertia of the new sectional shape of the floe about its centre of flotation. Strictly, the floe rolls about an axis which moves during motion. However, the locus traced out by this axis is only known for simple shapes, and it is usual to approximate by forcing the roll axis to coincide with that for rotation under a static moment. The results of the simulations indicate two expected results: firstly, the motions of long thin floes are not affected substantially by the presence of a keel, and, secondly, a keel can severely alter the wave-induced motions when the beam-to-draft ratio Π is small. In heave, the keel is equivalent to a slight increase in thickness or, alternatively, a decrease in the value of Π. In roll, the keel modifies the locus of the rolling axis and causes a significant change in n. For short floes of large draft the keel may well cause the floe to become unstable. The roll motions of a floe with a keel are highly damped, indicating that it is good at generating outgoing surface waves. The behaviour of a floe with a keel in sway is coupled to the corresponding roll motion. The presence of a keel creates the effect of increased thickness, so the first resonance becomes less although the sway minimum is unchanged. Undercuts appear to have very little effect on floe motion until they become so unrealistically large that the ice at the edges of the floe is liable to fracture. We may therefore assume that the motion of an undercut floe is very nearly the same as the equivalent floe of rectangular section.

Ice floes with sills are particularly difficult to model. Indeed, the existence proof of John (1950) for the solution to the integral equation is not applicable to any floating body whose underwater contour has points outside the point at which it cuts the water surface. This is not to say that solutions cannot be found, and Frank (1967[a]) has evaluated motion characteristics for several partially immersed bulbous cylinders of various forms. His results are particularly interesting and indicate that the damping coefficients of such a body vanish at certain periods, i.e. that it ceases to generate outgoing surface waves as it moves in the water. This effect had been observed experimentally before Frank’s theoretical study. Figure 7 shows the results of a simulation to compute heave, roll and sway for an ice floe with a significant sill projecting for 10% of the beam beneath the waterline of the floe. It appears that the sill can modify the motion considerably. In heave, the sill induces a resonance which makes the floe behave as though it had decreased beam-to-thickness ratio. A secondary peak also develops at a wavelength comparable to the beam of the floe, close to a point where damping due to heave is particularly small. Damping in roll is also very small for values of X>4b and this has important consequences for the roll RAO. The roll maximum is split into two distinct resonances with a minimum between them, where the body ceases totally to generate outgoing waves. Sway is equally complicated, and has a pair of resonant peaks also. The usual minimum in sway is displaced slightly towards longer waves. In all cases, the asymptotic limits of the motions remain unaffected by the sill. These results, although complicated, are in agreement with those of Frank (1967[a]) and show the significant differences which can occur when any form of underwater projection is introduced.

Fig. 7. Heave, roll and sway RAOs for an ice floe with sill and with beam-to-draft ratio of 11=10. Non-dimensional wavelength: λ/b, gain factor is the amplitude of motion of the floe/amplitude of forcing wave.

5. Bending Response

A facet of the behaviour of ice floes which has not yet been addressed is their ability to bend slightly as each ocean wave passes beneath. Bending is important as, if the strains are large enough, the floe will fracture and its character will change. As each floe contributes to the response of the MIZ as a whole, wave-induced fracture will be a principal factor in determining the distribution of floe size and other significant MIZ features. Although the numerical model described above was originally devised for the purpose of studying rigid bodies, the author has extended the model to allow some flexural response. This has been achieved using a finite-element scheme in order to avoid placing unnecessary restrictions on geometry. The method is reported in detail elsewhere (Squire 1981) and we shall only summarize the technique in this paper. We regard the solution as having two independent phases: first the calculation of the motions of a rigid body and the associated pressure beneath the floe, and then the modelling of the bending after the assumption of rigidity has been relaxed and the floe is allowed to bend on the pressure field. The material properties of sea ice may be specified though we have initially assumed elastic behaviour. Implicit in this method is the assumption that the bending of the floe does not interfere with the hydrodynamics. We believe this to be a good approximation since the natural modes of flexure lie well beyond the frequencies of typi cal ocean waves.

For simple rectangular floes, the flexural solution agrees well with that for a plate on an elastic foundation, although thick-plate theory or a full three-dimensional solution must be included when n is small. Simulations for more realistically shaped floes are still in their early stages but already it appears that modification to the basic shape can lead to zones of high stress intensity and hence to preferred regions for fracture.

6. “Real” Ice Floes

The precise nature of the underwater profile of each floe cannot be determined routinely for all types of sea ice unless an unrealistically thorough coring programme or under-ice dives are conducted. Impulse radar and other recent remote-sensing techniques look promising, but results have been inconsistent so far due to extreme variability in the physical properties of sea ice. It is therefore necessary to infer underwater shape from the upper surface character of each floe: viz. freeboard measurements. However, it is well known that the density of sea ice is not uniform through its thickness and will vary due to desalination as the floe ages. In fact, it is virtually impossible to predict the draft of a particular ice floe or its underwater shape from theoretical considerations alone, and one must turn to empirical results. Ackley and others (1976) provide a useful nomogram (Fig.8.) derived from data in the literature and from coring measurements collected during the pilot study in 1972 for the Arctic Ice Dynamics Joint Experiment. Given the effective freeboard (defined as the measured freeboard adjusted to take account of the lower density of the snow cover) at a point, Figure 8 enables us to calculate the thickness of the floe at that point and, after smoothing, to infer its underwater profile. In view of the results of the present study, it is essential that this process be carried out before the behaviour of individual ice floes in waves can be accurately modelled.

Fig. 8. Nomogram for calculating approximate ice thickness from effective freeboard. Shaded area indicates standard deviation between measured and predicted values (after Ackley and others 1976).

After the underwater shape of the ice floe has been found, the behaviour of a single floe in ocean waves may be evaluated. Each floe will respond to the incoming waves as discussed above, and will only allow some proportion of the incident energy to proceed into the pack, the rest being reflected back. The transmission and reflection properties of single floes will therefore determine the overall attenuation rate of ocean waves as they propagate into the MIZ. Using the simple two-dimensional scattering theory proposed by Wadhams (in press) and neglecting floe-floe interactions, the appropriate attenuation coefficients at all frequencies present may be found and hence the spectral character of the waves may be tracked into the pack. No satisfactory alternative to two-dimensional scattering theory exists at present although the author is now investigating the possibility of using the theory of wave propagation through random media as a more effective means of calculating ice-induced energy decay. It should be stressed, however, that the Wadhams theory gives quite acceptable agreement between theory and experiment for two-dimensional transects.