2.1. Solving the amplitude profile

Equation (2.5) represents a challenging nonlinear ordinary differential equation. Here, we aim to derive an analytical solution and to facilitate the mathematical treatment we assume the parameters $ C_d$ and $ \alpha$ to be spatially homogeneous (constant). The model can be extended for heterogeneous conditions but numerical solutions should be sought. A polynomial in $ a$ for the integral of the skin drag (2.4) can be found explicitly in closed form for $|v| \geqslant a\varOmega$ , see Appendix A. In effort to enable closed form analytical solutions we seek an approximate polynomial form using an asymptotic expansion of the integrand for $|v |\ll a\varOmega$ that is then extended into the region $|v|\lt a\varOmega$ (see Appendix A). We obtain

(2.7) \begin{equation} I = \int _0^{2\pi } | a\varOmega \sin (\phi ) -v|^3 {\rm d}\phi \begin{cases} \approx \dfrac {8}{3} (a\varOmega )^3 + 12 a\varOmega v^2 , & \text{if } |v |\lt a\varOmega ; \\[8pt] = 3\pi (a\varOmega )^2 |v| + 2\pi |v|^3 , & \text{if } |v| \geqslant a\varOmega . \end{cases} \end{equation}

Note that, in oceanic conditions, $c_g$ is two orders of magnitude lager than the wind- and current-induced drift, i.e. $\mathcal{O}(10)\,\textrm{m}\,\textrm{s}^{-1}$ versus $\mathcal{O}(0.1)\,\textrm{m}\,\textrm{s}^{-1}$ , therefore $|v|\ll c_g$ implying that, with $\alpha \text{ and }\varGamma$ positive, (2.7) is strictly positive, therefore the amplitude $ a$ decreases monotonically in (2.5). For $a_0 \varOmega \gt |v|$ , we define $\hat {x}^*$ to be the location into the ice domain such that $a(\hat {x}^*) \varOmega = |v|$ . Therefore, the condition $|v| \lessgtr a\varOmega$ corresponds to $\hat {x} \lessgtr \hat {x}^*$ in the physical domain. In contrast, for $a_0\varOmega \leqslant |v|$ , one obtains $\hat {x}^*=0$ .

Solving (2.5) with the integral approximation (2.7), the amplitude solutions are written in piecewise form as

(2.8) \begin{equation} a(\hat {x}) \approx a_A(\hat {x}) \, \textbf{1}_{\{\hat {x} \lt \hat {x}^*\}} + a_B(\hat {x}) \, \textbf{1}_{\{\hat {x} \geqslant \hat {x}^*\}}, \end{equation}

where $\textbf{1}$ is the indicator function, equal to $1$ if the condition is true and $0$ otherwise. The piecewise solution components result in

(2.9) \begin{align}& a_{A}(\hat {x}) = \begin{cases} \displaystyle \frac {3\sqrt {\Delta }}{16\varGamma \varOmega ^3} \tan \Biggl (-\frac {\sqrt {\Delta }}{2} \hat {x} + \arctan \Bigl (\frac {16\varGamma \varOmega ^3 a_0 + 3\alpha }{3\sqrt {\Delta }}\Bigr )\Biggr ) - \frac {3\alpha }{16\varGamma \varOmega ^3}, & \delta \gt 1, \\[12pt] \displaystyle \frac {32 \varGamma \varOmega ^3 a_0 - \left (16 \alpha \varGamma \varOmega ^3 a_0 + 3 \alpha ^2\right )\hat {x}} {32 \varGamma \varOmega ^3 + \left (32 \varGamma \varOmega ^3 a_0 + 6 \alpha \right )\hat {x}}, & \delta = 1, \\[12pt] \displaystyle \frac {r_-(a_0-r_+)\exp \left (\dfrac {8 \varGamma \varOmega ^3}{3}(r_- -r_+) \hat {x} \right ) - r_+(a_0-r_-)} {(a_0-r_+)\exp \left (\dfrac {8 \varGamma \varOmega ^3}{3}(r_- -r_+) \hat {x} \right ) - (a_0-r_-)}, & \delta \lt 1, \end{cases} \end{align}

