Inferring flow energy, space scales, and timescales: freely drifting vs. fixed-point observations

A novel method for the inference of spatiotemporal decomposition of oceanic surface flow variability is presented and its performance assessed in a synthetic idealized configuration with horizontally divergentless flow. Inference methodology is designed for observations of surface velocity. The ability of networks of surface drifters and moorings to infer the spatiotemporal scales of surface ocean flow variability is quantified. The sensitivity of inference performance for both types of platforms to the number of observations, geometrical configurations, and flow regimes is presented. As drifters simultaneously sample spatial and temporal variability, they are shown to be able to capture both spatial and temporal flow properties even when deployed in isolation. Moorings are particularly adept for the characterization of the flow's temporal variability and may also capture spatial scales provided they are deployed as arrays. In particular, we show that our method correctly identifies whether drifters are preferentially sampling spatial vs. temporal variability. Pending further developments, this method opens novel avenues for the analysis of existing datasets as well as the design of future experimental campaigns targeting the characterization of small-scale (e.g., <100 km) ocean variability.

1 Introduction

Characterizing oceanic surface motions in terms of their spatial and temporal scales is a recognized pathway toward the identification of the numerous processes that occur in the ocean as well as toward an improved understanding of their occurrences, life cycle, interactions, and impact on other components of the ocean variability (Ferrari and Wunsch, 2009). For example, Arbic et al. (2014) relied on horizontal wavenumber-frequency decompositions in order to quantify and rationalize the impact of ocean mesoscale turbulence on longer-term ocean variability in idealized, realistic numerical simulations and altimetric observations. At higher frequencies, wavenumber-frequency decomposition enables the separation of internal gravity waves and balanced motions that share similar spatial scales and are therefore entangled in instantaneous two-dimensional datasets (Torres et al., 2019; Jones et al., 2023). For example, using a wavenumber-frequency decomposition, Qiu et al. (2018) were able to quantify the so-called “transition scale” above which altimetric observations are dominated by balanced turbulence and below which smaller scales are dominated by internal gravity waves. These decompositions are easily performed with numerical simulation output and are provided on complete and regular spatial and temporal grids. However, the lack of observational knowledge of the high-frequency and small-scale distribution of energy is a recognized limitation of the validation of tide-resolving kilometer-resolution global- or basin-scale numerical models of the ocean circulation (Arbic et al., 2018; Yu et al., 2019b; Arbic et al., 2022).

The characterization of ocean variability in terms of spatial and temporal scales is also relevant from an operational perspective. The description of an ocean variable's autocorrelation properties is required to map sparse observations via optimal interpolation (Bretherton et al., 1976; Bretherton and McWilliams, 1980). For instance, estimation of surface currents heavily relies on the accurate mapping of altimetric observations which consist of narrow (order 5 to 10 km) geographically and temporally distant tracks (Pujol et al., 2016). The advent of wide-swath altimetric (Morrow et al., 2019; Fu et al., 2024) and upcoming current measuring satellite missions introduces novel challenges regarding the mapping of the observed variables and the separation of slower balanced motions and faster internal gravity waves. This has motivated the development of novel strategies for the separation of the signatures associated with both classes of motion. These strategies rely on a priori knowledge of the motions' spatial and temporal scales (Barth et al., 2014, 2021; Ubelmann et al., 2021, 2022).

