# Scattering of surface waves by ocean currents: the U2H map

Han Wang1<EMAIL_ADDRESS>Ana B. Villas Bôas2  Jacques Vanneste1
William R. Young3  1School of Mathematics and Maxwell Institute for
Mathematical Sciences, University of Edinburgh, EH9 3FD, UK 2Department of
Geophysics, Colorado School of Mines, Golden CO 80401, USA 3 Scripps
Institution of Oceanography, University of California at San Diego, La Jolla
CA 92093-0213, USA

###### Abstract

Ocean turbulence at meso- and submesocales affects the propagation of surface
waves through refraction and scattering, inducing spatial modulations in
significant wave height (SWH). We develop a theoretical framework that relates
these modulations to the current that induces them.

We exploit the asymptotic smallness of the ratio of typical current speed to
wave group speed to derive a linear map – the U2H map – between surface
current velocity and SWH anomaly. The U2H map is a convolution, non-local in
space, expressible as a product in Fourier space by a factor independent of
the magnitude of the wavenumber vector. Analytic expressions of the U2H map
show how the SWH responds differently to the vortical and divergent parts of
the current, and how the anisotropy of the wave spectrum is key to large
current-induced SWH anomalies.

We implement the U2H map numerically and test its predictions against
WAVEWATCH III numerical simulations for both idealised and realistic current
configurations.

##  1 Introduction

Surface gravity waves (SGWs) propagate through currents resulting from ocean
submesoscale turbulence and from the surface expression of tides and internal
gravity waves. Much work, both historical and recent, has focussed on the
effect of internal waves and tides on SGWs (e.g. Barber, 1949; Perry &
Schimke, 1965; Phillips, 1977; Osborne & Burch, 1980; Tolman, 1990; Hao &
Shen, 2020; Ho et al., 2023). Recognition of the role of submesoscale
turbulence in shaping the open-ocean surface wave field is comparatively
recent and relies on ocean observations and modeling (Romero et al., 2017;
Ardhuin et al., 2017; Romero et al., 2020; Villas Boâs et al., 2020). In this
work we develop a theoretical framework that can be used to understand the
effect of submesoscale turbulence on SGWs.

Refraction and scattering of open-ocean deepwater SGWs by eddies, fronts,
filaments and vortices results in fluctuations in significant wave height
(SWH) with length scales reflecting those of the underlying eddy field – see
figure 1. Fluctuations in SWH modulate SGW breaking and thus affect all
aspects of air-sea exchange (Cavaleri et al., 2012; Villas Bôas et al., 2019).
There are also implications for the formation of extreme waves and remote
sensing.

Scale separation between SGW wavelengths and the larger spatial scale of the
currents makes it possible to adopt a phase-averaged description, focusing on
the density of wave action $\mathcal{A}({\bm{x}},{\bm{k}},t)$ in the
position–wavevector $({\bm{x}},{\bm{k}})$-space. Action density satisfies a
transport equation which, because of its high dimensionality, poses analytic
and numerical challenges, even when linearised by neglecting wave–wave
interactions. As a result, numerous open questions remain about the relation
between the currents and the fluctuations in SWH they induce. Some of these
questions can be addressed by numerical solution of the transport equation
(Villas Boâs et al., 2020). These computations are costly and the results can
be difficult to interpret.

In §2 we develop an alternative approach that directly links SGW amplitude to
current, reducing computational costs and providing new insights. This
approach relies on the smallness of the ratio $\varepsilon=U/c$ of typical
current speed $U$ to SGW group speed $c$. SGWs with wavelengths greater than
10 m have group speed in excess of 2 m s-1. For these relatively long waves
$\varepsilon\ll 1$ holds in all but the most extreme ocean conditions. In a
steady-state scenario and in the absence of current, that is, for
$\varepsilon=0$, the action density $\mathcal{A}({\bm{x}},{\bm{k}})$ can be
taken as spatially uniform, $\bar{\mathcal{A}}({\bm{k}})$ say. For
$0<\varepsilon\ll 1$, the current induces a small, $O(\varepsilon)$
${\bm{x}}$-dependent modulation so that the current-perturbed action spectrum
is
$\mathcal{A}({\bm{x}},{\bm{k}})=\bar{\mathcal{A}}({\bm{k}})+a({\bm{x}},{\bm{k}})$.
The anomaly $a({\bm{x}},{\bm{k}})$ is linearly related to the current.
Similarly the anomaly of any measure of wave amplitude deduced from
$\mathcal{A}({\bm{x}},{\bm{k}})$ is linearly related to the current.

Figure 1: (a) Surface current speeds in a simulation of the California Current
system carried out with the MITgcm model configured as in Villas Boâs et al.
(2020). SWH anomaly $h_{s}$ (given as the relative anomaly $h_{s}/\bar{H}_{s}$
expressed in percent) computed using (b) WW3 and (c) the U2H map. The wave
spectrum is narrow banded in frequency (with periods around 10.3 s and
wavelength 165.5 m) and in direction (propagating primarily along the
$x$-axis, with spectral-width parameter $s=10$; see Appendix B for details).
Panels (a) and (c) can be produced from the notebook accessible at
https://shorturl.at/bef14.

We focus on SWH and obtain an explicit form for the linear map that relates
the SWH anomaly $h_{s}({\bm{x}})$ to the current velocity $\bm{U}({\bm{x}})$.
This map, which we term U2H map, turns out to be a convolution, best expressed
as

$\hat{h}_{s}=\hat{\bm{L}}\bcdot\hat{\bm{U}}$ (1)

in terms of the Fourier transforms $\hat{h}_{s}$ and $\hat{\bm{U}}$ of $h_{s}$
and $\bm{U}$. In section 3 we obtain several alternative and approximate
expressions for the transfer function $\hat{\bm{L}}$ that embodies the U2H
map. Depending on details of the currents $\bm{U}$ and the base-state action
spectrum $\bar{\mathcal{A}}({\bm{k}})$ one of these different expression of
$\hat{\bm{L}}$ may be most effective.

We show a complicated example in figure 1. Panel (a) shows the surface current
speed in a simulation of the California Current system. Panel (b) shows the
SWH anomaly $h_{s}$ computed using WAVEWATCH III (Wave Height, Water Depth,
and Current Hindcasting third generation wave model, hereafter WW3) which
solves the four-dimensional transport equation satisfied by the action
$\mathcal{A}({\bm{x}},{\bm{k}},t)$ (Tolman et al., 2009). Panel (c) shows the
result of applying the U2H map to the current in panel (a). The match between
the results of the (computationally expensive) WW3 computation in (b) and of
the (much cheaper) application of U2H in (c) is excellent. Throughout the
paper we assess the accuracy of the U2H map (1) by comparing its predictions
to numerical simulations using WW3.

In §4 we examine the SWH anomaly $h_{s}$ induced by simple flows such as
vortices. In §5 we consider $h_{s}$ produced by special wave spectra. We show,
in particular, that $h_{s}$ vanishes (to leading order in $\varepsilon$) for
isotropic wave spectra. The complementary limit is highly directional wave
spectra, characteristic of ocean swell. Swell produces strongly anisotropic
SGW anomalies aligned with the dominant direction of wave propagation, i.e.
streaks in SWH. The limit of highly directional wave spectra is delicate in
that it is non-uniform in the small parameter, $\delta$ say, characterising
the spectral width. The results of the present paper assume that
$\delta\gg\varepsilon$. For $\delta=O(\varepsilon)$, the SWH response is
nonlinear in the current and the assumptions leading to the U2H map break
down. An asymptotic form for the SWH in this case is derived in Wang et al.
(2023) under the additional assumption of a localised current.

##  2 Formulation

We start with the conservation equation

$\partial_{t}\mathcal{A}+\bnabla_{\bm{k}}\omega\bcdot\bnabla_{\bm{x}}\mathcal{A}-\bnabla_{\bm{x}}\omega\bcdot\bnabla_{\bm{k}}\mathcal{A}=0$
(2)

for the wave action density $\mathcal{A}({\bm{x}},{\bm{k}},t)$ in
position–wavevector space (e.g. Komen et al., 1996; Janssen, 2004). Here
$\omega$ is the absolute frequency of deep-water SGWs, related to the
intrinsic frequency $\sigma$ by

$\omega({\bm{x}},{\bm{k}})=\sigma(k)+\bm{k}\bcdot\bm{U}({\bm{x}}),\quad\textrm{where}\quad\sigma(k)=\sqrt{gk},$
(3)

$k=|{\bm{k}}|$, and $\bm{U}({\bm{x}})$ is the surface velocity the ocean
current, taken to be horizontal and time independent. The wave energy spectrum
is related to the action density according to

$\mathcal{E}({\bm{x}},{\bm{k}},t)=\sigma(k)\mathcal{A}({\bm{x}},{\bm{k}},t).$
(4)

With equipartition between kinetic and potential energy of deep-water SGW the
root mean square sea-surface displacement is related to the wave energy
spectrum by

$g\eta^{2}_{\mathit{rms}}=\int\mathcal{E}({\bm{x}},{\bm{k}},t)\,\mathrm{d}{\bm{k}}.$
(5)

Our focus is on the spatial distribution of wave energy, obtained by
integrating (4) in ${\bm{k}}$ and conventionally reported in terms of the SWH
defined as $H_{s}({\bm{x}},t)=4\eta_{rms}({\bm{x}},t)$. SWH can be obtained
from the action spectrum with

$H_{s}({\bm{x}},t)=4\left(g^{-1}\int\sigma(k)\mathcal{A}({\bm{x}},{\bm{k}},t)\,\mathrm{d}{\bm{k}}\right)^{1/2}\,.$
(6)

The action equation (2) relies on an assumption of spatial scale separation
between surface gravity waves and currents. We make the further assumption
that the typical current velocity is small compared with the typical group
speed of SGWs:

$\varepsilon\stackrel{{\scriptstyle\text{def}}}{{=}}U/c\ll 1\,.$ (7)

As preparation for a perturbation expansion we highlight that the Doppler
shift is a small correction to the intrinsic frequency by rewriting the wave
frequency (3) as

$\omega({\bm{x}},{\bm{k}})=\sigma(k)+\varepsilon\bm{k}\bcdot\bm{U}({\bm{x}})\,.$
(8)

We develop an approximation to the action equation (2) in dimensional
variables, hence the dimensionless parameter $\varepsilon$ in (8) should be
interpreted as a bookkeeping parameter that is ultimately set to $1$. We
expand the action as

$\mathcal{A}({\bm{x}},{\bm{k}},t)=\bar{\mathcal{A}}({\bm{k}})+\varepsilon
a({\bm{x}},{\bm{k}},t)+O(\varepsilon^{2}).$ (9)

Here we assume that the leading-order action $\bar{\mathcal{A}}({\bm{k}})$ is
independent of space and time. The presence of the currents leads to anomalies
at order $\varepsilon$ captured by $a({\bm{x}},{\bm{k}},t)$. We relate
$a({\bm{x}},{\bm{k}},t)$ to $\bm{U}({\bm{x}})$ by introducing (9) into (2) to
obtain

$\partial_{t}a+\bm{c}\bcdot\bnabla_{{\bm{x}}}a=(\bnabla_{{\bm{k}}}\bar{\mathcal{A}}\bcdot\bnabla_{\bm{x}})\bm{U}\bcdot{\bm{k}},$
(10)

where $\bm{c}(k)=\sqrt{g/4k^{3}}{\bm{k}}$ is the group velocity. We focus on
the steady-state response $a({\bm{x}},{\bm{k}})$, independent of $t$. This
satisfies (10) where the time derivative term is omitted. Causality is
enforced by adding a linear dissipation term to find

$(\bm{c}\bcdot\bnabla_{{\bm{x}}}+\mu)a=(\bnabla_{{\bm{k}}}\bar{\mathcal{A}}\bcdot\bnabla_{\bm{x}})\bm{U}\bcdot{\bm{k}}.$
(11)

The non-dissipative, causal solution is then obtained in the limit $\mu\to
0^{+}$.

We solve (10) in terms of the Fourier transforms

$\hat{a}(\bm{q},{\bm{k}})\stackrel{{\scriptstyle\text{def}}}{{=}}\int
a({\bm{x}},{\bm{k}})\,\mathrm{e}^{-\mathrm{i}\bm{q}\cdot{\bm{x}}}\,\mathrm{d}{\bm{x}}\quad\textrm{and}\quad\hat{\bm{U}}(\bm{q})\stackrel{{\scriptstyle\text{def}}}{{=}}\int\bm{U}({\bm{x}})\,\mathrm{e}^{-\mathrm{i}\bm{q}\cdot{\bm{x}}}\,\mathrm{d}{\bm{x}}.$
(12)

We emphasise the distinction between the newly introduced wavevector $\bm{q}$,
which captures spatial variations at the current scale, and the original
wavevector ${\bm{k}}$, which represents spatial variations of the wave phase.
Introducing (12) into (11) leads to

$\hat{a}(\bm{q},{\bm{k}})=\lim_{\mu\to
0^{+}}\frac{({\bm{k}}\bcdot\hat{\bm{U}})(\bm{q}\bcdot\bnabla_{{\bm{k}}}\bar{\mathcal{A}})}{\bm{c}\bcdot\bm{q}-\mathrm{i}\mu}.$
(13)

The limit $\mu\to 0^{+}$ above is taken in all expressions involving $\mu$ and
we proceed with lighter notation in which the limit is understood.

Our focus is on the SWH, which we expand as

$H_{s}({\bm{x}})=\bar{H}_{s}+\varepsilon h_{s}({\bm{x}})+O(\varepsilon^{2}).$
(14)

The leading order term $\bar{H}_{s}$ is a constant. The anomaly
$h_{s}({\bm{x}})$ is deduced from $a({\bm{x}},{\bm{k}})$ by Taylor expanding
(6):

$h_{s}({\bm{x}})=\frac{8}{g\bar{H}_{s}}\int\sigma(k)a({\bm{x}},{\bm{k}})\,\mathrm{d}{\bm{k}}\,.$
(15)

An analogous formula relates the Fourier transform $\hat{h}_{s}(\bm{q})$ of
$h_{s}({\bm{x}})$ to $\hat{a}(\bm{q},{\bm{k}})$. Substituting (13) into (15)
gives

$\hat{h}_{s}(\bm{q})=\hat{\bm{L}}(\bm{q})\bcdot\hat{\bm{U}}(\bm{q})\,,$ (16)

where

$\hat{\bm{L}}(\bm{q})=\frac{8}{g\bar{H}_{s}}\int\frac{\bm{q}\bcdot\bnabla_{{\bm{k}}}\bar{\mathcal{A}}}{\bm{c}\bcdot\bm{q}-\mathrm{i}\mu}\sigma{\bm{k}}\,\mathrm{d}{\bm{k}}.$
(17)