(2.10) \begin{align}& a_B(\hat {x}) = \left [ e^{-(6 \pi \varOmega ^2 |v|\varGamma + 2\alpha )(\hat {x}-\hat {x}^*)}\left (\frac {v^2}{\varOmega ^2}+ \frac {2\pi |v|^3 \varGamma }{3\pi \varOmega ^2 |v| \varGamma + \alpha }\right ) - \frac {2\pi |v|^3 \varGamma }{3\pi \varOmega ^2 |v| \varGamma + \alpha } \right ]^{\tfrac {1}{2}} . \end{align}

Here, $\Delta \equiv \alpha ^2 (\delta ^2-1)$ , $r_\pm \equiv -3(\alpha \pm \sqrt {-\Delta })/(16 \varGamma \varOmega ^3)$ and

(2.11) \begin{equation} \delta = \frac {8\sqrt {2} \varGamma \varOmega ^2 |v|}{\alpha }\,, \end{equation}

which controls the roots of the quadratic polynomial in $ a$ affecting the solution in the $ \hat {x} \lt \hat {x}^*$ region (see Appendix B). Note that the physical interpretation of $ \delta$ will be discussed in § 2.2.

The solution has a few interesting features. The amplitude profile is continuous in $\hat {x}$ (see Appendix B). The solution is nonlinear, i.e. drag effects and other mechanisms cannot be added as linear superposition. For $\hat {x}\lt \hat {x}^*$ , the different branches of the solution in (2.9) differ in form but possess similar qualitative behaviour, where we note that $\hat {x}^*$ depends on $\delta$ . For a wave with initial orbital velocity magnitude $U_0 \lt |v|$ we obtain the solution in region B, $a_B$ , applies to the entire domain $\hat {x}\geqslant 0$ .

A distinctive consequence of such solution is the existence of an extinction location

(2.12) \begin{equation} \hat {x}_{\textit{end}} = \hat {x}^* + \frac {8}{3 \sqrt {2} \pi \alpha \delta +8} \log \left (\frac {5}{2} + \frac {4 \sqrt {2}}{\pi \delta }\right )\,, \end{equation}

such that $a_B(\hat {x}=\hat {x}_{\textit{end}})=0$ . The extinction location is always in region B because it is where amplitude tends to zero and the condition $a\varOmega \lt |v|$ is satisfied. The extinction location is unique to dissipation due to drift-induced skin drag and defines the extent of the wave-affected region. It moves closer to the sea-ice edge for increasing $C_D$ coefficient and $\alpha$ , see solution in red, for $v=0.05$ m s^–1 and $\alpha =7\times 10^{-6}{ \text{ m}^{-1}}$ ( $\delta =3.4$ ), vs in blue, for $v=0.05$ m s⁻¹ and $\alpha =0$ ( $\delta =\infty$ ), in figure 1(a), while its dependence on drift is not trivial. Such extinction locations are not present in the solutions of Kohout et al. (Reference Kohout, Meylan and Plew2011) and Herman et al. (Reference Herman, Cheng and Shen2019), shown in dashed green and orange respectively in figure 1 a), for which $\alpha =0$ and $v=0$ .

In the special case of no drift, i.e. $v=0$ , the solution reduces to

(2.13) \begin{equation} a(\hat {x}) = \frac {3a_0 \alpha e^{-\alpha \hat {x}}}{3\alpha +8a_0\varGamma \varOmega ^3 -8a_0\varGamma \varOmega ^3 e^{-\alpha \hat {x}}}. \end{equation}

