Inferring ocean transport statistics with probabilistic neural networks
Using a probabilistic neural network and Lagrangian observations from the
Global Drifter Program, we model the single particle transition probability
density function (pdf) of ocean surface drifters. The transition pdf is
represented by a Gaussian mixture whose parameters (weights, means and
covariances) are continuous functions of latitude and longitude determined to
maximise the likelihood of observed drifter trajectories. This provides a
comprehensive description of drifter dynamics allowing for the simulation of
drifter trajectories and the estimation of a wealth of dynamical statistics
without the need to revisit the raw data. As examples, we compute global
estimates of mean displacements over four days and lateral diffusivity. We use
a probabilistic scoring rule to compare our model to commonly used transition
matrix models. Our model outperforms others globally and in three specific
regions. A drifter release experiment simulated using our model shows the
emergence of concentrated clusters in the subtropical gyres, in agreement with
previous studies on the formation of garbage patches. An advantage of the
neural network model is that it provides a continuous-in-space representation
and avoids the need to discretise space, overcoming the challenges of dealing
with nonuniform data. Our approach, which embraces data-driven probabilistic
modelling, is applicable to many other problems in fluid dynamics and
oceanography.
Martin T. Brolly
# Inferring ocean transport statistics with probabilistic neural networks

###### Abstract

Using a probabilistic neural network and Lagrangian observations from the Global Drifter Program, we model the single particle transition probability density function (pdf) of ocean surface drifters. The transition pdf is represented by a Gaussian mixture whose parameters (weights, means and covariances) are continuous functions of latitude and longitude determined to maximise the likelihood of observed drifter trajectories. This provides a comprehensive description of drifter dynamics allowing for the simulation of drifter trajectories and the estimation of a wealth of dynamical statistics without the need to revisit the raw data. As examples, we compute global estimates of mean displacements over four days and lateral diffusivity. We use a probabilistic scoring rule to compare our model to commonly used transition matrix models. Our model outperforms others globally and in three specific regions. A drifter release experiment simulated using our model shows the emergence of concentrated clusters in the subtropical gyres, in agreement with previous studies on the formation of garbage patches. An advantage of the neural network model is that it provides a continuous-in-space representation and avoids the need to discretise space, overcoming the challenges of dealing with nonuniform data. Our approach, which embraces data-driven probabilistic modelling, is applicable to many other problems in fluid dynamics and oceanography.

## 1 Introduction

The motion of turbulent fluids can be characterised usefully by dynamical statistics such as dispersion, energy spectra and velocity structure functions (e.g., Batchelor 1953, Monin & Yaglom 1971). In oceanography much effort has been directed towards inferring such statistics from observations (e.g., LaCasce 2008, van Sebille et al. 2018). In many cases, these inference tasks can be related to problems in conditional probability density estimation. For example, estimating single-particle dispersion is related to estimating the conditional density

\[p\left(\mathbf{X}(t+\tau)-\mathbf{X}(t)\mid\mathbf{X}(t),\,t,\,\tau\right), \tag{1}\]

where \(\mathbf{X}(t)\) is the position of a particle at time \(t\), in that the dispersion is the variance of this distribution. Similarly, the velocity structure functions are moments of the conditional
density

\[p\left(\mathbf{u}(\mathbf{x}_{1},\,t)-\mathbf{u}(\mathbf{x}_{2},\,t)\mid\mathbf{x}_{1},\, \mathbf{x}_{2},\,t\right), \tag{2}\]

where \(\mathbf{u}(\mathbf{x},\,t)\) is the fluid velocity at position \(\mathbf{x}\). By estimating full conditional densities like (1) and (2), it is possible to estimate simultaneously a number of related statistics. For instance, (1) describes entirely the single-particle displacement statistics, while (2) encodes velocity structure functions of all orders, providing two-point Eulerian velocity statistics. It is no surprise, then, that estimating these conditional densities accurately is a nontrivial task.

In this work we consider a particular tool for conditional density estimation, the mixture density network (MDN) (Bishop, 1994), and test its performance in learning fluid statistics from observations. MDNs are machine learning models, which combine artificial neural networks with probabilistic mixture models to represent conditional densities (Bishop, 2006). Their use has increased rapidly in recent years with applications in a variety of fields for a range of reduced order modelling and emulation tasks, including surrogate modelling of fluid flow (Maulik et al., 2020), parameterisation of subgrid momentum forcing in ocean models (Guillaumin and Zanna, 2021), emulation of complex stochastic models in epidemiology (Davis et al., 2020) and multi-scale models of chemical reaction networks (Bortolussi and Palmieri, 2018), and subgrid scale closures in large eddy simulations of turbulent combustion (Shin et al., 2021).

We focus on learning the single-particle transition density (1) in the ocean near-surface using Lagrangian trajectory data collected as part of the Global Drifter Program (Lumpkin and Centurioni, 2019). A model of the transition density provides, at every point in the ocean, a probabilistic forecast for drifter displacements from that location. We show that the MDN model outperforms existing stochastic models of drifter dynamics based on Ulam's method (Ulam, 1960; Froyland, 2001), as well as another simple benchmark model, and eliminates the difficulty of designing appropriate discretisations of space needed for such models.

From the transition density it is possible to derive estimates of a range of single-particle statistics. As examples, we provide maps of the mean displacement over four days as a function of initial position \(\mathbf{X}_{0}\), as well as the lateral diffusivity. The transition density produces highly non-Gaussian statistics in some regions. By calculating the Kullback-Leibler divergence between our full model and a simplified Gaussian model, we quantify and map non-Gaussianity in drifter displacements.

The MDN model also provides the basis for a discrete-time Markov process model of drifter dynamics, offering a continuous space alternative to Markov chain models which have been used in numerous studies (Maximenko et al., 2012; van Sebille et al., 2012; Miron et al., 2017, 2021). We perform a global simulation of drifters for a period of ten years with initial positions given on a uniform grid, and reproduce the 'garbage patches' in subtropical gyres seen in previous studies.

The article is structured as follows. In SS2 we discuss conditional density estimation and the estimation of conditional statistics. In SS3 we introduce MDNs. In SS4 we describe the MDN model of the single-particle transition density from drifter observations. We compare its performance with alternative models, present derived single-particle statistics and simulate the clustering of drifters in subtropical gyres. In SS5 we conclude and suggest further problems where MDNs may be a useful tool.
Conditional modelling

