## 1 Introduction  
The spatio-temporal high resolution monitoring of sea surface geophysical parameters (e.g., temperature, salinity, ocean colour) is of key interest for a variety of scientific fields [1; 2; 3]. Direct observations of these geophysical tracers are provided by satellite remote sensing observations and in-situ networks. However, due to sensors characteristics (e.g., space-time sampling, sensor type) and their sensitivity to the atmospheric conditions (e.g., rain, clouds), only partial, with potentially high missing data rates, and possibly noisy observations are available. As a consequence, providing high resolution gape free spatio-temporal fields, in both space and time, based on these observations have long been a crucial challenge that motivated the development of several spatio-temporal interpolation tools. Within the satellite ocean community, Optimal Interpolation (OI) is a standard technique used in several operational products [4; 5; 6; 7; 8; 9; 10]. Given a covariance model of spatio-temporal dynamics, the interpolated field results from a linear combination of the observations. In general, stationary covariance hypotheses are considered, which prove relevant for the reconstruction of horizontal scales above 100 km. Fine scale components in the other hand may hardly be retrieved with such approaches and a variety of research studies aim to improve the reconstruction of high-resolution components of spatio-temporal fields. Empirical Orthogonal Function (EOF) based interpolation is an other category widely used in geosciences [11; 12; 13]. It relies on a Singular Value Decomposition (SVD) to compute an EOF basis, the field is then reconstructed by projecting the observations on the EOF subspace until a convergence criterion is reached [14]. Unfortunately, dealing with high missing data rates decreases the encoded variability in the EOF components witch results in smoothing fine scale structures. Data assimilation is the state-of-the-art framework for the reconstruction of dynamical systems from partial observations based on a given numerical model [15; 16]. Statistical data assimilation schemes, especially ensemble Kalman filters, have become particularly popular due to their trade-off between computational efficiency and modeling flexibility. Unlike OI and EOF based techniques, these schemes explicitly rely on dynamical priors to address interpolation issues resulting in better representation of fine scale components. However, When dealing with sea surface dynamics, the analytical derivation of these priors involves simplifying assumptions which may not be satisfied by real observations [17]. By contrast, realistic analytical parameterizations may lead to highly computationally-demanding numerical models associated with modeling and inversion uncertainties [18], which may limit their relevance for an application of the interpolation of a single sea surface tracer. Recently, data-driven approaches [11; 19] have emerged as relevant alternatives to model-driven schemes. They take benefit from the increasing availability of remote sensing observations and simulation data to derive computationally efficient [20] dynamical priors. Analog methods are one of the first data-driven techniques developed within a data assimilation framework [19]. In our recent study [20; 21; 22], we proved the relevance of such data-driven approache when addressing the spatio-temporal interpolation of sea surface geophysical tracers. Combining analog data assimilation (AnDA) with a patch-based representation have shown great results with respect to the state-of-the-art OI and EOF-based schemes. However, the parametrization of the proposed framework involves tuning several parameters principally due to the data-driven formulation of the dynamical prior based on analog forecasting.