# Localised proper orthogonal decomposition

Dyhia Elhaddad, Igor Maingonnat, Gilles Tissot, Noé Lahaye

► To cite this version:

Dyhia Elhaddad, Igor Maingonnat, Gilles Tissot, Noé Lahaye. Localised proper orthogonal decomposition. Theoretical and Computational Fluid Dynamics, 2025, 39 (5), pp.40. 10.1007/s00162-025-00760-2. hal-05268871

HAL Id: hal-05268871

https://inria.hal.science/hal-05268871v1

Submitted on 19 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Localised proper orthogonal decomposition

Dyhia Elhaddad, Igor Maingonnat,
Gilles Tissot*, Noé Lahaye

**Abstract** Proper Orthogonal Decomposition (POD) extracts an optimal orthonormal basis for reconstructing a sequence of data. This technique is widely used in fluid mechanics to construct reduced-order models for physical analysis, modelling, estimation, and control. However, when the number of snapshots is small and the spatial domain is large, convergence of the basis with respect to the number of samples is compromised. This often results in a degradation of the quality of the resulting reduced-order models. In this study, we propose a localisation procedure to mitigate convergence artefacts in the context of limited data. Specifically, we adapt Schur product-based localisation – commonly used in operational data assimilation – to the POD framework. We develop an algorithm that operates directly on the snapshots rather than on the large-scale correlation tensor. In addition, we introduce a family of partition-of-unity localisation functions to preserve the reconstruction properties of the localised POD basis. By artificially increasing the rank of the correlation tensor, we demonstrate that the localised POD, computed from a limited number of snapshots, can approximate the POD basis obtained from a larger dataset. This is illustrated using a wave field scattered by a turbulent jet in a rotating shallow water system. We further show that using smooth, overlapping localisation functions enables physically meaningful reconstructions. The reconstruction capabilities of the method are also demonstrated using a realistic internal wave field in the Gulf Stream region, based on a high-resolution numerical simulation.

## 1 Introduction

Proper orthogonal decomposition (POD) is the most popular basis construction technique in fluid flow reduced-order modelling. From a set of snapshots, it extracts an orthonormal basis, optimal to reconstruct the data with respect to a given energy norm. This technique is known under several names, POD in turbulence (Lumley, 1967; Berkooz et al., 1993), empirical orthogonal functions (EOF) in geophysical fluid dynamics (Lorenz, 1956), Karhunen-Loève (KL) expansion in probability (Karhunen,

<footnote>Centre Inria de l'Université de Rennes, IRMAR - UMR CNRS 6625, Odyssey team, 263 av. Général Leclerc, 35042, Rennes, France</footnote>

<footnote>* coresp: Gilles.Tissot@inria.fr</footnote>
1946; Loève, M., 1955), or principal component analysis (PCA) in statistics (Jolliffe, 2005). It can be employed to construct Galerkin reduced-order models (Aubry et al., 1988), solve optimisation, estimation (Cordier et al., 2013) and control (Bergmann et al., 2005) problems. Sustained by the development of intensive experimental acquisition data, such as particle image velocimetry (PIV), and of high fidelity numerical simulations, the use of snapshot POD algorithm (Cordier and Bergmann, 2002) has become the most employed, since the number of data points along the space coordinates is most often significantly larger than the number of time samples.

The recent increase of numerical capabilities in direct and large eddy simulations allows to consider complex configurations, with very large numerical domains compared to the physical processes at play. When huge snapshots are considered, their number can become limited due to storage constraints or due to the difficulty to perform very long runs. In this context, POD modes can be difficult to converge, in the sense that when the basis is converged, increasing the number of samples should not change much the alignment of the modes, at least the first ones of interest. As a result, techniques of estimation or control can become less effective. Convergence issues have for instance been encountered in a spectral version of POD defined in the Fourier domain (Blanco et al., 2022). An objective of the present paper is to address this issue through localisation techniques. As an example of possible strategy, in the context of ocean dynamics reconstruction, Egbert and Erofeeva (2021) have performed a PCA of internal tidal waves using overlapping spatial patches, which can be viewed as a form of localisation. This works well, but the treatment of the overlapping regions is not obvious. In particular, a straightforward procedure is to perform independent reconstructions on overlapping patches and to interpolate the result to ensure smoothness of the reconstruction. However, this breaks the orthogonality of the overlapping modes and thereby the optimality of the reconstruction. Patch localisation have also been proposed on domains defined using clustering techniques (Sahyoun and Djouadi, 2013). We will propose in the present paper a more general framework, allowing in particular, even if it is localised, to define the problem globally.

In data assimilation, such as in weather forecast, ensemble methods are employed to estimate at best the state of a system knowing some sparse and noisy observations and an ensemble of realisations of the model. A key ingredient to perform good estimations is to provide background and observation covariance matrices, informing the degree of confidence to the model and measurements, respectively. In ensemble methods (Evensen, 1994), background covariance matrices are approximated by an empirical covariance provided by the ensemble members, whose number is constrained to be relatively small ($O(10) - O(100)$) due to limited computational resources. This results in spurious correlations between far space locations, induced by the rank deficiency of the empirical covariance matrix. To circumvent this issue, localisation (Hamill et al., 2001) consists in performing a Schur product (or Hadamard product, i.e. element by element) between the covariance matrix and a localisation matrix, to enforce these spurious correlations to vanish. The direct consequence is an artificial increase of the rank of the covariance matrix, leading to better estimation performances. The main idea of the present paper is to adapt this technique to improve the POD in the context of limited number of snapshots compared to the space dimension.

In ensemble methods, localisation is implemented based on a low-rank square-root approximation of the localisation matrix, as through a so-called $\alpha$ control variable (Lee et al., 2024), or through augmented ensembles (Farchi and Bocquet, 2019) as employed in operational conditions (Arbogast et al., 2017). The strategy proposed in the present
paper follows these lines. In order to obtain a simple algorithm we further propose to define directly a set of localisation functions of the state variables, associated with a low rank approximation of the localisation matrix.

In the context of reduced-order modelling, localisation strategies have been explored in several studies. Anderson et al. (2023) introduced a method based on non-overlapping local bases combined with sparsity-promoting techniques. Smetana and Taddei (2023) employed a decomposition inspired by the partition of unity (PU) discontinuous Galerkin method, in which basis functions associated with a given subdomain are allowed to be used outside their local region – unlike in our approach. Gastelen et al. (2024) proposed a PU-based strategy for defining kernels of modified inner products. While this approach is close to ours, it differs by resulting in a loss of orthogonality of the localised modes. Localisation has also been applied in the context of proper generalised decomposition by Ibáñez et al. (2019). In another line of work, local principal component analysis (PCA) on subdomains identified via clustering techniques has been investigated by Muñoz et al. (2023). In their approach, overlapping is implicit: at each spatial location, the best reconstruction among the clusters is retained. More recently, Eichberger et al. (2024) proposed a localisation strategy to enhance resolvent analysis (McKeon and Sharma, 2010). Given the strong connection between resolvent analysis and spectral POD (Towne et al., 2018; Cavalieri et al., 2019), the localised POD method introduced in the present work – applicable in both spectral and non-spectral contexts – can be viewed as a data-driven counterpart to the localised resolvent approach of Eichberger et al. (2024).

In the present paper, we propose an algorithm that defines a set of localisation functions acting directly on the snapshots, thereby avoiding the need to compute a Schur product on the full-scale covariance tensor, which is often too large to be stored. We establish the equivalence between this snapshot-based approach and the classical Schur-product-based localisation. In addition, we derive the constraint that this family of functions must satisfy to form a partition of unity. This property enables straightforward reconstructions using the localised POD basis. Beyond certain technical differences from previously proposed methods, a central contribution of the present work is to draw a connection between Schur-product-based localisation techniques used in data assimilation and localisation strategies for reduced-order modelling. The proposed method is characterised by the use of overlapping localisation functions, the extension of POD in its snapshot-based formulation, and the enforcement of global orthogonality and optimality (in a sense defined later) of the resulting localised basis.

