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

HAL
open science

# Ensemble Filter on dynamically driven RKHS: Application to a Multilayer Quasi-Geostrophic ocean model

Maël Jaouen, Benjamin Dufée, Etienne Mémin, Gilles Tissot

► To cite this version:

Maël Jaouen, Benjamin Dufée, Etienne Mémin, Gilles Tissot. Ensemble Filter on dynamically driven RKHS: Application to a Multilayer Quasi-Geostrophic ocean model. 2025. hal-04915412v2

HAL Id: hal-04915412

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

Preprint submitted on 4 Feb 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ensemble Filter on dynamically driven RKHS: Application to a Multilayer Quasi-Geostrophic ocean model

Maël Jaouen<fnref n="1" />, Benjamin Dufée<fnref n="2" />,<fnref n="3" />, Etienne Mémin<fnref n="1" />, Gilles Tissot<fnref n="1" />

## Key Points:

* An ensemble filter defined on a time evolving RKHS family driven by the system dynamics is presented.

* The corresponding filter enables extending classical Bayesian ensemble filters to efficiently solve a 4D Var assimilation problem.

* Assessed on a multilayer quasi-geostrophic simulation for the North Atlantic Ocean basin in an OSSE framework with SSH observations.

<footnote n="1">¹Univ Rennes, CNRS, Inria, IRMAR - UMR CNRS 6625, F-35000 Rennes, France</footnote>

<footnote n="2">²Univ Bretagne Sud, IUT Vannes, Department of Data Science, 8 rue Michel de Montaigne, 56000, Vannes, France</footnote>

<footnote n="3">³Laboratory of Mathematics in Bretagne Atlantique, Rue de Saint-Maudé, 56100 Lorient, France</footnote>
16

**Abstract**

17

18

19

20

21

22

23

Estimating the ocean’s dynamic state, particularly at the mesoscale, is challenging due to the non-linear nature of ocean processes. Traditional data assimilation methods like optimal interpolation (OI) and four-dimensional variational assimilation (4D-Var) improve ocean state estimates but face trade-offs: OI is computationally efficient but struggles with small-scale features, while 4D-Var is more accurate but computationally demanding. Ensemble methods offer a compromise, with ensemble Kalman filters assuming Gaussianity and particle filters limited by high-dimensional systems.

24

25

26

27

28

29

30

31

This article introduces a novel ensemble filter using a Reproducing Kernel Hilbert Spaces (RKHS) framework, embedding observables in time-evolving RKHSs to relax Gaussianity assumptions and provide computational efficiency comparable to 3D-Var while addressing 4D-Var problems. A key innovation is a tiled assimilation approach that enhances the ensemble's physical properties over time. Applied to a quasi-geostrophic model of the North Atlantic with synthetic sea surface height (SSH) data, this RKHS-based filter accurately assimilates observations with non-Gaussian noise, improving both the accuracy and efficiency of ocean state estimations in complex systems.

32

**Plain Language Summary**

33

34

35

36

37

Estimating the ocean’s dynamic state, such as eddies and currents, is challenging due to the complexity of ocean processes. This paper introduces a new data assimilation method, grounded in a rigorous theoretical framework, that uses an ensemble of simulations to improve these estimates. It builds new trajectories by combining ensemble members in a time-invariant way, enabling efficient assimilation of observations across space and time.

38

39

40

41

42

By updating the ensemble at the start of each observation period and simulating forward, the method enhances physical consistency. Tests with a North Atlantic model and synthetic satellite data demonstrate its ability to accurately estimate the ocean’s state, even with noisy, non-Gaussian observations, while improving both accuracy and efficiency compared to traditional methods.

43

**1 Introduction**

44

45

46

47

48

49

50

51

52

Data assimilation methods in ocean modeling aim to integrate sparse observations, both in time and space, of the ocean surface into numerical circulation models to produce accurate estimates of the ocean's state. In this context, a major challenge is the assimilation of sea surface height (SSH) observations, collected over time from orbiting satellite swaths, to accurately represent ocean mesoscale dynamics. The mesoscale, characterized by features such as eddies, fronts, and jets, is crucial for understanding and modeling global ocean circulation. To address this complex problem, several methods proposed by the data assimilation community (Carrassi et al., 2018; Evensen, 2009; Reich & Cotter, 2015) have been adapted and are now routinely used for various applications.

53

54

55

56

57

58

59

60

61

62

63

Optimal interpolation (OI) is a geostatistical method used to reconstruct a dense field of relevant oceanic variables from spatially sparse satellite observations by considering predetermined statistical properties (Le Traon et al., 2003). This method integrates observations collected at various times, assuming spatiotemporal correlation function between the observations and the estimated field. It is in particular employed to construct the Data Unification and Altimeter Combination System (DUACS) products (Taburet et al., 2019). OI's reliance on statistical interpolation to merge observations from varying times often results in an over-smooth degraded representation of small-scale structures. This limitation is particularly evident in the reconstruction of mesoscale features, where finer details are crucial for accurate modeling. Ubelmann et al. (2015) observed that optimal interpolation often misrepresents the transport over time of small-scale eddies. To address this issue, they
introduced dynamic interpolation, which offers a more accurate reconstruction of small-scale structure transport between two observations. Assuming an invertible dynamical system, the method involves running the model forward and backward between two observation times, and then reconstructing the sea surface height as the mean of the results from both simulations. Although this method better represents transport in the reconstructions, it still requires two space interpolations to obtain the complete fields needed to initialize the forward and backward simulations. Consequently, the effective resolution remains relatively low.

Variational methods aim to minimize a cost function that penalizes the distance between the observations and the estimated state. Among these, four-dimensional variational assimilation (4D-Var), assimilates spatially and temporally distributed observations. This method uses optimizations tools, such as adjoint methods (Le Dimet & Talagrand, 1986), which efficiently compute the gradient of the cost function with respect to its parameters by running the dynamics forward and then the tangent linear adjoint system backward. Despite the effectiveness of these methods, the associated computational costs can be prohibitively high, hardly parallelizable, and constructing a numerically exact tangent linear operator can be a challenging task. To address this, neural networks such as 4DVarNet (Fablet et al., 2021) can be use to efficiently solve it after a learning phase. Recently, 4DVarNet has been applied to the reconstruction of the sea surface currents from sea surface height (SSH) and sea surface temperature (SST) observations (Fablet et al., 2024), demonstrating reconstructions with a higher effective resolution than DUACS. Machine learning methods offer powerful tools for efficient data assimilation (Cheng et al., 2023). These approaches can be integrated with classical assimilation methods to improve future state predictions (Brajard et al., 2020; Farchi et al., 2021). It can also be used to interpolate partial observations to retrieve the ocean mesoscale dynamic state (Ouala et al., 2018; Manucharyan et al., 2021). In these methods, neural networks are employed to model the system's dynamics, thereby bypassing the need to solve the dynamical system equations directly. However, this class of approach opens the question how to assess, without the help of a dynamical model, the physical quality of the generated predictions. In particular, neural networks introduce generalization challenges, which may lead to poor performance when the model encounters new or previously unseen physical phenomena.

Another strategy, called nudging, introduces a correction term into the dynamics equations to fit better with the observations (Hoke & Anthes, 1976). Back-and-forth nudging (BFN) (Auroux & Blum, 2008), involves iteratively solving the forced dynamical equation forward and backward in time until the algorithm converges toward a trajectory. This method circumvents the need for minimization tools used in 4D-Var methods, but it requires uncontrolled modifications of the dynamical system, and hence of the physical prior we have on the system's dynamics. Le Guillou, Metref, et al. (2021) apply BFN to data from the Surface Water Ocean Topography (SWOT) mission to assimilate a quasi-geostrophic model, demonstrating the potential of SWOT data to enhance BFN performance. It can be remarked, that BFN incurs computational costs due to the necessity of running the model multiple times until convergence. This makes it impractical as it is for application to complex models. To address this limitation, Le Guillou, Lahaye, et al. (2021) proposed a hybrid approach that combines BFN with 4D-Var. In this approach, BFN is used to assimilate balanced motions in a quasi-geostrophic model, while 4D-Var is employed to assimilate internal tides in a shallow water model.

To manage computational complexity in high-dimensional ocean models, one effective approach is the use of ensemble Kalman filter (EnKF) methods, which operate locally in the phase space (Evensen, 2009; Currassi et al., 2018), in the sense that the solution is assumed to be well characterized by the span of an ensemble of realizations of the dynamical model relatively close from each other. The use of Kalman filters in nonlinear systems is justified by their ability to provide the best linear unbiased estimator within a Gaussian framework. EnKFs assimilate observations sequentially as they become available, assuming
117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

