Variability of surface gravity wave field over a realistic cyclonic eddy

Recent remote sensing measurements and numerical studies have shown that surface gravity waves interact strongly with small-scale open ocean currents. Through these interactions, the significant wave height, the wave frequency, and the wave direction are modified. In the present paper, we investigate the interactions of surface gravity waves with a large and isolated realistic cyclonic eddy. This eddy is subject to instabilities, leading to the generation of specific features at both the mesoscale and submesoscale ranges. We use the WAVEWATCH III numerical framework to force surface gravity waves in the eddy before and after its destabilization. In the wave simulations the source terms are deactivated, and waves are initialized with different wave intrinsic frequencies. The study of these simulations illustrates how waves respond to the numerous kinds of instabilities in the large cyclonic eddy from a few hundred to a few tens of kilometres. Our findings show that the spatial variability of the wave direction, the mean period, and the significant wave height is very sensitive to the presence of submesoscale structures resulting from the eddy destabilization. The intrinsic frequency of the incident waves is key in the change of the wave direction resulting from the current-induced refraction and in the location, from the boundary where waves are generated, of the maximum values of significant wave height. However, for a given current forcing, the maximum values of the significant wave height are similar regardless of the frequency of the incident waves. In this idealized study it has been shown that the spatial gradients of wave parameters are sharper for simulations forced with the destabilized eddy. Because the signature of currents on waves encodes important information of currents, our findings suggest that the vertical vorticity of the current could be statistically estimated from the significant wave height gradients down to a very fine spatial scale. Furthermore, this paper shows the necessity to include currents in parametric models of sea-state bias; using a coarse-resolution eddy field may severely underestimate the sea-state-induced noise in radar altimeter measurements.

1 Introduction

The ubiquity of mesoscale (10–100 km) and submesoscale (1–10 km) eddies, fronts, and filaments in the surface of the ocean, leads to a strong variability in the wave field generated by wind (waves): wave–current interactions result in a change of significant wave height (H_s), frequency, and direction (Phillips, 1977 and Mei, 1989).

From these modulations, it has been proven recently, thanks to both field measurements and numerical simulations, that the effects of currents on waves induce strong regional inhomogeneity of the wave field (Romero et al., 2017, 2020). In particular, Ardhuin et al. (2017) showed, using realistic numerical simulations, that the spatial variability of H_s is closely linked to surface kinetic energy (KE) at the mesoscale range. Quilfen et al. (2018), Quilfen and Chapron (2019), and Marechal and Ardhuin (2021) used high-resolution H_s measurements from altimetry and highlight the close link between the surface current gradients and the significant wave height gradients (∇H_s). Villas Bôas and Young (2020) showed numerically, in the absence of wave dissipation and wind momentum input, that the current-induced refraction is necessarily induced by the solenoidal component of the surface currents (vorticity). Finally, Villas Bôas et al. (2020), under the same assumptions, emphasized a relationship between the vertical vorticity of the flow and the ∇H_s.

Beyond a simple modification of the wave kinematics, the surface currents seem to increase the deep-water breaking wave probability and the related air–sea fluxes (Romero et al., 2017, 2020). The reader can refer to the instantaneous numerical outputs of Romero et al. (2020) and notify the local effect of the sharp current gradients on the simulated whitecap coverage (see Fig. 5d and i of the same reference). Wave breaking at the air–sea interface is the major source of momentum and heat exchanges between the atmosphere and the ocean (Cavaleri et al., 2012) or gas and sea spray production (Monahan et al., 1986; Bruch et al., 2021). Therefore, surface mesoscale and submesoscale currents have a significant impact on air–sea fluxes (momentum, gas, heat, sea-spray, etc.) through their interactions with the wave field.

In the ocean and particularly in western boundary currents, eddies are ubiquitous from the mesoscale to the submesoscale range (Chelton et al., 2007, 2011; Gula et al., 2015b; McWilliams, 2016; Rocha et al., 2016). The interactions between the eddy field and the waves are thus of primary importance for the global distribution of wave properties. In the present study, we numerically analyse the effects of an isolated realistic eddy on the wave properties (H_s, mean period, and direction). Former similar works have already been performed, but only for idealized eddy cases (Gaussian profiles; see Mapp et al., 1985; Mathiesen, 1987; White and Fornberg, 1998; Holthuijsen and Tolman, 1991; Gallet and Young, 2014), with a particular attention to the evolution of the wave direction. However, the structure of eddies in the ocean can strongly differ from textbook analytical idealized profiles (Le Vu et al., 2018; de Marez et al., 2019a), making the study of the interactions between the waves and eddy with a Gaussian shape an unrealistic framework. Indeed, the instabilities occurring in a large and isolated eddy result in the strong production of energy in the oceanic submesoscale range (Hua et al., 2013; de Marez et al., 2020b), which would interact strongly with waves. Furthermore, most of the previous studies solely focused on the refraction induced by an eddy without discussing the modulation of other wave parameters (H_s or mean wave period, Mapp et al., 1985; White and Fornberg, 1998; Gallet and Young, 2014). Here, our goal is to investigate the long-term mean effects of an isolated cyclonic eddy with a realistic shape (highly dynamical at the meso- and submesoscale range) on the wave properties. We demonstrate that wave field characteristics are strongly modified by the presence of the eddy and that those changes are even more significant for an eddy field with dynamics in the meso- and the submesoscale range.

This study highlights the importance of working with vortex fields with realistic spatial structures rather than with idealized eddies with a Gaussian shape. For example, in a real ocean, the resulting deviation of the waves from the great circle path due to eddy-induced refraction is certainly underestimated when eddies are considered Gaussian (Gallet and Young, 2014). Also, previous studies in eddy rings in the vicinity of the western boundary current, as in the Gulf Stream, highlighted spatial wave height gradients at the regional scale (Holthuijsen and Tolman, 1991). These spatial gradients would have been strongly underestimated due to the too coarse aspect of the vortex geometry (Gaussian shape). Also, the estimated ocean circulation from altimeter measurements is affected by noise correlated to the H_s (called sea state bias). Some proposed methods to remove the contribution of waves in altimeter measurements assume that the wave field is smooth at scales less than 90 km (Sandwell and Smith, 2005). However, the variability of H_s over a realistic eddy field pattern (more realistic than a Gaussian eddy) reveals very sharp wave parameter gradients of tens to hundreds of kilometres. In a current field, the assumption that the wave field is homogeneous at the mesoscale range is therefore not appropriate. Finally, the signature of currents on waves encodes important information that could be used to infer properties of the underlying current (Huang et al., 1972; Sheres et al., 1985, or more recently Villas Bôas et al., 2020). The last reference showed that sharp ∇H_s can be inverted to infer the vertical vorticity (ζ) field that has generated them. In the similar numerical framework of Villas Bôas et al. (2020), we will show that the statistic of the ζ field can be estimated by inverting the variability of the wave field induced by the eddy field. Reconstructing the ζ field would be relevant for a wide range of applications such as search and rescue, plastic debris monitoring, biological activities, or short-term wave forecast among others.

The paper is organized as follows: in the Sect. 2, we introduce the eddy structure based on the works of de Marez et al. (2020b) and the numerical framework WAVEWATCH III (The WAVEWATCH III® Development Group, 2016) without source terms. Results are presented and discussed in Sect. 3. The limits and the perspectives of this present work close the paper in Sect. 4.

2 Method

2.1 A cyclonic eddy from in situ measurements

To study the wave propagation through an eddy field, we used the current field from the simulation performed by de Marez et al. (2020b). In this study, authors performed idealized simulations, using the Coastal and Regional Ocean COmmunity model, CROCO (Shchepetkin and McWilliams, 2005), that solves the hydrostatic primitive equations for the velocity $\mathit{u}=(u,v,w)$, temperature T, and salinity S, using a full equation of state for seawater (Shchepetkin and McWilliams, 2011). The spatial resolutions are chosen to accurately resolve both the frontal dynamics and the forward energy cascade at the surface. The simulation is initialized with a composite cyclonic eddy as revealed by Argo floats in the northern Arabian Sea (details of the composite extraction are fully described in de Marez et al., 2019a). The eddy is intensified at the surface, but has a deep-reaching influence down to about 1000 m depth. Its initial horizontal shape corresponds to a shielded vorticity monopole: a positive core of vorticity and a shield of negative vorticity (Fig. 1c). Its radius, R=100 km, is large compared to the mean regional Rossby radius R_D (47 km, see Chelton et al., 1998). It is a mesoscale eddy. In the following, mentions of “submesoscale” refer to features and processes occurring at scales that are small compared to Rossby deformation radius (i.e. Bu>1 with $Bu=\frac{R_D^{\mathrm{2}}}{L^{\mathrm{2}}}$, is bigger than 1).