The theoretical lower bound on the error achievable by a reduced-order basis is formalised through the concept of Kolmogorov N-width (Arbes et al., 2025). While rapid convergence can be attained for certain problems (Maday et al., 2002), the convergence rate is known to deteriorate significantly in the case of advection-dominated flows (Olberger and Rave, 2016) and wave propagation problems (Greif and Urban, 2019). For wave propagation, this degradation is particularly pronounced when the wavelength is small compared to the size of the computational domain. In such cases, the flow exhibits strong spatial variability, while the physical structures of interest must still be represented globally due to long-range propagation effects. We evaluate the proposed localisation technique on a wave field governed by the rotating shallow water (RSW) equations, where waves propagate and are scattered by a fluctuating jet. The configuration is similar to that described in Maingonnat et al. (2025), but with a domain significantly larger than the wavelength. By restricting the analysis to a small number of snapshots, the setup highlights a regime where localisation is necessary. We
126 demonstrate that the proposed algorithm yields improved reconstructions at a fixed energy content, and enables a more accurate recovery of the first ideal POD mode (i.e. the mode computed from a large number of snapshots).

129 Section 2 introduces classical methods on which the localised POD is built. First, a brief presentation of the POD procedure is provided in section 2.1. Secondly classical localisation techniques employed in data assimilation are described in section 2.2. The proposed localised POD methodology is presented in section 3. Numerical results applied to a wave field emerging from a RSW system are given in section 4. In section 5, the localised POD is applied to the sea surface height signature of an internal wave field in the Gulf Stream region, using the realistic numerical simulation eNATL60. Some conclusions are drawn in section 6.

137 **2 Classical methods**

138 Localised POD is built by applying localisation techniques, classically employed in data assimilation, on the POD method. The aim of this section is to briefly present these two well established techniques.

141 2.1 Proper orthogonal decomposition

142 Starting from a space-time dependent flow variable $q \in L^2(\Omega, \mathbb{C}^d) \times \mathbb{R}^+$ defined on the bounded spatial domain $\Omega$, we define the spatial inner product between two functions $q_1$ and $q_2$ as

$$ (q_1, q_2)_\Omega = \int_{\Omega} q_1^*(x, t) W_q q_2(x, t) dx, \quad (1) $$

145 with $\cdot$* the conjugation operation and $W_q$ an Hermitian matrix. The associated norm corresponds to an energy $E = \frac{1}{2} \|q\|_{\Omega}^2$, and we define the temporal average

$$ \langle q \rangle_T = \frac{1}{T} \int_0^T q(x, t) dt. \quad (2) $$

147 We can note that variations of the definition of these innerproducts and average operators leads to variants of POD (as for instance in Iollo et al., 2000). The POD consists in searching, from a sequence of realisations of $q(x, t)$ through a time sequence $\{t_n, n = 1, \dots, N_t\}$, for the function $\Phi$ the most correlated with $q$. This translates in maximising the function

$$ \lambda_{\max} = \max_{\Phi \in L^2(\Omega)} \frac{\langle |(q, \Phi)|_{\Omega}^2 \rangle_T}{\|\Phi\|_{\Omega}^2}, \quad \text{s.t. } \|\Phi\|_{\Omega} = 1. \quad (3) $$

152 It can be shown (see for instance Cordier and Bergmann, 2002, for a detailed derivation) that maximising the problem (3) is equivalent to solve the Fredholm equation

$$ \int_{\Omega} \mathbf{R}(x, x') \Phi_i(x') dx' = \lambda_i \Phi_i(x), \quad (4) $$

