Introduction

Several observational studies suggest that long-period oscillations occur in ice shelves and ice tongues at periods well beyond those of the typical open ocean swell; for example, long-period oscillations were observed on an ice shelf during the “90° South” expedition (personal communication from M. Kristensen). Oscillations with 40–50 s period, for example, have been measured in the Erebus Glacier Tongue, McMurdo Sound, Antarctica. Further, there is some evidence to suggest that the waves are propagating (personal communication from W.U. Robinson), i.e. that they are not entirely due to free oscillations set up by some applied action which may only exist for a short time and is aperiodic. The questions which then remain are where floes the energy come from to drive these travelling waves, and how efficient is the coupling between the ocean and the ice shelf/tongue? In this paper we will examine the possibility of long-period surface gravity waves coupling between the ocean and one of those geometries, namely a massive ice shelf.

During the summer months in particular, meltwater from the ice shelf may overlie the denser sea water leading to a stratified water column and the possibility of internal gravity waves. For a given wave period, there are two travelling internal waves, one corresponding to an homogeneous water body with the average density and another in which the surface and internal waves are in antiphase. The wavelength of the latter type of internal wave is usually much less than the wavelength of the former. For example, for the periods and geometries that we will consider in this ptipcr, the wavelengths of the first type of wave are in the range 1 to 10 km, whereas the wavelength of the second type is less than 10 m, even with extreme stratification. Coupling to an ice shelf at such short wavelengths is unlikely to be significant. Hence, the short internal waves can be neglected as a forcing mechanism or, equivalently, the ocean can be regarded as being homogeneous.

The coupling question is investigated in this paper by using a thick elastic plate equation to study the transmission of wave energy into an ice shelf. The complete set of wave potentials is used, viz. travelling plus an infinite sum of evanescent waves in the open water, and travelling, damped-travelling, plus an infinite sum of evanescent waves in the ice-covered region. Unlike earlier approximate solutions, which omit the evanescent contribution (see Fox and Squire, Reference Fox and Squire1990, for references), this ensures continuity of the potential and the velocity over the depth, thereby solving the mathematical problem completely. Hence, the amplitude reflexion and transmission coefficients can be found to any desired precision. Further, the strain on the ice shelf induced by ocean waves or swell can be calculated from the velocity potential, once found.

By working backwards, the incident wave heights in the open water, required to induce flexural oscillations in the correct frequency band with the observed magnitude on the ice shelf, will then be found. There is little energy available in the open sea at the observed periods, and it is of interest to discover from where the input derives. Swells at such long periods are unlikely, for example, unless ice shelves and tongues react selectively to a band of long wave periods. This is an interesting possibility, as it suggests that existing accelcrometer-type instrumentation (as typically used in wave buoys) has insufficient resolution to detect such low frequency waves, whereas the band-pass filtering action of an ice shelf, combined with alternative curvature-based measurement techniques, makes detection possible.

Seiches, and harmonics of seiches, are other possible forcing mechanisms, as are nonlinear modulation, surf beats, edge waves, or wave packet formation, although these seem rather esoteric mechanisms.

Mathematical Model

We investigate the propagation of ocean waves into an ice shelf by solving the following mathematical model of the physical system. We will assume that the ocean is homogeneous and has constant depth, the incoming ocean waves have small amplitude and propagate normal to the edge of a massive (half plane) ice shelf of constant thickness. This model has no variation parallel to the edge of the ice shelf, and so is essentially two-dimensional. In cases where the thickness of the ice shelf is a small proportion of the ocean’s depth, we can think of the ice shelf as floating on the water. This geometry is shown in Figure 1.

Using the general linear theory for small amplitude waves, along with the usual assumption of irrotational flow, we define a velocity potential within the water, Φ(x, y, t), which satisfies Laplace’s equation at each instant in time. Initially, we will consider a single monochromatic ocean wave with time dependence e^iwt where the period of the wave is Τ = 2π/ω. As our model is linear, and invariant in time, it follows that Φ (x, y, t) = ϕ(x, y)e^iwt where the spatial dependence also satisfies Laplace’s equation; i.e.,

At the sea floor, the vertical component of velocity must vanish leading to the boundary condition

