## Abstract

Flexural oscillations of floating sea ice sheets induced by ocean waves travelling at the boundary between the ice and the water below can propagate great distances. But, by virtue of scattering, changes of ice thickness and other properties encountered during the journey affect their passage, notwithstanding attenuation arising from several other naturally occurring agencies. We describe here a two-dimensional model that can simulate wave scattering by long (approx. 50 km) stretches of inelastic sea ice, the goal being to replicate heterogeneity accurately while also assimilating supplementary processes that lead to energy loss in sea ice at scales that are amenable to experimental validation. In work concerned with scattering from solitary or juxtaposed stylized features in the sea ice canopy, reflection and transmission coefficients are commonly used to quantify scattering, but on this occasion, we use the attenuation coefficient as we consider that it provides a more helpful description when dealing with long sequences of adjoining scatterers. Results show that scattering and viscosity both induce exponential decay and we observe three distinct regimes: (i) low period, where scattering dominates, (ii) high period, where viscosity dominates, and (iii) a transition regime. Each regime’s period range depends on the sea ice properties including viscosity, which must be included for the correct identification of decay rate.

## 1. Introduction

Swell waves generated in the distant major oceans or created by intense polar storms whose paths encircle the North Polar Basin penetrate deep into the interior of the ice-covered Arctic Ocean as flexural-gravity waves, so named because both the bending of the sea ice canopy and the inertia of the water column beneath affect their dispersion. These waves are also called ice-coupled waves. The goal of the work described in this paper is to replicate the passage of such waves over substantial lengthscales, a particularly demanding task because sea ice is a highly heterogeneous material in regard to its physical topography. Over a few kilometres, one would typically encounter open and refrozen leads, cracks, pressure ridges and gradual or abrupt changes of thickness, for each of which the growth and development phase and the history of deformation may be attributed to meteorological and oceanographical factors. Furthermore, inhomogeneity in this situation is particularly awkward to deal with, as thickness variation, for example, markedly alters the way the floating sheet bends in response to the wave and, through this, the dispersion, so its effect is more significant than one would imagine at first sight.

It was Greenhill (1916) who first proposed a standard mathematical treatment for flexural-gravity waves. He conjectured that the ice behaves as a uniform thin elastic (Euler–Bernoulli) plate and employed linear wave theory, which requires the ice-coupled waves to be small in amplitude. These assumptions are, in fact, quite reasonable, as over small spatial scales the strain rates involved are accordant with perfect elastic deformation and the waves of the Arctic interior are almost always of low amplitude. Many papers have appeared since that have followed Greenhill’s lead and used the Euler–Bernoulli plate as the basic boundary condition for the sea ice; these are documented in two review papers (Squire *et al*. 1995; Squire 2007) along with the historical progression of advances towards more realistic sea ice that appears in the later review.

In the present paper, two interrelated notions are propounded and explored. The first relates to the incorporation of viscosity into the plate equation that describes sea ice bending (Squire & Fox 1992; Vaughan & Squire 2008*b*). In the current context, this introduces exponential decay into an otherwise sinusoidal, propagating ice-coupled wave train. It recognizes the inherent inaccuracies in the perfectly elastic Euler–Bernoulli plate boundary condition when waves are propagating over several wavelengths and includes viscosity to interfuse and absorb the many processes that remove energy from the propagating wave when fetch is significant. The inclusion of viscosity by Vaughan & Squire (2008*b*) was also partly motivated by the desire to reduce the effects of length resonance, which is an artefact of the construction of their model and the two dimensionality of the problem. They continued to use the Euler–Bernoulli plate equation but included an additional, velocity-dependent viscosity term (Robinson & Palmer 1990; Banks & Inman 1991).

The second notion concerns scattering. The leads, polynyas, cracks and pressure ridges that cause sea ice to be heterogeneous partially transmit and partially reflect the incoming waves. To a degree, this has been investigated in the past by simulating how isolated features such as these scatter the incoming wave train, for example, by cracks (Squire & Dixon 2000), leads (Chung & Linton 2005), ridges (Williams & Squire 2004) and by sequences of such features. However, natural sea ice is rarely this simple and this observation led Vaughan & Squire (2008*a*) to simulate scattering by sea ice transects sampled from data collected using upward-pointing sonar during submarine voyages (NSIDC 2006). Those authors assimilated measured natural variations by assuming that sea ice thickness was piecewise constant but their solution method did not allow Archimedean draught to be accommodated and transects were relatively short in length (≤1.2 km but dependent on the computer employed). Vaughan & Squire (2008*a*) confirmed that flexural-gravity waves decayed exponentially, but their short transects prevented wave periods of real interest from being examined. Kohout & Meylan (2008) report a model of the marginal ice zone in which scattering may be simulated for up to 4000 ice floes. Their model is not directly comparable with the one described here because they use free edge conditions, but some of their observations are relevant. For example, they note that the wave energy decays exponentially in their results and that decay rates are under predicted at high periods.

By using the multi-mode approximation of Bennetts *et al*. (2007), this paper takes a new path. Instead of modelling flexural-gravity wave scattering due to a very small region of continuously varying thickness, we accommodate very long stretches of sea ice that have piecewise constant thickness. To do this, we again use transects sampled from the NSIDC (2006) data. These data are presented as lists of ice draughts at discrete points, but measurement and processing limitations require that each draught should be considered as being a representative average value so that the reconstruction of an ice sheet draught profile necessitates making an assumption about the nature of the sea ice between draughts. By this means, we create a piecewise constant profile, an approximation motivated primarily by numerical expediency that is probably the best that can be done with the NSIDC dataset. In summary, we describe in this paper a method that allows the correct physical draught to be incorporated and long transects of up to 50 km to be considered.

While focused here on sea ice in nature, the model and methods developed in this paper may be equally useful in helping to interpret data from laboratory experiments. These could either be done in a freezable wave flume such as the Arctic Environmental Test Basin at the Hamburgische Schiffbau-Versuchsanstalt or a tank at the US Army Cold Regions Research and Engineering Laboratory, or they could exploit synthetic sea ice fabricated from a suitable plastic like those done at Iwate University in Japan by Sakai and co-workers (Kohout & Meylan 2008) or planned by the authors and colleagues at the Laboratoire de Mécanique des Fluides, Ecole Centrale de Nantes.

