Data-Driven Interpolation of Sea Level Anomalies Using Analog Data Assimilation

Abstract

From the recent developments of data-driven methods as a means to better exploit large-scale observation, simulation and reanalysis datasets for solving inverse problems, this study addresses the improvement of the reconstruction of higher-resolution Sea Level Anomaly (SLA) fields using analog strategies. This reconstruction is stated as an analog data assimilation issue, where the analog models rely on patch-based and Empirical Orthogonal Functions (EOF)-based representations to circumvent the curse of dimensionality. We implement an Observation System Simulation Experiment (OSSE) in the South China Sea. The reported results show the relevance of the proposed framework with a significant gain in terms of Root Mean Square Error (RMSE) for scales below 100 km. We further discuss the usefulness of the proposed analog model as a means to exploit high-resolution model simulations for the processing and analysis of current and future satellite-derived altimetric data with regard to conventional interpolation schemes, especially optimal interpolation.

1. Introduction

Over the last two decades, ocean remote sensing data has benefited from numerous remote earth observation missions. These satellites measured and transmitted data about several ocean properties, such as sea surface height, sea surface temperature, ocean color, ocean current, sea ice, etc. This has helped building big databases of valuable information and represents a major opportunity for the interplay of ideas between the ocean remote sensing community and the data science community. Exploring machine learning methods in general and non-parametric methods in particular is now feasible and is increasingly drawing the attention of many researchers [1,2].

More specifically, analog forecasting [3] which is among the earliest statistical methods explored in geoscience benefits from recent advances in data science. In short, analog forecasting is based on the assumption that the future state of a system can be predicted throughout the successors of past (or simulated) similar situations (called analogs). The amount of currently available remote sensing and simulation data offers analog methods a great opportunity to catch up their early promises. Several recent works involving applications of analog forecasting methods in geoscience fields contribute in the revival of these methods, recent applications comprise the prediction of soil moisture anomalies [4], the prediction of sea-ice anomalies [5], rainfall nowcasting [6], numerical weather prediction [7,8,9], etc. One may also cite methodological developments such as dynamically-adapted kernels [10] and novel parameter estimation schemes [11]. Importantly, analog strategies have recently been extended to address data assimilation issues within the so-called analog data assimilation (AnDA) [12,13], where the dynamical model is stated as an analog forecasting model and combined to state-of-the-art stochastic assimilation procedures such as Ensemble Kalman filters.

Producing time-continuous and gridded maps of Sea Surface Height (SSH) is a major challenge in ocean remote sensing with important consequences on several scientific fields from weather and climate forecasting to operational needs for fisheries management and marine operations (e.g., [14]). The reference gridded SSH product commonly used in the literature is distributed by Copernicus Marine Environment Monitoring Service (CMEMS) (formerly distributed by AVISO+). This product relies on the interpolation of irregularly-spaced along-track data using an Optimal Interpolation (OI) method [15,16]. While OI is relevant for the retrieval of horizontal scales of SSH fields up to ≈100 km, the prescribed covariance priors lead to smoothing out finer-scales. Typically, horizontal scales from a few tens of kilometers to ≈100 km may be poorly resolved by OI-derived SSH fields, while they may be partially revealed by along-track altimetric data. This has led to a variety of research studies to improve the reconstruction of the altimetric fields. One may cite both methodological alternatives to OI, for instance locally-adapted convolutional models [17] and variational assimilation schemes using model-driven dynamical priors [18], as well as studies exploring the synergy between different sea surface tracers, especially the synergy between SSH and SST (Sea Surface Temperature) fields and Surface Quasi-Geostrophic dynamics [17,19,20,21,22,23].

In this work, we build upon our recent advances in analog data assimilation and its application to high-dimensional fields. While the works in [12,13] presented the AnDA framework by combining the analog forecasting method and stochastic filtering, these works have only shown applications to geophysical toy models. It was not until the work in [24] that the AnDA methodology was applied to realistic high dimensional fields, namely, Sea Surface Temperature (SST). Dealing with the curse of dimensionality is a critical challenge, in [24] we have shown that the use of patch-based representations (a patch is a term used by the image processing community to refer to smaller image parts of a given global image [25]) combined with EOF-based representations (EOF stands for Empirical Orthogonal Function, a classic dimensionality reduction technique also known as Principal Component Analysis (PCA)) leads to a computationally-efficient interpolation of missing data in SST maps outperforming classical OI-based interpolation schemes. Another development in AnDA applied to high dimensional fields was the introduction of conditional and physically-derived operators [26], where the analog forecasting operators account for the theoretical studies relating to synergies between ocean variables (e.g., SSH and SST) and those highlighting the importance of inter-scale dependencies. In this paper, we make use of these previously developed methodologies and tools and apply the AnDA to Sea Level Anomaly fields. The contribution of this work is two-fold: (i) Confronting AnDA to the reconstruction of an ocean tracer with scarce observations compared to SST (due to the nature of the available altimeters); (ii) designing an Observation System Simulation Experiment (OSSE) based on numerical simulation data to build the archived datasets used for the analog search; (iii) Reconstructing Sea Level Anomaly (SLA) by using SST or large scale SLA as auxiliary variables embedded in the analog regression techniques as shown in Section 4.

Using OFES (Ocean General Circulation Model (OGCM) for the Earth Simulation) numerical simulations [27,28], we design an Observation System Simulation Experiment (OSSE) for a case-study in the South China Sea using real along-track sampling patterns of spaceborne altimeters. Several particular mesoscale variation patterns characterizing this region were studied in the literature, we refer the reader to [29] and references therein. We also note that our method is not region specific and can be applied to any region of interest . Using the resulting groundtruthed dataset, we perform a qualitative and quantitative evaluation of the proposed scheme, including comparisons to state-of-the-art schemes.