Fig. 1.

Schematic showing the semi-infinite sheet of sea ice of thickness h floating on water of depth H. Surface gravity waves travel from left to rigid and are reflected and transmitted at the sea-ice boundary. Note that the y-axis is vertically upwards with y = 0 at the sea floor and y = H at the surface; this is done to simplify subsequent algebra.

A boundary condition for the surface of the open sea can be established using the kinematic condition and Bernoulli’s theorem (Phillips, Reference Phillips1977) which in terms of ϕ is

where g is the acceleration due to gravity.

In the absence of viscous and other nonlinear terms, Equations (1) to (3) are the general equations governing the potential within an homogeneous constant-depth sea with a free surface. In this paper, we are concerned with long-period waves and later we will find that only small amplitudes are involved. Consequently, the nonlinear terms in a more complete formulation are negligible.

Thick-plate model of the ice shelf

The properties of the ice shelf will be modelled as a thick isotropic elastic plate. The most complete work to date on the problem of ocean coupling to floating ice sheets has modelled the ice as a thin elastic plate (Reference Evans, Davies and HobokenEvans and Davies, 1968; Reference Fox and SquireFox and Squire, 1990, Reference Fox and Squire1991; Reference Fox and SquireSquire and Fox, 1990). The thin-plate models neglect the effects of the finite thickness of the plate and, hence, neglect its rotary inertia and transverse shear deformation. An equation relating the transverse displacement, w, of a thick plate to the applied normal pressure, p, has been derived by Reference MindlinMindlin (1951) from the three-dimensional equations including the effects of rotary inertia and transverse shear with suitable assumptions such as the plate faces remaining parallel. This relation is

where ρ_i is the density of the ice shelf, h is the thickness of the ice shelf, μ, is a shear constant which is close to unity and, for a Young’s modulus Ε and Poisson’s ratio υ, G = Ε/2(1 + ν) is the shear modulus and L = Eh³/12(1 − ν² ) is the flexural rigidity of the plate. In order to determine the motion of the plate uniquely, boundary conditions must be applied at its edges. The natural boundary conditions can be found by comparing Equation (4) with the variational principle that the plate must satisfy (Hildebrand, Reference Hildebrand1965), which can in fact be derived from the differential equation. The two boundary conditions found this way are

where ∂/∂n and ∂/∂s denote differentiation normal and parallel to the boundary, respectively, and the Laplacian acts in the plane of the plate.

When the plate thickness, h, is small and when we are only concerned with waves with a long period, T, the terms in Equation (4) representing rotary inertia and transverse shear can be neglected, leaving the simpler thin-plate equation

The boundary conditions for this simpler case are given by neglecting the time derivative in Equation (5) and by Equation (6) as it stands. The case of ocean wave coupling to thin ice plates has been solved in detail by Fox and Squire (Reference Fox and Squire1990, Reference Fox and Squire1991) and Squire and Fox (Reference Fox and Squire1990).

Restriction of the thick-plate equations to the model

Assuming no cavitation between the ice and water, the thick-plate equation can be combined with the dynamic and kinematic conditions at the surface of the water to give a boundary condition for ϕ at that surface. When restricted to the two-dimensional monochromatic model, that condition can be written

where for brevity we have defined the constants

with ρ_w denoting the density of sea water and m = ρ_ih being the mass per unit area of the ice shelf. As an example of the relative magnitudes of the constants, for t = 50 s, H = 500 m and h = 200 m, we find L = 4.4 × 10¹⁵, A = −6.4 × 10⁷, Β = 7.1 × 10³, C = 1.7 × 10⁵ and D = −16. When the boundary conditions (Equations (5) and (6)) are similarly restricted they become

Equations (1), (2), (3), (8), (9) and (10) give a full statement of the mathematical model.

Modes of the system

A complete set of potentials can be found, each satisfying the differential Equation (1) and the bottom boundary condition (Equation (2)) and one of the surface conditions (Equations (3) and (8)), by separating variables on cither side of the line x = 0, 0 < y < H. These potentials are the modes of the system and any solution of the model can be represented as a combination of these modes. Each mode has the form ϕ = e^kxe^±iky. Pairs of these modes can be combined to give the equivalent set of modes

