Abstract
Lagrangian (following the flow) motion of fluid particles is connected to the Eulerian (spatio-temporally
defined) velocity field by an ordinary differential equation. Lagrangian motion causes transport (of fluid,
as well as of other quantities such as energy and pollutants), and is therefore highly relevant to climate,
weather, the environment, industrial mixers, etc. Velocities,
on the other hand, are often known in an Eulerian sense from observational, experimental or computational
fluid dynamics data. The uncertainties in this Eulerian data will lead to uncertainties in Lagrangian prediction.
This short article examines some recent approaches to ascribing uncertainties to Lagrangian-derived
information.
1. Lagrangian motion and uncertainty
At all spatial scales (microfluidic, laboratory, oceanic, atmospheric, Jovian), transport is governed
by where particles move, based on the forces acting on them. This movement strongly influences how chemical
species, pollutants, nutrients, heat energy, and biological organisms disperse. Forces give rise to a velocity
with the relationship given, for example, by the Navier-Stokes equations in the classical fluids context. Velocities
are typically unsteady (time-varying), and the rate of change of a particle's position is then instantaneously given by
the local velocity. The tracks followed by particles as a result of these velocities is called Lagrangian motion
in fluid mechanics. In many cases, we have the velocities from observations, experiments, computational
fluid dynamics simulations, or via some proxy (e.g., inferring Jupiter's atmospheric motion by performing optical
processing on videos, imputing oceanic velocities from satellite measurements of sea-surface heights by using
the so-called geostrophic approximation). These are Eulerian measurements, meaning that they are
given as functions of position and time. Since the information is moreover typically available on a spatial grid,
at discrete values of time, there are uncertainties in the velocity away from gridpoints. In addition to this
resolution uncertainty, there are uncertainties even at the gridpoints, due to measurement error, modelling (proxy) limitations,
and intermittent coverage (e.g., spatio-temporally varying cloud cover when taking satellite measurements). Given this, computed
Lagrangian trajectories will certainly inherit these errors.
While earlier work in fluid mechanics typically used Eulerian measurements to (incorrectly) infer transport, seminal work
by Haller and Yuan [16] highlighted the necessity of understanding Lagrangian trajectories instead. One aspect is to try to find coherent entities
which, while moving in time as is inevitable in an unsteady velocity field, retain aspects which might be deemed
`coherent.' Some examples of this could be the Gulf stream, an oceanic eddy, a moving hurricane, the Antarctic
circumpolar vortex (`ozone hole'), Jupiter's Great Red Spot, a blob of fluid in a mixing device which resists mixing
with its surroundings, and so on. There is considerable debate as to what `coherent' might mean, as is evidenced by
the multitude of
diagnostic methods which have—and continue to be—developed in this endeavor. In general nomenclature,
these have come to be known as `Lagrangian Coherent Structure' (LCS) methods, for which several recent reviews are available
[8, 13, 22, 24] (note that Haller uses the acronym LCS only for a specific class of methods
[14]; more generality will be assumed here). As a general rule, coherent structures are detected from LCS methods
based on numerical computations from the
velocity field, treated as if it had no uncertainty whatsoever. The spatial resolution issue is dealt with
in the natural way by interpolation between gridpoints, thereby implicitly using a deterministic and (sufficiently)
smooth velocity across the entire domain.
Recently, the issues of uncertainty is coming to the fore. In climate and weather prediction, unresolved
subgrid-scale uncertainties are sometimes referred to as the `stochastic parametrization problem' (typically
in these situations, the emphasis is not on tracking coherent structures but in determining the impact of
this lack of resolution on global properties such as atmospheric carbon-dioxide). From the Lagrangian
perspective, there have been some recent studies which examine the impact of uncertainty in some specific
ways [3, 7, 11]. Moreover, there has been an increased interest in
understanding diffusion in the context of fluid flows [15, 18]; this is related because diffusion can be thought of as
the macroscopic impact of trajectory uncertainties at the microscopic level.
For an unsteady flow in \({\mathbb R}^n\) (typically, \(n = 2\) or \(3\)) from time \(0\) to \(T\), purely deterministic trajectories solve
|
\( \frac{\mathrm{d} \vec{x}}{\mathrm{d} t} = \vec{u}(\vec{x},t) \, , \quad \mbox{also written as} \quad \mathrm{d} \vec{x}_t = \vec{u}(\vec{x}_t,t) \, \mathrm{d} t \, . \) |
(1) |
In determining the trajectories \(\vec{x}_t\) numerically, since \(\vec{u}\) is only available on a spatio-temporal grid, it is
necessary to interpolate \(\vec{u}\); this artificially extends \(u\) beyond the gridpoints, while imposing smoothness.
A model for including uncertainties in the velocity field would be to consider instead
|
\( \mathrm{d} \vec{y}_t = \vec{u}(\vec{y}_t,t) \, \mathrm{d} t + \varepsilon \, \mathrm{d} \vec{W}_t \) |
(2) |
in which \(y_t\) is the random trajectory, where now \(\varepsilon \, \mathrm{d} \vec{W}_t\) gives the uncertainty in location over
an infinitesimal time \(\mathrm{d} t\), due to \(n\)-dimensional Brownian motion \(\vec{W}_t\). As necessary in considering
the noise as a perturbative effect on (1), it is required that \(\left| \varepsilon \right|\) be small in some appropriate sense. A more general formulation
which allows for the uncertainties to be spatio-temporally dependent would be to replace \(\varepsilon \, \mathrm{d} \vec{W}_t\) with
\(\varepsilon \sigma(\vec{y}_t,t) \mathrm{d} \vec{W}_t\), where \(\sigma\) is an \(n \times n\) matrix [3, 7];
its \((\vec{y}_t,t)\)-dependence allows, for example, modeling of cloud cover, or certainty at gridpoints. Given a stochastic
model such as (2), along with spatio-temporal data \(\vec{u}\) on a grid, it is a straightforward exercise to evolve
trajectories \(\vec{y}_t\) using, for example, the Euler-Maruyama scheme [19]. For each initial condition, this
generates one realization of the stochasticity. By taking all initial conditions (consisting of a finite number of discrete points in any numerical situation),
one can then pretend that this generates a flow map from time
\(0\) to \(T\), which is the basic building block used in most LCS methods. Thus, various LCS techniques, such
as the Finite-Time Lyapunov Exponent (FTLE) can easily be applied to this map [12, 9]. Typically,
a function \(f\) on the space of initial conditions \(\vec{x}_0\) is generated by an LCS method, and based on its structure,
coherent structures are identified. In the case of the FTLE, \(f\) would be the stretching field: if an
infinitesimal \(n\)-dimensional sphere were positioned at the initial location \(\vec{x}_0\), it would get mapped to an ellipsoid under the flow, and
\(f\) is the ratio of the ellipsoid's longest axis to the original diameter. Indeed in this case, \(f\) is equal to the
operator norm of the gradient of the flow map, thereby illustrating that \(f\) would be in general derived
from the flow map in some way. \((n-1)\)-dimensional ridges of this stretching field (usually scaled using
a logarithm as befitting the `exponent' part of the definition, and divided by the time-of-flow \(T\)) are important flow
separators between coherent structures. Other LCS methods similarly use different diagnostic fields \(f\),
and interpret them in diverse ways such as using zero contours of \(f\) to define a region (see [8] for more detailed descriptions of many methods).
In doing this for one implementation of (2), it should be noted that
the map is random (more properly, there is a different random dynamical system that each initial condition is subject
to). Thus computing an \(f\) for one such map cannot give us uncertainty information in general. The solution,
then, would be to perform this numerous times—thereby using many realizations of the stochasticity—and compute
\(f\) for each such instance. A function \(\bar{f}\), which is \(f\) averaged over all implementations, can then
be developed. When using \(\bar{f}\) as a diagnostic method instead of \(f\), alternative sets associated
with LCSs are obtained: more diffused ridges of the FTLE field [12, 6] or of the
stable/unstable manifolds [7].
The question, though, is how an uncertainty can be assigned to any features obtained this way. As a computational
tool, the standard deviation (with respect to implementations) of \(f\) is the obvious way to proceed. Any such spatially-varying
uncertainty field shall be denoted \(s\). At each initial location \(\vec{x}_0\), \(s(\vec{x}_0)\) will provide an indicative value of the uncertainty of \(\bar{f}(\vec{x}_0)\). Put
another way, the true diagnostic field \(f(\vec{x}_0)\) could then be considered as
|
\( f(\vec{x}_0) = \bar{f}(\vec{x}_0) \pm s(\vec{x}_0) \quad {\mathrm{or}} \quad
f(\vec{x}_0) \in \left[ \bar{f}(\vec{x}_0) - s(\vec{x}_0), \bar{f}(\vec{x}_0) + s(\vec{x}_0) \right] \, . \) |
(3) |
While this process will numerically generate heuristic, spatially-dependent errors, it remains unsatisfying
and computationally expensive. Is it possible
to find errors in theoretical ways? A recent article [5] does just this, though it is confined to
two-dimensional flows. For the computation of any diagnostic
function \(f\), the flow map from time \(0\) to \(T\) is required. This flow map acquires an uncertainty, and determining
a quantifier for this is therefore necessary for almost any LCS method one considers. For the more general
spatio-temporally modulated noise \(\varepsilon \sigma(\vec{y}_t,t) \, \mathrm{d} \vec{W}_t\) situation, the theory in [5] concerns the
scaled deviation of the random trajectory of (2) from that of (1) by time \(T\), if they began
at the same initial condition \(\vec{x}_0\), i.e.,
|
\( \vec{z}(\vec{x}_0) := \frac{\vec{y}_T - \vec{x}_T}{\varepsilon} \quad {\mathrm{where}} \quad \vec{x}_0 = \vec{y}_0 \, . \) |
(4) |
Under some natural conditions, the expectation \(\mathbb{E} \left( \vec{z}(\vec{x}_0) \right)\) is shown to be zero
for all relevant \(\vec{x}_0\) in the limit as \(\varepsilon \rightarrow
0\) [5], indicating that to leading-order the expected final location due to (2) is precisely the deterministic one. Here, the expectation is with respect to many realizations of (2). A theoretical expression for the uncertainty is suggested as the stochastic sensitivity field \(S^2\) given by [5]
|
\( S^2(\vec{x}_0) := \lim_{\varepsilon \rightarrow 0} \sup_\theta {\mathbb{V}} \left( \vec{z}(\vec{x}_0) \cdot \hat{\vec{n}}_\theta \right) \, , \) |
(5) |
where \({\mathbb{V}}\) is the variance with respect to many realizations of (2), and \(\hat{\vec{n}}_\theta =
\left( \cos \theta \, \, \, \sin \theta \right)^\top\) is a unit vector pointing in the direction \(\theta
\in [-\pi/2, \pi/2)\). Thus, the variance
is computed for a projected direction \(\theta\), and then optimized over all directions in defining \(S^2\). This
can be expressed in terms of trajectories and the evolution of the flow map gradient of the
deterministic system (1), as well as the matrix \(\sigma\), as is given in Appendix A. Once \(S^2\) is known, it is noted that the lengthscale of the uncertainty is to leading-order
|
\( L(\vec{x}_0) = \varepsilon \sqrt{S^2(\vec{x}_0)} \, . \) |
(6) |
If \(\mathrm{d} \vec{W}_t\) is thought of in the standard sense (\(\mathrm{d} \vec{W}_t \sim \sqrt{\mathrm{d} t}\) is assumed in terms of units) in (2), a natural way to
think of \(\varepsilon\) is that
|
\( \varepsilon = \sqrt{ v_r \, L_r} \) |
(7) |
where \(L_r\) and \(v_r\) are resolution scales of the length (i.e., the grid spacing) and velocity respectively.
If so, (6) gives \(L(\vec{x}_0) = \sqrt{ v_r \, L_r \, S^2(\vec{x}_0)}\), enabling the computation of the
lengthscale uncertainty of a trajectory beginning at \(\vec{x}_0\), in terms of scales in the problem. (In this interpretation,
\(S^2\) has dimensions of time, and \(\vec{z}\) in (4) has dimensions of \({\mathrm{time}}^{-1/2}\).)
Having this information on Lagrangian uncertainty can be invaluable in evaluating uncertainties of various
LCS measures. The implications on the FTLE field have already been determined, and are
presented in another recent article [6]. There are unexpected implications—validated via stochastic
simulations as well—that the error increases when the spatial resolution scale \(L_r\) is decreased.
An alternative—but connected—way to understand uncertainties in Lagrangian locations is through diffusion.
Equations such as (2) are directly connected to a Fokker-Planck equation, which governs how densities
\(\rho(x,t)\) of trajectories evolve. The Brownian term in (2) contributes to a standard diffusive (Laplacian)
operator on \(\rho\) (it is more complicated if the spatio-temporally modulated form is considered instead). The
Fokker-Planck equation is thus an advection-diffusion equation for \(\rho\); indeed, precisely so if \(u\) is area-preserving
and \(\sigma\) is the identity. Because of this reason or otherwise, there have been several recent approaches
[10, 15, 18] to uncertainty by using an advection-diffusion partial differential equation
rather than the stochastic ordinary differential equation (2).
Thus theoretical methods for uncertainty quantification of Lagrangian trajectories is an
emerging area, for which the methods mentioned here may set the foundations.
In the remainder of this article, some examples related to the discussion above will be presented.
2. Rotating Arc Mixer (RAM)
Figure 1. Movie of the RAM flow protocol: (left) deterministic, and (right) stochastic with \(\varepsilon = 0.01\).
The Rotating Arc Mixer (RAM) is an industrial mixing device, where a thin layer of fluid lies within a circle, and
where the circular boundary has a number of arcs which can be moved independently [20, 21]. A three-dimensional
version of this has the circle extending outwards to form a cylinder, with the additional potential for the arcs to be in
different locations as one proceeds in the extended direction [23]. The idea is to use arc-activation in a
sequential way, usually to maximize mixing of fluid within the device. For the standard two-dimensional situation
with \(\vec{x} = (x,y)\),
there is an exact expression for the velocity within the circular domain \(x^2 + y^2 < 1\)
when one arc is activated
[17]. This explicit solution is derived using Stokes' flow with appropriate boundary conditions, and repeated activation
of different arcs can then be implemented by combining these solutions appropriately. The presence of chaotic
mixing in the RAM is well-established [21, 23]; hence it is an example in which laminar flow can be used to generate
significant mixing.
Figure 2. Final \(x\)-locations of the particles, shown as a contour plot over initial conditions: (left) deterministic, and
(right) average of the stochastic simulations with \(N = 500\).
In the simulations to be presented, a RAM with three equally spaced arcs, each subtending an angle of \(\pi/4\), is considered. The flow protocol is that the arcs are rotated counterclockwise
with angular frequency \(1\), in a sequential order which also goes counterclockwise. Four complete cycles of this
procedure constitutes the full flow. Fig. 1 shows on the left the impact of the flow on a scatter of fluid particles. The right figure, for comparison, is a stochastically perturbed version; specifically one realization of
(2) with \(\varepsilon = 0.01\). Here, each particle is subject to its own stochastic perturbation (as opposed
to the approach in `random dynamical systems,' in which one stochastic realization affects all initial conditions
equally [2]). The arcs are shown in blue, and when active, an arc is highlighted in magenta. Particles near an
active arc get pushed strongly counterclockwise, and the impacts are also felt (more weakly) in the interior regions.
The RAM system is highly sensitive at the points where the moving arcs end, because
the angular velocity on the circular boundary at these points has a discontinuity during the instances in which that
particular arc is active: it is \(1\) on the arc, and \(0\) just off it. While there is no discontinuity within the fluid
region, the intense gradients near to these points leads to high uncertainty. The deterministic simulation gives evidence
of this; particles apparently `disappear' because in the advection they are sometimes transported across the flow
boundary. The stochastic simulation, in which additional velocity uncertainties are imposed on top of the high-gradient
issue, is understandably affected substantially more by this issue than the deterministic case.
Figure 3. Uncertainty in final location of the RAM flow, expressed as a lengthscale: (left) using stochastic sensitivity theory (5), and (right) using a distance measure from the deterministic location, based on \(N = 500\) stochastic simulations.
The expectation is that if the final location is averaged over many stochastic simulations, the deterministic final
location is obtained. This is investigated in Fig. 2, where the final \(x\)-locations are averaged
from the \(N = 500\) simulations, and illustrated here as a contour plot as a function of the initial locations \(\vec{x}_0\).
Thus, the value of \(x_T\) is the indicative \(f\) function considered in this instance.
There is clear spatial correlation with \(x_T\) of the deterministic advection (shown on the left).
The blank areas in this figure correspond to particles for which information has been lost due to the effect described
in the previous paragraph. This effect can be reduced by choosing more refined grids locally and—in the case of
the stochastic simulations—by performing many more stochastic realizations.
The stochastic sensitivity (5) is a measure of the uncertainty in the final location. More specifically,
it can be converted to a lengthscale according to (6). By computing \(S^2\) using the expressions
and methodology outlined in [5] (and summarized in Appendix A), the lengthscale of the uncertainty in the final location can be computed. This is displayed in the left
figure of Fig. 3. Strong uncertainty `ridges' are seen to emanate from the ends of the arcs,
which is expected because of prior discussions. There is a nonuniform distribution of the uncertainty (which is
related to the presence of chaotic regions and islands in the flow; however, these concepts only make sense in
reality for infinitely many iterations of the flow). The figure on the right shows a proxy for the uncertainty in
final location as the distance from the deterministic to the averaged stochastic final location, based on the \(500\)
stochastic simulations. While the two concepts are not identical (compare with the definition (5)), it is
clear that the spatial structures of the uncertainties are somewhat correlated. The two figures in Fig. 3 thus provide a theoretical and a numerical estimate for the Lagrangian uncertainty,
as fields \(s(\vec{x}_0)\).
This example has been chosen to highlight how stiff velocity gradients can cause considerable difficulties, as evidenced here by `white'
regions in the plots where information is unavailable. In this case, this generates uncertainties such that even
deterministic advection—with an explicitly known velocity field—leads to problems. Basically, numerical uncertainties
in integration are enough to generate this issue. When performing stochastic simulations, in which the velocity
uncertainties are more explicitly included in the numerics, computations are affected in even larger regions. Most
examples will not have as extreme a velocity gradient as this, and in such cases computations can often be
performed without such clear evidence of them being wrong. Nevertheless, they will carry some errors. The
stochastic approaches described in this article can help access information on these.
3. Oceanic velocities
Figure 4. Movie of the evolution of the FTLE field (in units of \({\mathrm{day}}^{-1}\)) in the northern Atlantic. A time-of-flow of \(T = 20\) days is used throughout, and the FTLE field is at \(t_0\) (the number of days from 1 June 2008) in each frame of the movie.
Next, oceanic velocity data (Unidata Dataset Discovery V1.0) from the European Union's Copernicus Marine Environment Monitoring Service
is briefly analyzed. This dataset has zonal and meridional geostrophic velocities (on the ocean surface), derived
from multiple satellite measurements of sea-surface height of the global ocean. In this case, daily data
from 1 June 2008 for a 3380 day period is available for the region in the northern Atlantic bounded by longitudes
\(77.125^0\)W to \(58.875^0\)W and
latitudes \(30.125^0\)N to \(50.125^0\)N. (Note: \(x\) will be represented as longitude eastwards, and \(y\) as
latitude northwards; time will be measured in days and distances in degrees.) Gridpoints are separated by 0.25 degrees in both directions, and hence
the resolution lengthscale is \(L_r = 0.25 \, {\mathrm{deg}}\).
A particularly significant issue in this dataset is to be noted. While oceanic water is incompressible, this is with
respect to three-dimensional flow. Thus, the horizontal velocity components \(u\) (eastward) and \(v\)
(northward) of ocean surface velocities need not
necessarily obey (two-dimensional) incompressibility. In other words, it is not necessarily true that \(u_x + v_y = 0\),
and indeed this data set indicates strong violations of two-dimensional incompressibility, particularly so when close to
shorelines where strong downwelling may occur.
In this example, rather than viewing the final locations, a particular field associated with these locations will be
examined. This is the Finite-Time Lyapunov Exponent (FTLE) field, which essentially measures the stretching of
infinitesimal fluid parcels over a given time duration. More explicitly, this is given by
|
\( {\mathrm{FTLE}}(\vec{x}_0) = \frac{1}{T} \ln \tau(\vec{x}_0) \) |
(8) |
where the stretching field \(\tau\) intuitively positions a tiny disk (say diameter \(\delta\)) at \(\vec{x}_0\), allows it
to evolve with the deterministic flow (1) till time \(T\) to form an ellipse, and then takes the ratio of the longest
extended axis of this ellipse with \(\delta\). Computationally, \(\tau(\vec{x}_0)\) is exactly the operator norm of the
gradient of the flow map from time \(0\) to \(T\), and thus after advecting particles using (1), numerical
differentiations for estimating flow map derivatives are necessary. The logarithm and then division by the time-of-flow
in (8) converts the FTLE to an exponential rate value, in units of inverse time. It is common in oceanography
to use the Finite-Size Lyapunov Exponent (FSLE) instead of the FTLE; this uses the same formula as (8), but
stops the calculation once the stretching has gone beyond a specified threshold and then uses the \(T\) at that
instance in the formula. The FTLE is preferred from the dynamical systems perspective, because it ensures that all
initial conditions are subject to the identical dynamical system (namely, (1) flowing from time \(0\) to \(T\)),
which is not true for an FSLE computation. In any case, FTLEs/FSLEs are probably the most commonly used diagnostic
technique for evaluating regions of coherence.
Figure 5. FTLE computed from two stochastic simulations of the ocean data from 11 November 2013 for \(20\) days using
\(v_r = 0.02 \, {\mathrm{deg}} / {\mathrm{day}}\).
Thus, the FTLE on 1 June 2008 (the first day of data availability), if considering the flow for time \(T = 20\) days into
the future, will give a field. This will be indexed by stating that the initial day is \(t_0 = 0\). Similarly, it is possible to
compute an FTLE field on 2 June 2008, again by considering a flow over a time of \(T = 20\) days, but this will then
be indexed by \(t_0 = 1\). In general, the time \(t_0\) will be taken to be the number of days after 1 June 2008.
With this understanding, a movie of the evolution of the FTLE field as \(t_0\) changes (but the time-of-flow remains at \(20\)
in all cases), is shown in Fig. 4. A two-year period is shown in the movie, and the white region
is of course the northern American continental landmass. Large FTLE values indicate regions which will stretch considerably,
while low FTLE value regions will not. Hence, the usual argument goes, the former regions would be associated with high levels of flow/shear/mixing
(e.g., the flanges of the Gulf Stream), while
the latter regions have lower values (e.g., coherent eddies). It should be noted that these conclusions must be
considered relative to the flow duration (\(20\) days into the future) considered; choosing a different time-of-flow
will give information on that particular time-window. (FTLE/FSLE fields such as these are used often to reach
various conclusions about coherent regions of flows; `ridges' have received significant focus in analogy to
stable/unstable manifolds of infinite-time flows. However, stable and unstable manifolds cannot be theoretically
defined in general for finite-time flows such as this.)
Figure 6. FTLE computed from the ocean data from 11 November 2013 for \(20\) days using \(v_r = 0.02 \, {\mathrm{deg}} / {\mathrm{day}}\): (left) deterministic, and (right) average of \(200\) stochastic simulations.
The velocities which are used as input here have significant uncertainties: spatial resolution, measurement error,
modelling error (the geostrophic approximation is often used, and other processing also occurs). Thus, any computation
of eventual fluid parcel location will have errors. Put another way, there are errors in the flow map derived from the
Eulerian velocities. Since this flow map is often used in computing various diagnostics \(f\) of the flow (case in point: the
FTLE field), it is instructive to try to put some error-bars on any fields \(f\) computed. For \(f =\) FTLE, there is
recent work [6] which does exactly this, generating an error field \(s(\vec{x}_0)\) which depends on
estimates of the velocity uncertainty scale \(v_r\) and the spatial resolution scale \(L_r\) according to [6]
|
\( s(\vec{x}_0) = \frac{2 \sqrt{v_r \, S^2(\vec{x}_0)} }{e^{{\mathrm{FTLE}}(\vec{x}_0) \left| T \right|}
\sqrt{L_r} \, \left| T \right| } \, . \) |
(9) |
\(\mathrm{FTLE}(\vec{x}_0)\) is the FTLE field computed from (8) assuming that no uncertainties
are present.
The stochastic sensitivity field \(S^2\) enters into this because the uncertainty in the final location is what leads to
uncertainties in the FTLE field. In the general sense, it is noted that \(s \sim \sqrt{v_r / L_r}\) in terms of the two
uncertainty scales.
Figure 7. Error in FTLE computed from the ocean data from 11 November 2013 for \(20\) days using \(v_r = 0.02 \, {\mathrm{deg}} / {\mathrm{day}}\): (left) theoretical bound using (9), and (right) standard deviation of FTLE from \(200\) stochastic simulations.
To examine the impact of uncertainty, let us focus on the FTLE field on 21 November 2013, using once again a
\(20\)-day period into the future. The resolution lengthscale \(L_r = 0.25\) deg for this data set, and a velocity uncertainty
of \(v_r = 0.02 \, {\mathrm{deg}} / {\mathrm{day}}\) will be assumed; this is very
small considering that the velocities are typically of the order \(1 \, {\mathrm{deg}} / {\mathrm{day}}\). Thus from
(7), \(\varepsilon = 0.0707 \, {\mathrm{deg}} \, {\mathrm{day}}^{-1/2}\). Now, (2) can be implemented
for this duration of time, and pretending that the final flow map is accurate, the FTLE field can be computed from
(8). Fig. 5 displays the result of two such simulations. Some differences are to be
observed, as expected. Thus, even with about \(2\%\) uncertainty in the velocity field, there can be errors
in computed FTLEs. Next, following recent ideas that averaging stochastic simulations give a better approximation of
the truth [12, 9, 6], the averaged FTLE field for \(200\) such simulations is shown on the right
in Fig. 6. The left figure is of the FTLE field based on the deterministic data,
pretending that there are no uncertainties. On a gross level, the figures are very similar, indicating a level of robustness
of the deterministic FTLE prediction to this level of noise. However, there are differences in the two figures in
particular with respect to finer details. Clearly, the high-resolution detail from the deterministic FTLE field is
unreasonable to use when uncertainties are considered. The averaged version of many stochastic realizations, on
the other hand, smoothens out these finer (and unjustifiable) details.
Fig. 7 displays two different ways of attempting to capture the spatially-varying errors in the
deterministically computed FTLE field. The left figure is computed using (9). Note its scale in comparison
to Fig. 6. The error values are large, and a prime factor in this is that the data turns out
to not obey area-preservation, presumably due to upwelling and downwelling. When computing \(S^2\) (see equation
(10) in Appendix A) for the
usage of (9), the term \(\vec{\nabla} \cdot \vec{u}\) turns out to contribute significantly. The implications from the theory are that the deterministic FTLE field is therefore
highly uncertain. However, it should be noted that the computation of this FTLE error field is difficult in this instance
because of the irregular boundary (the coastline), which in particular renders the \(S^2\)-field computation
difficult. (Several levels of interpolation is needed; the errors that can arise in nonregularly scattered data points
is well-known.) The right figure is obtained from the \(200\) stochastic simulations.
The standard deviation (across the \(200\) values) is computed for the FTLE value at each point to generate this
field. It should be noted that in performing stochastic simulations here, the simple strategy of excising particles when
they cross to land has been done. If using sufficiently many stochastic simulations, there will be sufficient relevant
data left, in order to do this calculations effectively. However, boundaries are a vexing problem when using
stochastic simulations, since it is not always clear what the physically appropriate boundary conditions should be.
Another difficulty with using the standard deviation of stochastic simulations is the computational cost: many
simulations need to be performed. Theoretical estimates are therefore highly desirable. In
any case, the yellow regions (lower theoretical errors) in the left figure do appear to be in regions where the
standard derivations (right figure) are small. While both approaches (left and right) give methods for estimating
errors in the FTLE field, there is clearly scope for improvement of these ideas because of issues such as boundary
conditions, as well as improving computational efficiencies.
It should be noted that a conservative error scale for the velocity uncertainty was used here: \(v_r = 0.02
\, {\mathrm{deg}} / {\mathrm{day}}\)
corresponds to an error of about \(2 \%\) in the velocity measurements. The true errors are likely to be larger,
which means that understanding the uncertainties in eventual predictions is crucial.
4. Concluding remarks
This short article has highlighted the emerging issue that uncertainties in velocities are relevant when
inferring Lagrangian motion. Of the two examples shown, the first was mainly concerned about the uncertainty in
the final location, whereas the second focussed on a quantity derived from the location. Uncertainties play a
role in either situation. Moreover, the first
was related to a patented [20] industrial device, whereas the second was from observational data. It is for such
instances—not merely in toy models—that knowing the uncertainties has a strong practical value.
In general, the main point made in this article is that there needs to be more ways of ascribing uncertainties to any
predictions made using spatio-temporally available data. One approach—using stochasticity to model the uncertainty in the context of
the governing ordinary differential equation—is the prime focus here, and
in this case, glimpses of both numerical and theoretical methods have been provided. A closely related
approach is to use the Fokker-Planck equation [10, 15, 4] corresponding to the stochastic
ordinary differential equation
(2); this encapsulates the macroscopic impact of the microscopic Brownian motions, though not per se
capturing uncertainties in Lagrangian trajectories. A continuing development of uncertainty quantification in
this realm—building possibly on the methods discussed here or those from different viewpoints [1, 11]—is essential in interpreting Lagrangian results
correctly.
Acknowledgements
Partial support from the Australian Research Council via grant DP200101764 is gratefully acknowledged. I also express my thanks to Michel Speetjens for allowing me to mine his expertise in the RAM device, and to Irina Rypina for providing me with the Copernicus data.
A. Stochastic sensitivity
This section gives the formula for the stochastic sensitivity field associated with the stochastic
differential equation (2), as developed in [5]. The formula to be stated is only
valid for two-dimensional situations. Extensions to three dimensions are in progress.
Let \(\vec{F}_\xi^\eta\) be the flow map of the deterministic flow from time \(\xi\) to \(\eta\). Thus,
\(\vec{x}_T = \vec{F}_0^T(\vec{x}_0)\) is the location to which the particle beginning at \(\vec{x}_0\) at time \(0\)
gets mapped by the flow till time \(T\). For \(t \in [0,T]\), define the \(2 \times 2\) matrix function
|
\( H(\vec{x}_T,t) := \exp \left[ \int_t^{T} \! \! \left[ \vec{\nabla} \cdot \vec{u} \right]
\left( \vec{F}_{T}^\xi(\vec{x}_T),\xi \right) \, \mathrm{d} \xi
\right] \, \left( \begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array} \right) \, \vec{\nabla} \vec{F}_{T}^t(\vec{x}_T) \, , \) |
(10) |
where the gradient in the final term is with respect to locations \(\vec{x}_T\). In terms of the components
of \(H\), let
|
\( \begin{eqnarray*}
P(\vec{x}_T) & = & \left| \frac{1}{2} \sum_{i=1}^2 \sum_{j = 1}^2 \int_{0}^{T} H_{ij}^2(\vec{x}_T,t) \, \mathrm{d} t \, \right| \, , \\
L(\vec{x}_T) & = & \frac{1}{2} \int_{0}^{T} \left[ \sum_{i=1}^2 H_{i2}^2(\vec{x}_T,t) - \sum_{i=1}^2 H_{i1}^2(\vec{x}_T,t) \right] \mathrm{d} t \, , \\
M(\vec{x}_T) & = & \int_{0}^{T} \! \! \sum_{i=1}^2 \left[ H_{i1}(\vec{x}_T,t) \, H_{i2}(\vec{x}_T,t) \right] \! \mathrm{d} t \, ,\quad {\mathrm{and}} \\
N(\vec{x}_T) & = & \sqrt{L^2(\vec{x}_T) + M^2(\vec{x}_T)} \, .
\end{eqnarray*} \) |
|
Each of these quantities can be computed quite easily for a given final location \(\vec{x}_T\). The complication
is that the flow map gradient needs to be evaluated along backwards trajectories starting at \(\vec{x}_T\), until
they reach \(\vec{x}_0\). At each instance in time, the gradient can be evaluated in the standard way (e.g., as is
normally done in finite-time Lyapunov exponent computations). Now, the stochastic sensitivity—as defined
in (5)—at \(\vec{x}_0\) is then given by [5]
|
\( S^2(\vec{x}_0) = P(\vec{x}_T) + N(\vec{x}_T) \, . \) |
(11) |
A field across all \(\vec{x}_0\) can thus be constructed numerically. The proof of this result is based in stochastic calculus
methods (It\^{o} lemma and isometry), coupled with some techniques more familiar to classical ordinary differential
equations [5]. Moreover, a more general version associated with the spatio-temporally varying
stochasticity \(\varepsilon \sigma(\vec{y}_t,t) \mathrm{d} \vec{W}_t\) is also available [5]. Numerically applying this technique does
have its difficulties, because one needs to map the values of \(S^2\) (which is computed at \(\vec{x}_T\) by
this process) to \(\vec{x}_0 = \vec{F}_T^0(\vec{x}_T)\) to generate a field at the time instance \(0\). If a
regular grid were chosen for the distribution of \(\vec{x}_T\), since the \(\vec{x}_0\)s will not be regularly
distributed, it will be necessary to extend \(S^2\) from these scattered points to the
relevant domain. Depending on the scatter, this process can be error-prone (see [5] for more details).
References
[1] F. Antown, D. Dagičević, and G. Froyland, Optimal linear responses for Markov chains and stochastically
perturbed dynamical systems, J. Statis. Phys., 170 (2018), pp. 1051-1087.
[2] L. Arnold, Random Dynamical Systems, Springer-Verlag, Berlin, 1998.
[3] S. Balasuriya, Stochastic uncertainty of advected curves in finite-time unsteady
flows, Phys. Rev. E, 95 (2017), 062201.
[4] S. Balasuriya, Stochastic approaches to Lagrangian coherent structures, Adv. Stud. Pure. Math., accepted, 2019.
[5] S. Balasuriya, Stochastic sensitivity: a computable uncertainty measure for unsteady
flows, SIAM Review, accepted, 2020.
[6] S. Balasuriya, Uncertainty in finite-time Lyapunov exponent calculations, J. Comp. Dyn., under review, 2020.
[7] S. Balasuriya and G. Gottwald, Estimating stable and unstable sets and their role as transport
barriers in stochastic flows, Phys. Rev. E, 98 (2018), 013106.
[8] S. Balasuriya, N. Ouellette, and I. Rypina, Generalized Lagrangian coherent structures, Physica D, 372 (2018), pp. 31-51.
[9] F. Balibrea-Iniesta, C. Lopesino, S. Wiggins, and A. Mancho, Lagrangian descriptors for stochastic differential equations: a tool
for revealing the phase portrait of stochastic dynamical systems, Intern. J. Bifurc. Chaos, 26 (2016), 1630036.
[10] A. Denner, O. Junge, and D. Matthes, Computing coherent sets using the Fokker-Planck equation, J. Comp. Dyn., 3 (2016), pp. 163-177.
[11] G. Froyland, C. González-Tokman, and A. Quas, Stochastic stability of Lyapunov exponents and Oseledets
splitting for semi-invertible matrix cocycles, Commun. Pure Appl. Math., 68 (2015), pp. 2052-2081.
[12] H. Guo, W. He, T. Peterka, H.-W. Shen, S. Collis, and J. Helmus, Finite-time Lyapunov exponents and Lagrangian coherent structures
in uncertain unsteady flows, IEEE Trans. Visual. Comp. Graphics, 22 (2016), pp. 1672-1682.
[13] A. Hadjighasem, M. Farazmand, D. Blazevski, G. Froyland, and G. Haller, A critical comparison of Lagrangian methods for coherent structure
detection, Chaos, 27 (2017), 053104.
[14] G. Haller, Lagrangian Coherent Structures, Annu. Rev. Fluid Mech., 47 (2015). pp. 137-162.
[15] G. Haller, D. Karrasch, and F. Kogelbauer, Material barriers to diffusive and stochastic transport, Proc. Nat. Acad. Sci., 115/37 (2018), pp. 9074-9079.
[16] G. Haller and G.-C. Yuan, Lagrangian coherent structures and mixing in two-dimensional turbulence, Phys. D, 147 (2000), pp. 352-370.
[17] T.-Y. Hwu, D.-L. Young, and Y.-Y. Chen, Chaotic advection for Stokes flow in circular cavity, J. Engin. Mech., 123 (1997), pp. 774-782.
[18] D. Karrasch and J. Keller, A geometric heat-flow theory of Lagrangian coherent structures, arXiv:1608.05598, under review, 2018.
[19] P. Kloeden, E. Platen, and H. Schurz, Numerical solution of SDE through computer experiments, Universitext, Springer-Verlag, Berlin, 1994.
[20] G. Metcalfe and M. Rudman, Fluid mixer using viscous drag, U.S. Patent, 7,121,714, 2006.
[21] G. Metcalfe, M. Rudman, A. Bryden, L. Graham, and R. Hamilton, Composing chaos: An experimental and numerical study of an open duct
mixing flow, AIChE J., 52 (2006), pp. 9-28.
[22] S. Shadden, Lagrangian coherent structures, In R. Grigoriev, editor, Transport and mixing in laminar flows:
from microfluidics to oceanic currents, Wiley, 2011.
[23] M. Speetjens, E. Demissie, G. Metcalfe, and H. Clercx, Lagrangian transport characteristics of a class of three-dimensional
inline mixing flows with fluid inertia, Phys. Fluids, 26 (2014), 11360.
[24] M. Speetjens, G. Metcalfe, and M. Rudman, Lagrangian transport and chaotic advection in three-dimensional
laminar flows, arXiv:1904.07580, under review, 2019.