Taking $C_d=0$ , i.e. no drag, we recover the exponential attenuation at a rate $\alpha$ , whereas for $\alpha \rightarrow 0$ , the solution converges to the one of Herman et al. (Reference Herman, Cheng and Shen2019). Therefore, (2.13) is an extension Herman et al. (Reference Herman, Cheng and Shen2019) that accounts for additional mechanisms. As a result, our solution predicts lower amplitudes compared with Herman et al. (Reference Herman, Cheng and Shen2019), see figure 1(a).

Figure 1.

Sample amplitude (a) and corresponding attenuation (b) profiles according to our model and formulations by Kohout et al. (Reference Kohout, Meylan and Plew2011) and Herman et al. (Reference Herman, Cheng and Shen2019). Dashed lines denote $x^*$ , which denotes the separation between regions A and B. Amplitude is normalised by the initial amplitude and distance by $x_{\textit{end}}$ for $\alpha =0$ . Attenuation rate is normalised by $\alpha$ . Solutions are shown for $a_0=1$ m, $C_d=0.05$ , $v=0.053\,\mathrm{m\, s^{-1}}$ , $\omega = 0.63\,\mathrm{s^{-1}}$ , $\alpha =7.2\times 10^{-6} \mathrm{\,m^{-1}}$ .

Figure 1(a) displays the normalised amplitude as a function of the distance into the MIZ normalised by the value $ x_{\textit{end}}$ computed for the non-zero drift and exponential attenuation case ( $v=0.053\,\mathrm{m\, s^{-1}}$ and $\alpha =0$ ; the case displayed in blue). The attenuation is stronger when drift is present, as shown in red ( $v\gt 0$ and $\alpha =7\times 10^{-6}$ ), compared with the solution in purple with $v=0$ m s⁻¹ and $\alpha =7\times 10^{-6}$ . The increased attenuation due to other mechanisms also manifests when drift is present, see solutions in blue vs in red.

2.2. Predicting the attenuation rate

To allow for better comparison with attenuation results for non-drifting sea ice as they appear in the literature (Kohout et al. Reference Kohout, Meylan and Plew2011; Herman et al. Reference Herman, Cheng and Shen2019) and to better highlight the attenuation characteristics, we compute the effective amplitude attenuation rate in the moving frame in the ice-covered region defined as

(2.14) \begin{equation} \alpha _{\textit{eff}} \equiv -\frac {1}{a} \frac {{\rm d}a}{{\rm d}\hat {x}}= \alpha \, + \frac {\varGamma }{a^2} I \, . \end{equation}

Using the integral estimation (2.7) we obtain

(2.15) \begin{equation} \alpha _{\textit{eff}}(\hat {x}) \begin{cases} \approx \alpha +\dfrac {8}{3} a \varOmega ^3 \varGamma + \dfrac {12 \varOmega v_{}^2 \varGamma }{a} , & \text{if } \hat {x}\lt \hat {x}^*; \\[10pt] = \alpha +3\pi \varOmega ^2 |v| \varGamma + \dfrac {2\pi |v|^3 \varGamma }{a^2} , & \text{if } \hat {x}\geqslant \hat {x}^*. \end{cases} \end{equation}

Note that, by extending asymptotic solution of $I$ for $|v|\lt a\varOmega$ to $|v|=a\varOmega$ , the decay rate inherits a jump discontinuity at $\hat {x}^*$ , shown in figure 1(b), with the suitability of the integral approximation justified by the discontinuity magnitude being only 7 %. The nonlinearity between skin drag and other exponential decay mechanisms leads to the lower bound

(2.16) \begin{equation} \alpha _{\textit{eff}} \geqslant \alpha +8 \sqrt {2} \varGamma \varOmega ^2 |v| = \alpha (1+\delta ) \,, \end{equation}