The boundary condition at the sea floor rules out the modes with sin(ky) y-dependence, leaving only those with cos(ky) y-dependence. Only certain values of k give rise to modes which also satisfy one of the boundary conditions at the sea surface.

In the region x < 0 the modes must satisfy the open sea condition (Equation(3)) which demands that k satisfies the dispersion equation

This equation has two pure imaginary roots ± ik_T (k_T > 0), corresponding to travelling waves propagating in the positive and negative x directions, and infinitely many real roots ±k_n (n = 1,2,…), where

nπ/H corresponding to evanescent modes.

In the ice-covered region, x > 0, we will write the modes as ϕ = e^kx cos(ky), replacing k by k to distinguish between the two regions. The wavenumber k must satisfy the more complicated dispersion equation

For physically relevant geometries and most periods this dispersion equation also has two imaginary roots ±ik _T(k _T = 0) corresponding to travelling waves, four complex roots occurring as plus and minus a complex conjugate pair ±k_D and

(k_D has positive real and imaginary parts) corresponding to damped travelling waves, and infinitely many real roots ±κ_n (n = 1,2,…), where in general (n - 1/2)π/H < k_n < (n+1/2)π/H, each giving an evanescent mode.

The evanescent modes have no phase term and, hence, carry no energy but are vital to the coupling of energy between the open sea and the ice shelf. Because of their e^kx (e^kx for beneath the ice) x-dependence, the physically relevant evanescent modes decay away from the ice edge. For large magnitudes of x, only the travelling waves are present, whereas near the ice edge the evanescent modes have significant amplitude. Thus, in the stationary case that we are considering, of a single period wave incident on an ice shelf, each evanescent mode traps a fixed amount of energy near the boundary between the ice shelf and the open sea. Each evanescent mode does not propagate any energy, but facilitates the transfer of energy from the incident wave to the reflected and transmitted waves.

We will consider only the roots of the dispersion equations giving rise to bounded modes and relevant to the problem of an ocean wave being reflected from, and transmitted into, the ice covered region. In particular, we do not allow any wave propagating from large positive x assuming that the transmitted wave is dissipated at, or before, the end of the ice shelf so that no significant wave is reflected there. We may use the relevant modes to write any bounded potential in the open sea region, x < 0, as

where I and R are the complex coefficients of the incident and reflected travelling waves, respectively, and {a_n } is the set of coefficients of the bounded evanescent modes. Similarly, any relevant potential in the ice covered region, x > 0, can be written

where T, b₊ , b₋ and {b_n } are the complex coefficients of the transmitted wave, the damped travelling waves propagating in the positive and negative x directions, and the infinity of evanescent modes, respectively. The potential in this region must also satisfy the ice-edge conditions (Equations (9) and (10)).

In reality, an ice sheet is not of infinite extent and the assumption that no energy is reflected from its far end may be questioned as the travelling wave suffers no attenuation in our model. With perfect reflection it would be a simple matter to include an additional term in Equation (14) representing the returning wave. An alternative approach, due to Gui and Squire (Reference Gui and Squire1989), assumes the ice to behave as a fixed-free beam and solves for the normal modes of the system when the beam is acted upon by a (white) wind pressure spectrum or by a realistic ocean wave spectrum.

Method of solution

A non-trivial solution is found by specifying unit incident amplitude and then matching the potentials ϕ_o and ϕ_i on the line x = 0. The natural quantities to equate are the potential and its derivative normal to the matching boundary, i.e. ∂ϕ/∂x. This can be achieved by minimising the error integral

evaluated at x = 0. The constant α is a weighting term included to improve convergence which we set empirically to 10 in most cases. The solution is thus found by minimising ε, subject to the conditions of unit incident amplitude and the boundary conditions at the ice edge given in Equations (9) and (10). The first of these conditions restricts ϕ_o , so that I = 1, whereas the second and third give rise to two linear equations relating the coefficients of ϕ_i in Equation (14). In order to make the calculations possible, the infinite sums in Equations (13) and (14) are terminated at some finite N; that is, only the first Ν evanescent modes are included. The matching is performed for a sequence of increasing Ν until the solution has converged to its final value. An independent check is made on the resulting potentials to ensure that they represent a system in energy balance (Fox and Squire, Reference Fox and Squire1990). Using the property that the error, ε, is a quadratic form in the real and imaginary parts of the coefficients I, R, {a_n }, T, b₊ , b₋ and {b_n }, the minimisation is performed, with respect to these coefficients, using a preconditioned linear conjugate gradient algorithm (Gill and others, Reference Gill, Murray and Wright1981) restricted to the sub-space in which the three extra conditions are satisfied.