Eq. (16) shows that $h_{s}$ is obtained from $\bm{U}$ via a linear map – the
U2H map. This map is naturally expressed in terms of Fourier transforms, with
the complex vector $\hat{\bm{L}}(\bm{q})$ acting as a transfer function. It is
clear from (17) that $\hat{\bm{L}}$ depends on the wavevector $\bm{q}$ only
through its orientation: introducing the polar representation

$\bm{q}=q\,\bm{e}_{\bm{q}}\stackrel{{\scriptstyle\text{def}}}{{=}}q\begin{pmatrix}\cos\varphi\\\
\sin\varphi\end{pmatrix}\,,$ (18)

with $\bm{e}_{\bm{q}}$ the unit vector in the direction of $\bm{q}$ and
$-\pi<\varphi\leq\pi$, we can rewrite (17) as

$\hat{\bm{L}}(\bm{q})=\hat{\bm{L}}(\varphi)=\frac{8}{g\bar{H}_{s}}\int\frac{\bm{e}_{\bm{q}}\bcdot\bnabla_{{\bm{k}}}\bar{\mathcal{A}}}{\bm{c}\bcdot\bm{e}_{\bm{q}}-\mathrm{i}\mu}\sigma{\bm{k}}\,\mathrm{d}{\bm{k}}.$
(19)

It follows from the reality of $h_{s}({\bm{x}})$ that
$\hat{h}_{s}(-\bm{q})=\hat{h}_{s}^{*}(\bm{q})$ and hence

$\hat{\bm{L}}(-\varphi)=\hat{\bm{L}}^{*}(\varphi)\,.$ (20)

Eq. (19) implies that, in physical space, $h_{s}({\bm{x}})$ is expressed as a
convolution of $\bm{U}({\bm{x}})$ with a kernel $\bm{L}({\bm{x}})$ – the
inverse Fourier transform of $\hat{\bm{L}}(\bm{q})$ – that is homogeneous of
degree $-2$, that is, $\bm{L}(\lambda{\bm{x}})=\lambda^{-2}\bm{L}({\bm{x}})$.
Linear maps of this type are known as (two-dimensional) Calderón–Zygmund
transforms (e.g. Stein, 1970, Ch. 2). While the right-hand side of (19)
appears ambiguous for $\bm{q}=0$ (since $\varphi$ is then not defined), we
simply take $\hat{h}_{s}(0)=0$, corresponding to the vanishing of the spatial
mean of $h_{s}({\bm{x}})$, consistent with the definition of $h_{s}({\bm{x}})$
as an anomaly.

From (16) and (19) we draw the important conclusion that patterns in $h_{s}$
have scales set by the scales of $\bm{U}$ (not vorticity). But the angular
dependence of $h_{s}$ depends on the wave spectrum. In the next sections we
refine this conclusion by obtaining alternative and approximate expressions
for $\hat{\bm{L}}(\bm{q})$.

##  3 The U2H map

###  3.1 Alternative forms of $\hat{\bm{L}}(\varphi)$

A useful expression for $\hat{\bm{L}}(\varphi)$ is obtained by substituting
the identity
$(\bm{q}\bcdot\bm{p}){\bm{k}}=({\bm{k}}\bcdot\bm{q})\bm{p}-({\bm{k}}^{\perp}\bcdot\bm{p})\bm{q}^{\perp}$
with $\bm{p}=\bnabla_{{\bm{k}}}\bar{\mathcal{A}}$ in (17) to obtain

$\hat{\bm{L}}(\varphi)=\frac{8}{g\bar{H}_{s}}\left(\int\frac{\sigma{\bm{k}}\bcdot\bm{q}}{\bm{c}\bcdot\bm{q}-\mathrm{i}\mu}\bnabla_{\bm{k}}\bar{\mathcal{A}}\,\mathrm{d}{\bm{k}}-\int\frac{\sigma{\bm{k}}^{\perp}\bcdot\bnabla_{\bm{k}}\bar{\mathcal{A}}}{\bm{c}\bcdot\bm{q}-\mathrm{i}\mu}\,\mathrm{d}{\bm{k}}\,\bm{q}^{\perp}\right).$
(21)

Noting that $\bm{c}=\sigma{\bm{k}}/(2k^{2})$, and that the multiplication of
$\mu$ by a positive factor is irrelevant, we rewrite (21) as

$\hat{\bm{L}}(\varphi)=\frac{16}{g\bar{H}_{s}}\left(\int\frac{k^{2}{\bm{k}}\bcdot\bm{q}}{{\bm{k}}\bcdot\bm{q}-\mathrm{i}\mu}\bnabla_{\bm{k}}\bar{\mathcal{A}}\,\mathrm{d}{\bm{k}}-\int\frac{k^{2}{\bm{k}}^{\perp}\bcdot\bnabla_{\bm{k}}\bar{\mathcal{A}}}{{\bm{k}}\bcdot\bm{q}-\mathrm{i}\mu}\,\mathrm{d}{\bm{k}}\,\bm{q}^{\perp}\right).$
(22)

Now, $\mu$ can be safely set to $0$ in the first integral which, on
integrating by parts, reduces to

$\int k^{2}\bnabla_{\bm{k}}\bar{\mathcal{A}}\,\mathrm{d}{\bm{k}}=-2\bm{P},$
(23)

where

$\bm{P}\stackrel{{\scriptstyle\text{def}}}{{=}}\int\\!\bar{\mathcal{A}}({\bm{k}}){\bm{k}}\,\mathrm{d}{\bm{k}}$
(24)

is the wave momentum. For the second integral we use the polar representation

${\bm{k}}=k\begin{pmatrix}\cos\theta\\\ \sin\theta\end{pmatrix}$ (25)

of the SGW wavevector. Noting that
$\bnabla_{\bm{k}}\bar{\mathcal{A}}=\partial_{k}\bar{\mathcal{A}}\,{\bm{k}}/k+\partial_{\theta}\bar{\mathcal{A}}\,{\bm{k}}^{\perp}/k^{2}$
and that ${\bm{k}}\bcdot\bm{q}=kq\cos(\theta-\varphi)$ we obtain