The in situ characterization of ocean variability at mesoscale to submesoscale (e.g., <100 km, <10 d) has been a central objective for a number of ambitious experimental efforts over the last decade: LatMIX (Shcherbina et al., 2015), Carthe Consortium (Poje et al., 2014; D'Asaro et al., 2018), OSMOSIS (Buckingham et al., 2016; Yu et al., 2019a), and SMODE (Farrar et al., 2020). Estimation of the time–space decomposition of upper-ocean variability has resulted from the dense dedicated mooring deployments of OSMOSIS and further highlighted difficulties associated with the Doppler shifting of small-scale structures when observed from fixed platforms (Callies et al., 2020). Such experiments incur significant financial and environmental costs; therefore any optimization in the experimental design and/or improved data analysis strategies are advantageous. Drifters are cheap and experimentally light platforms for the spatial and temporal characterization of ocean variability but require adequate inference methodologies. This study presents one such methodological development.

The characterization of horizontal and temporal variability of oceanic surface motions from observations represents a challenge that depends on the class of motions of interest, the quantity and nature of observations available, and the lack of a methodology that is both sufficiently versatile to the differing observation platforms and mathematically coherent. Fixed-point platforms provide information that is horizontally localized over potentially extended time periods with fine temporal resolution. Such data are conducive to temporal decomposition (Polzin and Lvov, 2011). The tracking of surface and subsurface drifting platforms provides ocean current observations, which are also amenable to temporal decomposition (Lumpkin et al., 2017), albeit representing the Lagrangian particle and thereby aliasing certain spatial characteristics. At daily to monthly timescales, drifters have enabled characterization of mesoscale eddy variability via inspection of surface current autocorrelation or spectral properties (Zhang et al., 2001; Lumpkin et al., 2002; Veneziani et al., 2004; Sykulski et al., 2016) or rotary wavelet decomposition (Lilly and Gascard, 2006; Lilly et al., 2011). The Global Drifter Program has collected surface current information worldwide for ∼30 years. Recently, the advent of GPS and wider-bandwidth satellite communications has enabled high-frequency sampling of surface drifter positions and a generation of surface drifter velocity datasets with global hourly coverage (Elipot et al., 2016). Over the last decade, global descriptions of the ocean surface high-frequency variability have emerged (Elipot et al., 2010, 2016; Yu et al., 2019b; Arbic et al., 2022). These descriptions are timely to validate recent kilometer-scale tide-resolving basin-scale numerical simulations (Arbic et al., 2018).

Satellite observations are well posed to characterize surface ocean spatial variability. The constellation of conventional nadir altimeters provide maps of sea level and surface currents that resolve larger mesoscale motions (Ballarotta et al., 2019). However, spatial and temporal gaps between nadir altimeters impose limitations on the resolvable spatial and temporal resolutions (Ballarotta et al., 2023). Consequently, there are multiple spatial and temporal characterizations of ocean variability that combine altimetry with other in situ datasets, e.g., moorings, XBTs (expendable bathythermographs), and tomography (Zang and Wunsch, 2001; Wunsch, 2010; Wortham and Wunsch, 2014). For smaller spatial scales, ship-based measurements of tracers and currents have informed the estimation of spatial scales of ocean variability (Callies and Ferrari, 2013), but such measurements potentially entangle spatial and temporal contributions to an unclear extent. Drifters are thought to offer promising data for the description of smaller mesoscale and submesoscale variability (Balwada et al., 2016, 2021). Dedicated experiments with deployments of a large number of surface drifters such as the experiment conducted by the Carthe Consortium have provided useful datasets to demonstrate small-scale ocean variability despite also highlighting potential biases associated with the horizontally divergent character of the flow at these scales (Poje et al., 2017; Pearson et al., 2019, 2020; Wang and Bühler, 2021).

Here we present a new method for the spatial and temporal characterization of oceanic surface flow variability. To test the method we consider an idealized configuration of ocean variability, whose properties and synthetic generation are described in Sect. 2.1. The novel method for the inference of the flow properties is described in Sect. 2.3. The inference is then applied to several scenarios of observations in order to explore the performance of the inference relative to the number of observations (Sect. 3.2), to platform spatial separation (Sect. 3.1), and to flow regime (Sect. 3.4). The results are discussed and conclusions drawn in Sect. 4.

2 Method

2.1 Idealized ocean surface flow design

We consider a two-dimensional and time-variable flow, described by the sum of rotational and divergent contributions:

$$\begin{array}{}\text{(1)}& {\displaystyle }u& {\displaystyle }=-{\partial }_y\mathit{\psi }+{\partial }_x\mathit{\varphi },\text{(2)}& {\displaystyle }v& {\displaystyle }={\partial }_x\mathit{\psi }+{\partial }_y\mathit{\varphi },\end{array}$$

where u and v are the zonal (toward positive x) and meridional (toward positive y) velocities in the respective directions x and y, ψ is the streamfunction, ϕ is the velocity potential, and ∂_x and ∂_y are the partial derivatives in x and y. We can describe the second-order behavior of ψ and ϕ, equivalently, by either their covariance functions or spectral densities. For general random fields a and b, defined over x, we define the stationary covariance function as $C_{ab}\left(\mathit{\tau }\right)=\langle a\left({\mathbf{x}}_{\mathrm{0}}\right),b({\mathbf{x}}_{\mathrm{0}}+\mathit{\tau })\rangle$, where the inner product is given as the covariance inner product $\langle a,b\rangle =\mathrm{E}\left[\right(a-\mathrm{E}\left[a\right]\left)\right(b-\mathrm{E}\left[b\right]\left)\right]$. Here, the boldface x₀ and τ denote a location and distance in x, respectively, in space and time. As stationarity is assumed, covariance is defined only as a function of τ. We define the corresponding spectral density as S_ab(ω), where the boldface ω represents a location in wavenumber and frequency space. As shown by Wiener–Khinchin's theorem, the covariance function and the spectral density are Fourier pairs, so that

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{C}_{ab}\left(\mathit{\tau }\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{S}_{ab}\left(\mathit{\omega }\right)\exp \left(\mathrm{2}\mathit{\pi }i\mathit{\omega }\mathit{\tau }\right)\mathrm{d}\mathit{\omega },\text{and}\\ {\displaystyle }& {\displaystyle }{S}_{ab}\left(\mathit{\omega }\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{C}_{ab}\left(\mathit{\tau }\right)\exp (-\mathrm{2}\mathit{\pi }i\mathit{\omega }\mathit{\tau })\mathrm{d}\mathit{\tau }.\end{array}\end{array}$$

Given an assumed parameterization of C_ψψ, C_ϕϕ, and C_ψϕ, the horizontal velocity auto- and cross-covariances are thus

$$\begin{array}{}\text{(4)}& {\displaystyle }{C}_{uu}& {\displaystyle }=-{\partial }_{yy}{C}_{\mathit{\psi }\mathit{\psi }}-{\partial }_{xx}{C}_{\mathit{\varphi }\mathit{\varphi }}+{\partial }_{xy}{C}_{\mathit{\varphi }\mathit{\psi }}+{\partial }_{xy}{C}_{\mathit{\psi }\mathit{\varphi }},\text{(5)}& {\displaystyle }{C}_{vv}& {\displaystyle }=-{\partial }_{xx}{C}_{\mathit{\psi }\mathit{\psi }}-{\partial }_{yy}{C}_{\mathit{\varphi }\mathit{\varphi }}-{\partial }_{xy}{C}_{\mathit{\varphi }\mathit{\psi }}-{\partial }_{xy}{C}_{\mathit{\psi }\mathit{\varphi }},\text{(6)}& {\displaystyle }{C}_{uv}& {\displaystyle }={\partial }_{xy}{C}_{\mathit{\psi }\mathit{\psi }}-{\partial }_{xy}{C}_{\mathit{\varphi }\mathit{\varphi }}+{\partial }_{yy}{C}_{\mathit{\psi }\mathit{\varphi }}-{\partial }_{xx}{C}_{\mathit{\varphi }\mathit{\psi }}.\end{array}$$

Similarly, given the spectral densities S_ψψ, S_ϕϕ, and S_ψϕ, we define the power and cross-power spectral densities of the horizontal velocities as

$$\begin{array}{}\text{(7)}& {\displaystyle }{S}_{uu}& {\displaystyle }=l^{\mathrm{2}}{S}_{\mathit{\psi }\mathit{\psi }}+k^{\mathrm{2}}{S}_{\mathit{\varphi }\mathit{\varphi }}-kl(S_{\mathit{\psi }\mathit{\varphi }}+S_{\mathit{\varphi }\mathit{\psi }}),\text{(8)}& {\displaystyle }{S}_{vv}& {\displaystyle }=k^{\mathrm{2}}{S}_{\mathit{\psi }\mathit{\psi }}+l^{\mathrm{2}}{S}_{\mathit{\varphi }\mathit{\varphi }}+kl(S_{\mathit{\psi }\mathit{\varphi }}+S_{\mathit{\varphi }\mathit{\psi }}),\text{(9)}& {\displaystyle }{S}_{uv}& {\displaystyle }=kl(S_{\mathit{\varphi }\mathit{\varphi }}-S_{\mathit{\psi }\mathit{\psi }})-k^{\mathrm{2}}{S}_{\mathit{\psi }\mathit{\varphi }}+l^{\mathrm{2}}{S}_{\mathit{\varphi }\mathit{\psi }},\end{array}$$

where k and l are horizontal wavenumbers. For our numerical experiment, we derive a purely rotational flow by setting ϕ=0 and, so, simply, $u=-{\partial }_y\mathit{\psi }$ and $v={\partial }_x\mathit{\psi }$. This leads to the covariance functions $C_{uu}=-{\partial }_{yy}{C}_{\mathit{\psi }\mathit{\psi }}$, $C_{vv}=-{\partial }_{xx}{C}_{\mathit{\psi }\mathit{\psi }}$, and $C_{uv}={\partial }_{xy}{C}_{\mathit{\psi }\mathit{\psi }}$ and spectral densities S_uu=l²S_ψψ, S_vv=k²S_ψψ, and $S_{uv}=-kl{S}_{\mathit{\psi }\mathit{\psi }}$.

To parameterize the flow we seek either a covariance function or spectral density that satisfies the physical requirements of the streamfunction ψ; namely, we require a log-linear decay in the high frequency or wavenumber of the spectral density. A good candidate for this is the isotropic Matérn covariance function (Rasmussen and Williams, 2005) with auto-covariance function and power spectral density (PSD),

$${\displaystyle }\begin{array}{rl}{\displaystyle }C\left(\mathit{\tau }\right)& {\displaystyle }={\displaystyle \frac{{\mathrm{2}}^{\mathrm{1}-\mathit{\nu }}}{\mathrm{\Gamma }\left(\mathit{\nu }\right)}}\left(\mathit{\lambda }{∥\mathit{\tau }∥}_{\mathrm{2}}{)}^{\mathit{\nu }}{\mathcal{K}}_{\mathit{\nu }}\right(\mathit{\lambda }{∥\mathit{\tau }∥}_{\mathrm{2}}),\text{and}\\ {\displaystyle }S\left(\mathit{\omega }\right)& {\displaystyle }={\displaystyle \frac{c_{\mathit{\nu }}}{{\left({∥\mathit{\omega }∥}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\lambda }}^{\mathrm{2}}\right)}^{\mathit{\nu }+\mathrm{D}/\mathrm{2}}}},\text{where}\\ {\displaystyle }& {\displaystyle }{c}_{\mathit{\nu }}={\displaystyle \frac{{\mathrm{2}}^{\mathrm{D}}{\mathit{\pi }}^{\mathrm{D}/\mathrm{2}}{\mathit{\lambda }}^{\mathrm{2}\mathit{\nu }}\mathrm{\Gamma }(\mathit{\nu }+\mathrm{D}/\mathrm{2})}{\mathrm{\Gamma }\left(\mathit{\nu }\right)}},\end{array}$$

where $\Vert \cdot {\Vert }_{\mathrm{2}}$ denotes the Euclidean norm/distance, D is the dimension of τ and ω, Γ(⋅) denotes the Gamma function, and 𝒦_ν is the modified Bessel function of the second kind of order ν≥0. For positive half-integers of ν, i.e., $\mathit{\nu }=p-\mathrm{1}/\mathrm{2}$, where $p\in {\mathbb{N}}^{+}$, 𝒦_ν has an analytical expression, otherwise it must be numerically calculated. We assume ψ to follow a separable Matérn process in space (D=2) and time (D=1), so that $C_{\mathit{\psi }\mathit{\psi }}\left(\mathit{\tau }\right)={\mathrm{\Psi }}^{\mathrm{2}}{C}_{\text{ss}}\left({\mathit{\tau }}_{\mathrm{d}}\right)\cdot C_{\text{tt}}\left({\mathit{\tau }}_{\mathrm{t}}\right)$, where Ψ is the standard deviation of the streamfunction, $\mathit{\tau }=[{\mathit{\tau }}_{\mathrm{d}},{\mathit{\tau }}_{\mathrm{t}}]$, where τ_d represents the isotropic distance in space and τ_t represents the time lag, and both C_ss(τ_d) and C_tt(τ_t) are specified as correlation functions, that is, $C_{\text{ss}}\left(\mathrm{0}\right)=C_{\text{tt}}\left(\mathrm{0}\right)=\mathrm{1}$. For the kernel defined over space C_ss(τ_d), we define the slope and decorrelation parameters ν_s and λ_s, respectively. For the kernel defined over time C_tt(τ_t), we define the slope and decorrelation parameters ν_t and λ_t, respectively. This separability assumption is a concession on realism, which substantially eases the computational cost of the flow generation step and is not expected to affect our evaluation of the inference performance (Wortham and Wunsch, 2014; De Marez et al., 2023). The covariance functions with respect to u and v are thus

$$\begin{array}{}\text{(10)}& {\displaystyle }{C}_{uu}\left(\mathit{\tau }\right)& {\displaystyle }=-{\mathrm{\Psi }}^{\mathrm{2}}{C}_{\text{tt}}\left({\mathit{\tau }}_{\mathrm{t}}\right)\cdot {\displaystyle \frac{y^{\mathrm{2}}{C}_{\text{ss}}^{\prime \prime }\left({\mathit{\tau }}_{\mathrm{d}}\right)+x^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}^{-\mathrm{1}}{C}_{\text{ss}}^{\prime }\left({\mathit{\tau }}_{\mathrm{d}}\right)}{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{2}}}},\text{(11)}& {\displaystyle }{C}_{vv}\left(\mathit{\tau }\right)& {\displaystyle }=-{\mathrm{\Psi }}^{\mathrm{2}}{C}_{\text{tt}}\left({\mathit{\tau }}_{\mathrm{t}}\right)\cdot {\displaystyle \frac{x^{\mathrm{2}}{C}_{\text{ss}}^{\prime \prime }\left({\mathit{\tau }}_{\mathrm{d}}\right)+y^{\mathrm{2}}{\mathit{\tau }}_{\mathrm{d}}^{-\mathrm{1}}{C}_{\text{ss}}^{\prime }\left({\mathit{\tau }}_{\mathrm{d}}\right)}{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{2}}}},\text{(12)}& {\displaystyle }{C}_{uv}\left(\mathit{\tau }\right)& {\displaystyle }={\mathrm{\Psi }}^{\mathrm{2}}{C}_{\text{tt}}\left({\mathit{\tau }}_{\mathrm{t}}\right)\cdot {\displaystyle \frac{xy\left(C_{\text{ss}}^{\prime \prime }\right({\mathit{\tau }}_{\mathrm{d}})-{\mathit{\tau }}_{\mathrm{d}}^{-\mathrm{1}}{C}_{\text{ss}}^{\prime }({\mathit{\tau }}_{\mathrm{d}}\left)\right)}{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{2}}}},\end{array}$$