Numerical Results

The geometries and constants

The following results are based on two water depths, H, of 1000 m and 500 m, and an ice-shelf thickness, h, of 200 m. The shelf thickness was chosen to approximate the average dimension of the Erebus Glacier Tongue which tapers from about 300 m at its hinge line to near 100 m at its terminus (Holdsworth and Holdsworth, Reference Holdsworth and Holdsworth.1978). The 500 m water depth was chosen to give a shallow water geometry analogous to the water depth in Erebus Bay, beneath the glacier tongue, which is seldom more than 400 m.

The density of sea water is taken to be ρ_w = 1025 kg m⁻³. Constants relating to the physical properties of the ice shelf are as follows. The Young’s modulus is Ε = 6 × 10⁹ N m⁻², Poisson’s ratio is ν = 0.3 and the density of the ice is taken to be ρ_i = 922.5 kg m⁻³. The value of the shear constant μ is chosen to ensure that the velocity of plane waves predicted by the thick-plate equation corresponds to the equivalent velocity predicted by the full three-dimensional theory (Mindlin, Reference Mindlin1951). This implies that μ satisfies the relations

For the value of υ above we find that μ = 0.9274.

Comparison of the thick-plate and thin-plate models

It was stated earlier that the thick-plate equation reduces to the simpler thin-plate equation in the thin ice and long wave period limits. The ranges of thickness and period over which this simplifying assumption can be made can be determined by comparing the roots of the dispersion equations resulting from use of the two plate equations. We will compare these formulations in terms of the travelling wavelengths for the geometry given above, over a range of periods.

The travelling wavelengths for the ice-covered region and the region of open sea, given by 2π/k_T and 2π/k_T respectively, are shown in Figure 2 for wave periods from 1 to 60s. Both the thick- and thin-plate wavelengths are shown and it can be seen that these coincide for long periods and diverge for periods less than, say, 10 s. In each case of thick-plate, thin-plate and open sea, the wavelengths for the 500 m depth are less than the wavelengths for 1000 m water depth at the same period. This is most easily seen for periods above 20 s. It is interesting to note that whereas 1000 m water is effectively deep

Fig. 2.

Wavelengths in metres plotted against period for travelling waves in the open sea and ice-covered region as given by both, the thick- and thin-plate models. Two geometries are represented, 200 m thick ice on water of 1000 m and 500 m depth. The wavelengths for the two depths coincide for short-period waves and are longer for the greater depth at long periods.

with respect to the ocean wave for periods less than 20 s, 1000 m of water is deep with respect to the ice-coupled wave only for the smallest period shown. Also note that 500 m of water is not deep with respect to t he ice-coupled wave for any period shown and is shallow for all waves of period 40 s and above. It is also noteworthy, and perhaps contrary to expectation, that the thick-plate model gives a wavelength closer to the open-sea wavelength than does the thin-plate model. This suggests that the thick plate is softer than the thin plate, and the effect may be due to the greater number of degrees of freedom associated with the thick plate (personal communication from K. Hutter).

The ratio of travelling wavenumbers for the thick- and thin-plate models is plotted in Figure 3 for 1000 m water depth and this shows more clearly the agreement and divergence of the two values. The analogous curve for 500 m water depth is essentially identical. Note that for periods above 10 s the two values differ by less than 2%, falling to less than 1% difference for periods above 20 s. The difference rises to 10% for periods less than 2 s, but for such short periods both models will predict the total reflection of an incoming wave.