While the aim of regression is to model \(\mathbb{E}[\mathbf{Y}\mid\mathbf{X}]\), where \(\mathbf{X}\) and \(\mathbf{Y}\) are random variables, conditional modelling (or conditional density estimation, CDE) is the task of inferring the full conditional probability density \(p(\mathbf{Y}\mid\mathbf{X})\)1. By modelling conditional densities, rather than just conditional means, we incorporate information about the variability of \(\mathbf{Y}\mid\mathbf{X}\); more than this, conditional models can capture skewness, excess kurtosis and multimodality. This comprehensive description of conditional statistics is valuable in applications where single point-estimates are insufficient due to inherent variability, and where there is interest in non-Gaussian statistics, including those associated with rare events. Conditional models can be used in two ways: (i) as stochastic surrogate models (or emulators), and (ii) as a tool for estimating conditional statistics.

Footnote 1: We restrict attention to the case of continuous random variables.

Parametric conditional models (such as MDNs) assume that, for each possible value of \(\mathbf{X}\), the distribution of \(\mathbf{Y}\mid\mathbf{X}\) belongs to a certain family of parametric distributions, i.e.

\[p(\mathbf{Y}\mid\mathbf{X})=\rho(\mathbf{Y}\,;\,\mathbf{\theta}(\mathbf{X})), \tag{3}\]

where \(\rho(\,\cdot\,;\,\mathbf{\theta})\) is the probability density corresponding to a family of distributions parameterised by \(\mathbf{\theta}\). In this case, not only must the form of \(\rho\) be chosen, but the dependence on the conditioned variable must also be modelled by some representation of \(\mathbf{\theta}(\mathbf{X})\).

### Estimating conditional statistics

Given data \(\{\mathbf{X}_{i},\mathbf{Y}_{i}\}\), a standard approach to estimating conditional statistics \(\mathbb{E}[\mathbf{f}(\mathbf{Y})\mid\mathbf{X}]\) is to first discretise (or 'bin') in \(\mathbf{X}\) and produce local estimates \(\mathbb{E}[\overline{\mathbf{f}(\mathbf{Y})}](\tilde{\mathbf{X}})\) for each value of the discretised variable \(\tilde{\mathbf{X}}\), typically by Monte Carlo estimation, such that

\[\mathbb{E}[\overline{\mathbf{f}(\mathbf{Y})}](\tilde{\mathbf{X}}):=\frac{\sum_{i}\mathbf{f}( \mathbf{Y}_{i})\,\,\mathbb{I}_{B}(\mathbf{X}_{i})}{\sum_{i}\,\mathbb{I}_{B}(\mathbf{X}_{i })}, \tag{4}\]

where \(\mathbb{I}_{B}\) is the indicator function of \(B\), the set of values of \(\mathbf{X}\) whose discretised value is \(\tilde{\mathbf{X}}\). For estimates (4) to be useful, one must design a suitable discretisation of the domain of \(\mathbf{X}\), which balances the need to choose a fine enough discretisation to resolve details in \(\mathbf{X}\) with the need to take sufficiently large bins to have enough data for these estimates to have reasonably small variance. This can be especially challenging when data is sparse, or when the density of data is highly inhomogeneous.

Conditional modelling offers an alternative approach wherein one first constructs a model of the conditional density, as in (3), that is continuous in both \(\mathbf{X}\) and \(\mathbf{Y}\), then computes estimates

\[\mathbb{E}_{\mathcal{M}}[\mathbf{f}(\mathbf{Y})\mid\mathbf{X}]:=\int\mathbf{f}(\mathbf{Y})\,\, \rho(\mathbf{Y}\,;\,\mathbf{\theta}(\mathbf{X}))\,\mathrm{d}\mathbf{Y} \tag{5}\]

for as many statistics as desired at any value of \(\mathbf{X}\) in the domain, without the need to revisit the raw data. In some cases the expectations \(\mathbb{E}_{\mathcal{M}}\) can be calculated using a closed-form expression. Where no such expression is known, the expectation can be computed by numerical integration or a Monte Carlo method. Since these calculations rely only on evaluating the modelled conditional density, or sampling from it, they are not limited by sparsity of data. Also, for a given \(\mathbf{X}^{*}\), estimates of the form (4) are informed only by observations in the same
bin as \(\mathbf{X}^{*}\), whereas in a conditional model, all observations are used to fit \(\rho(\mathbf{Y}\,;\,\mathbf{\theta}(\mathbf{X}^{*}))\). The schematic in figure 1 contrasts the standard approach and the conditional modelling approaches.

## 3 Mixture density networks

A mixture density network (Bishop 1994, 2006) is a conditional model where an artificial neural network is employed to represent the function \(\mathbf{\theta}(\mathbf{X})\) in (3) and the parametric form \(\rho(\,\cdot\,;\,\mathbf{\theta})\) corresponds to a mixture distribution. The density of a general mixture distribution is

\[\rho(\,\cdot\,;\,\,\mathbf{\theta})=\sum_{i=1}^{N_{\mathrm{c}}}\alpha_{i}\ \rho_{i}(\,\cdot\,;\,\,\mathbf{\theta}_{i}), \tag{6}\]

where \(N_{\mathrm{c}}\) is the number of components in the mixture, the \(i^{\mathrm{th}}\) component has density \(\rho_{i}(\,\cdot\,;\,\,\mathbf{\theta}_{i})\) with parameters \(\mathbf{\theta}_{i}\), \(\mathbf{\theta}=[(\alpha_{1},\ \mathbf{\theta}_{1}),\ \cdots,\ (\alpha_{N_{\mathrm{c}}},\ \mathbf{\theta}_{N_{\mathrm{c}}})]\) and the \(\alpha_{i}\) are component weights subject to the constraint

\[\sum_{i=1}^{N_{\mathrm{c}}}\alpha_{i}=1. \tag{7}\]

Commonly, the component densities \(\rho_{i}\) are chosen from the same family and, in particular, Gaussian, but components can be chosen differently. In the Gaussian case, the \(\mathbf{\theta}_{i}\) are conditional means and covariances.

The neural network representation of \(\mathbf{\theta}(\,\cdot\,)\) is itself parametric with parameters \(\mathbf{w}\); hence, MDNs model \(p(\mathbf{Y}\,\mid\,\mathbf{X})\) with \(\rho(\mathbf{Y}\,;\,\mathbf{\theta}(\mathbf{X}\,;\,\mathbf{w}))\). The network can have any architecture,

Figure 1: Estimating conditional statistics by a standard approach versus by first constructing a model for the conditional density.
but that of a multilayer perceptron (Rumelhart et al., 1985) (also known as a fully connected multilayer feedforward neural network) with nonlinear activation functions is common -- in this case \(\mathbf{w}\) consists of the weights and biases.

A natural loss function for conditional models, which quantifies how well they fit data, is the negative (conditional) log likelihood of observations \(\mathcal{D}=\{\mathbf{X}_{i},\mathbf{Y}_{i}\}\) under the model. In MDNs this is