154 with the correlation tensor $(\mathbf{R}(x, x'))^{ij} = \langle q^i(x) q^{j*}(x') \rangle_T$, with $(i, j) \in [1, \dots, d]$ the indices associated with the flow components.
In practice when the spatial dimension is discretised on $N_x$ degrees of freedoms, it is common that $N_x \gg N_t$. In this case, $\mathbf{R}$ is of very large dimension and may be impossible to store. The snapshot POD algorithm is then well suited. By defining the vector $\mathbf{q}_k \in \mathbb{C}^N$, with $N = N_x \times d$, associated with the spatial discretisation of $\mathbf{q}(x, t_k)$, we assemble the matrix of snapshots

$$Q = [\mathbf{q}_1, \dots, \mathbf{q}_{N_t}]. \qquad (5)$$

Instead of approaching $\mathbf{R}$ by $\frac{1}{N_t} \mathbf{Q} \mathbf{Q}^*$ with $\cdot^*$ the transpose-conjugate operation, we consider instead the correlation matrix

$$C = \frac{1}{N_t} \mathbf{Q}^* \mathbf{W} \mathbf{Q}, \qquad (6)$$

with $\mathbf{W}$ the matrix containing $\mathbf{W}_q$ times the weight quadrature coefficients associated with the space integral in the inner product (1). The POD modes can then be computed by solving the eigenvalue problem

$$Ca_i^* = \lambda_i a_i^*, \qquad (7)$$

by enforcing the scaling of the temporal modes $\langle a_i^2(t) \rangle_T = \lambda_i$, with $(a_i)^k = a_i(t_k)$. The spatial modes can then be computed as follows

$$\Phi_i(\mathbf{x}) = \frac{1}{\lambda_i} \langle \mathbf{q}(\mathbf{x}, t) a_i^*(t) \rangle_T. \qquad (8)$$

This procedure recovers the POD eigen-elements $(\lambda_i, \Phi_i)$ of (4), and the snapshots can be reconstructed as such

$$\mathbf{q}(\mathbf{x}, t) = \sum_{i=1}^{N_t} \Phi_i(\mathbf{x}) a_i(t). \qquad (9)$$

Truncating this modal decomposition allows to obtain a reduced-order representation of the solution. In the next section, we present the concept of localisation as employed in the framework of data assimilation.

## 2.2 Localisation in data assimilation

The aim of data assimilation is to combine heterogeneous sources of informations to estimate at best the state of the system, and eventually perform some forecasts. We seek the initial conditions – and eventually some parameters – of a model, that matches at best some observations $\mathcal{Y}$. Distinctions are often made between statistical and variational methods (Carrassi et al., 2018). For the latter, a cost function is typically sought to be minimised (here written for a 3D-Var problem)

$$\mathcal{J}(\mathbf{q}) = \frac{1}{2} \|\mathbb{H}(\mathbf{q}) - \mathcal{Y}\|_{R_o^{-1}}^2 + \frac{1}{2} \|\mathbf{q} - \mathbf{q}^b\|_{R_b^{-1}}^2, \qquad (10)$$

where $\mathbb{H}$ is an observation operator mapping a state to its representation in the observation space. A regularisation term measuring the distance to a background state $\mathbf{q}^b$ aims at penalising deviations from an a priori given by the model. The relative confidence in these sources of information are controlled by the covariance matrices $R_o$ and $R_b$, through the Mahalanobis norms $\|\cdot\|_{R_o^{-1}}$ and $\|\cdot\|_{R_b^{-1}}$, respectively. While the observation covariance matrix $R_o$ is representative of measurements uncertainties which
186

187

could be estimated, here we will discuss the definition of the background covariance matrix $R_b$, and show the analogy with the POD problem.

188

189

190

191

192

193

194

195

196

197

198

199

A wide class of data assimilation methods exist in the literature, which can differ especially by various treatments of the time component, such as in 4D-Var (Le Dimet and Talagrand, 1986) or Kalman filters. It can be noted that Kalman filters, which are sequential statistical techniques, can be connected to a cost functional of the type of (10). For practical reasons, ensemble methods (Evensen, 1994) are often considered. Indeed, it provides a sampling of states, giving a proxy for the probability density functions, and saves the derivation and implementation of an adjoint equation. In this context, the background covariance tensor $R_b$ can be approximated by its empirical covariance from the background ensemble. Due to the limited number of ensemble members the empirical covariance matrix is low rank and spurious correlations arise between distant locations, where the flow variables are supposed to be uncorrelated. Indeed, the empirical covariance matrix is

$$R_b = Q'Q'^* \in \mathbb{R}^{N \times N}, \qquad (11)$$

200

201

202

where $Q'$ contains in columns the normalised anomaly ensemble members $\{\frac{1}{\sqrt{p-1}}(q_k - \bar{q})\}$, $k = 1, \dots, p$}, with $p$ the number of ensemble members and $\bar{q}$ the ensemble average. Its rank is $p-1$ at most.

203

204

205

206

To circumvent this issue, one may define a localisation matrix $\rho(x, x') \in \mathbb{R}^{N \times N}$, such that $\rho \mapsto 0$ when $\|x - x'\| \gg r$, with $r$ a localisation radius (eventually space-dependent). Gaspari-Cohn functions (Gaspari and Cohn, 1999) are often considered. A localised covariance matrix is then defined (Houtekamer and Mitchell, 2001) as

$$\tilde{R}_b = \rho \circ R_b, \qquad (12)$$

207

208

where $\circ$ denotes the Schur product, which is simply an element-wise multiplication. Properties of this product are given in appendix A.

209

210

As in POD, this covariance tensor is too large to be computed or stored. The localisation matrix is then approximated through a low-rank square-root decomposition

$$\rho \approx \rho^{1/2} \rho^{*/2}, \qquad (13)$$

211

212

213

214

with $\rho^{1/2} \in \mathbb{R}^{N_x \times N_\rho}$, and $N_\rho$ the number of localisation functions considered in the decomposition. Based on such decomposition, several techniques have been proposed to obtain an efficient algorithm (see for instance Farchi and Bocquet, 2019, comparing several algorithms). In section 3, we adapt these concepts to the POD problem.

215

## 3 Localised POD

216

217

218

219

220

221

222

223

224

For a small number of snapshots, POD modes can show some strong convergence artefacts, which can be critical for applications in reduced-order modelling. In the present paper, we assume that there is no additional available snapshots, and we aim at mitigating these artefacts using localisation. While previous methods have focused on directly applying partition of unity techniques in the Galerkin procedure (Ibáñez et al., 2019; Smetana and Taddei, 2023), here, we propose a strategy wherein we localise the covariance tensor $R$, and rewrite the problem in order to obtain a localised equivalent of the snapshot POD algorithm. Finally, we add a partition of unity constraint, in order to be able to reconstruct a snapshot, similarly as in (9).
225 First, we localise the correlation tensor and decompose the localisation matrix as
226 in (13):

$$
\tilde{R} = \rho \circ R = \left( \sum_j \rho_j^{1/2} \rho_j^{*/2} \right) \circ (QQ^*), \quad (14)
$$

227 with $\rho^{1/2} = [\rho_1^{1/2}, \dots, \rho_{N_\rho}^{1/2}]$. Here, we specify a set of localisation functions $\rho_j^{1/2}$, which
228 define implicitly a localisation matrix $\rho = \rho^{1/2}\rho^{*/2}$. As a consequence, it is a strict
229 equality. By using the Schur product calculation rules (appendix A), we obtain

$$
\tilde{R} = \sum_j D_{\rho_j^{1/2}} Q Q^* D_{\rho_j^{1/2}}^* = \sum_j \left( D_{\rho_j^{1/2}} Q \right) \left( D_{\rho_j^{1/2}} Q \right)^*, \quad (15)
$$

230 with $D_{\rho_j^{1/2}} = \operatorname{diag}(\rho_j^{1/2})$. We define $\tilde{Q}_j = D_{\rho_j^{1/2}} Q = \{\tilde{q}_{1,j} \dots \tilde{q}_{N_t,j}\}$, and obtain

$$
\tilde{R} = \sum_j \tilde{Q}_j \tilde{Q}_j^*. \tag{16}
$$

231 It can be seen that $\tilde{q}_{i,j}(\boldsymbol{x}) = \rho_j^{1/2}(\boldsymbol{x})q_i(\boldsymbol{x})$ is the $i^{th}$ snapshot multiplied by the $j^{th}$
232 localisation function. In this way, the localised covariance matrix can be viewed as
233 the empirical covariance matrix computed using the $N_t \times N_\rho$ synthetic snapshots $\tilde{q}_{i,j}$.
234 As a consequence, we artificially increase the rank of the covariance matrix, which is
235 now $N_t \times N_\rho$, at most. The proposed method can thus be viewed as a kind of data
236 augmentation technique.
237

237 In equation (14), the POD problem is artificially modified. In order to ensure that
238 the POD problem is still able to reconstruct the data, we add a reconstruction con-
239 straint, stating that the localisation functions $\rho_j^{1/2}$ constitute a partition of unity

$$
\sum_{j=1}^{N_{\rho}} \rho_j^{1/2}(\mathbf{x}) = 1. \tag{17}
$$

240 Indeed, doing this, we can verify that

$$
\[
\mathbf{q}(\mathbf{x}, t_i) = \mathbf{q}(\mathbf{x}, t_i) \underbrace{\sum_{j=1}^{N_{\rho}} \rho_j^{1/2}(\mathbf{x})}_{1} = \sum_{j=1}^{N_{\rho}} \mathbf{q}(\mathbf{x}, t_i) \rho_j^{1/2}(\mathbf{x}) = \sum_{j=1}^{N_{\rho}} \tilde{\mathbf{q}}_{i,j}(\mathbf{x}). \tag{18}
\]
$$

241 This property will simplify further reconstructions.

242 We now simply need to apply the snapshot POD algorithm to $\tilde{q}_{i,j}$. We define the
243 localised correlation matrix

$$
C_{ij,mn}^{L} = \frac{1}{N_t} \int_{\Omega} \tilde{q}_{ij}^{*}(\mathbf{x}) \tilde{q}_{mn}(\mathbf{x}) \, d\mathbf{x}, \quad (19)
$$

244 which belongs to $\mathbb{C}^{N_t N_\rho \times N_t N_\rho}$ due to the rank inflation, and compute the temporal
245 localised POD modes

$$
C^L a_k^{L*} = \lambda_k^L a_k^{L*}. \tag{20}
$$
246 Here, $a_k^L$ gathers the temporal coefficients $a_k^L(t_i, j)$ (time $t_i$) of the k-th mode for each localisation function (index $j$). It is assumed to be normalised such that

$$ \frac{1}{N_t} \sum_{ij} |a_k^L(t_i, j)|^2 = \lambda_k^L. \qquad (21) $$

248 The POD modes can then be obtained as such

$$ \Phi_k^L(\mathbf{x}) = \frac{1}{N_t \lambda_k} \sum_{ij} \tilde{q}_{ij}(\mathbf{x}) a_k^{L*}(t_i, j). \qquad (22) $$

Finally, due to the reconstruction property (18) the snapshots can be reconstructed as such:

$$ q(\mathbf{x}, t_i) = \sum_{j=1}^{N_\rho} \sum_{k=1}^{N_t \times N_\rho} \Phi_k^L(\mathbf{x}) a_k^L(t_i, j). \qquad (23) $$

In classical POD, linear constraints such as homogeneous boundary conditions or divergence-free property are inherited by the modes – and then by the reconstructed fields – through linear combinations. A caveat of the proposed method lies in the fact that the pre-multiplication of the snapshots $q$ by the localisation functions $\rho_j^{1/2}$ may break these properties, which become only approximately respected. The larger the gradients of $\rho_j^{1/2}$, the stronger may be the effect. In addition, in appendix B, we detail how the localisation affects the optimality of the decomposition.

Localising the correlation tensor associated with the PU constraint represents the main methodological contribution of this work. The demonstrated equivalence between Schur-product localisation and snapshot-based localisation makes the computation tractable on the one hand, while on the other hand it establishes a direct connection between PU-based localisation (employed in reduced-order modelling) and the underlying manipulation of the correlation tensor (in line with Schur-product localisation commonly employed in data assimilation). The main objective is to circumvent convergence issues often encountered in POD (e.g. Blanco et al., 2022), while preserving a physical structure of the POD modes. These aspects are explored numerically in section 4.

## 4 Rotating shallow water equations

### 4.1 Numerical simulation and database

The localisation is experimented on POD of an inertia-gravity wave field interacting non-linearly with a turbulent jet in quasi-geostrophic balance. It is the solution of a rotating shallow water (RSW) model (Zeitlin, 2018), which reads in nondimensional form:

$$ \partial_t v_t + R_o(v_t \cdot \nabla) v_t + f v_t^{\perp} + \nabla h_t = 0, \qquad (24a) $$

$$ \partial_t h_t + \nabla \cdot ((B_u + R_o h_t) v_t) = 0, \qquad (24b) $$

270 with the zonal ($x$) and meridional ($y$) coordinates and $\mathbf{x} = (x, y)^T$. Here, the zonal velocity, meridional velocity and sea-surface anomaly are denoted $u_t$, $v_t$ and $h_t$, respectively; and $v_t^{\perp} = (-v_t, u_t)^T$. The Rossby number $R_o = U/(Lf)$ characterises the
typical flow velocity $U$ compared to the Coriolis frequency $f$ and its dominant length-scale $L$ (width of the jet). The Burger number $B_u = R_d^2/L^2$ measures the Rossby deformation radius $R_d = \sqrt{gH}/f$ compared to the typical lengthscale, with $g$ the gravitational acceleration and $H$ the water depth. The numerical simulation, implemented with the Dedalus library (Burns et al., 2020), is performed in a configuration similar to Maingonnat et al. (2025). In contrast to their configuration, the domain size ($L_x = 40$, $L_y = 60$) is large compared to the wavelength ($\lambda_y \approx 1$), which requires localisation techniques. The Rossby and Burger numbers are $R_o = 0.2$ and $B_u = 1$, respectively, and the (nondimensional) Coriolis parameter $f = 1 + \beta y$ with $\beta = 0.03$. The domain is discretised with Fourier modes associated with a grid of size ($N_x, N_y$) = (256, 384). Sponge zones in the North and South domain allow us to break homogeneity in $y$ and to enforce an incident plane wave propagating northward. To extract the wave field, a complex demodulation of $q_t = (u_t, v_t, h_t)^T$ at the angular frequency $\omega_d = 2$ (corresponding to the frequency of the incoming wave) is performed:

$$q(x, t) = \langle q_t e^{-i\omega_d t} \rangle, \qquad (25)$$

with $\langle \cdot \rangle$ denotes a low-pass filter, with a cut-off frequency $\omega_c = 0.63$. The variable $q(x, t)$ represents then the slowly varying envelop of the wave field $q_t^{(w)} = 2\Re(q(x, t)e^{i\omega_d t})$. The real part of a snapshot of $q = (u, v, h)^T$ is displayed in figures 1(d) to 1(f), besides a visualisation of the jet flow $q_t^{(j)} = \langle q_t \rangle$ in figures 1(a) to 1(c). It shows the wave propagating northward, which is almost plane in the South region, is perturbed by the jet (centered around $y = -1$) and becomes disorganised in the North region.

POD and its localised versions will be applied on $q(x, t)$, i.e. on the complex demodulated wave field (this approach was introduced as broadband POD in Maingonnat et al. (2025)). An innerproduct and associated weight matrix $W$ are defined such that the square of the induced norm corresponds to the linearised total energy. To highlight the interest of localisation, a reduced number of snapshots $N_t = 66$, taken every $\Delta t = 0.3$, will be considered. As will be confirmed below, this number is not enough to converge the POD. Besides, another POD with a larger number of snapshots $N_t^g = 371$ (while keeping the same time interval between snapshots) will provide us a reference basis.

## 4.2 Localisation functions

Four different localisation functions are considered. In order to build easily partition of unity functions on a rectangular domain, we define finite-element-like polynomials associated with $p \times p$ rectangular elements of equal size. They are built by considering 1D Lagrange interpolating polynomial functions $\phi_i^x(x)$ and $\phi_j^y(y)$ of degree $d_g$ in the $x$ and $y$ directions, respectively, and constructing the associated 2D function $\phi_k(x, y) = \phi_i^x(x)\phi_j^y(y)$, with the range of $k$ covering all $(i, j)$ combinations. As in the assembling procedure in finite elements, continuity is enforced (except for $Q_0$, see below), by gathering in the same shape localisation functions all polynomials sharing the same non-null degree of freedom, up to a reindexing of the functions.

$Q_0$: first, piece-wise constant elements with $p = 3$ denoted $Q_0$ are considered. This is equivalent to perform a naïve localisation by splitting the domain into non-overlapping patches. Due to the lack of overlap, the correlation matrix $C_{ij,mn}^L$ in
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Real part of a snapshot of the low-pass filtered jet component ((a) to (c)) and complex-demodulated wave component ((d) to (f)) field.</figcaption></figure>

equation (19) is block-diagonal and the methods boils down to performing independent PODs on each element. As a consequence, such a procedure is computationally efficient, but it may lead to discontinuities of the reconstructed field at the interfaces between patches.

$Q_1$: linear polynomials are considered with $p = 3$. These simple functions have the advantage to provide an overlap between localisation functions and to be easily implemented.

$Q_2^{p \times p}$: quadratic polynomials are also considered, with $p = 2$ which lead to a similar number of localisation functions as $Q_1$, and with $p = 3$ which keeps the number of elements fixed. It can be noted that with $p = 3$, a large number of localisation functions are generated, and increasing more the degree of the polynomials may lead to large computational/memory costs.

It can be noticed that increasing the polynomial order with a fixed element size increases the number of localisation functions, increases the overlap, and results then in smaller localisation function shapes. This corresponds to a smaller effective localisation radius, in the sense that if we build the matrix $\rho$ equation (13), it would result in a sharper function. In the following, we will say by abuse that a smaller localisation radius "localise more". The localised bases can be sorted by "intensity" of localisation as follows: $Q_0, Q_1, Q_2^{2 \times 2}$ and $Q_2^{3 \times 3}$. This ordering is supported by the following results. Due to the fast increase of number of localisation functions, only $p = 2$ and $p = 3$ are considered.
<figure><img src="image_3.png" /><figcaption>Fig. 2 POD spectrum, and relative information criteria (RIC). Comparison between non-localised POD (black) and localised POD using $Q_0$ (red), $Q_1$ (blue), $Q_2^{2 \times 2}$ (green) and $Q_2^{3 \times 3}$ (pink). In the normalised RIC (c), the abscissa is normalised by the total number of modes.</figcaption></figure>

In the present section, we focus on finite-element functions, but partition of unity can be easily achieved (see for instance Babuška and Melenk, 1997) by considering a set of functions $w_k(x)$ – here written in 1D – and by renormalising them as $\phi_k(x) = w_k(x)/\sum_j w_j(x)$. In section 5, Gaussian functions for $w_k(x)$ will be considered.

## 4.3 Localised POD

Figure 2 compares the POD spectrum built with the small amount of data. Non-localised and localised versions are compared. Figure 2(a) shows the POD eigenvalues. As expected, the higher the degree $d_g$ – or number of localisation functions –, the slower the decay of the eigenspectrum. Since the rank of the covariance matrix has been artificially increased, the total number of modes increases. In all these curves, we have $\sum_i \lambda_i$ constant, equal to the total energy contained in the snapshots. This explains a reduction of the amplitude of the first modes as the localisation effect is more pronounced (higher degree or more elements).

The relative information criteria (RIC), defined as

$$ \text{RIC}(i) = \frac{\sum_{j=1}^{i} \lambda_j}{\sum_{j=1}^{N_t} \lambda_j}, \qquad (26) $$

is a standard criterion measuring the energy contained by the projection of the snapshots on the $i$ first modes. By orthonormality of the basis, it is a monotonic quantity increasing progressively until $\text{RIC}(N_t) = 1$. This is displayed in figure 2(b). As expected, the more the localisation, the slower the increase of RIC. This visualisation
may be inadequate to represent the ability of the subspace to span the data. We have thus complemented the plot by figure 2(c) where the abscissa is normalised by the number of available modes. Compared to figure 2(b), it exhibits a reverse trend, and the RIC of the non-localised POD is very progressive until the total number of modes is reached. This shape is typical of a lack of data to converge the POD modes. Indeed, the RIC is only informative of the energy contained in the available snapshots, and such a profile suggests that a new snapshot will add additional information, i.e. will not be fully spanned by the previous snapshots. On the contrary, localisation sharpens the profile, by showing a fast increase, then a plateau. This indicates that the rank increase by localisation is enough to encompass the available snapshots, which in turn suggests that it may possibly represent unseen additional snapshots and that a reduced-order representation is possible. On the other hand, a too sharp profile indicates that the rank has been increased too much, which is at the price of generating a very large number of modes. This is incompatible with model reduction, and such a behaviour can be detected by a too slow increase in figure 2(b). A compromise is then expected.

The first three POD modes are shown in figures 3 and 4. The reference modes (figures 3(a) to 3(c)) clearly show waves deflected eastward for mode 1, and westward for mode 3. Mode 2 is a stationary wave, combining correlated eastward and westward deflections. This physical behaviour has been reported in Maingonnat et al. (2025) for smaller domains. The non-localised modes (figures 3(d) to 3(f)) exhibit no clear pattern, confirming their poor convergence. The $Q_0$ localisation results in modes with non-zero spatial support within a given patch (or element) and zero elsewhere. This is expected due to the block-diagonal structure of the correlation matrix (19), and corresponds to POD modes computed independently in each patch. These modes reveal deflection structures, but the lack of overlap between patches hinders the ability to perform smooth reconstructions. Modes associated with $Q_1$ localisation (figures 3(j) to 3(l)) recover a global structure, and are remarkably able to reproduce deflection and stationary patterns of the first two modes. Increasing the localisation effect using $Q_2^{2 \times 2}$, and more markedly with $Q_2^{3 \times 3}$, leads to very (and possibly excessively) localised modes. These modes are well suited as reconstruction bases but fail to capture a physically meaningful global structure. As with $Q_0$, deflection patterns are reproduced, but here the modes exhibit spatial overlap.

## 4.4 Reconstructions

In this section we evaluate the ability of the various POD bases to reconstruct a snapshot which does not belong to the small dataset used to build the POD basis, but belongs to the large dataset used to build the reference POD basis. The new dataset is issued from the same simulation, but at a time far enough to ensure decorrelation with the training set. To perform fair comparisons, the reconstructions are performed at fixed RIC: RIC = 0.6 in figure 5 and RIC = 0.9 in figure 6. For each reconstruction, normalised root mean square error of the snapshot is provided in the captions:

$$ \text{RMSE} = \frac{\|q^{\text{proj}} - q\|_{\Omega}}{\|q\|_{\Omega}}, \qquad (27) $$

where $q$ is the snapshot and $q^{\text{proj}}$ is its projection onto the considered basis.

For RIC = 0.6, reconstruction using the reference POD (Figures 5(b)) allows to represent the main interference patterns. However, the reconstruction in figure 5(c) is
<figure><img src="image_4.png" /><figcaption>(a) POD reference mode 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) POD reference mode 2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) POD reference mode 3.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) POD non-loc. mode 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) POD non-loc. mode 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) POD non-loc. mode 3.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(g) POD $Q_0$ mode 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(h) POD $Q_0$ mode 2.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(i) POD $Q_0$ mode 3.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(j) POD $Q_1$ mode 1.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(k) POD $Q_1$ mode 2.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(l) POD $Q_1$ mode 3.</figcaption></figure>