A similar pattern of percentage difference, between the roots of the thin- and thick-plate dispersion equations, occurs for the complex root k_D and the first real root k₁ . The second real root, k₂ shows less dependence on which equation is used and the third and subsequent real roots (k_n , n = 3, 4,…) are essentially equal for the two equations over the entire range of periods. Thus, we expect that the simpler thin-plate model should correctly predict the thick-plate behaviour for periods of 20 s and above.

More usual as a measure of relative plate thickness is the ratio of plate thickness to the wavelength of the

Fig. 3.

Ratio of ice-coupled travelling wavelengths as predicted by the thick-plate and thin-plate models plotted against period for 200 m thick ice on water of 1000 m depth. The analogous curve for 500 m water depth is essentially identical.

travelling wave. This ratio is shown in Figure 4 for the same range of periods used above and for both water depths. For periods below approximately 8 s the plate thickness is greater than one tenth of the wavelength, and it is for that range of periods that the wavelengths differ significantly. This observation agrees with that of Mindlin (Reference Mindlin1951), who reached the same conclusion on the basis of the wave velocities.

It is also worth noting here that a scale analysis of Equation (4) may be carried out, choosing physically reasonable ranges for its various parameters, and for the period and travelling wavelength. When this is done, it is found that the most significant terms in Equation (4) are precisely those of Equation (7), thus confirming the conclusions reached above. It was, however, necessary

Fig. 4.

Plate thickness divided by travelling wavelength as given by the thick-plate model plotted against period for 200 m thick ice on water of 1000 m depth and 500 m depth.

Fig. 5.

Reflection and transmission coefficients plotted against period for 200 m thickice on water of 1000 m and 500 m depth found using the thick-plate model.

to go through the full procedure because of the presence of the evanescent wave modes. These modes have wavelengths which approach zero as n → ∞ and, hence, could be altered markedly by the additional terms of Equation (4). In the event they were not.

Reflection and transmission characteristics

The reflection and transmission coefficients for the amplitude of surface displacement, denoted R and Τ respectively, can be calculated from the coefficients of potential I, R and Τ using the following relationships

These coefficients are plotted over the range of periods 10 to 140 s in Figure 5. Note that the transition from total reflection to total transmission occurs around a wave period of 60s. Incoming ocean waves with periods less than 20 s are completely reflected for both water depths, whereas those with periods somewhat longer than 140 s are fully transmitted in the 1000 m water depth case. Full transmission in the 500 m water depth case occurs at still longer periods. Note that the transmission coefficient is generally smaller when the water is 500 m deep as opposed to 1000 m deep. However, the reflection coefficient for 500 m deep water is appreciably larger than that for 1000 m deep water only for periods above 50 s, whereas that relationship is reversed for periods below 50 s. Essentially identical results are predicted when the thin-plate equations are employed.

Ocean wave-induced strain in an ice shelf

Measurements of flexural gravity waves at the surface of an ice shelf are more usually made in terms of the surface strain rather than the transverse displacement. The strain in a plate is simply half its thickness divided by the radius of curvature which, for small amplitude waves, reduces to

Fig. 6.

Magnitude of strain at the surface of an ice shelf plotted against distance from the edge of the ice in metres for wave periods of 40, 50, 60, 70, 80, 90 and 100s. The geometry is 200 m thick ice on water of 1000 m depth and the incoming wave amplitude is 1 m. The thick-plate model was used.

Thus, the strain is the product of the amplitude of a wave and the curvature produced by a unit amplitude wave of the same period. For a fixed wave height in the open sea, increasing period will increase the transmission coefficient and, hence, the amplitude at that period in the ice. At the same time, however, increasing period will decrease the curvature of a unit height wave. It is clear, therefore, that the maximum strain at the surface of the ice shelf will be generated, for fixed amplitude ocean waves, when the period is in the range between total reflection and total transmission.

The magnitude of strain generated by 1 m amplitude ocean waves with periods of 40 to 100 s, in 10 s intervals, is plotted in 6 and 7 for water depths of 1000

Fig. 7.

Magnitude of strain for the same geometry and periods as Figure 6 except with 500 m deep water.

