A transport equation for flexural-gravity wave propagation under a sea ice cover of variable thickness

Wave Motion 88 (2019) 153–166

Contents lists available at ScienceDirect

Wave Motion journal homepage: www.elsevier.com/locate/wamot

A transport equation for flexural-gravity wave propagation under a sea ice cover of variable thickness ∗

J.E.M. Mosig , F. Montiel, V.A. Squire Department of Mathematics and Statistics, University of Otago, Dunedin, New Zealand

highlights • Derivation of a transport equation for 1D waves under ice. • Derivation of attenuation coefficient in terms of statistical ice cover properties. • Resulting attenuation coefficient appears to overestimate actual attenuation rates.

article

info

Article history: Received 1 July 2018 Received in revised form 11 March 2019 Accepted 18 March 2019 Available online 21 March 2019 Keywords: Sea ice Ocean waves Transport equation Random thickness

a b s t r a c t In the polar regions, ocean wave propagation is affected by the presence of sea ice. In particular, the waves are observed to attenuate exponentially due to scattering and dissipative effects. Here we derive a transport equation and the associated attenuation coefficient for linear ocean wave packets in one horizontal dimension due to scattering only, assuming that the ocean is covered with ice of spatially varying thickness. This thickness variation is assumed to occur at length scales that are comparable to the wavelength but short compared to the observation scale over which the attenuation is measured. We use a multiple scale expansion and a Wigner transform to arrive at the transport equation. We find that the resulting attenuation coefficient generally overestimates the attenuation that we expect to see in a real ice-covered ocean. We argue that this is likely to be due to the one-dimensional nature of our model. © 2019 Published by Elsevier B.V.

1. Introduction The Earth’s atmosphere and oceans are two dynamical systems that are tightly coupled through the exchange of heat, matter and momentum at their mutual interface [1]. In the polar regions, this coupling is regulated by a sea ice cover that grows or shrinks with time due to changing external factors. One of these factors is the effect of ocean waves, which propagate through the ice-infested seas, break the sea ice up, or cause ice floes to move and potentially collide. Concurrently, the presence of ice floes affects the propagation of waves, which is the focus of the present study. Sea ice can affect ocean wave propagation in two distinct ways. First, waves can scatter from ice floes or inhomogeneities of the ice cover. In this case the wave energy is conserved (although the ice temporarily holds some of the energy arising from, e.g., elastic bending). Second, waves can dissipate energy. This can occur due to heat generation via turbulent flow at the ice edges or rough underside, or energy can be transferred to permanent structural changes in the ice. In the extreme case, this results in floe breakup. Incidentally, wave induced floe breakup can lead to a positive ∗ Corresponding author. E-mail addresses: [email protected] (J.E.M. Mosig), [email protected] (F. Montiel), [email protected] (V.A. Squire). https://doi.org/10.1016/j.wavemoti.2019.03.010 0165-2125/© 2019 Published by Elsevier B.V.

154

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

feedback loop, since waves can travel further when the floes are smaller [2]. In reality both scattering and dissipative processes affect wave propagation simultaneously. The above processes have been modelled time and time again with different approaches. Most of these approaches can be associated with one of three categories. First, by floe models that represent individual ice floes and their interaction with waves. Floe models can be two or three dimensional and involve a single floe [3], or an array of floes [4] (see also [5] for a review). Second, by effective medium models that parametrize the effect of the ice cover with the rheological parameters of a postulated ‘‘effective’’ material. Several such models have been proposed, with a variable number of free tuning parameters [6,7]. The third model category comprises transport equation models, where the effect of ice floes is directly integrated into the transport equation that is typically used to describe the ocean wave spectrum on large scales [8]. An extensive overview of these three modelling approaches is reported in [9]. In the present article, we derive a new transport equation model for waves propagating in a highly concentrated sea ice cover, such as pack ice in the Arctic Basin during winter. In this case, the ice may be thought of as a continuous cover of spatially varying thickness. Squire et al. [10] consider a similar situation by modelling ocean wave propagation over a long transect of a sea ice cover of variable thickness, where the thickness profile was taken from submarine sonar measurements. Using a floe model for wave propagation under an ice cover of variable thickness developed by Bennetts et al. [11] and Vaughan et al. [12], Squire et al. find that the attenuation coefficient depends nearly exponentially on the wave period. While Squire et al. [10] investigate the scattering of phase-resolved plane waves by a thickness profile with large variations, e.g. including pressure ridges and polynyas, our model approximates the transport of ocean wave packets at the water/ice interface when the ice cover’s thickness varies on spatial scales that are similar to the packet’s associated wavelength, while the magnitude of the thickness variations is small. More precisely, the goal of this work is to derive an analytic expression for the rate of attenuation of energy that is stored in a wave packet or wave train in one horizontal dimension. Since the wave modes that make up the packet scatter from the sea ice cover, the packet can lose energy to modes that, e.g., travel in the opposite direction (or a different angle if in two horizontal dimensions), while the total energy of the system is conserved. The attenuation coefficient that we seek is central to a large amount of literature on ocean wave/sea ice interaction, as it is required for wave forecasting models such as WAVEWATCH III. The key equation in the latter model is the transport equation that describes the time-evolution of an ensemble of wave packets in phase space (position × wavenumber), as they interact with the random ice cover. The transport equation is also known as the action balance equation, the radiative transfer equation, the kinetic equation and in some restricted context it may also be referred to as an energy balance equation or a Boltzmann equation. In the most general (one-dimensional) case it takes the form

∂t W + ∂k Ω ∂x W − ∂x Ω ∂k W =

∑

Si ,

(1)

i

where the packet frequency Ω = Ω (x, k) depends on the wavenumber k and position x, and the Si terms represent various source terms that affect the evolution of the wave action density W = W (t , x, k) with time t. This wave action density is closely related to the wave energy density, as we will explain shortly. Eq. (1) is rather generic, as it appears in a variety of contexts ranging from statistical mechanics, via light scattering from dust particles [13,14], to ocean wave propagation [15–17]. For a general class of wave equations, Bretherton [18] and Bretherton and Garret [19] show that the wave action W (t , x, k) is conserved for a wave train that propagates in a medium with properties that vary slowly in space. Specifically,

( ) ∂t W + ∂x ∂k Ω W = 0,

(2)

where ∂k Ω is the group velocity of the wave train. While Bretherton and Garret arrive at this result using variational principles, Andrews and McIntyre [20] derive the same general result from considerations about the energy–momentum tensor of the system. The wave action density W (t , x, k) in Eq. (2) equals the wavenumber-resolved wave energy density E(t , x, k) divided by the intrinsic frequency ω ˜ (x, k), which is the wave’s frequency measured in a frame of reference in which the mean state of the medium is locally in equilibrium and at rest. When there is no net flow in the system, then ω˜ = ω is just the wave frequency. When, instead of a single slowly varying wave train, a linear superposition of slowly varying wave trains is considered, then equation (2) generalizes to the transport equation (1) without source terms (Si = 0 for all i), as is explained at length by Komen et al. [17], and Willebrand [21]. In other words, when the medium is inhomogeneous (on large scales), each single wave train satisfies the conservation equation (2) but the ensemble-average of all possible wave trains satisfies the transport equation (1) with the source terms removed. The transport equation is well known in the theory of non-linear ocean waves, where the Si contain terms for wave– wave and wave–wind interactions. Early works on the transport equation for ocean waves [22–25,15] only consider spatially homogeneous systems, where the transport equation is free of spatial gradients of Ω and thus assumes the simple form