$\hat{\bm{L}}(\varphi)=-\frac{16}{g\bar{H}_{s}}\left(2\bm{P}+\int\frac{k\partial_{\theta}\bar{\mathcal{A}}}{\cos(\theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}{\bm{k}}\,\bm{e}_{\bm{q}}^{\perp}\right)\,,$
(26)

where

$\bm{e}_{\bm{q}}^{\perp}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{pmatrix}-\sin\varphi\\\
\ \ \cos\varphi\end{pmatrix}$ (27)

is the unit vector perpendicular to the wavevector $\bm{q}$.

Starting from (26) we can derive an explicit expression for the transfer
function $\hat{\bm{L}}(\varphi)$ as a Fourier series in $\varphi$. We rewrite
(26) as

$\hat{\bm{L}}(\varphi)=-\frac{16}{g\bar{H}_{s}}\left(2\bm{P}+\,\bm{e}_{\bm{q}}^{\perp}\,\partial_{\varphi}\int_{0}^{2\pi}\frac{\mathcal{P}(\theta)}{\cos(\theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}\theta\right)\,,$
(28)

where

$\mathcal{P}(\theta)\stackrel{{\scriptstyle\text{def}}}{{=}}\int_{0}^{\infty}\\!\\!\\!\bar{\mathcal{A}}(k,\theta)\,k^{2}\mathrm{d}k\,.$
(29)

The wave momentum in (24) is

$\bm{P}=\oint\mathcal{P}(\theta)\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\mathrm{d}\theta\,.$ (30)

With the results above $\hat{\bm{L}}(\varphi)$ depends on the leading-order
action spectrum $\bar{\mathcal{A}}({\bm{k}})$ only through the function
$\mathcal{P}(\theta)$. This function can be expanded in Fourier series as

$\mathcal{P}(\theta)=\sum_{n=-\infty}^{\infty}p_{n}\,\mathrm{e}^{n\mathrm{i}\theta},\quad\textrm{with}\quad
p_{n}=\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!\\!\mathcal{P}(\theta)\,\mathrm{e}^{-n\mathrm{i}\theta}\,\mathrm{d}\theta\,.$
(31)

Computations detailed in Appendix A express the integral in (28) (as $\mu\to
0^{+}$) in terms of the $p_{n}$. This puts the transfer function into the form

$\hat{\bm{L}}(\varphi)=\frac{16}{g\bar{H}_{s}}\left(\bm{e}_{\bm{q}}^{\perp}\,\sum_{n=-\infty}^{\infty}n(-\mathrm{i})^{|n|}\,2\pi
p_{n}\,\mathrm{e}^{n\mathrm{i}\varphi}-2\bm{P}\right)\,,$ (32)

where the wave momentum is

$\bm{P}=\begin{pmatrix}+\mathop{\mathrm{Re}}2\pi p_{1}\\\
-\mathop{\mathrm{Im}}2\pi p_{1}\end{pmatrix}.$ (33)

Equation (32) provides the transfer function $\hat{\bm{L}}(\varphi)$ in a
readily computable form.

###  3.2 Contributions of the divergent and vortical parts of current

The 2D Helmholtz decomposition of $\bm{U}$ into divergent and vortical parts
is

$\bm{U}=\underbrace{\bnabla\phi}_{\bm{U}_{\phi}}+\underbrace{\bnabla^{\perp}\psi}_{\bm{U}_{\psi}}\,,$
(34)

where $\phi$ and $\psi$ are the potential and streamfunction, and
$\bnabla^{\perp}=(-\p_{y},\p_{x})$. The corresponding decomposition of the
Fourier transform $\hat{\bm{U}}$ is

$\hat{\bm{U}}(\bm{q})=\mathrm{i}q\hat{\phi}(\bm{q})\bm{e}_{\bm{q}}+\mathrm{i}q\hat{\psi}(\bm{q})\bm{e}_{\bm{q}}^{\perp}.$
(35)

In view of (16), we can separate the contributions of the divergent and
vortical parts of the currents by expressing the transfer function
$\hat{\bm{L}}(\bm{q})$ in terms of its components along $\bm{e}_{\bm{q}}$ and
$\bm{e}_{\bm{q}}^{\perp}$. Projecting (32) on $\bm{e}_{\bm{q}}$ and
$\bm{e}_{\bm{q}}^{\perp}$ gives

$\hat{\bm{L}}(\varphi)=\hat{{L}}_{\parallel}\,\bm{e}_{\bm{q}}+\hat{{L}}_{\perp}\,\bm{e}_{\bm{q}}^{\perp},$
(36)

where

$\hat{{L}}_{\parallel}(\varphi)=-\frac{32}{g\bar{H}_{s}}\bm{P}\bcdot\bm{e}_{\bm{q}}$
(37)

and

$\hat{{L}}_{\perp}(\varphi)=\frac{16}{g\bar{H}_{s}}\left(\sum_{n=-\infty}^{\infty}n(-\mathrm{i})^{|n|}2\pi
p_{n}\,\mathrm{e}^{n\mathrm{i}\varphi}-2\bm{P}\bcdot\bm{e}_{\bm{q}}^{\perp}\right).$
(38)

Note that, because $\bm{e}_{-\bm{q}}=-\bm{e}_{\bm{q}}$, the symmetry property
(20) implies that
$\hat{{L}}_{\parallel}(-\varphi)=-\hat{{L}}_{\parallel}^{*}(\varphi)$ and
$\hat{{L}}_{\perp}(-\varphi)=-\hat{{L}}_{\perp}^{*}(\varphi)$.

Combining the contributions proportional to $p_{\pm 1}$ (stemming from $n=\pm
1$ in the series and from $2\bm{P}\bcdot\bm{e}_{\bm{q}}^{\perp}$), we can
rewrite (38) as

$\hat{{L}}_{\perp}(\varphi)=\frac{16}{g\bar{H}_{s}}\sum_{n=-\infty}^{\infty}n(-\mathrm{i})^{|n|}2\pi\tilde{p}_{n}\,\mathrm{e}^{\mathrm{i}n\varphi},$
(39)

where

$\tilde{p}_{n}=\begin{cases}2p_{\pm 1}\,,\qquad&\text{if $n=\pm 1$;}\\\
p_{n}\,,\qquad&\text{if $n\not=\pm 1$.}\end{cases}$ (40)

With the form (36) for $\hat{\bm{L}}(\varphi)$ and the Helmholtz decomposition
(35), the linear map (16) becomes

$\hat{h}_{s}(\bm{q})=\mathrm{i}q\hat{{L}}_{\parallel}(\varphi)\hat{\phi}(\bm{q})+\mathrm{i}q\hat{{L}}_{\perp}(\varphi)\hat{\psi}(\bm{q})\,.$
(41)

$\hat{{L}}_{\parallel}$ and $\hat{{L}}_{\perp}$ control the dependence of
$h_{s}$ on, respectively, the divergent and vortical parts of the current. In
general, $\hat{{L}}_{\parallel},\,\hat{{L}}_{\perp}\not=0$, and both the
divergent and vortical parts of the current induce modulations in SWH.
However, we show in §5.3 that for highly directional SGW spectra
$\hat{{L}}_{\perp}\gg\hat{{L}}_{\parallel}$, i.e. the vortical part of the
current is dominant.

##  4 Application to specific currents

Given the wave spectrum $\bar{\mathcal{A}}({\bm{k}})$ and the current
$\bm{U}({\bm{x}})$, the U2H map $\hat{h}_{s}=\hat{\bm{L}}\bcdot\hat{\bm{U}}$,
with $\hat{\bm{L}}$ in (32) or (36)–(38), enables the computation of $h_{s}$.
In this section, we describe this computation, first using a numerical
procedure suitable for arbitrary currents, and then considering particular
cases of purely divergent and purely vortical currents for which analytic
results can be obtained.

Throughout the section we use the separable (in $k$ and $\theta$) action
spectrum detailed in appendix B for illustration. The wavenumber dependence is
defined by a truncated Gaussian in $\sigma(k)$, with peak period $10.3$ s
corresponding to a wavelength of 166 m and a group speed of 8 m s-1. The
angular dependence $D(\theta)$ follows the model of Longuet-Higgins et al.
(1963, LHCS hereafter)

$D(\theta)\propto\cos^{2s}(\theta/2).$ (42)

The parameter $s$ controls the directional spread with large values, say
$s\gtrsim 10$, corresponding to swell-like sea states.

###  4.1 Numerical implementation for arbitrary current

We have implemented the numerical computation of $h_{s}$ using the expression
for $\hat{\bm{L}}(\varphi)$ in (32). The velocity field $\bm{U}({\bm{x}})$ is
discretised on a regular grid and its Fourier transform $\hat{\bm{U}}(\bm{q})$
is obtained on the dual Fourier grid by Fast Fourier transform (FFT). The
transfer function $\hat{\bm{L}}(\varphi)$ is readily computed for the values
of $\varphi$ corresponding to grid-values of $\bm{q}$: for the LHCS spectrum
with integer $s$, which we assume for simplicity, the coefficients $p_{n}$ in
(31) vanish for $n>p$, reducing the series in (32) to a finite sum, and they
take a simple explicit form for $n\leq p$ (see appendix B) The inverse FFT of
the product
$\hat{h}_{s}(\bm{q})=\hat{\bm{L}}(\varphi)\bcdot\hat{\bm{U}}(\bm{q})$ yields
$h_{s}({\bm{x}})$ on the spatial grid.

This is the implementation we use for panel (c) of figure 1. We extracted the
snapshot of surface currents in a $600\times 600$ km domain shown in panel (a)
from a MITgcm simulation. To prevent numerical artefacts due to the non-
periodicity of the currents, we carried out the computation of $h_{s}$ in a
larger domain, zero-padding $\bm{U}$ in 150 km wide strips at the periphery of
the domain. A Jupyter Notebook of this implementation is available on
https://shorturl.at/bef14, where users can customize the input currents and
incoming wave spectrum. The resulting $h_{s}$ is in good agreement with the
WW3 simulation shown in panel (b). The setup of the simulation is as described
in Wang et al. (2023) except that the forcing imposes the spectrum
$\bar{\mathcal{A}}({\bm{k}})$ on the entire boundary i.e. waves enter the
rectangular domain from all four sides. This improved formulation ensures that
the wave spectrum in the absence of currents is uniform even for spectra with
broad directional spread. (In Wang et al. (2023) waves entered the
computational domain only from the western boundary. Even in the complete
absence of currents the resulting steady-state solution decreased with $x$ as
‘wave shadows’ from the northern and southern boundaries encroach into the
centre of the domain.)

We now consider idealised scenarios to gain insight into the dependence of
$h_{s}$ on $\bm{U}$.

###  4.2 Divergent current

For a purely divergent current, with $\bm{U}_{\psi}=0$, (41) reduces to

$\hat{h}_{s}(\bm{q})=-\frac{32}{g\bar{H}_{s}}\mathrm{i}q\hat{\phi}(\bm{q})\,\bm{P}\bcdot\bm{e}_{\bm{q}}\,.$
(43)

Since $\mathrm{i}q\bm{e}_{\bm{q}}=\mathrm{i}\bm{q}$, the inverse Fourier
transform of (43) is

$h_{s}=-\frac{32}{g\bar{H}_{s}}\bm{U}_{\phi}\bcdot\bm{P}\,.$ (44)

Thus the SHW anomaly that arises in response to a divergent current is
proportional to the component of the current’s velocity along the wave
momentum. In particular, the response is local and vanishes where the current
vanishes.

Figure 2: SWH anomaly $h_{s}$ for the divergent flow with Gaussian divergence
(45) with characteristic radius $r_{v}=25$ km (indicated by the dashed circle)
and maximum speed $0.8$ m s-1. The results of WW3 simulations (left) are
compared with the U2H prediction (44) (right) for three values of the
parameter $s$ characterising the directionality of the wave spectrum.

We illustrate (44) with a simple axisymmetric, divergent current whose
divergence is the Gaussian

$\bnabla\bcdot\bm{U}=\nabla^{2}\phi=\frac{\kappa}{2\pi
r_{\mathit{v}}^{2}}\,\mathrm{e}^{-r^{2}/2r_{\mathit{v}}^{2}},$ (45)

where $r_{\mathit{v}}=25\,\text{km}$ is the characteristic radius and $\kappa$
is the area flux, set such that the maximum current speed
$U_{m}=\sqrt{U^{2}+V^{2}}$ is $0.8$ m s-1.

Figure 2 compares the U2H prediction (44) for this current with results for
WW3 simulations for three values of the directionality parameter $s$. Figure 2
confirms the validity of the U2H prediction and the local nature of the SHW
response. This response to divergent currents has a spatial structure
independent of the directional spread of wave energy i.e. $h_{s}$ in (44)
depends only on $\bm{P}$. This striking result is in sharp contrast with the
response to vortical currents as we show next.

###  4.3 Vortical current

For a purely vortical currents, $\bm{U}_{\phi}=0$ in (34) and the U2H map in
(41) is determined by the scalar transfer function
$\hat{{L}}_{\perp}(\varphi)$, which is explicitly computed from the series in
(39). As a demonstration, we consider a Gaussian vortex, with zero divergence
and vorticity in physical and Fourier space given by

$\zeta({\bm{x}})=\frac{\kappa}{2\pi
r_{\mathit{v}}^{2}}\mathrm{e}^{-r^{2}/(2r_{\mathit{v}}^{2})}\quad\textrm{and}\quad\hat{\zeta}(\bm{q})=\kappa\mathrm{e}^{-r_{\mathit{v}}^{2}q^{2}/2},$
(46)

where $\kappa$ is the circulation. We take advantage of the axisymmetry of
this flow to carry out the Fourier inversion leading to $h_{s}({\bm{x}})$
analytically. Calculations detailed in Appendix D yield the explicit
expression

$h_{s}({\bm{x}})=-\frac{16\mathrm{i}}{g\bar{H}_{s}}\frac{\kappa}{r_{\mathit{v}}}\sqrt{\frac{\pi}{2}}\mathrm{e}^{-r^{2}/4r_{v}^{2}}\sum_{n=-\infty}^{\infty}n\,\tilde{p}_{n}\,I_{|n|/2}(r^{2}/4r_{v}^{2})\,\mathrm{e}^{\mathrm{i}n\nu},$
(47)

where ${\bm{x}}=r(\cos\nu,\sin\nu)$ and the $I_{|n|}$ are modified Bessel
functions. The coefficients $\tilde{p}_{n}$ depend only on the wave spectrum
and are related to the already obtained $p_{n}$ according to (40). Formula
(47) has the advantage over the general implementation of the U2H map
described in §4.1 in that it gives $h_{s}({\bm{x}})$ at any location without
the need for entire computational domains in both physical and Fourier
domains.

Figure 3 compares the SWH anomaly (47) with that obtained in WW3 simulations.
The parameters $r_{v}=25$ km and $U_{m}=0.8$ m s-1 are the same as those of
the divergent flow in §4.2. The SWH response in figure 3 is very different
from the response to divergent currents in figure 2. The SWH anomaly in figure
3 extends beyond the vortex. For the swell-like case $s=10$ there is a wake-
like feature decaying slowly in the direction of wave propagation. This
physically important limiting case is discussed in §5.3 and in Wang et al.
(2023).

Figure 3: Same as figure 2 but for the vortical Gaussian vortex with vorticity
(46).

##  5 Particular wave spectra

In this section we examine the role of the wave spectrum
$\bar{\mathcal{A}}({\bm{k}})$ in shaping the SWH anomaly by considering
special and limiting cases.

###  5.1 Isotropic wave spectrum

Figure 4: SWH anomaly $h_{s}$ computed using WW3 for an isotropic wave
spectrum ($s=0$ in the LHCS model) with the Gaussian vortical currrent (46).
Left panel: $h_{s}$ for $U_{m}=0.8$ m s-1, with the dashed circle indicating
the vortex radius. Right panel: cross section of $h_{s}$ at $x=0$ (slicing
through the centre of the vortex) for $U_{m}=1.6$ m s-1 (blue solid line) and
$U_{m}=0.8$ m s-1 (yellow solid line). The yellow dashed line is obtained by
multiplying $h_{s}$ for $U_{m}=0.8$ m s-1 by 4.

Eq. (26) shows that $\hat{\bm{L}}(\varphi)=0$ if the SGW action (or energy)
spectrum is isotropic since $\partial_{\theta}\bar{\mathcal{A}}({\bm{k}})=0$
and $\bm{P}=0$ as a consequence. Thus, if the wave spectrum is isotropic,
currents do not induce modulations of the SWHs (at the order we consider).

To verify this, we ran WW3 simulations with the isotropic wave spectrum
obtained by setting $s=0$ in the LHCS model of Appendix B and the currents of
either the MITgcm simulation of figure 1 or the Gaussian vortex of 3. The SWH
anomaly $h_{s}$ in both cases (MITgcm outcome not shown) is at most 3%, much
smaller than found for anisotropic spectra. The small but non-zero $h_{s}$ for
isotropic spectra can be attributed to quadratic effects not captured by the
linear U2H map. We confirmed this by checking that increasing the velocity of
the Gaussian vortex by an overall factor of 2 (setting $U_{m}=1.6$ m s-1
instead of 0.8 m s-1) increases $h_{s}$ by a factor 4, as shown in Figure 4.

###  5.2 Mildly anisotropic wave spectrum

The U2H map is particularly simple for the spectrum

$\bar{\mathcal{A}}({\bm{k}})=\mathcal{A}_{0}(k)+\mathcal{A}_{c}(k)\cos\theta,$
(48)

e.g., as in the LHCS spectrum with $s=0$ used for the first row of figures 2
and 3. For the action spectrum in (48) the wave momentum (24) can be written
as

$\bm{P}=|\bm{P}|\begin{pmatrix}1\\\
0\end{pmatrix}=|\bm{P}|\cos\varphi\,\bm{e}_{\bm{q}}-|\bm{P}|\sin\varphi\,\bm{e}_{\bm{q}}^{\perp},$
(49)

where

$|\bm{P}|=\pi\int\\!\\!\mathcal{A}_{c}(k)k^{2}\mathrm{d}k.$ (50)

The function $\mathcal{P}(\theta)$ defined in (29) is then

$\mathcal{P}(\theta)=\int\\!\\!\mathcal{A}_{0}(k)k^{2}\mathrm{d}k+\frac{|\bm{P}|}{\pi}\cos\theta.$
(51)

and $2\pi p_{1}=2\pi p_{-1}=|\bm{P}|$. The transfer function in (32) reduces
to

$\hat{\bm{L}}(\varphi)=\frac{32}{g\bar{H}_{s}}|\bm{P}|\left(2\sin\varphi\,\bm{e}_{\bm{q}}^{\perp}-\cos\varphi\,\bm{e}_{\bm{q}}\right)\,.$
(52)

With the Helmholtz decomposition (35), the U2H map is

$\displaystyle\hat{h}_{s}(\bm{q})$
$\displaystyle=\hat{\bm{L}}(\bm{q})\bcdot\hat{\bm{U}}(\bm{q})\,,$ (53)
$\displaystyle=\frac{32}{g\bar{H}_{s}}\left(2|\bm{P}|\,\mathrm{i}q\sin\varphi\,\hat{\psi}-|\bm{P}|\,\mathrm{i}q\cos\varphi\,\hat{\phi}\right)\,.$
(54)

The inverse Fourier transform can be taken by inspection and the result
written as

$h_{s}=-\frac{32}{gH_{s}}\left(2\bm{U}_{\psi}\bcdot\bm{P}+\bm{U}_{\phi}\bcdot\bm{P}\right)\,.$
(55)

In this simple case, only the component of current along $\bm{P}$ produces a
SWH anomaly which turns out to be local, vanishing where the current vanishes.
In (55) the vortical part of the current, $\bm{U}_{\psi}$, is twice as
effective as the divergent part, $\bm{U}_{\phi}$. The divergent contribution
in (55) is identical to that in (44) for arbitrary wave spectra. This example,
which applies to the first row of figures 2 and 3, shows that the response to
divergent currents is not always negligible relative to the vortical response.
This is in contrast with Villas Boâs et al. (2020)’s suggestion that only the
vortical part of the current affects $h_{s}$. The next section however shows
that the imprint of the vortical part of the current is much larger than that
of the divergent part for highly directional spectra.

###  5.3 Highly directional wave spectrum

We conclude above that $h_{s}$ is small for an isotropic wave spectrum. It is
of interest to examine the opposite limit of a highly directional wave
spectrum. This limit corresponds to an action spectrum of the form

$\bar{\mathcal{A}}(k,\theta)=\delta^{-1}\bar{\mathcal{A}}(k,\Theta),\quad\textrm{where}\quad\Theta=\theta/\delta,$
(56)

with $\delta\ll 1$ the relevant small parameter, and we assume that $\theta=0$
is the primary propagation direction. The prefactor $\delta^{-1}$ ensures that
the action spectrum integrated over $\theta$ is $O(1)$. Correspondingly,

$\mathcal{P}(\theta)=\delta^{-1}\mathcal{P}(\Theta).$ (57)

For simplicity, we abuse notation by using the same symbols
$\bar{\mathcal{A}}$ and $\mathcal{P}$ on both sides of (56)–(57),
distinguishing them by their arguments. We treat $\delta$ as a bookkeeping
parameter that can be set to $1$ in the final result.

Taking (28) as a starting point, we obtain an asymptotic approximation to
$\hat{\bm{L}}(\varphi)$ in appendix C. There we show that the dominant
contribution to $\hat{\bm{L}}(\varphi)$ comes from the integral term and is
large in small regions around $\varphi=\pm\pi/2$. In terms of the rescaled
variable

$\Phi_{\pm}=(\varphi\mp\pi/2)/\delta=O(1),$ (58)

the leading-order approximation to the transfer function in these regions is

$\hat{\bm{L}}(\varphi)\sim\hat{{L}}_{\perp}(\varphi)\begin{pmatrix}\mp 1\\\
0\end{pmatrix}\sim\frac{16}{g\bar{H}_{s}\delta^{2}}\,\partial_{\Phi_{\pm}}\int_{-\infty}^{\infty}\frac{\mathcal{P}(\Theta)}{\Theta-\Phi_{\pm}\mp\mathrm{i}\mu}\,\mathrm{d}\Theta\,\begin{pmatrix}1\\\
0\end{pmatrix}.$ (59)

Thus $\hat{\bm{L}}(\varphi)$ is dominant and $O(\delta^{-2})$ in narrow,
$O(\delta)$, sectors around $\varphi=\pm\pi/2$. We conclude that:

  1. (i)

for typical $\hat{\bm{U}}(\bm{q})$, patterns of $h_{s}$ take the form of
structures elongated in the direction of propagation of the waves, with an
$O(\delta)$ aspect ratio i.e. streaks;

  2. (ii)

$\hat{{L}}_{\perp}=O(\delta^{-2})\hat{{L}}_{\parallel}$, since
$\hat{{L}}_{\parallel}=O(1)$ (see (76)). According to (41), this implies that
the divergent-free, vortical part of the velocity field $\bm{U}$ has an
asymptotically larger impact on $h_{s}$ than the potential part;

  3. (iii)

highly directional waves produce SWH anomalies larger by a factor
$\delta^{-1}$ than those induced for spectra with $O(1)$ directional spread.
(This estimate we account for both the factor $\delta^{-2}$ in (76) and the
$O(\delta)$ width of the support of $\hat{{L}}_{\perp}(\varphi)$ implied by
(58)).

We illustrate the asymptotic result (76) by considering the limit $s\to\infty$
of the LHCS spectrum. Taking this limit in (70) gives

$\mathcal{P}(\Theta)=\frac{\alpha}{\sqrt{2\pi}}\mathrm{e}^{-\Theta^{2}/2},\quad\textrm{where}\quad\delta=\sqrt{2/s}$
(60)

and $\alpha$ a constant determined by the dependence of the spectrum on $k$.
Using (60) and the Sokhotski–Plemelj theorem, we rewrite the integral term in
(76) as

$\displaystyle\int_{-\infty}^{\infty}\frac{\mathcal{P}(\Theta)}{\Theta-\Phi_{\pm}\mp\mathrm{i}\mu}\,\mathrm{d}\Theta$
$\displaystyle=\frac{\alpha}{\sqrt{2\pi}}\left(\pm\mathrm{i}\pi\,\mathrm{e}^{-\Phi_{\pm}^{2}/2}+\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int_{-\infty}^{\infty}\frac{\mathrm{e}^{-\Theta^{2}/2}}{\Theta-\Phi_{\pm}}\,\mathrm{d}\Theta\right)$
$\displaystyle=\alpha\left(\pm\mathrm{i}\sqrt{\pi/2}\,\mathrm{e}^{-\Phi_{\pm}^{2}/2}-\sqrt{2}\,\mathrm{daw}(\Phi_{\pm}/\sqrt{2})\right),$
(61)

where
$\mathchoice{{\vbox{\hbox{$\textstyle-$}}\kern-4.86108pt}}{{\vbox{\hbox{$\scriptstyle-$}}\kern-3.25pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-2.29166pt}}{{\vbox{\hbox{$\scriptscriptstyle-$}}\kern-1.875pt}}\\!\int$
denotes the Cauchy principal value and $\mathrm{daw}(\cdot)$ denotes the
Dawson function (DLMF, 2023). Using that
$\mathrm{daw}^{\prime}(x)=1-2x\,\mathrm{daw}(x)$ (DLMF, 2023) we can evaluate
the right-hand side of (76) to find

$\hat{{L}}_{\perp}(\varphi)=\frac{16\alpha}{g\bar{H}_{s}\delta^{2}}\begin{cases}\mathrm{i}\sqrt{\pi/2}\,\Phi_{+}\mathrm{e}^{-\Phi_{+}^{2}/2}+1-\sqrt{2}\Phi_{+}\,\mathrm{daw}(\Phi_{+}/\sqrt{2})&\textrm{for}\
\ 0\leq\varphi<\pi\\\
\mathrm{i}\sqrt{\pi/2}\,\Phi_{-}\mathrm{e}^{-\Phi_{-}^{2}/2}-1+\sqrt{2}\Phi_{-}\,\mathrm{daw}(\Phi_{-}/\sqrt{2})&\textrm{for}\
\ -\pi<\varphi<0\end{cases}$ (62)

Figure 5: Magnitudes of the transfer functions
$\hat{{L}}_{\perp}(\varphi)/H_{s}$ (left) and $\hat{{L}}_{\parallel}(\varphi)$
(right) associated with the vortical and potential part of the current as
functions of $\varphi$ for the LHCS spectrum with directionality parameter
$s=1,\,10$ and $15$. Magnitudes are normalized by $H_{s}$. The exact values
computed from (37)–(38) are shown by the dashed lines; the solid lines in the
left panel show the large-$s$ approximation (62) for
$\hat{{L}}_{\perp}(\varphi)$. We take advantage of the symmetry (20) to show
only the range $\varphi\in[0,\pi]$.

Figure 5 compares the asymptotic approximation (62) of $\hat{{L}}_{\perp}$
with the exact values obtained from (38) for $s=1$, $10$ and $15$. It shows
the asymptotic approximation to be reasonably accurate for $s=10$. We have
checked that the error scales as $O(\delta^{2})$. The figure also shows
$\hat{{L}}_{\parallel}$ in (37) to confirm that
$\hat{{L}}_{\perp}\gg\hat{{L}}_{\parallel}$, and hence that vortical part of
the current dominates over the divergent part, for $s\gg 1$.

Figure 6: SHW $h_{s}$ for the wave spectrum (with $s=10$) and MITgcm current
of figure 1 computed with the full U2H map (left panel, identical to panel (c)
of figure 1), and with its asymptotic approximation for large $s$ (right
panel). This figure can be produced from the notebook accessible at
https://shorturl.at/bef14.

As an application of (62), in figure 6 we compare the predictions of the U2H
map for the MITgcm simulation current of figure 1 computed with the exact
$\hat{\bm{L}}$ and with the asymptotic approximation (62). The match is very
good, even though for $s=10$, $\delta\approx 0.45$ is only marginally small.
The code applying the expression (62) is available on the Jupyter Notebook
https://shorturl.at/bef14, where readers can also experiment with different
choices of the parameter $s$ to observe how the agreements get better/worse
with larger/small $s$.

We conclude by connecting the results of this section with those of Wang et
al. (2023). They focus on the regime $\delta\ll 1$ and on localised currents.
Using matched asymptotics, they obtain an asymptotic expression for the total
SWH $H_{s}=\bar{H}_{s}+h_{s}$ in the presence of currents. This expression
holds without the linearity assumption $h_{s}\ll\bar{H}_{s}$ that underpins
the U2H map. Specifically, they consider the distinguished limit
$\delta=O(\varepsilon)$ which leads to $h_{s}/\bar{H}_{s}=O(1)$. Wang et al.
(2023) give a simplified form valid when $\varepsilon\ll\delta\ll 1$. In
appendix C we show that U2H in the approximation (76) reduces to this form for
localised currents.

##  6 Discussion and conclusion

In the oceanographic regime with $U/c\ll 1$ the effect of currents on an
underlying spatially uniform action spectrum $\bar{\mathcal{A}}({\bm{k}})$ can
be determined by solving the linear problem in (10) and (11) for the anomaly
in action density $a({\bm{x}},{\bm{k}})$. We have focussed on extraction of
the anomaly in SWH $h_{s}({\bm{x}})$ via the weighted ${\bm{k}}$-integral of
$a({\bm{x}},{\bm{k}})$ in (15). The results are in good agreement with
numerical solutions of WW3. There is a significant generalization of this
procedure: given $a({\bm{x}},{\bm{k}})$, other important SGW properties, such
as the current-induced anomaly in the Stokes velocity, are only a
${\bm{k}}$-integral away.

Linearity of the partial differential equation (10) implies that there is
linear map from $\bm{U}$ to $h_{s}$. Can this U2H map be inverted to produce
an H2U map? Not in general: U2H maps a vector field to a scalar field, so
cannot be expected to have an inverse. The non-invertibility of U2H is
illustrated by considering the mildly anisotropic spectrum of §5.2: (54)
implies that $h_{s}=0$ for a velocity field with potential and streamfunction
satisfying $2\psi_{y}-\phi_{x}=0$, demonstrating the non-uniqueness of
$\bm{U}$ for a given $h_{s}$. However, observations of $h_{s}$ provide partial
information about $\bm{U}$. For swell-like waves, in particular, §5.3 shows
that $h_{s}$ can be approximated by a linear operator acting on the vortical
component of the current; in this case, we can infer vorticity from $h_{s}$.

An important qualitative result emerging quickly from the analysis is that the
transfer function, $\hat{\bm{L}}$ in (19), does not depend on the magnitude
$q$ of the current wavenumber $\bm{q}$, but only on its direction $\varphi$.
This implies that the spatial scale of variations in $h_{s}$ is set by those
of the current $\bm{U}$, e.g. power-law $\bm{U}$-spectra result in power-law
$h_{s}$-spectra with the same slope (Villas Boâs et al., 2020). Using (41) we
are now exploring the ramifications of this result.

[Funding]JV and HW are supported by the UK Natural Environment Research
Council (grant NE/W002876/1). ABVB is supported by NASA award 80NSSC23K0979
through the International Ocean Vector Winds Science Team and NASA award
80NSSC24K0411 through the S-MODE Science Team. WRY is supported by the
National Science Foundation award 2048583.

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The customizable Jupyter Notebook that computes
the U2H map from $\bm{U}$ to $h_{s}$ under an LHC incoming wave spectrum is
available on https://shorturl.at/qrxEJ. The configuration files used to
produce the outputs from WW3 used in this work are available at
https://github.com/biavillasboas/U2H.

[Author ORCID]H. Wang, https://orcid.org/0000-0002-5841-5474. A. B. Villas
Bôas, https://orcid.org/0000-0001-6767-6556. W.R Young,
https://orcid.org/0000-0002-1842-3197. J. Vanneste,
https://orcid.org/0000-0002-0319-589X.

##  Appendix A Computable form of the transfer function $\hat{\bm{L}}$

Introducing (31) into the integral in (28) and changing integration variable
from $\theta$ to $\theta+\varphi$ gives

$\int_{0}^{2\pi}\frac{\mathcal{P}(\theta)}{\cos(\theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}\theta=\sum_{n=-\infty}^{\infty}\left(\int_{0}^{2\pi}\frac{\mathrm{e}^{n\mathrm{i}\theta}}{\cos\theta-\mathrm{i}\mu}\,\mathrm{d}\theta\right)\,p_{n}\,\mathrm{e}^{n\mathrm{i}\varphi}.$
(63)

The integrals on the right-hand side can be evaluated explicitly using contour
integration in the complex plane. Consider first $n\geq 0$ and let
$z=\mathrm{e}^{\mathrm{i}\theta}$ so that

$\int_{0}^{2\pi}\frac{\mathrm{e}^{n\mathrm{i}\theta}}{\cos\theta-\mathrm{i}\mu}\,\mathrm{d}\theta=-2\mathrm{i}\oint\frac{z^{n}}{z^{2}-2\mathrm{i}\mu
z+1}\,\mathrm{d}z.$ (64)

We can now use the residue theorem, noting that the integrand has to the two
poles $z_{\pm}=\pm\mathrm{i}+\mathrm{i}\mu+O(\mu^{2})$ and that only $z_{-}$
is enclosed by the integration contour since $\mu\to 0^{+}$. This leads to

$\lim_{\mu\to
0^{+}}\int_{0}^{2\pi}\frac{\mathrm{e}^{n\mathrm{i}\theta}}{\cos\theta-\mathrm{i}\mu}\,\mathrm{d}\theta=2\pi(-\mathrm{i})^{n-1}\quad\textrm{for}\
\ n\geq 0.$ (65)

For $n<0$, we can let $\theta\mapsto-\theta$ on the left-hand side of (64) to
conclude that (65) holds with $n\mapsto-n$. Hence, for any $n$,

$\lim_{\mu\to
0^{+}}\int_{0}^{2\pi}\frac{\mathrm{e}^{n\mathrm{i}\theta}}{\cos\theta-\mathrm{i}\mu}\,\mathrm{d}\theta=2\pi(-\mathrm{i})^{|n|-1},$
(66)

and (63) reduces in the limit $\mu\to 0^{+}$ to

$\lim_{\mu\to
0^{+}}\int_{0}^{2\pi}\frac{\mathcal{P}(\theta)}{\cos(\theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}\theta=2\pi\sum_{n=-\infty}^{\infty}(-\mathrm{i})^{|n|-1}p_{n}\,\mathrm{e}^{n\mathrm{i}\varphi}.$
(67)

Taking the derivative in $\varphi$ lead to the coefficient in front of
$\bm{e}_{\bm{q}}^{\perp}$ in (32).

The computation of the wave momentum $\bm{P}$ in terms of the Fourier
coefficients of the action spectrum is straightforward: using (24), (25), (29)
and (31) we find

$\displaystyle\bm{P}$
$\displaystyle=\int_{0}^{\infty}\\!\\!\int_{0}^{2\pi}\mathcal{A}(k,\theta)k^{2}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\,\mathrm{d}k\mathrm{d}\theta=\int_{0}^{2\pi}\mathcal{P}(\theta)\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\,\mathrm{d}\theta$
$\displaystyle=\sum_{n=-\infty}^{\infty}\int_{0}^{2\pi}\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\,\mathrm{e}^{n\mathrm{i}\theta}\,\mathrm{d}\theta\,p_{n}=\begin{pmatrix}+\mathop{\mathrm{Re}}2\pi
p_{1}\\\ -\mathop{\mathrm{Im}}2\pi p_{1}\end{pmatrix}.$ (68)

Substituting (67)–(68) into (26) yields (32).

##  Appendix B The LHCS spectrum

We use the spectrum proposed by Longuet-Higgins et al. (1963). This takes the
separable form

$\bar{\mathcal{A}}(k,\theta)=f(k)\times\underbrace{\frac{\Gamma(s+1)}{2\sqrt{\pi}\Gamma(s+1/2)}\cos^{2s}(\theta/2)}_{D(\theta)}\,.$
(69)

In (69) the parameter $s\geq 0$ controls the directional spread around the
primary direction of wave propagation conventionally taken along the positive
$x$-axis. The wavenumber function $f(k)$ is chosen so that the frequency
spectrum is a truncated Gaussian which we take centred at $\sigma=0.61$ rad
s-1 corresponding to a peak period of $10.3$ s, with standard deviation of
$0.040$ rad s-1. (See appendix A of Wang et al. (2023) for the complete form
of the spectrum.) The amplitude of the Gaussian, proportional to
$\bar{H}_{s}$, is unimportant in this paper since we take advantage of
linearity to report all numerical results in terms of the relative SWH anomaly
$h_{s}/\bar{H}_{s}$. We vary the sole remaining parameter, $s$, to assess the
impact of wave directionality on $h_{s}$.

With (69), the function $\mathcal{P}(\theta)$ in (29) becomes

$\mathcal{P}(\theta)=\alpha
D(\theta),\quad\textrm{where}\quad\alpha\stackrel{{\scriptstyle\text{def}}}{{=}}\int_{0}^{\infty}\\!\\!\\!f(k)k^{2}\mathrm{d}k.$
(70)

The Fourier coefficients of $D(\theta)$, denoted $p_{n}$ in (31), are real and
satisfy $p_{-n}=p_{n}$. When $s$ is an integer $\Gamma(s+1/2)$ can be
expressed in terms of factorials and then these coefficients are

$p_{n}=\frac{\alpha}{2\pi}\begin{cases}\displaystyle{\frac{(s!)^{2}}{(s+n)!(s-n)!}}&\textrm{for}\
\ |n|\leq s\\\ 0&\textrm{for}\ \ |n|>s.\end{cases},$ (71)

The wave momentum is

$\bm{P}=\frac{\alpha s}{s+1}\begin{pmatrix}1\\\ 0\end{pmatrix}.$ (72)

The spectrum (48) is a particular case of (69) with $s=1$.

##  Appendix C U2H map for highly directional spectra

We derive an asymptotic approximation for $\hat{\bm{L}}(\varphi)$ in (28) for
a spectrum (56) in the limit $\delta\to 0$. We first note that

$\bm{P}=\int_{-\pi/\delta}^{\pi/\delta}\mathcal{P}(\Theta)\begin{pmatrix}\cos(\delta\Theta)\\\
\sin(\delta\Theta)\end{pmatrix}\,\mathrm{d}\Theta=2\pi
p_{0}\begin{pmatrix}1\\\ 0\end{pmatrix}+O(\delta).$ (73)

Thus the first term in (28) makes an $O(1)$ contribution to
$\hat{\bm{L}}(\varphi)$. To evaluate the second term, we approximate the
integral involved as

$\displaystyle\int\frac{\mathcal{P}(\Theta)}{\cos(\delta\Theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}\Theta$
$\displaystyle=\frac{1}{\cos\varphi}\int\mathcal{P}(\Theta)\,\mathrm{d}\Theta+O(\delta)$
$\displaystyle=\frac{2\pi p_{0}}{\cos\varphi}+O(\delta).$ (74)

This is also $O(1)$ except near $\varphi=\pm\pi/2$. There we use the rescaled
variables $\Phi_{\pm}=(\varphi\mp\pi/2)/\delta$ to write

$\displaystyle\int\frac{\mathcal{P}(\Theta)}{\cos(\delta\Theta-\varphi)-\mathrm{i}\mu}\,\mathrm{d}\Theta$
$\displaystyle=\pm\int\frac{\mathcal{P}(\Theta)}{\sin(\delta\Theta-\delta\Phi_{\pm})\mp\mathrm{i}\mu}\,\mathrm{d}\Theta$
$\displaystyle=\pm\delta^{-1}\int\frac{\mathcal{P}(\Theta)}{\Theta-\Phi_{\pm}\mp\mathrm{i}\mu}\,\mathrm{d}\Theta+O(\delta).$
(75)

Introducing (73), (74) and (75) into (28), we find

$\hat{\bm{L}}(\varphi)\sim\frac{16}{g\bar{H}_{s}}\times\begin{cases}-\displaystyle{\delta^{-2}\partial_{\Phi_{+}}\int_{-\infty}^{\infty}\frac{\mathcal{P}(\Theta)}{\Theta-\Phi_{+}-\mathrm{i}\mu}\,\mathrm{d}\Theta\,\,\bm{e}_{\bm{q}}^{\perp}+O(1)}&\textrm{for}\
\ |\varphi-\pi/2|=O(\delta)\\\
\displaystyle{\delta^{-2}\partial_{\Phi_{-}}\int_{-\infty}^{\infty}\frac{\mathcal{P}(\Theta)}{\Theta-\Phi_{-}+\mathrm{i}\mu}\,\mathrm{d}\Theta\,\,\bm{e}_{\bm{q}}^{\perp}+O(1)}&\textrm{for}\
\ |\varphi+\pi/2|=O(\delta)\\\ \displaystyle{-2\pi
p_{0}\left(\frac{\sin\varphi}{\cos^{2}\varphi}\,\bm{e}_{\bm{q}}^{\perp}+\begin{pmatrix}2\\\
0\end{pmatrix}\right)}&\textrm{for}\ \ |\varphi\mp\pi/2|=O(1)\end{cases}.$
(76)

Under the assumption of a localised current, we now derive an approximation
for $h_{s}({\bm{x}})$ by introducing (76) in the inverse Fourier transform

$h_{s}({\bm{x}})=\frac{1}{(2\pi)^{2}}\int\hat{\bm{L}}(\varphi)\bcdot\hat{\bm{U}}(\bm{q})\mathrm{e}^{\mathrm{i}\bm{q}\bcdot{\bm{x}}}\,\mathrm{d}\bm{q}$
(77)

of (16). We compute explicitly the contribution of the boundary layer at
$\varphi=\pi/2$; the contribution of the other boundary layer is its complex
conjugate, denoted by $\mathrm{c.c.}$, and is added. Using the polar
representation of $\bm{q}$ and the approximation

$\bm{q}=q(-\delta\Phi,1)+O(\delta^{2}),$ (78)

where we denote $\Phi_{+}$ by $\Phi$ for simplicity, we obtain

$\displaystyle h_{s}({\bm{x}})\sim$
$\displaystyle\frac{16\delta^{-1}}{(2\pi)^{2}g\bar{H}_{s}}\int_{0}^{\infty}\hat{U}(0,q)q\,\mathrm{e}^{\mathrm{i}qy}\,\mathrm{d}q\int\left(\frac{\mathrm{d}}{\mathrm{d}\Phi}\int\frac{\mathcal{P}(\Theta)}{\Theta-\Phi-\mathrm{i}\mu}\,\mathrm{d}\Theta\right)\,\mathrm{e}^{-\mathrm{i}\delta
qx\Phi}\,\mathrm{d}\Phi$ $\displaystyle+\mathrm{c.c.},$ (79)

Here we implicitly assume that $x=O(\delta^{-1})$ since this turns out to be
the range of $x$ for which $h_{s}$ is the largest. Integrating by parts in
$\Phi$ and swapping the order of the $\Theta$ and $\Phi$ integrations produces

$h_{s}({\bm{x}})\sim\frac{16\mathrm{i}x}{(2\pi)^{2}g\bar{H}_{s}}\int_{0}^{\infty}\hat{U}(0,q)q^{2}\mathrm{e}^{\mathrm{i}qy}\,\mathrm{d}q\int\mathcal{P}(\Theta)\left(\int\frac{\mathrm{d}\Phi}{\Theta-\Phi-\mathrm{i}\mu}\,\mathrm{e}^{-\mathrm{i}q\delta
x\Phi}\right)\,\mathrm{d}\Theta+\mathrm{c.c.}$ (80)

The $\Phi$-integral can be evaluated as a contour integral. For $x\geq 0$, we
close the contour in the lower half plane where the simple pole
$\Phi=\Theta-\mathrm{i}\mu$ is located. In the limit $\mu\to 0^{+}$ the
integral becomes $2\pi\mathrm{i}\mathrm{e}^{-\mathrm{i}q\delta x\Theta}$,
leading to

$h_{s}({\bm{x}})\sim-\frac{8x}{\pi
g\bar{H}_{s}}\int_{0}^{\infty}\hat{U}(0,q)q^{2}\mathrm{e}^{\mathrm{i}qy}\,\mathrm{d}q\int\mathcal{P}(\Theta)\mathrm{e}^{-\mathrm{i}q\delta
x\Theta}\,\mathrm{d}\Theta+\mathrm{c.c.}.$ (81)

Recognising the integral with respect to $\Theta$ as the Fourier transform
$\hat{\mathcal{P}}(q\delta x)$ of $\mathcal{P}(\Theta)$, we obtain the compact
expression

$h_{s}({\bm{x}})\sim-\frac{8x}{\pi
g\bar{H}_{s}}\int_{0}^{\infty}\hat{U}(0,q)\hat{\mathcal{P}}(q\delta
x)q^{2}\mathrm{e}^{\mathrm{i}qy}\,\mathrm{d}q+\mathrm{c.c.}\quad\textrm{for}\
\ x\geq 0.$ (82)

For $x<0$, the integration contour can be closed in the upper half
$\Phi$-plane, where there are no poles, hence $h_{s}({\bm{x}})=0$ for $x<0$.
The reality conditions $\hat{U}^{*}(0,q)=\hat{U}(0,-q)$ and
$\hat{\mathcal{P}}^{*}(q\delta x)=\hat{\mathcal{P}}(-q\delta x)$ can then be
used to reduce the result to:

$h_{s}({\bm{x}})\sim-\frac{8x}{\pi
g\bar{H}_{s}}\int_{-\infty}^{\infty}\hat{U}(0,q)\hat{\mathcal{P}}(q\delta
x)q^{2}\mathrm{e}^{\mathrm{i}qy}\,\mathrm{d}q\quad\textrm{for}\ \ x\geq 0,$
(83)

and $h_{s}({\bm{x}})=0$ for $x<0$.

We now show that (83) is equivalent to the result obtained by Wang et al.
(2023) in the linear regime $\varepsilon\ll\delta\ll 1$. In the notation of
the present paper, their result reads

$h_{s}({\bm{x}})=-\frac{\bar{H}_{s}}{2}\int
D^{\prime}(\Theta)\Delta(y-x\Theta)\mathrm{d}\Theta\quad\textrm{for}\ \ x\geq
0,$ (84)

where $D(\Theta)$ is the (normalised) angular distribution of the spectrum and

$\Delta(y)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{c_{\star}}\int_{-\infty}^{\infty}\\!\\!\\!\zeta(x,y)\,\mathrm{d}x,$
(85)

with $c_{\star}$ the peak group speed. The function $\Delta(y)$ can be written
in terms of the Fourier transform of the current as

$\Delta(y)=-\frac{\mathrm{i}}{2\pi c_{\star}}\int
q_{2}\hat{U}(0,q_{2})\mathrm{e}^{\mathrm{i}q_{2}y}\,\mathrm{d}q_{2}.$ (86)

Introducing (86) into (84)

$h_{s}=-\frac{\bar{H}_{s}x}{4\pi
c_{\star}}\int\hat{U}(0,q)\hat{D}(qx)q^{2}\mathrm{e}^{\mathrm{i}qy}\mathrm{d}q\quad\textrm{for}\
\ x\geq 0.$ (87)

The result in (87) is shown to be identical to (83) by noting that

$\displaystyle\bar{H}_{s}^{2}D(\theta)$
$\displaystyle=16g^{-1}\iint\sigma(k)\mathcal{\bar{\mathcal{A}}}k\mathrm{d}k\mathrm{d}\theta\,D(\theta)=32g^{-1}\iint
c(k)\mathcal{\bar{\mathcal{A}}}k^{2}\mathrm{d}k\mathrm{d}\theta$ (88)
$\displaystyle\approx 32g^{-1}c_{\star}\mathcal{P}(\theta).$ (89)

Here we use the definition (29) of $\mathcal{P}(\theta)$ in the case of a
separable spectrum and we approximate $c(k)\approx c_{\star}$ as appropriate
for a highly directional spectrum.

##  Appendix D Axisymmetric vortex

For a purely vortical, axisymmetric flow, with $\hat{\phi}(q)=0$ and
$\hat{\psi}(q)=-q^{-2}\hat{\zeta}(q)$, combining (39) into (41) yields

$h_{s}(\bm{q})=-\frac{16\mathrm{i}}{g\bar{H}_{s}}\sum_{n=-\infty}^{\infty}n(-\mathrm{i})^{|n|}\,2\pi\tilde{p}_{n}\mathrm{e}^{\mathrm{i}n\varphi}\,q^{-1}\hat{\zeta}(q)\,.$
(90)

The inverse Fourier transform is

$h_{s}({\bm{x}})=-\frac{16\mathrm{i}}{g\bar{H}_{s}}\sum_{n=-\infty}^{\infty}n(-\mathrm{i})^{|n|}\,\frac{2\pi\tilde{p}_{n}}{2\pi}\,\int_{0}^{\infty}\\!\\!\\!\hat{\zeta}(q)\int_{0}^{2\pi}\mathrm{e}^{\mathrm{i}qr\cos(\varphi-\nu)+\mathrm{i}n\varphi}\frac{\mathrm{d}\varphi}{2\pi}\,\mathrm{d}q\,\,.$
(91)

In the $\varphi$-integral above, $\nu$ is the polar angle in physical space.
The $\varphi$-integral can be reduced to Bessel functions

$\displaystyle\int_{0}^{2\pi}\mathrm{e}^{\mathrm{i}qr\cos(\varphi-\nu)+\mathrm{i}n\varphi}\frac{\mathrm{d}\varphi}{2\pi}$
$\displaystyle=\int_{0}^{2\pi}\mathrm{e}^{\mathrm{i}qr\cos\alpha}\cos(n\alpha)\,\frac{\mathrm{d}\alpha}{2\pi}\,\mathrm{e}^{\mathrm{i}n\nu}$
(92) $\displaystyle=\mathrm{i}^{|n|}J_{|n|}(qr)\mathrm{e}^{\mathrm{i}n\nu}\,.$
(93)

Substituting (93) into (91) we obtain

$h_{s}({\bm{x}})=-\frac{16\mathrm{i}}{g\bar{H}_{s}}\sum_{n=-\infty}^{\infty}n\tilde{p}_{n}\,\int_{0}^{\infty}\\!\\!\\!\hat{\zeta}(q)J_{|n|}(qr)\,\mathrm{d}q\,\mathrm{e}^{\mathrm{i}n\nu},$
(94)

valid for any vortical axisymmetric flow. For the Gaussian vortex in (46) the
$q$-integral reduces to

$\displaystyle\int_{0}^{\infty}\\!\\!\\!J_{|n|}(qr)\hat{\zeta}(q)\,\mathrm{d}q$
$\displaystyle=\kappa\int_{0}^{\infty}\\!\\!\\!J_{|n|}(qr)\mathrm{e}^{-r_{\mathit{v}}^{2}q^{2}/2}\,\mathrm{d}q$
(95)
$\displaystyle=\frac{\kappa}{r_{v}}\sqrt{\frac{\pi}{2}}\mathrm{e}^{-r^{2}/4r_{v}^{2}}I_{|n|/2}(r^{2}/4r_{v}^{2}),$
(96)

where $I_{|n|/2}$ is the modified Bessel function. This simplifies
$h_{s}({\bm{x}})$ in (94) to (47).

## References

  * Ardhuin et al. (2017) Ardhuin, Fabrice, Gille, Sarah T, Menemenlis, Dimitris, Rocha, Cesar B, Rascle, Nicolas, Chapron, Bertrand, Gula, Jonathan & Molemaker, Jeroen 2017 Small-scale open ocean currents have large effects on wind wave heights. Journal of Geophysical Research: Oceans 122 (6), 4500–4517.
  * Barber (1949) Barber, Norman F 1949 The behaviour of waves on tidal streams. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 198 (1052), 81–93.
  * Cavaleri et al. (2012) Cavaleri, L, Fox-Kemper, B & Hemer, M 2012 Wind waves in the coupled climate system. Bulletin of the American Meteorological Society 93 (11), 1651–1661.
  * DLMF (2023) DLMF 2023 NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.1.11 of 2023-09-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Hao & Shen (2020) Hao, Xuanting & Shen, Lian 2020 Direct simulation of surface roughness signature of internal wave with deterministic energy-conservative model. Journal of Fluid Mechanics 891, R3.
  * Ho et al. (2023) Ho, Allison, Merrifield, Sophia & Pizzo, Nick 2023 Wave–tide interaction for a strongly modulated wave field. Journal of Physical Oceanography 53 (3), 915–927.
  * Janssen (2004) Janssen, P 2004 The Interaction of Ocean Waves and Wind. Cambridge University Press.
  * Komen et al. (1996) Komen, Gerbrand J, Cavaleri, Luigi, Donelan, Mark, Hasselmann, Klaus, Hasselmann, S. & Janssen, P. A. E. M. 1996 Dynamics and Modelling of Ocean Waves. Cambridge University Press.
  * Longuet-Higgins et al. (1963) Longuet-Higgins, Michael S, Cartwright, David E & Smith, ND 1963 Observations of the directional spectrum of sea waves using the motions of a floating buoy. In Ocean Wave Spectra (ed. Stephan EC), pp. 111–136. Englewood Cliffs, New Jersey: Prentice-Hall.
  * Osborne & Burch (1980) Osborne, AR & Burch, TL 1980 Internal solitons in the Andaman Sea. Science 208 (4443), 451–460.
  * Perry & Schimke (1965) Perry, Richard B & Schimke, Gerald R 1965 Large-amplitude internal waves observed off the northwest coast of Sumatra. Journal of Geophysical Research 70 (10), 2319–2324.
  * Phillips (1977) Phillips, O M 1977 The Dynamics of the Upper Ocean. Cambridge University Press.
  * Romero et al. (2020) Romero, Leonel, Hypolite, Delphine & McWilliams, James C 2020 Submesoscale current effects on surface waves. Ocean Modelling p. 101662.
  * Romero et al. (2017) Romero, Leonel, Lenain, Luc & Melville, W Kendall 2017 Observations of surface wave–current interaction. Journal of Physical Oceanography 47 (3), 615–632.
  * Stein (1970) Stein, Elias M 1970 Singular integrals and differentiability properties of functions. Princeton university press.
  * Tolman (1990) Tolman, Hendrik L 1990 The influence of unsteady depths and currents of tides on wind-wave propagation in shelf seas. Journal of Physical Oceanography 20 (8), 1166–1174.
  * Tolman et al. (2009) Tolman, Hendrik L & others 2009 User manual and system documentation of wavewatch iii tm version 3.14. Technical note, MMAB Contribution 276, 220.
  * Villas Boâs et al. (2020) Villas Boâs, Ana B., Cornuelle, Bruce D., Mazloff, Matthew R., Gille, Sarah T. & Ardhuin, Fabrice 2020 Wave-current interactions at meso and submesoscales: Insights from idealized numerical simulations. Journal of Physical Oceanography 50, 3483–3500.
  * Villas Bôas et al. (2019) Villas Bôas, Ana B & others 2019 Integrated observations of global surface winds, currents, and waves: Requirements and challenges for the next decade. Frontiers in Marine Science 6, 425\.
  * Wang et al. (2023) Wang, H., Villas Boâs, A. B., Young, W. R. & Vanneste, J. 2023 Scattering of swell by currents. J. Fluid. Mech. 975, A1.


# Method of Exact Solutions Code Verification of a Superelastic Constitutive
Model in a Commercial Finite Element Solver

Kenneth I. Aycock Division of Applied Mechanics, Office of Science and
Engineering Laboratories, Center for Devices and Radiological Health, United
States Food and Drug Administration, 10903 New Hampshire Avenue, Silver
Spring, MD 20993 Nuno Rebelo Nuno Rebelo Associates, LLC, Fremont, CA 94539
Brent A. Craven*

###### Abstract

The superelastic constitutive model implemented in the commercial finite
element code ABAQUS is verified using the method of exact solutions (MES). An
analytical solution for uniaxial strain is first developed under a set of
simplifying assumptions including von Mises-like transformation surfaces,
symmetric transformation behavior, and monotonic loading. Numerical
simulations are then performed, and simulation predictions are compared to the
exact analytical solutions. Results reveal the superelasticity model agrees
with the analytical solution to within one ten-thousandth of a percent
(0.0001%) or less for stress and strain quantities of interest when using
displacement-driven boundary conditions. Full derivation of the analytical
solution is provided in an Appendix, and simulation input files and post-
processing scripts are provided as supplemental material.

##  1 Introduction

Superelastic nickel titanium (nitinol) alloys are commonly used in medical
devices such as guidewires, dental arches, and self-expanding peripheral
stents, stent grafts, heart valve frames, and inferior vena cava filters.
Because of nitinol’s unique material behavior and the complex geometry of most
nitinol devices, engineers and scientists often use physics-based
computational modeling and simulation to predict device mechanics and fatigue
safety factors as part of non-clinical bench performance testing. As described
in ASME V&V40-2018 [1], model predictions relied on for decision making should
be accompanied by verification and validation (V&V) evidence demonstrating
simulation credibility commensurate with the risk associated with the intended
model use. However, rigorous code verification evidence for medical device
simulations is often omitted (e.g., see Figure 1 in [2]), in part due to the
lack of detailed examples to facilitate these studies. Herein, we aim to
provide such an example for superelastic nitinol.

In previous work, we demonstrated gold-standard method of manufactured
solutions (MMS) code verification of the commercial finite element software
ABAQUS for various linear and nonlinear elastostatics problems [3]. However,
direct MMS verification of the superelastic model commonly used to simulate
nitinol was not possible due to the lack of a closed-form representation of
the underlying rate- and history-dependent constitutive equations. An approach
recommended in the literature for rigorously verifying similarly complex,
plasticity-based constitutive models is to perform method of exact solutions
(MES) verification on an affine deformation problem with prescribed strains or
displacements [4]. In this study, we perform MES code verification of the
superelastic constitutive model in ABAQUS.

##  2 Methods

###  2.1 Constitutive model summary

Superelastic constitutive behavior was first implemented in ABAQUS/Standard as
a user-material (UMAT) by Rebelo et al [5, 6] in 2000. In brief, the model is
based on the work of Auricchio and Taylor [7, 8] and leverages generalized
plasticity theory to model the dependency of the material stiffness on the
current stress state. More specifically, the model uses a mixture-based
approach to simulate the stress-induced solid-solid phase transformation
between cubic (B2) austenite and monoclinic (B19’) martensite, tracked by the
martensite fraction parameter $\zeta$. Additional details on the constitutive
model and the associated transformation flow rule are provided in the
Appendix. A notional stress-strain response and associated input parameters
are summarized in Figure 1 and Table 1, respectively.

Figure 1: Notional uniaxial stress-strain curve for superelasticity UMAT in ABAQUS/Standard. Table 1: Input parameters for superelastic UMAT in ABAQUS/Standard. parameter |  description
---|---
$E_{a}$ |  Young’s modulus, austenite
$\nu_{a}$ |  Poisson’s ratio, austenite
$E_{m}$ |  Young’s modulus, martensite
$\nu_{m}$ |  Poisson’s ratio, martensite
$\epsilon^{L}$ |  transformation strain
$\left(\frac{\partial\sigma}{\partial T}\right)_{L}$ |  change in transformation stresses with temperature, loading
$\sigma_{L}^{S}$ |  start of transformation, loading
$\sigma_{L}^{E}$ |  end of transformation, loading
$T_{0}$ |  reference temperature
$\left(\frac{\partial\sigma}{\partial T}\right)_{U}$ |  change in transformation stresses with temperature, unloading
$\sigma_{U}^{S}$ |  start of transformation, unloading
$\sigma_{U}^{E}$ |  end of transformation, unloading
$\sigma_{cL}^{S}$ |  start of transformation, compression
$\epsilon_{V}^{L}$ |  volumetric transformation strain

###  2.2 Simplifying assumptions

The superelastic constitutive model in ABAQUS uses pressure-dependent
Drucker–Prager-like transformation surfaces and cubic transformation equations
to define nonlinear hardening during phase transformation between austenite
and martensite. As such, a general analytical solution to the associated rate
equations is not easily obtained and to our knowledge has not been derived.
Here, we instead derive an analytical solution for linear transformation
behavior. Because of the way the nonlinear transformation equations are
defined in ABAQUS, the linear and nonlinear transformation solutions should be
equivalent at the beginning, mid-point, and end of both the loading (austenite
$\rightarrow$ martensite) and unloading (martensite $\rightarrow$ austenite)
transformations under the following assumptions (see Figure 2 and Table 2):

  * •

symmetric transformation behavior in tension and compression, i.e.
$\sigma_{L}^{S}=\sigma_{cL}^{S}$ (constitutive model becomes von Mises-like
rather than Drucker–Prager-like)

  * •

constant temperature (isothermal)

  * •

equal elastic moduli for austenite and martensite, i.e. $E_{a}=E_{m}$

  * •

pseudo-plasticity/superelasticity behavior only (i.e., the model is not
superelastic-plastic)

  * •

monotonic and proportional (i.e., radial) loading.

Figure 2: Illustration comparing the linear and nonlinear transformation solutions. Markers indicate locations at the beginning (1,5), mid-point (2,6), and end (3,7) of the loading and unloading transformation plateaus where the solutions are equivalent, and $\xi$ indicates the martensite fraction. The cubic shape of the nonlinear curve is exaggerated to facilitate visualization. Table 2: Pairs of Mises equivalent stresses $q$ and martensite fractions $\xi$ defining verification points where the linear and nonlinear transformation solutions are equivalent. Note each combination of $q$ and $\xi$ defines a unique location in Figure 2. |  | description | $q$ | $\xi$
---|---|---|---|---
loading | 1 | beginning of transformation | $\sigma_{L}^{S}$ | 0.0
| 2 | mid-point of transformation | $\frac{\sigma_{L}^{S}+\sigma_{L}^{E}}{2}$ | 0.5
| 3 | end of transformation | $\sigma_{L}^{E}$ | 1.0
| 4 | end of loading (pure martensite) | $q_{\textrm{max}}$ | 1.0
unloading | 5 | beginning of transformation | $\sigma_{U}^{S}$ | 1.0
| 6 | mid-point of transformation | $\frac{\sigma_{U}^{S}+\sigma_{U}^{E}}{2}$ | 0.5
| 7 | end of transformation | $\sigma_{U}^{E}$ | 0.0

With these assumptions, we derive an exact analytical solution for the
uniaxial strain of a single cubic element undergoing linear transformation and
monotonic loading.

###  2.3 Problem description

Uniaxial strain of a unit cube is considered, i.e.,

$\displaystyle\epsilon_{11}$ $\displaystyle=f(t)$ (1)
$\displaystyle\epsilon_{22}$
$\displaystyle=\epsilon_{33}=\epsilon_{12}=\epsilon_{13}=\epsilon_{23}=0$ (2)

where $\epsilon_{ij}$ are the components of the logarithmic strain tensor and
$t$ is the simulation pseudo-time for the elastostatic analysis. The cube has
an initial side length $L_{0}=1$ and final side length in the direction of the
applied strain of $L_{F}=L_{0}+u$, where $u$ is the applied displacement
(Figure 3).

Figure 3: Schematic showing uniaxial strain of a unit cube with starting side
length $L_{0}$ and final length $L_{F}=L_{0}+u$.

An unloading analysis is additionally performed to return the single element
to the reference configuration.

During the loading and unloading analyses, the current length $l$ is defined
as a linear ramping between $L_{0}$ and $L_{F}$,

$l=\begin{cases}L_{0}+u\,t&\textrm{loading}\\\
L_{0}+u\left(1-t\right)&\textrm{unloading}\end{cases}$ (3)

with pseudo-time $t$ in the range

$0<t<1\,.$ (4)

We can also quantify the deformation using the non-dimensional stretch,

$\lambda=\frac{l}{L_{0}}=\begin{cases}1+\frac{u\,t}{L_{0}}&\textrm{loading}\\\
1+\frac{u\left(1-t\right)}{L_{0}}&\textrm{unloading.}\end{cases}$ (5)

###  2.4 Analytical solution

The exact analytical solution is derived in detail in the Appendix. In
summary, given

  * •

Mises stress $q$ and martensite fraction $\xi$ at the verification points in
Table 2,

  * •

shear and bulk moduli $G$ and $K$,

  * •

assumptions from Section 2.2,

the analytical solutions for pseudo-time and the remaining field variables are
summarized in Table 3.

Table 3: Analytical solutions for key quantities of interest at verification points $\left(q,\,\xi\right)$ listed in Table 2 and illustrated in Figure 2. description | parameter | expression
---|---|---
total logarithmic strain | $\epsilon$ | $\frac{q}{2G}+\frac{3}{2}\epsilon^{L}\xi$
hydrostatic pressure | $p$ | $-K\epsilon$
deviatoric stress tensor | $\mathbf{S}$ | $\frac{2}{3}q\,\mathbf{n}$
Cauchy stress tensor | $\boldsymbol{\sigma}$ | $\mathbf{S}-p\mathbf{I}$
stretch | $\lambda$ | $\exp\left(\epsilon\right)$
simulation pseudo-time | $t$ | $\begin{cases}\frac{L_{0}(\lambda-1)}{u}&\text{loading}\\\ 1-\frac{L_{0}(\lambda-1)}{u}&\text{unloading}\end{cases}$

###  2.5 Numerical simulations

Single-element verification simulations were performed in ABAQUS/Standard
versions R2016x and 2022 using a single core of an Intel Xeon E5-4627 v4
processor. Commonly used continuum element types C3D8, C3D8R, C3D8I, C3D20,
and C3D20R are investigated in both software versions. Material constants are
prescribed as defined in Table 4. The input files for all cases are provided
in supplemental material.

ABAQUS solves nonlinear mechanics problems in an incremental fashion for each
simulation “step” or load sequence. For the present verification exercise,
multiple steps are defined to facilitate extraction of field outputs at the
precise pseudo-times where the analytical and numerical solutions agree if the
constitutive code is correctly implemented (Tables 2 and 3). Displacement-
controlled boundary conditions are prescribed to enforce the uniaxial strain
conditions described in Eqn. 1–3. Increment sizes are defined such that each
step consists of 100 increments. Numerical calculations of Mises stress $q$,
uniaxial strain $\epsilon_{11}$, and stress components $\sigma_{11}$,
$\sigma_{22}$, and $\sigma_{33}$ are extracted at each pseudo-time point using
an ABAQUS/Python script. Note these results could also be obtained from ABAQUS
printed output provided such output is requested at appropriate simulation
pseudo-times. For elements with multiple integration points, quantities of
interest are averaged across all integration points. Percentage error
magnitudes between the numerical and analytical solutions are calculated as
$|\phi_{\textrm{numerical}}-\phi_{\textrm{analytical}}|/\phi_{\textrm{analytical}}\times
100$ for each quantity of interest $\phi$.

A second set of single-element verification simulations using traction
boundary conditions was also performed. In brief, the moving displacement
boundary condition used above was replaced with a pressure boundary condition,
with the prescribed pressure equal to the negative of the $\sigma_{11}$
component of the Cauchy stress tensor from the analytical solution (Table 3).

Table 4: Input variables for single-element verification problem. description | parameter | value | units
---|---|---|---
initial length | $L_{0}$ | 1 | mm
maximum displacement | $u$ | 0.1 | mm
shear modulus | $G$ | 19,000 | MPa
bulk modulus | $K$ | 42,000 | MPa
Young’s modulus | $E_{\left\\{a,m\right\\}}$ | $\frac{9KG}{3K+G}$ | MPa
Poisson’s ratio | $\nu_{\left\\{a,m\right\\}}$ | $\frac{3K-2G}{2(3K+G)}$ | –
transformation strain |  $\epsilon^{L}=\epsilon_{V}^{L}$* | 0.05 | –
change in transformation stress with temperature | $\left(\frac{\partial\sigma}{\partial T}\right)_{\\{L,U\\}}$ | 0 | $\frac{\textrm{MPa}}{{{}^{\circ}}C}$
start of transformation, loading | $\sigma_{L}^{S}=\sigma_{cL}^{S}$ | 370 | MPa
end of transformation, loading | $\sigma_{L}^{E}$ | 410 | MPa
reference temperature | $T_{0}$ | 0 | °C
start of transformation, unloading | $\sigma_{U}^{S}$ | 160 | MPa
end of transformation, unloading | $\sigma_{U}^{E}$ | 120 | MPa
*The input of a desired volumetric transformation strain $\epsilon_{V}^{L}$ in ABAQUS other than $\epsilon^{L}$ toggles
the use of a non-associated flow rule. Because the model herein simplifies to
Mises-equivalent
transformation, the volumetric transformation strain is automatically zero.

##  3 Results

Each single-element simulation was completed within approximately 30 seconds
with no clear influence of element type or solver version on simulation time
(25.7 $\pm$ 1.3 sec and 35.6 $\pm$ 3.4 sec for displacement- and traction-
controlled simulations, respectively). Simulations converged at each pseudo-
time increment without the need for decreasing the increment size.
Qualitatively, the numerical and analytical solutions are identical across all
element types and software versions considered (Figure 4). The loading and
unloading plateaus exhibit relatively large stiffnesses under uniaxial strain
conditions compared to those associated with uniaxial stress conditions
(Figure 4). Indeed, the $\sigma_{11}$ Cauchy stresses reach values nearly an
order of magnitude larger under uniaxial strain conditions.

Quantitatively, for the displacement-controlled simulations, the results are
relatively insensitive to the continuum element type or software version used,
and extracted stress and strain measures are equivalent to within six to eight
significant figures (see Supplemental Materials). One exception is the final
equilibrium state after unloading where small residual stresses and strains
are observed (Table 5). Comparing the analytical and numerical results,
maximum errors in stress and strain are on the order of one ten-thousandth of
a percent (Table 5). Although differences in the percent error magnitudes are
observed among the various solver versions and element types investigated, the
differences are relatively small. Accordingly, only the largest error
magnitudes are reported in Table 5 for brevity.

Results are similar for traction-controlled conditions, although percent error
magnitudes increase to approximately 1e-03 at the end of transformation in
both loading and unloading. One exception is encountered using C3D8I elements,
where percent error magnitudes reach approximately 1% at the end of
transformation in unloading, and finite residual stresses and strains are
observed at the end of unloading. For full details, see the supplemental
material.

Figure 4: Comparison of the numerical and analytical results under uniaxial strain conditions. The corresponding superelastic stress-strain curve under uniaxial stress is also shown for reference. Table 5: Numerical quantities of interest and percent error magnitudes calculated at the verification time points using displacement-controlled boundary conditions. The final row shows the largest deviations from the anticipated zero-stress and zero-strain condition across the simulations performed. The final three columns likewise report the maximum percent error magnitudes (% err.) observed across all simulations comparing the numerical and analytical solutions. Percent error magnitudes are omitted for the last row given division by zero. Stresses are in units of megapascals (MPa). For full results using each solver version and element type, see supplemental material.  step | $t$ | $q$ | $\epsilon_{11}$ | $\sigma_{11}$ | $\sigma_{\\{22,33\\}}$ |  $\epsilon_{11}$ % err. |  $\sigma_{11}$ % err. |  $\sigma_{\\{22,33\\}}$ % err.
---|---|---|---|---|---|---|---|---
1 | 0.09784 | 370 | 0.009737 | 655.6 | 285.6 | 3.806e-06 | 3.267e-06 | 3.187e-06
2 | 0.4892 | 390 | 0.04776 | 2266 | 1876 | 5.419e-07 | 6.237e-06 | 1.027e-06
3 | 0.8958 | 410 | 0.08579 | 3876 | 3466 | 2.121e-06 | 4.420e-07 | 4.942e-07
4 | 1 | 771.8 | 0.09531 | 4518 | 3746 | 6.165e-06 | 7.726e-06 | 3.672e-06
5 | 0.1757 | 160 | 0.07921 | 3434 | 3274 | 1.168e-06 | 4.990e-07 | 5.234e-07
6 | 0.5796 | 140 | 0.04118 | 1823 | 1683 | 3.009e-06 | 5.874e-06 | 6.362e-06
7 | 0.9684 | 120 | 0.003158 | 212.6 | 92.63 | 1.078e-04 | 1.088e-04 | 1.179e-04
8 | 1 | 3.796e-04 | 1.000e-08 | 6.739e-04 | 2.943e-04 | N/A | N/A | N/A

##  4 Discussion and Conclusions

Using the method of exact solutions, we demonstrate verification of the
superelastic constitutive model implemented in the ABAQUS/Standard implicit
finite element solver. Specifically, uniaxial strain conditions are used to
facilitate derivation of a closed-form solution for monotonic loading and
unloading through the full range of transformation behavior. Although the
uniaxial conditions are relatively simple to implement, they generate a
nontrivial stress state that differs from the uniaxial stress condition
typically used for model calibration. The verification exercise is performed
by extracting quantities of interest at specific simulation increments where
the numerical and analytical solutions should theoretically agree if the
superelastic model is properly implemented. Simulation results reveal maximum
errors in quantities of interest are on the order of one ten-thousandth of a
percent, providing evidence that the model is indeed correctly implemented.
The results are quantitatively similar for all solver versions and continuum
element types investigated.

Note that, while we observe relatively small error magnitudes, the errors
exceed machine precision for the double-precision floating-point operations
performed by ABAQUS/Standard. A potential explanation is numerical
incrementation. Each ABAQUS solution is computed in 800 increments, and there
is a potential for errors to be generated at each increment. The results do
not indicate obvious trends in error magnitude with increasing pseudo-time
increment during loading. Accordingly, errors generated through incrementation
are either negligible compared to the overall observed error magnitudes, or
they are offset by subsequent increments such that they do not accumulate in
simulation pseudo-time. In contrast, error magnitudes do increase i) during
unloading when using C3D8I elements in the 2022 solver under displacement
control, and ii) at the end of transformation for all element types and solver
versions under traction control (see supplemental material), possibly due to
incrementation error. The largest error magnitudes are observed at the end of
unloading when using the C3D8I elements under traction control. Based on
further investigation, this latter observation is uniquely associated with a
combination of C3D8I elements, traction boundary conditions, and the
superelastic model and is believed to be caused by the particular formulation
of the incompatible mode element type.

A few limitations should be noted. First, although the fully integrated and
simplified closed-form solution provided in Table 3 is convenient for
performing the verification exercise, the solution is limited to strictly
monotonic loading. Alternatively, the solution could be generalized to
consider the activation and evolution of the martensite fraction $\xi$ under
other loading paths, for example, unloading from the midpoint on the upper
plateau (i.e., between verification points 2 and 6 in Figure 2). Second, a
number of simplifying assumptions were used to derive the closed-form
analytical solution. Consideration of tension-compression asymmetry,
differences in austenite and martensite elastic moduli, and extensions of the
superelasticity model such as the superelastic-plastic implementation could be
investigated in future work. Third, the verification tests considered herein
only address three-dimensional continuum elements using the implicit solver in
ABAQUS. Verification using other element (e.g., beam and tetrahedral) and
solver types (e.g., ABAQUS/Explicit) could be investigated in future work.

In conclusion, method of exact solutions code verification of the
superelasticity model in ABAQUS was successful under uniaxial strain
conditions. Code verification evidence like that generated by this study,
alongside solution verification, experimental validation, and uncertainty
quantification evidence, are useful to support the credibility of
computational models, especially when model predictions are used to inform
high-risk decision making. To facilitate reproducibility of this study using
other hardware systems or software versions, and adaptation of the approach to
other rate-based constitutive models, the full derivation of the analytical
solution is provided in the Appendix, and simulation input files and post-
processing scripts are provided as supplemental material.

## Conflict of interest statement

One of the authors (NR) was formerly an employee of Dassault Systèmes Simulia,
makers of ABAQUS.

## Acknowledgments

We thank Snehal S. Shetye and Andrew P. Baumann (U.S. FDA) for reviewing the
manuscript. This study was funded by the U.S. FDA Center for Devices and
Radiological Health (CDRH) Critical Path program. The research was supported
in part by an appointment to the Research Participation Program at the U.S.
FDA administered by the Oak Ridge Institute for Science and Education through
an interagency agreement between the U.S. Department of Energy and FDA. The
findings and conclusions in this article have not been formally disseminated
by the U.S. FDA and should not be construed to represent any agency
determination or policy. The mention of commercial products, their sources, or
their use in connection with material reported herein is not to be construed
as either an actual or implied endorsement of such products by the Department
of Health and Human Services.

##  Appendix A Appendix

In the following, standard typeface symbols are scalars (e.g., $\phi$),
boldface symbols denote second-order tenors (e.g., $\mathbf{F}$ or
$\boldsymbol{\sigma}$), and blackboard bold characters denote fourth-order
tensors (e.g., $\mathbb{A}$). Additionally, the overdot operator (e.g.,
$\dot{A}$) indicates a time derivative.

###  A.1 General kinematic equations

Begin with the problem setup provided in Section 2.3 of the manuscript. The
deformation gradient tensor is

$\mathbf{F}=\mathbf{I}+\nabla\mathbf{u}$ (6)

or, in matrix form,

$\mathbf{F}=\begin{bmatrix}\lambda&0&0\\\ 0&1&0\\\ 0&0&1\end{bmatrix}$ (7)

where $\mathbf{I}$ is the second-order identity tensor and
$\lambda=\frac{l}{L_{0}}$ is the stretch. For finite-strain problems, ABAQUS
solves the problem incrementally and calculates the strain by integrating the
strain increments. The resulting strain measure is thus the logarithmic or
“true” strain

$\boldsymbol{\epsilon}=\ln\mathbf{V}$ (8)

where $\ln$ is the principal matrix logarithm and

$\mathbf{V}=\sqrt{{\mathbf{F}\mathbf{F}}^{T}}$ (9)

is the left Cauchy stretch tensor.

Since $\mathbf{F}$ is diagonal here, $\mathbf{F}=\mathbf{F}^{T}$ and
$\mathbf{V}=\mathbf{F}$. Therefore, the logarithmic strain tensor is simply
$\ln\mathbf{F}$ or

$\boldsymbol{\epsilon}=\begin{bmatrix}\ln\lambda&0&0\\\ 0&0&0\\\
0&0&0\end{bmatrix}\,.$ (10)

Since the only non-zero strain component for this uniaxial strain problem is
$\epsilon_{11}$, let

$\displaystyle\epsilon$ $\displaystyle=\epsilon_{11}$
$\displaystyle=\ln\lambda$ (11)

to simplify notation.

The volumetric strain $\epsilon_{V}$ is next defined as

$\displaystyle\epsilon_{V}$ $\displaystyle=\mathrm{tr}\,\boldsymbol{\epsilon}$
$\displaystyle=\epsilon_{11}+\epsilon_{22}+\epsilon_{33}$
$\displaystyle=\epsilon\,,$ (12)

which represents a measure of volume change or dilation, where
$\textrm{tr}\,\mathbf{A}=\mathbf{A}:\mathbf{I}$ is the trace operator for a
second-order tensor and $\left(:\right)$ is the double-inner product operator
$\mathbf{A}:\mathbf{B}=A_{ij}B_{ij}$. Note that the volumetric strain is
sometimes defined in other literature as
$\epsilon_{V}=\frac{1}{3}\textrm{tr}\,\boldsymbol{\epsilon}$, which represents
a measure of mean normal strain.

The strain rate tensor

$\dot{\boldsymbol{\epsilon}}=\begin{bmatrix}\dot{\epsilon}&0&0\\\ 0&0&0\\\
0&0&0\end{bmatrix}$ (13)

is obtained by taking the time derivative of $\boldsymbol{\epsilon}$. The
deviatoric strain rate tensor $\dot{\mathbf{e}}$ is

$\dot{\mathbf{e}}=\boldsymbol{\dot{\epsilon}}-\frac{1}{3}\,\textrm{tr}\,\dot{\boldsymbol{\epsilon}}\,\mathbf{I}\,,$
(14)

or, in matrix form,

$\dot{\mathbf{e}}=\begin{bmatrix}\frac{2}{3}\dot{\epsilon}&0&0\\\
0&-\frac{1}{3}\dot{\epsilon}&0\\\
0&0&-\frac{1}{3}\dot{\epsilon}\end{bmatrix}\,.$ (15)

The scalar equivalent strain rate $\dot{e}$ is

$\displaystyle\dot{e}$
$\displaystyle=\sqrt{\frac{2}{3}\dot{\mathbf{e}}:\dot{\mathbf{e}}}$
$\displaystyle=\frac{2}{3}\dot{\epsilon}\,.$ (16)

###  A.2 Equations for linear, monotonic transformation behavior

The superelastic constitutive model in ABAQUS [5, 6] is based on the work of
Aurrichio and Taylor [7, 8] and leverages generalized plasticity theory to
model the dependency of the material stiffness on the current stress state
(see Online SIMULIA User Assistance 2022 >Abaqus >Materials >Elastic
Mechanical Properties >Superelasticity). The constitutive model uses the
additive strain rate decomposition

$\dot{\boldsymbol{\epsilon}}^{e}=\dot{\boldsymbol{\epsilon}}-\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}\,,$
(17)

where $\dot{\boldsymbol{\epsilon}}^{e}$ is the elastic strain rate tensor and
$\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}$ is the transformation strain rate
tensor. The Cauchy stress rate tensor is then

$\dot{\boldsymbol{\sigma}}=\mathbb{D}:\dot{\boldsymbol{\epsilon}}^{e}$ (18)

where $\mathbb{D}$ is the fourth-order elasticity or stiffness tensor and
$\left(:\right)$ is the operator $\mathbb{A}:\mathbf{B}=A_{ijkl}B_{kl}$. For
an isotropic material,

$\mathbb{D}=2G\,\left(\mathbb{I}-\frac{1}{3}\,\mathbf{I}\otimes\mathbf{I}\right)+K\,\mathbf{I}\otimes\mathbf{I}\,,$
(19)

where $G$ is the shear modulus, $K$ is the bulk modulus,
$\left(\mathbf{I}\otimes\mathbf{I}\right)_{ijkl}=\delta_{ij}\delta_{kl}$,
$\left(\mathbb{I}\right)_{ijkl}=\frac{1}{2}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)$,
and $\delta_{ij}$ is the Kronecker delta function
$\delta_{ij}=\begin{cases}0&i\neq j\\\ 1&i=j\end{cases}$ [9]. Substituting
Eqn. 19 into Eqn. 18 and simplifying, the Cauchy stress rate can be written in
Hooke’s law form as