de Marez et al. (2020b) observed that the eddy is unstable with respect to a mixed barotropic–baroclinic instability. The latter deforms the eddy, which eventually evolves into a tripole after about 4 months of simulation. Sharp fronts are subsequently generated in the surface mixed layer at the edge of the tripole. These fronts then become unstable, and this generates submesoscale cyclones and filaments. Near these fronts, diapycnal mixing occurs, causing the potential vorticity to change sign locally and symmetric instability to develop in the core of the cyclonic eddy. Despite the instabilities, the eddy is not destroyed and remains a large-scale coherent structure for 1 year of simulation. A full description of instability processes can be found in de Marez et al. (2020b). Snapshots of the current velocity and vorticity of the fully developed eddy field after 210 d of simulation are represented in Fig. 1b, d respectively. The main core of the cyclone is surrounded by filaments, submesoscale eddies, and fronts that lead to sharp vorticity gradients. This vorticity field is far from the idealized representation of eddies often considered in the literature and is closer to reality (see e.g. Fig. 1 in Lévy et al., 2018, for an example of a realistic turbulent field above mesoscale eddies).

For the purpose of the present study, we consider the surface velocity fields (the simulated level closest to the ocean surface) from the simulation outputs described above. We use the initial state that represents the eddy before instabilities occur (Fig. 1a) and the state after 210 d of simulation, in which submesoscale features have been generated by the spontaneous instability of the eddy (Fig. 1b). At 210 d all instabilities have occurred (mixed barotropic–baroclinic instabilities). After 210 d, the eddy field starts to dissipate, making some small-scale features disappear (de Marez et al., 2020b). Note that the use of strictly 2D surface current is an approximation of what happens in nature. In reality, waves feel the effects of an “average current”, i.e. averaged over the top few metres of the water column. The maximum depth of the current where waves can interact depends on the wavelength of the waves (Kirby and Chen, 1989).

Figure 1 Surface current intensity and direction for the initial and Gaussian eddy (a) and after 210 d of destabilization (b). Their associated normalized relative vorticity ($\mathit{\zeta }={\partial }_{X}V-{\partial }_{Y}U$) are given in (c) and (d). The Coriolis parameter is kept constant in the simulations: $f_{\mathrm{0}}=\mathrm{5.2}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹. The original zonal and meridional velocities simulated in de Marez et al. (2020b) have been multiplied by 2 here.

2.2 The wave model

To describe the dynamics of waves over the eddy described above, we use the WAVEWATCH III numerical framework (The WAVEWATCH III® Development Group, 2016) forced with both the initial state (Fig. 1a, c) and the fully developed state of the eddy (Fig. 1b, d). The model integrates the wave action equation

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\partial }_{t}N(\mathit{\sigma },\mathit{\theta })+& \mathrm{\nabla }.\left(\dot{X}N\right(\mathit{\sigma },\mathit{\theta }\left)\right)+{\partial }_k\left(\dot{k}N\right(\mathit{\sigma },\mathit{\theta }\left)\right)\\ +& {\partial }_{\mathit{\theta }}\left(\dot{\mathit{\theta }}N\right(\mathit{\sigma },\mathit{\theta }\left)\right)=S,\end{array}\end{array}$$

where N(σ,θ) is the wave action spectrum ($N(\mathit{\sigma },\mathit{\theta })=\frac{E(\mathit{\sigma },\mathit{\theta })}{\mathit{\sigma }}$, with E(σ,θ) the two-dimensional wave energy spectrum), θ the direction of wave propagation, σ the wave intrinsic frequency equal to $\sqrt{gk}$ in deep water (where water depth is largely greater than wave wavelength; here k is the wavenumber and is a scalar), and g is the gravitational acceleration. $\dot{X}$ is the wave action advection velocity (equal to the sum of the wave group velocity and the surface current velocity), and $\dot{k}$ and $\dot{\mathit{\theta }}$ are the wave advection velocities in the spectral space. The expressions of $\dot{k}$ and $\dot{\mathit{\theta }}$ are developed from wave ray equations (Eq. 3) and are fully given in Phillips (1977), Benetazzo et al. (2013), and Ardhuin et al. (2017). The right-hand side of Eq. (1) is the sum of the source terms describing the wind energy input, the dissipation due to the wave breaking and bottom friction, and the nonlinear energy exchange between waves.

The dispersion of the waves is described by the dispersion relationship that links σ and k. In a current field, it is necessary to consider a moving frame of reference. The waves' dispersion relationship is thus impacted because the current induces a Doppler shift on the wave frequency (Eq. 2),

$$\begin{array}{}\text{(2)}& \mathit{\omega }=\mathit{\sigma }+\mathit{k}.\mathit{u}.\end{array}$$

The wave ray equation is also modified:

$$\begin{array}{}\text{(3)}& {\partial }_t\mathit{k}={\partial }_x\mathit{\omega }.\end{array}$$

ω is the absolute frequency, k the wavenumber vector, and u the surface current vector. Bold characters refer to vector notation all along this paper.

For this study we consider waves already well developed, far from their generation areas, propagating in the current field without any source term (no dissipation, no nonlinear exchanges between waves, and no wind input; i.e. the right-hand side of Eq. 1 is equal to 0). The aim of the current study is to investigate, in a very idealized case, how long wave properties can be modified by an isolated eddy field more realistic than a Gaussian eddy. In a more realistic framework, the wave steepness modified by the current or due to nonlinear wave–wave interactions would lead to local wave breaking as observed in Romero et al. (2017). Also, wind input would generate higher-frequency waves, which will also interact with the eddy field.

Throughout this paper we discuss the evolution of the significant wave height (H_s) and the mean wave period weighted on the low-frequency part of the wave spectrum ($T_{\text{m}\mathrm{0},-\mathrm{1}}$), known as “bulk” quantities. We called them “bulk” because they are integrated over the wave energy spectrum, E(σ,θ). They are defined as