where primes denote derivatives with respect to the horizontal distance τ_d.

2.2 Synthetic flow data generation

The amplitude of the streamfunction Ψ is related to the flow standard deviation U via $\mathrm{\Psi }=U{\mathit{\lambda }}_{\mathrm{s}}\sqrt{({\mathit{\nu }}_{\mathrm{s}}-\mathrm{1})/{\mathit{\nu }}_{\mathrm{s}}}$. The reference flow simulation is defined such as to be representative of moderately energetic mesoscale turbulence with U=0.1 m s⁻¹, λ_s=100 km, and λ_t=5 d (Fig. 1) (Ferrari and Wunsch, 2009). Matérn slope parameters are chosen to be ${\mathit{\nu }}_{\mathrm{t}}=\mathrm{1}/\mathrm{2}$ and ${\mathit{\nu }}_{\mathrm{s}}=\mathrm{5}/\mathrm{2}$, leading to a −2 temporal spectrum slope and a spatial isotropic spectral slope of −6. While the temporal spectral slope fits expectations, the spatial spectral slope is steeper by a value of 1 (20 %) compared to the value typical of quasi-geostrophic turbulence (Callies and Ferrari, 2013; Wortham and Wunsch, 2014). This concession to realism was made because it yields an analytical form for the Matérn covariance function, which alleviates the computational cost of the inference substantially. We reparameterize the covariance functions by Ψ=γλ_s, where γ is interpreted as the amplitude parameter of the horizontal velocity process; as well as interpretability, this has some computational benefits.