\[\mathcal{L}(\mathbf{w}\,;\,\mathcal{D})=\sum_{i}-\log\rho\left(\mathbf{Y}_{i}\,;\,\bm {\theta}(\mathbf{X}_{i}\,;\,\mathbf{w})\right). \tag{8}\]

Training an MDN then amounts to finding optimal values for the neural network's parameters

\[\mathbf{w}^{*}=\operatorname*{arg\,min}_{\mathbf{w}}\;\mathcal{L}(\mathbf{w}\,;\,\mathcal{ D}). \tag{9}\]

Minimising the negative log likelihood is equivalent to maximising the log likelihood of training data, also referred to as the log score in probabilistic forecasting (Bernardo, 1979; Gneiting and Raftery, 2007; Brocker and Smith, 2007). Maximum likelihood estimation in this context differs from the more familiar setting of fitting an unconditional model for \(p(\mathbf{Y})\) given observed data \(\{\mathbf{Y}_{i}\}\) -- here, there is generically only one observed value of \(\mathbf{Y}\mid\mathbf{X}\) corresponding to each observed value of \(\mathbf{X}\), and for most values of \(\mathbf{X}\) there are no observations at all. It is clear, then, that, for each value of \(\mathbf{X}\), we are certainly not in the large-data regime that would allow one to invoke asymptotic properties of maximum likelihood estimates. The quality of parametric conditional models (3) depends critically on how well \(\mathbf{\theta}(\mathbf{X}\,;\,\mathbf{w}^{*})\) represents how the distribution of \(\mathbf{Y}\mid\mathbf{X}\) varies with \(\mathbf{X}\). In particular, since MDNs employ a neural network to model \(\mathbf{\theta}(\mathbf{X})\), and neural networks are highly flexible models, it is common for MDNs to exhibit poor generalisation unless regularisation techniques are used. In the following section we employ a widely used regularisation technique known as early stopping (see e.g. Prechelt (2012)), wherein a small proportion of training data (referred to as the test set) are not used to inform steps in the optimisation scheme, but are instead used to track the evolution of an estimate of the model's generalisation error (the value of the loss function evaluated on data outside the training set). The guiding heuristic is that it is typical for the generalisation error of neural networks to reach a minimum as training progresses before increasing due to overfitting -- early stopping is a strategy where one terminates model training when the generalisation error is believed to have reached this minimum. Details of our implementation are given in the following section.

## 4 Application to single-particle statistics of the ocean near-surface

In this section we present an MDN model of the single-particle transition density (1) of ocean surface drifting buoys (drifters). The model's parameters are inferred from trajectory data collected as part of the Global Drifter Program (Lumpkin and Centurioni, 2019).

### Data

We use the Global Drifter Program quality-controlled 6-hour interpolated dataset, which includes positions (latitude and longitude) and sea-surface temperatures. Drifter velocity estimates are also provided, though these are obtained by simple finite-differencing of position
in time and subject to error. The raw measurements are treated according to the procedure of Hansen & Poulain (1996), which involves the removal of suspected spurious values and interpolation to regular 6-hour intervals. The interpolation method, which is a form of kriging (Hansen & Herman 1989), assumes contamination by an uncorrelated zero-mean noise and makes assumptions about the structure functions of the discretised position process. We leave as a caveat to our results that this preprocessing of the data could be questioned and proceed taking the interpolated data as our ground-truth. Only position observations are used in our modelling. Figure 2 shows how many observed displacements are recorded per squared kilometre in each \(1^{\circ}\) latitude \(\times 1^{\circ}\) longitude square. These data were recorded between 1989 and 2021 and include a total of 23893 drriter trajectories. We split the data in two parts, by selecting approximately half (11946) of the drriter trajectories at random to use for creating the model and set the remaining data aside for validation. The overall dataset contains over 18 million observations of 6-hour displacements.

In section 4.3 we perform a model comparison. Skill scores are computed for the full training and validation datasets with global coverage, as well as for three restricted regions, labelled \(A\), \(B\), and \(C\), shown in figure 3, having extents \(20-50^{\circ}\) W, \(30-50^{\circ}\) N; \(145-175^{\circ}\) E, \(20-40^{\circ}\) N; and \(110-130^{\circ}\) W, \(10^{\circ}\) S-\(10^{\circ}\) N.

### Model

The transition density is not modelled in the most general form. Instead, we (i) consider, at first, a fixed value of the time-lag \(\tau\), so that the transition density may be written

\[p(\mathbf{X}_{n+1}\mid\mathbf{X}_{n}), \tag{10}\]

where \(\mathbf{X}_{n}=\mathbf{X}(t_{0}+n\tau)\), and (ii) assume the process \(\mathbf{X}(t)\) is time-homogeneous, such that (1) is independent of the initial time \(t\), (10) is independent of \(n\) and

\[p(\Delta\mathbf{X}\mid\mathbf{X}_{0}) \tag{11}\]

Figure 2: Count of drriter observations per squared kilometre.
\(\Delta\mathbf{X}\) is the displacement of a drifter from its position at the previous timestep, denoted \(\mathbf{X}_{0}\). By assuming time-homogeneity we neglect the effects of seasonality and low-frequency variability in ocean dynamics. If, additionally, an assumption of Markovianity is made, then (10) is enough to construct a discrete-time Markov process model (\(\mathbf{X}_{n}\)) for drifter position (Pavliotis, 2014). For a Markov assumption to be accurate, the discretisation timescale \(\tau\) must be chosen appropriately. We choose a timescale of 4 days on the basis that the Lagrangian velocity decorrelation time (or integral timescale) at the surface was previously estimated from drifters to be approximately 2-3 days in all four ocean basins (Rupolo, 2007).

We model (11) using an MDN \(-\) see the schematic in figure 4. The model takes as input \(\mathbf{X}_{0}\), given in longitude-latitude coordinates, and its output is a Gaussian mixture distribution with \(N_{c}=32\) mixture components modelling \(\Delta\mathbf{X}\mid\mathbf{X}_{0}\), also in degrees of longitude and latitude from \(\mathbf{X}_{0}\). The neural network part of the model thus encodes

\[\mathbf{\theta}(\cdot)=\left\{\alpha_{i}(\cdot),\ \mathbf{\mu}_{i}(\cdot),\ \mathbf{C}_{i}( \cdot)\right\}_{i=1}^{N_{c}} \tag{12}\]

such that