For an isolated thickness change or a short sequence of commensurate features, scattering of ice-coupled waves has customarily been quantified by means of the (asymptotic) reflection and transmission coefficients on the near and far side of the transect being modelled. These parameters necessitate the presence of long, uniform sea ice on either side of the inhomogeneity. Sea ice in the real world is not so obligingly constructed, and this makes identification of these parameters in field measurements ambiguous. Vaughan & Squire (2008*b*) investigated using wave decay as a surrogate for wave scattering, but were again limited by the model that they employed. They also observed that the wave elevation envelopes varied greatly between any two transects and sought a method to characterize decay via the statistical properties of the transect. This is also the approach used here.

In the following section, a brief description of the mathematical theory underpinning this study is provided. Following this, in §3, model results are validated against other independent models for simple test cases, thereby providing compelling evidence that the model operates correctly and accurately. Finally, the effect of a long stretch of natural heterogeneous sea ice on flexural gravity wave trains at different periods is considered in §4, where results are presented for NSIDC Arctic ice draught data.

## 2. Basic theory

This section begins with a discussion of the equations that govern the problem (§2*a*). The geometry of the sea ice is introduced (§2*b*), along with an appropriate terminology and the conditions it imposes. The solution requires that the roots of the dispersion relation be found. In this respect, our method is no different from the majority of existing methods, but here the presence of viscosity and the need for computational efficiency require special consideration, which is discussed in §2*c*. The solution method is described (§2*d*), followed by a technical note on how performance may be improved (§2*e*). Here, the main body of results are simulations of scattering and damping in real sea ice, with transects sampled from the files of draught data that are archived by the National Snow and Ice Data Center (NSIDC 2006). Issues pertaining to the use of these data are described in §2*f*.

### (a) Governing equations

The situation we envisage is that of a long-crested flexural-gravity wave travelling through a large expanse of sea ice and being attenuated by the viscosity in the ice. In addition to attenuation due to viscosity, the travelling wave is scattered by variations in the thickness and draught of the ice. The mathematical model that we design to represent this situation is restricted so that these geometrical variations may only occur in one horizontal dimension, which we define through the use of the Cartesian coordinate *x*, and the incident wave propagates in this direction. Consequently, there is no variation perpendicular to *x* in the horizontal plane and it is only necessary to consider a two-dimensional cross section parallel to the *x*-axis in order to investigate the scattering process. The vertical plane is defined by the coordinate *z*, which is directed upwards and has its origin, *z*=0, set to coincide with the equilibrium fluid surface in the absence of ice cover.

The two-dimensional geometry that we consider is of an infinite extent, that is, , and a fluid domain of finite depth occupies this entire interval. We suppose the seabed to be fixed and impermeable and, because it is considered to be sufficiently deep (approx. 1 km) that its undulations may be neglected in the scattering process, to be flat and lying in the plane *z*=−*h*. An ice sheet rests on the surface of the fluid, with the undisturbed interface between the fluid and ice denoted by *z*=−*d*(*x*) and the thickness of the plate described by the function *D*(*x*).

We will seek the motion of the fluid, which is assumed to be inviscid, incompressible, homogeneous and in irrotational motion, through a velocity potential from which the velocity field may be recovered via . By assuming a harmonic time dependence of given angular frequency *ω*, we write , where *g*≈9.81 m s^{−2} is the acceleration due to gravity. The unknown function *ϕ* is the (reduced) velocity potential, which satisfies Laplace’s equation within the fluid domain, that is,
2.1a
and on the bed it must also obey the no-flow condition
2.1b

It is standard practice to model sea ice as a thin-elastic plate, and this allows us to deduce all of the stresses it experiences from knowledge of the small-scale flexural oscillations induced by fluid motion at the fluid–ice interface. This moving boundary is denoted *z*=−*d*(*x*)+Re{*W*(*x*)e^{−iωt}}, where *W* is a displacement function that must be calculated during the solution process. Assuming that the amplitudes of the oscillations are sufficiently small for linear theory to be applied, the displacement function is coupled to the velocity potential through the interfacial conditions
2.2a
2.2b
which are applied at the linearized boundary *z*=−*d*. Here, a prime denotes differentiation with respect to *x*, *κ*=*ω*^{2}/*g* is a frequency parameter, *μ*=*μ*(*x*)=*ρ*_{i}*D*(*x*)/*ρ*_{w} is the scaled mass of the ice and *β*=*β*(*x*)=*E**D*(*x*)^{3}/12(1−ν^{2})*ρ*_{w}*g* is the scaled flexural rigidity of the ice. Additional quantities have been introduced in the definitions of *μ* and *β* and these are: the density of the ice *ρ*_{i}=922.5 kg m^{−3}; the density of the fluid *ρ*_{w}=1025 kg m^{−3}; Poisson’s ratio ν=0.3; and Young’s modulus for sea ice *E*=5×10^{9} Pa.

The derivation of equation (2.2a) is based on a combination of the linearized Bernoulli equation and thin-plate theory to describe the oscillations owing to imbalances between the pressure exerted at its upper and lower surfaces caused by fluid motion. Equation (2.2b) is a kinematic condition that assumes that contact is maintained between the fluid and the underside of the ice during motion.

The term −i*γ**W* in equation (2.2a) is a modification that adds artificial viscosity (Robinson & Palmer 1990; Banks & Inman 1991; Vaughan & Squire 2008*a*), which we consider to provide a first-order parametrization of a number of dissipative mechanisms. Here, the quantity *γ*=*ω**Γ*/*ρ*_{w}*g*, where *Γ* has units of dynamic viscosity per unit length (Pa s m^{−1}). Any value *Γ*≥0 is acceptable, but for ice-coupled waves it seems that *Γ*=4×10^{3} Pa s m^{−1} is a reasonable upper limit (Squire & Fox 1992). During our investigation, we will take *Γ* to be an arbitrary non-negative constant and, as part of our numerical results, we will look at the effects of changing the level of viscosity in the ice sheet.

### (b) The geometry of the sea ice

In §4, we will use our model with the NSIDC Arctic ice draught data (NSIDC 2006). The bottom of sea ice is typically rough and it is therefore reasonable to use this dataset to represent the thickness and draught of the ice as being piecewise constant. Although it is possible to accommodate the thickness and draught of the ice plate as continuous functions in our solution, taking advantage of this constraint constitutes a considerable computational gain. Furthermore, we note that the errors introduced by approximating a smooth surface by a piecewise constant one have been investigated by Vaughan & Squire (2006), and it was found that the computational efficiency of the method far outweighs the negligible errors introduced by the approximation.