**Fig. 3** Meridional component *v* of the 3 first POD modes. Comparison between the reference POD built on a large dataset ((a) to (c)), the non-localised POD built on a small dataset ((d) to (f)) and localised with $Q_0$ ((g) to (i)) and $Q_1$ ((j) to (l)) localisation functions.
<figure><img src="image_16.png" /><figcaption>**Fig. 4** Meridional velocity component, $v$, of the 3 first POD modes (Figure 3 continued). Localised POD with $Q_2^{2 \times 2}$ ((a) to (c)) and $Q_2^{3 \times 3}$ ((d) to (f)) localisation functions.</figcaption></figure>

worst, which shows that the small dataset is not enough to span a snapshot which does not belong to the training dataset. Employing localisation improves the reconstruction, and the RMSE decreases gradually as the localisation is stronger, until approaching the value associated with the reference POD. As expected, $Q_0$ reconstructions, which correspond to a naïve patch reconstruction, suffers from discontinuities.

For $RIC = 0.9$, reconstructions are accurate for all bases, which is a consequence of the high RIC value, but differences of performances can be observed between the various bases. The same trend than for $RIC = 0.6$ is observed. Moreover, even for these high RIC, $Q_0$ still suffers from discontinuities artefacts.

The RMSE are reported on table 1, and complemented with $RIC = 0.3$ and $RIC = 1$. The gradual improvement of RMSE as more localisation is employed is confirmed for all RIC. The reconstruction RMSE belong systematically between the reference POD and non-localised POD. This behaviour shows the ability of the localised basis to span snapshots which have not been seen before, and then to mitigate the rank deficiency associated with the small amount of data. Especially, for $RIC = 1$, all modes are kept, and even for this case, the non-localised POD basis lead to strong RMSE (62% of the snapshot RMS). It quantifies the projection residual of the reconstructed snapshot on the space spanned by the training dataset. In this case, localisation drastically reduces the RMSE.