both the state estimation error and observation error to be Gaussian. They approximate
the state estimate and its uncertainty using an ensemble of realizations, where the analyzed
ensemble members are constructed as linear combinations of the forecast ensemble. However,
this superposition principle lacks theoretical justification for highly nonlinear systems and
it has been shown that EnKF do not converge toward the true filtering distribution (Le
Gland et al., 2011). For a general review on the asymptotic results of the EnKF in a
Gaussian context, interested readers may refer to (A. N. Bishop & Del Moral, 2023). Despite
this limitation, EnKFs have been successfully applied in various oceanographic studies,
demonstrating their efficiency in estimating the ocean state by assimilating SST or SSH data
(Leeuwenburgh et al., 2005; Wan et al., 2010). A significant limitation of classical EnKFs
is their sequential nature, which leads to correct the model at the final observation time of
a given time window (and not at the window's initial time). To address this, smoothing
methods based on forward-backward assimilation of trajectories from data collected over the
whole time window have been developed. Among these, ensemble extension of the 4DVar
techniques (4DEnVar) have been proposed (Lorenc, 2003; Liu et al., 2008; Desroziers et al.,
2014). Such techniques replace the exact numerical tangent linear operator by an ensemble
approximation, leading to much cheaper and practical methods (Yang et al., 2015). Other
methods use time-constant weights to reconstruct trajectories, and formulate a cost function
of this weight vector (Hunt et al., 2004; Bocquet & Sakov, 2014). Connections between these
methods and 4DEnVar can be established. However, the assumption of constant weights
over time lacks clear theoretical justification in a classical setting.

138

139

140

141

142

143

144

145

Alternatively, particle filters (Doucet et al., 2001; Gordon et al., 1993; Van Leeuwen et al., 2019) relax the Gaussianity and linearity assumptions of Kalman filters by representing the state distribution with an ensemble of particles, each associated with a weight. Unlike EnKFs, particle filters do not alter the ensemble's physical properties by creating members as linear combinations of existing members. However, particle filters are subject to the curse of dimensionality, meaning that in high-dimensional state spaces, a huge number of particles is required to accurately represent the state distribution. This makes particle filters computationally expensive and challenging to apply to high-dimensional systems.

146

147

148

149

150

151

152

153

154

155

156

157

Reproducing Kernel Hilbert Spaces (RKHS) have been employed in data assimilation
methods to address system non-linearity and to relax the inherent assumptions of Kalman
filters. Mauran et al. (2023) propose a reformulation of the Ensemble Transform Kalman
Filter (ETKF) (C. H. Bishop et al., 2001; Hunt et al., 2007), operating within a RKHS
framework that uses feature maps (i.e. the mapping from the phase space onto the RKHS)
defined from the ensemble state space and observations. Gottwald and Reich (2021) intro-
duced an ensemble Kalman filter designed to simultaneously retrieve the state variable and
learn the dynamical propagator, using random feature map bases. Both approaches aim
to determine weight vectors whose size scales with the state space, making them compu-
tationally expensive for high-dimensional systems. This issue of increasing computational
cost with system dimensionality is also observed in other kernel-based assimilation methods
(Luo, 2019; Zhang et al., 2022).

158

159

160

161

162

163

RKHS-based assimilation methods share similarities with the decomposition method introduced by Ubelmann et al. (2021), where the state variable of the dynamical system is decomposed into a reduced space, and optimal interpolation is performed directly on the reduced state vector to define the analyzed state. In RKHS approaches, the reduced state vector corresponds to the weight vector in the RKHS, and the ensemble filter operates directly in this space.

164

165

166

167

168

169

The Koopman operator (Koopman, 1931) is an infinite-dimensional linear operator
that governs the evolution of observables in a dynamical system. It enables the analysis
of complex, nonlinear systems using the tools of linear operator theory, offering a powerful
framework for understanding and predicting dynamical behavior. Since the seminal work of
Mezić (Mezić, 2005), numerous methods have been developed to approximate the spectral
decomposition of the Koopman operator, enabling the identification of Koopman modes
170

171

172

173

174

175

and eigenfunctions (Budišić et al., 2012; Williams, Kevrekidis, & Rowley, 2015; Kutz et al., 2016). These techniques have been extensively applied in various domains, particularly in fluid mechanics (Rowley et al., 2009; Schmid, 2010; Tu et al., 2014), where they offer insights into coherent structures and modal dynamics. Moreover, some advancements have integrated Koopman operator theory with kernel methods to enhance data-driven approximations (Williams, Rowley, & Kevrekidis, 2015; Das & Giannakis, 2020).

176

177

178

179

180

181

Such a decomposition has also been used for analyzing ocean circulation or climate from satellite data (Hogg et al., 2020; Navarra et al., 2021; Zhen et al., 2022). In particular, this has been performed by relying on the spectral representation of the time series autocovariance (Ghil et al., 2002; Kondrashov et al., 2020). It has been proven in (Zhen et al., 2022) that the leading eigenpairs of the autocovariance matrix correspond to Koopman eigenpairs for ergodic systems with a finite invariant measure.

182

183

184

185

186

187

188

189

190

Let us point out that the study (Navarra et al., 2021), which relies on the approach of (Klus et al., 2020), embeds the estimation of the Koopman spectral decomposition within a given RKHS. This is performed under the restrictive assumption that the RKHS is invariant under the Koopman operator, which would require defining the RKHS kernel from the system's dynamics rather than using a fixed kernel. As shown in the following, we will consider such a dynamical construction, leading to a family of RKHS, together with an ensemble of runs instead of long time series, and a compactification technique applied to the unbounded infinitesimal generator of the Koopman operator, which also results in dealing with a family of RKHS (Das & Giannakis, 2020).

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

The method proposed in this study relies on the theoretical framework developed in (Dufée et al., 2024), which differs from usual RKHS based estimation. In this setting the dynamical system is embedded within a family of time-varying reproducing kernel Hilbert spaces (RKHS). This framework enables establishing an evolution law for the RKHS bases across time. The RKHS are transported along time by the Koopman operator, defining an isometry between the kernel evaluations over time. This approach enables the use of embarrassingly simple data assimilation schemes, reconstructing the trajectory as a time-invariant linear combination of the ensemble members. These schemes rely on a fully justified superposition principle allowing to compare the trajectories at initial time, leading consequently to effectively solve a 4D-Var problem with a time series of observations at the computational cost of a three-dimensional variational (3D-Var) methods, in which time evolution is not considered. Compared to the other RKHS schemes proposed in the literature our method works on a small dimensional vector of function evaluations corresponding to the ensemble dimension. The small dimension of these evaluation vectors allows for direct matrix inversion in the weight space making our approach applicable to high-dimensional systems. Similar to particle filters, each ensemble member in the RKHS are assigned with (a signed) weight, representing the state distribution.

208

209

210

211

This method has been applied to a multilayer quasi-geostrophic system in the North Atlantic basin, using synthetic observations generated from real-world satellite tracks, following an observing system simulation experiment (OSSE) methodology. We show that, within this framework, the proposed method leads to improved reconstruction resolution.

212

213

214

215

The paper is structured as follows: section 2 provides a brief overview of the RKHSs, the Koopman operator and how with the combination of these two tools a RKHS family transported by the dynamics is defined. Section 3 focuses on assimilation methods. Section 4 introduces the experimental framework. Finally, section 5 describes the results obtained.

216

## 2 Theoretical framework

217

218

219

This section introduces the RKHS theory, the Koopman operator, and the framework established by Dufée et al. (2024), which defines a family of RKHSs, built from the transport by the dynamical system of an initial kernel at the time instant $t_0$.
## 2.1 RKHS

A Reproducing Kernel Hilbert Space (RKHS), $\mathcal{H}$, is a Hilbert space of functions $f: E \to \mathbb{C}$, where $E$ is a non-empty set, equipped with an inner product $\langle \cdot, \cdot \rangle_{\mathcal{H}}$. A map $k: E, E \to \mathbb{C}$ is defined as being a reproducing kernel of $\mathcal{H}$ if it has the two following properties:

* membership property: $\forall x \in E, k(\cdot, x) \in \mathcal{H}$,

* reproducing property: $\forall x \in E, \forall f \in \mathcal{H}, \langle f, k(\cdot, x) \rangle_{\mathcal{H}} = f(x)$.

The reproducing property implies that the kernel satisfies $k(x, y) = \overline{k(y, x)}$, where $\cdot$ denotes the complex conjugate. The RKHS is fully and uniquely determined by its kernel $k$ (Aronszajn, 1950). The functions $k(\cdot, x)$ for $x \in E$ are usually referred to as feature maps. It can be shown that their span is dense in $(\mathcal{H}, \|\cdot\|_{\mathcal{H}})$. The integral operator associated with the kernel $k$, is defined as $\mathcal{L}_k: L^2(E, \nu) \to L^2(E, \nu)$, such that for any $f \in L^2(E, \nu)$,

$$ (\mathcal{L}_k f)(x) = \int_E k(x, y) f(y) \nu(dy), \quad (1) $$

where $L^2$ denotes the space of square-integrable functions, and $\nu$ is a measure on $E$. The inner product in $L^2(E, \nu)$ for $f$ and $g$ in $L^2(E, \nu)$ is given by

$$ \langle f, g \rangle_{L^2(E, \nu)} = \int_E f(y) \overline{g(y)} \nu(dy). \quad (2) $$

Consequently, the inner product between a function $f \in L^2(E, \nu)$ and a kernel function $k(\cdot, x)$ becomes:

$$ \langle f, k(\cdot, x) \rangle_{L^2(E, \nu)} = \int_E f(y) \overline{k(y, x)} \nu(dy) = \int_E k(x, y) f(y) \nu(dy) = (\mathcal{L}_k f)(x). \quad (3) $$

The range of $\mathcal{L}_k$ is assumed to be dense in $L^2(E, \nu)$, implying that its eigenvalues are strictly positive. By Mercer's theorem as $k$ is a bounded, continuous, symmetric, positive semi-definite kernel, $\mathcal{L}_k$ eigenfunctions $\{\varphi_i\}_{i \in \mathbb{N}}$ and eigenvalues $\{\mu_i\}_{i \in \mathbb{N}}$ define the RKHS $\mathcal{H}$ as

$$ \mathcal{H} = \left\{ f \in L^2(E, \nu), f = \sum_{i=1}^{\infty} a_i \varphi_i \quad \text{s.t.} \quad \sum_i \frac{|a_i|^2}{\mu_i} < \infty \right\}. \quad (4) $$

This last expression highlights the fact that $\mathcal{H}$ is a space of functions with fast decaying coefficients within $L^2(E, \nu)$.

## 2.2 Dynamical system

Let consider a dynamical system defined by a continuous flow $\Phi_t$, operating on functions $\mathbf{X} \in L^2(\mathbb{R}^+ \times \Omega_x, \mathbb{R}^d)$ included in an invariant phase space manifold, with $\Omega_x$ representing the spatial support. The state of the system at time $t \in \mathbb{R}^+$, denoted by $\mathbf{X}(t) = \mathbf{X}_t$ is given by $\mathbf{X}_t = \Phi_t(\mathbf{X}_0)$, where $\mathbf{X}_0$ is the initial condition. In differential form the flow is given by the system:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{X}_t}{\partial t} = \mathcal{M}(\mathbf{X}_t), \\ \mathbf{X}(0) = \mathbf{X}_0 \end{array} \right. , \quad (5) $$

where $\mathbf{X}_t$ belongs to $\Omega$ a compact (invariant) subset of the space of square-integrable functions $L^2(\Omega_x, \mathbb{R}^d)$. The nonlinear operator $\mathcal{M}$ is assumed to be differentiable with continuous derivatives. An observable of this system is defined as a square-integrable function that takes as input a function from $\Omega$ and outputs a real valued vector. In the following, we assume that the dynamical system is invertible, and that there exists a measure $\nu$ preserved by this dynamical system.
**Remark 1** Let us point out that in our data assimilation setting, the invertibility condition will be required only on the assimilation time window. Existence of invariant measures for continuous system is ubiquitous.

The Koopman operator $U_t$ is an infinite-dimensional linear operator (Koopman, 1931), that acts on the space of square-integrable functions $L_\mathbb{C}^2(\Omega, \nu)$, mapping it into itself. It propagates the observables over time, such that for every observable $f$ and every $\mathbf{X}$ in $\Omega$, we have $U_t f(\mathbf{X}) = f(\Phi_t(\mathbf{X}))$. The Koopman operator is a powerful tool used to study non linear dynamical systems. It has found extensive applications across various scientific and engineering disciplines (Brunton et al., 2022; Eisner et al., 2015; Mezić, 2021).

## 2.3 Dynamical system embedded in a RKHS

With all the previous ingredients, we are now able to present the framework proposed by Dufée et al. (2024), on which our assimilation schemes are built.

Let $\Omega_0$ denote the ensemble of initial conditions, and define $\Omega_t := \Phi_t(\Omega_0)$ as the ensemble of states (or snapshots) at time $t \in \mathbb{R}^+$. Since the system is invertible, for any snapshot $\mathbf{X}_t \in \Omega_t$, there exists a unique initial condition $\mathbf{X}_0 \in \Omega_0$ associated with it. Moreover in the present framework, we consider that for each state $\mathbf{X}_t$, the current time $t$ is known, in such a way it can be uniquely associated with its initial condition $\mathbf{X}_0$. Given an initial kernel $k_0: \Omega_0 \times \Omega_0$ it is then possible to define a positive Hermitian kernel $k_t: \Omega_t \times \Omega_t$ such that,

$$ \forall \mathbf{X}_t, \mathbf{Y}_t \in \Omega_t, k_t(\mathbf{X}_t, \mathbf{Y}_t) = k_0(\Phi_t^{-1}(\mathbf{X}_t), \Phi_t^{-1}(\mathbf{Y}_t)). \quad (6) $$

Let $(\mathcal{H}_0, \langle \cdot, \cdot \rangle_{\mathcal{H}_0})$ and $(\mathcal{H}_t, \langle \cdot, \cdot \rangle_{\mathcal{H}_t})$ represent the RKHSs associated with $k_0$ and $k_t$, respectively. A RKHS family $\mathcal{W} = (\mathcal{H}_t)_{t \ge 0}$ can be constructed by considering the RKHSs along time.

The operator $\mathcal{U}_t: \mathcal{H}_0 \to \mathcal{H}_t$ is defined as:

$$ \forall \mathbf{X}_0 \in \Omega_0, \mathcal{U}_t k_0(\cdot, \mathbf{X}_0) = k_t(\cdot, \Phi_t(\mathbf{X}_0)). \quad (7) $$

This operator propagates the features maps over time. It can be seen as a useful tool to define the kernel $k_t$ at any time instant. Operator $\mathcal{U}_t$ is closely related to the adjoint of the Koopman operator $U_t$, known as the Perron Frobenius operator. It follows from equation (6) that $\mathcal{U}_t$ defines an isometry between $\mathcal{H}_0$ and $\mathcal{H}_t$.

The theory presented here operates on infinite ensembles of functions, given that both $\Omega$ and $\mathcal{H}$ are infinite dimensional spaces. However, in practice, we are constrained to consider only a finite number of trajectories, denoted as $\{\mathbf{X}^{(i)}, i = 1, \dots, p\}$. Figure 1 illustrates the transported feature maps associated with the ensemble members, which will be employed to solve the data assimilation problem in section 3.2.

In ensemble data assimilation methods, a superposition principle is most of the time assumed, which in theory is not applicable for non-linear dynamical systems. However, this principle can be rigorously justified within the RKHS family described previously (Dufée et al., 2024). Considering an observable $g_x: \Omega_t \to \mathbb{R}^d$ such that for $t \ge 0$, $\forall \mathbf{X}_t \in \Omega_t, \forall x \in \Omega_x$, $g_x(\mathbf{X}_t) = \mathbf{X}_t(x)$. Assuming $g_x$ is smooth enough to belong to $\mathcal{H}_t$, the kernel reproducing property allows the following decomposition:

$$ g_x(\mathbf{X}_t) = \langle g_x(\cdot), k_t(\cdot, \mathbf{X}_t) \rangle_{\mathcal{H}_t}. \quad (8) $$

Since $k_t(\cdot, \mathbf{X}_t)$ is included in $\mathcal{H}_t$ it admits a decomposition such that

$$ k_t(\cdot, \mathbf{X}_t) = \sum_{i=1}^{p} \beta_i k_t(\cdot, \mathbf{X}_t^{(i)}), \quad (9) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Feature maps transported over time from $t_0$ to $t_k$ by the dynamical system, with the observations collected between $t_0$ and $t_k$ represented at the center of the figure in the form of satellite swath.</figcaption></figure>