and 500 m, respectively, as a function of distance from the edge of the ice shelf. Each curve shows some characteristic features. In each case there is no strain at the ice edge: indeed, this is the boundary condition in Equation (10). After a transition, the magnitude of strain reaches a far field value, due solely to the transmitted wave. The strain in the intermediate range is affected by the strain due to the evanescent and damped travelling modes as well as the travelling wave. For some wave periods, e.g. 50 s, the magnitude of strain reaches a peak, at an intermediate distance, which is greater than the far field value.

Note that the maximum magnitude of strain is exhibited for 60 and 70s wave periods with 1000 and 500 m deep water, respectively. Indeed, for 1000 m deep water, the 60 s period generates the greatest strain for distances up to 3000 m from the edge and close to the greatest strain for greater distances. For 500 m deep water, waves with a 70 s period generate close to the greatest strain for distances up to 2000 m from the edge and the greatest strain for greater distances. This implies that the maximum coupling between ocean waves, or swell, and strain in the ice shelf occurs at 60 s for water depths of 1000 m and at 70s for depths of 500 m, at least for 200 m thick ice. It is interesting to note that the strain generated in 200 m thick ice, at each period, is greater for the shallower water geometry despite the fact that the transmission coefficient is lesser for that case. This is evidently due to the corresponding decrease in wavelength and, hence, increase in curvature of an ice-coupled wave, of a given period, when the water depth is decreased. This more than compensates for the decrease in transmission.

Figures 6 and 7 show that the strain in an ice shelf should be a sensitive detector of ocean waves with periods of 50–80 s. There have been some measurements of long-period oscillations in an ice tongue against which we can calibrate the strains predicted in Figure 7. Whereas the geometry of an ice tongue docs not conform exactly to the model analyzed, we use these measurements in the absence of more appropriate data. The authors would like to obtain more testing data against which the theory can be compared and would welcome input on this point. Propagating oscillations with periods of 40–50 s and with a strain amplitude of approximately 5 × 10⁻⁸ were observed in the Erebus Glacier Tongue in the austral summer of 1989 (personal communication from W.H. Robinson). As a rough calculation, we can use the coupling calculated above to infer the size of seas required to drive such an oscillation. With reference to the curve for 50 s in Figure 7, we see that a 50 s ocean wave would only need an amplitude of about 1 mm to generate such a strain. Given the difficulty of measuring long-period ocean waves, it is not improbable that long-period ocean waves of this amplitude could exist, even though none have been reported. The greatest strain observed in the Erebus Glacier Tongue occurred the year before the observation above and was of the order of 5 × 10⁻⁷ strain (personal communication from W.H. Robinson). Again on the basis of our crude comparison, this would imply that 50 s seas rarely go above 1 cm in amplitude.

The Gui and Squire (1980) analysis mentioned earlier reaches similar conclusions. It is likely that reality lies somewhere between the two extremes, with the strain field due to the ice-tongue oscillation being partitioned into a propagating component travelling up the tongue together with a standing vibration.

Conclusions

The results of this paper have implications for the range of periods within which we might expect to measure strain on an ice shelf, and the suitability of the thin-plate model for calculating that strain. Some conclusions are summarised below:

• • The simpler thin-plate model coincides with the thick-plate model for periods above 20s. Below that period both models predict total reflection of incoming waves. Hence, in terms of the reflection and transmission coefficients, the two models give essentially identical results for all periods and for the two geometries considered. This conclusion may well hold in general, that is, for a range of ice thicknesses and water depths.

• • For the geometries modelled, the maximum coupling between ocean waves and strain in the ice occurs for a 60 to 70s wave period. This is nearly true for all distances less than far field and in the far field. Further, the ice shelf acts as a band-pass filter, reacting most strongly to periods around 60 to 70s. It is likely all geometries react similarly, with a band-pass response to incoming wave periods.

• • Ocean wave driving of long wave period strain in an ice shelf is a plausible mechanism. The comparisons made with existing data indicate that mean wave-heights of 1 mm at a 50 s period would account for the observed strains.

Acknowledgements

The authors are indebted to the New Zealand Foundation for Research, Science, and Technology, the University of Otago and the New Zealand University Grants Committee for their continued support. Colin Fox was a University Grants Committee Post-Doctoral Research Fellow when this work was carried out. The helpful comments of the reviewers and editor are acknowledged.