The remainder of the paper is organized as follows: Section 2 presents the different datasets used in this paper to design an OSSE, Section 3 gives insights on the classical methods used for mapping SLA from along track data, Section 4 introduces the proposed analog data assimilation model. Experimental settings are detailed in Section 5 and results for the considered OSSE are shown in Section 6. Section 7 further discusses the key aspects of this work.

2. Data: OFES (OGCM for the Earth Simulator)

An Observation System Simulation Experiment (OSSE) based on numerical simulations is considered to assess the relevance of the proposed analog assimilation framework. Our OSSE uses these numerical simulations as a groundtruthed dataset from which simulated along-track data are produced. We describe further the data preparation setup in the following sections.

2.1. Model Simulation Data

The Ocean General Circulation Model (OGCM) for the Earth Simulator (OFES) is considered in this study as the true state of the ocean. The simulation data is described in [27,28]. The coverage of the model is ${75}^{\circ }$S–${75}^{\circ }$N with a horizontal resolution of $1/{10}^{\circ }$. 34 years (1979–2012) of 3-daily simulation of SSH maps are considered, we proceed to a subtraction of a temporal mean to obtain SLA fields. In this study, our region of interest is located in the South China Sea (${105}^{\circ }$E to ${117}^{\circ }$E, $5^{\circ }$N to ${25}^{\circ }$N). This dataset is split into a training dataset corresponding to the first 33 years (4017 SLA maps) and a test dataset corresponding to the last year of the time series (122 SLA maps).

2.2. Along Track Data

We consider a realistic situation with a high rate of along track data. More precisely we use along-track data positions registered in 2014 where four satellites (Jason2, Cryosat2, Saral/AltiKa, HY-2A) were operating. Data is distributed by Copernicus Marine and Environment Monitoring Service (CMEMS).

From the reference 3-daily SLA dataset and real along-track data positions, we generate simulated along-track data from the sampling of a reference SLA field: More precisely, for a given along-track point, we sample the closest position of the $1/{10}^{\circ }$ regular model grid at the closest time step of the 3-daily model time series. As we consider a 3-daily assimilation time step (see Section 2.1 for details), we create a 3-daily pseudo-observation field, to be fed directly to the assimilation model. For a given time t, we combine all along-track positions for times $t-1$, t and $t+1$ to create an along-track pseudo-observation field at time t. We denote by $s3dAT$ the simulated 3-daily time series of along-track pseudo-observation fields. An example of these fields is given in Figure 1.

3. Problem Statement and Related Work

3.1. Data Assimilation and Optimal Interpolation

Data assimilation consists in estimating the true state of a physical variable $\mathbf{x}\left(t\right)$ at a specific time t, by combining (i) equations governing the dynamics of the variable, (ii) available observations $\mathbf{y}(1,\dots ,T)$ measuring the variable and (iii) a background or first guess on its initial state ${\mathbf{x}}^b$. The estimated state is generally called the analyzed state and noted by ${\mathbf{x}}^a$. Data assimilation is a typical example of inverse problems, and similar formulations are known to the statistical signal processing community through optimal control and estimation theory [30]. We adopt here the unified notation of [31] and formulate the problem as a stochastic system in the following:

$$\{\begin{array}{}\hfill \mathrm{(1)}& \mathbf{x}\left(t\right)=\mathcal{M}\left(\mathbf{x}\right(t-1\left)\right)+\eta \left(t\right),\hfill \hfill \mathrm{(2)}& \mathbf{y}\left(t\right)=\mathcal{H}\left(\mathbf{x}\right(t\left)\right)+ϵ\left(t\right).\hfill \end{array}$$

Equation (1) represents the dynamical model governing the evolution of state $\mathbf{x}$ through time, while $\eta$ is a Gaussian centered noise of covariance $\mathbf{Q}$ that models the process error. Equation (2) explains the relationship between the observation $\mathbf{y}\left(t\right)$ and the state to be estimated $\mathbf{x}\left(t\right)$ through the operator $\mathcal{H}$. The uncertainty of the observation model is represented by the $ϵ$ error, considered here to be Gaussian centered and of covariance $\mathbf{R}$. We assume that $ϵ$ and $\eta$ are independent and that $\mathbf{Q}$ and $\mathbf{R}$ are known. Two main approaches are generally considered for the mathematical resolution of the system (1) and (2), namely, variational data assimilation and stochastic data assimilation. They differ in the way they infer the analyzed state ${\mathbf{x}}^a$, the first is based on the minimization of a certain cost function while the latter aims to obtain an optimal a posteriori estimate. We encourage the reader to consider the book of [32] for detailed insights on the various aspects and methods of data assimilation.

A popular data assimilation algorithm that is largely used in the literature to grid sea level anomalies from along-track data is called Optimal Interpolation (OI) (e.g., [15,33]), this algorithm is also the method adopted in CMEMS altimetry product. Optimal Interpolation (OI) aims at finding the Best Linear Unbiased Estimator (BLUE) of a field $\mathbf{x}$ given irregularly sampled observations $\mathbf{y}$ in space and time and a background prior ${\mathbf{x}}^b$. The multivariate OI equations were derived in [34] for meteorology and numerous applications in oceanography have been reported since the early work of [16]. Supposing that the background state ${\mathbf{x}}^b$ has covariance $\mathbf{B}$, and the observation operator is linear $\mathcal{H}=\mathbf{H}$, the analyzed state ${\mathbf{x}}^a$ and the analyzed error covariance ${\mathbf{P}}^a$ can be calculated using the following OI set of equations:

$$\begin{array}{}\mathrm{(3)}& \mathbf{K}=\mathbf{B}\mathbf{H}{(\mathbf{R}+\mathbf{H}\mathbf{B}{\mathbf{H}}^T)}^{-1}\mathrm{called}\mathrm{the}\mathrm{Kalman}\mathrm{gain},\hfill \mathrm{(4)}& {\mathbf{x}}^a={\mathbf{x}}^b+\mathbf{K}(\mathbf{y}-\mathbf{H}{\mathbf{x}}^b),\hfill \mathrm{(5)}& {\mathbf{P}}^a=(\mathbf{I}-\mathbf{K}\mathbf{H})\mathbf{B}.\hfill \end{array}$$

It is worth mentioning that [35] showed that OI is closely related to the 3D-Var variational data assimilation algorithm which obtains ${\mathbf{x}}^a$ by minimizing the following cost function:

$$J\left(\mathbf{x}\right)={(\mathbf{x}-{\mathbf{x}}^b)}^T{\mathbf{B}}^{-1}(\mathbf{x}-{\mathbf{x}}^b)+{(\mathbf{y}-\mathbf{H}\mathbf{x})}^T{\mathbf{R}}^{-1}(\mathbf{y}-\mathbf{H}\mathbf{x}).$$

(6)

An important limitation of OI is that the Gaussian-like covariance priors lead to smoothing out the small-scale information (e.g., mesoscale eddies), more specifically, this is a limitation due to the use of a static climatological B matrix. For satellite-derived altimetry fields, this usually results in over-smoothing altimetry fields for structures below ≈100 km [18]. This limitation may be even more critical in the context of future high-resolution altimetry missions, which supports the development of new OI-based methods (e.g., multi-scale OI schemes as in [36]) or alternatives as addressed by our work.

3.2. Analog Data Assimilation

Endorsed by the recent development in data-driven methods and data storage capacities, the Analog Data Assimilation (AnDA) was introduced as an alternative to classical model-driven data assimilation under one or more of the following situations [13]:

• The model is inconsistent with observations.
• The cost of the model integration is high computationally.
• (mandatory) The availability of a represenative (large) dataset of the dynamics of the state variables to be estimated. These datasets are hereinafter called catalogs and denoted by $\mathcal{C}$. The catalog is organized in a two-column dictionary where each state of the system is associated with its successor in time, forming a set of couples $({\mathcal{A}}_i,{\mathcal{S}}_i)$ where ${\mathcal{A}}_i$ is called the analog and ${\mathcal{S}}_i$ its successor.

Given the considerations above, AnDA resorts to evaluating filtering, respectively smoothing, posterior likelihood, i.e., the distribution of the state to be estimated $\mathbf{x}\left(t\right)$ at time t, given past and current observations $\mathbf{y}(1,\dots ,t)$, respectively given all the available observation $\mathbf{y}(1,\dots ,T)$. This evaluation relies on the following state-space model:

$$\{\begin{array}{}\mathrm{(7)}& \mathbf{x}\left(t\right)=\mathcal{F}\left(\mathbf{x}\right(t-1\left)\right)+\eta \left(t\right),\hfill \mathrm{(8)}& \mathbf{y}\left(t\right)=\mathcal{H}\left(\mathbf{x}\right(t\left)\right)+ϵ\left(t\right).\hfill \end{array}$$

The difference between AnDA and classical data assimilation resides in the transition model Equation (7). The counterpart of a model-driven operator $\mathcal{M}$ of Equation (1) is here the operator $\mathcal{F}$ which refers to the considered data-driven operator, so called, the analog forecasting operator. This operator makes use of the available catalog $\mathcal{C}$ and assumes that the state forecast can be inferred from similar situations in the past.

Provided the definitions of the analogs and successors given above, the derivation of this operator resorts to characterizing the transition distribution i.e., $p\left(\mathbf{x}\right(t\left)\right|\mathbf{x}(t-1))$. Following [13], a Gaussian conditional distribution is adopted:

$$p\left(\mathbf{x}\right(t\left)\right|\mathbf{x}(t-1))=\mathcal{N}(\mu \left(\mathbf{x}\right(t-1\left)\right),\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right)),$$

(9)

where $\mathcal{N}$ is a Gaussian distribution of mean $\mu \left(\mathbf{x}\right(t-1\left)\right)$ and covariance $\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right)$. These parameters of the Gaussian distribution are calculated using the result of a K nearest neighbors search. The K nearest neighbors (analogs) ${\mathcal{A}}_{k\in (1,\dots ,K)}$ of state $\mathbf{x}(t-1)$ and their successors ${\mathcal{S}}_{k\in (1,\dots ,K)}$, along with a weight associated to each pair $({\mathcal{A}}_k,{\mathcal{S}}_k)$ are used to calculate $\mu \left(\mathbf{x}\right(t-1\left)\right)$ and $\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right)$ as we will show in the next paragraph, the forecast state $\mathbf{x}\left(t\right)$ is then sampled from $\mathcal{N}\left(\mu \right(\mathbf{x}(t-1)),\mathsf{\Sigma }(\mathbf{x}(t-1)\left)\right)$. The weights are defined using a Gaussian kernel ${\mathcal{K}}_G$.

$${\mathcal{K}}_G\left(u,v\right)=exp\left(-\frac{{\parallel u-v\parallel }^2}{\sigma }\right).$$

(10)

Scale parameter $\sigma$ is locally-adapted to the median value of the K distances $\parallel x(t-1)-{\mathcal{A}}_k{\parallel }^2$ to the K analogs. Other types of kernels might be considered (e.g., [4,10]), investigating kernel choice is out of the scope of this paper.

The mean and the covariance of the transition distribution might be calculated following several strategies. We consider in this work the three analog forecasting operators introduced in AnDA [13], more details can be found in Appendix A:

• Locally-constant operator: Mean $\mu \left(\mathbf{x}\right(t-1\left)\right)$ and covariance $\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right))$ are given by the weighted mean and covariance of the K successors ${\mathcal{S}}_{k\in (1,\dots ,K)}$.
• Locally-incremental operator: Here, the increments between the K analogs and their corresponding successors are calculated ${\mathcal{S}}_{k\in (1,\dots ,K)}-{\mathcal{A}}_{k\in (1,\dots ,K)}$. The weighted mean of the K increments is then added to the $\mathbf{x}(t-1)$ to obtain $\mu \left(\mathbf{x}\right(t-1\left)\right)$. While $\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right))$ results in the weighted covariance of these differences.
• Locally-linear operator: A weighted least-square estimation of the linear regression of the state at time t given the state at time $t-1$ is performed based on the K pairs $({\mathcal{A}}_k,{\mathcal{S}}_k)$. The parameters of the linear regression are then applied to state $\mathbf{x}(t-1)$ to obtain $\mu \left(\mathbf{x}\right(t-1\left)\right)$. Covariance $\mathsf{\Sigma }\left(\mathbf{x}\right(t-1\left)\right))$ is represented by the covariance of the residuals of the fitted weighted linear regression.

We may state clearly the key difference between the AnDA and reduced-rank Kalman filters and the OI method. It lies in the fact that the AnDA introduces a dynamical operator and not a prescribed space-time covariance model, and this dynamical operator is state-dependent and globally non-linear. The proposed analog forecasting operator can be seen as a state-dependent linear Gaussian operator, meaning that it is locally Gaussian and linear at each time step with a parameterization that depends on the current state, such that globally the dynamical operator is non-linear and non-Gaussian. A special case where AnDA is equivalent to OI is when all the elements of the catalog are considered as neighbors of any state vector. This case comes to assume that the dynamical operator is linear and state-independent. It is obviously of low interest due to the computational burden resulting from using all the catalog.

The application of the AnDA framework faces the curse of dimensionality i.e., the search of analogs is highly affected by the dimensionality of the problem and can fail at finding good analogs for state vector dimensions above 20 [13]. As proposed in [24], the extension of AnDA models to high-dimensional fields may then rely on turning the global assimilation issue into a series of lower-dimensional ones. We consider here an approach similar to [24] using a patch-based and EOF-based representation of the two-dimensional (2D) fields, i.e., the 2D fields are decomposed into a set of overlapping patches, each patch being projected onto an EOF space. Analog strategies then apply to patch-level time series in the EOF space.

Overall, as detailed in the following section, the proposed analog data assimilation model for SLA fields relies on three key components: A patch-based representation of the SLA fields, the selection of a kernel to retrieve analogs and the specification of a patch-level analog forecasting operator.

4. Analog Reconstruction for Altimeter-Derived SLA

4.1. Patch-Based State-Space Formulation

As stated above, OI may be considered as an efficient model-based method to recover large-scale structures of SLA fields. Following [24], this suggests considering the following two-scale additive decomposition:

$$X=\overline{X}+dX+\xi ,$$

(11)

where $\overline{X}$ is the large-scale component of the SLA field, typically issued from an optimal interpolation, $dX$ the fine-scale component of the SLA field we aim to reconstruct and $\xi$ is the remaining unresolved scales.

The reconstruction of field $dX$ involves a patch-based and EOF-based representation. It consists in regarding field $dX$ as a set of $P\times P$ overlapping patches (e.g., $2^{\circ }\times 2^{\circ }$), an example of patch locations is shown in Figure 2. This set of patches is referred to as $\mathcal{P}$, and we denote by ${\mathcal{P}}_s$ the patch centered on position s. After building a catalog ${\mathcal{C}}_{\mathcal{P}}$ of patches from the available dataset of residual fields $X-\overline{X}$ (see Section 3.2), we proceed to an EOF decomposition of each patch in the catalog. The reconstruction of field $dX({\mathcal{P}}_s,t)$ at time t is then stated as the analog assimilation of the coefficients of the EOF decomposition in the EOF space given an observation series in the patch space. Formally, $dX({\mathcal{P}}_s,t)$ decomposes as a linear combination of a number $N_E$ of EOF basis functions:

$$dX({\mathcal{P}}_s,t)=\sum _{k=1}^{N_E}{\alpha }_k(s,t)EO{F}_k,$$

(12)

with $EO{F}_k$ the $k^{th}$ EOF basis and ${\alpha }_k(s,t)$ the corresponding coefficient for patch ${\mathcal{P}}_s$ at time t. Let us denote by $\mathsf{\Phi }({\mathcal{P}}_s,t)$ the vector of the $N_E$ coefficients ${\alpha }_k(s,t)$: $\mathsf{\Phi }({\mathcal{P}}_s,t)=\{{\alpha }_1(s,t),\dots ,{\alpha }_{N_E}(s,t)\}$. This vector represents the projection of $dX({\mathcal{P}}_s,t)$ in the lower-dimensional EOF space.

4.2. Patch-Based Analog Dynamical Models

Given the considered patch-based representation of field $dX$, the proposed patch-based analog assimilation scheme involves a dynamical model stated in the EOF space. Formally, Equation (9) leads to the following Gaussian conditional distribution in the EOF space:

$$p\left(\mathsf{\Phi }({\mathcal{P}}_s,t)\right|\mathsf{\Phi }({\mathcal{P}}_s,t-1))=\mathcal{N}(\mu \left(\mathsf{\Phi }({\mathcal{P}}_s,t-1)\right),\mathsf{\Sigma }\left(\mathsf{\Phi }({\mathcal{P}}_s,t-1)\right)),$$

(13)

We consider the three analog forecasting operators presented in Section 3.2, namely, the locally-constant, the locally-incremental and the locally-linear. The calculation of the weights associated to each analog-successor pair relies on a Gaussian kernel ${\mathcal{K}}_G$ (Equation (10)). The search for analogs in the $N_E$-dimensional patch space (in practice, $N_E$ ranges from 5 to 20) ensures a better accuracy in the retrieval of relevant analogs compared to a direct search in the high-dimensional space of state $dX$. It also reduces the computational complexity of the proposed scheme.