where $\beta$ is a vector of length $p$. At the continuous level, the sum involved in (9) is infinite, however with the finite number of ensemble members considered, it reduces to a finite sum approximation. By the inner product bilinearity and the kernel reproducing property it follows that:

$$g_x(\mathbf{X}_t) = \sum_{i=1}^{p} \beta_i g_x(\mathbf{X}_t^{(i)}). \qquad (10)$$

As this holds for every $x \in \Omega_x$, it can be concluded that:

$$\mathbf{X}_t = \sum_{i=1}^{p} \beta_i \mathbf{X}_t^{(i)}, \qquad (11)$$

at condition that all the observables constituting the state vector – or at least the variables we aim at estimating – belong to $\mathcal{H}_t$. This superposition principle allows the construction of new trajectories as linear combinations of previously computed ones, with the linear combination being constant over time due to equation (6). As further described in the following sections this property leads to the setting of very appealing data assimilation methods.

**Remark 2** *From a machine learning perspective, the RKHS family can be viewed as latent spaces defined from the dynamics, where the weights in the latent spaces are equal to those in the physical space, and characterize a given trajectory.*

In the next section we provide a brief recap on ensemble data assimilation techniques.

# 3 Assimilation methods

## 3.1 Classical methods

Although the previous setting is expressed in full generality within functional spaces, for the sake of simplicity, we will describe the assimilation techniques in the following sections using finite-dimensional spaces.

Assimilation methods aim to provide the best estimate of the state of a dynamical system from observations. The observations $\mathbf{Y} \in \mathbb{R}^m$ are assumed to be related to the state of the system $\mathbf{X} \in \mathbb{R}^n$, as follows,

$$\mathbf{Y} = H(\mathbf{X}) + \epsilon, \qquad (12)$$

where $H: \mathbb{R}^m \to \mathbb{R}^n$, with usually $n \gg m$, denotes the observation operator, $\epsilon \in \mathbb{R}^m$ represents an additive centred Gaussian observation noise with covariance matrix $\mathbf{R} \in \mathbb{R}^{m \times m}$.
329

330

331

332

333

334

Classical assimilation methods aim to solve the ill-posed problem $H(\mathbf{X}) \simeq \mathbf{Y}$, seeking to find the state that best fits the observations. For sake of simplicity, $H$ will be assumed to be linear in the following, and thus represented through matrix $\mathbb{H} \in \mathbb{R}^{m \times n}$ such that $H(\mathbf{X}) = \mathbb{H}\mathbf{X}$. To have a well-posed problem and therefore a unique solution, one searches for the closest solution to the forecasted state $\mathbf{X}^f$. Consequently, assimilation methods generally seek to minimize a cost function as follows

335

$$J(\mathbf{X}) = \frac{1}{2} \|\mathbb{H}\mathbf{X} - \mathbf{Y}\|_{\mathbf{R}^{-1}}^2 + \frac{1}{2} \|\mathbf{X} - \mathbf{X}^f\|_{\mathbf{P}^{-1}}^2, \quad (13)$$

336

337

where $\mathbf{P} \in \mathbb{R}^{n \times n}$ is the state variable covariance matrix, and as usual $\|\bullet\|_{\mathbf{A}}$ stands for the norm $(\bullet^{\mathrm{T}}\mathbf{A}\bullet)$, with upperscript, $T$, denoting the transpose operator.

338

339

Cancellation of the cost function gradient gives the following expression of the optimal analyzed state

340

$$\mathbf{X}^a = -(\nabla^2 J(\mathbf{X}))^{-1} (\nabla J(0)) = (\mathbb{H}^T \mathbf{R}^{-1} \mathbb{H} + \mathbf{P}^{-1})^{-1} (\mathbb{H}^T \mathbf{R}^{-1} \mathbf{Y} + \mathbf{P}^{-1} \mathbf{X}^f). \quad (14)$$

341

Using the Sherman-Morrison-Woodbury formula, the previous equation can be rewritten,

342

$$\mathbf{X}^a = \mathbf{X}^f - \mathbf{G} (\mathbb{H}\mathbf{X}^f - \mathbf{Y}), \quad (15)$$

343

with the Kalman gain expressed as

344

$$\mathbf{G} = \mathbf{P}\mathbb{H}^T (\mathbf{R} + \mathbb{H}\mathbf{P}\mathbb{H}^T)^{-1}. \quad (16)$$

345

### 3.1.1 Ensemble Square root Filter (ESRF)

346

347

348

349

Ensemble methods consider $p$ realizations (or members) in the state space, whose mean and covariance matrix represent the estimated state and covariance. Note that in high-dimensional problems, such as those encountered in ocean dynamics, we have $p \ll m \ll n$, with typical values of $p \sim 50$, $m \sim 10^3$, $n \sim 10^6$.

350

351

352

353

354

355

356

The ensemble matrix $\mathbf{X}$ is defined as a $n \times p$ matrix. Each columns $i$ of $\mathbf{X}$ contains a member of the ensemble noted $\mathbf{X}^{(i)}$, with $1 \le i \le p$. We briefly present in the following the ensemble square root filter (Whitaker & Hamill, 2002), which will serve as a reference in this study. Its name comes from the fact that the state covariance matrix is expressed through its square root, the anomaly matrix $\mathbf{A} = [(\mathbf{X}^{(1)} - \bar{\mathbf{X}}), (\mathbf{X}^{(2)} - \bar{\mathbf{X}}), \dots, (\mathbf{X}^{(p)} - \bar{\mathbf{X}})]$, where $\bar{\mathbf{X}}$ denotes the ensemble mean, which represents the estimated state. The empirical covariance matrix is then expressed as:

357

$$\mathbf{P} = \frac{1}{p-1} \mathbf{A} \mathbf{A}^T. \quad (17)$$

358

With this decomposition, $\mathbf{A} \in \mathbb{R}^{n \times p}$ is stored instead of $\mathbf{P} \in \mathbb{R}^{n \times n}$.

359

360

361

The analyzed covariance matrix is expressed as the inverse of the Hessian of cost function $J$ defined in equation (13). Using the Sherman-Morrison-Woodbury formula the covariance matrix update reads

362

$$\mathbf{P}^a = (\mathbf{I} - \mathbf{G}\mathbb{H})\mathbf{P}^f, \quad (18)$$

363

364

365

where $\mathbf{G}$ is the gain matrix expressed in equation (16). Exponents $\bullet^f$ and $\bullet^a$ stand for forecast and analyzed, respectively. Substituting the gain into the covariance matrix formula leads to,

366

$$\frac{1}{p-1} \mathbf{A}^a (\mathbf{A}^a)^T = \mathbf{A}^f \left( \mathbf{I} - \frac{1}{p-1} (\mathbb{H}\mathbf{A}^f)^T \left( \mathbf{R} + \mathbb{H} \frac{1}{p-1} \mathbf{A}^f (\mathbf{A}^f)^T \mathbb{H}^T \right)^{-1} (\mathbb{H}\mathbf{A}^f) \right) \frac{1}{p-1} (\mathbf{A}^f)^T. \quad (19)$$
367

368

Applying the Sherman–Morrison–Woodbury formula the previous equation can be written as

$$ \mathbf{A}^a (\mathbf{A}^a)^T = \mathbf{A}^f \left( \mathbf{I} + \frac{1}{p-1} (\mathbb{H}\mathbf{A}^f)^T \mathbf{R}^{-1} (\mathbb{H}\mathbf{A}^f) \right)^{-1} (\mathbf{A}^f)^T . \quad (20) $$

369

370

Using the square-root matrix $\mathbf{S} \in \mathbb{R}^{p \times p}$ defined as

$$ \mathbf{S} = \left( \mathbf{I} + \frac{1}{p-1} (\mathbb{H}\mathbf{A}^f)^T \mathbf{R}^{-1} (\mathbb{H}\mathbf{A}^f) \right)^{-1/2} . \quad (21) $$

371

372

The analyzed anomaly is expressed from the forecasted anomaly, as:

$$ \mathbf{A}^a = \mathbf{A}^f \mathbf{S} . \quad (22) $$

373

374

The analyzed ensemble mean is defined using equation (15),

$$ \overline{\mathbf{X}}^a = \overline{\mathbf{X}}^f - \mathbf{G} (\overline{\mathbb{H}\mathbf{X}}^f - \mathbf{Y}) , \quad (23) $$

375

376

377

and the analyzed ensemble members are obtained by adding to the analyzed mean (23) the analyzed anomaly (22),

$$ \mathbf{X}^a = \overline{\mathbf{X}}^a + \mathbf{A}^a \quad (24) $$

$$ = \overline{\mathbf{X}}^a + \mathbf{A}^f \mathbf{S} . \quad (25) $$

378

379

380

It is important to note that (25) defines the update of the ensemble in the state space. In ocean models, as indicated previously the state is typically very large. To address this issue it is preferable to work in the weight space.

381

### 3.1.2 ESRF in the weight space

382

383

384

385

386

387

388

We now express directly the ESRF in the weight space. To that end, let $\beta^f$ and $\beta^a$ denote the $p \times p$ weight matrices representing the forecast and analyzed ensembles, such that $\mathbf{X}^f = \mathbf{X}^f \beta^f$ and $\mathbf{X}^a = \mathbf{X}^f \beta^a$. The $i^{th}$ column of the weight matrix contains the coefficients to construct the $i^{th}$ member of the ensemble from the forecasted ensemble members. It appears clearly from its definition, that $\beta^f$ is the identity matrix. This formulation is similar to decompositions found in the literature, such as in Lorenc (2003), Hunt et al. (2004), and Liu et al. (2008).

389

Injecting $\mathbf{P}^f = \mathbf{A}^f (\mathbf{A}^f)^T$ in the gain definition leads to

$$ \mathbf{G} = \frac{1}{p-1} \mathbf{A}^f (\mathbf{A}^f)^T \mathbb{H}^T \left( \mathbf{R} + \mathbb{H} \frac{1}{p-1} \mathbf{A}^f (\mathbf{A}^f)^T \mathbb{H}^T \right)^{-1} . \quad (26) $$

390

391

392

Writing $\mathbf{A} = \mathbf{X} (\mathbf{I} - \mathbf{1}/p)$, where $\mathbf{1}$ denotes the $p \times p$ matrix filled with one, the previous equation reads then

$$ \mathbf{G} = \frac{1}{p-1} \mathbf{X}^f (\mathbf{I} - \mathbf{1}/p) (\mathbf{I} - \mathbf{1}/p)^T (\mathbf{X}^f)^T \mathbb{H}^T \left( \mathbf{R} + \mathbb{H} \frac{1}{p-1} \mathbf{A}^f (\mathbf{A}^f)^T \mathbb{H}^T \right)^{-1} . \quad (27) $$

393

394

395

A gain matrix in the weight space, denoted $\mathbf{G}_w \in \mathbb{R}^{p \times m}$, can thus be defined such that $\mathbf{G} = \mathbf{X}^f \mathbf{G}_w$, with,

$$ \mathbf{G}_w = \frac{1}{p-1} (\mathbb{H}\mathbf{A})^T \left( \mathbf{R} + \frac{1}{p-1} \mathbb{H}\mathbf{A}^f (\mathbb{H}\mathbf{A}^f)^T \right)^{-1} . \quad (28) $$

396

397

398

Considering the mean weight vectors $\overline{\beta}^f$ and $\overline{\beta}^a$, such that $\overline{\mathbf{X}}^f = \mathbf{X}^f \overline{\beta}^f$ and $\overline{\mathbf{X}}^a = \mathbf{X}^f \overline{\beta}^a$, which correspond to the column averages of the associated weight matrices. As $\beta^f$ is the
identity, $\overline{\beta^f}$ is a $p$ length vector filled with $\frac{1}{p}$. The equation (23) can thus be rewritten in the weight space:

$$ \overline{\beta^a} = \overline{\beta^f} - G_w (\mathbb{H}X^f\overline{\beta^f} - Y). \qquad (29) $$

As $X^a = X^f\beta^a$ and $A^f = X^f(\beta^f)'$, where $\cdot'$ denotes the right matrix multiplication of the operand by the matrix $I - 1/p$, equation (25) can be rewritten

$$ X^f \beta^a = X^f \overline{\beta^a} + X^f (\beta^f)' S, \qquad (30) $$

which yields the weight update

$$ \beta^a = \overline{\beta^a} + (\beta^f)' S. \qquad (31) $$

## 3.2 Ensemble filter on RKHS family

The formalism proposed by Dufée et al. (2024), as detailed in section 2, has the advantage to be directly formulated in terms of piece-wise in time, linear combination of trajectories. It naturally enables the estimation of an optimal trajectory from time-series of observations. More precisely, this ensemble filter (EnF) defined on the RKHS family, operates on the weights, within the RKHS family, assuming the observation operator is well defined on the RKHS family through the smooth evaluation operator assumption (8). In the following, we will assume without loss of generality that the kernel is real-valued. The filter, referred to as RKHS-EnF in the following, is built from a 4DVar-like cost function, accounting for time-series of observations, and incorporating a regularization term based on the norm of the kernel function in the initial RKHS $\mathcal{H}_0$:

$$ J(\mathbf{X}) = \frac{1}{2} \sum_{j=0}^{J} \left\| \mathbb{H}_{t_j} \mathbf{X}_{t_j} - \mathbf{Y}_{t_j} \right\|_{\mathbf{R}_{t_j}^{-1}}^2 + \frac{1}{2} \left\| k_0(\cdot, \mathbf{X}_0 - \overline{\mathbf{X}}_0^f) \right\|_{\mathcal{H}_0}^2, \qquad (32) $$

where $\mathbb{H}_{t_j}, \mathbf{X}_{t_j}, \mathbf{Y}_{t_j}$ stand respectively for the observation matrix, the ensemble matrix and the observation vector at time $t_j$. Given $\beta$ the weight vector such that $\mathbf{X}_{t_j} = \sum_{i=1}^p \beta_i \mathbf{X}_{t_j}^{(i)}$, and using equations (9), the second term of the right-hand side can be expressed as

$$ \|k_0(\cdot, \mathbf{X}_0 - \overline{\mathbf{X}}_0)\|_{\mathcal{H}_0}^2 = \left\| \sum_{i=1}^p (\beta_i - \overline{\beta}_i^f) k_0(\cdot, \mathbf{X}^{(i)}) \right\|_{\mathcal{H}_0}^2. \qquad (33) $$

By the kernel reproducing property, we obtain

$$ \|k_0(\cdot, \mathbf{X}_0 - \overline{\mathbf{X}}_0)\|_{\mathcal{H}_0}^2 = (\beta - \overline{\beta}^f)^T \mathbf{K} (\beta - \overline{\beta}^f), \qquad (34) $$

where $\mathbf{K}$ is the kernel matrix that contains the evaluation of the initial kernel between every member's pair of the ensemble, $\mathbf{K}_{i,j} = k_0(\mathbf{X}_0^{(i)}, \mathbf{X}_0^{(j)})$, for $1 \le i, j \le p$. The second member of the cost function penalizes hence the weight anomaly in the RKHS.

The previous cost function can thus be rewritten as a function of the weight vector in the RKHS, thereby formulating the assimilation problem directly in the weight space:

$$ J_{\text{RKHS}}(\boldsymbol{\beta}) = \frac{1}{2} \sum_{j=0}^{J} \left\| \sum_{i=1}^{p} \beta_i \mathbb{H}_{t_j} \mathbf{X}_{t_j}^{(i)} - \mathbf{Y}_{t_j} \right\|_{\mathbf{R}_{t_j}^{-1}}^2 + \frac{1}{2} \left\| \boldsymbol{\beta} - \overline{\boldsymbol{\beta}^f} \right\|_{\mathbf{K}}^2. \qquad (35) $$

Let us define $\mathbf{V}$ such that $\mathbb{H}_t \mathbf{X}_t = \mathbf{V}_t \mathbf{k}$ and $\mathbf{V} = (\mathbf{V}_{t_j})_{j \in [1, J]}$ where $\mathbf{k} = (k_0(\mathbf{X}_0, \mathbf{X}^{(i)}))_{i \in [1, p]}$.

By the reproducing property, we have $\mathbf{k} = \mathbf{K}\boldsymbol{\beta}$. Denoting $\mathbf{Y} = (\mathbf{Y}_{t_j})_{j \in [1, J]}$ and $\tilde{\mathbf{R}} = \mathbf{I}_J \otimes (\mathbf{R}_{t_j})_{j \in [0, J]}$, where $\otimes$ represents the Kronecker product, the cost function can be finally compactly expressed as:

$$ J_{\text{RKHS}}(\boldsymbol{\beta}) = \frac{1}{2} (\mathbf{V}\mathbf{K}\boldsymbol{\beta} - \mathbf{Y})^T \tilde{\mathbf{R}}^{-1} (\mathbf{V}\mathbf{K}\boldsymbol{\beta} - \mathbf{Y}) + \frac{1}{2} (\boldsymbol{\beta} - \overline{\boldsymbol{\beta}^f})^T \mathbf{K} (\boldsymbol{\beta} - \overline{\boldsymbol{\beta}^f}) . \qquad (36) $$
**Remark 3** In equation (35), a linear observation operator is considered through the matrix $\mathbb{H}_{t_j}$. However, the expression (36) remains valid for non-linear observation operators, at the condition that the observables associated with these observations belong to $\mathcal{H}_t$. This is the reason of introducing explicitly the notation $\mathbb{H}_t X_t = V_t K \beta$. This linear combination representation is a direct consequence of the reproducing property. The fact that these coefficients are constant in time ensues from the framework proposed by Dufée et al. (2024) presented in section 2.3.

### 3.2.1 Updating the kernel matrix

The inverse of the kernel matrix $K \in \mathbb{R}^{p \times p}$ can be interpreted as the weight covariance matrix, $P_w^f = (K)^{-1}$. We notice that due to the RKHS property expressed in equation (4), $K$ is invertible. To simplify the notation, let us denote $B^f = K^{-1/2}$ as the square root matrix of $P_w^f$, so that $P_w^f = B^f (B^f)^T$. Similarly, we denote the analyzed covariance square root matrix as $B^a = (P_w^a)^{1/2}$. In practice, the square root of the covariance matrix is computed using an eigenvalue decomposition. Analogous to the ESRF, we can then define the square root matrix $S_{\text{RKHS}}$,

$$B^a = B^f S_{\text{RKHS}}, \qquad (37)$$

the EnF gain,

$$G_{\text{RKHS}} = P_w^f (VK)^T (\tilde{R} + (VK) P_w^f (VK)^T)^{-1}, \qquad (38)$$

and the covariance matrix, which is updated according to:

$$P_w^a = (I - G_{\text{RKHS}}(VK)) P_w^f. \qquad (39)$$

Substituting the gain into this equation yields:

$$B^a (B^a)^T = B^f \left( I - (VKB^f)^T (\tilde{R} + (VK)B^f (B^f)^T (VK)^T)^{-1} (VKB^f) \right) (B^f)^T. \qquad (40)$$

Upon applying the Sherman–Morrison–Woodbury formula, the above equation can be rewritten as:

$$B^a (B^a)^T = B^f \left( I + (VKB^f)^T \tilde{R}^{-1} (VKB^f) \right)^{-1} (B^f)^T. \qquad (41)$$

As $B^a (B^a)^T = B^f S_{\text{RKHS}} (S_{\text{RKHS}})^T (B^f)^T$, the RKHS EnF square root matrix is finally identified as:

$$S_{\text{RKHS}} = (I + (VKB^f)^T \tilde{R}^{-1} (VKB^f))^{-1/2}. \qquad (42)$$

The matrix $S_{\text{RKHS}}$ is then used to update the covariance square root matrix, from which the analyzed covariance matrix can be derived. Both the gain $G_{\text{RKHS}}$ and transition matrix $S_{\text{RKHS}}$ are defined for assimilating an observation time series, with these matrices being weighted by the inverse of the kernel matrix.

### 3.2.2 Updating the weights

It is important to note that unlike the ESRF, the RKHS-based approach does not directly provide a formula for updating the ensemble weight anomalies. The matrix $S_{\text{RKHS}}$ does not update the anomalies in the same direct manner. However, since the anomaly can be treated as a centered random vector with a non-singular covariance matrix, it is possible to apply a whitening transformation to this vector. This is achieved by multiplying it by the square root of its covariance matrix inverse, $K^{1/2}$, resulting in a decorrelated vector (Koivunen & Kostinski, 1999; Friedman, 1987). To transform this vector into one that aligns with the analyzed covariance matrix $P_w^a$, a subsequent coloring transformation is applied. This transformation, being the inverse of the whitening process, involves multiplying the vector by the square root of the analyzed covariance matrix $(P_w^a)^{1/2}$. Consequently, the updated weight anomaly can be expressed as follows:
$$ (\beta^a)' = (\beta^f)' (\mathbf{P}_w^f)^{-1/2} (\mathbf{P}_w^a)^{1/2} \\ = (\beta^f)' \mathbf{K}^{1/2} (\mathbf{B}^f \mathbf{S}_{\text{RKHS}}) . \qquad (43) $$

Given that $\mathbf{B}^f = \mathbf{K}^{-1/2}$, it follows that the previous equation can be simplified to

$$ (\beta^a)' = (\beta^f)' \mathbf{S}_{\text{RKHS}} . \qquad (44) $$

The ensemble mean can be then updated in a manner consistent with classical ensemble Kalman filters:

$$ \overline{\beta}^a = \overline{\beta}^f - \mathbf{G}_{\text{RKHS}} (\mathbf{V}\mathbf{K}\overline{\beta}^f - \mathbf{Y}), \qquad (45) $$

and the ensemble weights are updated as:

$$ \beta^a = \overline{\beta}^a + (\beta^f)' \mathbf{S}_{\text{RKHS}}, \qquad (46) $$

similarly as in (31).

A similar approach can be considered in the ESRF context to express the ensemble update. Given that $\mathbf{A}^a = \mathbf{A}^f \mathbf{S}$, the transition matrix reads $\mathbf{S} = (\mathbf{A}^f)^{-1} \mathbf{A}^a$. Within the ESRF framework, where the anomaly is assumed to be the square root of the covariance matrix, the update of the anomaly can also be interpreted as a whitening transformation followed by a coloring transformation. In the ESRF setting, the change of variable $(\mathbf{P}_w^f)^{-1/2} \beta$ allows to work with variables decorrelated from each other – through the diagonalisation of the covariance matrix. In the RKHS-based case, $(\mathbf{P}_w^f)^{-1/2} = \mathbf{K}^{1/2}$, and the equivalent change of variable procedure constitutes an embedding in the RKHS.

### 3.3 Link between ESRF and EnF on a RKHS family

The ESRF can be viewed as a specific instance of RKHS-EnF, given an appropriately selected initial kernel. To establish the connection between these two techniques, we derive in this section the ESRF equations from the RKHS EnF equations by selecting a suitable initial kernel.

Let us consider a RKHS EnF which assimilates observations directly at the time they are available, as a consequence $\mathbf{V} = \mathbf{V}_t$ and $\tilde{\mathbf{R}} = \mathbf{R}_t$ for a fixed time instant $t$. As $\mathbf{k} = \mathbf{K}\beta$, by the definition of $\mathbf{V}$, we have

$$ \mathbf{V}_t \mathbf{K} \beta = \mathbb{H} \mathbf{X}_t, \qquad (47) $$

and it follows that $\mathbf{V}_t \mathbf{K} \beta^f = \mathbb{H} \mathbf{X}_t^f$. However, as seen previously $\beta^f$ is equal to the identity matrix. This represents the forecast members' weights expressed in their own basis, therefore $\mathbf{V}_t \mathbf{K} = \mathbb{H} \mathbf{X}_t^f$. For simplicity, the time index $t$ is omitted in the subsequent discussion since all observations are assumed to be taken simultaneously. Assuming the kernel is chosen such that $(\mathbf{K})^{-1} = \frac{1}{p-1} (\mathbf{I} - \mathbf{1}/p) \mathbf{I} (\mathbf{I} - \mathbf{1}/p)$, the square root of $\mathbf{K}$ is given by $\mathbf{B}^f = \frac{1}{\sqrt{p-1}} (\mathbf{I} - \mathbf{1}/p)$, and the RKHS-EnF gain reads

$$ \mathbf{G}_{\text{RKHS}} = \frac{1}{p-1} (\mathbf{I} - \mathbf{1}/p) \mathbf{I} (\mathbf{I} - \mathbf{1}/p) (\mathbb{H}\mathbf{X}^f)^T \left( \tilde{\mathbf{R}} + \frac{1}{p-1} (\mathbb{H}\mathbf{X}^f) (\mathbf{I} - \mathbf{1}/p) \mathbf{I} (\mathbf{I} - \mathbf{1}/p) (\mathbb{H}\mathbf{X}^f)^T \right)^{-1} . \qquad (48) $$

It is worth noting that the kernel $\mathbf{K}$ defined here has a Mercer spectrum $\mu_i = 1$ for $0 \le i \le p$, and can be interpreted as the zero-correlation limit of a Gaussian, tending towards a Dirac measure, which is not a RKHS reproducing kernel (as it is unbounded). As $\mathbf{X}^f (\mathbf{I} - \mathbf{1}/p) = \mathbf{A}^f$ and $(\mathbf{I} - \mathbf{1}/p)^2 = (\mathbf{I} - \mathbf{1}/p)$, we have