$$\begin{array}{}\text{(4)}& H_{\text{s}}=\mathrm{4}\sqrt{\underset{\mathit{\theta }=\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{{\mathit{\sigma }}_{min}}{\overset{{\mathit{\sigma }}_{max}}{\int }}E(\mathit{\sigma },\mathit{\theta })d\mathit{\sigma }d\mathit{\theta }}\end{array}$$

and

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{T}_{\text{m}\mathrm{0},-\mathrm{1}}=& {\displaystyle \frac{\mathrm{1}}{{\int }_{\mathit{\theta }=\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{{\mathit{\sigma }}_{min}}^{{\mathit{\sigma }}_{max}}E(\mathit{\sigma },\mathit{\theta })d\mathit{\sigma }d\mathit{\theta }}}\underset{\mathit{\theta }=\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{{\mathit{\sigma }}_{min}}{\overset{{\mathit{\sigma }}_{max}}{\int }}\\ \\ & {\mathit{\sigma }}^{-\mathrm{1}}E(\mathit{\sigma },\mathit{\theta })d\mathit{\sigma }d\mathit{\theta }.\end{array}\end{array}$$

The evolution of the wave peak direction (θ_p, θ where E(σ,θ) is maximum) is also studied. The performance of the wave model used here has already been discussed in boundary current systems such as in the Gulf Stream, the Drake Passage, and the Agulhas Current, especially concerning the H_s estimation (Ardhuin et al., 2017; Marechal and Ardhuin, 2021). In those previous studies, wind forcing, wave dissipation, and nonlinear wave–wave interactions have been taken into account.

We initialize simulations with waves that are propagating from the left boundary of a 500 km × 500 km Cartesian domain, with a resolution of 500 m in both X and Y. The right boundary is open. The initialization is done with narrow-banded wave spectra that are Gaussian in frequency and centred at different peak frequencies, f_p=0.1428, 0.097, and 0.0602 Hz. The frequencies have been chosen to correspond to the mean periods used in the work of Villas Bôas et al. (2020) (7, 10.3, and 16.6 s). The wave spectra have a frequency spread of 0.03 Hz around the peak frequency, and an initial H_s equals 1 m. Waves are generated at the left boundary hourly from the spectra described above. The initial direction of waves is 270^∘. The direction convention follows the meteorological convention such that 270^∘ waves are coming from the left and are propagating toward the right boundary parallel to the x axis. The wave model global time step is 12 s, the spatial advection time step is 4 s, and the spectral time step is 1 s. The model provides outputs every 15 min. Wave spectra are computed at each grid point, discretized into 32 frequencies and 48 directions. Fine directional resolution is required for a better description of wave refraction, especially in strong rotational currents (Ardhuin et al., 2017; Marechal and Ardhuin, 2021). The surface current forcing fields are from the de Marez et al. (2020b) simulation outputs. In one case we consider the initial shape of the cyclonic eddy (Fig. 1a, c). In the other case, we consider the fully developed state of the cyclonic eddy (Fig. 1b, d). The variation timescale of the current is much longer (𝒪(1) week) than the waves (𝒪(1) min); thus the current is assumed to be stationary during one wave train propagation. The simulations forced with the initial eddy are similar to the former works performed for Gaussian eddies (Mathiesen, 1987; Holthuijsen and Tolman, 1991; White and Fornberg, 1998; Gallet and Young, 2014).

The eddy described in the previous section and in de Marez et al. (2020b) is an averaged composite eddy reconstructed from measurements in the Arabian Sea (de Marez et al., 2019a). The method of reconstruction tends to an underestimation of the eddy intensity, which is why both the zonal and meridional velocities have been multiplied by 2 in order to increase the potential effects of currents on wave properties. The eddy is staying geophysically realistic (current velocity remains around 1 m s⁻¹ and normalized vorticity lower than 2; see Fig. 1). Those values are comparable with surface vorticity measured in the first hundred metres of the Arabian Sea (de Marez et al., 2020a) and in other current regimes as in the western boundary currents (Gula et al., 2015a; Tedesco et al., 2019). Although the eddy field represented in Fig. 1b, d is from an averaged composite eddy (solely estimated using in situ data), it has been considered, in this study, to be realistic because it differs from an analytical vortex. Also, the fully developed eddy has been compared with altimeter and drifter data in the Arabian Sea region where it has been estimated. The cyclonic eddy was coherent with those measurements (see Figs. 12, 13, and 14 of de Marez et al., 2019a).

3 Results

The frequency sensibility of the incident waves is studied in both the initial and the fully developed eddies. Waves are dispersive in deep water; their group and their energy propagate at the group velocity (C_g). For T_p=7 s ($T_{\text{p}}=\frac{\mathrm{1}}{f_{\text{p}}}$), T_p=10.3, and T_p=16.6 s, group velocities are 11, 16, and 26 m s⁻¹. To reach X=X₀ (a given value of X), short waves take more time than long waves. As waves are generated continuously from the left boundary, a stationary state is reached. The wave field reaches the stationary state after 10 h, 9 h, and 8 h of simulation for initializations of T_p=7, T_p=10.3, and T_p=16.6 s incident waves respectively. In Figs. 2, 3, and 4, the fields are taken once the stationary state is reached. Surface currents modulate the wave amplitude, the wave frequency, and the wave direction; the variability of these wave properties is highlighted through the H_s, $T_{\text{m}\mathrm{0},-\mathrm{1}}$, and θ_p variables. Other aspects of waves’ variability, e.g. directional spread or mean direction, are not described here.

Figure 2 Significant wave height (H_s) fields for (a, d) T_p=7 s, (b, e) 10.3 s, and (c, f) 16.6 s incident waves. Without current forcing, the entire domain is equal to the initial H_s (1 m). The first row (a, b, c) shows fields for simulations forced with the initial eddy (Fig. 1a, c); the second row (d, e, f) shows the same fields but for simulations forced with the fully developed eddy (Fig. 1b, d). Panel (g) shows H_s along X=300 km.

3.1 Modulation of wave parameters

3.1.1 Significant wave height

Surface currents induce a strong regional H_s variability, especially in a highly solenoidal field (Ardhuin et al., 2017; Villas Bôas et al., 2020; Marechal and Ardhuin, 2021). The presence of the vortex induces strong significant wave height gradients ($\left|\mathrm{\nabla }{H}_{\text{s}}\right|=\sqrt{({\partial }_X{H}_{\text{s}}{)}^{\mathrm{2}}+({\partial }_Y{H}_{\text{s}}{)}^{\mathrm{2}}}$, noted ∇H_s hereinafter), inside and outside the eddy (Fig. 2). Simulations forced with the initial eddy (Fig. 2a, b, c) show coherent alternate sign H_s structures along lines of constant X. An important lens shape dipole of H_s increase and decrease is noticeable in the field. H_s reaches a maximum of 1.63 m at X=308 km, 1.62 m at X=324 km, and 1.57 m at X=340 km for simulations initialized with T_p=7, T_p=10.3, T_p=16.6 s waves respectively. A transect at X=300 km is given for every initialization in Fig. 2g. Two maximums are visible, the main one at Y=310 km (H_s∼1.6 m) and a secondary one at Y=125 km (H_s∼1.2 m). Two minimums are visible, one at Y=200 km (H_s=0.8 m) and a secondary one near Y=380 km (H_s=0.85 m). One can see that, at Y=200 km (300 km), shorter incident waves result in lower (higher) H_s. Globally, H_s follows the current vorticity signal (Fig. 1c). The enhanced H_s areas are associated with the boundary of the inner eddy core (ζ>0) and the vorticity ring (ζ<0) that surrounds the eddy core. Where waves are propagating against the current, H_s is enhanced, which agrees with wave–eddy interactions simulated in realistic fields (see Fig. 1 of Ardhuin et al., 2017, and Fig. 6 of Romero et al., 2020).

Simulations forced with the fully developed eddy show a stronger spatial inhomogeneity in the wave field (Fig. 2d, e, f). The initial H_s is more scattered (mostly in the X direction due to the initial direction of the incident wave packet) than in the initial eddy. As noticed for simulations forced with the initial eddy (Fig. 2a, b, c), the H_s field matches pretty well with the current forcing (Fig. 1b, d); in other words where surface current gradients are important, strong ∇H_s values are noticed. H_s is mostly modulated by the fully developed eddy core. The modulation of H_s by the fully developed eddy core occurs for X>50 km, which is more upstream than the H_s modulations induced by the initial eddy. Let us note that ∇H_s values are apparent in the submesoscale eddies that have emerged spontaneously all around the eddy core. In the submesoscale eddy field, the wave field shows alternate sign fluctuations of H_s, with globally the same intensity regardless of the period of incident waves. It is explicitly shown in Fig. 2g at Y<180 and Y>350 km for every initialization. In the same transect, at Y=200 km, as previously, for shorter incident waves, the H_s values are lower. However at Y=300 km, the ∇H_s values are almost identical regardless of the periods of the incident waves, and the ∇H_s values along Y are strongly sharper for simulations forced with the fully developed eddy with higher extreme values. One can see that ∇H_s values are important downstream of the eddy field. The horizontal size of H_s patches (intensified or decreased H_s structures) is comparable to the width of the eddy (Fig. 2a–f). Finally one can see that for all simulations, the signatures of the eddy in the H_s field are not totally symmetric with respect to the Y axis, whereas the current fields used to force the wave model seemed to be so.

The intensity and the patterns of ∇H_s are very sensitive to the underlying current: the more turbulent the vortex, the sharper the ∇H_s (Fig. 2). Figure 2g shows that, at X=300 km, the (minimum) maximum values of H_s are (lower) higher for the fully developed eddy but are very similar regardless of the periods of the incident waves. For the two current forcings and all initializations of the model, we computed the 95th percentile of the H_s values, the maximum value of H_s, and the distance from the left boundary where the maximum value of H_s is located. The results are given in Table 1. Regardless of the periods of the incident waves, the 95th percentile of H_s is similar for the two current forcings and varies between 1.18 and 1.24 m, with a maximum of the 95th percentile of H_s for simulation initialized with 10.3 and 16.6 s and forced with the fully developed eddy. The maximum values of H_s are higher for the simulations forced with the fully developed eddy. Finally, the shorter the incident waves, the closer the maximum H_s values are to the left boundary, with a minimum distance for the simulation forced with the fully developed eddy and initialized with T_p = 7 s.

Table 1 The 95th percentile of the significant wave height (H_s), the maximum value of H_s, and the distance from the left boundary where the maximum value of H_s is located. T_p is the peak period of the incident waves.

3.1.2 Peak direction

The effect of currents on wave direction can be captured to the first order by the θ_p field. Waves turn in the current field due to the refraction induced by the vorticity of the flow (Kenyon, 1971; Dysthe, 2001). Waves turn toward Y=0 km (θ_p increase) in the bottom part of the domain and toward Y=500 km (θ_p decreases) in the upper part (Fig. 3). When waves pass through the eddy, θ_p changes due to the vorticity field, at X=125 km for the initial eddy (Fig. 3a, b, c), and slightly upstream, at X=79 km, for the fully developed eddy (Fig. 3d, e, f). Patterns shown in Fig. 3 are similar to the ∇H_s patterns shown in Fig. 2 with a large-scale dipole for simulations forced with the initial eddy and both large- and small-scale signal gradients for simulations forced with the fully developed eddy. The peak direction gradient intensity ($\left|\mathrm{\nabla }{\mathit{\theta }}_{\text{p}}\right|=\sqrt{({\partial }_X{\mathit{\theta }}_{\text{p}}{)}^{\mathrm{2}}+({\partial }_Y{\mathit{\theta }}_{\text{p}}{)}^{\mathrm{2}}}$, noted ∇θ_p hereinafter) depends on both the period of the incident waves and the underlying vorticity field (Dysthe, 2001; Kenyon, 1971). ∇θ_p is stronger for simulations initialized with T_p=7 s (Fig. 3a, d) than for simulations initialized with T_p=10.3 and 16.6 s, with the sharpest gradients for simulation forced with the fully developed eddy (Fig. 3d). In this simulation waves can be deviated by 30^∘ with respect to the initial direction of the waves. The result corroborates the findings of Villas Bôas et al. (2020) where authors forced wave model with synthetic surface currents inverted from the kinetic energy spectrum (with a random phase) with different spectral slopes. The more turbulent the current is, the more the waves are refracted. Very long wave trains (T_p=16.6 s) hardly reach a deviation of wave direction higher than 10^∘, in both the fully developed and the initial eddies. Finally, one can see that θ_p differs downstream of the eddy with respect to the initial direction (270^∘). Downstream of the eddy field, the wave fields retain the effects of the remote currents that passed through.

Figure 3 Peak direction (θ_p) fields for (a, d) T_p=7 s, (b, e) 10.3 s, and (c, f) 16.6 s incident waves. Without current forcing, the entire domain is equal to the initial θ_p (270^∘). The first row (a, b, c) shows θ_p for simulations forced with the initial eddy (Fig. 1a, c); the second row (d, e, f) shows the same fields but for simulations forced with the fully developed eddy (Fig. 1b, d). The narrow yellow bands in the left part of every panel are spurious; they marked the boundary where waves are generated at the left boundary.

3.1.3 Mean wave period

The surface currents have an effect on the wave frequency (Phillips, 1977). Due to the conservation of the absolute frequency in a current field (ω in Eq. 2), the intrinsic frequency (σ) is modified, which subsequently changes the $T_{\text{m}\mathrm{0},-\mathrm{1}}$ (Eq. 5). Wave simulations are initialized with different wave peak frequencies, so $T_{\text{m}\mathrm{0},-\mathrm{1}}$ is directly impacted. The different initializations of the wave field justify the representation of the relative difference of $T_{\text{m}\mathrm{0},-\mathrm{1}}$ ($\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$) rather than the raw outputs. $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ is the difference between the outputs of simulations performed with and without surface current forcing. The results are given in Fig. 4. At first glance, the spatial variability is more striking for simulations forced with the fully developed eddy with patterns similar to the H_s fields (Fig. 2). For the fully developed eddy, $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ exceeds 3 s in the eddy core for X between 200 and 400 km. For the initial eddy, for all the initializations, $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ does not exceed 2 s (Fig. 4a, b, c). Similarly to the H_s, the $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ does not much depend on the period of the incident waves, or at least, not as much as the θ_p fields studied above. Slight differences are, however, noticeable for simulations forced with the fully developed eddy. It is not clear if there is a link between the wave period of the incident waves and the slight differences in $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ shown in Fig. 4g both in the main eddy structure and in the submesoscale eddies. Indeed, $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ values are stronger for long incident waves (T_p=16.6 s) in the submesoscale eddies, whereas we see the opposite in the core of the fully developed eddy (T_p=10.3 s).

$\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ is positive where waves and current are aligned and negative where waves and current are opposed. This change of $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ is because the current induces a Doppler shift on the wave frequency (Eq. 2) and the absolute frequency is conserved. If we focus on the maximum of $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$, at Y=200 km, wavelengths increase to about 153 m and H_s values decrease by ∼0.65 m. Where waves and current are opposite, we see that H_s values are enhanced (Fig. 2) and wave wavelengths are shortened and vice versa. It is due to the conservation of wave action (D_tN=0, Eq. 1). One can see that stripe structures induced by the refraction (Fig. 3) are also visible through the mean wave period fields.

We recall that the change of H_s induced by current is due to a superposition of processes. Indeed, in current field, in the absence of wind, the regional ∇H_s results mainly from the current-induced refraction and the advection of wave action by the current and the group speed (Ardhuin et al., 2017). The current-induced changes in the wave frequency can also increase the H_s (see introduction of Benetazzo et al., 2013). Note that current refracts waves such that waves and current can become aligned (or opposite). So refraction can lead to a change of mean wave period downstream from the refraction areas in the same manner that refraction induces a non-local change of H_s.

Figure 4 Mean wave period difference ($\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$) between simulations forced with and without current ($\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}=T_{\text{m}\mathrm{0},-\mathrm{1}}\left(\text{curr}\right)-T_{\text{m}\mathrm{0},-\mathrm{1}}(\text{Nocurr}$)). Panels (a, d) show $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ fields initialized with T_p=7 s waves. Panels (b, e) show $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ fields initialized with T_p=10.3 s. Panels (c, f) show $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ fields initialized with T_p=16.6 s. The first row (a, b, c) shows fields for simulations forced with the initial eddy (Fig. 1a, c); the second row (d, e, f) shows the same fields but for simulations forced with the fully developed eddy (Fig. 1b, d). Panel (g) shows $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ along X=300 km.

For all the variables studied here (Figs. 2, 3, 4), waves are continuously generated from the left boundary, a solitary incident wave train strongly affects the results presented above; for instance the non-local effect of refraction on the wave field is strongly less pronounced (not shown).

3.2 Ray tracing

Knowing that the wave action (N(σ,θ)) is conserved along the wave trajectory in the current field (Bretherton and Garrett, 1968), we show in this section, from a ray-tracing framework, that waves respond very differently to the two eddy fields. In the present study, the isolated vortex refracts the waves and changes the wave frequency, which leads to a strong inhomogeneity in both the H_s and $T_{\text{m}\mathrm{0},-\mathrm{1}}$ fields (Figs. 2, 4). The current-induced refraction is highlighted, here, thanks to Monte Carlo ray-tracing simulations, as performed in literature for different current regimes (White and Fornberg, 1998; Ardhuin et al., 2012; Gallet and Young, 2014; Kudryavtsev et al., 2017; Villas Bôas and Young, 2020). This ray-tracing method is used in order to follow the conservation of the wave action in the current field. For the ray tracing, we assume that the surface current is stationary $\left(\frac{\left|\mathit{u}\right|}{C_{\text{g}}}\ll \mathrm{1}\right)$ and that incident waves are monochromatic. In the real ocean, the wave field is a superposition of wave trains with specific directions and frequencies; thus ray tracing is only a very simplified view of how the direction of the waves is modified by the presence of current.

For the ray-tracing model calculations, the initial direction is 270^∘ (waves are coming from the left boundary), and the initial frequencies are the same as the ones discussed above (T_p=7, 10.3, and 16.6 s peak periods). We see that the current-induced refraction is sensitive to both the nature of the underlying current and the frequency (or wavelength) of the incident waves (Fig. 5). The radius of curvature of wave rays is larger where the current field is highly rotational (Fig. 5d, e, f) and when the ray-tracing simulations are initialized with T_p=7 s waves (Fig. 5a, d). This confirms the works of Kenyon (1971) and Dysthe (2001). In the initial eddy, the wave train is refracted by both the eddy's edge (toward the lower part of the domain) and the core of the eddy (toward the upper part of the domain; Fig. 5a, b, c). It leads to two wave ray focalization areas downstream of the initial eddy. These focalization areas, or caustics, are slightly shifted toward the right boundary when the incident waves are longer. Indeed, the caustic in the upper part of the Fig. 5a, b, c appears at X=330, X=370, and X=445 km respectively. The locations of caustic formation appear further downstream of the eddy than the location of the maximum values of H_s (Table 1). However, the positions of the caustics are proportional with respect to the location of the maximum values of H_s; i.e. the shorter the incident waves, the closer the caustic to the left boundary.

In the fully developed eddy field, both mesoscale and submesoscale eddies refract waves. In comparison with the initial eddy, one can see that the number of caustics increases in the fully developed eddy, with a maximum of caustics for T_p=7 s incident waves (Fig. 5d). Even if isolated submesoscale eddies have a vorticity comparable with the eddy core ($\frac{\mathit{\zeta }}{f_{\mathrm{0}}}\sim \mathrm{1.5}$), they do not refract waves as much as the centre structure. Indeed, if we look at the southernmost submesoscale eddy, we see that one wave ray deviates about 30 km from the left boundary to the right boundary, whereas one wave ray at the centre of the domain deviates more than 200 km. So, the shape of vorticity patterns is key in the intensity of the refraction. One can note that the ray convergent areas are also localized almost where the maximum values of H_s are spotted (Fig. 2), especially at the edge of the positive vorticity core. The main caustic at Y=300 km is slightly shifted toward the right boundary for longer incident waves, which is also qualitatively consistent with results shown in Table 1.

In a strong rotational current field, the change of H_s is mostly driven by refraction induced by mesoscale and submesoscale currents (Irvine and Tilley, 1988; Ardhuin et al., 2017; Romero et al., 2020). This has been confirmed in additional simulations where the refraction has been deactivated, showing maximum values of H_s not exceeding 1.36 m in the core of the eddy field (not shown). With realistic numerical studies in strong current fields, Ardhuin et al. (2012) and Kudryavtsev et al. (2017) qualitatively showed the link between caustics and areas where H_s is enhanced. Here, the ray-tracing model highlights the current-induced refraction but it can also explain how the surface currents induce H_s variability. If we assume that one ray is carrying a certain amount of wave action with a certain value of H_s (here 1 m), caustic locations can be assimilated to areas of wave action accumulation and, subsequently, assimilated to areas of increases in H_s. If we consider an infinite number of rays, the expected H_s at caustic locations is infinite. However, in a real ocean, because the wave action is distributed in a range of frequencies and directions, these H_s enhancement are limited. In the fully developed eddy, there are more caustics than in the initial eddy due to the submesoscale eddies, which could explain why the H_s fields presented in Fig. 2d, e, f show more ∇H_s structures. Also, this partially explains why the extreme values of H_s are very slightly higher for simulations initialized with short waves (Table 1).

Figure 5 (a, b, c) Ray tracing for waves travelling over the initial eddy with T_p=7 s, (a) 10.3 s (b), and 16.6 s (c) peak period. Panels (d, e, f) show the same ray tracing but for waves travelling over the fully developed eddy. The vorticity fields are given in the background.

The strong vorticity field, for both the initial and the fully developed cyclonic eddies, induces wave ray scattering, which can reach a deviation of several hundred kilometres in comparison with simulations without background current. This deviation is more important for short wave incidence (Fig. 5a, d). In the ocean, the strong wave scattering can be responsible for the space-time bias in the forecast of waves' arrival (Gallet and Young, 2014; Smit and Janssen, 2019).

The present ray-tracing simulation shows that refraction has a local effect on wave direction; strong ray deviations appear where surface gradients are strong. However, refraction effects on wave parameters are non-local. Indeed, the sharp ∇H_s areas seem to be associated with wave ray caustics and can appear both inside and outside the eddies (Figs. 2, 5). In other words, strong ∇H_s is not necessarily located where strong surface currents are spotted.

3.3 Is it possible to reconstruct ∇U via the measurement of the ∇H_s?

We have seen that the current-induced refraction and ∇H_s are driven by the underlying turbulence induced by the presence of the cyclonic eddy. Villas Bôas et al. (2020) and Marechal and Ardhuin (2021) showed that at scales between 200 and ∼10 km, ∇H_s is associated with the nature of the underlying current (structure and intensity). The current intensity gradients ($\left|\mathrm{\nabla }\mathit{u}\right|$=$\sqrt{({\partial }_{X}U{)}^{\mathrm{2}}+({\partial }_{Y}U{)}^{\mathrm{2}}}$ with $U=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}$, noted ∇U hereinafter), and more specifically the vorticity of the flow (ζ), induces refraction, resulting in ∇H_s patterns correlated to the vorticity patterns (Figs. 1, 2). Note that both ∇U and ∇H_s are scalars. Assuming that the group speed of waves is much bigger than the intensity of the current velocity and that the dominant balance in the conservation of wave action (Eq. 1) is between wave action advection and refraction, Villas Bôas et al. (2020) proposed a scaling between the root mean square (rms) of the vorticity and ∇H_s (see Eq. 15 of the same reference). We propose writing the scaling as a function of the wave steepness (k〈H_s〉) knowing that $C_{\text{g}}\propto \frac{\mathit{\sigma }}{k}$. It yields the following expression:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{\nabla }{H}_{s_{\text{rms}}}\mathit{\sigma }}{{\text{Slope}}_{\text{KE}}\langle H_{\text{s}}\rangle k}}\propto \mathrm{\nabla }{U}_{\text{rms}},\end{array}$$