\[\begin{split}\dot{p}(\Delta\mathbf{X}\mid\mathbf{X}_{0})=\sum_{i=1}^{N_{ c}}\alpha_{i}(\mathbf{X}_{0})\det\left(2\pi\mathbf{C}_{i}\left(\mathbf{X}_{0}\right) \right)^{-\frac{1}{2}}\\ \times\exp\left[-\frac{1}{2}\left(\Delta\mathbf{X}-\mathbf{\mu}_{i} \left(\mathbf{X}_{0}\right)\right)^{\mathrm{T}}\mathbf{C}_{i}^{-1}\left(\mathbf{X}_{0} \right)\ \left(\Delta\mathbf{X}-\mathbf{\mu}_{i}\left(\mathbf{X}_{0}\right)\right)\right],\end{split} \tag{13}\]

where \(\mathbf{\mu}_{i}\) and \(\mathbf{C}_{i}\) are the mean vector and covariance matrix of mixture component \(i\). The number of mixture components is a hyperparameter which could be optimised. We chose \(N_{c}=32\) on the basis that 32 component mixtures were found to be sufficiently expressive in trial experiments with MDNs.

The architecture chosen for the neural network is the standard multilayer perceptron, with six hidden (i.e. interior) layers. The first four hidden layers have 256 neurons and the remaining two have 512. The activation function \(\tanh(x)\) is applied to each of the hidden

Figure 3: Regions considered for model comparison in section 4.3.
layers. Thus, the activity of hidden layer \(i\), is

\[\mathbf{h}_{i}=\tanh\left(W_{i}\mathbf{h}_{i-1}+\mathbf{b}_{i}\right) \tag{14}\]

for \(i>1\), and

\[\mathbf{h}_{1}=\tanh\left(W_{1}\mathbf{X}_{0}+\mathbf{b}_{1}\right). \tag{15}\]

Here, \(\mathbf{W}_{i}\in\mathbb{R}^{d_{i}\times d_{i-1}}\) and \(\mathbf{b}_{i}\in\mathbb{R}^{d_{i}}\) are the weight and bias parameters corresponding to the \(i^{\text{th}}\) layer, having \(d_{i}\) neurons. Note that \(\mathbf{w}=\{\mathbf{W}_{i},\mathbf{b}_{i}\}\). The final layer has custom activation functions designed to enforce the natural constraints on the components of \(\mathbf{\theta}\). In particular, the softmax activation function \(a_{\text{sm}}(\mathbf{x})=\exp(\mathbf{x})\,/\,\sum_{i}\exp(x_{i})\) is applied to the neural network outputs which correspond to the mixture component weights, \(\mathbf{\alpha}\), to ensure that these are positive and satisfy the constraint (7). Each covariance matrix \(\mathbf{C}_{i}\) is represented by the components of a lower triangular Cholesky factor \(-\) positivity of the diagonal elements is enforced by taking an exponential. When \(N_{c}=32\) we have \(\dim(\mathbf{\theta})=192\), and the total number of neural network parameters, i.e. weights and biases, is \(\dim(\mathbf{w})=690,880\). We train the model by minimising the negative log likelihood loss function (8) using the Adam algorithm (Kingma & Ba, 2015). We note that the number of widths of hidden layers are further hyperparameters which we have chosen after experimentation with test problems. We do not attempt to find optimal values for these in this work.

As is common in machine learning, we standardise the data before training (LeCun et al., 2012), that is we transform both the input data, \(\{\mathbf{X}_{0i}\}\), and output data, \(\{\Delta\mathbf{X}_{i}\}\), separately, by subtracting the mean of the training data and dividing each component by its standard deviation in the training data, so that each component of the transformed data has zero mean and unit variance. While theoretical justifications for this practice are lacking or unsatisfactory, we found that it did improve noticeably the numerical stability of the optimisation procedure. In any case, the transformation that we apply is invertible, although care must be taken to correctly invert the rescaling of the transition density. For example, if we denote the standardised variables by \(\widehat{\mathbf{X}_{0}}\) and \(\widehat{\Delta\mathbf{X}}\), then the model approximates \(p(\widehat{\Delta\mathbf{X}}\mid\widehat{\mathbf{X}_{0}})\), and we can recover the transition density with the correct units as

\[p(\Delta\mathbf{X}\mid\mathbf{X}_{0})=\frac{p\left(\widehat{\Delta\mathbf{X}}\mid\widehat {\mathbf{X}_{0}}\right)}{\widehat{\text{std}}(\Delta X)\,\widehat{\text{std}}( \Delta Y)}, \tag{16}\]

where \(\widehat{\text{std}}(\,\cdot\,)\) denotes the sample standard deviation among the training data.

Figure 4: Schematic of the MDN model of the single-particle transition density of drifters.
One aspect of neural networks that is particularly relevant to the problem at hand, is that they struggle to represent periodic functions (Liu et al., 2020). Given that we operate in longitude-latitude coordinates, a model of the transition density ought to be periodic in longitude. However, since the neural network model receives the initial position \(\mathbf{X}_{0}\) as simply a vector in \(\mathbb{R}^{2}\), the concept of a spherical domain is not built in to the representation. Indeed, the MDN model produces discontinuities in \(p(\Delta\mathbf{X}\mid\mathbf{X}_{0})\) at the dateline due to model error on either side. To improve continuity at the dateline we employ a crude technique, wherein we replicate the data twice, once shifted by \(360^{\circ}\) longitude west, and once shifted by \(360^{\circ}\) east.

The model is implemented in Python using TensorFlow (TensorFlow Developers, 2021) and TensorFlow Probability (TensorFlow Probability Developers, 2021) and trained using four NVIDIA Tesla V100 16GB GPUs in parallel. Of the 50% of data used to construct the model, 90%, again chosen randomly, was used to inform iterations of the optimisation procedure, and 10% was used for early stopping -- we refer to these portions of the data as the training and test sets, respectively. Training took approximately 90 minutes. The evolution during training of the loss function on training and test sets is shown in figure 5 as a function of epoch. An epoch is the number of iterations taken for all data to be used once in the Adam algorithm. The stopping criterion used for the optimisation, an example of early stopping, was that the test loss had not decreased since 50 epochs previous.

### Model evaluation and comparison

Since the MDN model is probabilistic, its performance should be assessed using skill scores for probabilistic forecasts, as opposed to performance metrics commonly used for deterministic models, such as the mean squared error. As discussed above, minimising the negative log likelihood is equivalent to maximising the log score, since this is exactly the log likelihood. The log score has attractive properties, namely strict propriety (Brocker and Smith, 2007) and locality (Du, 2021). Indeed it is the only smooth local strictly proper scoring rule for continuous variables up to affine transformation (Bernardo, 1979). A scoring rule is strictly proper