Table 2 provides the number of modes in each basis for several RIC. This shows the penalty in terms of basis size to obtain an alignment improvement. As in figure 2(b), the higher the RIC or the localisation, the larger the number of modes.
<figure><img src="image_17.png" /><figcaption>(a) Original snapshot.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) POD reference.
RMSE = 0.5</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) POD non-localised.
RMSE = 0.87</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) POD $Q_0$.
RMSE = 0.8</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(e) POD $Q_1$.
RMSE = 0.7</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(f) POD $Q_2^{2 \times 2}$.
RMSE = 0.68</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(g) POD $Q_2^{3 \times 3}$.
RMSE = 0.64</figcaption></figure>

**Fig. 5** Reconstructions of the incoherent fluctuation of the meridional velocity component *v* of a snapshot outside the training set using POD modes with different localisations, at fixed RIC = 0.6. The root mean square error (RMSE) is reported in the subcaptions.

<table><caption>**Table 1** Relative root-mean-square error of the reconstructed snapshot.</caption><thead><tr><th>RIC</th><th>POD ref.</th><th>POD non-loc.</th><th>Q<sub>0</sub></th><th>Q<sub>1</sub></th><th>Q<sub>2</sub><sup>2×2</sup></th><th>Q<sub>2</sub><sup>3×3</sup></th></tr></thead><tbody><tr><td>0.3</td><td>0.76</td><td>0.92</td><td>0.89</td><td>0.87</td><td>0.86</td><td>0.85</td></tr><tr><td>0.6</td><td>0.50</td><td>0.87</td><td>0.80</td><td>0.70</td><td>0.68</td><td>0.64</td></tr><tr><td>0.9</td><td>0.30</td><td>0.68</td><td>0.53</td><td>0.46</td><td>0.42</td><td>0.37</td></tr><tr><td>1.0</td><td>0</td><td>0.62</td><td>0.23</td><td>0.12</td><td>0.092</td><td>0.055</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>**Fig. 6** Reconstructions of the incoherent fluctuation of the meridional velocity component *v* of a snapshot outside the training set using POD modes with different localisation, at fixed RIC = 0.9. The root mean square error (RMSE) is reported in the subcaptions.</figcaption></figure>