We therefore divide the *x*-axis into *s* sections *S*_{j}={*x*:*x*_{j−1}<*x*<*x*_{j}} (*j*=1,…,*s*), where the values *x*_{j} (*j*=0,…,*s*) are arbitrary but fixed. Within each section, the geometry of the ice is uniform and we write
2.3
to denote the changing values. The long, yet finite, interval occupied by the sections, namely , is referred to as a transect, and it is bordered by two additional sections, both of a semi-infinite extent. These occupy the intervals and and have uniform thicknesses *D*_{0} and *D*_{s+1}, respectively.

For each *x*_{j} that defines the position of one of the geometrical discontinuities in the ice that have just been introduced, governing equations 2.1 and (2.2) cannot be formed and we must, instead, apply alternative criteria to define the fluid and plate motions. First, at all points, the continuity of the fluid’s pressure and velocity must hold, and therefore
2.4a
where 〈〈·〉〉 denotes the jump in the included quantity across *x* over the fluid depth. The ice must also obey certain dynamic conditions that are applied through the displacement function. These are
2.4b
where 〈·〉 denotes the jump in the included quantity at a point *x*. These conditions represent the continuity of the plate’s displacement, slope of displacement, bending moment and shearing stress, respectively. When a draught changes between sections, we also have the condition of no-horizontal flow through the protruding edge, that is,
2.4c
where , and a superscript (±) is used on a horizontal point hereinafter to denote that a quantity is evaluated from the limit in which *d*=*d*_{±}.

The problem now only needs conditions that describe the way in which waves radiate away from the transect, far into the surrounding semi-infinite intervals of ice-covered fluid. These radiation conditions are
2.5
where *k*_{0} is the primary wavenumber that will be defined shortly. In the far-field *x*≪*x*_{0}, we have a partial-standing wave that consists of the incident wave and a reflected wave of amplitude *R*, whereas in the far-field *x*≫*x*_{s}, only a transmitted wave of amplitude *T* is present. An expression analogous to equation (2.5) holds for the displacement function *W*.

Our motivation is to investigate large-scale transects, and this necessitates the use of a method of a low numerical cost. A suitable approach is to restrict the vertical motion of the velocity potential *ϕ* to a finite-dimensional space spanned by the vertical modes (*n*=0,…,*N*), where *N* is a finite number chosen to give suitable accuracy. We therefore write
2.6
and seek the functions *ϕ*_{n} that control the horizontal motion. The wavenumbers *k*_{n} that define the vertical modes on which this approximation is based are discussed in the following section.

### (c) The dispersion relation

In order to define approximation (2.6), it is necessary to calculate the wavenumbers *k*_{n}, which are defined as the roots *k* of the ice-covered dispersion relation
2.7
where *H*=*h*−*d* is the fluid depth. We note that if *k* is a root of equation (2.7) then −*k* is also a root and, because they define the same vertical modes, we may restrict our attention to the upper half of the complex plane, including the positive real axis.

It is clear that the properties of the dispersion relation depend on the ice thickness, through the mass and rigidity, *μ* and *β*, respectively, and its draught, through the fluid depth *H*. We are therefore required to calculate a separate set of wavenumbers in each section *S*_{j} ( *j*=0,…,*s*+1), and we write *k*_{n}=*k*_{n}(*S*_{j})≡*k*_{n,j} (*n*=0,…,*N*) when it is pertinent to distinguish this fact. Consequently, we have partitioned our approximation into the various intervals defined by the sections and this will have important ramifications in the solution procedure.

#### (i) Zero viscosity

For *γ*=0, the primary root, *k*_{0}, is positive and real and the corresponding vertical mode, *χ*_{0}, supports propagating waves. This is always the first mode used, which ensures that the form of the radiation conditions for the full-linear problem (2.5) can be represented exactly in the approximation. Specifically, this means that in the far-field
2.8
and all of the subsequent functions *ϕ*_{n}∼0 (*n*=1,…,*N*) for *x*≪*x*_{0} and *x*≫*x*_{s}. We note that the only approximation in equation (2.8) of the full-linear solution is through the reflection and transmission coefficients *R* and *T*.

On the imaginary axis, there are an infinite number of purely imaginary roots *k*_{n} (*n*=1,…,), ordered so that |*k*_{n}|<|*k*_{n+1}|, and the modes that they define support increasingly evanescent waves. These motions are generated when a plane wave is scattered by a jump in the geometry between sections, and if the distance between the sections is relatively short, their interactions will be important.

The dispersion relation (2.7) also has two complex roots, *k*_{−1} and *k*_{−2}. They typically have both real and imaginary components and, if so, they are reflections of one another in the imaginary axis, , and their modes support oscillatory-evanescent waves. Under certain circumstances (set out in Bennetts 2007), these roots instead sit on the imaginary axis and produce modes that support evanescent waves.

Bennetts *et al*. (2007) proved that the modes defined by the roots *k*_{−1} and *k*_{−2} could be written as linear combinations of the modes *χ*_{n} (*n*=0,…,). The space spanned by the set {*χ*_{n}:*n*=0,…,} is therefore unaffected if it is augmented by the modes given by the complex roots and, for this reason, these modes are omitted from the expansion on which we base our approximation (2.6). Instead, the waves that they support are re-distributed among the included modes, as will be outlined shortly.

The roots of equation (2.7) are found using the Newton–Raphson method, with appropriately chosen starting values. In deep water, and the dispersion relation reduces to a quintic for which the roots may be found readily, providing good starting values for the finite depth case when seeking the primary root. The procedure for finding the imaginary roots is discussed by Bennetts *et al*. (2007). Essentially, it can be shown that there is an interval within which a root will lie, and knowledge of this interval allows starting values to be selected.

#### (ii) Non-zero viscosity

When viscosity is introduced through choosing *γ*>0, all of the above roots become perturbed so that the real root and imaginary roots are shifted into the first quadrant of the complex plane and we retain the notation used for each root by tracing it back to the *γ*=0 case. In particular, as the primary wavenumber *k*_{0} now possesses an imaginary component, the waves that it defines will decay as they propagate. Similarly, the evanescent waves, defined by *k*_{n} (*n*=1,…,*N*), will oscillate as they decay. For most values of viscosity that are relevant to the physical problem that we are interested in, the distribution of the roots is only a simple shift away from the *γ*=0 case discussed earlier, and it is straightforward to associate these roots with their no-viscosity counterparts. The behaviour of the roots as functions of *γ* is not always trivial but a discussion of the technical issues that arise here is not commensurate with the objectives of the current work.