Figure 5: Evolution of the training and test loss in the MDN model during optimisation. The loss shown is the mean negative log likelihood per datapoint (i.e. a normalised form of (8)) in terms of the standardised variables \(\mathbf{\widetilde{X}}_{0}\) and \(\widehat{\Delta\mathbf{X}}\).
if its expectation (with respect to data) is maximised uniquely by the correct/perfect model (assuming it exists). A scoring rule is local if it is a function only of the value of the forecast probability distribution evaluated at the observed data, and does not depend for example on other features of the forecast distribution, such as its shape. For validation purposes we can compute the log score on our validation data set. However, while the value of the log score can be easily interpreted in the case of forecasts of discrete/categorical variables, its value in the case of continuous variables is not immediately meaningful, since it refers to probability density, which has dimensions inverse to the area of its support, meaning that the scale of the log score is problem-dependent. On the other hand, the log score can be more easily interpreted when used as a relative score between models -- in particular, the mean difference of log scores between models reflects the average additional probability the first model places on observed outcomes compared to the other model, measured in units of information, nats (or shannons when \(\log_{2}\) is used in the definition of the score). The difference of log scores is invariant under smooth transformations of the forecast variable (Du, 2021); this means, in particular, that differences in log scores are unaffected by a change of units. Thus, in order to evaluate the MDN model we compare it with alternative models. We describe two alternative models, one used extensively in the literature, and one proposed here as a simple but reasonable alternative. We also compare with a simplified version of the MDN model, which features only one mixture component, i.e. for which \(N_{c}=1\). The log score of all models is computed on both training and validation data to assess relative performance.

#### 4.3.1 Transition matrix model

Previous work (Maximenko et al., 2012; van Sebille et al., 2012; Miron et al., 2017, 2021), modelled drifter dynamics with a discrete-time Markov chain using Ulam's method (Ulam, 1960; Froyland, 2001). This requires to discretise space into bins \(\{B_{i}\}\) and estimate the transition matrix

\[P_{ij}=\mathbb{P}(\mathbf{X}_{n+1}\in B_{j}\mid\mathbf{X}_{n}\in B_{i}), \tag{17}\]

which is the discrete analogue of the transition density (10). Indeed the primary difference between a Markov chain model and our Markov process model is that ours is continuous in space. The elements of the transition matrix are usually estimated by the standard approach sketched in figure 1, where we have \(\mathbf{Y}=\mathbf{X}_{n+1}\) and \(f(\mathbf{Y})=1_{B_{j}}(\mathbf{X}_{n+1})\) -- this corresponds to the maximum-likelihood estimate for each \(P_{ij}\) and, hence, maximises the log score on the training dataset. Note that the transition matrix can be used to construct a corresponding transition density which is piecewise constant on gridcells in \(\mathbf{X}_{n}\) and \(\mathbf{X}_{n+1}\) via

\[p(\Delta\mathbf{X}\mid\mathbf{X}_{0})=\frac{P_{ij}}{A(B_{j})},\quad\text{ when }\mathbf{X}_{0}\in B_{i},\ \mathbf{X}_{0}+\Delta\mathbf{X}\in B_{j}, \tag{18}\]

where \(A(B_{j})\) is the area2 of \(B_{j}\). This is important for allowing comparison with models which are continuous in space.

Footnote 2: For consistency with the transition density as given by the MDN model, these areas must be calculated in terms of the same variables, i.e. degrees longitude by degrees latitude.

An advantage of Markov chain models is that analysis of their long time behaviour is straightforward -- the left and right eigenvectors of the transition matrix can be studied to identify almost-invariant sets, as in Miron et al. (2017). This has been called the eigenvector method (Froyland et al., 2014). The extension of this analysis to the continuous-space setting
using our model, which we leave for future work, requires the calculation of eigenfunctions of the relevant Perron-Frobenius operator \(\mathcal{P}\), which acts on probability density functions to evolve them forward in time, such that

\[p(\mathbf{X}_{n+1}) =\mathcal{P}(p(\mathbf{X}_{n})) \tag{19}\] \[:=\int_{\Omega}p(\mathbf{X}_{n})\,p(\mathbf{X}_{n+1}\mid\mathbf{X}_{n})\, \mathrm{d}\mathbf{X}_{n}. \tag{20}\]

Alternatively, it is worth noting that the MDN model can be used to construct a transition matrix, by numerical integration of the transition density, that is by computing numerically

\[P_{ij}=\int_{B_{i}}\int_{B_{j}}p(\mathbf{X}_{n+1}\mid\mathbf{X}_{n})\,\mathrm{d}\mathbf{X}_{ n+1}\,\mathrm{d}\mathbf{X}_{n}. \tag{21}\]

In figure 6 we show the log transition density \(\log p(\Delta\mathbf{X}\mid\mathbf{X}_{0})\) derived from the transition matrix model via (18) for two different initial positions \(\mathbf{X}_{0}\). The first is located within the core of the Gulf Stream at \(34.85^{\circ}\) N, \(74.50^{\circ}\) W, and the second is just outside the Gulf stream at \(33.67^{\circ}\) N, \(72.55^{\circ}\) W. Notice that in each case the support of the density is the set of grid cells to which transitions were observed in the training data. In other words, transitions to other grid cells have probability zero under the model. We return to this point in section 4.3.3.

#### 4.3.2 Gaussian transitions with gridded parameters (GTGP)

A simple model for the transition density (11) is that, given initial positions \(\mathbf{X}_{0}\), transitions are conditionally Gaussian with conditional mean and covariance given by functions of \(\mathbf{X}_{0}\) which are piecewise constant on grid cells, i.e.

\[\Delta\mathbf{X}\mid\mathbf{X}_{0}\sim\mathcal{N}\left(\mathbf{\mu}\left(\mathbf{X}_{0}\right),\,\mathbf{C}\left(\mathbf{X}_{0}\right)\right), \tag{22}\]

with \(\mathbf{\mu}\left(\mathbf{X}_{0}\right)\) and \(\mathbf{C}\left(\mathbf{X}_{0}\right)\) piecewise constant in \(\mathbf{X}_{0}\). The parameters \(\mathbf{\mu}\) and \(\mathbf{C}\) are estimated by sorting the observations into bins and computing the sample mean and sample covariance for each bin. The sample mean is the maximum likelihood estimate of \(\mathbf{\mu}\), while the sample

Figure 6: Maps of the log transition probability density function, \(\log p(\Delta\mathbf{X}\mid\mathbf{X}_{0})\), for initial positions, \(\mathbf{X}_{0}\), (a) in the Gulf Stream (\(34.85^{\circ}\) N, \(74.50^{\circ}\) W), and (b) adjacent to the Gulf Stream (\(33.67^{\circ}\) N, \(72.55^{\circ}\) W), derived from the transition matrix model with \(\tau=4\) days via (18). Yellow dots indicate \(\mathbf{X}_{0}\).
covariance differs from the maximum likelihood estimate of \(\mathbf{C}\) only by a factor of \(\frac{N-1}{N}\approx 1\), where \(N\) is the number of training data in the given bin. Hence, the parameter estimates used are very close to those which maximise the log score on training data.

