# Coupled estimation of incoherent internal tides and turbulent motions via statistical modal decomposition

Igor Maingonnat, Gilles Tissot, Noé Lahaye

<figure><img src="image_1.png" /></figure>

## To cite this version:

Igor Maingonnat, Gilles Tissot, Noé Lahaye. Coupled estimation of incoherent internal tides and turbulent motions via statistical modal decomposition. Ocean Science, 2025, 21 (2), pp.807-827. 10.5194/egusphere-2024-1483. hal-04866439v2

HAL Id: hal-04866439

https://inria.hal.science/hal-04866439v2

Submitted on 29 Apr 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coupled estimation of internal tides and turbulent motions via statistical modal decomposition

Igor Maingonnat¹, Gilles Tissot¹*, and Noé Lahaye¹*

¹INRIA Rennes Bretagne Atlantique, IRMAR – UMR CNRS 6625, av. General Leclerc, 35042 Rennes, France

*These authors contributed equally to this work.

**Correspondence:** Igor Maingonnat (igor.maingonnat@inria.fr)

**Abstract.** We present a data-driven modal decomposition method that extracts the part of an incoherent internal tidal wave that correlates with the proper orthogonal decomposition (POD) of a turbulent mesoscale flow. This method exploits the *a priori* knowledge that the incoherent internal tide arises from interactions between an incident wave and the turbulent flow, and exploits the corresponding statistical correlation between the two types of motions. The method is presented and tested in an idealised framework based on the rotating shallow water model, where we provide a physical interpretation for the decomposition method based on theoretical considerations. Using idealised simulations with a plane wave propagating through a zonal turbulent jet, we first propose the use of the modal decomposition method as a data analysis technique to understand how the wave is scattered by the flow. In a second step, we construct an estimation algorithm capable of separating the entangled contributions of the wave and mesoscale motions from a single sea surface height snapshot. This algorithm, which consists of estimating the POD coefficients of the turbulent flow – shared by the wave and jet modes – is particularly suitable for configurations where the jet contribution to the SSH is larger than that of the wave.

**Keywords:** Internal tide scattering, data-driven methods, estimation.

## 1 Introduction

Internal tides (IT) are internal waves generated by the interaction of the barotropic tide with the irregular topography, such as ridges or continental slopes, that propagate mainly at tidal frequencies. They are ubiquitous in the ocean and play a crucial role in vertical mixing and energy transport, especially in the deep ocean (Munk and Wunsch, 1998; Vic et al., 2019). Propagating over large distances, they encounter regions with energetic mesoscale turbulence and lose their fixed phase relationship with the astronomical forcing through non-linear interactions with this turbulence. The resulting incoherent internal tide field (often called «non phase-locked» or «non-stationary» internal tide in the literature), highly unpredictable, complicates, for example, our ability to disentangle internal tides and low-frequency turbulent signals from satellite data (Richman et al., 2012). Due to their large coverage and overall impact on the ocean (Zaron, 2017; Nelson et al., 2019), there is a need to understand how incoherent waves propagate, which includes developing algorithms for estimating and separating the surface signature of internal tides and the mesoscale flow from observational data.
In the physical oceanography community, various studies have examined the impact of a low-frequency turbulent flow on internal waves, using various methodological approaches. Ponte and Klein (2015) demonstrated, using idealized numerical experiments featuring a plane wave propagating through a zonal baroclinic jet, that the increase of energy of the turbulent jet enhances the loss of coherence of the internal tide field. Savva and Vanneste (2018) described, based on kinetic transport theory, how random quasigeostrophic and barotropic flows impact energy exchanges in a scattered plane wave. Dunphy et al. (2017); Kelly et al. (2016); Rainville and Pinkel (2006), among others, studied the propagation of IT using reduced order models derived from a Galerkin projection onto a basis of vertical modes. Ward and Dewar (2010) applied a multiscale analysis in the rotating shallow water model (RSW) to study the scattering of inertia-gravity waves (IGW) by the low-frequency flow. The authors showed that waves that resonantly interact with the jet transfer their energy to waves of equal wavelength and that this transfer is intensified for short-wavelength and strong turbulent flows. A consequence of the non-linear interactions between the IT and the low-frequency flow is a spectral broadening, often referred to as «cusps» (Colosi and Munk, 2006; Zaron, 2022), which is a direct translation of the loss of coherence of the signal. Indeed, the quadratic non-linear term becomes, in time frequency-space, a convolution operator between the wave and the mesoscale fluctuations. The spectrum of coherent IT, which corresponds to a finite Dirac sum, is thus broadened by the mesoscale flow, giving rise to an incoherent broadband spectrum. Spectral broadening is a general characteristic of scattered wave fields, widely studied in other disciplines (e.g. in aeroacoustics, Campos, 1978; Clair and Gabard, 2016).

Recently, several methodological approaches have been developed to address the issue of separation of IT and Balanced Motions (BM), motivated in particular by the launch of the SWOT mission (Fu et al., 2024). Due to the presence of an incoherent component, similar wavelengths and a large time lag between two swaths, simple filtering or harmonic decomposition strategies are often ineffective. Some methods are based on physical approaches, such as in Ponte et al. (2017), who assumed a weak signature of IT on surface density fields and considered potential vorticity dynamics to identify the BM from the observations. In the realm of data-assimilation methods, Le Guillou et al. (2021) developed a coupled iterative approach based on a 4D-Var algorithm for IT and back-and-forth nudging for BM. In the range of data-driven methods, there has been a recent focus on deep-learning approaches to disentangle low and high frequency signals (Wang et al., 2022; Gao et al., 2024). Gao et al. (2024) address in particular the difficulty caused by the long revisit time of altimeters by developing a method based on single SSH snapshots. In closer connection to the methodology that will be presented in this paper, Egbert and Erofeeva (2021) and Tchilibou et al. (2024) proposed methods based on a Proper Orthogonal Decomposition (POD) of the IT sea surface height (SSH) field, computed from a HYCOM-based realistic simulation in the former, and from a filtered daily SWOT swath in the latter.

In this paper, we introduce a data-driven modal decomposition of a wave field scattered by a turbulent mesoscale jet. We focus on methods based on the POD (Berkooz et al., 2003; Lumley, 1967). The POD (often called empirical orthogonal functions – EOF – in the geophysical fluid dynamics community) is a modal decomposition method designed to extract recurrent phenomena from flow data (Long et al., 2021). We propose two different methods, namely the broadband POD (BBPOD) and an algorithm based on the extended POD (EPOD; Boree, 2003). Both methods rely on the physical consideration that the scattering of the wave field, corresponding to the time-evolution of its complex amplitude, is driven by the mesoscale dynamics
and therefore occurs at a similar timescale (Bühler, 2014). With this consideration in mind, the BBPOD algorithm corresponds to a POD of the complex demodulated variables and extracts the most energetic modes of variability of a high-frequency component associated with a broadband spectrum. It shares similarities with the spectral POD algorithm (Welch, 1967; Towne et al., 2018), as will be further discussed in the manuscript. The second method applies the extended POD method to the complex wave amplitude in order to extract spatial modes that are correlated with the POD modes of the jet, thereby providing a decomposition of the flow that takes into account the coupling between the mesoscale dynamics and the wave.

Then, we take advantage of this EPOD-based decomposition method to build an algorithm for estimating the mesoscale and wave fields – including velocities – from a single SSH observation. In particular, we address configurations where the jet component dominates the flow, in which case the weak internal tide signal can still be estimated by the coupling induced by their correlation. The algorithm is tested on idealised simulations of a one-layer RSW model representing such a configuration.

The plan of this paper is as follows. We begin in Sect. 2 by describing the RSW model on which this study is based. Derivation of equations for the ansatz of the complex amplitudes is then given in order to complement the data-driven methods with a physical interpretation. In Sect. 3, the BBPOD and EPOD-based methods are first presented (Sect. 3.1 and 3.2), before describing the algorithm to disentangle observations (Sect. 3.3). Finally, Sect. 4 gathers the numerical results computed from five idealised simulations, including a data analysis of the variability of the wave field (Sect. 4.2), and results of estimates of the mesoscale jet and the internal tide (Sect. 4.3).

## 2 Non-linear interactions between mesoscale turbulence and high frequency waves

We present in this section the (dimensionless) one-layer rotating shallow water model (RSW), which is an adequate idealistic model to analyze non-linear interactions between high-frequency waves and a low-frequency turbulent flow (e.g. Vallis, 2006; Ward and Dewar, 2010). A set of asymptotic equations is then derived for the complex wave amplitude and the evolution of the incoherent part. We shall mention that, since internal tides interacting with mesoscale turbulence are the main target of this study, the waves will be often referred to as internal tide (IT) and the low-frequency turbulent flow as mesoscale flow throughout the paper.

### 2.1 One-layer rotating shallow water model

For this specific analysis, the RSW equations are non-dimensionalized as follows.

$$t^{\ddagger} = f_0^{-1}t; (x^{\ddagger}, y^{\ddagger}) = (lx, ly); \mathbf{v}^{\ddagger} = (u^{\ddagger}, v^{\ddagger}) = (Uu, Uv); h^{\ddagger} = H_o \frac{R_o}{B_u} h,$$

where the superscript .$\ddagger$ refers to the dimensional variables. The parameter $f_0$ is the Coriolis frequency, $v$ is the horizontal velocity, and $h$ is the SSH, with a layer thickness at rest $H_o$. The dimensionless parameters are the Rossby number $R_o = U/f_0l$ and the Burger number $B_u = R_d^2/l^2$, where $R_d = \sqrt{gH_o}/f_0$ is the Rossby deformation radius. The characteristic timescale is taken here as the inertial time $T = f_0^{-1}$, which is well suited for studying IT propagation as it is the lower frequency bound of the internal wave spectrum. The reference length scale $l$ and the reference current $U$ are chosen to be of the typical jet thickness
and velocity, respectively. The Coriolis frequency follows the beta-plane approximation, which writes (in dimensionless form) $f(y) = 1 + \beta y$, where $\beta = R_d/(R_T\sqrt{B_u})$, with $R_T$ the Earth radius. Including the $\beta$-term is necessary to maintain a zonal jet structure (Vallis, 2006). Finally, we consider a monochromatic wave forcing term $f_{q,\omega} = (f_{h,\omega}, f_{v,\omega})^T = \Re(\tilde{f}_{q,\omega}(x,y)e^{i\omega t})$ (which will be specified in Sect. 4.1).

Introducing this non-dimensionalisation into the RSW model leads to the following set of equations, defined on $\mathbb{R}^+ \times \Omega$ where $\Omega \subset \mathbb{R}^2$ is the bounded spatial domain:

$$ \partial_t h + B_u \mathrm{div}(\mathbf{v}) = -R_o [(\mathbf{v} \cdot \nabla)h + h \mathrm{div}(\mathbf{v})] + f_{h,\omega} \quad (1a) $$

$$ \partial_t \mathbf{v} + (1 + \beta y) \mathbf{v}^{\perp} + \nabla h = -R_o (\mathbf{v} \cdot \nabla) \mathbf{v} + \mathbf{f}_{v,\omega} \quad (1b) $$

where $\mathbf{v}^{\perp} = (-v, u)^T$.

## 2.2 Complex wave amplitude ansatz

The total wave field $\mathbf{q}_{\omega} = (u_{\omega}, v_{\omega}, h_{\omega})^T$ is expressed as:

$$ \mathbf{q}_{\omega}(t) = \Re(\tilde{\mathbf{q}}_{\omega}(t)e^{i\omega t}), \quad (2) $$

where $\tilde{\mathbf{q}}_{\omega}(t)$ is the complex wave amplitude vector, and $\omega$ the dominant wave frequency (the same as in the forcing term). It is assumed that the complex wave amplitude is slowly varying (in time), which reflects that the wave is scattered by the low-frequency mesoscale flow, i.e. we assume scale separation (in time) between the mesoscale flow and the IT (see Sect. 2.3). This ansatz is borrowed from high-frequency asymptotic methods, such as ray tracing or WKB methods (c.f. Bühler, 2014), and also corresponds to the dominant term of a scattered wave in standard multiscale analysis (e.g. Ward and Dewar, 2010; Reznik et al., 2001, in the framework of the RSW model).

We further decompose the wave field into a coherent and incoherent component. The coherent part is obtained by applying a time average $\langle \cdot \rangle = \lim_{T \to \infty} \frac{1}{T} \int_0^T \cdot dt$ to the complex wave amplitude:

$$ \mathbf{q}_{\mathrm{coh},\omega}(t) = \Re(\langle \tilde{\mathbf{q}}_{\omega} \rangle e^{i\omega t}), \quad (3) $$

while the incoherent part is the residual:

$$ \mathbf{q}_{\mathrm{incoh},\omega}(t) = \mathbf{q}_{\omega}(t) - \mathbf{q}_{\mathrm{coh},\omega}(t) = \Re(\tilde{\mathbf{q}}'_{\omega}(t)e^{i\omega t}), \quad (4) $$

where $\tilde{\mathbf{q}}'_{\omega}(t) = \tilde{\mathbf{q}}_{\omega}(t) - \langle \tilde{\mathbf{q}}_{\omega} \rangle$ is the incoherent complex amplitude. From the definition (Eq. (3)), we recover the fact that the coherent wave is phase-locked to a monochromatic forcing with frequency $\omega$.

The timescale separation between the mesoscale and the wave allows extracting the complex wave amplitude from the data $\mathbf{q} = (u, v, h)^T$ by a complex demodulation technique:

$$ \tilde{\mathbf{q}}_{\omega} = 2\overline{\mathbf{q}e^{-i\omega t}}, \quad (5) $$

where $\tilde{\mathbf{q}}$ is a low-pass filter (in time) selecting the broadband spectrum of the wave.
## 2.3 Equations for the complex wave amplitude

From the RSW model (Eqs. (1)), we derive idealised equations for the coherent and incoherent complex amplitudes in order to provide a theoretical framework for the interpretation and justification of the extended POD method presented in Sect. 3.2. Yet, a reader more interested in the estimation algorithm could skip this section and go directly to the description of the extended POD technique (Sect. 3.2).

### 2.3.1 Mesoscale / wave separation

In order to simplify forthcoming developments, we assume that the wave field is of small amplitude. This hypothesis allows separating the equations into one for the mesoscale flow and one for the waves (possible relaxation of this hypothesis will be discussed below). The state vector is decomposed as $q = q_{jet} + \epsilon q_{\omega}$, where $q_{jet}$ is the slow mesoscale component and $\epsilon$ is the small parameter of the perturbation expansion. Likewise, the forcing term is also of order $O(\epsilon)$.

At leading order, one obtains the following equation for the mesoscale flow (which is not the focus of the study):

$$ \partial_t h_{jet} + B_u \mathrm{div}(v_{jet}) = -R_o \left[ (v_{jet} \cdot \nabla)h_{jet} + h_{jet} \mathrm{div}(v_{jet}) \right] \quad (6a) $$

$$ \partial_t v_{jet} + (1 + \beta y) v_{jet}^{\perp} + \nabla h_{jet} = -R_o (v_{jet} \cdot \nabla) v_{jet}. \quad (6b) $$

The equation for the wave is obtained at order $\epsilon^1$:

$$ \partial_t h_{\omega} + B_u \mathrm{div}(v_{\omega}) = -R_o \left[ (v_{jet} \cdot \nabla)h_{\omega} + (v_{\omega} \cdot \nabla)h_{jet} + h_{jet} \mathrm{div}(v_{\omega}) + h_{\omega} \mathrm{div}(v_{jet}) \right] + f_{h,\omega} \quad (7a) $$

$$ \partial_t v_{\omega} + (1 + \beta y) v_{\omega}^{\perp} + \nabla h_{\omega} = -R_o \left[ (v_{jet} \cdot \nabla) v_{\omega} + (v_{\omega} \cdot \nabla) v_{jet} \right] + f_{v,\omega}. \quad (7b) $$

For conciseness, we rewrite the above system of equations in terms of linear and bilinear operators:

$$ \partial_t q_{\omega} + L(q_{\omega}) = -R_o B(q_{jet}, q_{\omega}) + f_{q,\omega}, \quad (8) $$

with,

$$ L(q_{\omega}) = \begin{pmatrix} B_u \mathrm{div}(v_{\omega}) \\ (1 + \beta y) v_{\omega}^{\perp} + \nabla h_{\omega} \end{pmatrix}, \quad B(q_{jet}, q_{\omega}) = \begin{pmatrix} (v_{jet} \cdot \nabla)h_{\omega} + (v_{\omega} \cdot \nabla)h_{jet} + h_{jet} \mathrm{div}(v_{\omega}) + h_{\omega} \mathrm{div}(v_{jet}) \\ (v_{jet} \cdot \nabla)v_{\omega} + (v_{\omega} \cdot \nabla)v_{jet} \end{pmatrix}. $$

$L$ is a linear operator, and $B$ is a symmetric bilinear map associated to the non-linear interactions with the jet. Finally, we derive an equation for the wave amplitude by injecting the ansatz Eq. (2) into Eq. (8), then performing a complex demodulation operation Eq. (5), with a filter eliminating dynamics at frequencies higher than typical mesoscale frequencies, (e.g. frequencies at $2\omega$ resulting from complex demodulation). This repeats the hypothesis (Sect. 2.2) that the frequency $\omega$ of IT is much higher than the frequencies of mesoscale flows. We thus obtain:

$$ \partial_t \tilde{q}_{\omega} + i\omega \tilde{q}_{\omega} + L(\tilde{q}_{\omega}) = -R_o B(q_{jet}, \tilde{q}_{\omega}) + \tilde{f}_{q,\omega}. \quad (9) $$
A Reynolds decomposition of the jet component into a mean and a fluctuating part $q_{jet} = \langle q_{jet} \rangle + q'_{jet}$, and of the wave into its coherent and incoherent part $\tilde{q}_{\omega} = \langle \tilde{q}_{\omega} \rangle + q'_{\omega}$, gives the following equation:

$$ \partial_t \tilde{q}'_\omega + i\omega(\langle \tilde{q}_\omega \rangle + \tilde{q}'_\omega) + L(\langle \tilde{q}_\omega \rangle + \tilde{q}'_\omega) + R_o B(\langle q_{jet} \rangle, \langle \tilde{q}_\omega \rangle + \tilde{q}'_\omega) = -R_o B(q'_{jet}, \langle \tilde{q}_\omega \rangle + \tilde{q}'_\omega) + \tilde{f}_{q,\omega}. \quad (10) $$

Taking the time average of Eq. (10) leads to an equation for the coherent wave amplitude:

$$ i\omega\langle\tilde{q}_{\omega}\rangle + L(\langle\tilde{q}_{\omega}\rangle) + R_o B(\langle q_{jet}\rangle, \langle\tilde{q}_{\omega}\rangle) = -R_o\langle B(q'_{jet}, \tilde{q}'_{\omega})\rangle + \tilde{f}_{q,\omega} \quad (11a) $$

$$ \langle \tilde{q}_{\omega} \rangle = R[\tilde{f}_{q,\omega} - R_o \langle B(q'_{jet}, \tilde{q}'_{\omega}) \rangle], \quad (11b) $$

where $R = [i\omega I + L + B(\langle q_{jet} \rangle, \cdot)]^{-1}$ is the resolvent operator, which allows us to write the equation in a more compact form. The resolvent operator can be viewed as the integral operator associated with the Green's function (Cavalieri et al., 2019), and is well defined if $-i\omega$ is not an eigenvalue of $L + B(\langle q_{jet} \rangle, \cdot)$. This equation means that the coherent wave field is the sum of the linear response to the forcing $\tilde{f}_{q,\omega}$, and of the averaged interaction between the jet fluctuations and the incoherent wave field.

Finally, subtracting Eq. (11a) from Eq. (10) leads to the equation for the incoherent wave amplitude:

$$ \partial_t \tilde{q}'_\omega + i\omega \tilde{q}'_\omega + L(\tilde{q}'_\omega) + R_o B(\langle q_{jet} \rangle, \tilde{q}'_\omega) = -R_o B(q'_{jet}, \langle \tilde{q}_\omega \rangle) - R_o B(q'_{jet}, \tilde{q}'_\omega) + R_o \langle B(q'_{jet}, \tilde{q}'_\omega) \rangle \quad (12a) $$

$$ R \partial_t \tilde{q}'_\omega(t) + \tilde{q}'_\omega(t) = \underbrace{R_o R[\langle B(q'_{jet}(t), \tilde{q}'_\omega(t)) \rangle - B(q'_{jet}(t), \tilde{q}'_\omega(t))]}_{\text{Multiple scattering}} - \underbrace{R_o R B(q'_{jet}(t), \langle \tilde{q}_\omega \rangle)}_{\text{single scattering}} \quad (12b) $$

In this equation, the non-linear interactions are decomposed into a single scattering term, which is the interaction between the coherent wave $\langle \tilde{q}_\omega \rangle$ and the jet fluctuations, and a multiple scattering term, which is the interaction with the jet fluctuations and the incoherent wave component.

If we now neglect the slow variations in the complex amplitude, assuming that they evolve over long time scales relative to the wave period $2\pi/\omega$, we obtain the reduced equation:

$$ \tilde{q}'_\omega(t) = \underbrace{R_o R[\langle B(q'_{jet}(t), \tilde{q}'_\omega(t)) \rangle - B(q'_{jet}(t), \tilde{q}'_\omega(t))]}_{\text{Multiple scattering}} - \underbrace{R_o R B(q'_{jet}(t), \langle \tilde{q}_\omega \rangle)}_{\text{single scattering}} \quad (13) $$

It can be noted that this time scale separation can also be interpreted as assuming that the resolvent operator is approximately constant over the frequency band considered for the broadband scattered wave, since the bilinear term is broadband but $R$ is defined at the frequency $\omega$. This approximation consists in considering a jet «frozen» with respect to the propagation time of the wave to travel the domain. It becomes limiting in large domains where the propagation time cannot be neglected compared to typical times of the balanced flow. It is the limit of the use of instantaneous correlations between the jet and the wave in the EPOD method Sect. 3 and is referred to as the local scattering hypothesis.

The above developments can be justified without assuming a small wave amplitude, but asymptotically from the assumption of separation of time scales, using perturbation theory (chapter 4 of Sutherland (2010)). This involves introducing the change
of variable $\tilde{q}_{\omega}(\epsilon t)$ and retaining the dominant term of a multiscale decomposition of the wave. As a consequence, Reynolds stresses remain at order 0, and generalised Reynolds stresses associated with the rise of higher harmonics arise in Eq. (7b) at order 1. A similar separation of equations has been generalized to high amplitude waves by multiple-scale analysis in Dewar and Killworth (1995) in a 1-layer quasi-geostrophic model, in Reznik et al. (2001); Thomas (2016) for RSW models, and in Sect. 2 of Thomas (2023) for a primitive equations model in view of investigating the impact of internal waves on the balanced flow. It can be noticed that this is not a critical issue since we focus on the impact of the jet flow on the wave (and not the inverse) and because the identified modes are extracted by modal decomposition from the non-linear simulation, as it will be presented in Sect. 3.

## 3 Methods

This section details the data-driven methods, derived from POD, that are adapted to decompose a wave scattered by a turbulent flow. These methods are the Broadband POD (BBPOD – Sect. 3.1) and the Extended POD (EPOD– Sect. 3.2). We then propose an algorithm that uses the EPOD decomposition to disentangle the mesoscale flow and the wave field from SSH observations (Sect. 3.3).

### 3.1 The Broadband POD method

The algorithm that we call Broadband POD (BBPOD) consists of performing a POD on a complex demodulated (wave) field in order to capture its most energetic modes of variability from time series data. This algorithm enables capturing the finite- width frequency-band dynamics with a single basis. The proposed algorithm can be related to the spectral proper orthogonal decomposition (Towne et al., 2018; Schmidt and Colonius, 2020), as detailed in Appendix A. It is also similar to a POD on wavelet-transformed (in time) variables, in the sense that the complex demodulation provides, like the wavelet transform, a temporal description of the signal at a given frequency $\omega$. Such a method has been applied in the context of ocean internal waves by Wang et al. (2000); Pairaud and Auclair (2005), who named the method “wavelet EOF” and performed the transformation either along the vertical direction or the time coordinate.

We consider a set of data $q$ containing a broadband peak centered around a frequency $\omega$. The algorithm to build a BBPOD basis at frequency $\omega$ consists of the following steps. We first compute the complex demodulation of $q$ at the frequency $\omega$, $\tilde{q}_{\omega}$ (Eq. (5)), with the choice of an appropriate filter to capture the broadband structure of the data. We assume that the slowly varying wave amplitude is statistically stationary (at least first and second order), which is a requirement for the POD technique (Towne et al., 2018). We next compute the space auto-correlation tensor of the complex amplitudes:

$$C(x, y, x', y') = \langle \tilde{q}_{\omega}(x, y, t) \otimes \tilde{q}_{\omega}^{*}(x', y', t) \rangle, \quad (14)$$

where $\otimes$ is the dyadic product, product of the components $(q_i q_j)_{i,j}$, and the superscript $*$ denotes the transpose-conjugate operation.
205

Before solving the POD problem, we define an inner product representative of the quadratic energy $E$ of the model Eq. (1), encoded with a positive definite matrix $W_E$:

$$ \|q\|^2_{W_E} = (q, W_E q)_{L^2(\Omega)} = \frac{1}{2} \int_{\Omega} (u^2 + v^2) dxdy + \frac{1}{2B_u} \int_{\Omega} h^2 dxdy. \quad (15a) $$

This norm corresponds to the kinetic plus potential energy for a small perturbation (e.g. Vallis, 2006). The BBPOD modes $(\psi_{n,\omega})_n$ are then defined as the solution to the Fredholm equation:

$$ \int_{\Omega} C(x, y, x', y') W_E(x', y') \psi_{n,\omega}(x', y') dx'dy' = \lambda_{n,\omega} \psi_{n,\omega}(x, y), \quad (16) $$

210

with non-negative eigenvalues $\lambda_{n,\omega}$. The BBPOD modes form an orthonormal basis of square integrable functions (in space) respectively to the inner product defined in Eq. (15). The complex demodulated field can be expressed by the decomposition:

$$ \tilde{q}_{\omega}(t, x, y) = \sum_{n=0}^{\infty} a_{n,\omega}(t) \psi_{n,\omega}(x, y), \quad (17) $$

where $a_{n,\omega}(t) = \int_{\Omega} \psi_{n,\omega}^*(x,y) W_E \tilde{q}_{\omega}(t,x,y) dxdy$ is the $n$-th projection coefficient. The modes are also decorrelated from each other and are optimal to express the quadratic mean energy at frequency $\omega$, calculated as:

215

$$ \langle \|\tilde{q}_{\omega}\|^2_{W_E} \rangle = \sum_{n=0}^{\infty} \lambda_{n,\omega}. \quad (18) $$

In practice, the method of snapshots (Schmidt and Colonius, 2020) is considered for numerical implementation. In the literature, POD usually applies to a zero-mean process, and the mean is subtracted beforehand from the data. However, provided the mean field is a solution to the Fredholm equation above (Eq. (16)), the procedure remains unchanged if we consider the total field or only the fluctuations. Since the mean quantities are relevant both to examining the wave scattering and to the estimation algorithm presented in Sect. 4.3, we keep the total field in our application. The reader can refer to the literature for more details on the theoretical backgrounds of POD methods (Towne et al., 2018; Schmidt and Colonius, 2020; Berkooz et al., 2003).

220

## 3.2 The extended POD method

In this section, we present the Extended POD method (Boree, 2003), which consists of extracting the component of a signal that is correlated with a POD mode of a second signal. Here, we apply this method to the context of wave-current interactions by expressing a decomposition of the wave correlated to a POD decomposition of the mesoscale flow. There is indeed a natural correlation between these two components, as the variability of the wave field is driven by the mesoscale fluctuations via the bilinear term $B(q_{jet}, \tilde{q}_{\omega})$ in Eq. (12a).

225

Given a POD decomposition of the mesoscale flow, the $n$-th extended POD mode of the wave is defined by:

230

$$ \chi_n(\tilde{q}_{\omega}, q_{jet}) = \frac{\langle \tilde{q}_{\omega}(t, x, y) a_n(t) \rangle}{\lambda_n}, \quad (19) $$
where $a_n(t) \in \mathbb{R}$ is the projection coefficient of the jet component onto its POD basis at time $t$ (for the inner product defined in Eq. (15)) and $\lambda_n$ is the associated eigenvalue. In order to lighten the notation in the following, we will drop the second argument $q_{jet}$ and denote by $\chi_n(\tilde{q}_\omega)$ the $n$-th EPOD mode. As the POD modes of the jet are decorrelated from one another, the $n$-th EPOD contribution $a_n\chi_n(\tilde{q}_\omega)$ is the part of the wave correlated to the $n$-th POD mode of the jet, but completely decorrelated from the other POD modes of the jet. It provides a decomposition of the wave component $\tilde{q}_{c,\omega}$ that is correlated with the $N$-order truncation POD decomposition of the turbulent jet:

$$ \tilde{q}_{c,\omega}(t,x,y) = \sum_{n=0}^{N} a_n(t) \chi_n(\tilde{q}_\omega)(x,y). \qquad (20) $$

This decomposition filters out all wave contributions decorrelated from the jet on $\Omega$. These may come from outside the domain or from other sources of variability, such as variations in stratification (Zilberman et al., 2011). We will show now that the EPOD modes of the wave and the POD modes of the jet are not only correlated with each other but also dynamically linked through the resolvent operator. They are the response to the non-linear interactions between the wave and the jet through scattering, thus conferring an easier interpretation of this correlation-based technique. Let us consider Eq. (13) for the incoherent complex amplitude in a regime where the multiple scattering terms can be neglected, e.g. if the interaction zone remains of limited extent (Olbers, 1981). Taking the time average of this equation multiplied by a POD coefficient of the fluctuations $a'_n$ associated with mode $\psi'_n = (\psi'_u, \psi'_v, \psi'_h)^T$, one obtains the following equation:

$$ \frac{\langle a'_n \tilde{q}'_\omega \rangle}{\lambda'_n} = -R_o RB \left( \frac{\langle a'_n q'_{jet} \rangle}{\lambda'_n}, \langle \tilde{q}_\omega \rangle \right), $$

that is:

$$ \chi_n(q'_\omega) = -R_o RB(\psi'_n, \langle \tilde{q}_\omega \rangle). \qquad (21) $$

This equation indicates that the EPOD modes of the incoherent complex amplitude are the instantaneous response to the interaction between the coherent wave and the POD modes of the jet fluctuations. We point out that the extended POD modes actually extract multiple scattering interactions from data, but the link with the POD modes of the turbulent fluctuations is then no longer direct. This correspondence has been successfully tested numerically in Maingonnat (2024) on a configuration similar to the one studied in the present paper. A dynamical link between POD and EPOD has also been exploited in the context of wall-bounded turbulent flows (Karban et al., 2022) and turbulent jet flows (Karban et al., 2023).

Finally, we would like to point out that the numerical calculation of jet POD modes and EPOD modes can be carried out in a single operation. This is done by calculating a POD using the snapshot method on the «extended» vector $(q_{jet}, \tilde{q}_\omega)^T$ and considering the weight matrix $\begin{pmatrix} W_E & 0 \\ 0 & 0 \end{pmatrix}$, i.e. without considering any weight on the wave.
## 3.3 Estimation algorithm

We now propose a data-driven method that uses the EPOD formalism to estimate and distinguish the mesoscale and internal tide field from a snapshot of SSH observation. We consider potentially sparse observations: $Y(t) = h_{obs}(t) = 1_{\Omega_{obs}} h(t)$, where $\Omega_{obs}$ is the observation mask.

260

The proposed method, described in Algorithm 1, consists of two distinct stages, a training phase (1) and the estimation per se (2). In the former, the POD modes of the jet $\psi(q_{jet})$ and the extended POD modes of the wave field $\chi_n(\tilde{q}_\omega)$ are computed over an initial – “learning” – time window. As in the previous sections (Sect. 3.1 and 3.2), the training dataset is time resolved such that separation of IT and mesoscale is possible by time-filtering. The estimation consists in extracting the projection coefficients of the mesoscale part, which link mesoscale modes and IT EPOD modes from observations issued from test set snapshots. We shall stress that the test snapshots do not belong to the training window, but are taken at sufficiently far time instants. Furthermore, they can be isolated, in the sense that they are not time resolved. One may notice that the minimisation method is a simple least-squared regression Eq. (23), and no regularisation terms are considered for simplicity reasons. We rely indeed on the «rigid» structure of the problem conferred by the low-rank modal decomposition and the small ratio between the number of parameters and the size of the observation space. A regularisation term could be straightforwardly added.

265

A first advantage of this method is that the EPOD decomposition allows the wave to be estimated in configurations where the jet dominates in amplitude. These configurations can be hard for estimating the wave field due to the low signal-to-noise ratio and also because a strong jet can lead to a very incoherent wave, like in the Gulf Stream region (Zaron, 2017). A second advantage is that the estimation relies on single snapshots, with no assumption about time sampling and/or correlation. Finally, this method also reconstructs – by correlation – the velocity components from observations of SSH only.

270

275

# 4 Numerical results

The first part of the results (Sect. 4.2) concerns the study of the BBPOD and EPOD modes from five idealised numerical simulations. We next focus on the estimates of the IT and BM fields (Sect. 4.3).

280

## 4.1 Numerical configuration

Five numerical simulations of the RSW equations Eqs. (1), featuring a plane wave interacting with a zonal jet, have been performed. The parameters that vary from one simulation to another are the temporal frequency of the incoming wave, its direction of propagation, and the Rossby number of the turbulent jet – see Table 1. The simulations are labelled W1, which is the reference simulation, then W2 to W5. The Burger number is $B_u = 1$ for all simulations (units for space coordinates is therefore the Rossby radius of deformation) and $\beta = 0.05$. The wave amplitude is small such that wave-wave interactions are negligible. For W1-W4, the jet has approximately the same spectrum and the same energy level, and only the impact of different wave parameters is studied for these runs.

285
Algorithm 1 Coupled estimation of IT and mesoscale fields

1. **Training stage:** Compute the POD modes of the jet and the EPOD modes of the wave field over an initial time window, and construct the observation operator as the superposition of the wave and jet components.

(a) Compute the jet POD modes $\psi_n(q_{jet})$

(b) Compute the wave EPOD modes: $\chi_n(\tilde{q}_\omega) = (\chi_n(\tilde{q}_\omega)_v, \chi_n(\tilde{q}_\omega)_h)^T$ (consisting of the velocity component and SSH).

(c) Construct the observation operator $\mathbb{H}$ using N POD/EPOD modes:

$$ \mathbb{H}(q(t)) = (\hat{h}_{jet} + \hat{h}_{\omega}) \mathbf{1}_{\{(x,y) \in \Omega_{obs}\}} \quad (22a) $$

$$ = \Re[\sum_{n=0}^{N} a_n(t) (\psi_n(q_{jet})_h + \chi_n(\tilde{h}_\omega)e^{i\omega t})] \mathbf{1}_{\{(x,y) \in \Omega_{obs}\}}, \quad (22b) $$

Here, the jet component is approximated by its truncated POD decomposition:

$$ \hat{h}_{jet} = \sum_{n=0}^{N} a_n(t) \psi_n(q_{jet})_h, $$

where the hat symbol $\hat{\cdot}$ denotes reconstructed fields. The wave is expressed by the truncated EPOD decomposition of its complex amplitude (see Eq. (20)):

$$ \hat{h}_{\omega} = \Re[\tilde{h}_{c,\omega} e^{i\omega t}] = \Re[\sum_{n=0}^{N} a_n(t) \chi_n(\tilde{h}_{\omega}) e^{i\omega t}]. $$

2. **Estimation stage:** minimize a cost function to estimate the mesoscale and wave fields from single snapshots of SSH at a time t

(a) Solve the minimisation problem for the observation error to obtain the jet POD coefficients $a_n$:

$$ \mathbf{a}^{LSQ}(t) = \min_{(a_0(t), \dots, a_N(t))} \|\mathbb{H}(q(t)) - \mathcal{Y}(t)\|_{L^2(\Omega_{obs})}^2. \quad (23) $$

(b) Reconstruct the estimated mesoscale and wave (via its complex amplitude) fields:

$$ \hat{\tilde{q}}_{\omega}(t) = \sum_{n=0}^{N} a_{n}^{LSQ}(t) \chi_{n}(\tilde{q}_{\omega}), \quad (24a) $$

$$ \hat{\tilde{q}}_{jet}(t) = \sum_{n=0}^{N} a_{n}^{LSQ}(t) \psi_{n}(\tilde{q}_{jet}). \quad (24b) $$

An estimate of the true wave field can be recovered with the fast wave part $e^{i\omega t}$ following Eq. (2).
<table><caption>**Table 1.** Parameters of the different simulations.</caption><thead><tr><th colspan="6">Numerical simulations</th></tr><tr><th>Parameters</th><th>W1</th><th>W2</th><th>W3</th><th>W4</th><th>W5</th></tr></thead><tbody><tr><td>Frequency ω</td><td>2f<sub>0</sub></td><td>3f<sub>0</sub></td><td>2f<sub>0</sub></td><td>2f<sub>0</sub></td><td>2f<sub>0</sub></td></tr><tr><td>Mode number m<sub>x</sub></td><td>0</td><td>0</td><td>1</td><td>-1</td><td>0</td></tr><tr><td>Rossby number R<sub>o</sub></td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.35</td></tr></tbody></table>

The equations have been discretised using a spectral method in space and a Runge-Kutta time scheme, with the open-source code Dedalus (Burns et al., 2020). The domain $\Omega$ is a rectangular domain of size $[0, 20] \times [-20, 20]$, discretised on a $128 \times 256$ grid for W1-W4, and on a $256 \times 1024$ grid for W5, which requires a higher resolution. To simulate open boundaries, «sponge layers» are considered, in which the various fields are damped at both edges of the domain in the $y$ direction (grey area in Fig 2). These regions are a numerical artifice that delimit the «physical» domain where $|y| < 12$. In particular, $\beta y$ is linear in the physical domain, and periodicity of the Coriolis term is ensured in the sponge region through a smooth recovery function, thus allowing the use of a Fourier basis in a non-homogeneous direction. As these regions are not physical, they are not included in the different POD calculations.

All simulations are initialized with an eastward zonal jet at geostrophic equilibrium with a small perturbation superimposed to trigger its destabilisation. An eastward wind forcing term, which is constant in time and follows a Gaussian function in $y$, and a radiative damping term of the form $\alpha h$ (where $\alpha = 10^{-3}$) are added to the system in order to maintain the balanced current in a statistically stationary state (e.g. Brunet and Vautard, 1996). To ensure numerical stability, a small hyperviscosity diffusion term $-R_o/R_e \frac{1}{1+R_o h} \nabla^4$ is also added in the momentum equations, with $R_e = 2 \cdot 10^5$, following (Ochoa et al., 2011). Furthermore, the wave field is generated (and re-absorbed at the North) in the sponge layers by nudging toward the incoming plane wave solution – see the visualisation given in Fig. 2.

The model is first run without wave forcing. Once stationarity is reached (after $4000 f_0^{-1}$, which corresponds approximately to 450 days at mid-latitude), we activate the generation of a plane wave in the sponge layers at the bottom of the domain, which propagates along $y$ before interacting with the jet. The properties of the wave forcing follow the dispersion and polarisation relations computed from the harmonic solutions of Eqs. (1) linearised around a state at rest. This second phase of simulation is $4000 f_0^{-1}$ long for W1-W4 and $8000 f_0^{-1}$ for W5. Snapshots are saved every 1/10th of a wave period, enabling extraction of the wave field. The mesoscale part is then extracted (offline) using a low-pass filter (fourth-order Butterworth) with a cutoff frequency of $2\pi/10 f_0^{-1}$ (corresponding approximately to a period of 3 days at mid-latitude), and the wave field is extracted by complex demodulation at frequency $\omega$ (and using the same low-pass filter).

Figure 1 shows the power spectrum of the SSH field, computed using the Welch method with a Hann window of size $160 f_0^{-1}$ (512 snapshots) with 50% of overlap. It clearly exhibits 2 broadband spectral peaks: one, around $\omega = 0$, associated with the jet and another, around $\omega = 2$, associated with the wave, and highlights the timescale separation between both flow components that we discussed in Sect. 2. We point out that sub-mesoscale contributions are also extracted by complex demodulation, but
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Power spectrum of the SSH field in lin-log scale, in the center of the domain $x = 10, y = 0$, for W1. The orange dotted line indicates the filter cut-off frequency.</figcaption></figure>

these remain negligible. One may also notice that the SSH contribution of the wave is only 1 – 2% of that of the low-frequency turbulent flow, as well as the presence of a weak super-harmonic signal (at $\omega = 4f_0$), which we do not treat in this study.

A snapshot of the vorticity field and the complex amplitude of the SSH of the northward propagating wave for the reference simulation is displayed in Fig. 2. It shows that the incident wave in the south of the domain is almost plane, while more disorganised patterns are visible in the north of the domain, which is the signature of the loss of coherence caused by interaction with the mesoscale flow.

## 4.2 BBPOD and EPOD analysis

In this section, we show and discuss how the BBPOD and EPOD methods perform in extracting the variability of the wave field, and provide some information on the underlying mechanisms. The modes are computed in the physical domain, discarding the sponge regions (for $|y| > 12$).

### 4.2.1 BBPOD analysis of the wave field

Figure 3 shows the SSH component of the first three weighted BBPOD modes ($\sqrt{\lambda_{n,\omega}}\psi_{n,\omega}$) of the wave field for every simulation. The first dominant mode (the most energetic), in the first column, corresponds to harmonic waves propagating in the same direction as the incident wave. It corresponds to the coherent part of the wave, as confirmed by the time series of the corresponding projection, which is nearly constant in time (see Fig. 4b). We should mention that the fact that the coherent wave is captured by the first BBPOD mode is not guaranteed a priori. Here, it occurs because the energy is integrated over the whole domain, including the lower half of the domain where wave propagation is essentially coherent. We can also note a slight damping of the amplitude of the coherent mode for the W2 and W5 simulations, where the wave is subject to greater
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Snapshots of the reference simulation W1 with $\omega = 2f_0, m_x = 0$: (a) total vorticity field and (b) complex amplitude of the wave SSH contribution. The sponge layer is included in this visualisation for $|y| > 12$, with the wave forcing region to the south and the relaxation layer to the north and south.</figcaption></figure>

interactions. This damping can be interpreted as resulting from the wave-mesoscale correction (last right-hand side term in Eq. (11a)).

The sub-optimal modes represent contributions from the incoherent wave, i.e. they are a residue of the coherent part, which is decorrelated from the wave forcing. The first two dominant incoherent modes exhibit nearly plane waves deflected by the jet in the upper part of the domain in several directions (where we can also note weak reflections by the jet in the lower half of the domain).

These $x$-wise structures of these waves are very close to Fourier modes associated with a pair of wavenumbers $m_x = k_x L_x / (2\pi)$ that are, from top to bottom, $(-2, 2)$, $(-2, 2)$, $(3, -1)$, $(-3, 1)$. This is a consequence of the statistical homogeneity of our configuration (periodic domain, zonal jet), i.e. the statistics do not depend on $x$. This statistical property ensures that the BBPOD modes converge to Fourier modes in $x$ as the number of snapshots of the correlation matrix increases (Schmidt and Colonius, 2020).

Calculating the eigenvalues of the correlation matrix gives the average energy contribution of each mode (Eq. (18)). We show (Fig. 4c) the cumulative energy of the incoherent part as a function of the number of modes. This corresponds to the energy contained in the first $N$ incoherent/sub-optimal modes, which is equal to $\sum_{n=1}^{N} \lambda_{n,\omega}$, whereas $\lambda_{0,\omega}$ refers to the energy
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The three dominant weighted BBPOD mode of SSH (left to right), calculated as $\sqrt{\lambda_{n,\omega}}\psi_{n,\omega}$ for the runs W1 to W5 (top to bottom).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** In (a) the three leading POD coefficient of the jet: $a_n = (q_{jet}, \psi_n)_{W_E}$. In (b) the three leading BBPOD coefficient of the wave $a_{n,\omega} = (\tilde{q}_{\omega}, \psi_{n,\omega})_{W_E}$. They are computed from the reference run W1, and the respective mean values and root-mean squared errors are shown. In (c) modal cumulative incoherent energy $\sum_{n=1}^N \lambda_{n,\omega}$ of the different runs.</figcaption></figure>

of the coherent mode. The figure shows first of all, that the W5 simulation represents a more energetic incoherent wave than in the other simulations. This can be explained by the fact that the Rossby number is higher, and therefore the contribution of the non-linear terms between the wave and the jet responsible for the loss of coherence is a priori also higher (Eq. (13)). Secondly, the figure shows that the incoherent waves from the W1-W4 simulations have a comparable energy level, as the jet is approximately the same.

### 4.2.2 EPOD analysis

Figure 5 shows the first three dominant jet POD modes and the associated wave EPOD modes for simulations W1 to W4. Akin to the dominant wave BBPOD mode shown previously, the first POD mode of the jet is approximately the mean field and has nearly constant projection coefficients (Fig. 4b). The two sub-optimal modes are meanders in phase quadrature. The first EPOD mode of the wave, which is correlated with the mean field, corresponds to the coherent part and has a spatial structure very similar to the first BBPOD mode. The second and third EPOD modes represent the part correlated with the meanders of the jet and feature standing wave patterns along $x$, which are also in phase quadrature, suggesting a slow zonal propagation velocity following the jet. This is pronounced in simulations where the wave crosses the jet perpendicularly (W1, W2) with 4 nodes in the domain. For the W5 high Rossby number simulation (shown in Appendix C), the jet's sub-optimal modes also represent meanders, and the wave EPOD modes are standing wave patterns, both with doubled wavelengths compared to runs W1-W4.

The correspondence between the wavelengths of the jet meanders and those of the EPOD modes confirms the analytical link between these structures that we discussed in Sect. 3.2 (see Eq. (21)). The dominant wave EPOD modes appear to be the response to interactions between the coherent part and the jet POD modes. The validity of this equivalence suggests that the multiple interaction term is weak for all simulations, which reflects that the mesoscale fluctuations are confined $y$-wise, such that the scattered wave rapidly leaves this zone, leaving no room for further interaction with the fluctuating jet.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Three leading POD modes of the jet (a) and the associated EPOD modes of the wave field (b) for simulations W1-W4. They are weighted by the square-root of the respective POD eigenvalue of the jet $\sqrt{\lambda_n}$. The first three eigenvalues are $\lambda_0 = 22.13, \lambda_1 = 6.09, \lambda_2 = 5.59$ for W1 and $\lambda_0 = 22.65, \lambda_1 = 6.77, \lambda_2 = 6.55$ for W2.</figcaption></figure>
The temporal evolution of the BBPOD wave modes seems to be related to the EPOD/POD behaviour – see Fig. 4: modes 1 are approximately in phase, and modes 2 approximately in phase opposition. Indeed, the standing waves of EPOD modes correspond to the superposition of the upward and downward deflections extracted by the BBPOD method. More precisely, BBPOD extracts these two deviations separately, since they are decorrelated and associated with opposite sign wave numbers, whereas the EPOD method extracts the superposition of these deviations in a single mode by correlation with the jet. It can be noticed in figure 4b, that there is a slight phase shift from an exact phase opposition, which corresponds to the slow zonal propagation of the standing wave. By summing these two first modes, we can conclude that the contribution of these meandering modes to the incoherent wave represents between 40 – 50% of the incoherent energy for W1-W4 and 30% for W5 (Fig. 4c).

370

375

In conclusion, we have shown that the proposed statistical methods – BBPOD and EPOD – allow us to extract and interpret a wave field scattered by mesoscale turbulence. In particular, they allows us to quantify the energetic contribution of the dominant modes of variability of the wave field in connection with the variability of the mesoscale flow. The EPOD shows that the meanders of the mesoscale jet generate an incoherent wave in the form of a standing wave, made up of deviations in directions determined by the primary interactions.

380

## 4.3 Estimation

The second part of this study is dedicated to the problem of estimation, implementing the algorithm 1 introduced in the previous section. Most of the results are shown for the W1 reference simulation. The time series are divided into a learning window and a test window, which duration (for the W1 run) is ≈ 1300 wave periods – 16 months – and ≈ 200 wave periods (50 snapshots separated by 4 wave periods), respectively. For the sake of simplicity, the estimate of the wave is based on the estimate of its complex amplitude, as the total wave can be reconstructed without losing accuracy (using Eq. (2)).

385

### 4.3.1 Estimation from a full SSH observation

Figure 6 shows an estimate of the SSH component of the wave and mesoscale fields, calculated from a complete SSH observation using 30 POD/EPOD modes. Qualitatively, the estimation is in good agreement with the reference, and the error does not exhibit any particular structure or bias.

390

The accuracy of the estimate (with respect to the number of modes considered) is quantified via the $L^2$ norm of the error, averaged over the entire test window, for the complex wave amplitude: $\langle\|\hat{q}_\omega - \tilde{q}_\omega\|_{L^2}^2\rangle/\langle\|\tilde{q}_\omega\|_{L^2}^2\rangle$, and for the mesoscale part: $\langle\|\hat{q}_{jet} - q_{jet}\|_{L^2}^2\rangle/\langle\|q_{jet}\|_{L^2}^2\rangle$ (Figure 7). The estimation error is compared with the averaged projection error performed when the jet is projected on its POD basis and the wave on its BBPOD basis. Indeed, the projection error measures the ability of the basis learnt from the training window to span the observations from the test window: it indicates the maximum energy that can be captured by $N$ modes (this is an inherent property of the POD) and thereby provides a minimal bound for the estimation error. Also shown is the EPOD projection error, which corresponds to the error between the true complex amplitude field and a decomposition over $N$ extended modes, wherein the coefficients are calculated by projection of the jet on its POD basis.

395
<figure><img src="image_8.png" /><figcaption>Figure 6. Estimate of the SSH contribution of the wave (a) and of the jet (b) for 30 modes from one snapshot, for W1.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Time-averaged $L^2$ norm errors (blue) as a function of mode numbers included, together with the projection error of the corresponding POD decompositions, for the turbulent jet (a) and the wave field (b). For the wave, the EPOD decomposition Eq. (20) is also plotted (dotted line). The results correspond to the W1 run.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

We see in Fig. 7 that the estimation error decreases over the first 30 modes for each of the fields. Using 30 modes, 80% of the mesoscale energy is captured. For the wave (Fig. 7a), the estimate captures around 73% of the total energy and 63% of the incoherent energy (associated with the norm defined in Eq. (15)). The total coherent energy of the test window corresponds to the energy captured by the projection of the wavefield onto the first BBPOD mode, which is 46%. The total incoherent energy corresponds to the residual, and the incoherent estimate corresponds to the estimate of the sub-optimal EPOD modes. The error reduction decreases as the number of modes increases, which is very fast until $N = 5$, and then decreases progressively. This indicates a cost/accuracy trade-off, usual in model reduction, which leaves room for tuning as a function of targeted model performances. Furthermore, we observe that the two error estimates (for the mesoscale and the wave) have a similar shape, with a clear break in slope from the 3rd mode – mode 2 and 3 corresponding to the meandering of the jet and the associated response in the form of a standing wave for the wave (c,f, Sect. 4.2), which reflects that the dynamics of both motions are coupled.

For the mesoscale (Fig. 7a), since the POD basis is optimal and the SSH field is mainly dominated by the turbulent jet, the estimation and projection errors collapse. This suggests that the jet coefficients – including for the velocity component, which are included in the error norm – are perfectly estimated. This shows that the SSH is sufficient to estimate the jet flow, which is consistent with the fact that the jet is almost in geostrophic balance. For the wave estimate (Fig. 7b), one sees a difference between the estimation error and the BBPOD projection error. Already for the first mode (corresponding to the coherent part), there is approximately 8% of the total energy that is missing in the estimate. There are several possible explanations for this: firstly, the wave estimate only contains the part that is correlated with the jet and not the part that is decorrelated. Also missing is the fraction of the wave field correlated with the residual of the jet estimate over the 30 POD modes (which represents 20% of the jet’s energy). Secondly, non-stationarity effects in the flow can partially degrade the estimation of correlations between
<table><caption>**Table 2.** Highest average energy captured over the first 30 modes, for wave and jet estimation, for each W1-5 simulation. For W1, this value can also be seen in Fig. 7, which is reached considering the 30 first modes.</caption><thead><tr><th></th><th>W1</th><th>W2</th><th>W3</th><th>W4</th><th>W5</th></tr></thead><tbody><tr><td>Total wave field</td><td>73%</td><td>64%</td><td>82%</td><td>88%</td><td>75%</td></tr><tr><td>Incoherent part</td><td>63%</td><td>21%</td><td>44%</td><td>63%</td><td>50%</td></tr><tr><td>Jet</td><td>79%</td><td>73%</td><td>83%</td><td>88%</td><td>82%</td></tr></tbody></table>

the wave and the jet. Despite the limitations that have been mentioned, we see that for a total observation of SSH, one can have confidence in a large number of modes to estimate the wave and jet fields.

The wave estimation error presented above is calculated for each simulation W1–W5 in table 2. The numbers shown correspond to the minimum of the error as a function of the number of modes for the first 30 modes (corresponding to the maximum energy captured). With the exception of the W2 simulation, which represents a wave of higher frequency compared to the other runs, the method captures more than 73% of the total wave energy for each of the simulations and between 44 and 63% of the incoherent energy. Moreover, the wave estimate does not diverge as the number of modes increases for these simulations. However, the method seems less effective at capturing a higher frequency wave by correlation in the W2 simulation. The best estimate is reached in the 12th mode and diverges very slightly. Unfortunately, we do not have a clear explanation for this. We can simply point out that at high frequencies ($\omega = 3f_0$ being greater than $M_4 \approx 2.8f_0$ at mid-latitudes), the scattering of waves is more complex and more intense (Ward and Dewar, 2010). The estimate of the jet remains fairly robust, since its amplitude remains high for all the simulations. We can also confirm that, by construction of the method, the better the estimate of the jet, the better the estimate of the total wave.

### 4.3.2 Impact of a partial observation

In the final part of the study, we evaluate the sensitivity of the method to the spatial coverage of observations. Figure 8 shows an example of a subsampled observation, consisting of 4 SSH vertical bands that are approximately homogeneously spaced in the domain and one Rossby radius wide.

Figure 9 shows an estimate over 12 EPOD/POD modes, from the same snapshot as in the full observation case, of the complex amplitude of the wave and the jet. Since fewer modes have been taken into account for the estimation, the wave field is smoother, but the error is still qualitatively quite small. For the jet, the error is more pronounced with higher amplitudes than in the case of total observation. The choice of 12 modes corresponds approximately to the minimum error, as shown below.

We calculate Fig. 10 the jet and wave estimation errors in $L^2$ norm averaged temporally over the whole test window. We evaluate the sensitivity of the error by varying the observation coverage from 4 bands (covering 20% of the domain) to 2 bands (10% of the domain). With 4 vertical bands, the algorithm is still able to estimate the first 3 modes with little difference from total observation. These modes correspond to larger jet structures and therefore seem to be less affected by the degradation of the spatial sampling. From 3 modes, the estimate reaches a plateau for the wave and decreases slightly for the jet. At best,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Sparse SSH snapshot observation from the W1 run, covering 20% of the full domain.</figcaption></figure>

around 60% of the jet’s energy is captured and around 65% for the wave. For 3 vertical bands, the error deteriorates slightly, but the curves remain similar to the previous case. Finally, the case of 2 bands is more pathological and shows that when the observation space is too small and the number of parameters to be minimised is too large (here > 11 modes), the error increases. In addition, as in the case of total observations, we note that the coupling of the projection coefficients between the jet and the wave means that the error curves for the jet and the wave necessarily follow a similar trajectory, suggesting that a good estimate of the jet is required in order to have a good estimate of the wave. The increase of the error could be mitigated by adding regularising terms to the algorithm, e.g. penalizing the amplitude of high modes.

## 4.4 Discussion

The EPOD-based estimation method presented above relies on the correlation between the mesoscale fluctuations and the scattered wave field to estimate the latter in a regime where the SSH is dominated by the mesoscale, which is probably the main advantage of the method. To emphasise this, an attempt is made to estimate and separate the mesoscale and wave dynamics from single SSH snapshots using a simpler approach, which is detailed in Appendix B. This involves considering BBPOD modes (with a specific set of projection coefficients to be determined by optimisation) instead of EPOD modes to express the wave. This approach fails completely in reconstructing the wave field (Fig. B1). Although making these two bases independent could provide a better projection space (which here is optimal in terms of energy content), the estimation fails because *i*) it subjects the minimisation problem to overfitting, *ii*) it discards the statistical links between the two bases supported by their dynamical relation Eq. (21) and *iii*) the wave is poorly observed due to a large relative amplitude difference between the SSH contribution of the jet and the wave. Thus, this comparison highlights the benefit of explicitly accounting for the coupling between the mesoscale and the wave via extended modes.
<figure><img src="image_12.png" /><figcaption>Figure 9. Estimate of the SSH contribution of the wave (a) and of the jet (b) for 12 modes from one sparse snapshot.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Time-averaged $L^2$ norm errors for the jet (a) and the wave (b) for sparse observations. Three types of spatial domain coverage are considered (10% – 15% – 20%).</figcaption></figure>