The process for finding the primary root is as for the *γ*=0 case, but the presence of the viscosity means the deep-water dispersion relation has a complex coefficient. The technique employed to find the imaginary roots when *γ*=0 does not work for *γ*>0. Instead, we locate the corresponding root for *γ*=0 and use this as a starting value with *γ* increased, repeating the process until *k*_{n} is found for the desired value of *γ*. In practice, 15 intermediary steps prove to be enough to find the roots for even large values of *γ* at short periods.

#### (iii) Kheysin’s idea for finding ice thickness

It is pertinent at this point to make a brief discussion of Kheysin’s idea. While it lies outside the theme of the work, this discussion is useful in illustrating the frequency response of the viscous plate model that we have employed.

In figure 1, the primary root of the dispersion relation is plotted for many values of the period with *Γ*=30 Pa s m^{−1} and an ice thickness of 1 m. The imaginary part of the root is precisely the decay rate that the viscosity engenders. Through further numerical tests, we have found that for other values of viscosity the curve is shifted up or down but retains the same shape. Note that the peak at 9.2 s suggests that it is the most damped period: it shifts slightly as viscosity is varied but remains between 9 and 10 s.

Kheysin’s idea is that certain wavenumbers propagate with less attenuation (Nagurny *et al.* 1994), these wavenumbers being *k*_{res}, when the phase speed for a wave in open water matches that for a free ice plate, and *k*_{min}, when the phase speed for an ice-coupled wave has a minimum. Wadhams & Doble (2009) note that no mechanism has been described that explains these prefered wavenumbers. The two wavenumbers required by Kheysin’s idea are plotted in figure 1, but clearly there are no corresponding local minima in the primary roots, so Kheysin’s idea will not work for a system using the plate model that we use.

The idea rests on the assumption that these special wavenumbers will experience less attenuation and so will remain in wave spectra for longer, but the justification for this assumption is not apparent in the literature. In addition to scattering, mechanisms for attenuation in ice sheets include viscosity of the ice and water and turbulence. Vaughan & Squire (2007, 2008*a*) did not observe such a response in simulations of scattering by midlength transects. The viscous plate model that we have used provides a basic parametrization that introduces attenuation, but it would be a tenuous logic that inferred physical meaning and attributed the parameter to one of the aforementioned mechanisms.

Interestingly, in their examination of field data for the attenuation of ocean waves, Wadhams *et al.* (1988) presented results that are broadly reproduced here in figure 1, albeit serendipitously. In particular, their results show a maximum in the attenuation coefficient at periods of 7–10 s. They called this effect ‘rollover’ and suggested explanations for the effect (short-wave generation in polynyas and leads or by nonlinear transfer of energy from long waves). It is also possible that the rollover can be explained as the frequency response in a viscous system: strain rates in the ice are low for long waves so damping is minor, while for short waves oscillations are too fast for damping to be significant, and for intermediate waves there must be a maximum.

### (d) Solution method

In order to generate an approximation from expression (2.6), we make use of a functional whose stationary point coincides with the solution to the full-linear problem outlined in §2*a*. By restricting the space over which the stationary point is sought to functions of the form of the approximation given in equation (2.6) alone, the vertical motion is averaged out, and a set of coupled ODEs is produced in the horizontal plane that defines the functions *ϕ*_{n} (*n*=0,…,*N*). Despite it having no vertical dependence of its own, an approximation of the displacement function, *w*(*x*)≈*W*(*x*), is also obtained, through its relation to the velocity potential.

This method was devised by Bennetts *et al*. (2007) as a means of solving for transects in which the ice thickness varies continuously over a finite interval owing to variations on its upper and lower surfaces. However, it is equally applicable in the current piecewise uniform context. The use of vertical averaging still clearly defines it from an eigenfunction-matching method, for example, and has been shown to enhance the efficiency of the approximation (see Bennetts 2007).

The subspace over which we search for our approximation may be expanded to any size, and thus our approximation improved by choosing larger values of the dimension *N*. Incrementing the number of evanescent modes used in the approximation therefore creates a sequence of increasingly accurate approximations. As such, the full-linear solution may be obtained to any desired level of accuracy by using this method. The choice to base the approximation on the natural eigenfunctions of the problem, *χ*_{n}, has been shown to produce high accuracy at a low computational cost (see Bennetts 2007), and this is of fundamental importance to our goal of solving for long transects.

For each section of piecewise uniform geometry, the ODEs that govern the approximation may be solved by standard methods (for details, see Bennetts *et al.* 2007). We therefore have analytic forms available for the horizontal components of the velocity potential, which are most conveniently defined using the vector notation ** Φ**=(

*ϕ*

_{0},…,

*ϕ*

_{N})

^{T}, and we write

**(**

*Φ**x*)≡

*Φ*_{j}(

*x*) for

*x*∈

*S*

_{j}, where 2.9 in which the wavenumbers used in the approximation have been collected into the matrix

*K*

_{j}=diag{

*k*

_{0,j},…,

*k*

_{N,j}}, and a matrix of corresponding waves is given by e

^{±iKjx}=diag{e

^{±ik0,jx},…,e

^{±ikN,jx}}.

The role of the scalars and the vectors **v**_{u}=**v**_{u}(*S*_{j})≡**v**_{u,j} for *u*∈*U*={−1,−2} is in approximating the motions supported by the complex modes, as well as that of the evanescent waves absent from the approximation, and it is these terms that distinguish our vertically averaged approximation from a direct eigenfunction expansion. The wavenumbers are defined as the roots of the quartic equation
2.10a
that lie in the upper-half complex plane, and their weightings **v**_{u} by
2.10b
where the diagonal matrices *C* and *S* are and .

Here, *I* is the identity matrix of dimension *N*+1, and **1**=(1,…,1)^{T} and **0**=(0,…,0)^{T} are vectors of length *N*+1. Also introduced in equations (2.10a) and (2.10b) is the matrix *A* whose entries are .

The amplitudes of the rightward travelling waves are contained in the scalars *r*_{u}≡*r*_{u,j} and vectors **r**≡**r**_{j}, and similarly for the leftward travelling waves in *l*_{u}≡*l*_{u,j} and **l**≡**l**_{j} (the latter is distinct from the vector **1**). Note that we have extended the definition of the direction of a wave so that it includes the evanescent waves that decay in that direction. Furthermore, as the amplitudes differ between sections, we can clearly observe the partitioning of the solution introduced by approximation (2.6).