Figure 7 shows the mean of displacements from a GTGP model with a regular \(1^{\circ}\times 1^{\circ}\) longitude-latitude grid and \(\tau=4\) days, as a function of initial position.

#### 4.3.3 Model scores

We compute skill scores for the full training and validation datasets with global coverage, and for regions \(A\), \(B\), and \(C\). For both the transition matrix and GTGP models it is necessary to choose a spatial discretisation; herein we consider only square latitude-longitude grids,

Figure 7: Mean of displacements from the GTGP model, with \(\tau=4\) days, as a function of initial position.
so that the only parameter to be chosen is the grid cell side length. This choice affects their performance. If a relatively high resolution discretisation is used, the models attain relatively high scores in training, but generalise poorly, as reflected in poor scores on validation data. In the case of the GTGP model, an issue arises when validation data falls in grid cells not visited by drifters in the training set, since sample means and covariances cannot be estimated in bins where data is absent. As a simple solution, we set the value of \(\mathbf{\mu}\) and \(\mathbf{C}\) on unvisited grid cells equal to a global (or regional) estimate. A flaw of the transition matrix model, with the transition matrix estimated as discussed above, is that validation data can have zero probability under the model, and hence achieve a log score of minus infinity. This situation is avoided by taking sufficiently large grid cells, but this leads to exceptionally low scores. On the other hand, a validation score can be computed with smaller grid cells if one is prepared to simply discard validation data which have zero probability under the model. This seems overly generous, as the transition matrix model will be scored increasingly highly as the grid cell size is reduced to zero and an increasing number of the validation data are neglected. As a compromise, we fix the grid cell size for the transition matrix model to \(1^{\circ}\times 1^{\circ}\), the resolution used in some previous studies (van Sebille et al., 2012) where the transition matrix model was used, and discard validation data with zero probability -- the proportion of validation data discarded was 7% globally, and 2%, 9% and 4% in regions \(A\), \(B\) and \(C\), respectively. For the GTGP model the grid cell size was optimised to maximise validation scores using grid search cross validation: a common procedure which amounts to trying a range of values of a model hyperparameter (in this case the grid cell size) and choosing the value which optimises the validation score. The optimal grid cell size found ranged from \(1.1^{\circ}\) in region A to \(5^{\circ}\) globally. The scores are presented in table 1. In all regions the MDN models outperform the alternatives, with the 32-component model achieving slightly higher scores than the single-component model. Note that the scores reported happen to be negative -- this is not by convention, but instead reflects that log probability densities are often negative. A higher score is a better score.

### Results

Once trained, the model can be used in at least two ways: (i) to derive estimates of single-particle displacement statistics, and (ii) to simulate drifter trajectories. However, we first ex

\begin{table}
\begin{tabular}{l l|c c c c}  & & TM & GTGP & MDN1 & MDN32 \\ \hline Global: & Training & \(-1.61\) & \(-1.20\) & \(-1.07\) & \(-1.02\) \\  & Validation & \(-1.74\) & \(-1.30\) & \(-1.14\) & \(-1.11\) \\ \(A\): & Training & \(-1.78\) & \(-1.25\) & \(-1.35\) & \(-1.30\) \\  & Validation & \(-1.89\) & \(-1.42\) & \(-1.40\) & \(-1.35\) \\ \(B\): & Training & \(-1.95\) & \(-1.90\) & \(-1.91\) & \(-1.85\) \\  & Validation & \(-2.16\) & \(-2.01\) & \(-2.01\) & \(-1.96\) \\ \(C\): & Training & \(-1.93\) & \(-1.59\) & \(-1.63\) & \(-1.56\) \\  & Validation & \(-2.00\) & \(-1.66\) & \(-1.61\) & \(-1.57\) \\ \end{tabular}
\end{table}
Table 1: Training and validation scores for the transition matrix and GTGP models, as well as the single-component MDN and full 32-component MDN models, in each of the regions considered (see the map in figure 3). The scores are the mean log score (i.e. mean log likelihood) per datapoint calculated in terms of the variables \(\mathbf{X}_{0}\) and \(\Delta\mathbf{X}\) in their original degrees longitude/latitude units.
amine the transition density directly. In figure 8 we show the log transition density \(\log p(\Delta\mathbf{X}\mid\mathbf{X}_{0})\) for two different initial positions \(\mathbf{X}_{0}\). In the first case, where \(\mathbf{X}_{0}\) is located within the core of the Gulf Stream at \(34.85^{\circ}\) N, \(74.50^{\circ}\) W, the transition density is strongly non-Gaussian, with contours extending roughly to the south and northeast, showing the influence of the Gulf Stream on drifters. In the second case, where \(\mathbf{X}_{0}\) is just outside the Gulf stream at \(33.67^{\circ}\) N, \(72.55^{\circ}\) W, the transition density is closer to Gaussian.

In order to quantify the extent to which the transition density deviates from being Gaussian, and how this varies from one region of the ocean to another, we computed the Kullback-Leibler (KL) divergence3 of the single-component MDN model, which is Gaussian, from the full 32-component model as a function of initial position. The result is shown in figure 9. Note that, since a closed-form expression for the KL divergence between two Gaussian mixtures is not known (Cui & Datcu, 2015), we provide simple Monte Carlo estimates based on 5000 samples at each of the vertices of a \(1^{\circ}\times 1^{\circ}\) grid. Where the KL divergence is zero, the two models agree exactly, indicating that displacements are Gaussian. The larger the KL divergence is, the greater the disagreement between the models, and the further from Gaussian the full model is. As a point of reference for interpreting the magnitude of the KL divergence, note that the if \(Z_{0}\sim\mathcal{N}(m_{0},\,1)\) and \(Z_{1}\sim\mathcal{N}(m_{1},\,1)\), then, writing their pdfs as \(p_{0}\) and \(p_{1}\), \(D_{\text{KL}}(p_{1}\parallel p_{0})=(m_{1}-m_{0})^{2}\). Non-Gaussianity of displacements is likely due primarily to inhomogeneity of ocean velocities -- drifters can explore a range of flow statistics as they move, and the convolved effects of these are reflected in observed displacements. An alternative explanation is that the underlying velocity field is non-Gaussian -- evidence of non-Gaussian velocities in the North Atlantic has been presented by Bracco et al. (2000) and LaCasce (2005) on the basis of observations from both subsurface current meters and subsurface floats.