∂t W + ∂k Ω ∂x W =

∑ i

Si .

(3)

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

155

In this homogeneous case, the source terms Si describe, for example, interactions with the atmosphere, whitecapping and non-linear wave–wave interactions. The origin of these terms is outside the scope of the present article. Many proposals have been made for a source term that represents the effect of sea ice on the waves. Multiple scattering floe models [26,27,4,2] as well as effective medium models [6,7] can all be adapted to generate a source term in the transport equation. This is how some of these models are implemented in WAVEWATCH III and other wave/ice simulations. In 1989, Masson and Leblond [28] took Hasselmann’s transport equation (3) with a non-linear wave–wave interaction term [15] and introduced a new source term into its right hand side, which models the effect of a sparse field of sea ice floes. Each of these ice floes interacts with the wave field by scattering and absorbing wave energy density. Masson and Leblond [28] assume that each frequency mode of the wave field interacts with the floes independently from other frequencies and that the scattered wave field is incoherent, due to the assumed random and sparse distribution of the floes. They represent each individual floe by a rigid floating cylinder and then compute the ensemble averaged scattered wave field using a Foldy–Twersky integral equation [28,14]. Finally, they integrate this scattered wave field into the transport equation in terms of a transfer function, which maps the directional wave action density spectrum before a scattering event to the spectrum after the event and they then solve the transport equation using a time stepping method. Masson and Leblond’s work was applied and slightly modified by Perrie and Hu [29,30]. Meylan et al. [31] derive a source term for the scattering by sea ice that is similar to the one by Masson and Leblond [28], but Meylan et al. use a much simpler ad hoc argument. They start from the transport equation without source terms, as stated by Phillips [32], who just quoted Bretherton and Garret [19]. Then, Meylan et al. [31] add the source term Sice = −β (t , x, θ ) W (t , x, θ ) +

2π

∫

˜ σ (t , x, θ, θ˜ ) W (t , x, θ˜ ) dθ,

(4)

0

where β describes the overall re-distribution of W (t , x, θ ) (total cross-section) as well as dissipation associated with the wave–ice interaction, and σ is the differential scattering cross-section of a single ice floe. Expression (4) is written in terms of the direction of wave propagation θ instead of the wave vector k, because the open water dispersion relation can be used to recover the magnitude of the wavenumber from the frequency (see, e.g. Eq. (3) in [28]). Meylan et al. [31] also upgrade the rigid cylinder model used by Masson and Leblond [28] for the individual ice floes to a more realistic thin elastic plate model. Meylan and Masson [8] subsequently showed, once an error has been corrected in [31], that the two theories by Meylan et al. [31] and Masson and Leblond [28] are equivalent. In contrast to [8,31,28], we consider here a continuous ice cover with spatially varying thickness. This thickness variation is assumed to occur at length scales that are comparable to the wavelength but short compared to the observation scale over which the attenuation is measured. These conditions lead us to a transport equation of the form (3) with a scattering term that accounts for the wave–ice interaction. Note, that there is no ∂x Ω term — just as in the Meylan–Masson model, because we assume that there is no large scale variation in thickness. A similar approach was taken by Rupprecht [33,34], who also considers wave propagation under a rough floating elastic beam and uses a multi-scale expansion to derive an analytic expression for the attenuation coefficient (also see [35]). Rupprecht finds that this multi-scale expansion technique strongly over-estimates the average attenuation of individual wave fields, as it estimates the attenuation of the average wave field, which is different. This difference was also studied by Bennetts et al. [36] in the case of rough sea floor topography instead of sea ice thickness. In contrast to the works of Rupprecht and Bennetts et al., here we take the ensemble average over all thickness profiles after constructing the Wigner function (see Section 2), which describes wave packets and is not resolving the phase of individual wave modes. Consequently, we circumvent the problem of differentiating between individual and average wave fields. The starting point of our calculation in Section 2 is the article by Ryzhik et al. [37], who show how to construct a transport equation for general wave equations using a Wigner transform. As an example, Ryzhik et al. derive the transport equation that corresponds to the linear Schrödinger equation, which is well known in quantum mechanics and optics. The Schrödinger equation also appears in the context of ocean wave propagation. It can be derived by means of a multiple scale expansion [38]. In Section 3 we derive the Schrödinger equation for the amplitude envelope of wave packets that travel under a sea ice cover with spatially varying thickness, thereby expressing the terms of the previously derived transport equation in terms of measurable properties of the ice cover. In Section 4, we compute the predicted attenuation coefficient and demonstrate that it is consistent with the behaviour of solutions to the Schrödinger equation. Subsequently, we briefly study the steady state solution in Section 5. Finally, we conclude this article in Section 6 with a discussion of our findings and suggestions for future work. 2. Reworking of Ryzhik et al. [37] The theory developed by Ryzhik et al. [37] applies to general wave equations (acoustic, electromagnetic, elastic, etc.) under the assumptions that (a) typical wavelengths are short compared to macroscopic features of the medium, (b) the correlation lengths of the inhomogeneities are comparable to wavelengths and (c) the fluctuations of the inhomogeneities are weak. Note that the second condition allows for a strong interaction between the wave and the inhomogeneities. Amongst other types of wave equations, Ryzhik et al. [37] study the Schrödinger equation

( ) iw 2 ∂x A + i V (x) A = 0. ∂t + Cg ∂x A − 2

(5)

156

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

This Schrödinger equation arises in multiple contexts, such as quantum mechanics and optics. It describes the time evolution of a wave field A = A(t , x) as it propagates with group speed Cg , disperses (controlled by the factor w ) and interacts with a potential field V (x). Here we have added parameters Cg and w , for more generality (Ryzhik et al. set Cg = 0 and w = 1). Assuming that the propagation distance is long compared to the wavelength and that the propagation time is large compared to the period, Ryzhik et al. [37] introduce large-scale coordinates x2 = ε 2 x and t2 = ε 2 t, where ε ≪ 1 is a dimensionless scaling parameter. Moreover, they introduce the scaled amplitude Aε (t2 , x2 ) = A(t2 /ε 2 , x2 /ε 2 ) and assume that the potential is a random function with zero mean that varies on the small O(ε 0 ) scale. Specifically, V (x) = ε V1 (x2 /ε 2 ). In the scaled coordinates, the Schrödinger equation reads

( ) i ε 2 ∂t2 + Cg ∂x2 Aε (t2 , x2 ) − ε 4 w ∂x22 Aε (t2 , x2 ) + i ε V1 (x2 /ε2 ) Aε (t2 , x2 ) = 0. 2

(6)

Note that our ε 2 assumes the role of Ryzhik et al.’s ε and, in the context of quantum mechanics, ε 2 is typically denoted by h ¯ . (In quantum mechanics h¯ has units of energy times time, however, whereas ε is dimensionless.) The key tool that Ryzhik et al. [37] use to convert a wave equation into a transport equation is the scaled Wigner transform, which takes the wave amplitude Aε (t , x) defined in space–time and converts it into a Wigner function W ε (t , x, k) in phase space, defined by