At the interface between two sections, the approximate velocity potential defined by equations (2.6) and (2.9) must be linked. It is not, however, possible to directly apply the fluid continuities obeyed by the full-linear solution (2.4a) to the approximate potential owing to the discontinuities in the vertical modes over these points. Instead, we have the jump conditions
2.11
where
2.12
is a matrix that is evaluated from either side of the interface between adjacent sections and has entries consisting of integrals over the fluid depth of the vertical modes from that limit, multiplied by the vertical modes from the side with the shallower fluid depth. These jump conditions are derived through the variational principle (see Bennetts *et al.* 2007), and they are approximate versions of the conditions (2.4a) and (2.4c) that allow for discontinuities in the velocity potential but ensure that these conditions are regained as the approximation converges to the full-linear solution.

An expression equivalent to equation (2.9) may be found for the approximate displacement. This is written as *w*(*x*)≡*w*_{j}(*x*) for *x*∈*S*_{j} ( *j*=0,…,*s*+1), where
2.13
As the displacement is only indirectly approximated, the dynamical conditions it must obey between sections, given in equation (2.4b), are retained in the approximation.

It is possible to relate the amplitudes in adjacent sections by imposing jump conditions (2.11) and continuity conditions
2.14
to expressions (2.9) and (2.13). Straightforward manipulations can then be used to write these relations in terms of the transfer matrices ( *j*=0,…,*s*) that give the waves on the right-hand side of an intersection in terms of those on the left-hand side, with
2.15
where the vector amalgamates all of the rightward travelling waves in a section and likewise for the leftward travelling waves in .

Although transfer matrices may be used to obtain the unknown amplitudes in an iterative manner (see Bennetts *et al.* 2009), which enables the solution for a transect appended by an additional section to be obtained from that of the original transect, our primary intention is to investigate the effects of changes in the properties of the ice in long transects of fixed length. The most computationally efficient way of achieving this was found to be by solving the system of equations provided by relations (2.15) for all sections simultaneously. Specifically, from equation (2.15), we may derive the block tri-diagonal system
2.16a
2.16b
where the matrices are of size (*N*+3)-square, and are defined by
2.17
and
2.18
This system of equations is to be solved for the unknown vectors and ( *j*=0,…,*s*), with the two remaining vectors, and , deduced from the radiation conditions (2.8) to be **I**_{1}=(1,0,…,0)^{T} and **0**, respectively. Note that, here and in the following, we have redefined **0** to be of length *N*+3 (previously it was shorter) and the identity matrix *I* to have size *N*+3. That is due to the absence of leftward travelling waves as *x*≫*x*_{s} and is used in the second part of equation (2.16a), whereas contains the amplitude of the incident wave and forces the system through its appearance in the former part of equation (2.16a).

By inverting the tri-diagonal system (2.16a) and (2.16b), we obtain all of the unknown amplitudes and the problem is therefore solved. In particular, the reflection and transmission coefficients, *R* and *T*, are given by the first entries of the vectors and , respectively.

### (e) Dividing the transect

For a transect composed of *s* sections, the direct solution of the tri-diagonal system (2.16a) and (2.16b) would require the inversion of a 2(*N*+3)(*s*+1)-square matrix. However, we may restrict the size of the matrices that must be inverted, and hence the computational cost, by dividing the transect into a number of smaller partial transects and reconstructing the solution for the full transect from the solutions of the partial transects.

Therefore, let us suppose the transect *S* to be divided into a finite number, *p* say, of partial transects *P*_{j} ( *j*=1,…,*p*), where *P*_{j} consists of the scattering problem posed by the *p*_{j} discontinuities located at the points {*x*_{cj−1},…,*x*_{cj−1}} with 0=*c*_{0}<*c*_{1}<⋯<*c*_{p}=*s*+1 being fixed but arbitrary integers. For each partial transect *P*_{j}, we may associate the amplitudes and to be the waves incident upon it, and likewise and to be the waves reflected by it. Note that, as the incident waves are generated between sections, in general the amplitudes of the incident evanescent waves here are non-zero. It is possible to calculate the amplitudes within the partial transect *P*_{j} and those scattered by it as a linear combination of these incident amplitudes by solving the restricted tri-diagonal system
2.19
for *j*=*c*_{j−1},…,*c*_{j}−1 and forming a forcing array from the terms and . We may therefore write the amplitudes and (*i*=*c*_{j−1}+1,…,*c*_{j}) in terms of the known matrices and , with
2.20
which only requires inversion of the 2(*N*+3)*p*_{j}-square matrix that is presented by the above tri-diagonal system for the *j*th partial transect.

By solving for each of the *p* partial transects, we find a representation for all of the amplitudes within the full transect *S* through equation (2.20). In particular, we can find expressions for all of the waves that are incident on the partial transects, which are also outgoing waves for a neighbouring partial transect, that is
2.21
for *j*=0,…,*p*. Recalling that and , it is clear that relations (2.21) pose a further tri-diagonal system of dimension 2(*N*+3)*p*, which must be solved for the unknown amplitudes and ( *j*=1,…,*p*). The remaining amplitudes can then be explicitly calculated through equation (2.20) from knowledge of the amplitudes that force the partial transect to which they belong. This completes the solution process.

With *p*=1 or *p*=*s*+1, the method we have just described of dividing a transect collapses to the direct inversion of the tri-diagonal system for the full transect, which is given in equations (2.16a) and (2.16b). For all intermediate values of *p*, however, the two methods are distinct; there is a balance between the computational saving made by solving smaller systems of equations versus the increased computational cost of having to solve multiple systems and also the cost of reconstructing the full transect through system (2.21). In §3*c*(i), we make a numerical investigation of the value of *p* that optimizes computational savings.

Note that the division of the transect could be nested so, for example, we could divide a 1000 m transect into 10 partial transects each 100 m long that comprise 10 partial–partial transects which are 10 m long. Of course, there is no reason to stop at the three layers of this example. In this work, we have not employed nested divisions. The ultimate constraint on transect length is the number of scatterers: the more the scatterers, the larger the memory required to hold arrays such as the amplitudes of the waves within each section. This required memory is the main limit on the maximum size of the problems we may solve.

### (f) Real ice draught data