$\displaystyle\dot{\boldsymbol{\sigma}}$
$\displaystyle=2G\left(\dot{\boldsymbol{\epsilon}}^{e}-\frac{1}{3}\textrm{tr}\,\dot{\boldsymbol{\epsilon}}^{e}\mathbf{I}\right)+K\,\textrm{tr}\,\dot{\boldsymbol{\epsilon}}^{e}\mathbf{I}$
$\displaystyle=2G\dot{\mathbf{e}}^{e}+K\,\textrm{tr}\,\dot{\boldsymbol{\epsilon}}^{e}\mathbf{I}\,.$
(20)

In ABAQUS, the flow rule describing the transformation strain rate for
superelastic materials is

$\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}=\epsilon^{L}\,\dot{\xi}\,\frac{\partial
G^{\textrm{tr}}}{\partial\boldsymbol{\sigma}}\,,$ (21)

where $\epsilon^{L}$ is a material constant, $\xi$ is the martensite fraction,
and $G^{\textrm{tr}}$ is a Drucker–Prager type transformation potential

$G^{\textrm{tr}}=q-p\tan\psi\,.$ (22)

In the above, $p$ is hydrostatic pressure, $\psi$ is a scaling constant, and
$q$ is the von Mises stress

$q=\sqrt{\frac{3}{2}\mathbf{S}:\mathbf{S}}\,,$ (23)