where Slope_KE is the spectral slope of the kinetic energy spectrum (here equal to 3 for the fully developed eddy). Equation (6) shows that ∇H_s is a function of surface current gradients, wave steepness (〈H_s〉k), and wave incident frequency (σ). The motivation of this paragraph is to know if, from high-resolution wave-height measurements, the nature and the statistics of the flow can be estimated. Today's surface current measurements from sea-level anomalies can capture eddy with a shape similar to Fig. 1a, c (if their lifetime is sufficiently long according to the revisiting time of altimeters). However, eddies with a more realistic shape (Fig. 1b, d) are poorly captured (see Sect. 5.2 of de Marez et al., 2020b). If waves capture information about the current through their interaction with these currents, one can imagine that current signal can be inverted from wave measurements. This would be relevant for data assimilation in oceanic wave models among others.

Today, filtered altimeter data capture the wave height at a fine scale on the global scale (Dodet et al., 2020). The new spectrometer on board the CFOSAT satellite brings a new view of wave measurements from space through directional wave spectra measurements (Hauser et al., 2020). Combining the frequency–direction measurements of CFOSAT and altimeters and knowing the statistics of the surface current at global scale and so the term Slope_KE in Eq. (6), the rms of the current gradients could be estimated. Inverting the wave signal to retrieve surface current properties is not a new concept. To name a few, Rascle et al. (2014) showed that the images of sea surface roughness from synthetic aperture radars provide clear observations of meso- and submesoscale oceanic features due to the presence of waves. Dugan et al. (2001), thanks to the 3D wave spectrum (wavenumber–frequency), were able to estimate the current speed from the current-induced Doppler shift. Also, Yurovskaya et al. (2019) discussed the possibility of retrieving current from the phase shift spectrum between two successive band measurements provided by the Sentinel-2 satellite. However, all these strategies to infer current gradients are pretty much limited in space.