Another important extension of the current study is the possibility of exploiting auxiliary variables with the state vector $\mathsf{\Phi }$ in the analog forecasting models. Such variables may be considered in the search for analogs as well as regression variables in a locally-linear analog setting. Regarding the targeted application to the reconstruction of SSH fields and the proposed two-scale decomposition (Equation (11)), two types of auxiliary variables seem to be of interest: The low-resolution component $\overline{X}$ to take into account inter-scale relationship [17], and Sea Surface Temperature (SST) with respect to the widely acknowledged SST-SSH synergies [17,19,21]. We also apply patch-level EOF-based decompositions to include both types of variables in the considered analog forecasting models (Equation (13)).

4.3. Numerical Resolution

Given the proposed analog assimilation model, the proposed scheme first relies on the creation of patch-level catalogs from the training dataset. This step requires the computation of a training dataset of fine scale data $d{X}_{training}$, this is done by subtracting a large-scale component ${\overline{X}}_{training}$ from the original training dataset. Here, we consider the large-scale component of training data to be the result of a global (By global, we mean here an EOF decomposition over the entire case study region, by contrast to the patch-level decomposition considered in the analog assimilation setting.) EOF-based reconstruction using a number of EOF components that retains 95% of the dataset variance, which accounts for horizontal scales up to ≈100 km. This global EOF-based decomposition provides a computationally-efficient means for defining large-scale component $\overline{X}$. This EOF-based decomposition step is followed by the extraction of overlapping patches for all variables of interest, namely ${\overline{X}}_{training}$, $d{X}_{training}$ and potential auxiliary variables, and the identification of the EOF basis functions from the resulting raw patch datasets. This leads to the creation of a patch-level catalog ${\mathcal{C}}_{\mathcal{P}}$ from the EOF-based representations of each patch.

Given the patch-level catalog, the algorithm applied for the mapping SLA fields from along-track data, referred to PB-AnDA (for Patch-Based AnDA), is stated in Algorithm 1 and a sketch of the method is shown in Figure 3.

Algorithm 1 Patch-Based AnDA

1: Compute the large-scale component $\overline{X}$, here, we consider the result of optimal interpolation (OI) projected onto the global EOF basis functions. 2: Split the case study region into overlapping $P\times P$ patches, here, $20\times 20$ patches 3: For each patch position s, use the Analog Ensemble Kalman Smoother (AnEnKS) [13], for patch ${\mathcal{P}}_s$ of field $dX$. As stated in (13), the assimilation is performed in the EOF space, i.e., for EOF decomposition $\Phi ({\mathcal{P}}_s,t)$, using the operator derived from EOF-based reconstruction (12) and decomposition (11) as observation model $\mathcal{H}$ in (8) and the patch-level training catalog described in the previous section. The assimilation is sequential and is performed each 3-days. 4: Reconstruct fields $dX$ from the set of assimilated patches ${\left\{dX({\mathcal{P}}_s,·)\right\}}_s$. This relies on a spatial averaging over overlapping patches (here, a 5-pixel overlapping in both directions). 5: the reconstruction of fields X as $\overline{X}+dX$.

We may point out two important aspects in the implementation of the proposed patch-level AnDA setting:

• (step 3) In the analog forecasting setting, the search for analogs is restricted to patch exemplars in the catalog within a local spatial neighborhood (typically a patch-level 8-neighborhood), except for patches along the seashore for which the search for analogs is restricted to patch exemplars at the same location.
• (step 4) In practice, we do not apply the patch-level assimilation to all grid positions. Consequently, the spatial averaging may result in blocky artifacts. We then apply a patchwise EOF-based decomposition-reconstruction with a smaller patch-size (here, $17\times 17$ patches) to remove these blocky artifacts.

5. Experimental Setting

We detail below the parameter setting of the models evaluated in the reported experiments, including the proposed PB-AnDA scheme:

• PB-AnDA: We consider $20\times 20$ patches with 15-dimensional EOF decompositions ($N_E=15$), which typically accounts for 99% of the data variance for the considered dataset. The postprocessing step exploits $17\times 17$ patches and a 15-dimensional EOF decomposition. Regarding the parametrization of the AnEnKS procedure, we experimentally cross-validated the number of nearest neighbors K to 50, the number of ensemble members $n_{ensemble}$ to 100 and the observation covariance error in Equation (8) is considered to be diagonal $\mathbf{R}={\kappa }^2\mathbf{I}$ and $\kappa =0.001$m.
• Optimal Interpolation: We apply an Optimal Interpolation to the processed along-track data. It provides the low-resolution component for the proposed PB-AnDA model and a model-driven reference for evaluation purposes. The background field is a null field. We use a Gaussian covariance model with a spatial correlation length of 100 km and a temporal correlation length of 15 days (±5 timesteps since our data is 3-daily). These choices result from a cross-validation experiment.
• VE-DINEOF: We apply a second state-of-the-art interpolation scheme using a data-driven strategy solely based on EOF decompositions, namely VE-DINEOF [36]. Using an iterative reconstruction scheme, VE-DINEOF starts by filling the missing data with a first guess, here along-track pseudo-observation field for along-track data positions and $\overline{X}$ for missing data positions. For each iteration, the resulting field is projected on the most significant EOF components calculated from the clean catalog data, then missing data positions are updated using pixels from the reconstructed new field. We run this iterative process until convergence. To make this algorithm comparable to the proposed PB-AnDA setting, we perform the reconstruction for each patch position then regroup the results as in PB-AnDA.
• G-AnDA: With a view to evaluating the relevance of the patch-based decomposition, we also apply AnDA at the region scale, referred to as G-AnDA. It relies on an EOF-based decomposition of the detail field $dX$. We use 150 EOF components, which accounts for more than 99% of the total variance of the SSH dataset. From cross-validation experiments, the associated AnEnKS procedure relies on a locally-linear analog forecasting model with $K=500$ analogs, $n_{ensemble}=100$ ensemble members and a diagonal observation covariance similar to as in PB-AnDA.