Figure 1 Overview of the inference input data for the reference scenario: (a) streamfunction snapshot in color overlaid with eight drifter tracks and moorings used for the inference, (b) x velocity time series of the drifter identified by the black track in (a), and (c) x velocity time series at the mooring indicated by the black star in (a). In (b) and (c), black dots indicate the 2 d of sub-sampled data used in the inference.

With the previous choice of parameters, the streamfunction is generated over a 1000 km by 1000 km domain with 2 km grid spacing and over 100 d with hourly resolution (Fig. 1). This resolution is a factor of ∼50 times smaller than decorrelation, which is considered enough to resolve the synthesized variability and mitigate numerical interpolation errors in Lagrangian numerical simulations. Sizes of the spatial domain and the time series are ∼10 and ∼20 times larger than decorrelation scales, which ensures we are capturing multiple, effectively independent, realizations of the process.

The hourly synthetic flow is fed to the Parcels Python library configured with fourth-order Runge–Kutta time stepping and the default A-grid interpolation scheme in order to produce synthetic drifter trajectories (Delandmeter and van Sebille, 2019). Drifters are released initially at all flow grid points with the exception of a 20 km strip around all boundaries, amounting to a total of 9216 drifters for each simulation. Trajectories reaching domain boundaries are de-activated and not advected further in time and discarded from the list of observations that will be used for inference. The fraction of trajectories discarded was 52 % in the reference configuration. Drifter positions are stored at hourly resolution and velocities estimated from drifter positions with second-order finite differencing. An example of drifter trajectories is shown in Fig. 1. The flow amplitude averaged over time and space is about 1.8 % larger than that computed from drifter trajectories which reveals small turbophoresis, i.e., concentration of drifters in areas of lower energy (Freeland et al., 1975).

A non-dimensional parameter used to characterize flow is $\mathit{\alpha }=U{\mathit{\lambda }}_{\mathrm{t}}/{\mathit{\lambda }}_{\mathrm{s}}$. This parameter is expected to control the relative importance of spatial vs. temporal variability in the projection onto Lagrangian time series (Middleton, 1985). In the reference scenario, the value of α is 0.4, which is in the range of observed ocean values (Lumpkin et al., 2002). In order to obtain mooring and drifter time series with different α values, the synthetic flow is simply rescaled, and new Lagrangian trajectories are simulated with the rescaled flow.

2.3 Inference

Observed data y are composed of flow time series collected over time by N_p drifters or moorings to which a white noise n of standard deviation σ is added. The critical difference between drifter and mooring observations is that they are collected along drifter trajectories in the former case, i.e., u[x(t)]+n(t), where x(t) is a drifter trajectory, while they are collected at a fixed location in the latter one, i.e., $\mathbf{u}[\mathbf{x},t]+\mathbf{n}\left(t\right)$ where x is a mooring location.

We treat the collection of parameters $\mathbf{\Theta }=\mathit{\{}\mathit{\gamma },{\mathit{\lambda }}_{\mathrm{s}},{\mathit{\lambda }}_{\mathrm{t}},{\mathit{\sigma }}^{\mathrm{2}}\mathit{\}}$ as uncertain and unknown and probabilistically quantify this uncertainty. We treat Θ as a random variable and so naturally adopt the Bayesian paradigm of probability. Bayes' theorem states $\mathrm{p}(\mathbf{\Theta }\mid \mathbf{y})\propto \mathrm{p}(\mathbf{y}\mid \mathbf{\Theta })\mathrm{p}\left(\mathbf{\Theta }\right)$, where p(Θ∣y) is the posterior distribution, p(y∣Θ) is the likelihood, and p(Θ) is the prior distribution. The posterior is our target quantity and describes the probability distribution of Θ conditioned on the observed data. The likelihood is a probability distribution that assesses the probability of the data being generated, conditioned on some value of Θ. Finally, the prior represents our knowledge of Θ before we observe the data y; in this term we may include the results from previous analyses, bounds on values that Θ may take or any physically derived structure between the constituent parameters inside of Θ. Prior distributions are chosen here to be uniform between 0 and 10 times true parameter values.

Exact computation of p(Θ∣y) is analytically achievable for a small class of model problems; however, this is typically not so, and so p(Θ∣y) is computed numerically using Markov chain Monte Carlo (MCMC), as this is the gold standard in statistical computing. MCMC generates a dependent chain of draws from the posterior p(Θ∣y), such that subsets of Θ are visited proportionally to the posterior probability of the subsets. We show an example of this in Fig. 2 for the moored data reference scenario. The trace plots consist of 20 000 dependent samples from which we may derive summaries of the posterior distribution, p(Θ∣y), via standard Monte Carlo methods. For example, the marginal distributions of each parameter are represented by the histograms in the right-hand column of Fig. 2. MCMC is asymptotically exact in that the sampled draws converge to the exact posterior probability distribution. We generate samples using Metropolis–Hastings (MH), a well-known and accessible MCMC algorithm. Description and particulars are provided in the Appendix.

Figure 2 Trace plots of MCMC sampling for the flow parameters (left) and associated histograms (right) for a single inference based on eight drifter trajectories for the reference scenario. True parameter values are indicated by the black lines, while MAP location and values are indicated by thick gray lines.

As discussed above, we parameterize our model using the Matérn covariance function as it exemplifies a number of desirable physical characteristics. However, the derivative of the Matérn covariance function is difficult to obtain due to 𝒦_ν(⋅): analytical derivatives are only available at integer values of $\mathit{\nu }-\mathrm{1}/\mathrm{2}$, and numerical calculations of 𝒦_ν(⋅) are not available in any symbolic toolboxes that we are aware of. To mitigate the computational burden and enable the performance of ensemble of statistical experiments, we decided to fix the slope parameters ν_s and ν_t to half-integer values described in Sect. 2.2 and exclude these parameters from those inferred. This choice is relaxed in Sect. 3.5 in order to demonstrate that the inference of these parameters is possible as achieved in the purely temporal domain by Sykulski et al. (2016).