$$ \mathbf{G}_{\text{RKHS}} = \frac{1}{p-1} (\mathbb{H}\mathbf{A}^f)^T \left( \tilde{\mathbf{R}} + \frac{1}{p-1} (\mathbb{H}\mathbf{A}^f) (\mathbb{H}\mathbf{A}^f)^T \right)^{-1} = \mathbf{G}_w . \qquad (49) $$
Considering the matrix $S_{RKHS}$ under the same assumption, we get:

$$S_{RKHS} = \left( I + \left( \mathbb{H}X^f \frac{1}{\sqrt{p-1}} (I - 1/p) \right)^T \tilde{R}^{-1} \left( \mathbb{H}X^f \frac{1}{\sqrt{p-1}} (I - 1/p) \right) \right)^{-1/2}, \quad (50)$$

which simplifies to:

$$S_{RKHS} = \left( I + \frac{1}{p-1} (\mathbb{H}A^f)^T \tilde{R}^{-1} (\mathbb{H}A^f) \right)^{-1/2} = S. \quad (51)$$

Consequently, performing the assimilation step each time an observation is available, using a well-chosen kernel, provides an RKHS EnF gain and a square root matrix that are equivalent to those of a classical square root filter. Furthermore, since under these conditions $VK = \mathbb{H}X$, the two formulas used to update the ensemble are equivalent and the two filters are equivalent.

## 3.4 Ensemble update strategy

The embedding of the dynamics into the RKHS family ensures that the ensemble weights remain invariant over time. This invariance provides the reconstruction of the analyzed ensemble members at any time of the assimilation window, as a linear combination of the forecast ensemble. The weights are calculated with respect to the initial kernel, and consequently, in the method proposed herein, it is natural to construct the analyzed members at the beginning of the assimilation window. Due to the non-linear dynamics of the system and the finite number of ensemble members, the analyzed members reconstructed at the beginning of the assimilation window and subsequently propagated forward in time may differ at the end of the assimilation window from the analyzed members. It is noteworthy that, according to the theoretical framework, if an infinite number of ensemble members were available, for an invertible dynamics the reconstruction at the beginning of the assimilation window, when transported by the dynamical system, would be identical to the reconstruction at the end of the assimilation window. However, obviously due to the finite size of the ensemble members and due to noninvertible numerical dissipation process, this is only approximately true. The ability of a direct reconstruction of the initial condition at the beginning of the assimilation window is used to forecast updated trajectories of the ensemble members. This approach referred to as the tiled RKHS ensemble filter is described in Algorithm 1.

**Remark 4** The tiled RKHS EnF is consistent with the theoretical framework, as all estimations are re-expressed at the initial time $t_0$ through the isometry between the evolving RKHSs. This procedure has also the advantage of introducing a spin-up phase during the assimilation window, where spurious numerical artefacts induced by the finite size linear-combination reconstructions are filtered-out by the dynamical system, and are ultimately replaced by physically relevant features.

**Remark 5** The resulting tiled RKHS ensemble filter shares similarities with Bayesian smoothers and 4D-Var methods, with the key distinction that it avoids multiple forward-backward integrations of the dynamical system, as we infer directly the initial condition of the system (at the beginning of the assimilation window). In the proposed method, the model is integrated forward only once to propagate the analyzed ensemble to the next assimilation window. An intermediate method, which can be referred to as sequential RKHS EnF, can also be devised by reconstructing the ensemble at the end of the assimilation window, akin to classical ensemble filters. Although this approach does not fully leverage the theoretical framework's potential, it maintains the same computational cost as a classical sequential filter.
**Algorithm 1 Tiled RKHS EnF**

**Require:** Ensemble $X_0$, initial time $t_0$, end time $t_{end}$, time step $dt$, assimilation interval

$T_a$.

$X_t^f \leftarrow X_0$

$t \leftarrow t_0$

$t_a \leftarrow T_a$ the assimilation time horizon.

Save $X_t$ for future reconstruction.

Construct the kernel matrix $K$.

Initialize $VK$ and $Y$ as empty.

while $t < t_{end}$ do

$t \leftarrow t + dt$

Propagates the ensemble members in time to update $X_t^f$.

Append to $VK$ the observations of the ensemble members.

Append to $Y$ the reference observations.

if It is time to assimilate i.e. $t$ is a multiple of $t_a$ then

$\beta^a \leftarrow \text{assimilate}(VK, Y, K)$

▷ Equations (45) and (46).

Empty $VK$ and $Y$.

Apply the weights $\beta^a$ to $X_{t-t_a}^f$ saved at $t - t_a$, to obtain $X_{t-t_a}^a$.

Propagates $X_{t-t_a}^a$ on the time interval $[t - t_a, t]$, to obtain $X_t^a$.

$X_t^f \leftarrow X_t^a$

Save $X_t^f$ for future reconstruction.

Update the kernel matrix $K$.

end if

end while

# 4 Experimental framework

In the following we first describe the experimental setup we built to assess the data assimilation we propose.

## 4.1 Dynamical model

The data assimilation methods previously presented are compared using a North Atlantic Ocean model, covering the region between 9°N and 48°N in latitude and -98°E and -3.84°E in longitude with realistic coasts. This area represents a domain of 4329 km by 8658 km. The ocean bottom is assumed to be flat, and all boundaries are considered closed, even those at the equator and in the north, which are normally open boundaries. Due to the flat bottom, the Gulf Stream develops close to the North American coast. To achieve a well-positioned Gulf Stream, areas with a bathymetry shallower than -100 meters are treated as land. A realistic wind forcing is considered. A multi-layer quasi-geostrophic model implemented as the projection of a rotating shallow-water (Thiry et al., 2024) is considered. The rotating shallow water is simulated with the model proposed by Roullet and Gaillard (2022), while the QG projection uses the elliptic solver of Thiry et al. (2023). Three layers of heights $H = (400, 1100, 2600)$ are considered. The evolution of three prognostic variables are simulated: $u$, $v$ the two components of the horizontal velocity vector, and $h$ the height anomaly in each layers. These variables are propagated over time according to the following equations:

$$ \left\{ \begin{aligned} \partial_t u &= (\omega + f) v - \partial_x (p + k) \\ \partial_t v &= - (\omega + f) u - \partial_y (p + k) \\ \partial_t h &= -H (\partial_x u + \partial_y v) - \partial_x (uh) - \partial_y (vh) \end{aligned} \right. , \quad (52) $$

with $\omega = \partial_x v - \partial_y u$ the relative vorticity, $p$ the pressure vector, $k = (u^2 + v^2)/2$ the kinetic energy and $f_i = f_0 + \beta y$ for $0 \le i \le 3$ the Coriolis parameter under a beta plane assumption, with $f_0 = 9.375 \times 10^{-5} \text{ s}^{-1}$ the mean Coriolis parameter, and $\beta = 1.754 \times 10^{-11} \text{ m}^{-1}\text{s}^{-1}$
the Coriolis parameter gradient. The pressure values in layer $1 \le i \le 3$ are defined such that $p_i = \rho_1 \sum_{j=1}^{i} (g'_j \sum_{k=j}^{3} h_k)$ with $g_1 = g = 9.81 \text{ m} \cdot \text{s}^{-2}$ and $g'_i = g(\rho_i - \rho_{i-1})/\rho_1$ for $i = 2, 3$. $\rho$ is the density in the different layers.

The QG equations are then obtained by projecting this shallow water system using a projector involving an elliptic operator. The full details of this procedure as well as the numerical codes can be found in Thiry et al. (2024). The simulation is performed on a 512 × 1024 × 3 grid ($n_x = 512$, $n_y = 1024$, $n_l = 3$), resulting in a mesh size of approximately 8.46 × 8.46 kilometers square. The time step is set to 2000 seconds.

## 4.2 Ensemble generation

The initial ensemble $X_0$ is generated using a noise based on the spectral decomposition of random local fluctuations (Bauer et al., 2020; Brecht et al., 2021). This technique is applied to perturb the velocity field. Given a velocity field $u$ defined on a $n = n_l \times n_x \times n_y$ grid, a local spatial window $W_{i,j}$ of size $n_w \times n_w$ ($n_w$ being odd) in the horizontal direction is considered around each grid point $x_{i,j}$, with $1 < i < n_x$ and $1 < j < n_y$. Within each window, $n_o$ samples are randomly picked, where $n_o > 1$. The sampled values within the window $W_{i,j}$ are denoted as $w_k^{(i,j)}$, where $1 \le k \le n_o$, and these values represent $n_l$-dimensional vectors associated with the number of vertical layers. These values are centered by subtracting their mean as follows:

$$ (\mathbf{w}_k^{(i,j)})' = \mathbf{w}_k^{(i,j)} - \bar{\mathbf{w}}^{(i,j)} \quad \text{with} \quad \bar{\mathbf{w}}^{(i,j)} = \frac{1}{n_o} \sum_{\ell=1}^{n_o} \mathbf{w}_{\ell}^{(i,j)}. \qquad (53) $$

Hence a matrix $Z \in \mathbb{R}^{n \times n_o}(\mathbb{R})$ is constructed by organizing the centered values $(\mathbf{w}_k^{(i,j)})'$ for each grid point in the $n_l \times n_x \times n_y$ grid. A singular value decomposition (SVD) is performed $Z = U \Sigma V^T$, allowing to define perturbation modes $U_i$ from the columns of $U$ and the diagonal matrix of singular values $\Sigma_{jj} = \lambda_j$. The perturbation, $\{X_0^{(i)}, i = 1 \dots, p\}$ are generated from the $p$ eigenpairs of the empirical covariance matrix:

$$ X^{(i)} = \sum_{j=1}^{n_o} \lambda_j U_j \eta_j^{(i)}, \qquad (54) $$

where $\eta \in \mathbb{R}^{n_o}$ is a multivariate random vector composed of $n_o$ standard Gaussian variables and $\eta^{(i)}$ denotes the $i$th realization. The initial ensemble members are obtained by adding these random fields to the initial velocity component. In the experiments, 76 members have been generated using this method, with $n_w = 5$ and $n_o = 21$. One member is retained as the reference, while the remaining 75 members constitute the training set. The reference member is not included in the training ensemble.

## 4.3 Observations simulation

An Observing System Simulation Experiment (OSSE) was conducted using sea surface height (SSH) as the observed variable. The SSH corresponds to the sum of the layer height anomalies $h$ along the vertical. Synthetic observations were generated from the reference simulation by applying spatial masks derived from realistic data provided by the Copernicus along-track sea surface height database (DOI https://doi.org/10.48670/moi-00147). The initial observation date was set to January 1, 2022. Figure 2 illustrates the observation masks for two periods: one over 33 minutes and 20 seconds, corresponding to a single simulation time step, and another over 41 hours and 40 minutes. It can be shown that over 41h40min, the altimeters swaths constitute a reasonable coarse grid observational coverage of the basin. Usual OI methods, provide a dense field of observation, by interpolating these observation spatially assuming they are all acquired at the same time. In the technique we propose, the observations are treated as a time-series and are not interpolated.
<figure><img src="image_4.png" /><figcaption>Figure 2: Observation masks for the periods of 33 minutes and 20 seconds (left) and 41 hours and 40 minutes (right), both from January 1, 2022.</figcaption></figure>

Four distinct experiments were carried out. In the first two, observations were perturbed with additive independent white Gaussian noise. The noise standard deviations were set to 1% and 10% of the root mean square (RMS) of the reference member's SSH at the initial time, corresponding to standard deviations of $\sigma_o = 0.18$ cm and $\sigma_o = 1.8$ cm, respectively. Considering spatially uncorrelated white Gaussian noise is a strong idealization of real-world data noise. In operational setups, the noise is correlated, as raw sea surface height measurements are corrected to remove barotropic tidal effects and the predictable internal tide waves, compensate for satellite jitters, and account for atmospheric turbulence.

Due to this, two additional experiments with spatially correlated noise were conducted. The noise was generated by smoothing a field of white Gaussian noise using a Gaussian filter with a standard deviation of 250 km and a support length of 300 km. This approach produced noise patterns with spatial correlations and, notably, the resulting noise was non-Gaussian due to the relatively limited support area of the filter. The initial variance of the white Gaussian noise was calibrated so that the mean covariance square root of the resulting spatially correlated noise was approximately of 1% and 10% of the SSH RMS (0.18 cm and 1.8 cm respectively). Although this noise is more realistic than white Gaussian noise, it is still obviously idealized. In all experiments, the observation error covariance matrix $\tilde{R}$ used in the filters was assumed to be diagonal, with all observations having the same error variance $\sigma_o^2$. As a result, the filter was not informed of the correlations between observations, posing a challenge for the data assimilation method to remain robust despite this misinformation.

## 4.4 Filters Setting

### 4.4.1 Assimilation frequency

In the context of the RKHS EnF, the assimilation frequency can be adjusted to optimize the reconstruction quality. In the one hand, if the assimilation window is too short, the observation time-series might be too short to exploit all the potential benefit of the RKHS EnF. In the limit case where observations are sequentially assimilated at each time step and the kernel is defined as described in section 3.3, the RKHS EnF and ESRF filters are equivalent. However, an excessively long assimilation window may violate the invertibility assumption of the dynamical system, leading to poorer outcomes for the ensemble linear combination of trajectories over that period. Conversely, for the tiled version of the filter, a longer window can enhance the development of small-scale feature variability and assist in eliminating non-physical structures within the dynamical system.

In the following, for the untiled sequential version of the RKHS EnF, the assimilation window length is set to 8 hours and 20 minutes when the observations are perturbed by
664

665

666

667

668

669

670

671

white Gaussian noise, and to 2 hours, 46 minutes, and 40 seconds when perturbed by spatially correlated noise. These correspond to analysis updates every 15 and 5 time steps of the dynamics, respectively. The presence of spatially correlated noise degrades the filter's performance, necessitating more frequent corrections to the ensemble. For the tiled version, the assimilation window is set to 41 hours and 40 minutes, which corresponds to 75 time steps of the dynamics, regardless of the type of noise affecting the observations. This suggests a robustness in the tiled EnF's configuration, making it less sensitive to the noise characteristics. The ESRF is updated every time step (i.e. every 33 minutes and 20 seconds).

672

673

674

675

676

677

**Remark 6** It is important to note that a tiled version of the ESRF is theoretically not directly feasible. The possibility of using this strategy with the RKHS EnF arises from the stationary weights associated with the linear combination of trajectories and the time series of observations. Obtaining the assimilated state at the beginning of the assimilation windows requires either a forward-backward Bayesian smoother or the application of 4D VAR techniques.

678

### 4.4.2 Kernel selection

679

680

In the experiments, a centered Gaussian kernel, defined by the following expression, was employed:

681

$$ \forall \mathbf{X}, \mathbf{Y} \in \Omega_0, k_g(\mathbf{X}, \mathbf{Y}) = \exp\left(-\frac{\|\mathbf{u}_{\mathbf{X}} - \mathbf{u}_{\mathbf{Y}}\|_{\Omega_x}^2 + \|\mathbf{v}_{\mathbf{X}} - \mathbf{v}_{\mathbf{Y}}\|_{\Omega_x}^2}{l_G^2}\right), \quad (55) $$

682

683

684

685

686

687

where $(\mathbf{u}_{\mathbf{X}}, \mathbf{v}_{\mathbf{X}})$ and $(\mathbf{u}_{\mathbf{Y}}, \mathbf{v}_{\mathbf{Y}})$ denote the velocity vectors associated with states $\mathbf{X}$ and $\mathbf{Y}$, respectively, and $\|\cdot\|_{\Omega_x}$ denoting the $L^2$ norm over the domain $\Omega_x$. The kernel matrix associated with this kernel is denoted $K_G$. The length scale parameter $l_G$ is chosen to maintain a constant ratio between the smallest and largest Mercer eigenvalues of $K_G$ which is set to $10^{-2}$. The Gaussian kernel matrix is then normalized by multiplying it by $p-1$ and then centered by pre- and post-multiplying it by the projection operator $(\mathbf{I} - 1/p)$:

688

$$ K = \alpha(p-1)(\mathbf{I} - 1/p) K_G (\mathbf{I} - 1/p), \quad (56) $$

689

690

691

692

693

694

695

696

where $\alpha$ is a constant. This centering procedure is simply equivalent to a change of variable toward the anomaly of the coefficients $\beta \mapsto \beta'$. Consequently, this is equivalent to work with a full rank kernel matrix (i.e. the kernel and the associated RKHS is well defined), which operates on the weight anomalies. Specifically, $\alpha$ is set to 0.75 for uncorrelated Gaussian noise, 2 for correlated (non-Gaussian) low noise, and 5 for correlated (non-Gaussian) high noise. As a rule of thumb, for correlated noise, $\alpha$ is set to a value greater than 1. This adjustment compensates for the correlations not accounted for in the $\mathbf{R}$ matrix. The higher the correlation between observations, the more $\alpha$ should be increased.

697

698

699

Otherwise, since the rank of $(\mathbf{I} - 1/p)$ is $p-1$, the kernel matrix $K$ is not invertible. Thus, its Moore-Penrose pseudo-inverse must be used. Given that $(\mathbf{I} - 1/p)$ is a projection matrix, it equals its own pseudo-inverse, thus the Moore-Penrose inverse of $K$ is given by:

700

$$ K^{\dagger} = \frac{1}{\alpha \cdot (p-1)} (\mathbf{I} - 1/p) K_G^{-1} (\mathbf{I} - 1/p). \quad (57) $$

701

702

703

704

705

706

The square root of the kernel matrix pseudo-inverse is typically computed through SVD decomposition. The kernel pseudo-inverse and its square root are used to compute the gain (38) and the covariance square root (42). The centering of the kernel ensures consistency with classical square root filtering methods. The kernel used in section 3.3 to demonstrate the connection between the ESRF and the RKHS EnF is a specific case of the kernel presented here, with $l_G$ tending towards 0.
## 4.5 Diagnostics

To assess the reconstruction quality of each methods, we consider the reconstructed relative vorticity fields. Relative vorticity is a less smooth variable compared to velocity or water height, making it a more effective discriminator between the different reconstruction methods. Three metrics are considered to evaluate the effectiveness of the methods. In this subsection the vorticity of the ith ensemble member to be evaluated is denoted as $\omega^{(i)}$ and the reference is denoted as $\omega_{ref}$.

The root mean square error (RMSE) is a widely used metric for quantifying the discrepancy between the ensemble members and the reference state. It is defined as:

$$e_{\text{RMSE}} = \sqrt{\frac{1}{p} \sum_{i=1}^{p} \|\omega^{(i)} - \omega_{\text{ref}}\|_{\Omega_x}^2} \quad (58)$$

It penalizes the deviation of each ensemble member from the reference, providing a measure of the overall accuracy of the ensemble. To improve interpretability, the RMSE can be expressed as a percentage by normalizing it with respect to the norm of the reference state, as follows:

$$e_{\text{NRMSE}} = \frac{e_{\text{RMSE}}}{\|\omega_{\text{ref}}\|_{\Omega_x}} \quad (59)$$

This normalized RMSE offers a more intuitive understanding of the ensemble's mean performance. It is important to note that this measure reflects the mean accuracy of the individual ensemble members, rather than the accuracy of the ensemble mean, which is captured by the bias. The bias is defined as the distance between the reference and the ensemble mean estimation:

$$\omega_{\text{bias}} = \bar{\omega} - \omega_{\text{ref}} \quad (60)$$

By computing the norm of the bias vorticity field over the domain $\Omega_x$, a scalar value can be obtained to represent the bias error across the entire spatial domain:

$$e_{\text{bias}} = \|\omega_{\text{bias}}\|_{\Omega_x} \quad (61)$$

This indicator represents the accuracy of the ensemble, the lower it is the more accurate is the estimation. The standard deviation is expressed as

$$\sigma_e = \sqrt{\frac{1}{p} \sum_{i=1}^{p} \|\omega^{(i)} - \bar{\omega}\|_{\Omega_x}^2} \quad (62)$$

The standard deviation represents the dispersion of the ensemble. A small standard deviation indicates that the ensemble members are close to each others. The three metrics presented here are interconnected, considering the following inequality for $1 \le i \le p$,

$$\|\omega^{(i)} - \omega_{\text{ref}}\|_{\Omega_x}^2 \le \|\omega^{(i)} - \bar{\omega}\|_{\Omega_x}^2 + \|\bar{\omega} - \omega_{\text{ref}}\|_{\Omega_x}^2, \quad (63)$$

which gives with an empirical mean

$$e_{\text{RMSE}}^2 \le \sigma_e^2 + e_{\text{bias}}^2 \quad (64)$$

A low standard deviation may indicate a good result if the bias is also low, meaning that all ensemble members are close to the reference, resulting in a low RMSE. However, if the ensemble spread is too small, there is a risk that the reference may fall outside the ensemble's range, making it difficult for the method to adjust the ensemble members to track the reference in the future. This can lead to divergence over time. Therefore, a larger standard deviation within the ensemble can contribute to more robust predictions.

Striking a balance between standard deviation and bias, known as the bias-variance trade-off, is crucial. A standard deviation that is too low can lead to divergence, while one that is too high may indicate that individual ensemble members poorly represent the system. A wide spread can also impact the bias, as overly dispersed ensemble members may produce an over-smoothed mean, failing to capture small-scale structures accurately.
# 5 Results

In this section, we compare the performance of the data assimilation methods presented in this study over 100 days. Let us mention that the state space dimension is large, while the number of ensemble members is limited, making this scenario particularly challenging for ensemble methods. It is well known that in such cases, localization techniques are recommended to prevent filter divergence. However, we will not consider them here, as the focus is on the efficiency and robustness provided by the RKHS embedding, as well as its capacity to handle a 4D VAR problem at low cost using the tiled simulation strategy. While extending the proposed methods with localization techniques is possible, it lies beyond the scope of this paper and presents a clear avenue for future research.

As an initial comparison, Figure 3 shows the NRMSE for both the classical ESRF and the RKHS EnF. The NRMSE is averaged over five different reference datasets, each with its own initial ensemble, generated as described in Section 4.2. The shaded areas around the mean NRMSE curves represent the standard deviation across these experiments. A wide shaded area indicates high sensitivity of the results to both the initial reference (from which the observations are extracted) and the initial ensemble, while a narrow band highlights low variability in performance, suggesting greater robustness.

In the experiments with 1% Gaussian observation noise, the classical ESRF shows a constant increase in error and significant variability, indicating a lack of robustness. The RKHS EnF performs similarly, though with slightly lower variability. In contrast, the tiled RKHS EnF demonstrates remarkable long-term stability, maintaining consistent performance along time.

For higher, decorrelated noise (10% of the reference height field RMS), all methods experienced an increase in error. However, the tiled version of the RKHS EnF outperformed the ESRF, achieving lower NRMSE with smaller standard deviation across different experiments, demonstrating greater resilience to high noise levels. The ESRF performs better with uncorrelated noise compared to the RKHS EnF, but similarly to the lower noise level, it exhibits slightly higher robustness in terms of a tighter standard deviation.

For 1% spatially correlated non-Gaussian noise, the tiled RKHS EnF significantly outperforms both ESRF and RKHS EnF, maintaining consistently low error and much smaller standard deviation. This highlights its strong robustness against misrepresentation of the observation covariance, which was assumed diagonal and hence representative of independent, uncorrelated spatial noise. ESRF and RKHS EnF perform similarly in terms of NRMSE average and standard deviation for this noise type.

With higher spatially correlated noise (10%), all methods see a more substantial increase in NRMSE. However, both the RKHS EnF and its tiled version significantly outperform the ESRF, showing greater robustness (with a tighter standard deviation) and maintaining a constant NRMSE at the noise level up to day 50. In contrast, the ESRF shows increasing error from around day 15. This experiment reveals a clear separation in the filters' performance. Theoretically, ESRF is not well-suited to nonlinear, non-Gaussian systems, whereas the RKHS EnF appears better equipped to handle non-Gaussianity and correlated noise.

As noted earlier, the performance of these filters could likely be improved with the use of localization techniques. For the RKHS EnF, localization could be achieved by localizing the observables and with the consequence of introducing localized stationary weights (i.e. stationnary spatial vectors of weights) instead of global weights.

These experiments demonstrate that the tiled version of the RKHS EnF offers clear improvements. This approach consistently produces better results, showing less sensitivity to initial conditions (i.e., the initial ensemble) by showcasing a tighter standard deviation for all noise levels. Additionally, the RKHS EnF is notably less sensitive to errors in the
noise covariance modelling. The tiled RKHS EnF benefits from the theoretical properties of embedding the dynamical system in the dynamically evolving RKHS family, allowing it to formulate a 4D-Var assimilation as a single-iteration filtering process, similar to sequential Bayesian filters.

Figure 4 presents the analyzed ensemble mean of the upper-layer relative vorticity after 100 days, using observations perturbed by the 10% non-Gaussian, spatially correlated noise. This figure shows that all assimilation methods successfully reconstruct a vorticity map close to the reference vorticity field. At naked eye the estimations are almost indistinguishable from the reference. The enstrophy spectra from the upper layer, calculated within the rectangular region bounded by -72.25°E to -42.28°E longitude and 14.64°N to 39.47°N latitude, following the method described by Durran et al. (2017) are shown in Figure 5. This figure reveals that the enstrophy spectra for the three methods are nearly identical to the reference one, indicating accurate reproduction of the reference field energy across all scales. Enstrophy spectra for the other layers lead to the same conclusion (not shown). Importantly, the effective resolution of the reference data has not been compromised. This suggests that employing ensemble methods, without resorting to optimal interpolations or gridded data prior to assimilation, maintains the resolution integrity.

However, a clearer understanding of the differences between the methods is provided by the bias fields, which represent the discrepancies between the reference and analyzed fields. Figures 6, figure 7, and 8 illustrate the relative vorticity bias fields for the three assimilation methods across the upper, middle, and lower layers, respectively. Notably, the middle and lower layers are better reconstructed by all methods compared to the upper layer, likely due to the reduced variability and the fewer small-scale structures in these layers. To improve clarity, the color scale for the bias fields in these two layers has been refined. An analysis of the bias fields reveals similar trends across all three layers. The classical ESRF exhibits significant bias, particularly along the Gulf Stream, indicating poorer performance in capturing the dynamics of this region. In contrast, the RKHS EnF demonstrates improved results, with lower bias compared to the ESRF. The tiled RKHS EnF outperforms both methods, producing substantially lower biases. This suggests that the tiled method assimilates the data more effectively, leading to more accurate reconstructions. The study of the bias fields aligns with the conclusions drawn from the RMSE analysis done in the previous paragraph. This consistency is expected, given the inherent connection between these diagnostic variables.

To further quantify the errors, Figure 9 presents the enstrophy bias spectra for the upper layer. The spectra indicate that most errors occur at wavenumbers greater than $10^{-2}$, corresponding to structures smaller than 600 km. The bias error energy of the RKHS EnF is slightly lower than that of the ESRF, while the tiled RKHS EnF exhibits significantly reduced error energy. Notably, the improvement is substantial across all scales, from large-scale features to the eddy-diffusion cut-off, with the greatest enhancement around 300 km. This highlights the superior performance of the tiled method in capturing both large-scale and small-scale features in the vorticity field.

To finally complement these analysis results, we present in Figure 10 the RMSE obtained for the upper layer. As observed, the ESRF exhibits the highest RMSE, while the RKHS EnF shows a lower RMSE, and the Tiled RKHS EnF further improves upon it. The other layers, which follow the same trend of lower RMSE, are not shown. The spread (standard deviation) of the ensemble for the upper layer is also displayed in Figure 11. We observe that the spread is attenuated for the ESRF, suggesting a collapse of the ensemble. Both the RKHS EnF and Tiled RKHS EnF display some variation in spread. While the spread is weak for the Tiled RKHS EnF, the bias is also minimal, indicating a less problematic situation compared to the ESRF.
<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of the Normalized Root Mean Square Error (NRMSE) for the ESRF (blue), RKHS EnF (red), and Tiled RKHS EnF (green). The first two rows display results for experiments with additive white Gaussian noise, while the last two rows present results for experiments with non-Gaussian, spatially correlated noise. Rows one and three correspond to experiments with 1% standard deviation noise, and rows two and four to experiments with 10% standard deviation noise.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Reference upper layer relative vorticity field (top row) after 100 days, analyzed fields for the ESRF (second row), RKHS EnF (third row), and tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Upper layer enstrophy spectra after 100 days within the region defined by -72.25°E to -42.28°E longitude and 14.64°N to 39.47°N latitude, presented as a function of wavenumber. The reference spectra (black) is compared to the results from the ESRF (blue), RKHS EnF (orange), and tiled RKHS EnF (green), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Reference upper layer relative vorticity field (top row) after 100 days, bias field corresponding to it for the ESRF (second row), RKHS EnF (third row), tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Reference middle layer realtive vorticity field (top row) after 100 days, bias field corresponding to it for the ESRF (second row), RKHS EnF (third row), tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Reference lower layer real-time vorticity field (top row) after 100 days, bias field corresponding to it for the ESRF (second row), RKHS EnF (third row), and tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Upper layer enstrophy bias spectra after 100 days within the region bounded by -72.25°E to -42.28°E longitude and 14.64°N to 39.47°N latitude, plotted as a function of wavenumber. Results are shown for the ESRF (blue), EnF on an RKHS family (orange), and tiled EnF on an RKHS family (green), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Upper layer relative vorticity RMSE field after 100 days for the ESRF (top row), RKHS EnF (second row), tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Upper layer relative vorticity standard deviation field after 100 days for the ESRF (top row), RKHS EnF (second row), tiled RKHS EnF (bottom row), all corresponding to the experiment with 10% standard deviation spatially correlated noise.</figcaption></figure>
# 6 Conclusion and perspectives

This study introduces a novel approach for assimilating observations of a non-linear dynamical system by using a family W of reproducing kernel Hilbert spaces (RKHS). Each RKHS in this family is associated with kernels $k_t$, which are defined by transporting an initial kernel $k_0$ along the system's dynamics up to time $t$. This framework establishes an isometry between the evolving RKHS, fully justifying the superposition principle, allowing new system trajectories to be constructed as linear combinations of known ensemble members trajectories. Importantly, the isometry ensures these linear combinations remain invariant over time.

An ensemble filter is developed within this RKHS family. The theoretical properties of this framework naturally enable the assimilation of time series of observations, effectively accumulating data at different time instants before analyzing the entire trajectory within the assimilation window. This leads to more robust outcomes. By leveraging the framework's properties, the analyzed ensemble can be provided at the start of the assimilation window and propagated forward. This tiled version of the filter has demonstrated superior performance in reconstructing physical structures and improving accuracy in estimating the reference state.

The methods presented in this study have been applied to a quasi-geostrophic model of the North Atlantic Ocean and compared to a classical square root filter (ESRF), which can be viewed as a limit of the sequential version of the RKHS-based ensemble filter.

Improvements to the RKHS EnF could be made by exploring methods to optimize the choice of kernels. One such approach, as proposed by Akian et al. (2022), involves an algorithm that adjusts the kernel eigenvalues to minimize accuracy loss when data are removed from a training set. Investigating alternative kernels beyond the Gaussian could also be beneficial. For example, as noted by Stein (1999), the Gaussian kernel may be too smooth, whereas the Matérn kernel offers greater flexibility in terms of smoothness while maintaining a reasonable number of parameters to configure.

One significant obstacle to improving state reconstructions in ensemble-based filtering methods is the inherent limited number of ensemble members that can be computationally afforded. Localization methods, as discussed by Sakov and Bertino (2010), help address this limitation. Localization is introduced either by localizing the observation in the assimilation process or by localizing the covariance matrix. This second approach seems more difficult to implement with the RKHS EnF, which operates directly in the weight space. The first approach, which result in local filter is more appealing and doable with the introduction of localized observable and considering spatial vector of stationary weights. In the same spirit Bocquet (2016) introduces a localization strategy that reconstructs trajectories through linear combinations, where the localization operator is propagated over time. This class of methods could potentially be naturally adapted within the time-evolving RKHS framework to introduce dynamical localization process.

In addition let us mention that, applying a tiled filter with local analysis could mitigate imbalances from localized approximate reconstructions.

Another perspective is that the framework allows direct access to the system's Lyapunov times (Dufée et al., 2024), which indicate when the predictability horizon has been exceeded. This information is valuable for determining the optimal moments for assimilation, facilitating the renewal of the ensemble and creating an adaptive assimilation step filter.

Additionally, the initial ensemble significantly impacts the method's performance, making it crucial to select a well-chosen starting set. Generative methods can be employed to establish an initial ensemble that best represents the reference simulation's initial conditions, thereby enhancing the filtering method's effectiveness.
In this study, the reference dynamics were modeled using the same model and resolution as the training set elements. However, modeling a reference that differs from the training set ensemble poses challenges. A different resolution can lead to variations in energy levels between the reference and training ensemble members; for instance, higher-resolution reference models generally exhibit more energy than their coarser counterparts. To address this, variability can be introduced into the training ensemble through methods such as inflation or by employing a stochastic model with calibrated noise (Dufée et al., 2022).

Moreover, to model the reference differently from the training set, the quasi-geostrophic projection can be omitted in the reference model run. This approach would enable the reference to capture fast waves, which would not be reconstructed by the training ensemble. Incorporating multi-modal data could further enhance the assimilation process. For example, combining sea surface height data with sea surface temperature and integrating data from the Surface Water Ocean Topography (SWOT) mission alongside traditional altimeter data could improve the assimilation of real data, as discussed by (Le Guillou, Lahaye, et al., 2021).

## Open Research Section

The data and the code files necessary to reproduce the results presented in this paper are publicly available at: https://gitlab.inria.fr/mael.jaouen/RKHS_EnF/

## Acknowledgments

We acknowledge the support of the ERC EU project 856408- STUOD, the French National program LEFE (Les Enveloppes Fluides et l'Environnement), as well as the support of Inria, which enabled the research presented in this paper through Inria exploratory actions. The authors acknowledge the Pôle de Calcul et de Données Marines (PCDM) for providing DATARMOR computational resources. This study has been conducted using E.U. Copernicus Marine Service Information.

## References

Akian, J.-L., Bonnet, L., Owhadi, H., & Savin, E. (2022). Learning “best” kernels from data in gaussian process regression. with application to aerodynamics. *Journal of Computational Physics*, 470, 111595.

Aronszajn, N. (1950). Theory of reproducing kernels. *Transactions of the American mathematical society*, 68(3), 337–404.

Auroux, D., & Blum, J. (2008). A nudging-based data assimilation method: the back and forth nudging (bfn) algorithm. *Nonlinear Processes in Geophysics*, 15(2), 305–319.

Bauer, W., Chandramouli, P., Chapron, B., Li, L., & Mémin, E. (2020). Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: A stochastic approach. *Journal of Physical Oceanography*, 50(4), 983–1003.

Bishop, A. N., & Del Moral, P. (2023). On the mathematical theory of ensemble (linear-gaussian) Kalman–bucy filtering. *Mathematics of Control, Signals, and Systems*, 35(4), 835–903.

Bishop, C. H., Etherton, B. J., & Majumdar, S. J. (2001). Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects. *Monthly Weather Review*, 129(3), 420–436.

Bocquet, M. (2016). Localization and the iterative ensemble Kalman smoother. *Quarterly Journal of the Royal Meteorological Society*, 142(695), 1075-1089.

Bocquet, M., & Sakov, P. (2014, July). An iterative ensemble Kalman smoother. *Quarterly Journal of the Royal Meteorological Society*, 140(682), 1521-1535.

Brajard, J., Carrassi, A., Bocquet, M., & Bertino, L. (2020). Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observa-
tions: A case study with the lorenz 96 model. *Journal of Computational Science*, **44**, 101171.

949

950

Brecht, R., Li, L., Bauer, W., & Mémin, E. (2021). Rotating shallow water flow under location uncertainty with a structure-preserving discretization. *Journal of Advances in Modeling Earth Systems*, **13**(12), e2021MS002492.

951

952

953

Brunton, S. L., Budišić, M., Kaiser, E., & Kutz, J. N. (2022). Modern koopman theory for dynamical systems. *SIAM Review*, **64**(2), 229-340.

954

955

Budišić, M., Mohr, R., & Mezić, I. (2012, 12). Applied koopmanisma). *Chaos: An Inter-disciplinary Journal of Nonlinear Science*, **22**(4), 047510.

956

957

Carrassi, A., Bocquet, M., Bertino, L., & Evensen, G. (2018). Data assimilation in the geosciences: an overview of methods, issues, and perspectives. *WIREs Climate Change*, **9**.

958

959

960

Cheng, S., Quilodrán-Casas, C., Ouala, S., Farchi, A., Liu, C., Tandeo, P., ... Arcucci, R. (2023). Machine learning with data assimilation and uncertainty quantification for dynamical systems: A review. *IEEE/CAA Journal of Automatica Sinica*, **10**(6), 1361-1387.

961

962

963

964

Das, S., & Giannakis, D. (2020). Koopman spectra in reproducing kernel hilbert spaces. *Applied and Computational Harmonic Analysis*, **49**(2), 573-607.

965

966

Desroziers, G., Camino, J.-T., & Berre, L. (2014). 4DEnVar: link with 4D state formulation of variational assimilation and different possible implementations. *Quarterly Journal of the Royal Meteorological Society*, **140**(684), 2097-2110.

967

968

969

Doucet, A., Freitas, N. D., & Gordon, N. (2001). *Sequential Monte Carlo methods in practice*. Springer Verlag.

970

971

Dufée, B., Hug, B., Mémin, E., & Tissot, G. (2024). Ensemble forecasts in reproducing kernel Hilbert space family. *Physica D: Nonlinear Phenomena*, **459**, 134044.

972

973

Dufée, B., Mémin, E., & Crisan, D. (2022). Stochastic parametrization: An alternative to inflation in ensemble Kalman filters. *Quarterly Journal of the Royal Meteorological Society*, **148**(744), 1075-1091.

974

975

976

Durran, D., Weyn, J. A., & Menchaca, M. Q. (2017). Practical considerations for computing dimensional spectra from gridded data. *Monthly Weather Review*, **145**(9), 3901-3910.

977

978

Eisner, T., Farkas, B., Haase, M., & Nagel, R. (2015). *Operator theoretic aspects of ergodic theory*. Springer.

979

980

Evensen, G. (2009). *Data assimilation, the ensemble Kalman filter*. Springer Berlin, Heidelberg.

981

982

Fablet, R., Chapron, B., Drumetz, L., Mémin, E., Pannekoucke, O., & Rousseau, F. (2021). Learning variational data assimilation models and solvers. *Journal of Advances in Modeling Earth Systems*, **13**(10), e2021MS002572 (e2021MS002572 2021MS002572)

983

984

985

Fablet, R., Chapron, B., Le Sommer, J., & Sévellec, F. (2024). Inversion of sea surface currents from satellite-derived sst-ssh synergies with 4DVarNets. *Journal of Advances in Modeling Earth Systems*, **16**(6), e2023MS003609 (e2023MS003609 2023MS003609)

986

987

988

Farchi, A., Laloyaux, P., Bonavita, M., & Bocquet, M. (2021). Using machine learning to correct model error in data assimilation and forecast applications. *Quarterly Journal of the Royal Meteorological Society*, **147**(739), 3067-3084.

989

990

991

Friedman, J. H. (1987). Exploratory projection pursuit. *Journal of the American Statistical Association*, **82**(397), 249-266.

992

993

Ghil, M., Allen, M. R., Dettinger, M. D., Ide, K., Kondrashov, D., Mann, M. E., ... Yiou, P. (2002, 2025/02/04). Advanced spectral methods for climatic time series. *Reviews of Geophysics*, **40**(1), 3-1-3-41.

994

995

996

Gordon, N., Salmond, D., & Smith, A. (1993, April). Novel approach to nonlinear/non-gaussian bayesian state estimation. *IEE Proceedings F (Radar and Signal Processing)*, **140**, 107-113(6).

997

998

999

Gottwald, G. A., & Reich, S. (2021). Supervised learning from noisy observations: Combining machine-learning techniques with data assimilation. *Physica D: Nonlinear Phenomena*, **423**, 132911.

1000

1001

1002

Hogg, J., Fonoberova, M., & Mezić, I. (2020). Exponentially decaying modes and long-term

1003
prediction of sea ice concentration using koopman mode decomposition. *Scientific Reports*, 10(1), 16313.

Hoke, J. E., & Anthes, R. A. (1976). The initialization of numerical models by a dynamic-initialization technique. *Monthly Weather Review*, 104(12), 1551–1556.

Hunt, B. R., Kalnay, E., Kostelich, E., Ott, E., Patil, D., Sauer, T., ... Zimin, A. (2004). Four-dimensional ensemble Kalman filtering. *Tellus A: Dynamic Meteorology and Oceanography*, 56(4), 273–277.

Hunt, B. R., Kostelich, E. J., & Szunyogh, I. (2007). Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. *Physica D: Nonlinear Phenomena*, 230(1), 112-126. (Data Assimilation)

Klus, S., Schuster, I., & Muandet, K. (2020). Eigendecompositions of transfer operators in reproducing kernel hilbert spaces. *Journal of Nonlinear Science*, 30(1), 283–315.

Koivunen, A. C., & Kostinski, A. B. (1999). The feasibility of data whitening to improve performance of weather radar. *Journal of Applied Meteorology*, 38(6), 741–749.

Kondrashov, D., Ryzhov, E., & Berloff, P. (2020). Data-adaptive harmonic analysis of oceanic waves and turbulent flows. *Chaos*, 30(6), 061105. Retrieved from https://api.semanticscholar.org/CorpusID:220307281

Koopman, B. O. (1931). Hamiltonian systems and transformation in Hilbert space. *Proceedings of the National Academy of Sciences*, 17(5), 315–318.

Kutz, J. N., Fu, X., & Brunton, S. L. (2016). Multiresolution dynamic mode decomposition. *SIAM Journal on Applied Dynamical Systems*, 15(2), 713-735.

Le Gland, F., Monbet, V., & Tran, V. (2011). Large sample asymptotics for the ensemble Kalman filter. In D. Crisan & B. Rozovskii (Eds.), *Handbook on nonlinear filtering*. Oxford University Press.

Le Dimet, F.-X., & Talagrand, O. (1986). Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. *Tellus A*, 38A(2), 97-110.

Leeuwenburgh, O., Evensen, G., & Bertino, L. (2005). The impact of ensemble filter definition on the assimilation of temperature profiles in the tropical pacific. *Quarterly Journal of the Royal Meteorological Society*, 131(613), 3291-3300.

Le Guillou, F., Lahaye, N., Ubelmann, C., Metref, S., Cosme, E., Ponte, A., ... Vidard, A. (2021). Joint estimation of balanced motions and internal tides from future wide-swath altimetry. *Journal of Advances in Modeling Earth Systems*, 13(12), e2021MS002613. (e2021MS002613 2021MS002613)

Le Guillou, F., Metref, S., Cosme, E., Ubelmann, C., Ballarotta, M., Le Sommer, J., & Verron, J. (2021). Mapping altimetry in the forthcoming swot era by back-and-forth nudging a one-layer quasigeostrophic model. *Journal of Atmospheric and Oceanic Technology*, 38(4), 697–710.

Le Traon, P. Y., Faugère, Y., Hernandez, F., Dorandeu, J., Mertz, F., & Ablain, M. (2003). Can we merge GEOSAT follow-on with TOPEX/Poseidon and ERS-2 for an improved description of the ocean circulation? *Journal of Atmospheric and Oceanic Technology*, 20(6).

Liu, C., Xiao, Q., & Wang, B. (2008). An ensemble-based four-dimensional variational data assimilation scheme. part I: Technical formulation and preliminary test. *Monthly Weather Review*, 136(9), 3363–3373.

Lorenc, A. C. (2003). The potential of the ensemble Kalman filter for NWP—a comparison with 4D-Var. *Quarterly Journal of the Royal Meteorological Society*, 129(595), 3183-3203.

Luo, X. (2019, July). Ensemble-based kernel learning for a class of data assimilation problems with imperfect forward simulators. *PLOS ONE*, 14(7), e1–40.

Manucharyan, G. E., Siegelman, L., & Klein, P. (2021). A deep learning approach to spatiotemporal sea surface height interpolation and estimation of deep currents in geostrophic ocean turbulence. *Journal of Advances in Modeling Earth Systems*, 13(1), e2019MS001965. (e2019MS001965 2019MS001965)

Mauran, S., Mouysset, S., Simon, E., & Bertino, L. (2023). A kernel extension of the ensemble transform Kalman filter. In J. Mikyška, C. de Mulatier, M. Paszynski,
V. V. Krzhizhanovskaya, J. J. Dongarra, & P. M. Sloot (Eds.), *Computational science – iccs 2023* (pp. 438–452). Cham: Springer Nature Switzerland.

1059

1060

Mezić, I. (2005). Spectral properties of dynamical systems, model reduction and decompositions. *Nonlinear Dynamics, 41*(1-3), 309–325.

1061

1062

Mezić, I. (2021). Koopman operator, geometry, and learning of dynamical systems. *Not. Am. Math. Soc.*, **68**(7), 1087–1105.

1063

1064

Navarra, A., Tribbia, J., & Klus, S. (2021). Estimation of koopman transfer operators for the equatorial pacific sst. *Journal of the Atmospheric Sciences*, **78**(4), 1227–1244.

1065

1066

Ouala, S., Fablet, R., Herzet, C., Chapron, B., Pascual, A., Collard, F., & Gaultier, L. (2018). Neural network based Kalman filters for the spatio-temporal interpolation of satellite-derived sea surface temperature. *Remote Sensing*, **10**(12).

1067

1068

1069

Reich, S., & Cotter, C. (2015). *Probabilistic forecasting and bayesian data assimilation*. Cambridge University Press.

1070

1071

Roullet, G., & Gaillard, T. (2022). A fast monotone discretization of the rotating shallow water equations. *Journal of Advances in Modeling Earth Systems, 14*(2), e2021MS002663. (e2021MS002663 2021MS002663)

1072

1073

1074

Rowley, C. W., Mezić, I., Bagheri, S., Schlatter, P., & Henningson, D. S. (2009). Spectral analysis of nonlinear flows. *Journal of Fluid Mechanics*, **641**, 115–127.

1075

1076

Sakov, P., & Bertino, L. (2010). Relation between two common localisation methods for the enkf comput. *Computational Geosciences*, **15**, 225-237.

1077

1078

Schmid, P. J. (2010). Dynamic mode decomposition of numerical and experimental data. *Journal of Fluid Mechanics*, **656**, 5–28.

1079

1080

Stein, M. L. (1999). *Interpolation of spatial data*. Springer New York.

1081

Taburet, G., Sanchez-Roman, A., Ballarotta, M., Pujol, M.-I., Legeais, J.-F., Fournier, F., ... Dibarboure, G. (2019). Duacs dt2018: 25 years of reprocessed sea level altimetry products. *Ocean Science*, **15**(5), 1207–1224.

1082

1083

1084

Thiry, L., Li, L., Mémin, E., & Roullet, G. (2024). A unified formulation of quasi-geostrophic and shallow water equations via projection. *Journal of Advances in Modeling Earth Systems, 16*(10), e2024MS004510.

1085

1086

1087

Thiry, L., Li, L., Roullet, G., & Mémin, E. (2023). Mqgeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries. *EGUsphere*, **2023**, 1–25.

1088

1089

Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L., & Kutz, J. N. (2014). On dynamic mode decomposition: Theory and applications. *Journal of Computational Dynamics*, **1**(2), 391–421.

1090

1091

1092

Ubelmann, C., Dibarboure, G., Gaultier, L., Ponte, A., Ardhuin, F., Ballarotta, M., & Faugère, Y. (2021). Reconstructing ocean surface current combining altimetry and future spaceborne doppler data. *Journal of Geophysical Research: Oceans*, **126**(3), e2020JC016560. (e2020JC016560 2020JC016560)

1093

1094

1095

1096

Ubelmann, C., Klein, P., & Fu, L.-L. (2015). Dynamic interpolation of sea surface height and potential applications for future high-resolution altimetry mapping. *Journal of Atmospheric and Oceanic Technology*, **32**(1).

1097

1098

1099

Van Leeuwen, P. J., Künsch, H. R., Nerger, L., Potthast, R., & Reich, S. (2019). Particle filters for high-dimensional geoscience applications: A review. *Quarterly Journal of the Royal Meteorological Society*, **145**(723), 2335-2365.

1100

1101

1102

Wan, L., Bertino, L., & Zhu, J. (2010). Assimilating altimetry data into a hycom model of the pacific: Ensemble optimal interpolation versus ensemble Kalman filter. *Journal of Atmospheric and Oceanic Technology*, **27**(4), 753–765.

1103

1104

1105

Whitaker, J. S., & Hamill, T. M. (2002). Ensemble data assimilation without perturbed observations. *Monthly Weather Review*, **130**(7), 1913–1924.

1106

1107

Williams, M. O., Kevrekidis, I. G., & Rowley, C. W. (2015). A data-driven approximation of the koopman operator: Extending dynamic mode decomposition. *Journal of Nonlinear Science*, **25**(6), 1307–1346. Retrieved from https://doi.org/10.1007/s00332-015-9258-5 doi: 10.1007/s00332-015-9258-5

1108

1109

1110

1111

Williams, M. O., Rowley, C. W., & Kevrekidis, I. G. (2015). A kernel-based method for data-driven koopman spectral analysis. *Journal of Computational Dynamics*, **2**(2),

1112

1113
1114 247-265.

1115 Yang, Y., Robinson, C., Heitz, D., & Mémin, E. (2015). Enhanced ensemble-based 4DVar scheme for data assimilation. *Computers & Fluids*, 115, 201-210.

1116 Zhang, Z., Archibald, R., & Bao, F. (2022). A pde-based adaptive kernel method for solving optimal filtering problems. *Journal of Machine Learning for Modeling and Computing*, 3(3).

1117
1118
1119
1120
1121
1122
Zhen, Y., Chapron, B., Mémin, E., & Peng, L. (2022). Eigenvalues of autocovariance matrix: A practical method to identify the koopman eigenfrequencies. *Phys. Rev. E*, 105, 034205.