The patch-based experiments were run on Teralab infrastructure using a multi-core virtual machine (30 CPUs, 64G of RAM). We used the Python toolbox for patch-based analog data assimilation [24] (available at github.com/rfablet/PB_ANDA). Optimal Interpolation was implemented on Matlab using [36]. Throughout the experiments, two metrics are used to assess the performance of the considered interpolation methods: (i) Daily and mean Root Mean Square Error (RMSE) series between the reconstructed SLA fields X and the groundtruthed ones; (ii) daily and mean correlation coefficient between the fine-scale component $dX$ of the reconstructed SLA fields and of the groundtruthed ones. These two metrics allow a good evaluation on image reconstruction capabilities and are widely used in missing data interpolation literature [38,39].

6. Results

We evaluate the proposed PB-AnDA approach using the OSSE presented in Section 2. We perform a qualitative and quantitative comparison to state-of-the-art approaches. We first describe the parameter setting used for the PB-AnDA as well as benchmarked models, namely OI, an EOF-based approach [36] and a direct application of AnDA at the region level. We then report numerical experiments for noise-free and noisy observation data as well the relevance of auxiliary variables in the proposed PB-AnDA scheme.

6.1. SLA Reconstruction from Noise-Free Along-Track Data

We first performed an idealized noise-free experiment, where the along-track observations were not contaminated with noise. The interpolation performances for this experiment are illustrated in

Table 1. SLA Interpolation performance for a noise-free experiment: Root Mean Square Error (RMSE) (meters), correlation statistics for Optimal Interpolation (OI), VE-DINEOF, G-AnDA and PB-AnDA with regard to the groundtruthed SLA fields. The relative gain with regard to OI is also shown in percentage. See Section 5 for the corresponding parameter settings.

Criterion		RMSE	Correlation	$\frac{{\mathit{RMSE}}_{\mathit{OI}}-\mathit{RMSE}}{{\mathit{RMSE}}_{\mathit{OI}}}$
OI		0.026 ± 0.007	0.81 ± 0.08	-
VE-DINEOF		0.023 ± 0.007	0.85 ± 0.07	11.53%
G-AnDA		0.020 ± 0.006	0.89 ± 0.04	23.07%
PB-AnDA	Locally-constant	0.014 ± 0.005	0.95 ± 0.03	46.15%
	Locally-Increment	0.014 ± 0.005	0.95 ± 0.03	46.15%
	Locally-Linear	0.013 ± 0.005	0.96 ± 0.02	50.00%

. Our PB-AnDA algorithm significantly outperforms OI. More specifically, the locally-linear PB-AnDA results in the best reconstruction among the competing methods. We suggest that this improvement comes from the reconstruction of fine-scale features learned from the archived model simulation data. Figure 4a reports interpolated SSH fields and their gradient fields which further confirm our intuition. PB-AnDA interpolation shows an enhancement of the gradients and comes out with some fine-scale eddies that were smoothed out in OI and VE-DINEOF. This is also confirmed by the Fourier power spectrum of the interpolated SLA fields in Figure 4b.

6.2. SLA Reconstruction from Noisy Along-Track Data

We also evaluated the proposed approach for noisy along-track data. Here, we ran two experiments with an additive zero-mean Gaussian noise applied to the simulated along-track data. We considered a diagonal noise covariance of ${\gamma }^2\mathbf{I}$ where $\gamma =0.01$ m (Experiment A) and $\gamma =0.03$ m (Experiment B) which was more close to the instrumental error of conventional altimeters. Given the resulting noisy along-track dataset, we applied the same methods as for the noise-free case study.

We ran PB-AnDA using different values for $\kappa$. For Experiment A,

Table 2. Impact of standard variation of observation error $\kappa$ in AnDA interpolation performance using noisy along-track data ($\gamma$ = 0.01 m): RMSE (meters) of AnDA interpolation for different values of parameter $\kappa$. For the same dataset, OI RMSE is $0.039$.

$\kappa$	0.1	0.05	0.03	$0.01$	0.005	0.001	0.0001
$rms{e}_{PB-AnDA}$	0.035	0.030	0.028	$\mathbf{0.025}$	0.025	0.029	0.044

shows that the minimum is reached using the true value of the error $\kappa =0.01$ m. While for Experiment B,

Table 3. Impact of standard variation of observation error $\kappa$ in AnDA interpolation performance using noisy along-track data ($\gamma$ = 0.03 m): RMSE (meters) of AnDA interpolation for different values of parameter $\kappa$. For the same dataset, OI RMSE is $0.066$.

$\kappa$	0.1	0.05	0.03	$0.01$	0.005	0.001	0.0001
$rms{e}_{PB-AnDA}$	0.038	0.036	0.035	$\mathbf{0.0349}$	0.037	0.046	0.076

shows that the minimum is counter-intuitively reached again using value of the error $\kappa =0.01$ m with a negligible margin compared to the true value.

Our algorithm is then compared with the results of the application of the competing algorithms considered in this work. Results are shown in

Table 4. SLA Interpolation performance for noisy along-track data: RMSE (meters) and correlation statistics for OI, VE-DINEOF, G-AnDA and PB-AnDA w.r.t. the groundtruthed SLA fields. The relative gain with regard to OI is also shown in percentage. See Section 5 for the corresponding parameter settings.

	Criterion		RMSE	Correlation	$\frac{{\mathit{RMSE}}_{\mathit{OI}}-\mathit{RMSE}}{{\mathit{RMSE}}_{\mathit{OI}}}$