<figure><img src="image_14.png" /></figure>

Let us now discuss the limitations of the proposed approach. First, the EPOD decomposition of the wave only provides an estimate of the part of the wave that correlates with a truncated POD decomposition of the jet. As discussed above, this contribution is a priori important for an incoherent wave resulting from interaction with a turbulent flow, and if the number of modes included in the decomposition is sufficiently large. Secondly, the decomposition bases should be able to span future observations that have not been learnt from the data set. Some of the errors may therefore be due to errors in projecting the wave and jet observations onto their respective bases.

Finally, it can be expected that this approach would become inefficient in configurations where the signature of the IT becomes dominant, and in particular in regimes with very weak mesoscale flows, as the correlation between the mesoscale fluctuations and the wave field would be less pronounced. In this case, recourse to a BBPOD-based algorithm seems more appropriate, provided some improvements are achieved, such as the use of regularisation techniques in the minimisation algorithm to avoid overfitting. For example, Egbert and Erofeeva (2021) showed that a SPOD basis can be used to estimate part of the incoherent signal off the Amazon, using multiple snapshots. Their algorithm estimates the sea surface height (SSH) contribution of the internal tide by finding the IT projection coefficients that minimise the error with multiple observations.

Constructing an algorithm that takes advantage of each method depending on the dynamical regime is a straightforward development that is left for future work. In any case, more sophisticated minimisation algorithms that better constrain the temporal coefficients need to be implemented in order to maintain accuracy as observations degrade. This may include a regularisation term, e.g. that penalises higher modes amplitudes, depending on the choice of an appropriate norm and an initial guess on the solution. In a sense, the methods proposed by Egbert and Erofeeva (2021) or Tchilibou et al. (2024) are examples of constraints put on the temporal coefficient by using a harmonic function to link observations at different times.
# 5 Conclusions and perspectives