The noise chosen here (σ=0.01 m s⁻¹) is representative of oceanographic velocity observations. State-of-the-art surface drifters are for instance equipped with GPS and provide hourly position observations with 10 m accuracy (Poulain et al., 2022). Under the assumption that the noise is white, this leads to a noise standard deviation on velocity observations of about 0.004 m s⁻¹, which is less than half of what we considered. Acoustic Doppler current profilers (ADCPs), such as those deployed on moorings, exhibit typical accuracies in the 0.01 to 0.1 m s⁻¹ range (Klema et al., 2020). Our expectation is that the performance of the inference should eventually decay as the noise on velocity observations increases, and this should be verified in a future study.

2.4 Validation of the inference methodology

As the mooring and drifter data are simulated, we know the ground truth and so may validate the MCMC sampling methodology. We show this for two cases: first, we show the probabilistic parameter estimates from the reference flow (Sect. 2.1), and second, we compare the maximum a posteriori (MAP) estimates, i.e., $\widehat{\mathbf{\Theta }}:=\arg {max}_{\mathrm{\Theta }}\mathit{\left\{}\mathrm{p}\right(\mathbf{\Theta }\mid \mathbf{y}\left)\mathit{\right\}}$, of an 100-member ensemble with their true values. Examining a single scenario demonstrates the inherent uncertainty associated with a single experiment, whereas, inference across an ensemble looks at the variability that arises between data samples. In all cases, the data comprise a bivariate u and v time series collected either along eight trajectories (drifters) or at eight stationary locations (moorings), with 2 d temporal sampling over 100 d, amounting to 400 data points. The ensemble data are generated from the single spatiotemporal field with randomly sampled drifter tracks and mooring locations. Figure 2 shows the marginal posterior probability distributions of the single-member reference scenario. For all parameters, the true values lie well within the probability distribution. Note, σ² is not well resolved; this is because of the roughness of the Matérn process at the set sampling interval (see Fig. 1), confounds with the noise signal so that both processes may be viable in explaining the observed data. This is somewhat expected, and more detailed statistical diagnostics accompany the code in the Supplement. Figure 3 plots a histogram of the MAP values calculated from each ensemble member's MCMC chain against the true value. This shows the variability of the distributions about the true value over the ensemble. Again, all distributions are centered on the true values, and there exists some difficulty in observing σ² with precision. The precision of the inference will also be quantified by the difference between the third and first quartiles, which will be referred to as the interquartile width (IQW).

Figure 3 Distribution of parameters MAP values for the reference flow and reference observation scenario (scenario REF). True parameter values are represented by vertical black lines. First and third quartiles are vertical dashed gray lines and provide insight into the interquartile width (IQW).

2.5 Inference scenarios

This study reports on the performance of the inference method under several scenarios (summarized in Table 1):

• REF corresponds to the nominal configuration described in Sect. 2.4 with eight simultaneously deployed platforms.

• Under SEP[dx], when multiple platforms are simultaneously sampling the flow, the separation between platforms and more generally their geometrical distribution is expected to modulate the performance of the inference. To simplify the analysis, we restrict the configuration to two simultaneous observing platforms (e.g., two drifters or two moorings) and investigate the sensitivity of the inference performance to their separation (with 10 % tolerance). For drifters, the separation is the initial one between the two drifters.

• In the IND scenario, inference is performed by assuming time series from different platforms are independent from each other. Such a situation would occur if individual moorings/drifters were deployed at the same location, but at times sufficiently far apart, no correlation is expected across the velocity time series recorded by each platform. In effect this amounts to quantifying the ability of one platform at capturing flow parameters and investigating the sensitivity to the length of the time series.

• Under OPT[N_p], platforms are deployed in a spiral configuration that leads to separations that span the flow spatial decorrelation scale (see Sect. A2). The purpose of this experiment is to perform a simple experimental design optimization of the number of platforms deployed and of the choice between moorings and drifters.

• In REG[α], the amplitude of the flow is rescaled in order to explore different values of the flow regime parameter $\mathit{\alpha }=U{\mathit{\lambda }}_{\mathrm{t}}/{\mathit{\lambda }}_{\mathrm{s}}$. The amplitude of the noise is linearly scaled as a function of α in order to maintain a fixed signal-to-noise ratio. Inference is performed with a single platform.

• The NU scenario is similar to OPT[N_p] with N_p=8, with the exception that the spectral slope parameters ν_s and ν_t are also inferred.

Table 1 Inference scenarios. All other parameters are held constant across the scenarios.

3 Results

3.1 Platform separation sensitivity

Under scenario SEP[dx], estimations of the flow amplitude are comparable for observations from two moorings or two drifters and precise, with IQWs lower than 13 % of true amplitudes and no sensitivity to separation (Fig. 4a). We argue this follows from the fact that inferences are provided with velocity observations as inputs. Drifter inferences of the flow amplitude exhibit a 1 % to 3 % low bias, which is comparable to that associated with turbophoresis (Sect. 2.2).

Figure 4 Sensitivity of parameter MAP estimates to platform separation (in km) for the two platform configurations (scenario SEP[dx]). Lines represent the median, while shaded areas are bounded by first and third quartiles. True values are in black. When visible, gray shadings represent the no-go zone of the prior and inference parameter exploration.

Mooring spatial-scale estimates are sensitive to separation (Fig. 4b). After a modest decrease in performance of the inference with separation as measured by IQWs, the best inference is obtained for separation in the range of 40 to 80 km. For larger separations, the inference precision decreases, with IQW reaching values of about 0 % of true values at 300 km, i.e., 3 times the flow spatial scale. This loss of performance with separation reflects the loss of correlation between the flow measured by each mooring and thus the lack of information about spatial structure in the dataset. Drifters exhibit no clear sensitivity to separation for the spatial-scale estimate. This may first be explained by the substantial displacements of the drifters compared to the separations considered. A flow exponentially autocorrelated over 10 d with a standard deviation of 10 cm s⁻¹ leads to an absolute dispersion of (250 km)² (Gurarie et al., 2017). The natural ability of drifters to explore space and time and therefore constrain spatial scales (see Sect. 3.2) provides a second explanation.

Mooring and drifter inferences of the flow temporal scale both exhibit a modest high bias of 5 % to 10 % (Fig. 4c). As expected, drifters are overall less effective than moorings at estimating the flow temporal-scale parameter. IQWs associated with drifter inferences are systematically larger than those associated with moorings, which fluctuate around 30 % for moorings, compared with drifters, which increase with separation up to 60 %.