with the global minimum located in the region $\hat {x}\lt \hat {x}^*$ if $ a_0\varOmega \geqslant 3\sqrt {2} |v|/2$ , whereas for lower orbital velocity at the sea-ice edge the attenuation monotonically increases. This result allows for a physical interpretation of the dimensionless parameter $ \delta$ : for $\delta \gt 1$ wave attenuation is dominated by drift-induced skin drag and for $\delta \lt 1$ all other mechanisms dominate.

For drifting sea ice, the attenuation rate grows unbounded in the neighbourhood of the extinction location (see solutions in red and blue in figure 1(b)), i.e. $\alpha _{\textit{eff}}(\hat {x}) \sim [2(x_{\textit{end}} - \hat {x})]^{-1}$ , therefore the extinction location affects the magnitude of the attenuation rate in this region. When drift is absent, $\hat {x}^* \rightarrow \infty$ and the solution for $\hat {x}\lt \hat {x}^*$ applies to the entire domain. As a consequence of our model, skin drag dominates attenuation close to the sea-ice edge whereas other processes dominate deeper into the sea ice, where we retrieve the classical exponential attenuation, i.e. $\alpha _{\textit{eff}} \rightarrow \alpha$ for $\hat {x} \rightarrow \infty$ (see solution in purple in figure 1(b); the grey dashed line denotes $\alpha _{\textit{eff}} = \alpha$ ).

3.1. Individual transects

The two transects reported in Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025) are analysed in detail. Transect I, taken on 26 December 2019 at $85^\circ \mathrm{W},\; 69^\circ \mathrm{S}$ , is 61.5 km long and waves have a peak period of 15 s with amplitude $a^{1}_V = 0.45 \text{ m and attenuation},\alpha ^1_{\mathrm{V}} = 1.42 \times 10^{-5} \, \mathrm{m}^{-1}$ . Transect II, taken on 24 May 2019 at $88^\circ \mathrm{W},\; 68^\circ \mathrm{S}$ , is longer ( $x_{\textit{MIZ}}=137.5 \text{ km}$ ) and waves have a shorter peak period ( $T=$ 12 s) with amplitude $a^{1}_V = 0.55 \text{ m and attenuation},\alpha ^1_{{V}} = 2.43 \times 10^{-6} \, \mathrm{m}^{-1}$ .

For transect I the parameters $ v=0.26 \text{ m s}^{-1}, C_d = 6.0 \times 10^{-3}$ and $\alpha = 5.0\times 10^{-6}$ provide good qualitative and quantitative agreement in terms of amplitude decay (see figure 2(a)). Similar performances are achieved by the the non-drifting model, i.e. $v=0$ , and the exponential decay. However, the model with drift outperforms the others, reproducing the observed increase in effective attenuation rate along the transect (between $ x/x_{\textit{MIZ}} \approx 0.2 \text{ and } 0.6$ ; see figure 2(b)). Note that, for the chosen parameters, drift is stronger than the orbital velocity throughout the transect ( $\lvert v \rvert / U_0 \approx 1.2$ ), i.e. only the branch $x\gt \hat {x}^*$ of (2.15) is used. Remarkably, the model predicts an extinction location consistent with the observed wave-affected sea-ice extent, i.e. at $x/x_{\textit{MIZ}} \approx 1$ .

Figure 2.

Amplitude (a,c) and corresponding attenuation profiles (b,d) for transects A (top) and B (bottom). Measurements (red dots) are shown against model predictions. Distance is normalised with respect to the wave-affected sea-ice extent $x_{\textit{MIZ}}$ and amplitude/attenuation with respect to the measurement closest to the sea-ice edge as reported in Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025).