<table><caption>**Table 2** Number of modes required to obtain a given RIC.</caption><thead><tr><th>RIC</th><th>POD ref.</th><th>POD non-loc.</th><th>Q<sub>0</sub></th><th>Q<sub>1</sub></th><th>Q<sub>2</sub><sup>2×2</sup></th><th>Q<sub>2</sub><sup>3×3</sup></th></tr></thead><tbody><tr><td>0.3</td><td>7</td><td>3</td><td>11</td><td>11</td><td>14</td><td>22</td></tr><tr><td>0.6</td><td>24</td><td>9</td><td>37</td><td>38</td><td>49</td><td>76</td></tr><tr><td>0.9</td><td>100</td><td>28</td><td>129</td><td>138</td><td>176</td><td>264</td></tr><tr><td>1.0</td><td>371</td><td>66</td><td>602</td><td>1069</td><td>1645</td><td>3082</td></tr></tbody></table>
<table><caption>**Table 3** Alignment metric $\beta^L$ for various localisations obtained for a given RIC.</caption><thead><tr><th>RIC</th><th>POD ref.</th><th>POD non-loc.</th><th>Q<sub>0</sub></th><th>Q<sub>1</sub></th><th>Q<sub>2</sub><sup>2×2</sup></th><th>Q<sub>2</sub><sup>3×3</sup></th></tr></thead><tbody><tr><td>0.3</td><td>1.0</td><td>0.74</td><td>0.73</td><td>0.77</td><td>0.71</td><td>0.77</td></tr><tr><td>0.6</td><td>1.0</td><td>0.78</td><td>0.84</td><td>0.89</td><td>0.90</td><td>0.91</td></tr><tr><td>0.9</td><td>1.0</td><td>0.85</td><td>0.93</td><td>0.95</td><td>0.97</td><td>0.98</td></tr></tbody></table>

4.5 Subspace inflation

To provide a more quantitative analysis, we propose to compare an alignment criterion similar to that introduced in Cavalieri et al. (2013). We compare the alignment between the first reference POD mode $\Phi_1^{\text{ref}}$ using the large data set, with its projection onto the small data set (localised or not) POD basis. This aims at quantifying the ability of the POD basis to span unseen snapshots, whose first reference POD mode is the function the most aligned with the statistically stationary data. This is defined as

$$ \beta^L = \frac{|\langle \Phi_1^{\text{ref}}, \Phi_1^{\text{proj}} \rangle_{\Omega}|}{\|\Phi_1^{\text{proj}}\|_{\Omega}}, \quad \text{with} \quad \Phi_1^{\text{proj}} = \sum_{j=1}^{M} \langle \Phi_j^L, \Phi_1^{\text{ref}} \rangle_{\Omega} \Phi_j^L. \qquad (28) $$

The various definitions of the localised (or not) basis $\{\Phi_j^L, j = 1..., M\}$ will be compared. The truncation order $M$ is fixed based on the RIC criterion. The criterion $\beta^L$ is equal to 1 if the vectors are perfectly aligned, i.e. if the reference POD mode is fully spanned by the basis of $\Phi_i^L$. $\beta^L$ becomes null if the vectors are orthogonal.

The $\beta^L$ metrics are provided in table 3. This shows that, even for high RIC, the non-localised POD basis leads to low values of $\beta^L$ ($\beta^L = 0.85$ for RIC = 0.9). This confirms the non-convergence of the non-localised POD basis. As a sanity check, the reference POD basis has $\beta^L = 1$, since $\Phi_1^{\text{ref}}$ is a member of the basis. By construction, the higher the RIC, the higher $\beta^L$. Moreover as expected, the more the basis is localised, the better the alignment. This shows the ability of localisation to increase artificially the subspace, thus allowing to encompass $\Phi_1^{\text{ref}}$. An exception is made for RIC = 0.3 and $Q_2^{p \times p}$. This can be explained by a higher degree of the interpolating polynomials, suggesting a higher order of accuracy leading to a faster increase of $\beta^L$ with respect to the RIC. This translates in poor performances for low RIC.

Interpreting together tables 2 and 3 highlights the compromise between localisation and cost in terms of number of modes. In the present configuration, by eliminating $Q_0$ due to discontinuities issues, $Q_1$ seems to be a good compromise in terms of number of modes. It can be shown indeed that it meets a similar number of modes as the reference basis, while improving significantly $\beta^L$. The number of modes is slightly larger than the reference, but we have to recall that it has only been built using the small data set. This choice is corroborated by the physical properties of the modes shown in section 4.3 and the good reconstruction behaviour in section 4.4. Finally, it is worth mentioning that the number of modes of $Q_1$ is very similar to $Q_0$, but the quality of reconstruction is drastically improved, which highlights the benefit of considering overlapping localisation functions.
## 5 Realistic numerical simulation

### 5.1 Numerical simulation and database

The localised POD is applied to a high-resolution, realistic numerical simulation of the North Atlantic Ocean, namely the eNATL60 configuration (Brodeau et al., 2020). This simulation is based on the Nucleus for European Modelling of the Ocean (NEMO) model (Madec et al., 2008), which solves the three-dimensional primitive equations under the Boussinesq and hydrostatic approximations. The model domain spans the North Atlantic Ocean from approximately 5°N to 65°N, and includes the Mediterranean and Black Seas. The horizontal resolution is 1/60° (approximately 1.5 km at mid-latitudes), and the vertical coordinate is discretised using 300 levels.

For our analysis, we extract the sea surface height (SSH) from a subdomain encompassing the Gulf Stream region, delimited by the longitude range [73.2°W, 56.3°W] and the latitude range [31.2°N, 44.1°N]. The barotropic tide is removed from the SSH field, and the internal tide signal is isolated via complex demodulation at a central semi-diurnal frequency of $\omega = 1.415 \times 10^{-4} \text{rad s}^{-1}$, corresponding to a period of 12 h 20 min. A total of $N_t = 1000$ snapshots, sampled every 6 h over a time window of $T = 255$ days, are used for training. An additional snapshot, taken 366 days after the dataset start time, is employed to assess reconstruction performance outside the training interval [0, $T$].

Although the spatial structures observed in this simulation are similar to the shallow-water results discussed in section 4, the present case exhibits a substantially greater complexity due to several factors. First, as in the previous test case, the long-range propagation of short internal wave wavelengths leads to a slow convergence of the POD modes. Second, the dynamics are significantly richer, owing to the model's vertical resolution, realistic bathymetry, atmospheric wind forcing, and the broader spatial extent of the eNATL60 domain. Finally, the time window of 255 days remains insufficient to capture the statistical variability of the Gulf Stream, implying that the test snapshot lies significantly outside the subspace spanned by the training data. Thus, this setup provides a challenging benchmark to assess the capacity of localisation strategies to enhance the reduced-order space and to improve the reconstruction of out-of-sample states.

The configuration investigated here is similar to that addressed by Egbert and Erofeeva (2021), who examined a region off the Amazon shelf. Our method can be viewed as an extension of their work, albeit with some minor differences resulting from a more careful design of the localisation strategy.

### 5.2 Localisation functions

The localisation functions are defined consistently with the end of section 4.2, using Gaussian functions for

$$w_{ij}(x, y) = e^{-\frac{(x-x_i)^2}{l_x^2} - \frac{(y-y_j)^2}{l_y^2}},$$

with $(x, y)$ the curvilinear coordinates of the domain $[x_{\min}, x_{\max}] \times [y_{\min}, y_{\max}]$. We define then $\phi_{ij}(x, y) = w_{ij}(x, y) / \sum_{k,l} w_{kl}(x, y)$ to enforce the partition of unity property. The Gaussian widths are chosen $l_x = (x_{\max} - x_{\min}) / (N_{px} - 1)$ and $l_y = (y_{\max} - y_{\min}) / (N_{py} - 1)$, respectively. We choose $x_i = x_{\min} + il_x$ and $y_j = y_{\min} + jl_y$
<table><caption>**Table 4** Relative root-mean-square error of reconstruction of a training and test snapshot in the Gulf Stream using 400 modes. The number of localised snapshots $N_t \times N_\rho$ is reported in last line.</caption><thead><tr><th></th><th>Non loc.</th><th>loc. 3 × 3 Gaussian</th><th>loc. 3 × 3 Patches</th><th>loc. 10 × 10 Gaussian</th></tr></thead><tbody><tr><th>Training</th><td>0.0003</td><td>0.0075</td><td>0.26</td><td>0.25</td></tr><tr><th>Test</th><td>0.69</td><td>0.58</td><td>0.56</td><td>0.50</td></tr><tr><th>N<sub>t</sub> × N<sub>ρ</sub></th><td>1000</td><td>9000</td><td>9000</td><td>10000</td></tr></tbody></table>