Footnote 3: The KL divergence of \(p\) from \(q\), also known as the relative entropy, defined \(D_{\text{KL}}(q\parallel p)=\int q(x)\,\log\frac{q(x)}{p(x)}\,\text{dx}\), is a measure of the divergence of a probability density \(p\) from a reference probability density \(q\) – often interpreted as the amount of information lost when \(p\) is used to approximate \(q\).

As can be seen in figure 8, the model assigns nonzero probability to drifter displacements intersecting land. This is unavoidable given that the support of the assumed parametric form, that of a Gaussian mixture, extends to infinity; moreover, this may not be entirely spurious, given that some drifters do run aground. In 2012 Lumpkin et al. (2012) reevaluated drifter data to study the causes of drifter deaths. They concluded that approximately 27% of drifter

Figure 8: Maps of the log transition density, \(\log p(\Delta\mathbf{X}\mid\mathbf{X}_{0})\), for initial positions, \(\mathbf{X}_{0}\), (a) in the Gulf Stream (\(34.85^{\circ}\) N, \(74.50^{\circ}\) W), and (b) adjacent to the Gulf Stream (\(33.67^{\circ}\) N, \(72.55^{\circ}\) W), derived from the MDN model with \(\tau=4\) days. Yellow dots indicate \(\mathbf{X}_{0}\).
deaths were due to running aground, with a further 10% being picked up by humans, and the remainder failing due to internal faults. Outside of coastal regions this issue is unlikely to have a strong effect on the estimates of displacement statistics considered in section 4.4.1. The implications for dritter simulations are discussed further in section 4.4.2.

#### 4.4.1 Displacement statistics

In this section we present maps of single-particle statistics derived from the model. As a first example, we show the mean of displacements over the 4-day time increment of our model. We further provide global estimates of lateral diffusivity.

Figure 10 shows the mean of dritter displacements as a function of initial position. While the output of the model is in longitude-latitude coordinates \((\lambda,\phi)\), we apply a simple conversion to kilometres based on a local tangent-plane approximation

\[\Delta X =R\,\Delta\phi \tag{23a}\] \[\Delta Y =R\,\Delta\lambda\,\cos\phi_{0}, \tag{23b}\]

where \(R\) is the radius of the Earth at the equator. The imprint of several features of the surface dynamics, such as the western boundary currents and equatorial (counter) currents, is clear.

For the sake of comparison with previous work, we consider the estimation of lateral diffusivity from our model, though we emphasise that by modelling the full transition density, we provide a more accurate description of Lagrangian statistics than can be captured by the familiar advective-diffusive model of dispersion put forward by Davis (1987, 1991). The estimation of ocean diffusivity by various methods has been the subject of numerous papers (Oh et al. 2000, Zhurbas & Oh 2003, 2004, Klocker et al. 2012, Abernathey & Marshall 2013, Klocker & Abernathey 2014, Ying et al. 2019). The estimation of diffusivity from dritter displacements is straightforward only when there exists a suitable sampling time, which is larger than the time for dritter velocities to decorrelate, i.e. for dritter motion to become diffusive, and such

Figure 9: Kullback–Leibler divergence of the single-component MDN model from the full 32-component MDN model, as a function of initial position. Larger values indicate stronger deviations from Gaussianity in displacements.
Figure 10: Mean of displacements from the MDN model, with \(\tau=4\) days, as a function of initial position.
that the scale of drifter displacements over that time scale is small relative to spatial variations in the diffusivity. In this case a simple estimate of the lateral diffusivity tensor \(\mathbf{K}(\mathbf{x})\) is

\[\mathbf{K}(\mathbf{x})=\frac{1}{2\,\tau}\,\text{Cov}(\Delta\mathbf{X}\mid\mathbf{X}_{0}=\bm {x}), \tag{24}\]

where \(\tau\) is the suitably chosen time scale, and the conditional covariance is estimated by either one of the approaches sketched in figure 1. Unfortunately, such a time scale may not exist in the ocean, and, if it does exist, it likely varies in space, making its determination difficult. This challenge has been borne out in previous studies (LaCasce et al. 2014, Zhurbas et al. 2014). Oh et al. (2000) proposed a method to circumvent the issues created by inhomogeneity. They proposed to isolate the cross-flow component of the displacement covariance, identified by the minor principal component (the smaller eigenvalue of displacement covariance), and use this to provide a scalar estimate of diffusivity, since the cross-flow component is less affected by shear in the mean flow. In figure 11(a) we provide a similar estimate, derived from the MDN model with \(\tau=4\) days,

\[K(\mathbf{x})=\frac{1}{2\,\tau}\,\lambda_{2}(\mathbf{x}), \tag{25}\]

where \(\lambda_{2}(\mathbf{x})\) is the smallest eigenvalue of

\[\text{Cov}(\Delta\mathbf{X}\mid\mathbf{X}_{0}=\mathbf{x})=\sum_{i}\alpha_{i}\mathbf{C}_{i }+\sum_{i}\left[\alpha_{i}\left(\mathbf{\mu}_{i}-\sum_{i}\alpha_{i}\mathbf{\mu}_{i} \right)\left(\mathbf{\mu}_{i}-\sum_{i}\alpha_{i}\mathbf{\mu}_{i}\right)^{T}\right]. \tag{26}\]

The result agrees very well with estimates provided by Zhurbas & Oh (2004) for the Atlantic and Pacific oceans. Figure 11(b) shows the difference of estimates of the form (25) with \(\tau=14\) days and \(\tau=4\) days, respectively. In many areas, the diffusivity estimates are slightly amplified by taking a larger time lag \(\tau\), with greater differences visible in some particularly energetic regions; however, the effect is indeed much weaker than that observed with analogous along-flow diffusivity estimates derived from the largest eigenvalue of the displacement covariance matrix.

Leaving aside the challenges of estimating diffusivity from displacements, which are common to all methods, we highlight as this point some advantages to our approach. Using the MDN model, trained with maximum likelihood and effectively regularised by the use of early-stopping, removes the difficulty of tuning the resolution of bins. Instead, the effective resolution of our mean displacement and diffusivity estimates is set automatically by the resolution of the data, and is free to vary optimally in space. This allows us to produce at once global estimates, which resolve well-sampled flow features very well and are forgiving in regions where data is relatively sparse, with the exception of very high-latitude regions, where there is simply no data to constrain the model.

#### 4.4.2 Drifter simulations