$\gamma$ = 0.01 m	OI		0.039 ± 0.005	0.64 ± 0.09	-
	VE-DINEOF		0.035 ± 0.005	0.68 ± 0.09	10.25%
	G-AnDA		0.030 ± 0.005	0.78 ± 0.06	23.07%
	PB-AnDA	Locally constant	0.026 ± 0.005	0.82 ± 0.05	33.33%
		Increment	0.028 ± 0.005	0.81 ± 0.05	28.20%
		Local Linear	0.0245 ± 0.005	0.83 ± 0.05	37.17%
$\gamma$ = 0.03 m	OI		0.066 ± 0.006	0.41 ± 0.09	-
	VE-DINEOF		0.060 ± 0.006	0.45 ± 0.09	9.09%
	G-AnDA		0.039 ± 0.006	0.67 ± 0.09	40.90%
	PB-AnDA	Locally constant	0.035 ± 0.006	0.688 ± 0.064	46.96%
		Increment	0.036 ± 0.006	0.656 ± 0.07	45.45%
		Local Linear	0.032 ± 0.006	0.708 ± 0.063	51.51%

. PB-AnDA still outperforms OI in terms of RMSE and correlation statistics in both experiments. The locally-linear version of PB-AnDA depicts the best reconstruction performance. We report an example of the reconstruction in Figure 5. Similarly to the noise-free case study, PB-AnDA better recovers finer-scale structures in Figure 5a compared with OI, VE-DINEOF and G-AnDA. In Figure 5b, PB-AnDA also better reconstructs a larger-scale North-East structure, poorly sampled by along-track data and hence poorly interpolated by OI.

6.3. PB-AnDA Models with Auxiliary Variables

We further explore the flexibility of the analog setting to the use of additional geophysical variable information as explained in Section 4.2. Intuitively, we expect SLA fields to involve inter-scale dependencies as well as synergies with other tracers [19,40]. The use of auxiliary variables provide the means for evaluating such dependencies and their potential impact on reconstruction performance. We consider two auxiliary variables that are used in the locally-linear analog forecasting model (7): (i) To account for the relationship between the large-scale and fine-scale component, we may consider variable $\overline{X}$; (ii) considering potential SST-SSH synergies, we consider SST fields. Overall, we consider four parameterization of the regression variables used in PB-AnDA: The sole use of $dX$ (PB-AnDA-$dX$), the joint use of $dX$ and SST fields (PB-AnDA-$dX$+SST), the joint use of $dX$ and $\overline{X}$ (PB-AnDA-dX+$\overline{X}$), the joint use of $dX$ and the groudntruthed version of $\overline{X}$ denoted by ${\overline{X}}^{GT}$, (PB-AnDA-dX+${\overline{X}}^{GT}$). The later provides a lower-bound for the reconstruction performance, assuming the low-resolution component is perfectly estimated.

We report mean RMSE (meters) and correlation statistics for these four PB-AnDA parameterizations in

Table 5. PB-AnDA reconstruction performance using noisy along-track data for different choices of the regression variables in the locally-linear analog forecasting model: PB-AnDA-$dX$ using solely $dX$, PB-AnDA-$dX+$ SST (Sea Surface Temperature) using both $dX$ and SST, PB-AnDA-$dX+\overline{X}$ using both $dX$ and $\overline{X}$, and PB-AnDA-$dX+{\overline{X}}^{GT}$ using $dX$ and the true large-scale component ${\overline{X}}^{GT}$. The table shows the RMSE (meters) and correlation statistics.

	PB-AnDA Model	RMSE	Correlation
$\gamma$ = 0.01 m	PB-AnDA-$dX$	0.025 ± 0.005	0.83 ± 0.05
	PB-AnDA-$dX$ + SST	0.024 ± 0.005	0.83 ± 0.05
	PB-AnDA-$dX+\overline{X}$	0.023 ± 0.005	0.84 ± 0.05
	PB-AnDA-$dX+{\overline{X}}^{GT}$	0.021 ± 0.004	0.87 ± 0.04
$\gamma$ = 0.03 m	PB-AnDA-$dX$	0.032 ± 0.006	0.708 ± 0.06
	PB-AnDA-$dX$ + SST	0.031 ± 0.006	0.710 ± 0.06
	PB-AnDA-$dX+\overline{X}$	0.029 ± 0.006	0.717 ± 0.06
	PB-AnDA-$dX+{\overline{X}}^{GT}$	0.026 ± 0.005	0.730 ± 0.05

for the noisy case-study. Considering PB-AnDA-$dX$ as reference, these results show a very slight improvement when complementing $dX$ with SST information. Though limited, we report a greater improvement when adding the low-resolution component $\overline{X}$. Interestingly, a significantly greater improvement is obtained when adding the true low-resolution information. The mean results are in accordance with [17], which reported that large-scale SLA information was more informative than SST to improve the reconstruction of the SLA at finer scales. Though mean statistics over one year leads to rather limited improvement, daily RMSE time series (Figure 6) reveal that for some periods, for instance between day 130 and 150, relative improvements in terms of RMSE may reach 10% with the additional information brought by the large-scale component. In this respect, it may noted that PB-AnDA-dX+$\overline{X}$ always perform better than PB-AnDA-dX. An example of the reconstruction in reported in Figure 7.

7. Discussion