We make extensive use of ice draught data collected using upward pointing sonar deployed on submarines. These data, collated and archived by NSIDC (2006), are openly available for cruises from 1975 through to 2000. The data are presented as approximately 3700 files, each containing sea ice draughts spaced at 1 m intervals. However, only draughts are reported in these data, hence values of the ice thickness are unknown, and we desire a method for estimating these from the draughts. Mahoney *et al.* (2007) discuss such a method, noting that values of between 4.5 and 9 are reported in the literature for the ratio of keel depth to sail height, the latter value for an iceberg of constant density and porosity and the former for a pressure ridge having porous ice, snow and air pockets in the sail. For simplicity, we take the ratio to be 5 and assume that it will apply to all draughts, hence a thickness is 1.2 times the corresponding draught.

#### (i) Missing segments and the UK cruises

Problematic in the NSIDC data are missing segments, caused by issues such as equipment malfunction. Percival *et al.* (2008) addressed this deficiency by filling the gaps using a stochastic simulation scheme, but for 1 km averages. The problem with using their approach at the 1 m scale is that draughts no longer have a Gaussian distribution (D. B. Percival 2009, personal communication). To address this inadequacy, Percival developed a method for ‘Gaussianizing’ the distribution of draughts so that autocorrelation can be used. This approach is used here.

The NSIDC data contain files from cruises by UK Royal Navy and by US Navy submarines. Data in both sets were processed in a similar manner, but the description of gaps differs for the UK files. The details of how the gaps are recorded for the UK trips mean that it is not possible to determine where in the record they occurred. Because of this uncertainty those UK files that have gaps are not used here.

Two files from the SCICEX 1997 cruise contain unphysically large ridges (200 and 860 m), which are clearly errors, so these files are not used.

#### (ii) Block averaging

We use a process that we term block averaging in which a contiguous block of sections, all lying within a 2δ range, is replaced by a single long section having draught equal to the mean of the draughts in the block. So for example, a 100 m portion of a transect comprising 1 m sections, each having a draught less than δ from a mean of 2 m, can be replaced by a single, 100 m section with a 2 m draught. A steep-sided ridge is likely to be unaffected, unless the parameter δ is chosen too large. The net result is that regions of nearly uniform ice are replaced by long sections of uniform ice, greatly reducing the total number of sections and speeding up the solution process. Clearly there is an issue of balancing errors (these increase as δ is increased) against run-times (these reduce as δ is increased) and this is discussed in §3*c*(ii).

## 3. Model tests

Having described the details of the model developed for this work, we now present results that demonstrate its accuracy. In this section, we establish that the model converges and make comparisons with results from several independent models. The section concludes with an examination of the choice of certain operational parameters, namely, the number of modes *N*, block averaging parameter δ and the number of divisions *p*.

### (a) Demonstration of convergence

Here, there is just one convergence parameter, *N*. Convergence is achieved if a larger value of *N* produces a solution closer to the exact solution of the governing equations. In figure 2, for a 2 m thick, 1 m long pressure ridge in 1 m thick ice for which Archimedean draught is neglected, the reflection coefficient is plotted against *N*. Results are shown for several depths and for periods of 1 s (figure 2*a*) and 10 s (figure 2*b*).

As the number of modes is increased, the reflection coefficient converges, but the rate of convergence is strongly affected by water depth and by wave period, being lower for deeper water and shorter periods. This behaviour is typical of simulations, but the number of modes required to achieve an error of a predetermined size varies according to the transect.

### (b) Verification of the model

In the following, we compare the results from the current model against those from others in the literature. In all cases, we examine only piecewise constant transects and compare results with those of completely independent models that employ different solution methods. The ice geometries that we use here are a simple rectangular pressure ridge with Achimedean draught, a 1000 m transect (excluding Archimedean draught) of sampled NSIDC (2006) data and a 400 m transect of uniformly distributed random ice thickness in which viscosity is included.

The comparisons are satisfying and differ only in minor detail: differences that may be further reduced through using higher resolutions. That the comparisons are so good provides us with confidence in our model.

#### (i) Submergence and the model of Williams & Squire (2008)

Here, a comparison is made with results from Williams & Squire (2008), who describe a model in which Archimedean draught may be optionally accommodated. For a 1 m thick ice sheet punctuated by a 2 m thick ridge, and for two different ridge widths, results are shown in figure 3, in which the inclusion and non-inclusion of draught have been examined.

Apparent in figure 3 is that the two models agree well. For the 20 m wide ridge the minimum (caused by resonance; see Vaughan *et al.* 2007) gets progressively deeper as the number of modes is increased. The slightly different slopes of the curves at high periods is the remaining point of departure of our results from those Williams & Squire (2008) and does not improve with higher number of modes. For the 100 m wide ridge, results are similar, but there the positions of resonance minima are shifted. We have no explanation for the points of departure between the results but, being minor, we attribute no importance to them.

#### (ii) Midlength transects and the model of Vaughan & Squire (2008*a*)

Figure 4 shows a comparison of results between this work and those of Vaughan & Squire (2008*a*). In particular, a 1000 m long transect was sampled from ice draught data obtained using upward pointing sonar (NSIDC 2006). For both sets of results, the Archimedean draught is neglected. The transect used for this comparison, plotted in the centre of figure 4, exhibits a response typical of these transects. The results are in good, but not perfect, agreement. While the results presented here are close to being converged, we observe that different quantities converge at different rates for different periods; in particular, at low periods, convergence rates (with number of modes) are low. Near very large and abrupt changes in draught, more modes are required. The 15 modes used to obtain results from this model for figure 4, and the 100 modes used by Vaughan & Squire (2008*a*), do not ensure sufficient convergence for all transects at all periods. Notably, however, the current model is approximately 10 times faster than that of Vaughan & Squire (2008*a*).

#### (iii) Viscoelasticity and the model of Vaughan & Squire (2008*b*)

In figure 5, a comparison of representative results from this work and from Vaughan & Squire (2008*b*) is presented. The ice comprises a 400 m transect with thicknesses (between 1 and 3 m) drawn from a uniform random distribution and set in sea ice 1 m thick. Draught is neglected but viscosity is accommodated. The results are virtually indistinguishable and can be further improved by obtaining more highly converged solutions from each of the models.

### (c) Optimal model parameters

Here, we examine the choice of model parameters that provides optimal performance. The number of divisions (see §2*e*), the number of modes and the amount of blocking (see §2*f*) are the parameters of interest.

#### (i) The optimal number of divisions