with $i \in [0 : N_{px} - 1]$ and $j \in [0 : N_{py} - 1]$. We perform tests with $(N_{px}, N_{py}) = (3, 3)$ and $(10, 10)$, respectively. In addition, $Q_0$ functions with $p = 3$ are also considered.

In the case of $10 \times 10$ Gaussian functions, because it generates a large number $N_t \times N_\rho$ localised snapshots, we sub-sample the snapshots by taking only $N_t = 100$ (instead of 1000) snapshots over the same time window $[0, T]$.

## 5.3 Reconstructions and modal decomposition

The SSH of the slow flow component is displayed in figure 7(a), which shows the signature of the Gulf stream jet and of mesoscales. POD and localised POD are applied on the wave SSH anomaly. The figure 7 compares the test snapshot (fig. 7(b)) with reconstructions using 400 modes of several bases. Non-localised POD, $3 \times 3$ Gaussian (fig. 7(d)), $3 \times 3$ $Q_0$ (fig. 7(e)), and $10 \times 10$ Gaussian (fig. 7(f)) localisation functions are considered. It should be noted that a large number of modes (here 400) is required to achieve meaningful reconstructions. Despite this high number, the non-localised solution is degraded: small-scale structures and the overall variance are poorly represented. Localisation, particularly with $3 \times 3$ and more strikingly with $10 \times 10$ Gaussian functions enhances this overall variability. The small-scale structures, which result in an erroneous noisy texture signal in the non-localised POD, are effectively filtered with the $10 \times 10$ Gaussian localisation functions, resulting in much cleaner reconstructions. This highlights the regularising effect of localisation, which mitigates long-range spurious correlations. In addition, the $3 \times 3$ $Q_0$ functions lead to a discontinuous reconstruction, as expected and in agreement with the results presented in section 4.

Reconstruction relative root-mean-square errors are reported in Table 4, for a fixed truncation order of 400. The table presents reconstruction errors for both the test snapshot (corresponding to figure 7) and a snapshot from the training set, for the various localisation functions considered. As expected, increasing localisation leads to higher reconstruction errors on the training set at fixed truncation order. This is consistent with the optimality of the non-localised POD, which minimises reconstruction error over the training data. Conversely, the reconstruction error on the test snapshot decreases with increased localisation. This highlights the ability of the localised bases, at fixed order, to trade off reconstruction accuracy on the training set for improved generalisation to previously unseen data.

Concerning the numerical cost, the size $N_t \times N_\rho$ is reported in last line of table 4. Localisation increases computational cost by inflating the size of the correlation matrix. To limit this cost, we subsampled the training set when moving from the $3 \times 3$ to the $10 \times 10$ localisation functions, in order to maintain a similar overall computational load. Despite using ten times fewer training snapshots and the same truncation order, the $10 \times 10$ Gaussian localisation functions yield significantly better performance
<figure><img src="image_25.png" /><figcaption>**Fig. 7** Sea surface height reconstructions of the test snapshot using 400 modes of several POD bases. The reference (b) is the original test snapshot. The non-localised POD (c) is compared with localised bases using 3×3 Gaussian (d), 3×3 $Q_0$ (e), and 10×10 Gaussian (f) localisation functions.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 8** Eigenvalue spectrum and Relative Information Criteria for the Gulf Stream dataset. Comparison between non-localised and localised with 3 × 3 and 10 × 10 Gaussian functions and 3 × 3 patch functions.</figcaption></figure>

in reconstructing an unseen snapshot. This highlights the efficiency of localisation in improving generalisation, even under constrained data availability.

To assess the effect of truncation order, the POD eigenvalue spectra and RIC values are shown in figure 8. A similar trend to that observed in section 4 for the RSW equations can be identified. It is observed that the 3 × 3 patch-based $Q_0$ localisation leads to stronger localisation than the 3 × 3 Gaussian localisation functions, despite using the same number $N_ρ$. This is due to the absence of overlap in the $Q_0$ functions, unlike the Gaussian functions which have overlapping support. As a result, the effective size of the $Q_0$ functions is smaller than that of the Gaussian ones.

Finally, the spatial modes obtained from the non-localised POD and from localised POD using 3 × 3 and 10 × 10 Gaussian functions are examined. Localisation is observed to suppress long-range correlations, resulting in modes with reduced spatial support. However, this support remains broader than the width of the localisation functions themselves, due to the presence of correlations in the data that can propagate through the overlap of the localisation functions. The 10 × 10 Gaussian localisation is particularly effective in isolating physical structures, notably the internal wave beam generated at the U.S. continental slope. This wave beam interacts with flow structures associated with the Gulf Stream, as described in Bella et al. (2024). Capturing the full dynamical behaviour of this interaction requires several modes, highlighting the complexity of the associated processes. This example illustrates that an appropriate choice of localisation functions enables a modal decomposition that extracts physically meaningful structures. The typical width of the localisation functions may be guided by physical considerations, such as spatial correlation scales, prior physical knowledge, or characteristic length scales derived from the underlying dynamics.

## 6 Conclusion

We have proposed a localised version of POD. The aim of this procedure is to circumvent lack of convergence of the POD basis when the domain is large compared to the physical processes at play, and when a limited number of snapshots are available. This is a common situation when high-resolution numerical simulations are employed. Localised POD has many applications in reduced-order modelling, physical analysis,
<figure><img src="image_27.png" /><figcaption>(a) Non loc. $\Phi_1$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) Loc. $3 \times 3 \Phi_1$</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(c) Loc. $10 \times 10 \Phi_1$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(d) Non loc. $\Phi_2$</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(e) Loc. $3 \times 3 \Phi_2$</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(f) Loc. $10 \times 10 \Phi_2$</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(g) Non loc. $\Phi_3$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(h) Loc. $3 \times 3 \Phi_3$</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(i) Loc. $10 \times 10 \Phi_3$</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(j) Non loc. $\Phi_4$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(k) Loc. $3 \times 3 \Phi_4$</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(l) Loc. $10 \times 10 \Phi_4$</figcaption></figure>

**Fig. 9** Real part of SSH associated with the spatial POD and LPOD modes. Non-localised ((a)-(d)-(g)-(j)), $3 \times 3$ Gaussian ((b)-(e)-(h)-(k)) and $10 \times 10$ Gaussian ((c)-(f)-(i)-(l)) functions.
560 flow control and data assimilation. It may provide as well some consistent data-based comparisons with the localised resolvent analysis proposed by Eichberger et al. (2024).

561 The paradigm is slightly different from classical POD-based reduced-order modelling, whose aim is to reduce the size of the projection subspace. Here, we start with incomplete data leading to a poor basis quality. The strategy is first to increase the subspace through localisation in order to correctly represent the flow features, and in a second step to perform the truncation.

562 The localisation strategy proposed in this paper relies on generating from each snapshot several synthetic snapshots, and to ensure a reconstruction property of the basis by defining a set of localisation functions, which are a partition of unity. A contribution was to link this strategy with an equivalent Schur-product localisation of the covariance tensor, classically employed in data assimilation. We consider overlapping localisation functions to avoid discontinuous reconstructions, that would result from a naïve patch decomposition strategy.

563 We have tested the methodology on an inertia-gravity wave field interacting with a turbulent jet flow in a RSW simulation. Various definitions of localisation functions are compared, and we have demonstrated the ability of localisation to produce globally defined modes, which extract physically relevant structures. We have highlighted the ability of localisation to increase artificially the projection subspace to encompass unseen flow states. We provided some guidelines on the choice of the localisation functions as a trade-off between localisation and number of generated modes. The final choice is obviously configuration dependent, and an important characteristic is that the effective localisation radius should cover physical mechanisms to reproduce, in order to avoid breaking artificially their correlations.