In this study, we have proposed new data-driven statistical decompositions of a (internal) wave scattered by mesoscale turbulence. These decompositions are derived from the POD and describe the slow evolution of the complex amplitude of the wave, driven by non-linear interactions with the mesoscale flow. We first introduced the Broadband POD, which adapts the POD for a scattered wavefield with a broadband spectrum. Then, we proposed a decomposition of the wave field that extracts its fraction correlated with the jet, using the Extended POD (EPOD) method – which is probably the most important aspect of this study. We have highlighted a dynamic link that exists between the extended wave modes and the mesoscale POD modes, which holds under certain wave scattering regimes (dominant and localised primary interactions) and allows for a physical interpretation of the extended modes.

We have demonstrated, using idealised rotating shallow water simulations, the ability of both decompositions to analyse the scattering of a low-amplitude wave by a zonal jet. For different wave and mesoscale fields, the dominant BBPOD modes of the wave result from the interactions between the dominant mesoscale POD modes and the coherent part of the wave. A significant part of the incoherent wave is a standing wave generated by the meanders of the jet, determined by these primary interactions.

In the second part of this study, we addressed the issue of disentanglement of wave and mesoscale flow from a single SSH snapshot. Our method relies on a simple minimization algorithm wherein the SSH observation is expressed as the sum of a mesoscale contribution, decomposed on a POD basis, and a wave contribution correlated with the jet, decomposed on the Extended POD modes. This coupled POD/EPOD estimation algorithm allows us to estimate a total wavefield with its velocities, which is of very low amplitude (especially compared to the mesoscale contribution) from isolated SSH snapshots. This allows us to estimate the incoherent wave field (we captured between 44% and 63% of the incoherent wave energy in the idealised RSW simulations), which can potentially give us access to derived quantities such as incoherent wave energy fluxes.