Thanks to our numerical results, we will test the validity of Eq. (6) in the case of the fully developed eddy. The final aim is to know if the nature of the flow can be estimated by inverting fine-resolution H_s or σ (or k) measurements. For that, we propose plotting Eq. (6) for the mean state of both wave and current fields; i.e, replace $\mathrm{\nabla }{H}_{{\text{s}}_{\text{rms}}}$ and ∇U_rms with ∇H_s and ∇U respectively. The mean gradients of the right- and the left-hand sides of Eq. (6) are shown in Fig. 6. The two fields are plotted for the fully developed eddy case and for incident waves fixed at T_p=7 s. ∇H_s and ∇U have been projected along and perpendicular to the wave peak direction (Fig. 3) respectively.

Figure 6 (a) Surface current gradients (∇U) projected perpendicularly to the peak wave direction vector, i.e. the right-hand side of Eq. (6), and (b) normalized wave height gradients ($\frac{\mathrm{\nabla }{H}_{\text{s}}\mathit{\sigma }}{{\text{Slope}}_{\text{KE}}{H}_{\text{s}}k}$) projected along the peak wave direction vector, i.e. the left-hand side of Eq. (6). The two fields are for the fully developed eddy. Panel (b) shows the instantaneous field for simulation initialized with T_p=7 s waves.

