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

HAL
open science

# Coherent motions to predict Lagrangian trajectories

Ali Rahimi Khojasteh, Dominique Heitz

► To cite this version:

Ali Rahimi Khojasteh, Dominique Heitz. Coherent motions to predict Lagrangian trajectories. 2025. [https://hal.inrae.fr/hal-05232948v1](https://hal.inrae.fr/hal-05232948v1)

HAL Id: hal-05232948

[https://hal.inrae.fr/hal-05232948v1](https://hal.inrae.fr/hal-05232948v1)

Preprint submitted on 1 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Banner appropriate to article type will appear here in typeset article

1

# Coherent motions to predict Lagrangian trajectories

2

Ali Rahimi Khojasteh<sup>1,2†</sup>, and Dominique Heitz<sup>1</sup>

3

<sup>1</sup>INRAE, OPAALE, 17 avenue de Cucillé, 35044, Rennes, France

4

<sup>2</sup>Laboratory for Aero and Hydrodynamics, Delft University of Technology, 2628 CD Delft,

5

The Netherlands

6

(Received xx; revised xx; accepted xx)

7

Accurate prediction of Lagrangian trajectories in turbulent flow remains challenging due to limited temporal information in transport functions. This paper shows that even with sparse temporal observations, there might be enough information from surrounding coherent motions, sharing the same dynamics, to provide highly probable trajectories. The proposed *coherent predictor* is derived from the concept of Lagrangian coherent structures (LCSs), which are advective transport barriers that govern the cohesive motion of neighbouring particles. Coherent trajectories are quantified using a local segmentation with the finite-time Lyapunov exponents (FTLE). The coherent predictor incorporates information from the particle's position history and neighbouring coherent velocity and acceleration into a novel generic energy function to predict its trajectory. We assess our proposed approach using both three-dimensional (3D) synthetic and experimental data of the wake behind a smooth cylinder and two-dimensional (2D) homogeneous isotropic turbulent (HIT) flow. The coherent predictor is deemed generic due to its consistent behaviour regardless of flow dimensions, Reynolds number, and flow topology. Our results demonstrate that the optimal parameters of the proposed energy function can be modelled based on measurement uncertainties, resulting in enhanced prediction accuracy and reduced uncertainty compared to current methods. We reveal direct signatures of flow topology, including the cylinder leading edge boundary layer, sideward shear layers, and vortex formation structures, on the prediction error map. These topologies, which are fundamental structures in fluid dynamics, are marked by high Lagrangian gradients and 3D directional motions. These findings on coherent predictions hold great potential for various Lagrangian analyses in turbulence.

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

**Key words:** Authors should not enter keywords on the manuscript, as these must be chosen by the author during the online submission process and will then be added during the typesetting process (see Keyword PDF for the full list). Other classifications will be added at the same time.

32

**MSC Codes** (Optional) Please enter your MSC Codes here

† Email address for correspondence: a.r.khojasteh@tudelft.nl
33

# 1. Introduction

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

The transport of tracer particles, which forms Lagrangian trajectories, is a fundamental problem in fluid dynamics, and *Lagrangian particle tracking* (LPT) has become a crucial tool for quantifying the motion of these fluid particles. Applications of Lagrangian trajectories span various fields in physics and engineering, including the mixing process and diffusion properties of tracer particles (Viggiano et al. 2021), vortex shedding topology analysis (Gold et al. 2023), turbulent bubbly jet (Kim et al. 2022). All these applications require trajectory reconstruction that precisely represents Lagrangian motion before undertaking any further analysis. For decades, LPT studies have focused on reconstructing tracer particle trajectories by treating them as isolated signals, separate from their surroundings (Schröder & Schanz 2023). While the use of predictive models has expanded the spatial and temporal capabilities of reconstructing dense and complex Lagrangian trajectories (Schröder & Schanz 2023; Schanz et al. 2016), the remained challenging problem is to precisely predict motions of tracer particles as representatives of Lagrangian dynamics in turbulent and chaotic flows. Trajectory prediction refers to a function that receives tracer positions and estimates the positions of the next time step to reconstruct the Lagrangian trajectory. These techniques rely solely on mathematical models without incorporating physics-based constraints and keep treating particles as isolated signals. The main reason why Lagrangian trajectories are subject to error and uncertain reconstructions is the temporal resolution in experiments. Khojasteh et al. (2021) showed the temporal resolution could be over an order of magnitude larger than the Kolmogorov time scale due to experimental limitations. Consequently, we miss numerous small-scale motions before the next observation (i.e., camera image) is recorded, resulting in a loss of Lagrangian information between two observations. Thus, depending on position history alone is insufficient to fully recover the Lagrangian trajectory with certainty. Lagrangian prediction extends beyond LPT experiments where explicit time-stepping methods, such as the Euler and Runge Kutta schemes, are employed in ocean analyses to determine the transport of tracer particles based on Eulerian simulations (Sebille et al. 2018; Valdivieso Da Costa & Blanke 2004). Due to the high computational cost and complexities of performing *direct numerical simulation* (DNS), researchers often employ coarse models with sparse temporal resolutions such as Oceanic General Circulation Model for the Earth Simulator (OFES Masumoto et al. 2004). As Qin et al. (2014) demonstrated, this eventually can lead to a significant source of error in the temporal interpolation of simulated velocity fields. This is due to a lack of temporal information in the transport function, highlighting a common issue in Lagrangian studies across various fields. Yet, even with sparse temporal observations, there might be enough information from surroundings to provide highly probable trajectories. If successful particle position predictions are possible in turbulent flows, they can be obtained thanks to suitable coherent neighbours. This is the thrust of the paper.

71

72

73

74

75

76

77

78

79

80

81

Knowing that trajectories are characterised by *Lagrangian coherent structures* (LCS) acting as barriers between regions in which tracer particles have different kinematics, we propose a novel approach considering groups of coherent particles, rather than an isolated approach, to estimate Lagrangian particle positions. Numerous techniques have been developed to diagnose LCS barriers and understand fluid transport where particles are generally advected by the flow, from trajectory-based coherent segmentation to vector displacement-based techniques (Haller 2023; Mowlavi et al. 2022; Hadjighasem et al. 2017; Green et al. 2007). The concept of coherency in Lagrangian trajectories refers to regions with similar dynamics while the flow evolves. Haller & Yuan (2000) showed how to characterise these LCSs based on the deformation gradient tensor from noisy velocity fields obtained from experiments. This finding led to more investigations into applying nonlinear dynamical
system theory to fluid flows. Later on, Shadden et al. (2005) studied how passive finite-time Lyapunov exponents (FTLE) ridges are the same as LCS material lines and boundaries of coherent motions. FTLE is a quantitative metric that determines how to extract coherent structures in the Lagrangian framework. These coherent structures are objective and have the most robust representations of trajectories that act as skeletons of flow (MacMillan & Richter 2021). Objective means the material responses of their surface (or line) transport barriers are independent of the observer (Haller 2023, 2015). The further theory was introduced by Haller et al. (2020) to objectively quantify the transport of active vector fields such as momentum and vorticity. We propose a robust prediction of Lagrangian particle positions by leveraging LCS boundaries that control the transport of tracer particles. The use of LCS enables the identification of coherent and non-coherent neighbour trajectories and local geometric separatrices, which can be used to estimate the dynamics of trajectories in their coherent surroundings. Specifically, we utilise FTLE to quantify the boundaries between separated regions, as this method has been previously demonstrated in experimental studies (see, e.g., Eisma et al. 2021) and has shown to be effective in identifying coherent structures (Balasuriya et al. 2020; Hadjighasem et al. 2017).

We adopt a local segmentation approach instead of computing the global FTLE field to classify neighbours as coherent or non-coherent. The local FTLE approach offers a more targeted analysis of the particle's immediate surroundings, reduces computational and three-dimensional (3D) transport barrier quantification complexities, and provides a precise understanding of the coherent motions in the vicinity of particles. We introduce primary and secondary coherent neighbours in turbulent flows. Primary coherent neighbours travel along a similar path as the target particle within the same time interval. In contrast, secondary coherent neighbours exhibit a phase delay and are situated ahead in their history relative to the target particle. This secondary information, representing future knowledge for the target particle, provides valuable prior information for predicting Lagrangian trajectories. We, therefore, define a generic energy function called *coherent predictor* that considers the history of trajectories and local coherent motions for prediction with the least biased and low uncertainty level compared to the state-of-the-art. We show that the coherent predictor can be modelled as a function of measurement uncertainty regardless of the flow case. That means the proposed approach remains generic to the best of available test cases in the present study. By "generic," we mean that the behaviour of the coherent predictor remains independent of flow dimensions (2D or 3D), Reynolds number, and flow topology. The coherent predictor offers valuable insights into Lagrangian physics by examining groups of coherent motions, leading to a deeper understanding of how effectively they share the same dynamics even in low temporal resolutions. Precise prediction also allows for the reconstruction of longer trajectories instead of split tracklets throughout the entire observation event, which provides valuable temporal information for Lagrangian statistics analysis.

This paper is structured as follows. In § 2, we introduce physics-based coherent terms as an additional constraint to the mathematical prediction model in the form of an energy function, followed by local segmentation using the FTLE metric. Then § 3 presents the numerical and experimental flow configurations employed in this study, including DNS for two-dimensional (2D) homogeneous isotropic turbulent (HIT) flow at a Reynolds number of 3000, the DNS of the cylinder wake flow at a Reynolds number of 300, and both DNS and experimental investigations for the 3D wake flow behind a smooth cylinder at a Reynolds number of 3900. We discuss the generic energy function and its minimisation procedure in § 4.1 and present the sensitivity and confidence analyses of Lagrangian position estimation based on coherent motions in § 4.2 and § 4.3, respectively. We demonstrate the effectiveness of our proposed approach through a real particle tracking experiment in § 4.4, where we find a minimum
deviation from the optimised positions. Finally, we provide a summary of our conclusions in § 5.

132 § 5.

## 2. Lagrangian position estimation

In the Lagrangian framework, we focus on the study of fluid motion from the perspective of individual fluid particles as they move through space and time. Adopting this framework allows for the development of prediction functions that can capture the complex dynamics of particles within a fluid. The most simple prediction approach is using a polynomial function, suggested by Schanz et al. (2013), resulting in reasonable 3D trajectory reconstructions in simple flows (Schröder et al. 2015b; Schanz et al. 2014, 2013). Significant misprediction occurs in case of flow associated with complexities such as high turbulence level and high Reynolds number (Tan et al. 2020). In such conditions, even by increasing the polynomial order from three to ten, misprediction stays remained (Tan et al. 2020). As an alternative, Schröder et al. (2015b) introduced prediction by optimal temporal filtering of particle positions, such as the Wiener filter in Lagrangian tracking experiments. Since then, this concept has become widely adopted in the time-resolved particle tracking velocimetry (4D-PTV) studies (Tan et al. 2020; Schröder et al. 2015a). The Wiener filter showed robust behaviour in prediction with complex flows such as inside the turbulent boundary layer (Knopp et al. 2021). However, it still suffers in high-motion gradients because the temporal filter cancels gradients unless an additional constraint is defined. This implies that the prediction function sticks out from a lack of information to estimate the dynamics of particles. Unlike tracer particles, Lagrangian fibres provide broad details of their dynamics, including orientation and rotation rates (Alipour et al. 2021). This additional information on Lagrangian dynamics can facilitate the prediction process. For individual tracer particles, our knowledge is limited to their histories. Despite implementing filtering and smoothing schemes such as Shake-the-Box (STB) using Wiener filter (Schröder et al. 2015a), our prediction is limited by the history of the target particle, ignoring every Lagrangian trajectory is spatially and temporally coherent with a specific group of other tracer particles following the same behaviour. It is not quantitatively studied after what turbulence intensity level or velocity gradient the predictor function fails. This isolated prediction approach causes wrong trajectory reconstructions that are artificially created and fails to represent true Lagrangian dynamics. Therefore, any further analysis would lead to uncertain and inaccurate results, particularly in complex flow regions such as wake flows, impinging jets, and areas near solid boundaries (Khojasteh et al. 2021; Yang & Heitz 2021; Schanz et al. 2020).

The prediction function comprises a set of coefficients and independent variables. The coefficients are employed to predict the outcome of a dependent variable, which can be written in the form of an energy function. Therefore, the objective is to minimise the energy function to reach the optimal trajectory prediction at the next time step. Essentially, position histories of trajectories are the main input ingredients of the prediction process, as shown in figure 1.a.b. The predictor subsequently fits a smooth curve over the noisy history {$t_1, \dots, t_n$} of the target particle and estimates its possible position at time step $t_{n+1}$ (see figure 1.c). At this stage, we categorise recently available and the proposed prediction functions based on the input information they require into position-based and Lagrangian coherent energy functions, discussing the former in § 2.1 and the latter in § 2.2. The Lagrangian coherent energy function seeks a solution in which the estimated position satisfies neighbouring coherent motions (see figure 1.d). By incorporating Lagrangian physics into the prediction function, we demonstrate that even a weak signal resulting from the motion of just one coherent neighbour can enhance particle prediction accuracy in complex flow regions, thus facilitating a better understanding of the complex interactions and transport phenomena in fluid systems.
<figure><img src="image_2.png" /><figcaption>Figure 1: Particle position prediction scenario from $t_n$ to $t_{n+1}$. (a) Known particle positions from history starting from $t_{n-4}$ up to $t_n$. (b) Trajectory (dark blue line) obtained from filtered curve fitting of known particle positions. (c) Prediction based on extrapolating the fitted trajectory (golden dashed line) from $t_n$ to $t_{n+1}$. (d) Modified prediction (red dashed line) using velocity and acceleration information of coherent neighbours of the target particle at $t_n$.</figcaption></figure>

## 2.1. Position-based energy function

Two main position-based methods have been reported in the particle tracking algorithms: polynomial and Wiener filter predictors (Yang & Heitz 2021; Tan et al. 2020; Schanz et al. 2016). Both techniques rely solely on the history of the target particle. The polynomial predictor aims to minimise the least mean square of the history to find the optimal polynomial coefficients and then extrapolate the function with the same coefficients to the next time step. In the Wiener filter approach, we design and adjust filter parameters based on history and then shift the designed filter solution to the next time step. Both the Wiener filter and polynomial predictor functions require two parameters to set: the polynomial (or Wiener filter) order and the length of history to perform the Linear regression analysis. The regression window length (history) should be determined by the desired integral time. Then, we should select the order of the polynomial or Wiener filter. The polynomial coefficients must be determined optimally by minimising mean square error, such that the corresponding polynomial curve with an order of $l$ best fits the given positions. The Polynomial model can be written as

$$X_n = \sum_{j=0}^{l} a_j t_n^j, \qquad (2.1)$$

where {$a_0, \dots, a_l$} are unknown coefficients of the predictor function and $X_n$ is the polynomial position estimation at time step $n$. Therefore, the least square energy function is

$$\mathcal{J} = \frac{1}{n} \sum_{i=1}^{n} (\mathbf{X}_i - \mathbf{y}_i)^2, \qquad (2.2)$$
particle history

where $\mathbf{y}_i$ are $n$ known observed positions (i.e., history) of the last finite frames (in this study, $n \ge 4$). For simplicity, observation starts from the first time step for the rest of derivations :. $i = \{1, \dots, n\}$. The objective is to estimate unknown polynomial coefficients
200 $\hat{\mathbf{a}} = \{a_0, \dots, a_l\}$ by minimising the energy function in (2.2) as

$$
201
\left\{
\begin{array}{l}
\hat{\mathbf{a}} = \arg \min \mathcal{J}(\mathbf{X}, \mathbf{y}), \\
\mathbf{X} = \{X_1, \dots, X_n\}, \\
\mathbf{y} = \{y_1, \dots, y_n\}, \\
\hat{\mathbf{a}} = \{a_0, \dots, a_l\}.
\end{array}
\right.
\tag{2.3}
$$

202 In the finite impulse response (FIR) Wiener filter approach, as a short-term linear prediction
203 model, we first design a linear estimator (filter) for the history of the target particle. Consider
204 the signal $u_n$ (i.e., history) is given to a Wiener filter of order $l$ as

$$
205
X_n = \sum_{j=1}^{l} w_j u_n, \tag{2.4}
$$

206 where $w_j$ are the filter parameters and the filter output is indicated by $X_n$. Similar to the
207 polynomial solution in (2.3), the objective is to find $\{w_1, \dots, w_l\}$ filter parameters that
208 minimise a quadratic energy function with the mean square error (MSE). The resulting
209 Wiener filter in (2.4) is a linear minimum mean square error (LMMSE) estimator. Then we
210 predict the future signal value (i.e., particle position) with the designed filter parameters
211 at $t_{n+1}$. Predicting a signal from its past samples depends on the auto-correlation function
212 in (2.4), or equivalently the signal's bandwidth and power spectrum. Schanz et al. (2016)
213 showed that the Wiener filter is able to predict Lagrangian trajectories by employing the
214 mentioned auto-correlation functions. On the other hand, a Wiener filter can forecast the
215 amplitude of a signal in a short time using a linearly weighted combination of past samples.
216

2.2. *Lagrangian coherent energy function*

217 Particle trajectories in fluid flows are highly influenced by the coherent structures present in
218 the flow. To improve the prediction of particle positions, it is essential to consider the effects
219 of these coherent structures. Lagrangian coherent energy function incorporates information
220 from temporal and local spatial coherent motions to come up with extra constraints for the
221 position-based energy function in 2.2. Each Lagrangian trajectory carries sets of information,
222 including position, velocity, and acceleration. We impose the dynamics of each coherent
223 particle into the prediction function. Consequently, the energy function of each particle
224 will end up with the weighted average of local coherent velocity and acceleration values
225 ($\dot{\mathbf{y}}_{c,n}, \ddot{\mathbf{y}}_{c,n}$), where the subscript $c$ denotes coherent trajectory. The imposed coherent velocity
226 and acceleration terms create additional constraints to (2.2). Therefore, the modified weighted
227 energy function called coherent predictor can be written as

$$
228
\mathcal{J} = \frac{1}{n} \sum_{i=1}^{n} (\mathbf{X}_i - \mathbf{y}_i)^2 + \alpha_1 (\dot{\mathbf{X}}_n - \dot{\mathbf{y}}_{c,n})^2 + \alpha_2 (\ddot{\mathbf{X}}_n - \ddot{\mathbf{y}}_{c,n})^2, \quad (2.5)
$$

229 where $\alpha_1$ and $\alpha_2$ are dimensional velocity and acceleration weights, respectively. The
230 weighted energy function (2.5) contains three terms. The first term is the least mean square
231 minimisation problem of the polynomial predictor based on position history observations $\mathbf{y}_i$.
232 The second and third terms involve neighbouring coherent velocity $\dot{\mathbf{y}}_{c,n}$ and acceleration $\ddot{\mathbf{y}}_{c,n}$
233 observations. The energy function (2.5) adds up different squared quantities (i.e., differences
234 in position, velocity, and acceleration) with distinct dimensions. $\alpha_1$ and $\alpha_2$ values are directly
235 related to the dimensions of the particular problem considered. They will, therefore, change
236 according to the configurations studied. At this stage, it is difficult to determine the behaviour
237 of their optimal values, i.e., leading to the most accurate position prediction according to
238 the different situations. Our first objective is to reformulate the function by explaining it in
relation to the turbulent characteristic scales of the problem so that the weighting parameters become dimensionless. In a second step, in § 4.1, we propose to model their behaviour. To obtain a generic energy function, we sized the variables with the following scales

$$
\begin{align}
\mathbf{X}' &= \frac{\mathbf{X}}{D}, \quad \mathbf{y}' = \frac{\mathbf{y}}{D} \nonumber \\
\dot{\mathbf{X}}' &= \frac{\dot{\mathbf{X}}}{U}, \quad \dot{\mathbf{y}}' = \frac{\dot{\mathbf{y}}}{U} \tag{2.6} \\
\ddot{\mathbf{X}}' &= \frac{\ddot{\mathbf{X}}}{U^2}, \quad \ddot{\mathbf{y}}' = \frac{\ddot{\mathbf{y}}}{U^2}, \nonumber
\end{align}
$$

where *D* is an integral length scale, and *U* is the velocity reference. The integral length scale represents the size of the largest coherent structures in the flow, while *U* represents the characteristic velocity of the flow. These reference scales help to non-dimensionalise the energy function, enabling its application to a wide range of flow configurations. Therefore, the modified energy function becomes

$$
\mathcal{J} = D^2 \left( \frac{1}{n} \sum_{i=1}^{n} (\mathbf{X}'_i - \mathbf{y}_i)^2 \right) + \alpha_1 U^2 (\dot{\mathbf{X}}'_n - \dot{\mathbf{y}}'_{c,n})^2 + \alpha_2 \frac{U^4}{D^2} (\ddot{\mathbf{X}}'_n - \ddot{\mathbf{y}}'_{c,n})^2, \quad (2.7)
$$

which can be simplified as

$$
\mathcal{J}' = \underbrace{\frac{1}{n} \sum_{i=1}^{n} (\mathbf{X}'_i - \mathbf{y}'_i)^2}_{\text{particle history}} + \underbrace{\alpha'_1 (\dot{\mathbf{X}}'_n - \dot{\mathbf{y}}'_{c,n})^2}_{\text{coherent velocity}} + \underbrace{\alpha'_2 (\ddot{\mathbf{X}}'_n - \ddot{\mathbf{y}}'_{c,n})^2}_{\text{coherent acceleration}}, \quad (2.8)
$$

where

$$
\alpha'_{1} = \alpha_{1}(\frac{U}{D})^{2}, \quad \alpha'_{2} = \alpha_{2}(\frac{U}{D})^{4}, \quad \text{and} \quad \mathcal{J}' = \frac{\mathcal{J}}{D^{2}}. \qquad (2.9)
$$

Both non-dimensional $\alpha'_1$ and $\alpha'_2$ weights determine how much velocity and acceleration signals can constrain the energy function minimisation process, knowing that the position history has a weight equal to one in (2.8). If both parameters are set for any range below one, the history will have the most significant impact on the prediction function. The minimum of the energy function can be found by solving $\frac{\partial \mathcal{J}'}{\partial a} = 0$, where its mathematics is addressed in Appendix A. We can use a minimum of four-time step histories of particles to minimise the energy function and predict the next step. The solution is not only smooth on the history of the target particle but also satisfies the local coherent dynamics of the flow. In the worst-case scenario where there is no coherent neighbour information, the prediction function is just a simple polynomial predictor without additional constraints. To obtain two coherent $\dot{\mathbf{y}}_{c,n}$ and $\ddot{\mathbf{y}}_{c,n}$ observations, neighbour trajectories need to be identified as coherent or non-coherent with the target particle. In § 2.2.1, we discuss how to employ LCS metric locally for particle segmentation and coherent neighbour identifications over sparse trajectories.

2.2.1. Local segmentation by FTLE

This section explains segmenting the phase space into different coherent regions, knowing that every single particle is spatially and temporally coherent with a specific cluster of other particles following the same kinematics (see, e.g., Khojasteh et al. 2021). In the present study, our ultimate goal in utilising LCS is to locally determine coherent and non-coherent neighbour trajectories. A global study from Hadjighasem et al. (2017) compared twelve techniques to detect coherent Lagrangian structures in 2D flows. It was found that FTLE
<figure><img src="image_3.png" /><figcaption>Figure 2: Prediction using primary and secondary coherent particles. (a) Primary coherent. (b) Secondary coherent with one time step delay. (c) Secondary coherent with two time step delay.</figcaption></figure>

273 is a simple algorithm with suitable performance in capturing hyperbolic LCS. However, 274 the technique becomes unreliable in elliptic LCSs. In the present study, we are looking for 275 hyperbolic actions (i.e., local separation) between the target particle and its neighbours. 276

276 We adopt a unique approach to compute the FTLE metric for local segmentation. Instead 277 of calculating the global FTLE field, we define a local Eulerian frame around each target 278 particle (see figure 1). Our goal is to classify all neighbours within this area as coherent or 279 non-coherent concerning the target particle. The frame remains fixed during a series of time 280 steps, offering an Eulerian view of neighbourhood behaviour (see figure 2). The advantage 281 of computing local FTLE over a global computation is that it aligns with our objective of 282 finding coherent neighbours both temporally and spatially. There is no need to classify other 283 clusters of coherent motions when predicting the trajectory of the target particle. Such an 284 objective necessitates a local segmentation approach, allowing for a more targeted analysis 285 of the particle’s immediate surroundings, reducing computational and 3D transport barrier 286 quantification complexities and providing a precise understanding of the particle dynamics 287 in the vicinity of the target particle. We employ a local Eulerian frame since trajectories far 288 from the target particle are less likely to be coherent and would only add computational cost. 289 The maximum displacement of the target particle determines the size of the Eulerian frame 290 in each direction. If the 2D (or 3D) displacements are equal in all directions, the shape forms 291 a circle (or sphere) around the target particle. We then employ the FTLE metric to quantify 292 the stretching rate for each neighbouring trajectory and the target particle. 293

294 To compute FTLE, we assume that the local flow map corresponds to the particle trajectory 295 over the finite time from $t_0$ to $T$ (in this study, $(T-t_0) \ge 4\times dt$). We plot the effect of integration 296 time from 4 which is the minimum track length up to 100 in our analysis in Figure 3 for both 297 coherent and non-coherent neighbours. The FTLE value in both cases tend to drop as the 298 integration time increases, however, comparative difference between two plots remain the 299 same. This indicates using FTLE as a diagnostic tool to quantify coherent, and non-coherent 300 is still possible within smaller integration times.

Figure FTLE vs time
301

302

303

304

FTLE examines the deformation of the flow maps of the nearest neighbours over time. Our approach simplifies to computing the spatial displacements between the target particle and its neighbours over a finite time. Essentially, two neighbour trajectories can stretch, contract, or rotate. The flow map $\Phi_{t_0}^T$ of a single particle can be formulated as

305

$$ \Phi_{t_0}^T(\mathbf{x}) : \mathbf{x}(t_0) \rightarrow \mathbf{x}(T), \qquad (2.10) $$

306

307

308

309

310

where $\mathbf{x}(t_0)$ is the starting position of the interval time $T$, and $\mathbf{x}(T)$ is the final position. Let $\mathbf{x}_0^p$ and $\mathbf{x}_0^n$ denote the initial positions of a target and its neighbouring particle. The Lyapunov exponent is a term that highlights regions with the most stretching if the interval time is positive (known as forward). Differences between the flow maps of the target particle $\Phi_{t_0}^T(\mathbf{x}_0^p)$ and its neighbour $\Phi_{t_0}^T(\mathbf{x}_0^n)$ would result in a vector displacement as following

311

$$ \delta\mathbf{x}(t) = \Phi_{t_0}^T(\mathbf{x}_0^p) - \Phi_{t_0}^T(\mathbf{x}_0^n). \qquad (2.11) $$

312

313

314

315

This vector displacement contains transformations between the initial and final positions of two particles. In classic dynamical systems, if we add a perturbation to a particle's starting position $\mathbf{x}_0$, we can Taylor expand how the particle behaves forward in time ($T > 0$). So, (2.11) can be linearised by using the first term of the Taylor series of $\Phi_{t_0}^T(\mathbf{x}_0^n)$ expanded around $\mathbf{x}_0^p$

316

$$ \delta\mathbf{x}(t) \approx \frac{\partial\Phi_{t_0}^T(\mathbf{x}_0^p)}{\partial\mathbf{x}_0} \delta\mathbf{x}_0, \qquad (2.12) $$

317

318

where $\delta\mathbf{x}_0 = \mathbf{x}_0^p - \mathbf{x}_0^n$. The state-transition matrix $\partial\Phi_{t_0}^T(\mathbf{x}_0^p)/\partial\mathbf{x}_0$ is also known as deformation gradient tensor $\nabla\Phi_{t_0}^T(\mathbf{x}_0)$ (Shadden et al. 2005). The magnitude of $\delta\mathbf{x}(t)$ is

320

$$ \begin{aligned} |\delta\mathbf{x}(t)| &= \sqrt{\delta\mathbf{x}(t)\delta\mathbf{x}(t)} \\ &= \sqrt{(\nabla\Phi_{t_0}^T(\mathbf{x}_0)\delta\mathbf{x}_0)(\nabla\Phi_{t_0}^T(\mathbf{x}_0)\delta\mathbf{x}_0)} \\ &= \sqrt{\delta\mathbf{x}_0(\nabla\Phi_{t_0}^T(\mathbf{x}_0)\nabla\Phi_{t_0}^T(\mathbf{x}_0))\delta\mathbf{x}_0}, \end{aligned} \qquad (2.13) $$

321

where we define

322

$$ \Delta = \nabla\Phi_{t_0}^T(\mathbf{x}_0)\nabla\Phi_{t_0}^T(\mathbf{x}_0). \qquad (2.14) $$

323

324

325

326

327

328

329

$\Delta$ is a symmetric positive definite matrix, also known as the right Cauchy-Green deformation tensor (Shadden et al. 2005) with three real and positive eigenvalues in a 3D domain over finite time. The maximum eigenvalue (the largest singular value) of the Cauchy-Green tensor $\lambda_{max}(\Delta)$ shows the maximum amount that can be possibly stretched (expansion or separation) in finite time, representing how the flow field is sensitive to a perturbation. Furthermore, the eigenvector corresponding to $\lambda_{max}(\Delta)$ represents the direction of the separation. Ultimately, the FTLE value $\Lambda_{t_0}^T$ is defined by scaling the magnitude of the maximum displacement as

330

$$ \Lambda_{t_0}^T = \frac{1}{|T - t_0|} \log(\sqrt{\lambda_{max}(\Delta)}) = \frac{1}{|T - t_0|} \log\left(\frac{\delta x(T)}{\delta x(t_0)}\right). \qquad (2.15) $$

331

332

333

334

By using the Lagrangian framework instead of the Eulerian velocity field, one can directly measure the flow map, as the spatial derivative $\delta x(T)/\delta(t_0)$ can be directly obtained from the trajectories. Forward-FTLE corresponds to a positive interval time, where we assess forward flow maps. Conversely, a negative interval time leads to backward-FTLE values, indicating
<table><thead><tr><th>Case</th><th>2D HIT</th><th>3D wake flow</th><th>3D wake flow</th><th>3D wake flow</th></tr></thead><tbody><tr><td>Approach</td><td>DNS</td><td>DNS</td><td>DNS</td><td>4D-PTV</td></tr><tr><td>Reynolds</td><td>3000</td><td>3900</td><td>300</td><td>3900</td></tr><tr><td>Dimension</td><td>2π × 2π</td><td>4D × 2D × 2D</td><td>8D × 8D × 6D</td><td>24D × 14D × 4D</td></tr><tr><td>Grid size</td><td>256 × 256</td><td>308 × 328 × 87</td><td>145 × 145 × 108</td><td>-</td></tr><tr><td>Time step</td><td>0.01 2π/U</td><td>0.00075 D/U</td><td>0.005 D/U</td><td>307 µs</td></tr><tr><td>Camera resolution</td><td>-</td><td>-</td><td>-</td><td>1280 × 800</td></tr><tr><td>Number of particles</td><td>65536</td><td>150000</td><td>150000</td><td>≈ 12000</td></tr><tr><td>Number of snapshots</td><td>1000</td><td>1000</td><td>1000</td><td>4000</td></tr></tbody><tfoot><tr><td colspan="5">Table 1: Details of numerical (DNS) and experimental (EXP) test cases. The grey squares indicate the regions that were recorded. The figures are schematic and not to scale. The flow structures are adopted from Lehmkuhl et al. (2014).</td></tr></tfoot></table>

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

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

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

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

that particles stretching more in backward time accumulate material in forward time. Both backward and forward values provide essential information on the Lyapunov exponent field. Regions that stretch the most in forward and backward times determine stable and unstable manifolds over finite time. High FTLE field values indicate the presence of ridges that divide the local area into different clusters of coherent particles. Particles that start on one side of a ridge tend to stay on the same side, meaning that FTLE ridges act like invariant manifolds that particles cannot pass through. Lower FTLE values imply that neighbouring particles behave similarly, showing no signs of separation from the target particle over the finite time. This process segments the phase space into different coherent regions, allowing for the indexing of a group of neighbour particles as coherent or non-coherent with the target particle.

The coherent predictor indexes neighbour trajectories as coherent or non-coherent within the interval time from $t_1$ to $t_n$. In the present study, these trajectories are called primary spatial/temporal coherent neighbours. As shown in figure 2.a, primary coherent neighbours follow a similar path with the target particle during the same period (same phase). However, a secondary group of coherent particles with a phase delay is still possible. Secondary refers to neighbours whose history was coherent with the target particle with a phase delay. The history of the secondary coherent neighbours is priori knowledge for the target particle. Schematics of secondary coherent neighbours with one and two-time step phase delays are shown in figure 2.b.c. One and two-time step phase delays refer to groups of particles that were spatially located in the neighbourhood of the target particle at $t_{n-1}$ and $t_{n-2}$, respectively. Then the FTLE function determines secondary coherent neighbours between the target particle at $t_n$ and the secondary particles at $t_{n-\tau}$ where $\tau$ is the phase delay. In steady flow scenarios, secondary coherent particles show the exact motion of the target particle; however, as flow unsteadiness increases, the uncertainty of relying on information from secondary coherent neighbours also increases. Further study is needed to track coherent particle clusters in real experiments. For future work, the FTLE function could be replaced with techniques like *Coherent Structure Colouring* (CSC, Martins et al. 2021), which is beyond the scope of the present study. As a result of local segmentation, we can first identify coherent neighbour trajectories and then obtain the weighted averaged values of velocity $\dot{y}_{c,n}$ and acceleration $\ddot{y}_{c,n}$ over coherent neighbours based on their FTLE and inverse of their Euclidean distances to the target particle.
## 3. Numerical and experimental flow configurations

To validate and assess our proposed energy function in turbulent flows, we employed four numerical synthetic and experimental velocimetry measurements. Direct numerical simulations (DNS) were performed for 2D homogeneous isotropic turbulent (HIT) flow at a Reynolds number of 3000, as well as the 3D wake behind a circular cylinder at Reynolds numbers of 3900 and 300. The specifics of these numerical cases can be found in table 1. A variety of flow dynamics, including vortices, saddle points, and shear flows, are present in the 2D-HIT case, complicating the accurate prediction of flow motion (see figure 3. b.c.d). For the 3D wake behind the cylinder at a Reynolds number of 3900, we captured DNS snapshots of the formation region where unstable sideward shears collapse into the wake. This region, denoted by grey squares in table 1, is particularly challenging for identifying true trajectories due to the emergence of complex 3D directional trajectories with high gradients (Khojasteh et al. 2021; Schanz et al. 2020). We also examined the wake behind the cylinder at a lower Reynolds number of 300 to assess the effectiveness of the proposed approach across various Reynolds numbers. In addition to the synthetic analyses, particle tracking experiments of the wake behind a circular cylinder at a Reynolds number equal to 3900 were conducted in a wind tunnel. In the experiment, we captured the vortex-shedding behaviour after the formation region, where large-scale streamwise and spanwise coherent motions are significant. Having all these cases enabled us to evaluate the proposed approach across a large spectrum of complexities in the wake flow.

### 3.1. Direct numerical simulations (DNS)

#### 3.1.1. Homogeneous isotropic turbulence

We transported synthetic particles on 2D-HIT flow obtained from DNS at a Reynolds number equal to 3000. The DNS snapshots were obtained from the FLUID European project (Heitz et al. 2007). The Navier-Stokes equation was solved using incompressible condition ($\nabla \cdot \mathbf{u} = 0$). The vorticity and scalar equations are solved in Fourier space using de-aliased Fourier expansions in two directions with periodic boundary conditions. The time integration is the third-order Runge Kutta scheme. The square non-dimensional domain size of $2\pi \times 2\pi$ was discretised into $256 \times 256$ node elements with periodic boundary conditions in four boundaries. This means that a particle enters in one side of the domain if a particle leaves on the opposite side at the same mirrored position. The dimensional DNS time step was set at $0.01 D/U$, where $D$ is the integral scale and $U$ is the reference velocity. We collected particle transport of 10000 DNS time steps. Figure 3 shows instantaneous velocity and streamlines of the computational domain. We employed 2D-HIT synthetic data as the ground truth particle positions to examine how incorporating Lagrangian coherency can improve prediction accuracy and robustness.

#### 3.1.2. Wake behind a circular cylinder

To evaluate our novel position prediction approach in a 3D complex flow, we employed DNS simulations of the wake behind a smooth cylinder at Reynolds numbers equal to 3900 and 300 (based on the diameter $D$ of the cylinder and the free-stream velocity $U$) computed by an open-access code called *Incompact3d* (Laizet & Li 2011). In both cases, the dimensional DNS time step was set to $0.00075D/U$ and $0.005D/U$, respectively, and 150000 particles were transported every 10 DNS time steps using fourth-order Runge Kutta temporal and trilinear spatial schemes. Particle trajectories in the synthetic case are smooth and predictable when the synthetic temporal resolution is in the same order as the DNS time step due to the short travelling distance between two time steps (less than the Kolmogorov timescale).

For the Reynolds number 3900 case, the flow field was sampled at 1000 time steps for a
a)

<figure><img src="image_8.png" /><figcaption>Figure 3: One snapshot of 2D homogeneous isotropic turbulent (HIT) flow.</figcaption></figure>

domain of $4D \times 2D \times 2D$, with the cylinder located within the domain at (1D, 1D). Details of the simulation and synthetic transport of particles can be found in Khojasteh et al. (2022).

For the Reynolds number 300 case, a cylinder wake flow was generated from an initial

null velocity field within the reference domain of size $20D \times 12D \times 6D$, discretized into

$360 \times 217 \times 108$ nodes. Turbulent structures were induced by imposing an inflow velocity $u_{in}/U = 1$ at each DNS time step. The cylinder was modelled using an Immersed Boundary

Method (Parnaureau et al. 2008, 2004). We imposed free-slip conditions on the y-axis borders

and periodic conditions along the cylinder axis. The outflow was modelled by a first-order advection model. We sampled 1000 time steps of data within a domain of $8D \times 8D \times 6D$,

situated far downstream from the cylinder, which was located at (4D, 4D). These two cases

provide a comprehensive evaluation of the proposed approach, allowing us to examine its

performance under varying flow conditions and Reynolds numbers.

## 3.2. Time-resolved particle tracking velocimetry (4D-PTV)

We performed an experimental study of the cylinder wake flow at a Reynolds number equal to

$3900$ (the same as the synthetic data). Experiments were carried out in the INRAE low-speed

wind tunnel equipped with a centrifugal fan, a diffuser, a plenum chamber with honeycomb

and grids, a contraction section decreasing by four, and an area with transparent walls for

testing (Chandramouli et al. 2019). With the aid of hot wire anemometry, the velocity profile

at the wind tunnel entrance was checked to ensure uniformity. The free-stream turbulence

intensity level was found to be less than 0.2 %. The cross-section of the testing zone is square,

with a width of 28 cm and a length of 100 cm. It has a slightly tilted upper wall to reduce

longitudinal pressure gradients. The setup allowed for continuous flow velocity selection

between 1 and 8 m/s.

We designed an experimental setup for time-resolved volumetric measurement, as shown

in figure 4.a. Four CMOS SpeedSense DANTEC cameras with a resolution of $1280 \times 800$

pixels and a maximum frequency of 3.25 kHz are empowered. We equipped the cameras
<figure><img src="image_9.png" /><figcaption>Figure 4: Experiment setup design for the cylinder wake flow at Reynolds 3900. (a) Schematic of the camera setup design. (b) One snapshot of HFSB tracer particles passing downstream of the cylinder.</figcaption></figure>

with Nikon 105 mm lenses. To enable off-axis acquisitions, we mounted the lenses on Scheimpflug adapters (LaVision GmbH). The first two cameras are positioned in backward light scattering, while the second two cameras receive maximum intensity signal in forward scattering. The calibration error was lower than 0.06 pixel and reduced to 0.04 after the volume self-calibration in which each mm is equal to 6.28 pixels. The volume of interest was 280 mm × 160 mm × 46 mm starting from roughly 4D downstream of the cylinder. The aperture was set at 8 to achieve 46 mm depth of focus. We employed an LED system on the top and a mirror at the bottom of the test section to illuminate the volume of interest. The acquisition was long enough to observe dynamic evolutions of the von Kármán vortex streets downstream of the cylinder.

Helium Filled Soap Bubbles (HFSB) are key ingredients of performing large-scale volumetric experiments in air. Soap bubbles are the only tracers with a size significantly larger than 10 µm, leading to sufficient light scattering in the Large-scale volumetric measurements. Smaller particles, such as oil droplets, can create a very dense particle concentration and follow the flow accurately, but they scatter very little light. For this reason, volumetric measurements using these small seeding particles are restricted to small measurement volumes. Therefore, in the present experiment, the seeding particles were HFSB, resulting in desired intensity signal with appropriate particle size. In a similar experiment, Scarano *et al.* (2015) studied the application of using HFSB in the wake flow past a cylinder in a volume of 20 × 20 × 12 cm³ (4800 cm³). Based on a study carried out by Caridi *et al.* (2016), HFSB was determined to have several orders of magnitude higher intensity than fog droplets. Particle response time is a value that determines how particles follow the flow with fidelity. Time response is directly linked to the particle diameter and mass density discrepancy to the air in wind tunnel experiments. Due to this reason, large fog droplets (> 10 – 20 µm) do not follow the flow with enough fidelity because of their poor time-response value. However, the mass density discrepancy is close to zero since HFSB particles are filled with Helium (lighter than air). As a result, the particle time response for HFSB becomes small enough to follow the flow with fidelity. Accordingly, Scarano *et al.* (2015) reported that the HFSB time-response is maintained well below 100 µs, which means that particles should adequately follow the flow in low-speed experiments. So, large particles with favourable time-response values provide the ability to perform large-scale measurement volumes (Schneiders *et al.* 2016). However, during the experiment, we noticed that bubbles in wind tunnel experiments are limited by three primary factors: generation rate, lifetime, and image glare points. Due to these limitations, HFSB for large-scale volumetric measurements inside the wind tunnel
<figure><img src="image_10.png" /><figcaption>Figure 5: Contour plot of the position estimation RMS error with respect to $\alpha'_1$ and $\alpha'_2$ weighting terms: (a), 2D homogeneous isotropic turbulent flow at a Reynolds number equal to 3000; (b),(c), 3D wake behind a cylinder at Reynolds numbers equal to 3900 and 300, respectively. The synthetic trajectories were generated based on the 1st LPT challenge estimation errors.</figcaption></figure>

leads to low particle concentration (approximately 1 particle/cm³). One of the earliest studies of using HFSB reported ppp < 0.01 for the seeding density, which could only resolve the shedding large-scale structures and was unable to capture turbulent small scales (Scarano et al. 2015). In 2018, Gibeau & Ghaemi (2018) reached to 1.6 particle/cm³ with the idea of having a multi-wing seeding system. The impact of the multi-wing generator on the flow stream is found to be at most 1.9% of the turbulence intensity with a negligible deficit on the mean flow. Gibeau et al. (2020) demonstrated to achieve 0.02 ppp over a volume of 2000 cm³ using a full-scale HFSB generator with 48 nozzles.

In the present study, we placed 50 bubble generator nozzles with airfoil-shaped structures inside the wind tunnel chamber. Figure 4.a shows homogeneous distribution of tracer particles during the experiment. The nozzles were far upstream of the measurement section in the settling chamber to ensure a sufficient number of bubbles were created, and the flow field was not disturbed by the existence of nozzles. The bubble lifetime is very short (around 2–3 minutes) inside the wind tunnel, mainly because they explode by passing through honeycomb layers. To overcome this issue, we injected bubbles inside the chamber for a few minutes when the wind tunnel was off before starting the acquisition. We found that particles larger than three pixels create two glare points on two sides of the bubble when the illumination is LED. This requires more image treatments before running the 4D-PTV algorithm to avoid false particle reconstruction. However, the intensity of two glare points can diffuse and merge if the particle size is around two pixels. Therefore, we adjusted the camera magnification to reach particles with two pixel image sizes on average to surpass the glare point issue.

## 4. Coherent predictor results

We now apply the newly developed prediction approach, based on coherent velocity and acceleration priors as presented in § 2, to numerical and experimental flow configurations described in § 3. The predictor function involves two weighting parameters that need to be optimally tuned. In § 4.1, the behaviour of these parameters is analysed and modelled as a function of the uncertainties on the predicted position, coherent velocity and acceleration. In § 4.2, using the 2D HIT test case, we address how sensitive the proposed coherent predictor and the state-of-the-art approaches are to the change of particle concentration, noise level, and temporal resolution. We then quantify the proposed model’s instantaneous and time-averaged bias error in the 3D wake flow test case context, along with uncertainty quantification using
<figure><img src="image_11.png" /><figcaption>Figure 6: Contour plot of the position estimation RMS error as a function of observations uncertainty levels. White arrows show the changing direction of optimal solutions (i.e., minimum estimation error) by increasing each observation uncertainty.</figcaption></figure>

Monte Carlo simulation in § 4.3. Following this, we present the experimental demonstration of the coherent predictor in § 4.4.

## 4.1. Optimal weighting parameters

In this subsection, we aim to determine the optimal values for the weighting parameters $\alpha'$ by examining their influence on the prediction RMS error. The optimal solution minimises the energy function, resulting in the lowest estimation error. With the non-dimensional energy function (2.8) in hand (see appendix A), we seek to determine the most suitable $\alpha'$ values for optimal trajectory estimation. To achieve this, we conducted predictions across a range of $\alpha'$ weights. As a real experiment involve uncertainties and inaccuracies, it is necessary to incorporate uncertainties into the three observation terms within the energy function, specifically the positions $\mathbf{y}'_i$, the coherent velocity $\dot{\mathbf{y}}_{c,n}$, and the coherent acceleration $\ddot{\mathbf{y}}_{c,n}$.

Using data from the 1st LPT challenge (Sciacchitano et al. 2021; Leclaire et al. 2021) as a reference, we estimated the inaccuracies introduced into the predictor function. This challenge assessed the position estimation accuracy of six state-of-the-art time-resolved tracking algorithms, including our coherency-based track initialisation (Khojasteh et al. 2021), for particle densities from 0.05 to 0.2 ppp. We, therefore, converted averaged RMS dimensional errors into non-dimensionalised position $\epsilon_X/D$, velocity $\epsilon_X/U$, and acceleration $\epsilon_{\ddot{X}}D/U^2$ errors to create the synthetic data at each particle density. For example, the average of the reported RMS position error was 0.005 mm where the integral scale was $D = 10$ mm at the particle density of 0.12 ppp, which gives $\epsilon_X/D = 5.10^{-4}$. The averaged velocity and acceleration errors at the same particle density were found to be $\epsilon_{\ddot{X}}/U = 0.01$ and $\epsilon_{\ddot{X}}D/U^2 = 0.3$, respectively. The acceleration estimation has at least an order of magnitude higher error than other terms. For such a scenario, we can expect the optimal solution should have less $\alpha'_2$ weight than $\alpha'_1$ because of having more acceleration estimation error than other terms. We introduced these errors as input uncertainties into three ground truth Lagrangian trajectories of 2D-HIT and 3D cylinder wake flow cases and then computed the correlation of the final estimation error and $\alpha'$ weights. Data creation of these cases is addressed in § 3.

Next, we plotted the prediction RMS error across a range of $\alpha'$ weights to identify the optimal solution in figure 5. As acceleration errors exceeded other terms in both synthetic data sets, the minimum error corresponded to relatively lower $\alpha'_2$ values (i.e., $\alpha'_1/\alpha'_2 \gg 1$). Both optimal weighting parameters were significantly smaller than $\ll 1$, indicating that the particle history with unit weight is the most valuable signal in such scenarios. Comparing estimated position errors behaviour in 2D-HIT and 3D wake flow cases revealed that all
<figure><img src="image_12.png" /><figcaption>Figure 7: Prediction comparison for single trajectories in 2D-HIT: (a), Smooth trajectory; (b), Highly curved trajectory.</figcaption></figure>

three cases required similar weighting parameters with $\alpha'_1/\alpha'_2 \gg 1$ and small magnitudes $\alpha' \ll 1$. This implies that the proposed energy function remains generic in all three cases.

Its independent nature might be related to the supplementary details provided by coherent neighbours, simplifying prediction complexities by indicating the appropriate direction and acceleration of Lagrangian trajectories. The only clear similarity between all three cases was the identical input uncertainty level. Consequently, we can hypothesise that the optimal values for the weighting parameters are directly connected to the uncertainty levels of the input observations. To validate this hypothesis, we designed three further parametric test cases. In figure 6, we increased the uncertainty level of each parameter solely while other terms were fixed to assess how the optimal solution behaves concerning three position, velocity, and acceleration uncertainties. The optimal solution tends to linearly move away toward higher weighting magnitudes with increased position uncertainty (see figure 6.a). If we fix the position uncertainty, the optimal solution rotates with the same distance around the coordinate centre, depending on which parameter is increased (see figures 6.b and 6.c).

This means that we can model the correlation between two weighting terms. Magnitudes of both gains increase by increased position uncertainty, and the slope for $\alpha_2/\alpha_1$ is a function of relative velocity and acceleration uncertainties.

Assume $\sigma'_p$ is the dimensionless uncertainty of position, and $\sigma'_u$ and $\sigma'_a$ represent the dimensionless uncertainties of velocity and acceleration included in our observation terms.

We employed the same method for non-dimensionalisation as in (2.6), where $\sigma'_p = \sigma_p/D$, $\sigma'_u = \sigma_u/U$, and $\sigma'_a = \sigma_a D/U^2$. Based on figure 6, it can be said that $(\alpha'_1^2 + \alpha'_2^2) = c_1\sigma'_p$ where $c_1$ is a constant, and it can be said that $\alpha'_2/\alpha'_1 = c_2\sigma'_u/\sigma'_a$. Consequently, a model for estimating $\alpha'_1$ and $\alpha'_2$ can be formulated as
<table><caption>Table 2: Scenarios for 2D trajectory assessment.</caption><thead><tr><th>Case</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Particle concentration (ppp)</td><td>0.025</td><td>0.05</td><td>0.08</td><td>0.1</td><td>0.2</td></tr><tr><td>Noise ratio (%)</td><td>0</td><td>15</td><td>30</td><td>45</td><td>60</td></tr><tr><td>Temporal resolution (×dt<sub>DNS</sub>)</td><td>20</td><td>40</td><td>60</td><td>80</td><td>100</td></tr></tbody></table>

$$ \alpha'_{1} = \sqrt{\frac{c_{1}\sigma'_{p}}{1 + \left(\frac{c_{2}\sigma'_{u}}{\sigma'_{a}}\right)^{2}}} \qquad (4.1) $$

$$ \alpha'_{2} = \alpha'_{1} \left( c_{2} \frac{\sigma'_{u}}{\sigma'_{a}} \right). \qquad (4.2) $$

As discussed earlier, both $\alpha'$ weights must be less than one. This model offers an estimation for $\alpha'_1$ and $\alpha'_2$, dependent on $\sigma'_p$, $\sigma'_u$, and $\sigma'_a$. It is important to note, however, that the constants $c_1$ and $c_2$ should be determined through calibration using available synthetic data. This validates our hypothesis that the proposed energy function is correlated with introduced uncertainty levels. Therefore, an appropriate estimation of three uncertainty levels would provide enough information to set weighting parameters, regardless of the flow case (see figure 5). We demonstrate that the proposed energy function is generic, as it remains independent of flow dimensions (2D or 3D), Reynolds number, and flow topology, exhibiting the same estimated position RMS error pattern for different cases.

## 4.2. Sensitivity analyses

Having established the optimal parameters for our innovative physics-based precursor in § 4.1, we first aim to analyse its behaviour for different levels of complexity of trajectories. To perform this analysis, we used 2D synthetic Lagrangian trajectories in the HIT flow case and compared the results of the coherent predictor with those obtained with classical predictor models. As illustrated in figure 3, 2D-HIT carries a range of complex flow motions inside where numerous vortices interact, creating saddle points, strong shears, and vortical structures. Figure 7.a shows one example of two trajectories extracted from 2D-HIT. The first one is for a smooth motion, where all predictor functions can estimate the next position with a negligible bias error. However, the bias error can significantly increase without knowing the surrounding flow motions as soon as a particle starts rotating with both velocity and acceleration changes. As shown in figure 7.b, both lower and higher order polynomial functions, as well as the Wiener filter mispredict high variational dynamics of particles. The third-order polynomial predictor estimated the correct direction of the particle at time step $t_{n+1}$. However, the true direction estimation does not necessarily lead to a proper prediction. Misprediction in the third-order polynomial mainly resulted from the lack of correct acceleration estimation. For such a trajectory, the two other polynomial cases showed less accurate estimations. Unlike other methods, the coherent predictor has correct acceleration and direction estimations due to imposed additional constraints in its energy function. The velocity constraint governs the direction of the estimated position, while the acceleration constraint controls how far or near the prediction can go in the same direction when faced with acceleration variations.
<figure><img src="image_13.png" /><figcaption>Figure 8: Comparison performances of 2D Lagrangian reconstruction in terms of the fraction of true particle detection by changing each characteristic parameter.</figcaption></figure>

In the following, we evaluate the performance of the proposed prediction approach one step further with sequences of particle positions during their two-dimensional simulated trajectories. It allows not only to study the influence of parameters such as particle image concentration, noise level, and the temporal resolution of the image observations but also to compare the results with those obtained with state-of-the-art particle tracking methods. We compared a total of seven techniques: four open-access particle tracking techniques capable of reconstructing 2D Lagrangian trajectories along with three predictors in 2D-STB. The first technique, TracTrac (Heyman 2019), starts with the nearest neighbour approach and then searches for the best forward/backward match from $t_{n-1}$ to $t_{n+1}$. The key aspect of this technique is the use of a predictor function based on the velocity field. TracTrac is capable of reconstructing over 1000 tracks per second with 0.01 pixel resolution accuracy (Heyman 2019) and was evaluated with PIV challenge cases (Kähler et al. 2016). Part2Track (Janke et al. 2020), the second technique, is a 2D polynomial predictive tracking method. This technique can be considered a simplified 2D version of STB, featuring a well-organised and robust code in terms of particle density. We also used 2D Enhanced Track Initialisation (4BE-ETI Clark et al. 2019), representing four-frame-based techniques. 4BE-ETI examines all probabilities around the target particle and generates two consecutive predictions in the following four frames. Any particle close to the predicted position is considered for the prediction of the next frame. A solution is deemed reconstructed if a unique track is identified after four frames. Topology-based particle tracking (T-PT Patel et al. 2018) is the last algorithm for the present assessment. The aforementioned technique generates particle descriptors called feature vectors of nearest neighbour particles for each frame. Subsequently, for each particle, groups of nearest neighbours are stored in the descriptor by binning these particles into a gridded bin, assuming that particles tend to remain within the same bin in the next time step. Eventually, it performs an iterative wrapping scheme to reconstruct Lagrangian trajectories. This method has been evaluated for biological flow motions (Patel et al. 2018).

All the aforementioned techniques were assessed based on three characteristic parameters: particle concentration, temporal resolution and noise level, as detailed in table 2. The aim is to determine how these characteristic parameters cause false Lagrangian reconstruction. When particle density is low, a simple optimisation approach would lead to building true trajectories. However, as the density increases and more particles interact and move close to
<figure><img src="image_14.png" /><figcaption>Figure 9: Normal probability density function (pdf) of particle position estimation error in x direction of four predictor functions.</figcaption></figure>

<table><thead><tr><th>Method</th><th>Fit parameters</th><th>Energy function</th><th>Prediction</th></tr></thead><tbody><tr><td>a) DNS predictor</td><td>-</td><td>-</td><td>y<sub>n+1</sub> = X<sub>DNS</sub> · t<sub>n+1</sub></td></tr><tr><td>b) Polynomial predictor</td><td>X<sub>n</sub> = &sum;<sub>j=0</sub><sup>&ell;</sup> a<sub>j</sub> · t<sub>n</sub><sup>j</sup></td><td>J = <sup>1</sup>&frasl;<sub>n</sub> &sum;<sub>i=1</sub><sup>n</sup> (X<sub>i</sub> - y<sub>i</sub>)<sup>2</sup></td><td>y<sub>n+1</sub> = &sum;<sub>j=0</sub><sup>&ell;</sup> a<sub>j</sub> · t<sub>n+1</sub><sup>j</sup></td></tr><tr><td>c) Wiener filter</td><td>X<sub>n</sub> = &sum;<sub>j=1</sub><sup>&ell;</sup> w<sub>j</sub> · u<sub>n</sub></td><td>J = &sum;<sub>i=1</sub><sup>n</sup> (X<sub>i</sub> - y<sub>i</sub>)<sup>2</sup></td><td>y<sub>n+1</sub> = &sum;<sub>j=2</sub><sup>&ell;+1</sup> w<sub>j</sub> · u<sub>n</sub></td></tr><tr><td>d) Coherent predictor</td><td>J = <sup>1</sup>&frasl;<sub>n</sub> &sum;<sub>i=1</sub><sup>n</sup> [(a<sub>0</sub> + a<sub>1</sub>t<sub>i</sub> + a<sub>2</sub>t<sub>i</sub><sup>2</sup> + a<sub>3</sub>t<sub>i</sub><sup>3</sup>) - y<sub>i</sub>']<sup>2</sup> + a'<sub>1</sub> [(a<sub>1</sub> + 2a<sub>2</sub>t<sub>n</sub> + 3a<sub>3</sub>t<sub>n</sub><sup>2</sup>) - y'<sub>n</sub>']<sup>2</sup> + a'<sub>2</sub> [(2a<sub>2</sub> + 6a<sub>3</sub>t<sub>n</sub>) - y'<sub>n</sub>']<sup>2</sup></td><td>J' = <sup>1</sup>&frasl;<sub>n</sub> &sum;<sub>i=1</sub><sup>n</sup> (X'<sub>i</sub> - y<sub>i</sub}')<sup>2</sup> + a'<sub>1</sub> (X'<sub>n</sub> - y'<sub>c,n</sub>)<sup>2</sup> + a'<sub>2</sub> (X'<sub>n</sub> - y'<sub>c,n</sub>)<sup>2</sup></td><td>y<sub>n+1</sub> = &sum;<sub>j=1</sub><sup>n</sup> a<sub>j</sub> · t<sub>n+1</sub><sup>j</sup></td></tr></tbody><caption>Table 3: Prediction function formulations.</caption></table>

each other, a more sophisticated algorithm is required to detect true tracks. In this context, we define particle density based on the number of particles per pixel (ppp). Adding noise level creates a realistic situation in synthetic image studies since noises are inevitable in experiments due to their nature. Therefore, the algorithm's robustness over different noise levels would provide valuable information to determine which technique is appropriate for a particular experiment. The DNS time step is smaller than the smallest timescale of turbulence, which is the order of the Kolmogorov timescale. This means the transport of particles between two DNS time steps is smooth enough to reach true prediction even with linear extrapolation. However, the experiment time step is multiple times higher than the DNS time step, depending on the acquisition setup. To reach a realistic condition, five time steps starting from every 20 to 100 DNS time steps are considered. All techniques receive the same starting position reconstruction. We increased each characteristic parameter individually within five scenarios (see table 2). Finally, we compared the deviation of the final detected positions with the ground truth data. The trajectory is reconstructed correctly if the estimation deviation is in the same order as the reconstruction accuracy.

Figure 8 shows the fraction of true reconstructed tracks over the total number of particles. By increasing the noise ratio up to 60 % noisy reconstruction, Part2Track, T-PT, and coherent
<figure><img src="image_15.png" /><figcaption>Figure 10: Instantaneous position estimation error averaged in *z* direction. (a) DNS predictor. (b) Polynomial Predictor. (c) Wiener filter. (d) Coherent predictor. Also shown in Supplementary Movie 1.</figcaption></figure>

642 predictor tend to keep their robustness. Meanwhile, other techniques faced significant drops, 643 losing nearly half of the true trajectories. Among all techniques, TracTrac has the most 644 sensitivity to the noise ratio. Figure 8 shows that particle concentration has almost the same 645 impact on all techniques. Their performance decreases in the same order as the concentration 646 increases. By increasing the time step, we are not able to observe small-scale motions of 647 particles; therefore, less temporal information is available. The number of true particles drops 648 severely with increased temporal resolution. In all mentioned techniques, relying only on a 649 single particle as a single signal to find its true Lagrangian motion and losing information 650 in the middle causes more wrong trajectories. However, even a weak signal of coherent 651 particle behaviour would lead to correct direction and prediction. Results showed that when 652 adding spatial and temporal coherent information, the prediction function remained robust 653 for up to 85 % in all situations, while other techniques suffered from a lack of information. 654 In conclusion, our sensitivity analysis reveals the proposed prediction function's robust 655 performance across diverse scenarios. The 2D-HIT evaluation highlights the efficiency of 656 our approach in handling complex flow motions and vortical structures, even in situations 657 with low temporal resolution.

### 4.3. Confidence analyses

658 In this section, we aim to assess the reliability and accuracy of Lagrangian prediction methods 659 by comparing their performance against ground truth 3D wake flow trajectories. We evaluate 660 four schemes, including the DNS predictor, Wiener filter, polynomial predictor, and coherent 661 predictor (see table 3). In the following section, we analyse the position prediction error by 662 decomposing it into two general sources, the bias error and the measurement uncertainties. 663 We calculate the bias error from the deviation of the ground truth trajectories and the 664 predicted positions estimated from perfect observations. Hence, the bias error was caused by 665 a misprediction of flow dynamics. Then, to obtain the measurement uncertainty, we conduct 666 Monte Carlo simulations. 667
<figure><img src="image_16.png" /><figcaption>Figure 11: Schematic of the Monte Carlo uncertainty quantification (MC-UQ) algorithm with distributed input parameters.</figcaption></figure>

### 4.3.1. Bias error

The bias error increases due to the lack of accuracy in the prediction model. We defined the DNS predictor as a reference using the Euler integration method to transport particle positions by the ground truth DNS positions and velocities for every 20 DNS time steps. In such a scenario, we can estimate the error achieved for this sparse temporal resolution and perfect observations, which can be considered the minimum bias. Figure 9 shows the normal probability density function (pdf) of the predicted position errors in the x direction of four schemes. The bell curve distribution is obtained by spatiotemporal averaging of trajectory errors for all particle positions over 200 time steps (i.e., time-averaged). Position error in x direction shows that the deviations of the coherent predictor remain virtually below $\epsilon/D = 0.05$, where D is the cylinder diameter. On the contrary, a significant number of particles are mispredicted in both polynomial and Wiener filter techniques. Similar significant improvements by using the coherent predictor are observed in y and z directions. It is worth noting that the polynomial and Wiener filter predictors require nearly the same amount of computational time, while the coherent predictor demands four times more computational time for a single CPU core calculation. Figure 10 shows the instantaneous projected distribution of the bias error on the xy plane for each predictor function. Interestingly, the prediction error is highly correlated with the flow topology in all schemes.

In this classic wake flow case, the flow topology is captured by preserving the areas where motion gradients are produced. There are three such regions: the boundary layer, the two sideward shear layers at the edge of the formation region, and the wake. The boundary layer refers to the thin layer of fluid that forms along the surface of the cylinder. In this region, Lagrangian trajectories face strong deceleration from the free-stream velocity to zero at the cylinder's surface (due to the no-slip condition) as we move close to the cylinder's leading edge. A strong acceleration gradient also exists in the shear. Inside the wake, on the other hand, the complexity is different and arises from numerous small-scale motions and chaotic behaviour of Lagrangian trajectories. We observe direct signatures of these structures on the instantaneous error. Although the DNS predictor (see figure 10.a) uses known ground truth velocity information, the travelling distance between two observations is large enough to introduce minor errors, particularly in the near wake region where it is chaotic. As shown in figure 10.b, the third-order polynomial has the worst prediction error, which can be up to $\epsilon/D = 0.1$ around the cylinder leading edge and inside the wake
<figure><img src="image_17.png" /><figcaption>Figure 12: Probability distribution around the bias error obtained from the MC-UQ of four predictor functions. (a) Second-order polynomial. (b) Third-order polynomial predictor. (c) Wiener filter. (d) Coherent predictor. Blue dashed lines represent the averaged bias error over 10000 Lagrangian trajectories.</figcaption></figure>

region. The polynomial prediction error distribution is thoroughly shaped by the flow motion (i.e., topology), meaning that any variations inside the flow create a huge estimation error. Overall and local performance of the Wiener filter is better than the polynomial predictor. The Wiener filter succeeded in reducing the prediction error in most of the peak regions (see figure 10.b.c). The coherent predictor showed the best performance locally and globally compared to Wiener and polynomial predictors.

### 4.3.2. Monte Carlo uncertainty quantification (MC-UQ)

For each individual trajectory, we conducted a Monte Carlo simulation to quantify the uncertainty level of the prediction function. In Monte Carlo uncertainty quantification (MC-UQ), the parameter distributions of models are sampled randomly, followed by statistics calculated on the output model distribution (see Joint Committee for Guides in Metrology 2008). Figure 11 schematically shows the probability distribution of the predicted positions as a function of distributed input uncertainties. To start the MC-UQ simulation, we need to quantify the uncertainty level of the input parameters that are fed into the predictor function. In the classic 4D-PTV process (see § 3.2), the prediction function receives positions without trajectories from *iterative particle reconstruction* (IPR). The uncertainty level of IPR can be utilised as an input parameter for uncertainty quantification of the prediction function. This can be obtained from numerical or analytical IPR performance analyses reported by Wieneke (2012) and Jahn *et al.* (2021). However, these values are optimistic and might differ in practical conditions due to the IPR tuning parameters. Instead, similar to § 4.1, we average the estimation errors of all participants in the 1st LPT challenge (Sciacchitano *et al.* 2021; Leclaire *et al.* 2021) to mimic practical and generic uncertainty levels that are introduced to a predictor function. Then, we can quantify the output uncertainty level of predictor functions using MC-UQ with estimated input normal distributions.

MC-UQ process of four prediction functions, second-order polynomial, third-order polynomial, Wiener filter, and coherent predictor, for nearly 10000 trajectories, are shown in figure 12. We need to subtract the bias error from the predicted position error obtained
<figure><img src="image_18.png" /><figcaption>Figure 13: The cylinder wake flow at Reynolds 3900. (a) Experiment normal pdf results of particle position error in the x direction of three predictors. Each predictor is compared with the final optimised positions of STB Davis. (b) Side view of particle trajectories superimposed by vorticity iso-surfaces.</figcaption></figure>

728 from uncertain input parameters to decompose the impact of uncertainty with the impact of flow motion behaviour. MC-UQ requires nearly 10000 iterations per trajectory to achieve a smooth Gaussian distribution in the output. As a result, the total number of 100 million (10000 × 10000) predictions was computed for each predictor function. Figure 12.a.b shows that as the order of magnitudes in the polynomial predictor increases, the uncertainty level in position estimation rises while the bias error decreases. This suggests an inverse correlation between bias error and uncertainty level when using polynomial predictors. It can be concluded that the third-order polynomial maintains better overall accuracy and uncertainty performance compared to the second-order polynomial. The Wiener filter outperformed the third-order polynomial function in terms of both bias error reduction and narrower uncertainty distribution. The coherent predictor’s uncertainty level was found to be minimum and equivalent to that of the second-order polynomial. When considering both bias error and uncertainty level, the coherent predictor provided an optimal balance in comparison to the other predictors.
729
730
731
732
733
734
735
736
737
738
739
740
741

4.4. *Experimental demonstration*

742
743 To quantify the results of different schemes, we compared predictions with optimised positions obtained from STB Davis 10 (LaVision GmbH). Optimised positions refer to 744 the ones acquired after the minimisation process, where the camera image residuals reach 745 a minimum, commonly known as shaking. As a result of the experiment, STB successfully 746 constructed nearly 12000 particles. Lagrangian trajectories of the current experiment with 747 superimposed vorticity iso-surfaces are shown in figure 13.b. We successfully reconstructed 748 long trajectories with a duration of up to 300 time steps. The experimental results presented 749 in this study successfully captured the complex wake structures, such as braid vortices and 750 von Kármán vortex streets, which develop after the formation region behind the circular 751 cylinder. The presence of these flow structures highlights the rich complexity of the flow 752 field and emphasises the need for accurate prediction methods. Even in the instantaneous 753 snapshot of the experiment (refer to figure 4.b), coherent motions formed by passive transport 754 barriers (Haller 2023) are readily visible to the naked eye. This observation from raw images 755 reveals that neighbour coherent particles follow similar dynamics and collectively contribute 756 to the formation of large-scale flow patterns with respecting the passive transport barriers. 757 Therefore, incorporating the dynamics of these coherent neighbours can indeed contribute 758 to improved predictions of Lagrangian trajectories.
759

760 We compared three techniques, polynomial, Wiener filter, and coherent predictors, with
<figure><img src="image_19.png" /><figcaption>Figure 14: Correlation between the magnitude of acceleration and the estimation error.</figcaption></figure>

final optimised positions. The deviation of position estimated of each technique is shown in figure 13.a. The distribution shows that the coherent predictor has more accurate estimations within 1 pixel deviation from the optimised positions. Position estimations of the Wiener filter and coherent predictors stay below 2.5 pixels deviation for nearly all particles. On the contrary, the polynomial predictor has maximum deviation. It should be noted that the performance of the Wiener filter is biased in this comparison since the optimised positions (after shaking) are obtained Wiener predictor. We had to respect the original STB approach to obtain the optimised positions, which is why the Wiener filter performs better in this case than the synthetic case.

We can further assess the distribution of the estimation deviation as a function of Lagrangian acceleration. Figure 14 shows the correlation of predictor functions with acceleration. In agreement with the bell curve distribution in figure 13, the coherent predictor has less deviation than the other two techniques. Both polynomial and Wiener filter predictors show scattered correlations with respect to acceleration. The linear regression model for the polynomial predictor (the black dashed line in figure 14.a) tends to explode as Lagrangian acceleration increases. Wiener filter succeeded in slightly controlling the regression slope. Coherent predictor, on the other hand, tends to fit inside narrow upper and lower bounds along with the regression line. This implies that the coherent predictor has a strong correlation with Lagrangian acceleration and effectively controls error amplification as the acceleration increases.

## 5. Conclusions and outlook

In this study, we addressed the challenge of precise motion estimation of Lagrangian trajectories in fluid dynamics, particularly in the presence of temporal resolution limitations and high motion gradients. The main idea is to shift from traditional single-particle treatment to groups of coherent particles as they share the same dynamics. Inspired by the use of advective Lagrangian coherent structures, we developed a new generic energy function called the *coherent predictor* that takes into account the history of trajectories and local spatial and temporal coherent motions for less biased and lower uncertainty predictions compared to existing methods. We perform a local segmentation approach using FTLE to quantify primary and secondary coherent neighbours. Primary neighbours follow the target particle's path, while secondary neighbours exhibit a phase delay and are ahead in their history. Both offer valuable information: primary neighbours provide immediate insights, and secondary neighbours offer prior knowledge about the target particle's future.
794

795

796

797

798

799

800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

815

816

817

818

819

820

821

To assess the proposed approach, we performed the synthetic analysis of the 2D-HIT flow at a Reynolds number of 3000 and the wake behind a smooth cylinder at Reynolds numbers equal to 3900 and 300, with the same measurement input uncertainty levels. The optimal solution of the energy function for all synthetic cases showed similar weighting configurations. The similar behaviour in all test cases led us to the conclusion that modelling the weighting parameters based on the measurement uncertainty is a feasible investigation. We performed further parametric studies and quantified a model that receives the measurement uncertainties to design the optimal energy function. Therefore, the optimal solution of the energy function is generic and applicable to other experimental test cases when the input uncertainty levels are known. The proposed coherent predictor outperformed the recent predictor functions by having a lower bias error, particularly in complex regions. The polynomial predictor showed maximum deviation with the ground truth data (i.e., bias error). To quantify the output uncertainty level of the proposed approach, we performed Monte Carlo simulations. The coherent predictor showed narrow output uncertainty distribution compared to the Wiener filter and third-order polynomial predictors. We also performed the 4D-PTV experiment of the wake far behind a smooth cylinder at a Reynolds number of 3900. The predicted positions from the coherent predictor showed minimum deviation from the optimised positions compared with other predictor functions. It was found that the estimation error has a direct relation with Lagrangian acceleration. We found that the flow topology highly impacts the estimation error. In the synthetic case, we found visible signatures of three main topologies, the leading edge boundary layer, sideward shears, and the vortex formation zone on the prediction error. In the experiment, we capture further downstream structures where braid vortices and von Kármán vortex streets exhibit. These mentioned structures are featured by high acceleration and 3D directional motions. This independence in the cost function's behaviour might be attributed to the additional information provided by the coherent neighbours. In other words, coherent motions simplify the complexities associated with the prediction process by offering accurate direction and acceleration for Lagrangian trajectories.

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

We believe that coherent prediction holds great potential for a variety of Lagrangian analyses in turbulence. It can contribute to studies requiring accurate, complex and long trajectory reconstructions, such as studies in Lagrangian physics and statistics (see, e.g., Viggiano et al. 2021). Recent studies have focused on converting Lagrangian trajectories into Eulerian vector fields using data assimilation or interpolation techniques (Jeon et al. 2022; Godbersen & Schröder 2020). These studies typically search for a method to interpolate groups of trajectories within a cell to estimate Eulerian velocity. However, not all neighbour trajectories are coherent and interpolating them could lead to an over-smooth reconstruction through the transport barrier. The proposed approach can significantly enhance such research by interpolating only coherent velocity fields. In addition, coherent prediction can also be employed in the transport of synthetic particles from Eulerian numerical simulation with sparse temporal resolution (see, e.g., Sebille et al. 2018). Furthermore, knowing that tracer particles respect advective transport barriers, interested readers can implement other LCS stretching diagnostics (Haller 2023). Finally, local segmentation of Lagrangian trajectories is an intriguing topic that could be usefully explored in future research to effectively switch from a single particle treatment to groups of coherent particles.

838

839

840

841

842

843

**Acknowledgements** The authors gratefully acknowledge Sylvain Laizet from Imperial College London for the 3D wake DNS simulation. We also thank the LPT challenge committee, particularly Andrea Sciacchitano from TU Delft University, for the time-resolved LPT challenge results. Special thanks go to Philippe Georgeault and Johan Carlier, both from INRAE. Philippe Georgeault contributed to the experimental setup, and Johan Carlier provided the 2D DNS simulation data from the FLUID project.
**Funding** For the 3D wake DNS simulation, this work was supported by the EPSRC, providing computational time on the UK supercomputing facility ARCHER2 via the UK Turbulence Consortium (grant number EP/R029326/1). Additionally, the 2D HIT DNS simulation was funded by the FLUID European project under grant agreement ID 513663.

**Declaration of interests.** The authors report no conflict of interest.

**Data availability statement.** The numerical 3D data supporting the findings of this study are openly available in the Data INRAE repository at https://doi.org/10.15454/GLNRHK. The experimental data supporting the findings of this study are available from the corresponding author upon reasonable request.

**Author ORCIDs.** AR. Khojasteh, https://orcid.org/0000-0002-3545-8391; D. Heitz, https://orcid.org/0000-0001-6295-2822

**Author contributions.** AR. Khojasteh, conceptualisation, methodology, computation, and authoring the original draft; D. Heitz, conceptualisation, methodology, supervision, reviewing, and editing.

**Appendix A.**

This appendix describes how to minimise the non-dimensional energy function of the coherent predictor. We assume that the solution is a third-order polynomial with unknown $a_i$ coefficients as,

$$ \sum_{i=1}^{n} \left( \sum_{j=0}^{l} a_j \cdot t_i^j \right) = \sum_{i=1}^{n} y_i. \qquad (A 1) $$

We can expand the solution coefficients into the energy function as follows,

$$ \mathcal{J} = \frac{1}{n} \sum_{i=1}^{n} \left[ \left( a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3 \right) - \mathbf{y}'_i \right]^2 + \alpha'_1 \left[ \left( a_1 + 2a_2 t_i + 3a_3 t_i^2 \right) - \dot{\mathbf{y}}'_i \right]^2 + \alpha'_2 \left[ \left( 2a_2 + 6a_3 t_i \right) - \ddot{\mathbf{y}}'_i \right]^2. \qquad (A 2) $$

If all three terms of (2.5) have exact weights, the linear solution of the prediction function is in three sets of equations as

$$ \left\{ \begin{array}{l} \sum_{j=0}^{l} a_j \cdot t_i^j = \mathbf{y}'_i, \quad i = 1, \dots, n, \\ \\ \sum_{j=1}^{l} a_j \cdot (j) \cdot t_n^{j-1} = \dot{\mathbf{y}}'_{c,n}, \\ \\ \sum_{j=2}^{l} a_j \cdot (j) \cdot (j-1) \cdot t_n^{j-2} = \ddot{\mathbf{y}}'_{c,n}. \end{array} \right. \qquad (A 3) $$

The first sets are rows of particle position history for $n$ time step observations. The second and third sets are additional coherency based constraints. Therefore, the solution for the energy function in (A 3) is not only smooth on the history of the target particle but also satisfies local coherent dynamics of the flow. The minimized solution should satisfy the partial derivative of the energy function over $a_i$ coefficients as,

$$ \frac{\partial \mathcal{J}}{\partial a_i} = 0. \qquad (A 4) $$
874 Partials derivate of $a_0$ is given by,

$$ 
\begin{aligned}
\frac{\partial \mathcal{J}}{\partial a_0} &= \frac{2}{n} \sum_{i=1}^{n} \left[ (a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3) - \mathbf{y}_i \right] + 0 + 0 \\
&\qquad \sum_{i=1}^{n} \left[ a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3 \right] = \mathbf{y}_i.
\end{aligned}
\quad (A5) $$

875

Partial derivative of $a_1$ is,

$$ 
\begin{aligned}
\frac{\partial \mathcal{J}}{\partial a_1} &= \frac{2}{n} \sum_{i=1}^{n} \left[ (a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3) - \mathbf{y}_i \right] \cdot t_i \\
&\quad + 2 \cdot \left[ (a_1 + 2a_2 t_n + 3a_3 t_n^2) - \dot{\mathbf{y}}_{c,n} \right] + 0.
\end{aligned}
\quad (A6) $$

877

Partial derivative of $a_2$ is,

$$ 
\begin{aligned}
\frac{\partial \mathcal{J}}{\partial a_2} &= \frac{2}{n} \sum_{i=1}^{n} \left[ (a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3) - \mathbf{y}_i \right] \cdot t_i^2 \\
&\quad + 2 \cdot \left[ (a_1 + 2a_2 t_i + 3a_3 t_i^2) - \dot{\mathbf{y}}_i \right] \cdot (2t_i) \\
&\quad + 2 \cdot [(2a_2 + 6a_3 t_i) - \ddot{\mathbf{y}}_i] \cdot (2).
\end{aligned}
\quad (A7) $$

879

Partial derivative of $a_3$ is,

$$ 
\begin{aligned}
\frac{\partial \mathcal{J}}{\partial a_3} &= \frac{2}{n} \sum_{i=1}^{n} \left[ (a_0 + a_1 t_i + a_2 t_i^2 + a_3 t_i^3) - \mathbf{y}_i \right] \cdot t_i^3 \\
&\quad + 2 \cdot \left[ (a_1 + 2a_2 t_i + 3a_3 t_i^2) - \dot{\mathbf{y}}_i \right] \cdot (3t_i^2) \\
&\quad + 2 \cdot [(2a_2 + 6a_3 t_i) - \ddot{\mathbf{y}}_i] \cdot (6t_i).
\end{aligned}
\quad (A8) $$

881

882 To linearise the solution, we set each partial derivative to zero leading. A series of partial
derivatives can be written in a matrix form as,

883

$$ 
\begin{aligned}
& \left( \frac{1}{n} \sum_{i=1}^{n} \begin{bmatrix} t_i^6 & t_i^5 & t_i^4 & t_i^3 \\ t_i^5 & t_i^4 & t_i^3 & t_i^2 \\ t_i^4 & t_i^3 & t_i^2 & t_i \\ t_i^3 & t_i^2 & t_i & 1 \end{bmatrix} + \alpha_1' \begin{bmatrix} 9t_n^4 & 6t_n^3 & 3t_n^2 & 0 \\ 6t_n^3 & 4t_n^2 & 2t_n & 0 \\ 3t_n^2 & 2t_n & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} + \alpha_2' \begin{bmatrix} 36t_n^2 & 12t_n & 0 & 0 \\ 12t_n & 4 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \right) \begin{bmatrix} a_3 \\ a_2 \\ a_1 \\ a_0 \end{bmatrix} \\
&= \left( \frac{1}{n} \sum_{i=1}^{n} \begin{bmatrix} \mathbf{y}_i' t_i^3 \\ \mathbf{y}_i' t_i^2 \\ \mathbf{y}_i' t_i \\ \mathbf{y}_i \end{bmatrix} + \alpha_1' \begin{bmatrix} 3\dot{\mathbf{y}}_{c,n} t_n^2 \\ 2\dot{\mathbf{y}}_{c,n} t_n \\ \dot{\mathbf{y}}_{c,n} \\ 0 \end{bmatrix} + \alpha_2' \begin{bmatrix} 6\ddot{\mathbf{y}}_{c,n} t_n \\ 2\ddot{\mathbf{y}}_{c,n} \\ 0 \\ 0 \end{bmatrix} \right).
\end{aligned}
\quad (A9) $$

884

885 On the other hand, (A 9) is in the form of $AX = B$. So $a_i$ solutions can be solved by
$886 X = A^{-1}B$.

REFERENCES

887 ALIPOUR, MOBIN, DE PAOLI, MARCO, GHAEMI, SINA & SOLDATI, ALFREDO 2021 Long non-axisymmetric

888 fibres in turbulent channel flow. *J. Fluid Mech.* **916**, 1–32.
889 BALASURIYA, SANJEEVA, KALAMPATTEL, RAHUL & OUELLETTE, NICHOLAS T. 2020 Hyperbolic
890 neighbourhoods as organizers of finite-time exponential stretching. *J. Fluid Mech.* **807**, 509–545.

891 CARIDI, GIUSEPPE CARLO ALP, RAGNI, DANIELE, SCIACCHITANO, ANDREA & SCARANO, FULVIO 2016 HFSB-seeding for large-scale tomographic PIV in wind tunnels. *Exp. Fluids* **57** (12), 1–13.

892 CHANDRAMOULI, PRANAV, MEMIN, ETIENNE, HEITZ, DOMINIQUE & FIABANE, LIONEL 2019 Fast 3D flow reconstructions from 2D cross-plane observations. *Exp. Fluids* **60** (2), 1–27.

895 CLARK, ALICIA M., MACHICOANE, NATHANAEL & ALISEDA, ALBERTO 2019 A quantitative study of track initialization of the four-frame best estimate algorithm for three-dimensional Lagrangian particle tracking. *Meas. Sci. Technol.* **30**, 1–7.

898 EISMA, JERKE, WESTERWEEL, JERRY & VAN DE WATER, WILLEM 2021 Do coherent structures organize scalar mixing in a turbulent boundary layer? *J. Fluid Mech.* **929**, 1–19.

900 GIBEAU, BRADLEY & GHAEMI, SINA 2018 A modular, 3D-printed helium-filled soap bubble generator for large-scale volumetric flow measurements. *Exp. Fluids* **59** (12), 1–11.

902 GIBEAU, BRADLEY, GINGRAS, DREW & GHAEMI, SINA 2020 Evaluation of a full-scale helium-filled soap bubble generator. *Exp. Fluids* **61** (2), 1–18.

904 GODBERSEN, PHILIPP & SCHRÖDER, ANDREAS 2020 Functional binning: Improving convergence of Eulerian statistics from Lagrangian particle tracking. *Meas. Sci. Technol.* **31** (9), 1–13.

906 GOLD, THOMAS, REITERER, KEVIN, WORF, DOMINIK, KHOSRONEJAD, ALI, HABERSACK, HELMUT & SINDELAR, CHRISTINE 2023 Lagrangian particle tracking velocimetry investigation of vortex shedding topology for oscillating heavy spherical pendulums underwater. *J. Fluid Mech.* **960**, 1–14.

909 GREEN, MELISSA A., ROWLEY, CLARENCE & HALLER, GEORGE 2007 Detection of Lagrangian coherent structures in three-dimensional turbulence. *J. Fluid Mech.* **572**, 111–120.

911 HADJIGHASEM, ALIREZA, FARAZMAND, MOHAMMAD, BLAZEVSKI, DANIEL, FROYLAND, GARY & HALLER, GEORGE 2017 A critical comparison of Lagrangian methods for coherent structure detection. *Chaos* **27** (5), 1–27.

914 HALLER, GEORGE 2015 Lagrangian coherent structures. *Annu. Rev. Fluid Mech.* **47**, 137–162.

915 HALLER, GEORGE 2023 *Transport Barriers and Coherent Structures in Flow Data*. Cambridge University Press.

917 HALLER, GEORGE, KATSANOULIS, STERGIOS, HOLZNER, MARKUS, FROHNAPFEL, BETTINA & GATTI, DAVIDE 2020 Objective barriers to the transport of dynamically active vector fields. *J. Fluid Mech.* **905**, 1–52.

919 HALLER, GEORGE & YUAN, GUOCHENG 2000 Lagrangian coherent structures and mixing in two-dimensional turbulence. *Physica D* **147** (3-4), 352–370.

921 HEITZ, DOMINIQUE, CARLIER, JOHAN & ARROYO, GEORGES 2007 Deliverable 5.4 - Final report on the evaluation of the tasks of the workpackage 2 - FLUID Image analysis and Description (FLUID) - Project No. FP6-513663 - FET Open - Domain, Priority IST. *Tech. Rep.* Cemagref. (hal-03623191).

924 HEYMAN, JORIS 2019 TracTrac: A fast multi-object tracking algorithm for motion estimation. *Comput Geosci* **128**, 11–18.

926 JAHN, TOBIAS, SCHANZ, DANIEL & SCHRÖDER, ANDREAS 2021 Advanced iterative particle reconstruction for Lagrangian particle tracking. *Exp. Fluids* **62** (8), 1–24.

928 JANKE, THOMAS, SCHWARZE, RÜDIGER & BAUER, KATRIN 2020 Part2Track: A MATLAB package for double frame and time resolved Particle Tracking Velocimetry. *SoftwareX* **11**, 1–6.

930 JEON, YOUNG JIN, MÜLLER, MARKUS & MICHAELIS, DIRK 2022 Fine scale reconstruction (VIC#) by implementing additional constraints and coarse-grid approximation into VIC+. *Exp. Fluids* **63** (4), 1–24.

933 KÄHLER, CHRISTIAN J., ASTARITA, TOMMASO, VLACHOS, PAVLOS P., SAKAKIBARA, JUN, HAIN, RAINER, DISCETTI, STEFANO, LA FOY, RODERICK & CIERPKA, CHRISTIAN 2016 Main results of the 4th International PIV Challenge. *Exp. Fluids* **57** (6), 1–71.

936 KHOJASTEH, ALI RAHIMI, LAIZET, SYLVAIN, HEITZ, DOMINIQUE & YANG, YIN 2022 Lagrangian and Eulerian dataset of the wake downstream of a smooth cylinder at a Reynolds number equal to 3900. *Data Brief* **40**, 1–9.

939 KHOJASTEH, ALI RAHIMI, YANG, YIN, HEITZ, DOMINIQUE & LAIZET, SYLVAIN 2021 Lagrangian coherent track initialization. *Phys. Fluids* **33** (9), 095113.

941 KIM, DONG, SCHANZ, DANIEL, NOVARA, MATTEO, SEO, HYUNDUK, KIM, YOUNGWOO, SCHRÖDER, ANDREAS & KIM, KYUNG CHUN 2022 Experimental study of turbulent bubbly jet. Part 1. Simultaneous measurement of three-dimensional velocity fields of bubbles and water. *J. Fluid Mech.* **941**, 1–35.

945 KNOPP, TOBIAS, REUTHER, NICO, NOVARA, MATTEO, SCHANZ, DANIEL, SCHÜLEIN, ERICH, SCHRÖDER,
946 ANDREAS & KÄHLER, CHRISTIAN J. 2021 Experimental analysis of the log law at adverse pressure gradient. *J. Fluid Mech.* **918**, 1–32.

948 LAIZET, SYLVAIN & LI, NING 2011 Incompact3d: A powerful tool to tackle turbulence problems with up to O(105) computational cores. *Int J Numer Methods Fluids* **67**, 1735–1757.

950 LECLAIRE, BENJAMIN, MARY, IVAN, LIAUZUN, CÉDRIC, PÉRON, STÉPHANIE, SCIACCHITANO, ANDREA, SCHRÖDER, ANDREAS, CORNIC, PHILIPPE & CHAMPAGNAT, FRÉDÉRIC 2021 First challenge on Lagrangian Particle Tracking and Data Assimilation: datasets description and planned evolution to an open online benchmark. In *14th International Symposium on Particle Image Velocimetry*, pp. 1–2. Chicago, USA.

955 LEHMKUHL, ORIOL, RODRÍGUEZ, IVETTE, BORRELL, RICARD, CHIVA, JORGE & OLIVA, ASSENSI 2014 Unsteady forces on a circular cylinder at critical Reynolds numbers. *Phys. Fluids* **26** (12).

957 MACMILLAN, THEODORE & RICHTER, DAVID H. 2021 The most robust representations of flow trajectories are Lagrangian coherent structures. *J. Fluid Mech.* **927**, 1–16.

959 MARTINS, FLAVIO A C., SCIACCHITANO, ANDREA & RIVAL, DAVID E. 2021 Detection of vortical structures in sparse Lagrangian data using coherent-structure colouring. *Exp. Fluids* **62**, 1–15.

961 MASUMOTO, YUKIO, SASAKI, HIDEHARU, KAGIMOTO, TAKASHI, KOMORI, NOBUMASA, ISHIDA, AKIO, SASAI, YOSHIKAZU, MIYAMA, TORU, MOTOI, TATSUO, MITSUDERA, HUMIO, TAKAHASHI, KEIKO, SAKUMA, HIROFUMI & YAMAGATA, TOSHIO 2004 A Fifty-Year Eddy-Resolving Simulation of the World Ocean-Preliminary Outcomes of OFES (OGCM for the Earth Simulator). *J. earth simul.* **1**, 35–56.

965 JOINT COMMITTEE FOR GUIDES IN METROLOGY, JCGM 101 2008 Evaluation of measurement data-Supplement 1 to the "Guide to the expression of uncertainty in measurement"-Propagation of distributions using a Monte Carlo method. *Tech. Rep.*.

968 MOWLAVI, SAVIZ, SERRA, MATTIA, MAIORINO, ENRICO & MAHADEVAN, L 2022 Detecting Lagrangian coherent structures from sparse and noisy trajectory data. *J. Fluid Mech.* **948**, 1–34.

970 PARNAUDEAU, PHILIPPE, CARLIER, JOHAN, HEITZ, DOMINIQUE & LAMBALLAIS, ERIC 2008 Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900. *Phys. Fluids* **20** (8), 85101.

973 PARNAUDEAU, PHILIPPE, LAMBALLAIS, ERIC, HEITZ, DOMINIQUE & SILVESTRINI, JORGE H. 2004 Combination of the Immersed Boundary Method with Compact Schemes for DNS of Flows in Complex Geometry. In *Combination of the Immersed Boundary Method with Compact Schemes for DNS of Flows in Complex Geometry*, p. w. Springer, Dordrecht.

977 PATEL, MOHAK, LEGGETT, SUSAN E., LANDAUER, ALEXANDER K., WONG, IAN Y. & FRANCK, CHRISTIAN 2018 Rapid, topology-based particle tracking for high-resolution measurements of large complex 3D motion fields. *Sci. Rep.* **8**, 1–14.

980 QIN, XUERONG, VAN SEBILLE, ERIK & SEN GUPTA, ALEXANDER 2014 Quantification of errors induced by temporal resolution on Lagrangian particles in an eddy-resolving model. *Ocean Model.* **76**, 20–30.

982 SCARANO, FULVIO, GHAEMI, SINA, CARIDI, GIUSEPPE CARLO ALP, BOSBACH, JOHANNES, DIERKSHEIDE, UWE & SCIACCHITANO, ANDREA 2015 On the use of helium-filled soap bubbles for large-scale tomographic PIV in wind tunnel experiments. *Exp. Fluids* **56**, 1–12.

985 SCHANZ, DANIEL, GESEMANN, SEBASTIAN & SCHRÖDER, ANDREAS 2016 Shake-The-Box: Lagrangian particle tracking at high particle image densities. *Exp. Fluids* **57**, 1–27.

987 SCHANZ, DANIEL, NOVARA, MATTEO & SCHRÖDER, ANDREAS 2020 Shake-The-Box particle tracking with variable time-steps in flows with high velocity range (VT-STB). In *3rd Workshop and 1st Challenge on Data Assimilation & CFD Processing for PIV and Lagrangian Particle Tracking*, pp. 1–2. virtual online meeting.

991 SCHANZ, DANIEL, SCHRÖDER, ANDREAS & GESEMANN, SEBASTIAN 2014 'Shake The Box' - a 4D PTV algorithm: Accurate and ghostless reconstruction of Lagrangian tracks in densely seeded flows. In *17th International Symposium on Applications of Laser Techniques to Fluid Mechanics*, pp. 1–14. Lisbon, Portugal.

995 SCHANZ, DANIEL, SCHRÖDER, ANDREAS, GESEMANN, SEBASTIAN, MICHAELIS, DIRK & WIENEKE, BERNHARD 2013 'Shake The Box': A highly efficient and accurate Tomographic Particle Tracking Velocimetry (TOMO-PTV) method using prediction of particle positions. In *10th International Symposium on Particle Image Velocimetry*, pp. 1–13. Delft, The Netherlands.

999 SCHNEIDERS, JAN F.G., CARIDI, GIUSEPPE CARLO ALP, SCIACCHITANO, ANDREA & SCARANO, FULVIO 2016 Large-scale volumetric pressure from tomographic PTV with HFSB tracers. *Exp. Fluids* **57**, 1–8.

1001 SCHRÖDER, ANDREAS & SCHANZ, DANIEL 2023 3D Lagrangian Particle Tracking in Fluid Mechanics. *Annu. Rev. Fluid Mech.* **55**, 511–540.
1003 SCHRÖDER, ANDREAS, SCHANZ, DANIEL, GEISLER, REINHARD, GESEMANN, SEBASTIAN & WILLERT, CHRISTIAN E. 2015a Near-wall turbulence characterization using 4D-PTV Shake-The-Box. In *11th International Symposium on Particle Image Velocimetry*, pp. 1–14. California, USA.

1004
1005 SCHRÖDER, ANDREAS, SCHANZ, DANIEL, MICHAELIS, DIRK, CIERPKA, CHRISTIAN, SCHARNOWSKI, SVEN & KÄHLER, CHRISTIAN J. 2015b Advances of PIV and 4D-PTV "shake-The-Box" for Turbulent Flow Analysis -the Flow over Periodic Hills. *Flow Turbul Combust* **95** (2-3), 193–209.

1009 SCIACCHITANO, ANDREA, LECLAIRE, BENJAMIN & SCHRÖDER, ANDREAS 2021 Main results of the first Lagrangian Particle Tracking Challenge. In *14th International Symposium on Particle Image Velocimetry*,, vol. 1. Chicago, USA.

1012 SEBILLE, ERIK VAN, GRIFFIES, STEPHEN M., ABERNATHEY, RYAN, ADAMS, THOMAS P., BERLOFF, PAVEL, BIASTOCH, ARNE, BLANKE, BRUNO, CHASSIGNET, ERIC P., CHENG, YU, COTTER, COLIN J., DELEERSNIJDER, ERIC, DÖÖS, KRISTOFER, DRAKE, HENRI F., DRIJFHOUT, SYBREN, GARY, STEFAN F., HEEMINK, ARNOLD W., KJELLSSON, JOAKIM, KOSZALKA, INGA MONIKA, LANGE, MICHAEL, LIQUE, CAMILLE, MACGILCHRIST, GRAEME A., MARSH, ROBERT, MAYORGA ADAME, C. GABRIELA, MCADAM, RONAN, NENCIOLI, FRANCESCO, PARIS, CLAIRE B., PIGGOTT, MATTHEW D., POLTON, JEFF A., RÜHS, SIREN, SHAH, SYED H.A.M. A M, THOMAS, MATTHEW D., WANG, JINBO, WOLFRAM, PHILLIP J., ZANNA, LAURE & ZIKA, JAN D. 2018 Lagrangian ocean analysis: Fundamentals and practices. *Ocean Model.* **121**, 49–75.

1021 SHADDEN, SHAWN C., LEKIEN, FRANCOIS & MARSDEN, JERROLD E. 2005 Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. *Physica D* **212** (3-4), 271–304.

1024 TAN, SHIYONG, SALIBINDLA, ASHWANTH, MASUK, ASHIK ULLAH MOHAMMAD & NI, RUI 2020 Introducing OpenLPT: new method of removing ghost particles and high-concentration particle shadow tracking. *Exp. Fluids* **61** (2), 1–16.

1027 VALDIVIESO DA COSTA, MARIA & BLANKE, BRUNO 2004 Lagrangian methods for flow climatologies and trajectory error assessment. *Ocean Model.* **6** (3-4), 335–358.

1029 VIGGIANO, BIANCA, BASSET, THOMAS, SOLOVITZ, STEPHEN, BAROIS, THOMAS, GIBERT, MATHIEU, MORDANT, NICOLAS, CHEVILLARD, LAURENT, VOLK, ROMAIN, BOURGOIN, MICKAËL & CAL, RAÚL BAYOÁN 2021 Lagrangian diffusion properties of a free shear turbulent jet. *J. Fluid Mech.* **918**, 1–30.

1032 WIENEKE, BERNHARD 2012 Iterative reconstruction of volumetric particle distribution. *Meas. Sci. Technol.* **24** (2), 1–14.

1034 YANG, YIN & HEITZ, DOMINIQUE 2021 Kernelized Lagrangian particle tracking. *Exp. Fluids* **62** (12), 1–21.