where $\mathbf{S}$ is the deviatoric stress tensor

$\mathbf{S}=\boldsymbol{\sigma}-\frac{1}{3}\,\textrm{tr}\boldsymbol{\sigma}.$
(24)

As stated earlier, here we assume symmetric compression and tension behavior
(i.e., $\sigma_{L}^{S}=\sigma_{cL}^{S}$). Thus, $\psi=0$, and the
transformation potential takes on a von Mises form as simply

$G^{\textrm{tr}}=q\,.$ (25)

###  A.3 An aside: radial return and the direction tensor $\mathbf{n}$

Because the transformation behavior has been simplified as von Mises-type, the
deviatoric stress and deviatoric strain rate tensors point in the same
direction (in 6D space). The transformation potential is $q$, and the
transformation strain rate is proportional to (i.e., in the direction of) the
gradient of $q$ with respect to stress $\frac{\partial
q}{\partial\boldsymbol{\sigma}}$, which we denote as $\mathbf{n}$. Expanding
using Eqn. 23 and applying chain rule,

$\displaystyle\mathbf{n}$ $\displaystyle=\frac{\partial
q}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\frac{\partial\sqrt{\frac{3}{2}\mathbf{S}:\mathbf{S}}}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\frac{\partial\left(\mathbf{S}:\mathbf{S}\right)^{\frac{1}{2}}}{\partial\mathbf{S}}\cdot\frac{\partial\mathbf{S}}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\frac{1}{2}\left(\frac{3}{2}\mathbf{S}:\mathbf{S}\right)^{-\frac{1}{2}}\frac{3}{2}\left(2\mathbf{S}\right)\cdot\frac{\partial(\boldsymbol{\sigma}-\frac{1}{3}\,\textrm{tr}\boldsymbol{\sigma})}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\frac{3}{2}\frac{\mathbf{S}}{q}\cdot(\mathbb{I}-\frac{1}{3}\,\mathbf{I}\otimes\mathbf{I})$
$\displaystyle=\frac{3}{2}\frac{\mathbf{S}}{q}\,.$ (26)