In figure 6, the existence of an optimal number of divisions (*p*) is demonstrated. For this plot, *s*=10 000 random draughts were generated and the time required for a solution using just one processor core measured. The optimal speed-up is a factor of approximately 2 for *N*=1 and 45 for *N*=10, justifying its use.

The times required to construct the block tri-diagonal matrix and to perform an LU decomposition are examples of bottlenecks in the code. The two times vary in relative importance and in magnitude according to choice of *N*, *s* and *p*, but empirically quantifying the time in terms of these variables is very difficult. Such empirical functions would be useful in order to approximate the curves shown in figure 6 and to estimate the position of the minimum representing the optimal value for *p*. In lieu of a tool for providing accurate estimates, we assume that *p*=*s*/5 divisions is sufficiently close to the minimum.

#### (ii) Empirical rules for choosing simulation parameters

In §2*f*, the process of block averaging was described, a method by which the number of sections in a transect, and hence run-times, may be reduced. This method requires the choice of the parameter δ, so we seek the largest value of δ that can be employed while still having an acceptable error. Likewise, we are interested in knowing the lowest value of *N* that ensures a sufficiently converged solution.

To establish such values, a measure of the error in a solution is required. Recall that |*W*| is a function of *x* that may be found at an arbitrary number of points, so, for a large set of points, the error estimate we use here is
3.1
in which *W*_{c} is a highly converged solution. The overline denotes the arithmetic mean taken with respect to the set of points. This error estimate provides a better indicator of convergence across all periods than errors based on *R* or *T*, which are suited to mutually exclusive period ranges.

Figure 7 shows details of this error for the solution of a randomly selected 1000 m long transect, the draughts of which are shown in the central plot. The solution is found for a wide range of δ from 1 cm to 1 m and for *N* from 1 to 20. For each combination of δ and *N*, the solution is evaluated and the error found using equation (3.1).

Figure 7 contains plots of |*R*| and |*T*| for both the highly converged solution (which has *N*=20 and δ=1 cm) and dots for solutions (having different values for *N* and δ) for which *E*_{w}<0.02. The lowest plot shows the corresponding functions |*W*|. For both |*R*| and |*T*|, there is some spread in the results, but the range that these cover is typically small, exceptions occuring at local minima of both |*R*| and |*T*|. Note that the plots of |*W*| are all similar. The value *E*_{w}<0.02 is chosen, rather arbitrarily, to provide solutions that pass the unimpeachable ‘it looks okay’ test in those examples we examined in detail.

To ensure that a specific solution is sufficiently converged, we could run the model at increasing resolutions and judge the convergence from the changing error. This would ensure all results were quite accurate but it would be time consuming. Instead, we choose simulation parameters that provide accurate results for the majority of a set of simulations. For the 1000 1 km long transects sampled by Vaughan & Squire (2008*a*), but with draught included, we obtain results using a range of simulation parameters. From that data, we construct an empirical ‘rule of thumb’, presented in table 1. We use this rule to choose parameters in the following section.

## 4. Scattering and damping in long transects

With tests showing that our model operates accurately, we can now begin to investigate problems of geophysical importance. We begin with an example of wave elevation and decay rate for a single transect. Results for single transects are unlikely to provide useful conclusions, hence the task of parametrizing decay in terms of transect properties is examined.

### (a) Wave decay in very long transects

To identify the decay of wave elevation we find, using least squares, the best-fit exponential, *A**e*^{−αx}. Plots of |*W*| and *α* are shown in figure 8 for a 50 km transect. For non-zero viscosity, the imaginary part of the primary root at the median thickness is also plotted. This is the decay rate that would be observed for a transect of constant thickness.

Figure 8*a* shows that for a 15 s wave, |*W*| depends strongly on *Γ*. Plots of |*W*| appear as the superposition of an oscillating term and an exponential decay, with local variations introduced by partial resonance. The main effects of increasing viscosity are an increased decay and reduced magnitude of oscillations, as expected.

Figure 8*b* shows that, at very low periods, scattering, rather than viscosity, dominates and so the solutions for different viscosities are similar. At high periods, viscosity dominates and the solutions converge to the decay expected for a transect of constant thickness. The period at which the transition from viscous dominance occurs reduces as viscosity increases, but at those periods (10–25 s) of most significance for ice-coupled waves both viscosity and scattering are significant. Note the abrupt minima that appear, particularly for the zero viscosity case, which are caused by length resonance: some periods are passed with significantly less attenuation.

### (b) An ensemble of long transects

Vaughan & Squire (2008*a*) observed that wave elevation is strongly dependent on the details of each particular transect. They sampled transects from the NSIDC data, which they then binned according to either the median or the coefficient of variation of the draughts in the transect. The coefficient of variation is usually defined with the arithmetic mean as the denominator, but Vaughan & Squire (2008*a*) found that the results were sorted better with the coefficient defined by
4.1
Using sorted transects, they found the expected decay rate as a function of *c*_{v}, but only for periods less than 10 s.

We use a similar approach here, but with the 704 transects that have lengths between 45 and 55 km. The solution is found for each transect for periods at 1 s intervals between 5 and 35 s, using run parameters given by table 1. Finding solutions for such long transects at periods lower than 5 s is time consuming and is only possible on computers having in excess of 4 GB of RAM; hence, we have not included these periods. In addition, we find the solution for zero viscosity and for *Γ*=100 Pa s m^{−1}. Simulations were performed on a range of computers and required approximately 9000 cpu hours in total.

In the following, we distribute the 704 transects into 10 bins, according to their values of *c*_{v} or median draught *d*_{0}. These bins are chosen to have approximately the same number of members. We expect approximately 70 transects in each bin. By first sorting all of the transects according to the desired statistic (*c*_{v} or median) in ascending order, the first 70 transects are placed in the first bin, the second 70 in the second bin, etc. The contents of each bin are transects having similar statistical properties. In the following, we describe results drawn from these binned transects.

#### (i) The effect of length on calculated decay

In figure 9, the effect of included length on calculated decay is shown for a single transect. We have plotted the results for two values of viscosity at 35 s, the latter chosen because it demonstrates broadly representative behaviour. In figure 9*a* the decay rates are shown for included lengths of up to 10 km in order to emphasize the interesting detail. figure 9*b* shows the errors in the decay rates, where the error is given by
4.2
where *α*(*l*) and *α*(50) are decay rates calculated for included lengths of *l* km and 50 km, respectively.

The calculated decay rates show the expected behaviour. For included length that is too short, the decay is highly inaccurate and changes markedly with increased length. Periodic features occur because we are fitting a curve to what is essentially a series of varying length sampled from a sine wave.