Although this study was primarily motivated by the issue of estimating internal tides, the EPOD method has a fairly general formalism that could be applied to different types of waves scattered by a time-varying flow, provided we can extract from data the variable envelope of the scattered wave. This includes, for example, scattering of near-inertial waves propagating in the ocean, as studied in Danioux and Vanneste (2016); Young and Jelloul (1997) where the propagation of the slowly varying wave envelope is modelled here, the loss of coherence of barotropic tide in coastal areas, or the scattering of acoustic waves by turbulence (Thomas, 2017). Also, this method could be applied to waves of large amplitude, since they remain correlated with the turbulent flow, which is a straightforward extension of this study.

Finally, an extension of this method to the study of a scattered wave in three dimensions, e.g in Kafiabad et al. (2019), may be envisaged. Although EPOD modes can be defined in 3D, it seems appropriate to apply our method to a set of 2D RSW models, similar to the one in this study, after projecting a 3D model onto a vertical mode basis, which is widely used and relevant for the study of IT (Dunphy et al., 2017; Kelly et al., 2016).

Another important perspective of this study is to adapt these methods to more realistic cases. A few limitations need to be addressed in order to achieve this. First, the use of instantaneous correlation between the jet and the wave is mainly valid for estimating the wave locally. It is also necessary to develop techniques to improve the convergence of statistical estimates,
which is not achieved when the number of snapshots in the training time series is limited. One consequence of this convergence problem is that the proposed decompositions do not generate a space of sufficiently large dimension to estimate wave and jet observations. In particular, these techniques would improve the estimation of observations whose variability has not been sufficiently learnt from the reduced-size time series and that is not represented in the set of POD/EPOD modes. These techniques would also make it possible to extend the estimation of waves and currents to the scale of an ocean. Localisation techniques, which broadly consist of defining localised sub-domains to enable considering localised interactions and artificially increasing the number of samples, could address these issues (Hamill et al., 2001; Farchi and Bocquet, 2019). Another methodological challenge to address more realistic configurations is that the internal tide field contains several nearby frequencies. How to separate each coherent peak in the data and adequately describe their correlation with the turbulent flow remains to be understood. Therefore, we think that further work on idealized cases is required before applying this method to a realistic case.