{

W ε (t2 , x0 , x2 , k) = F Aε (t2 , x2 −

ε2 2

⋆

y) Aε (t2 , x2 +

ε2 2

}

(k),

y)

(7)

y

where A⋆ denotes the complex conjugate of A (we allow A to be complex to include information about phase) and F {·}y(k) symbolizes the Fourier transform, defined by

∫

∞

F {f (x)}x(k) = −∞

dx 2π

ei k x f (x) = fˆ (k).

(8)

The inverse Fourier transform is defined accordingly by

{

F −1 fˆ (k)

∞

∫

}

dk e−i k x fˆ (k) = f (x).

(x) = k

(9)

−∞

Note, that the Wigner function W ε (t2 , x0 , x2 , k) depends on both small x0 = x and large x2 = ε 2 x scales, since its evolution is influenced by the potential term V1 (x2 /ε 2 ) = V1 (x0 ). To find an evolution equation for W we write

{ ε2 ε2 ⋆ ∂t W ε = F ∂t Aε (x2 − y) Aε (x2 + y) 2 2 } ε2 ε2 ε ε⋆ + A (x2 − y) ∂t A (x2 + y) (k) 2

2

(10)

y

and then use the Schrödinger equation (6) and its complex conjugate to express the time derivatives ∂t Aε and ∂t Aε ⋆ . In (10) the t arguments are suppressed for brevity. For reasons that will become apparent in Section 3, we let the potential V1 depend on the wavenumber k and we write the potential as V1 (x0 ) = V1 (k, x0 ) =

∑

Ui (k) ∂xi 0 f (x0 /l),

(11)

i