Analog data assimilation can be regarded as a means to fuse ocean models and satellite-derived data. We regard this study as a proof-of-concept, which opens research avenues as well as new directions for operational oceanography. Our results advocate for complementary experiments at the global scale or in different ocean regions for a variety of dynamical situations with a view to further evaluating the relevance of the proposed analog assimilation framework. Such experiments should evaluate the sensitivity of the assimilation with respect to the size of the catalog. The scaling up to the global ocean also suggests investigating computationally-efficient implementation of the analog data assimilation. In this respect, the proposed patch-based framework intrinsically ensures high parallelization performance. From a methodological point of view, a relative weakness of the analog forecasting models (9) may be their low physical interpretation compared with physically-derived priors [18]. The combination of such physically-derived parameterizations to data-driven strategies appear to be a promising research direction. While we considered an OSSE to evaluate the proposed scheme, future work will investigate applications to real satellite-derived datasets, including the use of independent observation data such as surface drifters’ track data to further assess the performance of the proposed algorithm.

The analog method is at the heart of this work as it is appealing by its implementation ease and its intuitive strategy, but it must not be seen as the only data-driven method adapted to the framework we presented. As long as a data-driven forecasting operator can be derived, other data-driven methods can be investigated [41,42]. One promising path is the use of neural networks, as they sparked off a series of breakthroughs in other fields [42,43,44,45,46]. While neural network based approaches can lead to better performances due to their superior regressing capabilities, two clear advantages of adopting analog methods are the fact that they do not need a time consuming training phase and that analog methods are easy to understand and interpret compared to black-box approaches such as neural networks.

The Analog Data Assimilation method as used in this work relies on several hyperparameters, assumptions and design choices. These considerations are discussed in the following:

• In this work, we used all the available data for the creation of the catalog, we expect that a rich dataset is important in increasing the likelihood of finding good analogs, yet a more thorough study is needed to assess the impact of reducing the temporal resolution of the dataset versus reducing the amount of the total available data. To compensate the computational impact of using a large dataset for the search for analogs, we used the FLANN (Fast Library for Approximate Nearest Neighbors) library as in [26]. This method makes use of tree-indexing techniques and is suitable for this kind of high dimensional applications [47].
• While we used a conditional Gaussian distribution in Equation (9), another alternative is the use of a conditional multinomial distribution. This resorts to sampling one of the analogs’ successors as the forecast. Adopting this alternative would mean that we rely strongly on the archived catalog, so that forecasts of each ensemble member are actual elements of the catalog. This reduces the ability of the model to generate a rich variability of forecast scenes as done by the conditional Gaussian distribution.
• While we used a Gaussian kernel in Equation (10), other alternatives include cone kernels [10] which are more adapted to finding analogs in time series. The performance of our algorithm was slightly improved at the expense of more time execution and an additional hyperparameter to tune empirically, and we decided to prioritize simpler and efficient kernels. Regarding the scale parameter, the median was chosen due to its robustness to outliers.
• Although the AnEnKS was used in this work, a possible alternative is the use of an analog based Particle Smoother which drops Gaussian assumptions, however techniques based on particle filters need more ensemble members than their Ensemble Kalman Filter counterparts, thus causing a considerable increase in computational demands which was impractical for our application.
• We encourage the reader to refer to the discussion section in [24] for more insights about the rationale behind the use of patch-based representations and EOF-based dimensionality reduction.

Through the experiments conducted in this work, it was shown that the best performance was always reached using the locally-linear analog operator, which is in line with our previous findings [13,24]. An explanation for the superiority of this approach is that it better approximates locally the true dynamical model [48].

Beyond along-track altimeter data as considered in this study, future missions such as SWOT (NASA/CNES) promise an unprecedented coverage around the globe. More specifically, the large swath is expected to provide a large number of data, urging for the inspection of the potential improvements that this new mission will bring compared to classical along-track data. In the context of analog data assimilation, the interest of SWOT data may be two-fold. First, regarding observation model (8), SWOT mission will both significantly increase the number of available observation data and enable the definition of more complex observation models exploiting for instance velocity-based or vorticity-based criterion. Second, SWOT data might also be used to build representative patch-level catalogs of exemplars. Future work should investigate these two directions using simulated SWOT test-beds [49].

8. Materials and Methods

The experimental results presented in this work were obtained using the Python PB-AnDA toolbox made available by the authors at https://github.com/rfablet/PB_ANDA.

9. Conclusions

This work sheds light on the opportunities that data science methods are offering to improve altimetry in the era of big data. Assuming the availability of high-resolution numerical simulations, we show that Analog Data Assimilation (AnDA) can outperform the Optimal Interpolation method and retrieve smoothed out structures resulting from the sole use of OI both with idealized noise-free and more realistic noisy observations for the considered case study. Importantly, the reported experiments point out the relevance for combining OI for larger scales (above 100 km) whereas the proposed patch-based analog setting successfully applies to the finer-scale range below 100 km. This is in agreement with the recent application of the analog data assimilation to the reconstruction of cloud-free SST fields [24]. We also demonstrate that AnDA can embed complementary variables in a simple manner through the regression variables used in the locally-linear analog forecasting operator. In agreement with our recent analysis [17], we demonstrate that the additional use of large-scale SLA information may further improve the reconstruction performance for fine-scale structures.

We may state here the limitations of the present work and possible research avenues for the future. The experiments presented in this work were conducted on a numerical simulation derived dataset. A major future work direction would be then to apply the same procedure on real satellite-derived SLA contaminated with more complex noise models, then investigate the contribution of the use of numerical simulation datasets as catalogs. As combining multi-source datsets can also be challenging when using auxiliary variables relationships with SLA, an interesting experiment for example would be constructing a catalog with real SST observations combined with numerical simulation SLA datasets.

More efforts should be directed to assess the quality of the catalogs (spatio-temporal resolution, total years of measurements to consider, occurence of rare events, etc.). Besides, building a good catalog can represent an opportunity for the use of neural networks based methods, and confronting these powerful regressors to our method is a promising future step. We also note that PB-AnDA can be a relevant candidate for the interpolation of other geophysical variables (e.g., Sea Surface Salinity, Chlorophyll concentrations, etc.) under the condition that they verify the set of assumptions made in this work.