In this section we demonstrate the simulation of drifter trajectories using the MDN model as the basis for a discrete time Markov process model. In a discrete-time setting, assuming Markovianity means assuming that \(p(\mathbf{X}_{n+1}\mid\mathbf{X}_{n},\,\mathbf{X}_{n-1},\,\cdots)=p(\mathbf{X}_{n+1}\mid\mathbf{X} _{n})\). In this case, sampling trajectories amounts to repeatedly sampling displacements in sequence according to the transition density, since, given the current position, displacements are statistically independent of previous positions.
Figure 11: Global estimate of lateral diffusivity derived from the MDN model of the transition density.
A complication of simulating drifters in this way is that, for reasons discussed above, drifters can hit land. In this work we do not attempt to model the beaching of drifters, since it is not clear that the Global Drifter Program dataset contains sufficient reliable information \(-\) in particular, it remains a challenge to determine whether drifters have run aground or not (Lumpkin et al., 2012). To exclude the possibility of running aground in our drifter simulations we implement a simple rejection sampling scheme, wherein displacements sampled from the transition density which would bring a drifter on land are rejected, and a new displacement is sampled until a displacement which keeps the drifter in the ocean is drawn. This amounts to sampling according to the conditional density \(p(\Delta\mathbf{X}\mid\mathbf{X}_{0},\mathbf{X}_{0}+\Delta\mathbf{X}\not\in\text{land})\), and is equivalent to the standard practice when using transition matrix models of restricting the domain considered to the ocean and normalising probability estimates correspondingly. To determine whether a proposed new position is on land, we check intersection with a 110m-resolution land mask.

We simulated the evolution of a set of drifters initialised on the vertices of a \(2^{\circ}\times 2^{\circ}\) grid for a period of 10 years. Note that the evolution of each drifter is simulated independently. This means that multi-particle statistics that would characterise the joint evolution of drifters released simultaneously in the ocean are not represented and, in particular, that the current model is not appropriate for simulating the release of a cloud of tracer particles on short time scales; however, it can be expected to represent the behaviour of drifters or buoyant tracers over large spatial and temporal scales. Similar experiments, carried out by Maximenko et al. (2012) and van Sebille et al. (2012) using transition matrix models trained on Global Drifter Program data, studied the clustering of simulated drifters due to near-surface convergence and the formation of so-called garbage patches, including the North Pacific Garbage Patch (Moore et al., 2001) and others corresponding to the other subtropical ocean gyres. The simulations of van Sebille et al. (2012) showed a further cluster in the Barents Sea which formed only after several decades.

The results of our model simulation are largely in agreement with these previous studies. The distribution of the simulated drifters is shown in figure 12 at the beginning of the simulation and after one, three, and ten years of evolution under the MDN model. After one year the drifters have become relatively sparse in equatorial regions. After three years clusters in the subtropical gyres have begun to appear, and after ten years, these are very well defined. Smaller clusters are also seen to appear, notably in the North Sea and in the seas south of Papua, as well as in some high latitude regions including along the west coast of Greenland and off Antarctica around \(100-130^{\circ}\) E.. Validating these clusters, that is, assessing whether marine debris is likely to accumulate in these areas, is difficult, because in situ observations remain sparse (Ryan et al., 2009). It may be that the dynamics in these regions, which are poorly sampled by GDP drifters, are simply underresolved by the MDN model, leading to spurious convergence zones. We note, for example, in figure 10 that mean displacements do not appear to represent the detail of known currents in the southern portion of the North Sea, which is not visited by drifters in the GDP data (see figure 2). In general, as is true for any data-driven model, caution should be exercised when interpreting model outputs in regions where data is lacking.

## 5 Conclusions

This work demonstrates the use of conditional density estimation, and, in particular, stochastic neural networks, in a fluid dynamical problem, namely that of diagnosing single-particle
statistics from trajectory data. We show how such probabilistic models are useful both as emulators, and as an indirect means of estimating conditional statistics. By operating in the framework of probabilistic modelling we are able to appeal to the extensive literature on statistical inference, probabilistic forecasting, model comparison and validation, and thereby avoid ad hoc choices of loss functions and performance metrics. Our model is compared, using a probabilistic scoring rule, to alternative models, including a Markov chain model used extensively in the literature, and is shown to outperform these, both globally and in three specific regions.

By modelling the single-particle transition density of surface drifters, we gain estimates of a range of conditional statistics simultaneously, which capture the occurrence of strongly non-Gaussian statistics in some areas of the ocean. We provide global maps of mean displacement and lateral diffusivity, but emphasise that these examples provide only a limited summary of the information contained in the transition density; further statistics, including higher moments of displacements can readily be computed from our model. Interpreted as the basis for a discrete-time Markov process, our model is also used to simulate the evolution of a set of drifters seeded globally on a uniform grid, and shows the emergence of clusters of drifters in the subtropical gyres, in agreement with previous work on the formation of garbage patches.

The approach espoused in this work is equally applicable to other problems in fluid dynamics and oceanography. One example is the estimation of structure functions from either Eulerian or Lagrangian velocity data. Another is the estimation of multi-particle statistics, such as relative dispersion, via modelling of multi-particle transition densities. Yet another is the learning of stochastic parameterisations in climate/atmosphere/ocean models. Guillau

Figure 12: Histograms of simulated drifters initially and after one, three, and ten years of evolution under the MDN model, respectively.
min & Zanna (2021) made progress on the parameterisation of subgrid momentum forcing in an ocean model with a single-component MDN model, but the approach is applicable more broadly, e.g. to the parameterisation of subgrid transport.

In this work we have largely neglected the need to quantify uncertainty in model parameters and to incorporate prior knowledge in our modelling. These needs would be met by a Bayesian approach, where, instead of estimating parameters by maximum likelihood, we apply Bayesian inference to obtain posterior distributions on parameters, which account for prior knowledge. Indeed, all of the results presented herein would benefit from uncertainty quantification. In the case of conditional statistics, a Bayesian approach would, e.g., allow to identify where there is not enough data to inform reliable estimates of lateral diffusivity; and in general, incorporating prior knowledge may help to regularise our model of the transition density, so that, in the case of drifter simulations, spurious convergence zones can be avoided. The application of Bayesian inference to MDNs remains challenging, but we consider this an important future direction.

**Acknowledgements.** I am grateful to Jacques Vanneste, James Maddison and Aretha Teckentrup for their advice, input and overall support of this work. I also thank Dhruv Balwada for helpful discussions. Thanks are also due to the reviewers for their suggestions which have improved the manuscript.

**Funding**. The author was supported by the MAC-MIGS Centre for Doctoral Training under EPSRC grant EP/S023291/1. This work used the Cirrus UK National Tier-2 HPC Service at EPCC (www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1).

**Declaration of interests.** The author reports no conflict of interest.

**Data availability statement.** The code required to reproduce the results herein is available at doi.org/10.5281/zenodo.7737161, along with the trained MDN model and a Jupyter Notebook which demonstrates its use. The processed GDP data used and drifter simulation data are available at doi.org/10.7488/ds/3821.