3.2 Sensitivity to the number of independent platforms

Under scenario IND[N_p], single moorings (i.e., N_p=1) provide estimates of the flow amplitude, γ, and temporal decorrelation scale, λ_t, parameters that are precise, with IQW starting at about 16 % and 45 % of true values, respectively. Parameters γ and λ_t converge to true values as the number of independent moorings is increased (Fig. 5). For the maximum number of platforms considered, the IQW of the flow amplitude and temporal decorrelation has decreased to 4 % and 11 %, respectively. As expected from their inability to explore the spatial dimension, single moorings are however globally unable to capture the flow spatial scale, with IQW comparable to half the width of the parameter space allowed to be explored, i.e., [0,1000 km], which amounts to the prior uncertainty (that is, there is no resolution of uncertainty).

Figure 5 Sensitivity of parameter MAP estimates to the number of platforms (scenario IND[N_p]). Platforms are assumed independent from each other. Same representation as Fig. 4.

In comparison, drifters provide reasonable estimates of all three flow parameters (γ, λ_s, λ_t), with IQW starting at about 14 %, 92 %, and 95 % for one platform. These estimates converge toward truth as the number of platforms is increased, with IQW smaller than 16 % for all three parameters with 16 drifters. The ability of drifters to capture both spatial and temporal scales is explained by their ability to sample space and time simultaneously. MAP medians indicate mild biases with an underestimation of amplitude and overestimation of temporal decorrelation scale, which decrease as the number of drifters is increased. The amplitude low bias is about 7 % with a single drifter and reduces to about 1.4 % with 16 drifters, which is comparable to the turbophoresis bias (Sect. 2.2). The temporal decorrelation scale λ_t of drifters is always less accurate than that obtained with moorings, which we interpret as the price to pay for the simultaneous sampling of spatial and temporal variability.

3.3 Experimental design optimization

Optimizing an experimental design is a complex task that results from a compromise between scientific goals, a priori knowledge of the variables to be measured, financial and logistical constraints, and the need for redundancy, among other aspects. Scenario OPT[N_p] illustrates how one could identify what the minimum experimental design is, enabling an accurate estimation of flow properties.

Consistent with the results of the previous scenarios, no substantial bias is observed. IQW is used to quantify accuracy and therefore is the target variable to minimize to identify optimal design (Fig. 6). Apart from the one platform configuration, where the mooring is unable to estimate the spatial scale of variability, moorings and drifters present comparable sensitivities as a function of the number of platforms. The number of platforms required to reach a target IQW of 20 % of the true value for all parameters except for σ is 4 for both platforms (Fig. 6).

Figure 6 Sensitivity of parameter MAP estimates to the number of platforms (scenario OPT[N_p]). Same representation as Fig. 4.

In light of the low cost of drifters compared to moorings (factor of about 100 for deep sea applications), this result is particularly striking. However, we note the simplicity of the present exercise (idealized flow, constrained geometry of deployment; see Sect. A2) in light of past efforts on the matter (Bretherton and McWilliams, 1980; Barth and Wunsch, 1990). As stated in the preamble, optimizing for characterization of flow properties constitutes one consideration among many that may be taken in an experimental design optimization. Scientific goals may in general go well beyond the characterization of flow properties. If flow properties are suspected to evolve temporally, the use of drifters which are expected to eventually disperse will require multiple deployments in the area of interest unlike with moorings.

3.4 Flow regime sensitivity

We turn now to an investigation of the sensitivity of inferences to the flow parameter α (scenario REG[α]). We revert to the single-platform configuration in order to limit the exchange of information across platforms and the resulting constraint it brings for inference which may mask the α sensitivity. For comparison purposes we also perform a “time-only” inference of drifter velocity time series, which estimates flow amplitude, temporal decorrelation scale, and noise only and not the spatial decorrelation scale λ_s.

As anticipated from Sect. 3.2, inferences of flow amplitude and temporal decorrelation scale from mooring observations are relatively accurate, with IQW of about 15 % and 50 % of true values, respectively (Fig. 7a and c). The amplitude inference reflects the linear sensitivity to α. Spatial scales remain undetermined for all α values (Fig. 7b). This lack of sensitivity is expected due to exclusive sampling of temporal variability by a single mooring.

Figure 7 Sensitivity of parameter MAP estimates to flow regime α for the single-platform configuration (scenario REG[α]). Time-only drifter inference is in red in (a) (median MAP dashed) and (c) (quartiles and median). Same representation as Fig. 4 otherwise.

Inferences of flow amplitude from drifter observations are comparable to mooring inferences in terms of IQW albeit with a low bias of about 5 % (Fig. 7a). A comparable bias is observed on time-only inferences for small α values but is exacerbated for α larger than unity where it reaches about 35 % of the true amplitude (Fig. 7c). For large α, distortions of the temporal spectrum shape are likely affecting the overall performance of the time-only inferences which rely on the spectral distribution following that of a Matérn $\mathrm{1}/\mathrm{2}$ process.

For small α values (<0.2), inference of the flow spatial decorrelation scale from drifter observations are the worst, and the IQW is nearly comparable to those from mooring observations (Fig. 7b). Drifters indeed merely move over a flow timescale comparable to the spatial decorrelation scale in this flow regime, which has been historically coined a “fixed float” and can be effectively considered a mooring (Middleton, 1985; Lumpkin et al., 2002). Accordingly, when α<0.2, estimates of the flow amplitude γ and temporal decorrelation scale λ_t are comparable for moorings and for drifters whether with the standard inference or the “time-only” inference.

For larger values of α (e.g., >0.2), the precision of the flow spatial decorrelation-scale inference from drifter observations improves substantially with decreasing IQW (down to 50 % at α∼1). In contrast, estimates of the temporal decorrelation scale deteriorate, with a bias high of about 25 % and IQW width of about 120 %. At these values of α, the flow is in the so-called “frozen turbulence” regime, and drifters are in effect experiencing the spatial variability of the flow field (Middleton, 1985; Lumpkin et al., 2002). This is directly reflected in the estimate of the temporal scale obtained from the “time-only” inference, which monotonically decreases with α. The fact that the temporal scale from the space–time inference does not follow a similar trend is a testimony to the relevance of the latter method, which is able to identify that observations reflect a predominance of spatial variability and attribute reasonable space-scale and timescale estimates, albeit with moderate error and bias.

3.5 Spectral slope estimation