where we leave the coefficients Ui (k) unspecified for now and assume that f is a random function with zero mean, whose values are of order O(1). (The sum on the right hand side of (11) has a finite number of terms. This is the form that we will recover later when we construct the potential term in Section 3. Although this step is not necessary to arrive at the transport equation, it allows us to simplify some expressions later on. We can also express the potential term V1 (k, x0 ) by its Fourier transform Vˆ (k, p) = F {V1 (k, x0 )}x0 (p) =

∑

Ui (k) l (−i p l)i fˆ (p l).

(12)

i

With this, the evolution equation of W ε becomes (after simplification)

( ) ( ) ∂t2 + Cg ∂x2 W ε (x0 , x2 , k) + w k ∂x2 + ε−2 ∂x0 W ε (x0 , x2 , k) { } = ε −1 Lε W ε (x0 , x2 , k) ,

(13)

where, suppressing the x-dependence,

{

}

{

(

Lε W ε (k) = i F −1 Vˆ (k, p) W ε (k + p/2) − W ε (k − p/2)

)} p

(x0 ).

(14)

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

157

The details of this rather long but straightforward calculation can be found in [37] and [9]. We want to estimate the leading order behaviour of W ε when ε is small. Therefore, we seek the limit of Eq. (13) for ε → 0, by expanding W ε (x0 , x2 , k) = W (0) (x2 , k) + ε W (1) (x0 , x2 , k) + ε 2 W (2) (x0 , x2 , k) + O(ε 3 ).

(15)

(0)

Note that we assume that the lowest order term W does not depend on x0 , since we want to capture the large-scale effect of the small scale thickness perturbations in W (0) . (Neglecting the effect of small scale roughness on the short scale implies that we only consider the regime of weak scattering, which is reasonable for waves interacting in most types of sea ice covers.) More local effects would be described by W (1) , W (2) and so on. Substitution of (15) into (13) and collecting orders of ε yields 0=

1(

) w k ∂x0 W (1) + Lε {W (0) } ε ) ( ) ( + ∂τ W (0) + w k ∂x2 W (0) + ∂x0 W (2) + Lε {W (1) } + O(ε ),

(16)

where we use the co-moving derivative ∂τ = ∂t2 + Cg ∂x2 for brevity. Now we exploit the fact that W (0) does not depend on x0 and take the Fourier transform of the lowest order part of (16) to find an expression for the Fourier transform of W (1) (see [9] for details), i.e.

(

F W (1) (x0 , x2 , k)

{

Vˆ (k, p) W (0) (x2 , k + p/2) − W (0) (x2 , k − p/2)

} x0

(p) =

p w (k) k + i θ

) ,

(17)

where, in the spirit of [37], we have added a regularization parameter θ which will later be set to zero. With this, we can write the next order term as

) ( ∂τ W (0) + w k ∂x2 W (0) + ∂x0 W (2) { { ∑ 2 −1 −1 ˆ = −i l F F f (p l) fˆ (q l) (−i p l)m+n Um (k) m,n

(

Un (k − p/2)

W (0) (k) − W (0) (k − p) i θ + q w (k − p/2) (k − p/2)

− Un (k + p/2)

W (0) (k + p) − W (0) (k) i θ + q w (k + p/2) (k + p/2)

)}

} (x0 ),

(x0 ) q

p

where we have used the expanded form of the Fourier transformed potential term (12). We take the ensemble average, denoted by ⟨·⟩, over all possible realizations of the random function f . Following [37], we assume that ⟨∂x0 W (2) ⟩ = 0, i.e. the ensemble average of W (2) does not vary on the short x0 scale. Since we also constructed W (0) to be independent of x0 and f only varies on the x0 scale, the ensemble average has no effect on W (0) , i.e. ⟨W (0) ⟩ = W (0) . Thus,

∂τ W (0) + w k ∂x2 W (0) { { ∑ 2 −1 = −i l F F −1 ⟨fˆ (p l) fˆ (q l)⟩ (−i p l)m+n Um (k) i ,j

(

Un (k − p/2)

W (0) (k) − W (0) (k − p) i θ + q w (k − p/2) (k − p/2)

− Un (k + p/2)

W (0) (k + p) − W (0) (k) i θ + q w (k + p/2) (k + p/2)

)}

} (x0 ).

(x0 ) q

(18)

p

Finally, we assume that the random function f is spatially homogeneous and isotropic, such that its covariance function R(x) only depends on the absolute distance between any two points, i.e.

⟨f (x/l) f (y/l)⟩ = 2 π l R(|x − y|),

(19)

which implies that

ˆ l) δ (p l + q l), ⟨fˆ (p l) fˆ (q l)⟩ = R(p

(20)

ˆ is known as the power spectrum of f and δ (p) is the Dirac delta distribution. Substituting (20) into (18), where R(p) simplifying and taking the limit θ → 0 while using that θ lim = π δ (x), (21) θ →0 θ 2 + x2

158

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

yields the transport equation for W (0) :

∫ ( ) ∂t2 + Cg (k) ∂x2 W (0) + w (k) k ∂x2 W (0) = −Σ (k) W (0) (k) + σ (k, p) W (0) (p) dp,

(22)

where

σ (k, p) =8 π

∑

ˆ − k) l) (i p l)m+n R((p

m,n

( Um (k) Un

k+p 2

) ) ( ( ) k+p δ Cg (p2 − k2 ) ,

(23)

2

and

Σ (k) =

∫

σ (k, p) dp,

(24)

are the differential and total cross-sections, respectively. Eq. (22) has the form of (3) with Ω (k) = k Cg (k) and the source terms, given by the (23) and (24), resemble (4) in the Meylan and Masson model [8]. 3. Propagation of a wavepacket under an ice cover of variable thickness We now seek to construct the Schrödinger equation (5) in the context of water wave propagation under a sea ice cover of variable thickness. Consider a two-dimensional ocean of infinite depth and infinite extent in its one horizontal dimension. When in equilibrium, the water/ice interface is located at z = 0 and the sea floor is located at z = −H, where we consider the limit H → ∞, for simplicity. The water body is described by linear potential flow theory and the ice is modelled as a thin elastic beam of variable thickness that covers the water body. This setup is described by the aforementioned linearized boundary value problem 0 = ∂x2 Φ + ∂z2 Φ ,

−∞
(25)

0 = ∂z Φ ,

z → −∞

(26)

z=0

(27)

0 = ∂x2

+

(

3

Eh

12

E h3 12

)

∂x2 ∂z Φ

∂x4 ∂z Φ + ρi h ∂t2 ∂z Φ + g ρw ∂z Φ + ρw ∂t2 Φ ,

where g = 9.8 m s−2 is the acceleration due to gravity, E ≈ 108 Pa is the Young’s modulus of sea ice, and ρw = 1025 kg m−3 and ρi = 917 kg m−3 are, respectively, the average densities of sea water and sea ice. We are interested in the large-scale behaviour of wave packets that satisfy the above boundary value problem. To this end, we perform a multi-scale analysis similar to that of Mei et al. [38], Section 2.4.1, to bridge the local scale effects, i.e. thickness variations, to the large scale properties of the wave packet propagation. We introduce a dimensionless scaling parameter ε ≪ 1. Note that ε does not have any particular value; it is just postulated to exist and does not appear in the final equation (22). Intuitively, we may think of the three scales as the local scale where the ice cover appears rough, an intermediate scale and an observation scale. Later, in Eq. (31), we also relate ε to the amplitude of the thickness variations. The velocity potential Φ now depends on the multi-scale variables xn = ε n x and tn = ε n t, with n = 0, 1, 2, i.e.

Φ (t , x, z) = Φ (t0 , t1 , t2 , x0 , x1 , x2 , z).

(28)

The use of only three scales is an approximation and higher order scales could be introduced as well, but this is not necessary to derive the transport equation in Section 2. We consider plane waves that are locally time-harmonic with angular frequency ω and wavenumber k, but on larger spatial and temporal scales their amplitude envelope may vary, that is,

( ) Φ (t0 , t1 , t2 , x0 , x1 , x2 , z) = ψ0 + ε ψ1 + ε 2 ψ2 ei (k x0 −ω t0 ) ,

(29)

where

ψn = ψn (t1 , t2 , x1 , x2 , z),

n = 0, 1, 2.

(30)

Note that the ψn do not depend on x0 and t0 . Moreover, we allow the ice thickness to vary only on the short spatial scale x0 , i.e. h(x0 ) = h0 (1 + ε f (x0 /l)),

(31)

where h0 is the mean thickness, f is a random function with expectation value E[f ] = 0 and expected amplitude √ E[f 2 ] = 1, l is the correlation length scale and h0 ε assumes the role of the root-mean-square amplitude of h(x0 ).

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

159

Thus, ε describes the assumed separation scale of the solutions as well as the amplitude of the thickness variations. The notation is not a coincidence; we made these assumptions such that f and ε will turn out to be the same as in Section 2. Substituting the multi-scale expansion (29) into the Laplace equation (25) and collecting orders in ε gives O(ε 0 ) : Lψ0 = 0,

(32)

O(ε ) : Lψ1 = K1 ψ0 ,

(33)

1

O(ε ) : Lψ2 = K2 ψ0 + K1 ψ1 , 2

(34)

where L = ∂ − k , K1 = −2 i k ∂x1 and K2 = −2 i k ∂x2 + ∂ condition (27) are 2 z

2

2 x1 .

Similarly, the three orders of the surface boundary

O(ε 0 ) : Bψ0 = 0

(z = 0),

(35)

O(ε 1 ) : Bψ1 = (C1 + C˜1 )ψ0

(z = 0),

(36)

O(ε ) : Bψ2 = (C2 + C˜2 )ψ0 + (C1 + C˜1 )ψ1

(z = 0),

(37)

2

where the expressions for the operators B, C1 , C2 , C˜1 and C˜2 are given in the Appendix, noting that in the limiting case of an ice cover with constant thickness, the operators C˜1 and C˜2 vanish. Finally, the sea floor boundary condition is limz →−∞ ∂z ψi (z) = 0 for all i = 0, 1, 2. We can solve the homogeneous zeroth-order system by separating variables. We find, that g (38) ψ0 (t1 , t2 ; x1 , x2 ; z) = A(t1 , t2 ; x1 , x2 ) ek z ,

ω

satisfies both the zeroth order Laplace equation and the sea floor condition at z → −∞. The normalizing factor g /ω is chosen to give A the unit of metres. Substituting this general solution into the surface condition gives the local dispersion relation E 3 5 h k + g k ρw − (h0 k ρi + ρw ) ω2 = 0. (39) 12 0 The general solution to the O(ε 1 ) Laplace equation and the sea floor condition is i g ek z (2 k z − 1)

(40) ∂x1 A(t1 , t2 ; x1 , x2 ), 2kω where we have omitted the homogeneous part as it can be absorbed in the first order solution ψ0 . According to the Fredholm alternative theorem this solution exists if and only if

ψ1 (t1 , t2 ; x1 , x2 ; z) = −

∫

0

[

ψ0 K1 ψ0 dz = ψ0 ∂z ψ1 − ψ1 ∂z ψ0 −∞

]0

,

(41)

−∞

since L and B are both self-adjoint. Substituting in the general solutions and evaluating the integral shows that the solvability condition has the form

∂t1 A + Cg ∂x1 A + i V1 (x0 ) A = 0,

(42)

where Cg = dω/dk is the group velocity and V1 (x0 ) = U0 f (x0 /l) + U2 f ′′ (x0 /l),

(43)

is the scattering potential term that arises from the ice thickness variation. The coefficients U0 and U2 are U0 = U2 =

3 E h30 k5 /12 − h0 k ρi ω2 2 ω (h0 k ρi + ρw ) 3 E h30 k3 /12 2 ω l2 (h0 k ρi + ρw )

,

(44)

.

(45)

The procedure we have used to derive the solvability condition at O(ε 1 ) can be repeated to obtain a second-order solvability condition. The general solution to Eq. (34) can be written as

ψ2 (t1 , t2 ; x1 , x2 ; z) = −

i g ek z (2 k z − 1)

∂x2 A(t1 , t2 ; x1 , x2 ) 2kω g ek z (1 + 2 k z (k z − 1)) 2 − ∂x1 A(t1 , t2 ; x1 , x2 ). 4 k2 ω Similar to the first-order case, the solvability condition is

(46)

[ ]0 ψ0 (K2 ψ0 + K1 ψ1 ) dz = ψ0 ∂z ψ2 − ψ2 ∂z ψ0 .

(47)

∫

0

−∞

−∞

160

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

We evaluate the integral and simplify the resulting expression using the first order condition (42) and the dispersion relation (39), resulting in

∂t2 A + Cg ∂x2 A + ζ (x0 ) ∂x1 A −

i 2

w ∂x21 A + i V2 (x0 ) A = 0,

(48)

where we call w = d2 ω/dk2 the dispersion coefficient and ζ (x0 ) and V2 (x0 ) are complicated expressions that depend on x0 , which are given in the Appendix. We now seek to identify the two solvability conditions as parts of the Schrödinger equation

∂t A + c1 ∂x A −

i 2

c2 ∂x2 A + i V (x0 ) A = 0,

(49)

where V (x0 ) = ε V1 (x0 ) + ε 2 V2 (x0 ). Writing (49) in terms of the multi-scale variables t1 , t2 , x1 and x2 gives

(

0 = ε ∂t1 A + c1 ∂x1 A + i V1 (x0 ) A

)

( ) i + ε 2 ∂t2 A + c1 ∂x2 A − c2 ∂x21 A + i V2 (x0 ) A + O(ε 3 ).

(50) 2 Identifying c1 = Cg and c2 = w , we recognize the first- and second-order conditions (42) and (48), respectively, to be the first- and second-order parts of the Schrödinger equation (49), except for the term ζ (x0 ) in the second order condition. This means, that the amplitude envelope A does not exactly evolve according to a Schrödinger equation, but according to the augmented Schrödinger equation

( ) i ∂t A + Cg + ε ζ (x) ∂x A − w ∂x2 A + i ε V1 (x) A = 0,

(51) 2 with an x-dependent modification of the group velocity (Cg depends only on k, not on x). Note that we neglected the V2 term, as it is smaller than the V1 term by a factor ε . The same argument can be made for the ε ζ (x0 ) term, although this term changes the structure of the equation and may allow for solutions that do not exist in the standard Schrödinger equation. We demonstrate in Section 4 that this is not the case, however, and we can indeed neglect the ε ζ (x0 ) term. With the V2 and ζ terms neglected, we arrive precisely at our starting point (5) for the derivation of the transport equation in Section 2. (Note, that all position-dependent terms in Eq. (51) depend on x = x0 only.) Thus, substituting the expressions (44) and (45) for Ui (k) into (22) of Section 2, we obtain a transport equation for ocean waves that travel under a sea ice cover of spatially variable thickness. 4. Results 4.1. Attenuation coefficient We can now directly compute the total cross-section Σ (k) from the statistical properties of the ice. Whereas the differential cross section σ (k, p) symbolizes the rate at which wave energy associated with the wavenumber k is converted into wave energy associated with the wavenumber p, the total cross section indicates how much wave energy is redirected overall. In our one-dimensional setting, there are only two directions: forwards and backwards. The delta distribution in the differential cross section (23) picks out the two contributions p = k (forward to forward) and p = −k (forward to backward) to the integral in (22), but it also allows for more transitions given by the root(s) of w (k). The case p = −k does not contribute to the integral, since Un (0) = 0 for all n, but w (k) typically contributes with one root, which we can find numerically. Assuming a Gaussian power spectrum for the thickness variation, we can evaluate the total cross-section Σ (k) by solving the p-integral in Eq. (24). Since σ (k, p) contains a Dirac delta distribution, the integral turns into a sum over the roots of its argument. Specifically, we use the fact that for any smooth functions f and g,

∫

∞

f (q) δ g(q) dq =

(

−∞

)

∑ f (qi ) , |g ′ (qi )|

(52)

i

where the qi are the roots of g(q) and g ′ (q) is the first derivative of g evaluated at q. Here we define the attenuation rate as α = Σ (k)/Cg (k). A different notion of attenuation rate is discussed later, in Section 5. Fig. 1 shows the contours of the log10 of the attenuation rate α as a function of the wavelength λ = 2 π/k and the correlation length l for two different values for the Young’s modulus, E = 108 Pa and E = 104 Pa. (Sea ice is usually taken to have a Young’s modulus of about 109 Pa in the thin elastic plate approximation, but cracks in the ice effectively reduce the effective modulus.) In both cases we observe a transition from high to low attenuation rates as the wavelength increases. Moreover, as l increases, this transition happens at an increasingly slow rate, noting that for large correlation lengths, a change in l has little effect (the contours become parallel and vertical near the top of the figures). We also observe a sharp maximum at the wavelength where the group velocity vanishes, independent of l (dashed vertical lines). For E = 108 Pa this is at λ ≈ 76 m, while for E = 104 Pa it is at λ ≈ 7 m (not visible in the Figure). Less obvious is

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

161

Fig. 1. Logarithm of the predicted attenuation rate for a range of different wavelengths and correlation lengths. The left panel shows the result for E = 108 Pa and the right panel for E = 104 Pa. The vertical dashed lines indicate the wavelength λ0 at which Cg (2 π/λ0 ) is minimal.

Fig. 2. Spatial attenuation rate α as a function of wavelength λ, when λ and the correlation length l are equal. Red and blue curves show α for the Young’s modulus E = 108 Pa and E = 104 Pa, respectively. The maximum occurs where Cg (k) = 0, and the position of both extrema primarily depends on E. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the minimum that appears between λ = 100 m and 400 m of the left panel figure (where E = 108 Pa). This minimum always occurs at higher wavelengths than the maximum and it is also visible near the left edge of the right panel (where E = 104 Pa). Fig. 2 shows the attenuation rate α as a function of wavelength, when the correlation length is made equal to the wavelength. In essence, Fig. 2 shows α along the diagonal from the bottom left to the top right of each of the two panels of Fig. 1. We chose this case because we derived the transport equation for λ ≈ l. Both the maximum and minimum in the contour plot are also observed in Fig. 2. The Young’s modulus E affects the position of these extrema, while the asymptotic behaviour for large λ and l appears to be independent of E. Specifically, for large λ and l, the attenuation rate α follows a power law that is approximately α ∼ λ−3 . In Fig. 3 we plot the attenuation rate against the wave period T = 2 π/ω(k), where ω(k) is given by the local dispersion relation (39). We fix l = 45 m, which is the correlation length for the ice profile dataset [39] that was obtained between 1975 and 2000 by US and Royal Navy submarines in the Arctic Basin using upward pointing sonar. The same dataset was used by Squire et al. [10] for the numerical computation of the attenuation coefficient using the one-dimensional floe model by Bennetts et al. [11] and Vaughan et al. [12]. Their result is given in Fig. 3 as grey dashed line. (Their model also allows for a viscosity parameter, but the results in Fig. 3 were obtained with zero viscosity.) For the yellow, red and orange curves we chose the mean thickness to be h0 = 0.5 m, 1.0 m and 1.5 m, respectively, while E = 5 × 107 Pa. Coincidentally, these curves resemble the result by Squire et al. [10]. However, if we choose the same parameters as Squire et al. (h0 = 3.25 m, E = 5 × 109 Pa) then we obtain the black curve, which predicts attenuation that is orders of magnitude larger than that of Squire et al. Furthermore, for periods above 24 s our α becomes negative. This may be due to the fact that for large periods the wavelength is long compared to the correlation length l and our model was not constructed for this case. It also hints at why our model does not compare well with the Squire et al. simulation at

162

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

Fig. 3. Attenuation as a function of wave period. Yellow, red and orange curves show the attenuation coefficient when the mean thickness is h0 = 0.5 m, 1.0 m and 1.5 m, respectively, while E = 5 × 107 Pa. The black curve shows the attenuation coefficient for the parameters of the sea ice profile data that were used by Squire et al. [10] (h0 = 3.25 m, E = 5 × 109 Pa). The result of Squire et al.’s numerical computation using a phase-resolved floe model is given as grey dashed lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

lower periods; at these parameters, even for T = 10 s, the wavelength λ = 241 m is more than five times larger than the correlation length l, which again contradicts the assumptions of the model. Furthermore, the thickness profile that was used in Squire et al.’s study does not satisfy our requirement of small localized thickness perturbations. Instead, their ice terrain contains pressure ridges and leads. We also note that our model considers wave packets, while the model used by Squire et al. considers plane waves. We note that the over-estimated attenuation rate is not caused by phase cancellation, as in Rupprecht [33] and Bennetts et al. [36], since we applied the ensemble average to the Wigner function which describes wave packets and does not resolve the phase of individual wave modes. For large periods the asymptotic behaviour of the yellow, red and orange curves in Fig. 3 is a power law (approximately α ∼ T −8 ) in each case. This is different to the result of Squire et al. [10] (grey dashed line), who suggest that the behaviour is exponential. In our model, the mean thickness h0 appears to control the vertical offset of these aforementioned asymptotes, while the Young’s modulus E mostly affects the position of the extrema. We note that various power laws have been constructed as limiting cases of various continuum models by Meylan et al. [40]. Closest to the α ∼ T −8 that we found above is the Keller model for a thin ice cover and small Reynolds number for which α ∝ T −7 . The Keller model describes the ice cover as a viscous layer on the inviscid ocean. This situation can be formulated as a boundary value problem whose solutions allow only certain wave modes to travel and the imaginary part of the dominant wavenumber describes the attenuation (see [40] for details). No sea ice roughness is assumed in the Keller model and the viscous fluid layer is an effective material. Since the Keller model is very different from our approach, it is unlikely that the two power laws are related. The article by Meylan et al. [40] also contains a meta analysis of power laws that have been found experimentally. The best fit power law is given by α ∝ T −3.27 , which is far from the one we found above. 4.2. Consistency with Schrödinger equation The generally very large attenuation rates we found in the previous section surprise us. Therefore, to ensure our results are self-consistent, we implement a method-of-lines solver for the Schrödinger equation (6). At time t = 0, we specify the initial wavepacket Aε (x) = exp(−x2 k2 ),

(53)

where k is the wavenumber associated with this wavepacket. We must now give an explicit representation of the random function f , since we do not consider an ensemble average. Following Bennetts et al. [36], we construct the random function f with

√ f (x) = lim

N →∞

N 2 ∑

N

cos(ai x + bi ),

(54)

i=0

where the ai and bi are random variables, sampled from a normal distribution with mean 0 and variance 1 and a uniform distribution between 0 and 2 π , respectively. For numerical purposes we fix N = 200 to approximate the limit. We also

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

163

Fig. 4. Squared amplitude of a wave packet after t = 20 s of propagation according to the Schrödinger equation (without ζ (x0 ) term). The dashed vertical line is located at x = Cg t. Red: propagation with homogeneous ice cover. Blue, dashed: propagation with rough ice cover. The three panels show how packets with different associated wavelengths propagate (the wavelength λ is indicated). Large packets are least affected, while smaller ones are attenuated or disaggregated entirely. The associated attenuation coefficient is displayed in Fig. 1.

set ε = 0.01 and h0 = 1 m, l = 100 m and E = 108 Pa. To simulate a spatially infinite domain we impose periodic boundary conditions at x = ±400 m. Fig. 4 shows how the effect of the potential term V1 changes with the wavelength that is associated with the wave packet. This is a qualitative consistency check only, since we now consider one individual wave packet with one particular realization of the thickness profile, instead of an ensemble average. The 200 m wave packet (top panel) is barely effected by the ice roughness. The 180 m wave packet (middle panel) attenuates and seems to travel slightly faster than the solution with constant ice cover. Finally, the 140 m wave packet (bottom panel) completely disaggregates after 20 s. This is consistent with the predicted attenuation rates of Cg α = 1.2 × 10−5 s−1 , 4.7 × 10−3 s−1 and 6.5 × 10−1 s−1 , respectively. 4.3. Influence of the ζ (x0 ) term At the end of Section 3 we argued that the x0 dependent term ζ (x0 ), that modifies the group velocity in the Schrödinger equation by Cg → Cg +ε ζ (x0 ), can be neglected. We raised the concern, however, that an x0 -dependent group velocity may lead to solutions that behave differently from solutions to the standard Schrödinger equation (5). To alleviate this concern, we solve both the standard Schrödinger equation and the augmented Schrödinger equation with the same parameters and initial condition and then compare the two solutions. Specifically, we evolve the initial wave packet (53) with a wavelength λ = 2 π /k = 140 m over 80 s according to (6), with and without the ζ (x0 ) term. All other parameters are as before, i.e. E = 108 Pa and ε = 0.01. The result is shown in Fig. 5, which shows three solutions. First, the free solution, i.e. with a homogeneous ice cover (red). Second, the solution with a rough ice cover but no ζ (x0 ) term (blue). And finally, the solution including the ζ (x0 ) term (orange, short-dashed). While the solution with a rough ice cover is strongly attenuated compared to the free solution, the ζ (x0 ) term does not appear to have any significant effect, even in this extreme case of strong attenuation. We have tested this with other parameter configurations as well, drawing the same conclusion. Therefore, we believe that neglecting the ζ (x0 ) term is justified.

164

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

Fig. 5. Squared amplitude of a wave packet after t = 80 s of propagation with λ = 140 m (strong attenuation). Red: propagation without roughness. Blue: propagation with roughness but neglecting ζ (x0 ) term. Orange, short dashed: propagation with roughness and ζ (x0 ) included. Note that at x = ±400 m there is a periodic boundary. There is no significant difference between the augmented and standard solution, even in this extreme case. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

5. The steady state limit We construct the steady state limit of the transport equation by ignoring the time derivative term in Eq. (22). The resulting steady state equation is Cg (k) + w (k) k ∂x W (x, k) = −Σ (k) W (x, k) +

(

)

∫

σ (k, p) W (x, p) dp,

(55)

where we drop the index of x2 and the superscript of W (0) for brevity. To simplify the right hand side, we first write σ (k, p) (see Eq. (23)) as

σ (k, p) = r(k, p) δ (f (k, p)),

(56)

thereby explicitly writing the Dirac delta with argument f (k, p) = Cg 2 (p − k ). Let M be the set of roots of f (k, p). We can show that M always contains the root p = k. Therefore, the steady-state equation (55) simplifies to

( k+p )

∂x W (x, k) = a(k)

∑ pi ∈M′

2

2

r(k, pi (k)) W (x, pi (k)), |∂p f (k, pi (k))|

(57)

where we defined a(k) = (Cg (k) + k ω(k))−1 and we have used Eq. (52). Furthermore, we defined the set M′ = M\{k}. When we separate variables with W (x, k) = X (x)K (k), we find

∂x X (x) X (x)

=

∑ a(k) r(k, pi (k)) K (pi (k)) = µ, |∂p f (k, pi (k))| K (k) ′

(58)

pi ∈M

where µ is some real number that specifies the steady state attenuation rate of X (x) in the general solution X (x) = X0 eµ x ,

(59)

for some X0 . For any given µ, the k-dependent part of Eq. (58) imposes a constraint on the function K (k). Alternatively, for some µ the function K (k) may not exist. Analysing the nature of these complicated constraints is beyond the scope of the present study. It is therefore unclear whether the steady state model produces an exponentially decaying behaviour, as observed in Nature. 6. Discussion and conclusions In this article we have derived a transport equation for ocean wave packets that travel under an ice cover of spatially variable thickness. This thickness variation is assumed to be small in amplitude and the correlation length is assumed to be comparable to the wavelength that is associated with the packet. Furthermore, we only consider one horizontal dimension and linear wave theory for our derivation and we model the ice as a thin elastic plate. This work was inspired by ideas from Ryzhik et al. [37] to derive the transport equation from a Schrödinger equation in Section 2 and ideas presented in Mei et al. [38] that guided our derivation of the Schrödinger equation for a wave packet that propagates under an ice cover of variable thickness in Section 3.

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

165

The predictions of the transport equation are consistent with the behaviour of solutions to the Schrödinger equation, as we have demonstrated in Section 4.2. We also verified that the ζ (x0 ) term we neglected in the derivation has no significant impact on the wave propagation. The predicted attenuation rates (see Figs. 1–3), however, are unrealistically large, especially for short wavelengths. We suspect that this is due to the one-dimensional nature of our model. Future work may therefore focus on a two-dimensional generalization of our derivation. Acknowledgements The authors are appreciative of financial support from the University of Otago (New Zealand), ONR Departmental Research Initiative ‘Sea State and Boundary Layer Physics of the Emerging Arctic Ocean’ (award number N00014-1310279) and ONR-Global (United States) Grant Award N62909-17-1-2080. FM further acknowledges support from the New Zealand Ministry of Business, Innovation and Employment through the Deep South National Science Challenge (contract C01X1445). JEMM and VAS also thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme The Mathematics of Sea Ice Phenomena, supported by EPSRC (United Kingdom) grant number EP/K032208/1, where initial discussions on this paper took place — especially with Dr Luke Bennetts and Associate Professor Michael Meylan. The paper was prepared while JEMM was in receipt of a University of Otago Postgraduate Publishing Bursary. The source code necessary to reproduce the figures presented in this article is available from the authors upon request ([email protected]). Appendix The operators that appear in Eqs. (35)–(37) of Section 3 are B = −ρw ω2 +

(

E 12

)

h30 k4 + g ρw − h0 ρi ω2 ∂z ,

C1 = 2 i ρw ω ∂t1 + 2 i h0 ρi ω ∂t1 ∂z + 4 i C2 = i

E

h30 k3 ∂x2 ∂z +

3

E 2

E 12

(60)

h30 k3 ∂x1 ∂z ,

(61)

h30 k2 ∂x21 ∂z + 2 i h0 ω ρi ∂t2 ∂z

− h0 ρi ∂x21 ∂z + 2 i ω ρw ∂t2 − ρw ∂t21 , (( ) ) E 3 4 E 3 k2 ′′ C˜1 = −3 h0 k + h0 ρi ω2 f (x0 /l) + 3 h0 2 f (x0 /l) ∂z , 12 12 l ) ( ) ( E 3 2 ′ 2 3 3 C˜2 = f (x0 /l) i E h0 k ∂x1 ∂z + 2 i h0 ρi ω ∂t1 ∂z + f (x0 /l) h k ∂z 2 l2 0 ( ) ) ( E E + f (x0 /l)2 − h30 k4 ∂z + f ′′ (x0 /l) −i 2 h30 k ∂x1 ∂z 4 2l ( ) E 3 2 ′′ + f (x0 /l) f (x0 /l) h k ∂z . 2 0 2l

(62) (63)

(64)

Furthermore, the second order potential term V2 (x0 ) is V2 (x0 ) =

− −

E h30 k5 8 ω (h0 k ρi + ρw ) E h30 k3 4 l2

)2 1 ( V1 (x0 ) 2ω

f (x0 /l)2 −

ω (h0 k ρi + ρw ) E h30 k3

4 l2 ω (h0 k ρi + ρw )

)2

f ′ (x0 /l)

(

−

h0 k ρ i h0 k ρi + ρw

V1 (x0 ) f (x0 /l)

f (x0 /l) f ′′ (x0 /l),

(65)

and the term ζ (x0 ) is

(

9

ζ (x0 ) = h0 k l2 ( E h20 k4 − ρi ω (4 Cg k + ω))f (x0 /l) 4

− 2 l2 (2 Cg k (h0 k ρi + ρw ) + ω (h0 k ρi − ρw )) V1 (x0 ) ) ( ) 5 − E h30 k3 f ′′ (x0 /l) / 4 k l2 ω (h0 k ρi + ρw ) . 4

(66)

166

J.E.M. Mosig, F. Montiel and V.A. Squire / Wave Motion 88 (2019) 153–166

References [1] A.E. Gill, Atmosphere-ocean dynamics, in: International Geophysics, no. 30, Elsevier, 1982. [2] F. Montiel, V.A. Squire, Modelling wave-induced sea ice break-up in the marginal ice zone, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 473 (2206) (2017) 20170258, http://dx.doi.org/10.1098/rspa.2017.0258. [3] J. Mosig, F. Montiel, V. Squire, Water wave scattering from a mass loading ice floe of random length using generalised polynomial chaos, Wave Motion 70 (2017) 222–239, http://dx.doi.org/10.1016/j.wavemoti.2016.09.005. [4] F. Montiel, V.A. Squire, L.G. Bennetts, Attenuation and directional spreading of ocean wave spectra in the marginal ice zone, J. Fluid Mech. 790 (2016) 492–522, http://dx.doi.org/10.1017/jfm.2016.21. [5] V.A. Squire, Of ocean waves and sea-ice revisited, Cold Reg. Sci. & Technol. 49 (2) (2007) 110–133, http://dx.doi.org/10.1016/j.coldregions.2007. 04.007. [6] J.E.M. Mosig, F. Montiel, V.A. Squire, Comparison of viscoelastic-type models for ocean wave attenuation in ice-covered seas, J. Geophys. Res. 120 (9) (2015) 6072–6090, http://dx.doi.org/10.1002/2015JC010881. [7] R. Wang, H.H. Shen, Gravity waves propagating into an ice-covered ocean: A viscoelastic model, J. Geophys. Res. 115 (C6) (2010) http: //dx.doi.org/10.1029/2009JC005591. [8] M.H. Meylan, D. Masson, A linear Boltzmann equation to model wave scattering in the marginal ice zone, Ocean Model. 11 (3–4) (2006) 417–427, http://dx.doi.org/10.1016/j.ocemod.2004.12.008. [9] J.E.M. Mosig, Contemporary Wave–Ice Interaction Models (Ph.D. thesis), University of Otago, 2018, URL https://ourarchive.otago.ac.nz/handle/ 10523/7958. [10] V.A. Squire, G.L. Vaughan, L.G. Bennetts, Ocean surface wave evolvement in the Arctic Basin, Geophys. Res. Lett. 36 (22) (2009) L22502, http://dx.doi.org/10.1029/2009GL040676. [11] L.G. Bennetts, N.R.T. Biggs, D. Porter, A multi-mode approximation to wave scattering by ice sheets of varying thickness, J. Fluid Mech. 579 (2007) 413–443, http://dx.doi.org/10.1017/S002211200700537X. [12] G.L. Vaughan, L.G. Bennetts, V.A. Squire, The decay of flexural-gravity waves in long sea ice transects, Proc. R. Soc. A 465 (2009) 2785–2812, http://dx.doi.org/10.1098/rspa.2009.0187. [13] S. Chandrasekhar, Radiative Transfer, Clarendon Press, 1950. [14] A. Ishimaru, Wave Propagation and Scattering in Random Media, John Wiley & Sons, 1999. [15] K. Hasselmann, On the non-linear energy transfer in a gravity-wave spectrum Part 1. General theory, J. Fluid Mech. 12 (4) (1962) 481–500, http://dx.doi.org/10.1017/S0022112062000373. [16] P.A.E.M. Janssen, Long-time behaviour of a random inhomogeneous field of weakly nonlinear surface gravity waves, J. Fluid Mech. 133 (1983) 113–132, http://dx.doi.org/10.1017/S0022112083001810. [17] G.J. Komen, L. Cavaleri, M. Donelan, K. Hasselmann, S. Hasselmann, P.A.E.M. Janssen, Dynamics and Modelling of Ocean Waves, reprint ed., Cambridge University Press, 1996. [18] F.P. Bretherton, Propagation in slowly varying waveguides, Proc. R. Soc. A 302 (1471) (1968) 555–576, http://dx.doi.org/10.1098/rspa.1968.0035. [19] F.P. Bretherton, C.J.R. Garrett, Wavetrains in inhomogeneous moving media, Proc. R. Soc. A 302 (1471) (1968) 529–554, http://dx.doi.org/10. 1098/rspa.1968.0034. [20] D.G. Andrews, M.E. McIntyre, On wave-action and its relatives, J. Fluid Mech. 89 (4) (1978) 647–664, http://dx.doi.org/10.1017/ S0022112078002785. [21] J. Willebrand, Energy transport in a nonlinear and inhomogeneous random gravity wave field, J. Fluid Mech. 70 (1) (1975) 113–126, http://dx.doi.org/10.1017/S0022112075001929. [22] R. Gelci, H. Cazalé, J. Vassal, Utilization des diagrammes de propagation a la prevision energetique de la houle, Bulletin d’information du Comite d’Oceanographie et d’Etude des Cotes 8 (4) (1956) 160–197. [23] R. Gelci, H. Cazalé, J. Vassal, Prévision de la houle. {L}a méthode des densités spectroangulaires, Bulletin d’information du Comité d’Océanographie et d’Etude des Côtes 9 (1957) 416–435. [24] W. Pierson, G. Neumann, R. James, Practical Methods for Observing and Forecasting Ocean Waves by Means of Wave Spectra and Statistics, U. S. Hydrographic Office, 1955. [25] O.M. Phillips, On the generation of waves by turbulent wind, J. Fluid Mech. 2 (5) (1957) 417–445, http://dx.doi.org/10.1017/S0022112057000233. [26] A.L. Kohout, M.H. Meylan, An elastic plate model for wave attenuation and ice floe breaking in the marginal ice zone, J. Geophys. Res. 113 (C9) (2008) http://dx.doi.org/10.1029/2007JC004434. [27] L.G. Bennetts, M.A. Peter, V.A. Squire, M.H. Meylan, A three-dimensional model of wave attenuation in the marginal ice zone, J. Geophys. Res. 115 (C12) (2010) http://dx.doi.org/10.1029/2009JC005982. [28] D. Masson, P.H. Leblond, Spectral evolution of wind-generated surface gravity waves in a dispersed ice field, J. Fluid Mech. 202 (1989) 43–81, http://dx.doi.org/10.1017/S0022112089001096. [29] W. Perrie, Y. Hu, Air–ice–ocean momentum exchange. Part 1: Energy transfer between waves and ice floes, J. Phys. Oceanogr. 26 (9) (1996) 1705–1720, http://dx.doi.org/10.1175/1520-0485(1996)026<1705:AMEPTB>2.0.CO;2. [30] W. Perrie, Y. Hu, Air–ice–ocean momentum exchange. Part 2: Ice drift, J. Phys. Oceanogr. 27 (9) (1997) 1976–1996, http://dx.doi.org/10.1175/ 1520-0485(1997)027<1976:AIOMEP>2.0.CO;2. [31] M.H. Meylan, V.A. Squire, C. Fox, Toward realism in modeling ocean wave behavior in marginal ice zones, J. Geophys. Res. 102 (C10) (1997) 22981–22991, http://dx.doi.org/10.1029/97JC01453. [32] O. Phillips, The Dynamics of the Upper Ocean, second ed., Cambridge University Press, 1977. [33] S. Rupprecht, Wave Propagation on Multiple Scales along Rough Thin-Elastic Solids (Ph.D. thesis), Universität Augsburg, 2018. [34] S. Rupprecht, L.G. Bennetts, M.A. Peter, On the calculation of wave attenuation along rough strings using individual and effective fields, 85 57–66, http://dx.doi.org/10.1016/j.wavemoti.2018.10.007. [35] L.G. Bennetts, M.A. Peter, Approximations of wave propagation in one-dimensional multiple scattering problems with random characteristics, in: Proceedings of the 27th International Workshop on Water Waves and Floating Bodies, Copenhagen, 2012. [36] L.G. Bennetts, M.A. Peter, H. Chung, Absence of localisation in ocean wave interactions with a rough seabed in intermediate water depth, Quart. J. Mech. Appl. Math. 68 (1) (2015) 97–113, http://dx.doi.org/10.1093/qjmam/hbu024. [37] L. Ryzhik, G. Papanicolaou, J.B. Keller, Transport equations for elastic and other waves in random media, Wave Motion 24 (4) (1996) 327–370, http://dx.doi.org/10.1016/S0165-2125(96)00021-2. [38] C.C. Mei, M. Stiassnie, D.K.-P. Yue, Theory and Applications of Ocean Surface Waves, expanded ed., in: Advanced Series on Ocean Engineering, no. 23, World Scientific, 2005. [39] National Snow and Ice Data Center (comp.), Submarine Upward Looking Sonar Ice Draft Profile Data and Statistics, Version 1, G01360, 1998, http://dx.doi.org/10.7265/N54Q7RWK. [40] M.H. Meylan, L.G. Bennetts, J.E.M. Mosig, W.E. Rogers, M.J. Doble, M.A. Peter, Dispersion relations, power laws, and energy loss for waves in the marginal ice zone, J. Geophys. Res., http://dx.doi.org/10.1002/2018JC013776, URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/ 2018JC013776.