For transect II the parameters $v=0.49 \text{ m s}^{-1}, C_d = 2.8 \times 10^{-3}$ and $\alpha = 1.0\times 10^{-7}$ qualitatively reproduce the measured mild amplitude decay up to $ x/x_{\textit{MIZ}}\lt 0.6$ followed by a sharp decline (see figure 2 c). Other models, i.e. the exponential and non-drifting models, better reproduce the attenuation near to the ice edge, but fail to capture the sharp decline deeper into the sea ice. The predictive performance of the drifting model is clearly highlighted by the attenuation rate, shown in figure 2(d), where the steady attenuation rate near the ice edge followed by rapid growth for $x/x_{\textit{MIZ}} \gt 0.6$ are well captured qualitatively and quantitatively. Note that the observed attenuation rate fluctuates, resulting in negative values close to the edge and again at $x/x_{\textit{MIZ}}\approx 0.4$ , and the models cannot capture this behaviour. A jump discontinuity appears at $x/x_{\textit{MIZ}} \approx 0.1$ , consistent with the choice of $v$ resulting in $\lvert v \rvert / U_0 \approx 0.9$ and the presence of $\hat {x}^*$ within the domain. Unlike our formulation, all other models predict either a decreasing attenuation rate (no drift model) or constant attenuation (exponential model), as shown in figure 2(d).

In summary, transects I and II showcase two different regimes; in transect I drag attenuation is weak ( $\delta =0.6$ ) whereas in transect II it is the dominant mechanism ( $\delta =8.0$ ). Therefore, as expected, in transect I deviations from exponential attenuation are weaker than in transect II; however, also in transect II, the inclusion of exponential attenuation through $\alpha$ provides significantly improved predictions. When drag dominates, the extinction location $x_{\textit{end}}$ occurs earlier and the spike in attenuation is more prominent.

Our model is developed for homogeneous sea-ice conditions, implying constant $C_d$ along the transect. While this hypothesis is more likely to hold in the winter season (transect II) its applicability to transect I is more debatable. Moreover, direct comparison with individual observations is made difficult by inherent measurement uncertainty, including misalignenment of the transect with the wave direction Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025), lack of information on the local sea-ice properties and the possible presence of other forcing mechanisms in the wave transport equation, e.g. wind forcing. Nevertheless, we highlight that our simple model captures the observed feature of wave attenuation in the MIZ.

3.2. Averaged attenuation

Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025) report average attenuation profiles across seasons and sectors of the Antarctic MIZ. To reproduce the observed variability, we drive our model, i.e. (2.15), with $N_0=3\times 10^5$ random realisations with Gaussian distributed drift velocity (with mean $\mu _v$ and standard deviation $\sigma _v$ ), constant throughout the transect. Each realisation corresponds to an individual transect and $N$ denotes the number of transects at any given location according to our model, the decrease is due to individual transects reaching the extinction location. The simulations replicate the variability of drift around Antarctica (Kwok et al. Reference Kwok, Pang and Kacimi2017). The attenuation rate is computed for the wave component $T = 12\ \mathrm{s}$ (angular frequency $\varOmega = 0.52 \,\mathrm{s}^{-1}$ and wavelength $\lambda = 225$ m) and we set the wave amplitude at the sea-ice edge to $a_0=1\,\text{m}$ . The other parameters are $C_d = 0.02$ and $\alpha = 7.0 \times 10^{-6}\ \mathrm{m}^{-1}$ , to match the averaged measured attenuation rate at the sea-ice edge. While the drag coefficient is higher than the one chosen for the individual transects in the previous section, its value is within the range of observed values for the MIZ (McPhee Reference McPhee1979). Moreover, the chosen value allows us to quantitatively match the observed attenuation rate close to the sea-ice edge, as reported in Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025). Smaller values of drag would give rise to a similar qualitative behaviour. For the comparison, only $a\geqslant 0.05 \text{m}$ are averaged, in accordance with the methodology used by Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025) to exclude observations with high uncertainty.