The decay rates appear to be converging, but the errors tell a different story. Occasionally, because of the periodicity effect in the decay, marked minima occur. The convergence of the error is slow and a sufficiently converged value is not achieved for zero viscosity.

What figure 9 illustrates is that to calculate the decay rate accurately from the wave envelope requires a record having sufficient length. This observation is significant for geophysicists planning to measure decay.

In figure 10, contour plots of against included length and period are presented, where is the mean relative error () for the 704 transects. figure 10*a* shows that for zero viscosity a longer included length is required at higher periods to calculate the decay accurately. This is not surprising because for high periods the wavelength is long and scattering induces significant attenuation only over long distances.

For the viscosity *Γ*=10^{2} Pa s m^{−1}, decay rates can be found more accurately at high periods than in the zero viscosity case. This is because the presence of viscosity, which is much more significant than scattering at high periods, induces greater attenuation that is calculated more readily. At periods under 10 s, errors at the two different viscosities are very similar.

#### (ii) Examples of decay in long transects

In figure 11, the decay rates for the members of a single bin are shown, along with the median decay for the bin and bar graphs showing how the *α* for the transects are distributed. This plot shows that the effect of viscosity is to increase the decay rate and remove the spread at higher periods. The noise in the examples shown in figure 11*a* is a consequence of length resonance: this is removed by viscosity. For low periods, scattering dominates and the plots are similar irrespective of viscosity. A direct comparison with the decay rates calculated by Vaughan & Squire (2008*a*) is inadvisable because those authors did not assimilate draught, a significant factor for the low periods at which they obtained results.

#### (iii) Sorted transects

Here, we sort transects into bins according to values of the median draught (*d*_{0}) and also to values of *c*_{v}. Ten bins are used in each case and median decay rates for these bins are plotted in figure 12 for a selection of the 10 bins. In each case, we have chosen the 1st, 3rd, 5th, 7th and 10th bins.

The theme is the same for all bin medians and irrespective of viscosity: decay is high at low periods and decreases with period. Also apparent is that there are high- and low-period regimes separated by a transition between 10 and 15 s.

The low-period regime is where scattering dominates. The bin medians for *d*_{0} are rough, and while they are sorted here (the low *d*_{0} bin has the highest decay rate and the high *d*_{0} bin has the lowest decay rate), if the medians for more of the bins are plotted, it is apparent that the sorting is actually quite poor.

In the high-period regime, the median decay rates are very smooth and are sorted very well, albeit with a low spread. In this regime, the wave is behaving as though it was propagating through viscous ice of constant thickness *d*_{0}. At the very highest periods, the waves behave as though they are open-water swell waves.

Conversely, when the transects are binned by *c*_{v}, the bin medians are very nicely sorted in the low-period regime. There, scattering dominates and *c*_{v}, which characterizes the roughness of the transect rather well, determines typical response. In the high-period regime decay rate is independent of scattering (and hence of *c*_{v}), so values for *c*_{v} do not sort the bin medians.

## 5. Discussion and conclusions

A new mathematical model is described in this paper that may be used to investigate the decay of ice-coupled waves by scattering and viscous damping. The model domain is two dimensional and vertical and may be used with sea ice having variable thickness. The main set of results that are presented here were generated using transects sampled from data obtained using upward-pointing sonar during submarine voyages (NSIDC 2006). We consider that these transects represent sea ice having piecewise continuous thickness. Among the interesting features of this model is the inclusion of Archimedean draught and viscous damping. An important property of the model is that it may be used to calculate solutions for 50 km long heterogeneous transects, more than 10 times longer than the longest previously examined. Results from the model are compared with those from other models described in the literature. Solutions for simple cases provide evidence that the model accommodates viscosity and draught accurately.

All current models for simulating such problems include semi-infinite plates of constant thickness placed on either side of the transect, and this model is no different. This enables the imposition of radiation conditions, but implies geometries that do not match reality. In numerical studies of scattering, the reflection and transmission coefficients are used to parametrize the scattering process, but these quantities are firmly wedded to the semi-infinite plates in which they are defined. Hence, they are not suitable for real-world measurements. Here, to keep our results pertinent to the geophysics, we assert that the decay rate is a better parameter for quantifying scattering.

Scattering can induce local resonance, suggesting that regions within the transect can support partial standing waves of varying magnitude. The existence of these means that divining physical meaning from the wave amplitude for any single transect is extremely difficult. Our approach is to sort transects into sets having similarity. Two parameters are used for this binning; the median of the ice draughts (*d*_{0}) and the coefficient of variation (*c*_{v}) taken here to be the ratio of standard deviation of draughts to *d*_{0}. With transects sorted by one of these values, we can average across members of the bin in order to determine the typical response expected from a member of that bin.

The main findings of this paper are:

both scattering and viscous damping induce (approximately) exponential wave decay,

for low periods (typically less than 10 s) scattering is the dominant mechanism,

for high periods (typically higher than 20 s) viscosity dominates,

the transition between these two regimes varies according to viscosity and transect parameters,

the coefficient of variation (

*c*_{v}) for the transect characterizes the low-period regime, with a larger value for*c*_{v}causing greater decay, andthe median draught characterizes decay in the high-period regime, a thicker sheet having a larger decay for the same viscosity.

Recalling that only scattering is accommodated in their formulation, the significance of viscosity at high periods provides a potential explanation for the under-prediction of decay rates observed by Kohout & Meylan (2008) in their marginal ice zone model.

Also novel here are insights into the techniques of calculating decay. Not surprisingly, the minimum length of the sample of the wave elevation solution required to calculate the decay accurately is determined by period. High-period waves require a long sample and low-period waves require a much shorter sample. When viscosity is included, the length required to calculate decay for high-period waves reduces markedly.

The sophistication of the approach described in this paper and its success in replicating ice-coupled wave propagation in natural, heterogeneous sea ice cover is such that the authors believe that there is now a real prospect of assimilating waves into global circulation models along with satellite data products. This is one of the planned outcomes of a proposal that is currently under consideration.

## Acknowledgements

This work is supported by the Marsden Fund Council from Government funding administered by the Royal Society of New Zealand and by the University of Otago. We are grateful to Dr T. Williams for making results from his model available for comparison and to Drs L. Kavalieris and D. Percival for their advice on matters statistical.

## Footnotes

- Received April 8, 2009.
- Accepted May 28, 2009.

- © 2009 The Royal Society