Note the inner product of $\mathbf{n}$ with itself is

$\displaystyle\mathbf{n}:\mathbf{n}$
$\displaystyle=\frac{3}{2}\left(\frac{3}{2}\frac{\mathbf{S}:\mathbf{S}}{q^{2}}\right)$
(27) $\displaystyle=\frac{3}{2}\frac{q^{2}}{q^{2}}$
$\displaystyle=\frac{3}{2}\,.$

The double-inner product between the deviatoric stress tensor $\mathbf{S}$ and
the direction of the deviatoric strain $\mathbf{n}$ is also useful since

$\displaystyle\mathbf{S}:\mathbf{n}$
$\displaystyle=\frac{3}{2}\frac{\mathbf{S}:\mathbf{S}}{q}$
$\displaystyle=\frac{q^{2}}{q}$ $\displaystyle=q\,.$ (28)

Radial return algorithms project the trial stress back onto the yield (here
transformation) surface by scaling the stress radially with respect to the
hydrostatic axis $\sigma_{1}=\sigma_{2}=\sigma_{3}$, where $\sigma_{i}$ are
principal stresses. We assume radial return to be exact under the given
simplifications and approximations, specifically, equal transformation
stresses in compression and tension and thereby von Mises-like transformation
behavior, and proportional (radial) loading. Accordingly, the loading
direction coincides with the projection direction and the (pseudo)plastic
strain rate direction. As derived below,