564 Finally, we have demonstrated the ability of localised POD to reconstruct and extract relevant features for realistic applications, specifically focusing on the internal wave field in the Gulf Stream region. These data are particularly challenging for standard POD, due to slow convergence and complex spatial correlations. Localisation was shown to enhance the projection space available for reconstruction, while simultaneously removing spurious long-range correlations in the spatial modes.

## 590 A Calculation rules associated with the Schur product

591 For matrices $(A, B, C) \in \mathbb{C}^{N \times N}$, and vectors $(p, q) \in \mathbb{C}^N$ Schur product is defined as such:

$$ (A \circ B)_{ij} = A_{ij}B_{ij}. \qquad (29) $$

592 Useful properties are provided as follows:

593 $-A \circ B = B \circ A$

594 $-A \circ (B + C) = A \circ B + A \circ C$

595 $-(qp^*) \circ A = D_q AD_p^*$

596 with $D_p = \text{diag}(p)$.

## 597 B Optimality of the localised POD

598 In classical POD, the N-order reconstructed field $q_r = \sum_{i=1}^{N} \Phi_i a_i$ has the following energy

$$ \frac{1}{T} \int_{0}^{T} \|q_r\|_{\Omega}^{2} dt = \sum_{i=1}^{N} \lambda_{i}, \qquad (30) $$
599 which is maximised by definition of the decomposition.

600 Localisation has two effects on the reconstruction optimality. The first is to consider localised snapshots as samples, leading, using the normalisation constraint (21), to

$$ 
\begin{align} 
\sum_{i=1}^{N} \lambda_{i}^{L} &= \sum_{j=1}^{N_{\rho}} \frac{1}{T} \int_{0}^{T} \| \tilde{q}_{r}^{j} \|_{\Omega}^{2} dt \tag{31} \\ 
&\neq \frac{1}{T} \int_{0}^{T} \left\| \sum_{j=1}^{N_{\rho}} \tilde{q}_{r}^{j} \right\|_{\Omega}^{2} dt = \frac{1}{T} \int_{0}^{T} \| q_{r} \|_{\Omega}^{2} dt. 
\end{align} 
$$

602 with $\tilde{q}_r^j = \sum_{i=1}^N \Phi_i^L a_i^L(t, j)$. The inequality is due to the overlapping regions, and becomes an equality only for $\tilde{Q}_0$.

603 The second effect, is a penalty on modes with large spatial support caused by the artificial vanishing of space-correlations. This can be easily understood with the specific case $Q_0$. Here, patches are forced to be uncorrelated with a block-diagonal correlation matrix, and then the localised POD modes are constrained to belong to a unique patch. This constraint leads to a sub-optimal reconstruction compared to the classical POD modes. Overlapping localisation functions authorise the modes to have a larger spatial support, but the reconstruction is still suboptimal. We recall that this suboptimality is the relaxation allowing the rank inflation.

## 611 Declarations

### 612 Fundings

613 This work was supported by the French National program LEFE (Les Enveloppes Fluides et l'Environnement). This research received support from the French research funding agency (ANR) under the ModITO project (ANR-22-CE01-0006-01), and from the CNES through the project DIEGO.

### 617 Acknowledgements

618 The authors thank G. Gabard for fruitful discussions.

## 619 References

620 Anderson S, White C, Farhat C (2023) Space-local reduced-order bases for accelerating reduced-order models through sparsity. International Journal for Numerical Methods in Engineering 124(7):1646–1671, DOI https://doi.org/10.1002/nme.7179

621 Arbes F, Greif C, Urban K (2025) The Kolmogorov n-width for linear transport: exact representation and the influence of the data. Advances in Computational Mathematics 51(2):13

622 Arbogast E, Desroziers G, Berre L (2017) A parallel implementation of a 4DEnVar ensemble. Quarterly Journal of the Royal Meteorological Society 143(706):2073–2083

623 Aubry N, Holmes P, Lumley JL, Stone E (1988) The dynamics of coherent structures in the wall region of a turbulent boundary layer. Journal of Fluid Mechanics 192:115–173, DOI 10.1017/S0022112088001818

624 Babuška I, Melenk JM (1997) The partition of unity method. International Journal for Numerical Methods in Engineering 40(4):727–758, DOI https://doi.org/10.1002/(SICI)1097-0207(19970228)40:4<727::AID-NME86>3.0.CO;2-N

625 Bella A, Lahaye N, Tissot G (2024) Internal tide energy transfers induced by mesoscale circulation and topography across the north atlantic. Journal of Geophysical Research: Oceans 129(8):e2024JC020914, DOI https://doi.org/10.1029/2024JC020914, e2024JC020914 2024JC020914
637 Bergmann M, Cordier L, Brancher JP (2005) Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model. Physics of Fluids 17(9):097101, DOI 10.1063/1.2033624

638 Berkooz G, Holmes P, Lumley JL (1993) The proper orthogonal decomposition in the analysis of turbulent flows. Annual Review of Fluid Mechanics 25(1):539-575

639 Blanco DCP, Martini E, Sasaki K, Cavalieri AVG (2022) Improved convergence of the spectral proper orthogonal decomposition through time shifting. Journal of Fluid Mechanics 950:A9, DOI 10.1017/jfm.2022.790

640 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

641 Blanco DCP, Martini E, Sasaki K, Cavalieri AVG (2022) Improved convergence of the spectral proper orthogonal decomposition through time shifting. Journal of Fluid Mechanics 950:A9, DOI 10.1017/jfm.2022.790

642 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

643 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

644 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

645 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

646 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

647 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

648 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

649 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

650 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

651 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

652 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

653 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

654 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

655 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

656 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

657 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

658 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

659 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732

660 Brodeau L, Le Sommer J, Albert A (2020) ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations. DOI 10.5281/zenodo.4032732, URL https://doi.org/10.5281/zenodo.4032732
697 0870.1986.tb00459.x

698 Lee JCK, Amezcua J, Bannister RN (2024) Variable-dependent and selective multivariate localization for ensemble-variational data assimilation in the tropics. Monthly Weather Review 152(4):1097-1118

701 Loève, M (1955) Probability theory. Van Nostrand

702 Lorenz EN (1956) Empirical orthogonal functions and statistical weather prediction, vol 1. Massachusetts Institute of Technology, Department of Meteorology Cambridge

704 Lumley JL (1967) The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation pp 166-178

706 Maday Y, Patera AT, Turinici G (2002) Global a priori convergence theory for reduced-basis approximations of single-parameter symmetric coercive elliptic partial differential equations. Comptes Rendus Mathématique 335(3):289-294

709 Madec G, et al. (2008) NEMO ocean engine. Note du Pôle de modélisation. Institut Pierre-Simon Laplace (IPSL), France 27:1288-1619

711 Maingonnat I, Tissot G, Lahaye N (2025) Coupled estimation of internal tides and turbulent motions via statistical modal decomposition. Ocean Science 21(2):807-827, DOI 10.5194/os-21-807-2025

714 McKeon BJ, Sharma AS (2010) A critical-layer framework for turbulent pipe flow. Journal of Fluid Mechanics 658:336-382

716 Muñoz E, Dave H, D'Alessio G, Bontempi G, Parente A, Le Clainche S (2023) Extraction and analysis of flow features in planar synthetic jets using different machine learning techniques. Physics of Fluids 35(9):094107, DOI 10.1063/5.0163833

719 Olberger M, Rave S (2016) Reduced basis methods: Success, limitations and future challenges. In: Proceedings of the Conference Algoritmy, pp 1-12

721 Sahyoun S, Djouadi S (2013) Local proper orthogonal decomposition based on space vectors clustering. In: 3rd International Conference on Systems and Control, pp 665-670, DOI 10.1109/ICoSC.2013.6750930

724 Smetana K, Taddei T (2023) Localized model reduction for nonlinear elliptic partial differential equations: Localized training, partition of unity, and adaptive enrichment. SIAM Journal on Scientific Computing 45(3):A1300-A1331, DOI 10.1137/22M148402X

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

730 Zeitlin V (2018) Geophysical fluid dynamics: understanding (almost) everything with rotating shallow water models. Oxford University Press