*Code and data availability.* The subsampled time series of lowpassed filter and complex demodulated output of the RSW simulation W1 are provided at https://gitlab.inria.fr/imaingon/internal-tide-simulation.git, including codes to produce data and diagnostics.

## Appendix A: Equivalence between Broadband POD and SPOD

This appendix presents the algorithmic specificities of the common SPOD algorithm using the Welch method (Welch, 1967) to estimate the CSD, and its relation to the BBPOD algorithm, that is used in this study. For some well-chosen parameters, we show that the two algorithms are equivalent.

The following proof is performed with discrete variables, considering a signal $x_t$ with time spacing $\Delta t$ and $t_k = k\Delta t$.

*Proof.* The complex demodulation writes:

$$ \begin{aligned} \overline{x_t e^{-i\omega t_j}} &= \sum_{k=-m}^{m} b_k x_{j-k} e^{-i\omega t_{j-k}} \\ &= \sum_{k=j-m}^{j+m} b_{j-k} x_k e^{-i\omega t_k}, \end{aligned} \quad (A1) $$

where $(b_i)_{-m \le i \le m}$ are the discrete coefficients of the filter $\bar{\gamma}$.

Besides, the principle of the Welch method is to subdivide $x_t$ into possibly overlapping blocks of size $N$ with an overlap $N_o$. A Fast Fourier Transform is then performed on each windowed block to extract the Fourier component at the tidal frequency, denoted $X_\omega^l$ where $l$ is the block index. So,

$$ X_{\omega}^{l} = \sum_{k=-N/2}^{N/2} x_{k+l(N-N_o)} W_k e^{-i\omega t_k}, \quad (A2) $$

where $W_k$ is a window function defined on $[-N/2, N/2]$.
By changing variable $k' = k + l(N - N_o)$, it follows

$$545 \qquad X_{\omega}^{l} = e^{i\omega t_{l(N-N_o)}} \sum_{k'=l(N-N_o)-N/2}^{l(N-N_o)+N/2} x_{k'} W_{k'-l(N-N_o)} e^{-i\omega t_{k'}}. \qquad (A3)$$

Assuming that the window function is symmetric in the middle of each block (which is verified for most windows used in the literature), i.e. $W_k = W_{-k}$, Eq. (A3) gives

$$550 \qquad X_{\omega}^{l} = e^{i\omega t_{l(N-N_o)}} \sum_{k'=l(N-N_o)-N/2}^{l(N-N_o)+N/2} W_{l(N-N_o)-k'} x_{k'} e^{-i\omega t_{k'}}. \qquad (A4)$$

Finally, by choosing the window function as the filter coefficients, i.e. $W_k = b_k$ and $m = N/2$, relation Eq. (A1) yields:

$$555 \qquad X_{\omega}^{l} = e^{i\omega t_{l(N-N_o)}} \overline{x_t e^{-i\omega t}}_{l(N-N_o)}. \qquad (A5)$$

Consequently, up to a phase, the FFT of a block $l$ of size $N$ with overlap $N_o$ at $\omega$ corresponds to the complex demodulation of the signal at time $l(N - N_o)$. The phase shift cancels when computing the correlation over $N_b$ blocks:

$$560 \qquad \sum_{l=0}^{N_b} X_{\omega}^{l} X_{\omega}^{l*} = \sum_{l=0}^{N_b} \overline{x_t e^{-i\omega t}}_{l(N-N_o)} \overline{x_t e^{-i\omega t}}_{l(N-N_o)}. \qquad (A6)$$

Therefore, the Welch method computed with parameters $(N, N_o, N_b, W)$ is equivalent to computing the complex demodulation of the time series over $N_b$ snapshots sampled every $N - N_o$ and a filter chosen as the common window function $W$ (Hann, Hanning,...). $\square$

This simple proof highlights the fact that the window function in the Welch method is playing the role of the filter in the BBPOD method. Although the two methods are equivalent, they differ in their use and interpretation and also in the choice of parameters. BBPOD is designed to extract a broadband signal, the width of which is chosen via a filter using physical time scales. SPOD looks for an orthonormal basis for the Fourier coefficients at a given frequency, and the parameters of Welch's algorithm are chosen to minimise spectral leakage for a given amount of data. Techniques for reconstructing a broadband spectrum using several SPODs at each frequency of the spectrum are possible by reordering their energy contribution (Nekkanti and Schmidt 2021). On the other hand, the orthogonality property of the modes is lost, since the SPOD modes are no longer orthogonal at different frequencies.

## 565 Appendix B: Wave estimation using BBPOD instead of EPOD

A test has been performed in order to show the importance of informing the correlation between the jet contribution and the wave through EPOD in the estimation problem. Here, a «naive» algorithm is implemented, wherein the wave is decomposed on its optimal basis BBPOD (instead of EPOD) and the jet on its POD basis. Instead of considering the same coefficient $a_n(t)$
<figure><img src="image_15.png" /><figcaption>**Figure B1.** Estimation of the SSH contribution of the wave for the W1 run, computed from one snapshot, using 30 modes of the optimal basis of the jet and the wave vectors (POD for the jet, BBPOD for the wave).</figcaption></figure>

for the jet POD mode and the wave EPOD mode, two independent sets of coefficients – $a_n$ for the mesoscale and $b_n$ for the wave – are sought through the optimisation algorithm, which is given by:

$$ \min_{(a_n, b_n)} \left\| \sum_{n=0}^{30} a_n \psi_n(\mathbf{q}_{jet})_h + \Re[b_n(\psi_{n,\omega})_h e^{i\omega t}] - h \right\|^2. $$

As visible in Figure B1, to be compared with Fig. 6, the corresponding estimate is completely wrong. Possible reasons for this failure are provided in the body text.

## Appendix C: Effect of Rossby number on the modal decompositions

In this section, we show the POD modes of the jet (Fig. C1a) and the EPOD modes of the wave (Fig. C1b) for the W5 simulation. This simulation is characterised by a higher Rossby number. The consequence is a change in the jet dynamics and then in the scattered wave fields. This effect is visible in the results of the modal decompositions. The jet POD and wave EPOD modes show similar structures to W1-W4, but the wavelengths of modes 1 and 2 are twice as long.