$\dot{\mathbf{e}}^{tr}=\dot{e}^{tr}\,\mathbf{n}$ (29)

where $\dot{\mathbf{e}}^{tr}$and $\dot{e}^{tr}$ are the deviatoric
transformation strain rate tensor and equivalent scalar, respectively, and
$\mathbf{n}$ specifies the direction of the deviatoric transformation strain
rate (the normal direction to the transformation surface with the given
problem description). Following standard Mises plasticity arguments and given
proportional loading, $\dot{\mathbf{e}}$ and $\dot{\mathbf{e}}^{tr}$are
collinear. Accordingly, using Eqn. 16, the deviatoric strain rate tensor may
also be written

$\displaystyle\dot{\mathbf{e}}$ $\displaystyle=\dot{e}\,\mathbf{n}$
$\displaystyle=\frac{2}{3}\dot{\epsilon}\,\mathbf{n}$ (30)

and therefore, given Eqn. 15,

$\mathbf{n}=\begin{bmatrix}1&0&0\\\ 0&-\frac{1}{2}&0\\\
0&0&-\frac{1}{2}\end{bmatrix}\,.$ (31)

###  A.4 Equations for linear, monotonic transformation behavior (continued)

Continuing from Eqn. 21, using Eqn. 25 and the definition of $\mathbf{n}$, the
last term on the right-hand side becomes