Both terms of Eq. (6) are of the same order of magnitude (Fig. 6). ∇U shows rounded structures for the core of both the mesoscale and submesoscale eddies (Fig. 6a), whereas the normalized ∇H_s field shows more elongated horizontal structures aligned with the initial wave direction (270^∘). From X=0 to X=250 km, the normalized ∇H_s patterns are aligned with the direction of incident waves; downstream from X=250 km, patterns follow the rays trajectories shown in Fig. 5d. Albeit the two fields show a difference in shapes, the two eddy fields match both at the mesoscale (the central eddy) and at smaller scales (submesoscale eddies around the core of the ellipsoidal eddy) from X=0 to X=250 km. ∇U exhibits fronts at the boundary of the central eddy, which is also captured by the normalized ∇H_s field at Y=200 km. Inside the central ellipsoidal eddy (between Y=200 and Y=300 km), ∇U shows a smooth and homogeneous field which is captured in Fig. 6b only between Y=200 and 250 km. Readers can also see discrepancies in the areas between the central eddy and the submesoscales eddies, where sharp ∇H_s is shown for Y>300 km, whereas ∇U is very smooth. Downstream of the eddy, even if ∇U is null (Fig. 6a), normalized ∇H_s is very sharp (Fig. 6b).

The normalized ∇H_s shows similar structures to the surface current gradient in the first half of the domain, X between 0 and 250 km (Fig. 6). It is crucial to note that the current gradients estimated from the wave field variability are estimated without any information on the phase of the surface current features. The inversion of the ∇H_s to infer the underlying surface currents seems to be promising; however both the non-local effect of currents on waves and the initial incidence direction (resulting in a prevalent direction of ∇H_s patterns) show that the phase of current gradient is hardly reproduced in most of the domain.

In Fig. 7 we illustrate the scaling (Eq. 6) for all initializations (T_p=7, 10.3, and 16.6 s). The results presented in the figures are for the total current gradient (grey dots) and for the vorticity component of the flow (coloured dots). A point in Fig. 7 is the rms of a normalized ∇H_s and of a ∇U at fixed distance from the left boundary between X=79 and X=423 km (where current velocity is not null). The normalized $\mathrm{\nabla }{H}_{{\text{s}}_{\text{rms}}}$ and ζ_rms follow the first bisectrix of the plot, unlike the total current gradient (∇U). For the coloured dots, the spread around the first bisectrix is noticeable regardless of the intensity of the current gradient (or distance from the left boundary), with a maximum of spread at X<100 km (dark purple dots in Fig. 7). Villas Bôas et al. (2020) proved that ∇H_s is strongly proportional to the vorticity component of the flow (see their Fig. 12). In the present study, the fully developed eddy is strongly rotational; nevertheless the divergent component of the flow is not negligible ($\mathit{\delta }/f_{\mathrm{0}}\sim \mathrm{0.5}$, with δ the relative divergence of the flow). We wanted to show here the effect of the divergence on the proportionality between ∇H_s and ∇U.

Figure 7 Scatter plot of the normalized root mean square of significant wave height gradients as a function of the root-mean-square surface current gradients. Coloured points are the scatter plot for the vorticity component of the surface current gradients, and grey points show the full surface current gradient (diverging component and rotational component). One point corresponds to the root mean square of the two quantities for a constant X; the value of X is given as a colour scale. 〈H_s〉 is the average value of the significant wave height when simulations reach the stationary state. Panels (a, b, c) are for simulations forced with the fully developed eddy initialized with T_p=7, T_p=10.3, and T_p=16.6 s respectively.

A linear regression is performed between the normalized $\mathrm{\nabla }{H}_{{\text{s}}_{\text{rms}}}$ and ∇U_rms (and ζ_rms) in Fig. 7. For the vorticity field only, the slopes of the regression are equal to 0.72, 0.8, and 0.8 for simulations initialized with T_p=7, T_p=10.3, and T_p=16.6 s respectively. The R² score varies between 0.67 and 0.75 for every initialization. For the full-gradient field (vorticity and divergence of the current), the slopes of the linear regression are equal to 1.13, 1.20, and 1.17 for simulations initialized with T_p=7, T_p=10.3, and T_p=16.6 s with R² < 0 for every initialization. This negative R² score means that the linear regression fits the data very badly. These R² scores confirm the results of Villas Bôas et al. (2020) between X=79 and X=423 km, which states that, in the absence of source terms and for mature incident waves, the variability of the H_s field is mainly driven by the rotational contribution of the current.

Where oceanic eddy destabilizes spontaneously due to horizontal sheared current structures (barotropic instabilities) or vertical buoyancy gradients (baroclinic instabilities, mixed-layer instabilities), the resulting ocean surface shows specific ∇U features. Thanks to wave numerical experiments, we were able to observe ∇H_s structures which are similar to the structures of ∇U and more particularly to the vorticity component of ∇U. The amplitude of the two gradients is comparable. Knowing the wave incident direction and frequency, it seems promising to invert the wave signal to infer the underlying vorticity field and, perhaps, the instabilities that created such vorticity structures (according to the shape and the size of ∇H_s). Optical instruments have shown the robustness to retrieve the amplitude of the wave field and its associated directional spectrum at fine spatial resolution in a very wide swath (Kudryavtsev et al., 2017). The use of such an instrument seems to be a good candidate to capture very small-scale current features by inverting wave characteristics as shown in the fully developed eddy. Also, if the incident wave direction and frequency are known, the same work would be possible with H_s derived from altimeter measurements. Nevertheless there is one drawback, and not least, the non-local effects of current on H_s cause non-zero ∇H_s where current can be null.

Measuring surface currents from space has been a very challenging in past decades (Villas Bôas et al., 2019). Altimetry has proven its robustness to capture surface geostrophic current on a global scale by measuring the along-track sea-level anomaly from multiple altimeter missions. The effective resolution of the current products depends principally on the number of satellites. These resolutions have been calculated and show a mean effective resolution coarser than 200 km at mid-latitudes and coarser than 600 km in the equatorial band (Rio et al., 2014; Ballarotta et al., 2019). Even if mesoscale eddies are observable from space (Chelton et al., 2011), surface dynamics at smaller scales are not captured by present altimeter products. As an example, we can cite the small oceanic features in the fully developed eddy (see Sect. 5.2 of de Marez et al., 2020b). This reality has highlighted the necessity to measure surface currents at finer resolution, triggering the emergence of new satellite missions based on innovative measurement methods (Ardhuin et al., 2018; Morrow et al., 2019; Gommenginger et al., 2019; Wineteer et al., 2020). However, even without new current measurements, the wave measurements, which are available both on a global scale and at fine resolution, could be assimilated to current models to improve their accuracy, specifically for the intensity of the simulated current gradients. Nevertheless, additional works will be required to quantify the non-local effects of currents on H_s.