For the final experiment (NU), the assumption that spectral slopes are known is relaxed, and Matérn slope parameters ν_s and ν_t are inferred along with the other parameters, i.e., γ, λ_s, λ_t, and σ. The assumed prior distributions are uniform over [1,5] and [0,5] for ν_s and ν_t, which is larger than typical uncertainties in the ocean for these parameters. Estimating these parameters leads to a 45-fold increase in computing time, due to computation of the Bessel function 𝒦_ν, as discussed in Sect. 2.1.

The impact on flow parameter estimation is a modest increase of normalized IQW (Fig. 8) compared to OPT[N_p] with N_p=8 (Fig. 6). For instance, spatial and temporal decorrelation scales IQW estimated with mooring observations increase from 7 % and 14 % to 10 % and 18 % of true values, respectively. Inferences from drifter observations undergo comparable increases.

Figure 8 Distribution of parameter MAP values for the reference flow and reference observation scenario with inference of Matérn spectral slopes (scenario NU). True parameter values are represented by vertical black lines. First and third quartiles of the posterior distribution are vertical dashed gray lines and provide insight into IQW.

The inference of spatial and temporal Matérn slopes is successful with posterior distributions centered around their true values and IQW of less than 22 % of true values. Independent experiments with random platform deployments similar to REF lead to more contrasted results, with the temporal slope being effectively resolved but not the spatial slope (not shown). This is an indication that the estimation of Matérn slope parameters is more demanding on observation quality and information content. Pending improvements in the performance of the inference computation, these results present promising perspectives for the systematic inference of spectral slopes.

4 Conclusions

We have presented a novel Bayesian method to infer surface ocean circulation spectral parameters (e.g., amplitude and spatial and temporal decorrelation scales) from sparse observations of the flow. The intention was to quantify parameter uncertainty due to sampling and flow regimes. These results may guide the design and analysis of future field campaigns and open novel avenues for the analysis of existing datasets. We considered flow observation from two platforms typically employed in oceanography: moorings which provide fixed-point flow observations and drifters that provide along-flow flow observations. Inference based on both types of platforms provides flow characterization estimates that converge to true values as the number of observations is increased. The performance of the method was quantified in various observing configurations, which allowed us to highlight the pros and cons of each type of platform. As already recognized, moorings are well suited to characterizing temporal scales of variability and if deployed as appropriately spaced simultaneous networks can constrain flow spatial scales. Drifters naturally sample both space and time, and we showed they can simultaneously constrain and separate the flow's space scales and timescales even when deployed in isolation, which is the first demonstration to our knowledge. We also showed that the ability of drifter observations to characterize flow properties depends on a non-dimensional parameter that quantifies the relative magnitude of the spatial and temporal decorrelation scales. Given the relative low cost and low environmental impact associated with drifter deployments compared with moorings, we argue they provide a powerful and more sustainable means to characterize surface flow properties. Finally, the present inference method may be of more general relevance for the simultaneous space–time characterization of flow properties in other fluid mechanics configurations (Reneuve and Chevillard, 2020).

More developments are required in order to make this method applicable to realistic oceanographic configurations. First the method needs to be extended to flows that are composed of a superposition of processes commonly occurring in the ocean, e.g., internal waves and tides, near-inertial waves. Such an extension will present methodological challenges associated with the parameterization of the space–time variability associated with these processes. The assumption of space–time separability, which was imposed here by the selected method of flow field generation, may have to be relaxed in a realistic configuration (Wortham and Wunsch, 2014; De Marez et al., 2023). As long as correlations may be expressed in physical space, extension of the inference to non-separable cases is direct. The issue introduced by non-separable kernels is arguably the most difficult challenge that needs to be addressed. It may also be useful to generalize the inference method to simultaneously account for observations that are of diverse nature, for instance current observations from drifters, pressure from moorings, and sea level observations from satellite altimetry. Such an extension will require deriving the expected correlation between each of the variables concerned and will in any case depend on the process modeled. A first application of the method to real data may be with gridded altimetric sea level or current data (AVISO+, 2024). The smoothing applied to generate these products may allow alleviation of the complexity associated with high-frequency processes.

Moving to a more realistic flow configuration will require evolving the synthetic flow strategy. The present choice allowed us to generate flows with arbitrary spatiotemporal structure, including some flows that are unlikely to occur in the ocean, in order to enable a broad exploration of the inference performance. This approach could be pushed further with the superposition of multiple processes and non-separable kernels and will likely require leveraging spectral domain approaches. As highlighted in Sect. 3.5, there are some computational difficulties with estimating the spectral slope via the Matérn covariance function. Slope estimation in the spectral domain is simple as the slope appears in the PSD in an analytically tractable form (see Sykulski et al., 2016); however, for drifter-based inference, as we are interested in estimation of the Eulerian properties, we cannot use such Lagrangian spectral techniques. There are some recent results that resolve the computational burden imparted by the calculation of the Bessel function and its derivatives (Geoga et al., 2022). Regrettably, at the time of writing, code for this study's methodology is not widely available across coding platforms. We hope that this, or similar methodological advancements, may be included in future work that will focus on estimating more realistic flows. Using flows generated from dynamical models (quasi-geostrophic, primitive equations) may eventually be necessary to capture regimes of variability more closely, representative of the actual ocean dynamics, with more realistic representations of process life cycles. Finally, the advection of drifters could account for a stochastic component in order to represent inaccuracies associated with the dynamical system considered (Mínguez et al., 2012).

Applications of the inference method to realistic observation datasets (e.g., velocity observations from the Global Drifter Program – Lumpkin et al., 2017) would be computationally prevented in the present form due to the use of dense covariance arrays. Alleviating this constraint will require us to leverage sparsity in the inference inputs associated with observations that are distant in space and/or time compared to associated decorrelation scales. This represents the second-most-important challenge that will have to be faced for applications in realistic configurations. Data collected from regional campaigns may be more suitable in the short term.

Appendix A

A1 MCMC sampling

A1.1 Metropolis–Hastings algorithm

The Markovian property of MCMC implies that a sample Θ^[i] only depends on its previous sample Θ^[i−1]; the method by which Θ^[i] is generated from Θ^[i−1] distinguishes the various MCMC algorithms. All MCMC algorithms propose some Θ^[*] from Θ^[i−1] and with probability α either accept Θ^[*], in which case ${\mathbf{\Theta }}^{\left[i\right]}={\mathbf{\Theta }}^{[*]}$, or reject Θ^[*], in which case ${\mathbf{\Theta }}^{\left[i\right]}={\mathbf{\Theta }}^{[i-\mathrm{1}]}$. The Metropolis–Hastings (MH) algorithm, initially presented in Metropolis et al. (1953) and later extended by Hastings (1970), generates a proposal Θ^[*] from Θ^[i−1] using some user-specified proposal distribution $f({\mathbf{\Theta }}^{[*]}\mid {\mathbf{\Theta }}^{[i-\mathrm{1}]})$. Given a proposal Θ^[*], we accept the sample with probability r, where