$\displaystyle\frac{\partial G^{\textrm{tr}}}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\frac{\partial q}{\partial\boldsymbol{\sigma}}$
$\displaystyle=\mathbf{n}\,.$ (32)

The transformation strain rate from Eqn. 21 can thus be written

$\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}=\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\,,$
(33)

and the scalar equivalent transformation strain rate becomes

$\displaystyle\dot{\epsilon}^{\textrm{tr}}$
$\displaystyle=\sqrt{\frac{2}{3}\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}:\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}}$
$\displaystyle=\epsilon^{L}\,\dot{\xi}\,.$ (34)

Note that $\mathbf{n}$ is deviatoric and $\textrm{tr}\,\mathbf{n}=0$.
Therefore,

$\textrm{tr}\,\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}=0$ (35)

and

$\displaystyle\textrm{tr}\,\dot{\boldsymbol{\epsilon}}^{e}$
$\displaystyle=\textrm{tr}\,\left(\dot{\boldsymbol{\epsilon}}-\dot{\boldsymbol{\epsilon}}^{\textrm{tr}}\right)$
$\displaystyle=\textrm{tr}\,\dot{\boldsymbol{\epsilon}}$
$\displaystyle=\dot{\epsilon}_{V}\,,$ (36)

where $\dot{\epsilon}_{V}$ is the volumetric strain rate.

Using Eqns. 14, 17, and 33—36, the Cauchy stress rate from Eqn. 20 becomes

$\displaystyle\dot{\boldsymbol{\sigma}}$
$\displaystyle=2G\left[\left(\dot{\boldsymbol{\epsilon}}-\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\right)-\frac{1}{3}\textrm{tr}\,\left(\dot{\boldsymbol{\epsilon}}-\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\right)\mathbf{I}\right]+K\,\dot{\epsilon}_{V}\mathbf{I}$
$\displaystyle=2G\left[\left(\dot{\boldsymbol{\epsilon}}-\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\right)-\frac{1}{3}\textrm{tr}\,\left(\dot{\boldsymbol{\epsilon}}\right)\mathbf{I}\right]+K\,\dot{\epsilon}_{V}\mathbf{I}$
$\displaystyle=2G\left(\dot{\mathbf{e}}-\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\right)+K\,\dot{\epsilon}_{V}\,\mathbf{I}\,.$
(37)

The Cauchy stress rate can be split into deviatoric and hydrostatic components

$\dot{\boldsymbol{\sigma}}=\dot{\mathbf{S}}-\dot{p}\,\mathbf{I}\,,$ (38)

where

$\dot{\mathbf{S}}=2G\left(\dot{\mathbf{e}}-\epsilon^{L}\,\dot{\xi}\,\mathbf{n}\right)$
(39)

is the deviatoric stress rate and

$\dot{p}=-K\,\dot{\epsilon}_{V}$ (40)

is the hydrostatic pressure rate.

The martensite fraction rate $\dot{\xi}$ is a function of the equivalent
stress,

$\dot{\xi}=f(q)\,,$ (41)

and this function is called the transformation law (equivalent to a work-
hardening law in plasticity). In general, the martensite fraction rate must be
calculated for each increment, as its sign and magnitude depend on the current
stress state as well as the direction and magnitude of the stress rate.
However, for the linear hardening approximation and monotonic loading, we can
define the martensite fraction directly as a linear function of the von Mises
equivalent stress. For loading, we have

$\xi_{\textrm{L}}\left(q\right)=\begin{cases}0&q\leq\sigma_{L}^{S}\\\
\frac{q-\sigma_{L}^{S}}{\sigma_{L}^{E}-\sigma_{L}^{S}}&\sigma_{L}^{S}<q<\sigma_{L}^{E}\\\
1&q\geq\sigma_{L}^{E}\end{cases}\,,$ (42)

and for unloading,

$\xi_{\textrm{U}}\left(q\right)=\begin{cases}0&q\leq\sigma_{U}^{E}\\\
\frac{q-\sigma_{U}^{E}}{\sigma_{U}^{S}-\sigma_{U}^{E}}&\sigma_{U}^{E}<q<\sigma_{U}^{S}\\\
1&q\geq\sigma_{U}^{S}\end{cases}\,.$ (43)

For strictly monotonic proportional (constant $\mathbf{n}$) loading or
unloading, we can now integrate the rate tensors in Eqns. 38—40 from $t=0$ to
$t=t$ to obtain

$\boldsymbol{\sigma}=\mathbf{S}-p\,\mathbf{I}$ (44)
$\mathbf{S}=2G\,\left(\mathbf{e}-\epsilon^{L}\,\xi\,\mathbf{n}\right)$ (45)
$p=-K\,\epsilon_{V}\,.$ (46)

Similarly, substituting Eqn. 45 into Eqn. 28 and expanding using Eqns. 30 and
27, the equivalent stress $q$ becomes

$\displaystyle q$ $\displaystyle=\mathbf{S}:\mathbf{n}$
$\displaystyle=2G\left(\mathbf{e}:\mathbf{n}-\epsilon^{L}\,\xi\,\mathbf{n}:\mathbf{n}\right)$
$\displaystyle=2G\left(\frac{2}{3}\epsilon\,\mathbf{n}:\mathbf{n}-\epsilon^{L}\,\xi\,\mathbf{n}:\mathbf{n}\right)$
$\displaystyle=2G\left(\epsilon-\frac{3}{2}\epsilon^{L}\,\xi\right)$
$\displaystyle=2G\,\epsilon-3G\epsilon^{L}\,\xi\,.$ (47)

Note the integrated equations take the same form for both loading and
unloading. Pairings of $q$ and $\xi$, however, define specific locations on
the stress-strain curve such that the solutions for loading and unloading are
unique (see Figure 2 and Table 2).

###  A.5 Final analytical solution for uniaxial strain conditions

Given $q$ and $\xi$ (Table 2), solve Eqn. 47 for the total logarithmic strain
$\epsilon$,

$\boxed{\epsilon=\underbrace{\frac{q}{2G}}_{\text{elastic
strain}}+\underbrace{\frac{3}{2}\epsilon^{L}\xi}_{\text{transformation
strain}}}\,.$ (48)

Next, use Eqn. 26 to calculate the deviatoric stress tensor,

$\boxed{\mathbf{S}=\frac{2}{3}q\,\mathbf{n}}$ (49)

and Eqns. 46 and 12 to calculate the hydrostatic pressure,

$\boxed{p=-K\epsilon}\,.$ (50)

Use Eqn. 44 to calculate the Cauchy stress tensor,

$\boxed{\boldsymbol{\sigma}=\mathbf{S}-p\mathbf{I}}\,.$ (51)

Finally, using Eqns. 5 and 11, calculate the corresponding solution pseudo-
time $t$,

$\boxed{t=\begin{cases}\frac{L_{0}\left(\lambda-1\right)}{u}&\text{loading}\\\
1-\frac{L_{0}\left(\lambda-1\right)}{u}&\text{unloading}\end{cases}}$ (52)

where $\lambda=\exp\left(\epsilon\right)$.

## References

  * [1] ASME V&V40, “Assessing credibility of computational modeling through verification and validation: application to medical devices,” _The American Society of Mechanical Engineers_ , 2018.
  * [2] A. P. Baumann, T. Graf, J. H. Peck, A. E. Dmitriev, D. Coughlan, and J. C. Lotz, “Assessing the use of finite element analysis for mechanical performance evaluation of intervertebral body fusion devices,” _JOR SPINE_ , vol. 4, no. 1, jan 2021.
  * [3] K. I. Aycock, N. Rebelo, and B. A. Craven, “Method of manufactured solutions code verification of elastostatic solid mechanics problems in a commercial finite element solver,” _Computers & Structures_, vol. 229, p. 106175, 2020\.
  * [4] K. Kamojjala, R. Brannon, A. Sadeghirad, and J. Guilkey, “Verification tests in solid mechanics,” _Engineering with Computers_ , vol. 31, no. 2, pp. 193–213, 2015.
  * [5] N. Rebelo and M. Perry, “Finite element analysis for the design of nitinol medical devices,” _Minimally Invasive Therapy & Allied Technologies_, vol. 9, no. 2, pp. 75–80, 2000.
  * [6] N. Rebelo, M. Hsu, and H. Foadian, “Simulation of superelastic alloys behavior with ABAQUS,” in _SMST-2000: Proceedings of the International Conference on Shape Memory and Superelastic Technologies_ , 2001, pp. 457–469.
  * [7] F. Auricchio and R. L. Taylor, “Shape-memory alloys: modelling and numerical simulations of the finite-strain superelastic behavior,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 143, no. 1-2, pp. 175–194, 1997.
  * [8] F. Auricchio, R. L. Taylor, and J. Lubliner, “Shape-memory alloys: macromodelling and numerical simulations of the superelastic behavior,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 146, no. 3-4, pp. 281–312, 1997.
  * [9] E. A. de Souza Neto, D. Peric, and D. R. Owen, _Computational methods for plasticity: theory and applications_. John Wiley & Sons, 2011.