Those current gradients are crucial for a wide range of applications. To cite one example, at front location, where there is a clear contrast in the sea surface temperature field, strong exchanges between the upper ocean and lower atmosphere occur which affect the dynamics of the atmospheric boundary layer (Frenger et al., 2013).

3.4 Wave steepness and implications for satellite altimetry

Both H_s and $T_{\text{m}\mathrm{0},-\mathrm{1}}$ are strongly modulated by the presence of the large cyclonic eddy, which, consequently, modifies the wave steepness (μ). The more turbulent the eddy, the stronger the inhomogeneity in the H_s and $T_{\text{m}\mathrm{0},-\mathrm{1}}$ fields (Figs. 2, 4). The wave steepness is a key parameter for both recent parametric models of the sea-state bias (SSB) in remote sensing measurements (Badulin, 2014; Badulin et al., 2018) and the wave dynamics (wave growth, wind drag, wave breaking; Rapp and Melville, 1990; Song and Banner, 2002). Here, we quantify, still for an isolated eddy, the effect of wave–current interactions on the change of the wave steepness at the meso- and the submesoscale ranges. The aims of this section are to highlight the importance of the current effects in the modification of the wave steepness and to qualitatively discuss the implications of those modifications for remote sensing measurements.

μ is defined in terms of both H_s and wave period (see Eq. 3 of Badulin et al., 2018).

$$\begin{array}{}\text{(7)}& \mathit{\mu }={\displaystyle \frac{\mathit{\pi }{H}_{\text{s}}}{g{T}^{\mathrm{2}}}}\end{array}$$

Note that the wave steepness (μ) is a dimensionless parameters. We use the mean period ($T_{\text{m}\mathrm{0},-\mathrm{1}}$) to compute μ in Eq. (7).

From Eq. (7), we provide the modulations of wave steepness in both the initial and the fully developed eddies when the wave field reaches its stationary state (Fig. 8). The wave steepness is at a maximum where waves and current are opposed, X∼250, Y∼300 km (Fig. 8). The spatial gradients of μ ($\left|\mathrm{\nabla }\mathit{\mu }\right|=\sqrt{({\partial }_X\mathit{\mu }{)}^{\mathrm{2}}+({\partial }_Y\mathit{\mu }{)}^{\mathrm{2}}}$, noted ∇μ hereinafter) look more local than the ∇H_s (Fig. 2a, d). In the fully developed eddy, we can see very localized ∇μ at the location of submesoscale eddies. In these areas, the steepness can reach 0.75, which is equal to almost 75 % of the maximum steepness spotted in the eddy core. Where waves and current are aligned, the wave steepness is at a minimum and almost equal to 0 for the fully developed eddy, X>250, Y∼200 km. The maximum values of μ do not reach very high value (<1.2). In our simulations, the H_s value of incident waves is equal to 1 m; actually, in the ocean, H_s can be much larger, which would multiply μ, presented in Fig. 8, by a factor equal to the H_s of the incident waves. The reader can refer to Fig. 5 of Badulin et al. (2018) to have an idea of the mean values of μ measured by the Envisat altimeter on the global scale.

The wave steepness can be estimated on the global scale from altimeter data with different methods, physical (Badulin, 2014) or parametric (Gommenginger et al., 2003). The wave steepness is regionally modified by the presence of the current, in particular for high incident waves (higher than our initialization of 1 m). The fully developed eddy induces stronger ∇μ than eddy with a Gaussian shape. The presence of submesoscale eddies leads to the creation of local ∇μ (Fig. 8) at scales similar to the submesoscale eddies. As the fully developed eddy is more realistic than the initial eddy, the simulation presented here would help to better understand the quick change of μ measured by altimeters and better estimate the sea state bias (SSB) in altimeter measurements and provide, perhaps, certain bases for new parametric models of SSB in strong current fields. Even without discussing the contribution of small-scale current gradients, one can see the necessity of taking into account current in the estimation of SSB. Indeed, in the present operational SSB models, the wave field is considered homogeneous at the mesoscale range (Sandwell and Smith, 2005), whereas we see in our simulations that the wave field is strongly modified due to the wave–current interactions at the mesoscale range, in other words, in both the initial and the fully developed eddy. Finally, because the submesoscale currents induce significant changes in the μ field, future work on the effects of submesoscale flows on the SSB is needed.

Figure 8 Wave steepness multiplied by 100 computed from the mean state of the significant wave height and the mean period in the initial eddy (a) and in the fully developed eddy (b) for T_p=7 s incident waves.

3.5 Effects of broader-banded incident spectra and nonlinear wave–wave interactions on wave–current interactions

3.5.1 New model setup

In the previous analysis, the incident waves have been simulated with wave spectra Gaussian in frequency with a frequency spread (σ_f) equal to 0.03 Hz. For timescales much larger than the wave period and assuming that the surface elevation field is a Gaussian process, with negative and positive anomalies around the mean sea level, nonlinear wave–wave interactions lead to a change of the wave energy in the wave field (Hasselmann, 1962). Here we wanted to quantify the effects of nonlinear wave–wave interactions on both ∇H_s and $\mathrm{\nabla }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ in the eddy field. To study the cross-spectral energy flux between frequencies, we activate the nonlinear source term (S_nl). The right-hand side of Eq. (1) is thus not equal to 0 any more but to S_nl. Because simulations initialized with a very narrow-banded spectrum do not show a clear difference between simulations with and without S_nl (not shown), we extend the frequency spread of the initial wave spectra to σ_f=0.1 Hz.

For sufficiently steep waves, nonlinear wave–wave interactions redistribute wave energy between frequencies over the spectrum, which strongly modifies the shape of the spectrum (Komen et al., 1984). As ∇H_s is a function of the wave steepness (kH_s, Fig. 6), we expect that nonlinear wave–wave interactions would have an impact on the intensity of ∇H_s. Nonlinear wave–wave interactions are simulated using the discrete interaction approximation (Hasselmann et al., 1985). The wave simulation is run for a sufficiently long time to capture the long-term effects of nonlinear wave–wave interactions on the wave parameters. Wave simulation is performed only for 7 s incident waves over the fully developed eddy field. This section is a simple introduction of how both wave–wave interactions and wave–current interactions could induce inhomogeneity in the wave field, still in a very idealized framework.

More detailed studies will have to be conducted, such as with the activation of the other source terms. For instance, activating the wind input source term with a given wind field will have an effect in the high-frequency band of the wave spectrum (the development of a wind sea), which, subsequently, will interact with the current field. Also, the presence of both wind and current will modify the wind work at the surface of the ocean. This work is a function of the difference between the wind speed and the surface current speed. This relative wind will modulate the wind growth and therefore the wave height in the current field. It would be interesting to spatially scale the effect of this relative wind on the wave parameters (H_s, $T_{\text{m}\mathrm{0},-\mathrm{1}}$). Also, the wave dissipation source term will constrain the wave energy in the domain, especially in the areas where the wave steepness is very sharp (Fig. 8).

3.5.2 Results

For a given wave parameter (H_s or $T_{\text{m}\mathrm{0},-\mathrm{1}}$), the relative difference is computed between a simulation where the nonlinear source term is activated and deactivated (Eq. 8),

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }X={\displaystyle \frac{X_{S_{\text{nl}}}-X_{\text{no}{S}_{\text{nl}}}}{X_{\text{no}{S}_{\text{nl}}}}}\times \mathrm{100}.\end{array}$$

The nonlinear wave–wave interactions have a large effect on the spatial gradients of wave parameters studied before; H_s is globally enhanced whereas $T_{\text{m}\mathrm{0},-\mathrm{1}}$ is decreased (Fig. 9). These changes are more visible where waves and currents are aligned, X>250 at Y∼200 km. The spatial variability of the H_s can reach +80 % when S_nl is activated. It has been shown that at the same location, wave–current interactions alone showed a strong decrease in H_s (Fig. 2). One can also see that simulation with wave–wave interactions enhances the H_s in the submesoscale eddy field. Globally, we see that, taking into account the nonlinear wave–wave interactions, the H_s values increase where the currents decreased the H_s and vice versa. We cannot quantitatively compare Figs. 9a and 2d (7 s incident waves in the fully developed eddy), because the incident waves have a different spread in frequency.