$$\begin{array}{}\text{(A1)}& r=min\left(\mathrm{1},{\displaystyle \frac{\mathrm{p}({\mathbf{\Theta }}^{[*]}\mid \mathbf{y})f({\mathbf{\Theta }}^{[i-\mathrm{1}]}\mid {\mathbf{\Theta }}^{[*]})}{\mathrm{p}({\mathbf{\Theta }}^{[i-\mathrm{1}]}\mid \mathbf{y})f({\mathbf{\Theta }}^{[*]}\mid {\mathbf{\Theta }}^{[i-\mathrm{1}]})}}\right).\end{array}$$

Figure A1 (a) Illustration of an array of N_p=16 platforms for L=200 km and $\mathit{\beta }=\mathit{\pi }/\mathrm{3}$ used in OPT[N_p]. (b) Corresponding distribution of platform separations.

If the proposal density is symmetrical, that is, $f({\mathbf{\Theta }}^{[i-\mathrm{1}]}\mid {\mathbf{\Theta }}^{[*]})=f({\mathbf{\Theta }}^{[*]}\mid {\mathbf{\Theta }}^{[i-\mathrm{1}]})$, then Eq. (A1) reduces to the ratio of the posterior densities, and so the MH algorithm will always accept a proposed Θ^[*] that is more probable than Θ^[i−1]. The choice of $f(\cdot \mid \cdot )$ is critical to the success of the MH algorithm. If $f(\cdot \mid \cdot )$ is too wide, then the algorithm can become stuck for many iterations, thus generating very few unique proposals. Conversely, if $f(\cdot \mid \cdot )$ is too narrow, the algorithm will not effectively explore the parameter space; the sampled ${\mathbf{\Theta }}^{\left[\mathrm{1}\right]},\mathrm{\dots },{\mathbf{\Theta }}^{\left[n\right]}$ will be highly correlated; and again, few independent samples will be generated. One of the main drawbacks of the MH algorithm is that there are sampling parameters that need to be hand-tuned.

We parameterize $f(\cdot \mid \cdot )$ as a multivariate normal distribution with mean Θ^[i−1] and diagonal covariance matrix. A widely agreed upon rule of thumb to balance exploration and exploitation of the posterior distribution is an acceptance probability of ∼0.25. Accordingly, we set the standard deviations of the proposal distribution to be between 0.05 and 0.2 of the true parameter values, corresponding to situations where we have larger and lower instances of observed data. The reason for this is simple: as the number of observations increases, the uncertainty of our parameter values decreases, implying a tighter posterior distribution. Consequently, a tighter proposal distribution is required to achieve a comparable acceptance probability. Full validation results to guarantee fit and convergence of the MCMC estimation algorithm are presented alongside the code at https://github.com/apatlpo/nwastats (last access: 31 July 2024).

A1.2 Notes on alternative MCMC sampling algorithms

Modern MCMC algorithms have been dominated by gradient-based proposal methods, where a proposal Θ^[*] is generated by assessing the local topology surrounding Θ^[i−1]: this allows the algorithm to efficiently trade off notions of exploration and exploitation of the posterior distribution. Included in these algorithms are the popular Hamiltonian Monte Carlo techniques, such as those implemented in Stan (Carpenter et al., 2017), PyMC3 (Salvatier et al., 2016), and Pyro (Bingham et al., 2019); these implementations, as well as others such as GPJax (Pinder and Dodd, 2022), will typically use symbolic toolboxes to define the local topology of the posterior. Alternative MCMC algorithms should not affect the accuracy of the posterior estimation; but rather, they will differ in their sampling efficiency. This study is concerned with inference, and not operationalization, and so we choose the Metropolis–Hastings algorithm so as to avoid the issue of gradients at the cost of some hand tuning of the algorithm.

A2 Platform array design

For the experiment OPT[N_p], platforms are deployed at locations that aim to span a wide range of platform separations around some expectation of the spatial-scale decorrelation. For that purpose, locations were set along a spiral defined by its spatial footprint L, orientation β, and center (x_c,y_c) according to

$$\begin{array}{}\text{(A2)}& {\displaystyle }{x}_j+i{y}_j& {\displaystyle }=x_{\mathrm{c}}+i{y}_{\mathrm{c}}+r{\mathit{\theta }}_j\times e^{i({\mathit{\theta }}_j+\mathit{\beta })},\text{ with}\text{(A3)}& {\displaystyle }{\mathit{\theta }}_j=& {\displaystyle }j\times \mathit{\delta },\text{ and, }r=\left\{\begin{array}{ll}L/\mathit{\delta },& \text{if }{N}_{\mathrm{p}}=\mathrm{2}\\ L/\left(\mathrm{2}{N}_{\mathrm{p}}\mathit{\delta }\right),& \text{otherwise}\end{array}\right.\end{array}$$

where $\mathrm{0}\le j<N_{\mathrm{p}}$ is a platform digit identifier. We have made the choice $\mathit{\delta }=\mathit{\pi }/\mathrm{3}$. An illustration of such platform deployment is illustrated in Fig. A1a. The distribution of platform separations successfully spans the ensemble of length scales up to L (Fig. A1b).

Each draw in the OPT[N_p] ensemble experiment is based upon uniform random draws of the spiral center within the domain, of the spatial footprint L within [50 km,300 km], and of the orientation β within [0,2π].

Appendix B: Notations

Inferred parameters
γ	streamfunction amplitude to spatial- decorrelation-scale ratio
λs	spatial decorrelation scale
λt	temporal decorrelation scale
σ	noise standard deviation
νs	spatial slope parameter, only inferred in Sect. 3.5
νt	temporal slope parameter, only inferred in Sect. 3.5
Θ	vector composed of all inferred parameters
Other parameters
U	flow amplitude
Ψ	streamfunction amplitude
$\mathit{\alpha }=U{\mathit{\lambda }}_{\mathrm{t}}/{\mathit{\lambda }}_{\mathrm{s}}$	non-dimensional flow parameter
Np	number of observing platforms (e.g., drifters or moorings)
Variables
x,y	spatial coordinate or increment
t	temporal coordinate or increment
u,v	horizontal velocity field
ψ	streamfunction
ϕ	flow potential
Cab	cross-correlation between variables a and b
Sab	cross-spectrum between variables a and b
(k,l)	horizontal wavenumbers
𝒦ν	modified Bessel function of the second kind of order ν
Γ	Gamma function
y	observation vector