*Author contributions.* IM wrote the article, performed the numerical tests, and was implicated in the methodology of this study. GT and NL wrote and reviewed the article and both supervised the global validation of the results. They provided expertise on each section of the paper and contributed to the formal analysis and to the global methodology of this article.
<figure><img src="image_16.png" /><figcaption>**Figure C1.** Three leading POD modes of the jet (a) and the associated EPOD modes of the wave field (b) for simulation W5. They are weighted by the square-root of the respective POD eigenvalue $\sqrt{\lambda_n}$. The first three eigenvalues are $\lambda_0 = 8.86, \lambda_1 = 3.29, \lambda_2 = 3.21$.</figcaption></figure>

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* N. Lahaye had support from the French research funding agency under the ModITO project (ANR-22-CE01-0006-01) and from the TOSCA-ROSES SWOT project DIEGO. G. Tissot and N. Lahaye acknowledge support from the French National program LEFE (Les Enveloppes Fluides et l'Environnement).
References

Berkooz, G., Holmes, P., and Lumley, J.: The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows, Annual Review of Fluid Mechanics, 25, 539–575, https://doi.org/10.1146/annurev.fl.25.010193.002543, 2003.

Boree, J.: Extended proper orthogonal decomposition: A tool to analyse correlated events in turbulent flows, Experiments in Fluids, 35, 188–192, https://doi.org/10.1007/s00348-003-0656-3, 2003.

Brunet, G. and Vautard, R.: Empirical normal modes versus empirical orthogonal functions for statistical prediction, Journal of Atmospheric Sciences, 53, 3468–3489, https://doi.org/10.1175/1520-0469(1996)053<3468:ENMVEO>2.0.CO;2, 1996.

Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., and Brown, B. P.: Dedalus: A flexible framework for numerical simulations with spectral methods, Physical Review Research, 2, 023 068, https://doi.org/10.1103/PhysRevResearch.2.023068, 2020.

Bühler, O.: Waves and Mean Flows, Cambridge Monographs on Mechanics, Cambridge University Press, 2 edn., https://doi.org/10.1017/CBO9781107478701, 2014.

Campos, L. M. B. C.: The spectral broadening of sound by turbulent shear layers. Part 1. The transmission of sound through turbulent shear layers, Journal of Fluid Mechanics, 89, 723–749, https://doi.org/10.1017/S0022112078002827, 1978.

Cavalieri, A. V. G., Jordan, P., and Lesshafft, L.: Wave-packet models for jet dynamics and sound radiation, Applied Mechanics Reviews, 71, 2019.

Clair, V. and Gabard, G.: Numerical Investigation on the Spectral Broadening of Acoustic Waves by a Turbulent Layer, chap. 2, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.2016-2701, 2016.

Colosi, J. A. and Munk, W.: Tales of the Venerable Honolulu Tide Gauge, Journal of Physical Oceanography, 36, 967–996, https://doi.org/10.1175/JPO2876.1, 2006.

Danioux, E. and Vanneste, J.: Near-inertial-wave scattering by random flows, Physical Review Fluids, 1, 033 701, https://doi.org/10.1103/PhysRevFluids.1.033701, 2016.

Dewar, W. K. and Killworth, P. D.: Do fast gravity waves interact with geostrophic motions?, Deep Sea Research Part I: Oceanographic Research Papers, 42, 1063–1081, https://doi.org/10.1016/0967-0637(95)00040-D, 1995.

Dunphy, M., Ponte, A., Klein, P., and Le Gentil, S.: Low-mode internal tide propagation in a turbulent eddy field, Journal Of Physical Oceanography, 47, 649–665, https://doi.org/10.1175/JPO-D-16-0099.1, 2017.

Egbert, G. D. and Erofeeva, S. Y.: An Approach to Empirical Mapping of Incoherent Internal Tides With Altimetry Data, Geophysical Research Letters, 48, e2021GL095 863, https://doi.org/10.1029/2021GL095863, 2021.

Farchi, A. and Bocquet, M.: On the efficiency of covariance localisation of the ensemble Kalman filter using augmented ensembles, Frontiers in Applied Mathematics and Statistics, 5, 3, 2019.

Fu, L.-L., Pavelsky, T., Cretaux, J.-F., Morrow, R., Farrar, J. T., Vaze, P., Sengenes, P., Vinogradova-Shiffer, N., Sylvestre-Baron, A., Picot, N., and Dibarboure, G.: The Surface Water and Ocean Topography Mission: A Breakthrough in Radar Remote Sensing of the Ocean and Land Surface Water, Geophysical Research Letters, 51, e2023GL107 652, https://doi.org/10.1029/2023GL107652, 2024.

Gao, Z., Chapron, B., Ma, C., Fablet, R., Febvre, Q., Zhao, W., and Chen, G.: A Deep Learning Approach to Extract Balanced Motions From Sea Surface Height Snapshot, Geophysical Research Letters, 51, e2023GL106 623, https://doi.org/10.1029/2023GL106623, 2024.

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Monthly Weather Review, 129, 2776–2790, 2001.
Kafiabad, H. A., Savva, M. A. C., and Vanneste, J.: Diffusion of inertia-gravity waves by geostrophic turbulence, Journal of Fluid Mechanics, 869, R7, https://doi.org/10.1017/jfm.2019.300, 2019.

Karban, U., Martini, E., Cavalieri, A. V. G., Lesshafft, L., and Jordan, P.: Self-similar mechanisms in wall turbulence studied using resolvent analysis, Journal of Fluid Mechanics, 939, A36, https://doi.org/10.1017/jfm.2022.225, 2022.

Karban, U., Bugeat, B., Towne, A., Lesshafft, L., Agarwal, A., and Jordan, P.: An empirical model of noise sources in subsonic jets, Journal of Fluid Mechanics, 965, A18, https://doi.org/10.1017/jfm.2023.376, 2023.

Kelly, S. M., Lermusiaux, P. F. J., Duda, T. F., and Haley, P. J.: A Coupled-Mode Shallow-Water Model for Tidal Analysis: Internal Tide Reflection and Refraction by the Gulf Stream, Journal of Physical Oceanography, 46, 3661–3679, https://doi.org/10.1175/JPO-D-16-0018.1, 2016.

Le Guillou, F., Lahaye, N., Ubelmann, C., Metref, S., Cosme, E., Ponte, A., Le Sommer, J., Blayo, E., and Vidard, A.: Joint Estimation of Balanced Motions and Internal Tides From Future Wide-Swath Altimetry, Journal of Advances in Modeling Earth Systems, 13, e2021MS002613, https://doi.org/10.1029/2021MS002613, 2021.

Long, Y., Zhu, X.-H., Guo, X., Ji, F., and Li, Z.: Variations of the Kuroshio in the Luzon Strait Revealed by EOF Analysis of Repeated XBT Data and Sea-Level Anomalies, Journal of Geophysical Research: Oceans, 126, e2020JC016849, https://doi.org/10.1029/2020JC016849, 2021.

Lumley, J.: The Structure of Inhomogeneous Turbulent Flows, Atmospheric Turbulence and Radio Wave Propagation, pp. 166–177, 1967.

Maingonnat, I.: Compréhension et modélisation de mécanismes non-linéaires dans l'océan : les interactions entre ondes internes et écoulement, Ph.D. thesis, Université de Rennes, 2024.

Munk, W. and Wunsch, C.: Abyssal recipes II: energetics of tidal and wind mixing, Deep Sea Research Part I: Oceanographic Research Papers, 45, 1977–2010, https://doi.org/10.1016/S0967-0637(98)00070-3, 1998.

Nelson, A. D., Arbic, B. K., Zaron, E. D., Savage, A. C., Richman, J. G., Buijsman, M. C., and Shriver, J. F.: Toward Realistic Non-stationarity of Semidiurnal Baroclinic Tides in a Hydrodynamic Model, Journal of Geophysical Research: Oceans, 124, 6632–6642, https://doi.org/10.1029/2018JC014737, 2019.

Ochoa, J., Sheinenbaum, J., and Jiménez, A.: Lateral Friction in Reduced-Gravity Models: Parameterizations Consistent with Energy Dissipation and Conservation of Angular Momentum, Journal of Physical Oceanography, pp. 1894 – 1901, https://doi.org/10.1175/2011JPO4599.1, 2011.

Olbers, D. J.: A Formal Theory of Internal Wave Scattering with Applications to Ocean Fronts, Journal of Physical Oceanography, 11, 1078–1099, https://doi.org/10.1175/1520-0485(1981)011<1078:AFTOIW>2.0.CO;2, 1981.

Pairaud, I. and Auclair, F.: Combined wavelet and principal component analysis (WEOF) of a scale-oriented model of coastal ocean gravity waves, Dynamics of Atmospheres and Oceans, 40, 254–282, https://doi.org/10.1016/j.dynatmoce.2005.06.001, 2005.

Ponte, A. L. and Klein, P.: Incoherent signature of internal tides on sea level in idealized numerical simulations, Geophysical Research Letters, 42, 1520–1526, https://doi.org/10.1002/2014GL062583, 2015.

Ponte, A. L., Klein, P., Dunphy, M., and Le Gentil, S.: Low-mode internal tides and balanced dynamics disentanglement in altimetric observations: Synergy with surface density observations, Journal of Geophysical Research: Oceans, 122, 2143–2155, https://doi.org/10.1002/2016JC012214, 2017.

Rainville, L. and Pinkel, R.: Propagation of Low-Mode Internal Waves through the Ocean, Journal of Physical Oceanography, 36, 1220–1236, https://doi.org/10.1175/JPO2889.1, 2006.
660

Reznik, G. M., Zeitlin, V., and Ben Jelloul, M.: Nonlinear theory of geostrophic adjustment. Part 1. Rotating shallow-water model, Journal of Fluid Mechanics, 445, 93–120, https://doi.org/10.1017/S002211200100550X, 2001.

665

Richman, J. G., Arbic, B. K., Shriver, J. F., Metzger, E. J., and Wallcraft, A. J.: Inferring dynamics from the wavenumber spectra of an eddying global ocean model with embedded tides, Journal of Geophysical Research: Oceans, 117, https://doi.org/10.1029/2012JC008364, 2012.

Savva, M. and Vanneste, J.: Scattering of internal tides by barotropic quasigeostrophic flows, Journal of Fluid Mechanics, 856, 504–530, https://doi.org/10.1017/jfm.2018.694, 2018.

Schmidt, O. T. and Colonius, T.: Guide to Spectral Proper Orthogonal Decomposition, AIAA Journal, 58, 1023–1033, https://doi.org/10.2514/1.J058809, 2020.

Sutherland, B. R.: Internal Gravity Waves, Cambridge University Press, 2010.

670

Tchilibou, M., Carrere, L., Lyard, F., Ubelmann, C., Dibarboure, G., Zaron, E. D., and Arbic, B. K.: Internal tides off the Amazon shelf in the western tropical Atlantic: Analysis of SWOT Cal/Val Mission Data, EGUSphere, 2024, 1–23, https://doi.org/10.5194/egusphere-2024-1857, 2024.

Thomas, J.: Resonant fast-slow interactions and breakdown of quasi-geostrophy in rotating shallow water, Journal of Fluid Mechanics, 788, 492–520, https://doi.org/10.1017/jfm.2015.706, 2016.

675

Thomas, J.: New model for acoustic waves propagating through a vortical flow, Journal of Fluid Mechanics, 823, 658–674, https://doi.org/10.1017/jfm.2017.339, 2017.

Thomas, J.: Turbulent wave-balance exchanges in the ocean, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 479, 20220565, https://doi.org/10.1098/rspa.2022.0565, 2023.

Towne, A., Schmidt, O. T., and Colonius, T.: Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis, Journal of Fluid Mechanics, 847, 821–867, 2018.

680

Vallis, G. K.: Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation, Cambridge University Press, https://doi.org/10.1017/CBO9780511790447, 2006.

Vic, C., Naveira Garabato, A. C., Green, J. A. M., Waterhouse, A. F., Zhao, Z., Melet, A., de Lavergne, C., Buijsman, M. C., and Stephenson, G. R.: Deep-ocean mixing driven by small-scale internal tides, Nature Communications, 10, 2099, https://doi.org/10.1038/s41467-019-10149-5, 2019.

685

Wang, H., Grisouard, N., Salehipour, H., Nuz, A., Poon, M., and Ponte, A. L.: A Deep Learning Approach to Extract Internal Tides Scattered by Geostrophic Turbulence, Geophysical Research Letters, 49, https://doi.org/10.1029/2022GL099400, 2022.

Wang, J., Chern, C.-S., and Liu, A. K.: The Wavelet Empirical Orthogonal Function and Its Application Analysis of Internal Tides, Journal of Atmospheric and Oceanic Technology, 17, 1403–1420, https://doi.org/10.1175/1520-0426(2000)017<1403:TWEOFA>2.0.CO;2, 2000.

690

Ward, M. L. and Dewar, W. K.: Scattering of gravity waves by potential vorticity in a shallow-water fluid, Journal of Fluid Mechanics, 663, 478–506, https://doi.org/10.1017/S0022112010003721, 2010.

Welch, P.: The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, IEEE Transactions on Audio and Electroacoustics, 15, 70–73, https://doi.org/10.1109/TAU.1967.1161901, 1967.

Young, W. R. and Jelloul, M. B.: Propagation of near-inertial oscillations through a geostrophic flow, Journal of Marine Research, 55, 735–766, https://doi.org/10.1357/0022240973224283, 1997.

695

Zaron, E. D.: Mapping the nonstationary internal tide with satellite altimetry, Journal of Geophysical Research: Oceans, 122, 539–554, https://doi.org/10.1002/2016JC012487, 2017.
Zaron, E. D.: Baroclinic Tidal Cusps from Satellite Altimetry, Journal of Physical Oceanography, 52, 3123–3137, https://doi.org/10.1175/JPO-D-21-0155.1, 2022.

Zilberman, N. V., Merrifield, M. A., Carter, G. S., Luther, D. S., Levine, M. D., and Boyd, T. J.: Incoherent Nature of M2 Internal Tides at the Hawaiian Ridge, Journal of Physical Oceanography, 41, 2021–2036, https://doi.org/10.1175/JPO-D-10-05009.1, 2011.

700