Nonlinear wave–wave interactions also highlight a change in the $T_{\text{m}\mathrm{0},-\mathrm{1}}$ field. $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ shows the opposite spatial variation in ΔH_s. Indeed, where ΔH_s was (strongly) positive, $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ is (strongly) negative and vice versa. A transect at X=300 km shows the values of H_s and $T_{\text{m}\mathrm{0},-\mathrm{1}}$ along the vertical (Fig. 9c, d). One can see that ∇H_s is globally reduced due to nonlinear wave–wave interactions, especially in the core of the central eddy (Y between 200 and 350 km). At the location of submesoscale eddies, ∇H_s is also sharper for simulation without S_nl. The $\mathrm{\nabla }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ field shows a much more striking difference between simulations with and without nonlinear wave–wave interactions. The transect presented in Fig. 9d shows that $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ is also the most pronounced in the core of the eddy where $T_{\text{m}\mathrm{0},-\mathrm{1}}$ increases by 4 s with respect to the mean period at X=300 km (∼ 8 s) for simulation without source terms. The simulation with S_nl shows an increase in $T_{\text{m}\mathrm{0},-\mathrm{1}}$ values only by 2 s. Whether for H_s or $T_{\text{m}\mathrm{0},-\mathrm{1}}$, in the current field, wave–wave interactions have the tendency to smooth spatial gradients of the wave parameters driven by wave–current interactions. Here the choice of the parametrization of the nonlinear wave–wave interactions was arbitrary (Hasselmann et al., 1985).

It would be interesting to extend this study to other parameterizations of S_nl to better describe how nonlinear wave–wave processes modify regional wave parameter gradients in a strong current field. Also, because the wave–wave interactions modify the intensity of the ∇H_s, it would be interesting to again characterize the proportionality between ∇H_s and the vorticity of the flow (Eq. 6). In this new numerical framework, considering the nonlinear wave–wave interactions, the R² score between the rms ∇H_s and the rms vorticity of the flow drops from 0.67 (simulation without the source term) to 0.42.

This preliminary work on the effects of the source term S_nl on the wave field in a realistic eddy field has shown that wave–wave interactions modify the wave field in a current field with strong current gradients. Those nonlinear interactions led to a significant change of wave parameters in the whole domain, with the tendency to smooth wave parameter gradients. This work could be extended to other source terms such as S_in (describing processes of wave generation due to wind) or S_dis (describing a great number of processes of wave dissipation). As an example, we showed in Sect. 3.4 that the wave steepness (μ) is strongly modified due to the presence of the eddy field. These changes of μ could lead to an increase in the probability of breaking, subsequently leading to a strong dissipation of the wave energy. It would have large consequences on the potential inversion of the wave signal to estimate the statistic of the underlying current (Eq. 6).

Figure 9 Model difference between solutions with and without nonlinear wave–wave interactions. Panels (a) and (b) show the relative difference in percent of the significant wave height and the mean wave period. Panels (c) and (d) show a transect at X=300 km for simulations without (solid blue line) and with (solid red line) a nonlinear source term (S_nl) for H_s and $T_{\text{m}\mathrm{0},-\mathrm{1}}$ respectively.

4 Conclusion and perspectives

In this paper, we numerically studied the effect of an isolated composite cyclonic eddy on the wave properties. Fine-resolution wave simulations have been forced with a composite eddy reconstructed from in situ measurements in the Arabian Sea. The wave model has been forced on the one hand by an initial eddy field with a Gaussian shape and, on the other hand, by a fully developed eddy resulting from the destabilization within the composite eddy. Waves have been simulated by the use of a third-generation phase-averaged spectral model initialized with narrow wave spectra centred at different frequencies (T_p=7, T_p=10.3, and T_p=16.6 s). Although wave refraction by an oceanic vortex has already been studied in former papers (Mapp et al., 1985; White and Fornberg, 1998; Gallet and Young, 2014), this study complements studies performed in the past with (1) a description of the evolution of the wave bulk parameters (significant wave height and mean wave period) inside and outside the isolated vortex and (2) an investigation of how a fully developed eddy (which really occurs in a real ocean) modifies the wave field. Both wave dynamics and kinematics are changed by the presence of underlying currents. These changes are more pronounced where the underlying current has a stronger vorticity. We have shown that the current-induced refraction is stronger for short incident waves and for highly rotational flows. This is consistent with the studies of Kenyon (1971) and Dysthe (2001). As the eddies, dynamical at both the meso- and the submesoscale, are certainly more realistic than Gaussian eddies, former studies of interaction between waves and Gaussian eddy underestimate the current-induced refraction, the intensity of ∇H_s, and the wave steepness inside and in the vicinity of an isolated vortex. Those underestimations can have a large impact on the wave forecast but also on the source of noise induced by waves in the ocean level measurements by altimeters (sea-state bias, among others). Thanks to our simulations, we expect that relationships developed in the past for SSB models cannot be applied in strong current gradients. For instance, Tran et al. (2010) proposed combining altimeter measurements and wave simulations in order to develop a global sea-state bias model. Thanks to the sea-state measurements and period provided by wave model (only forced with wind), authors showed the possibility of significantly reducing the error budget in the SSB estimation. However, Tran et al. (2010) parameterized their wave model on a much too coarse grid ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$) without taking into account current forcing. As we proved here, short-scale currents induce large modifications of wave period at a regional scale (smaller than the wind scales). Indeed, in the current field, even in a very idealized eddy, $\mathrm{\Delta }{T}_{\text{m}\mathrm{0},-\mathrm{1}}$ oscillates within 1 s (Fig. 4a–c) and reaches ∼3 s for a more realistic eddy pattern (Fig. 4a–c). So it strongly affects the geometry of the ocean surface through the wave steepness. Redoing the same work of Tran et al. (2010) at finer resolution with current sufficiently resolved would improve their sea-state bias model at the regional scale.

Following the relationship introduced in Villas Bôas et al. (2020) based on the balance between wave action advection and current-induced refraction, the significant wave height gradients normalized by the incident wave frequency have been described as a function of the surface current gradients. Besides a good coherence in terms of magnitude between the two quantities, the structures of the normalized significant wave height gradients are very sensitive to the underlying surface current. This work was motivated by the idea of inverting wave measurements to infer current properties. We know that measurements of sea level anomaly from space are able to monitor geostrophic surface currents at a global scale with a wavelength resolution of several hundreds of kilometres in ice-free areas (Villas Bôas et al., 2019). The total surface dynamics at finer scales cannot be captured by altimeters, whereas a lot of oceanic processes occur at those scales (from 1–100 km). This paper has shown the possibility of inferring the rms of the vorticity of the eddy field from the inhomogeneity in the wave field, as proposed in Villas Bôas et al. (2020). Inferring vorticity patterns could allow the capture of the small-scale processes (vertical movements, mixing, shear flows, etc.) without measurement of surface currents. Nevertheless, this inversion could not work in the vicinity of a strong ∇U field because the wave field retains the effects of the remote ∇U encountered along their propagation. This results in regional inhomogeneities in the wave field, even at the location where current gradients are very weak. The wave inversion is, at best, only partial. So, the best solution to retrieve the current field at high resolution would be a direct measurement of surface currents from space as proposed in Ardhuin et al. (2018), Gommenginger et al. (2019), and Wineteer et al. (2020).

Finally, because the wave–current coupled system is much too complex, much more than the one proposed here, the assumptions proposed in this paper are hardly satisfied in nature and the potential effect of the non-linear wave–wave interactions probably as well. An example is the assumption that submesoscale currents are stationary during wave propagation.

In the present paper, we studied the effect of a turbulent eddy on wave parameters by assuming the underlying current as barotropic in the first metres of the water column. In reality, both the initial and fully developed eddies are strongly sheared along the vertical, particularly in the first 500 m (see Fig. 2 of de Marez et al., 2020b). It is certain that this vertical shear induces a change in the wave dispersion as described by Kirby and Chen (1989) and therefore, would modify the wave parameters. Also, because the geometry of surface oceanic features is strongly modified due to the presence of waves (Hypolite et al., 2021), another relevant study would be to study the deformation of the eddy field due to the waves.