For drift $\mu _v=0.22$ m s⁻¹ ( $\mu _v/U_0=0.42$ ) and standard deviation $\sigma _v=0.03$ m s⁻¹ (figure 3 a) the attenuation only slightly increases up to $x/\lambda \approx 550$ , consistent with the fact that simulations are in weak drag regime ( $\delta =0.42$ ). From $x/\lambda \approx 550$ to $x/\lambda \approx 700$ the increase is more noticeable, up to $\approx 3$ times the attenuation rate at the sea-ice edge, in qualitative agreement with simulations for transect I (cf. figure 2 b), before stabilising farther in the sea-ice domain. The spread of individual profiles, shaded in grey, also grows when the attenuation rate increases. Note that the increase in attenuation corresponds to a reduction in available measurements of amplitude over threshold (shown in red) and, similarly, individual amplitude profiles reaching the extinction location (distribution shaded in blue).

Figure 3.

Averaged attenuation rates (in black; for negative velocity in green) and interquartile range (shaded in grey). On the right axis, number of measurements (in red) and distribution of the extinction location (shaded in blue; not to scale). Insets show the corresponding drift velocity distribution: (a) $\mu _v=0.22$ m s $^{-1}$ $\sigma _v={}0.03$ m s $^{-1}$ ; (b) $\mu _v=0.22$ m s $^{-1}$ $\sigma _v=0.01$ m s $^{-1}$ ; (c) $\mu _v=0.29$ m s $^{-1}$ $\sigma _v=0.03$ m s $^{-1}$ .

With the same average drift ( $\mu _v=0.22$ m s^–1) but lower variability (from $\sigma _v=0.03$ m s^–1 to $\sigma _v=0.01$ m s^–1; figure 3 b) the same qualitative behaviour is observed. The increase in attenuation is sharper, i.e. it starts deeper in the sea ice (at ${\approx} 650$ wavelengths) and reaches higher values ( $\alpha _{\textit{eff}}\approx 5.5 \times 10^{-4}$ vs $\alpha _{\textit{eff}}\approx 3.5 \times 10^{-4}$ ).

The increase in drift ( $\mu _v=0.29$ m s^–1, $\mu _v/U_0=0.55$ and $\sigma _v=0.03$ m s^–1; figure 3 c) results in higher $\delta =0.55$ , still in the weak drag regime, and moves the increase in attenuation rate closer to the sea-ice edge, from $x/\lambda \approx 450$ to $x/\lambda \approx 550$ . A slightly higher attenuation is also achieved compared with the previous cases, up to $\alpha _{\textit{eff}}\approx 6\times 10^{-4}$ , after which it decreases slightly. For smaller drag coefficients results the increase in attenuation occurs further into the ice field. The number of observations decreases sharply due to individual transects reaching their extinction location (cf. figure 1 b), and only less than 0.1 % remain after 800 wavelengths, also explaining the larger variance in effective attenuation.

In figure 3 the mean attenuation rates for negative mean velocity are also reported in the dashed green line. The direction of the drift has only a marginal effect on the effective attenuation profile, i.e. the growth of the attenuation rate is deeper into the sea ice (by ${\approx}10$ wavelengths) but does not affect its magnitude, with a similar effect on the available measurements and extinction location distributions. This is expected for oceanic conditions, as $\alpha \approx \alpha _{\textit{exp}}$ and $\varGamma \approx C_d/( g c_g 2\pi )$ , hence the dependence of direction of the drift is negligible for operational purposes.

Figure 3(c) is remarkably similar to the measured averaged attenuation in Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025), i.e. a linear increase followed by slight decline, also matching the range of observed attenuation rates. Unlike in the measurements, in the model all individual transects are available from the edge but the decrease of available observations in Voermans et al. (Reference Voermans, Fraser, Brouwer, Meylan, Liu and Babanin2025) deeper into the sea ice is also captured. It can be argued that the presence of drift limits the extent of the wave-affected area, and its extent in our model ( $\approx 185$ km; figure 3(c)) is consistent with a MIZ width of the order of 200 km, as reported in observations (Brouwer et al. Reference Brouwer2022) and model simulations (Day et al. Reference Day, Bennetts, O’Farrell, Alberello and Montiel2024).