A bridge between invariant dynamical structures and uncertainty quantification

We develop a new quantifier for forward time uncertainty for trajectories that are solutions of models generated from data sets. Our uncertainty quantifier is defined on the phase space in which the trajectories evolve and we show that it has a rich structure that is directly related to phase space structures from dynamical systems theory, such as hyperbolic trajectories and their stable and unstable manifolds. We apply our approach to an ocean data set, as well as standard benchmark models from deterministic dynamical systems theory. A significant application of our results, is that they allow a quantitative comparison of the transport performance described from different ocean data sets. This is particularly interesting nowadays when a wide variety of sources are available, since our methodology provides avenues for assessing the effective use of these data sets in a variety of situations.


Introduction
Uncertainty Quantification (UQ) searches for a quantitative characterization of uncertainties in computational models of real-world applications. Uncertainties associated to these models have been studied for many years. Depending on the context they can arise as a result of noisy or incomplete data, errors introduced by numerical models, or by the lack of complete understanding of the governing physical processes.
In geophysical contexts, uncertainty is a topic of much interest because it is inherent to the equations describing the motion of fluids such as the ocean or the atmosphere. The motivation of this paper is trying to acquire a deeper understanding of uncertainty quantification, that allows a better characterization of its presence in ocean models. Large amount of oceanic data are becoming now available. For instance, the Copernicus Marine Environment Monitoring Service (CMEMS) provides regular and systematic information about the physical state and dynamics of the ocean for the global ocean and the European regional seas. The data cover the current and future state of variables and the provision of retrospective data records (re-analysis). Other oceanic services include ocean currents supplied by altimeter satellites, like AVISO; the HYbrid Coordinate Ocean Model (HYCOM), a consortium, which is a multi-institutional effort sponsored by the National Ocean Partnership Program (NOPP), as part of the U. S. Global Ocean Data Assimilation Experiment (GODAE), to develop and evaluate a data-assimilative hybrid coordinate ocean models, etc.
Ocean models are built on partial differential equations, which are only an approximation of reality. Indeed, typically ocean models approach ocean motion by the Reynolds-averaged Navier-Stokes equations using the hydrostatic and Boussinesq assumptions. In this context inaccurate or inadequate models lead to structural uncertainty. Even if these equations were perfectly accurate models, they contain parameters, such as viscosity or diffusivity, which are not known precisely and cause parameter uncertainty. Boundary conditions on input variables, lead to parametric variability. For instance input variables such as velocities, experience at the surface atmospheric winds, which impose forcings on the velocities, and these atmospheric forcings are not precisely known.
Eventually these models are not solved exactly, but numerically through discretization methods and numerical approximations, which are subjected to numerical errors, and causes algorithmic uncertainty. For all these reasons, the predicted ocean currents are also uncertain. Uncertainties in the solutions, system outputs, due to the uncertainties in the system inputs is referred to as forward uncertainty quantification [1].
In this work, we are interested in evaluating the reliability of the outputs, ocean currents, regarding the transport they produce. This is a perspective slightly different to the classical one now described, in which ocean models adequacy is judged just against the velocity fields. Transport in the ocean surface is produced by fluids parcels that follow trajectories x(t) that evolve according to the dynamical system: where the position is described in longitude (λ) and latitude (φ) coordinates, that is, x = (λ, φ), and v represents the velocity field, which has two components determined by the zonal (u) and meridional (v) velocities. In longitude/latitude coordinates, the dynamical system in Eq.(1) can be rewritten as: where R is the Earth's radius. The system (1), which is a general expression encompassing the specific problem dealt with in (2), is a nonlinear non-autonomous dynamical system in 2D. Uncertainties in the velocities v(x, t), or more generally in the vector field, produce uncertainties in the solutions x(t). In general, as illustrated in figure 1, the exact model connecting two successive observations is not known. Only the initial observation, x 0 , and the final state, x * , are measurable. They are presented in red color in the figure. The evolution of initial conditions, X 0 , in a neighbourhood close to the initial observation, x 0 , towards a final state, X * , predicted by a model is expressed in the pink color in the same figure. The uncertainty of the model in representing the observations can be expressed in a number of ways. It may be defined by an absolute error E, which measures distance, in a certain metrics, between the final observation x * , which is considered a target and the computed prediction X * . This work proposes measures for this error and links the proposed uncertainty quantifier with dynamical objects present in the model (1). This is done without making any mathematical assumption about the observations. x 0 . The final observed state x * at time t * is referred to as the "target". The evolution law for these observations is unknown, but is approached by a model, which in our setting involves the velocities v(x, t) of the system (1). The evolution according to a model, of a neighbourhood of points, X 0 , close to the initial observation, x 0 , is illustrated with the pink color.
This paper is structured as follows. In section 2, we will discuss and develop an approach to uncertainty quantification recently taken for an ocean study case in [2]. We will show evidence of connections between the introduced definitions and the dynamical objects of system (1). In section 3, we will provide formal results to show how, for specific simple examples of the system (1), the definition of uncertainty quantification given in section 2, is able to highlight these dynamical objects. These results support one of the main findings of this article which is that stable invariant manifolds of hyperbolic trajectories provide a structure to forward uncertainty quantification. We will see in detail what is meant by this statement. Section 4 presents a discussion with further examples that illustrate the findings of this work and links uncertainty quantification to Lagrangian Descriptors and other Lagrangian indicators found in the literature. Finally, in section 5, we will provide the conclusions.

Transport uncertainty quantification in ocean models
In this section, we propose uncertainty quantification metrics. We work in the context of transport in ocean models, starting the discussion from the results presented in [2], which considers the performance of very high resolution tools for the monitoring and assessment of environmental hazards in coastal areas. In particular, this work, following the scheme presented in figure 1, quantifies uncertainties of several high resolution hydrodynamic models in the area of Gran Canaria, by measuring an error between the observations and predictions for the evolution of a diesel fuel spill event, well documented by port authorities and tracked with very high resolution remote sensing products. The pollution event was produced after the collision of the passenger ferry 'Volcán de Tamasite' with the Nelson Mandela dike in La Luz Port on April 2017. After the crash, supply pipes along the dike were broken and diesel fuel poured into the sea. SAR Sentinel 1 images were available in the area approximately one day and a half after the accident from which it was possible to identify the spill. For that period ocean currents in the area were available from two sources. One was the Copernicus Marine Service model for the Iberian-Biscay-Irish region (CMEMS IBI-PHY, IBI hereafter) and other was a very high resolution model of Puerto de la Luz setup by Puertos del Estado currently implemented and running operationally in different Spanish Port Authorities within the SAMOA port forecast system [3]. In [2], for each model, uncertainties are quantified by measuring an error with respect to a target, a 'ground truth', which should be recovered by the model. In particular, following the scheme proposed in figure  domain. This domain is divided into sub-domains and in each one a small circular blob with radius r = 3.5 · 10 −3 ( • ), is evolved in a time interval from the initial time t 0 . This is represented just for one of the sub-domains. The blob distorts while it approaches the target observed spill. A way to measure the 'proximity' between the blob and the target spill at the final time, t * , is to compute the distance between the centroid of the evolved blob, (c m ), and that of the observed spill, (c * g ). Mathematically this can be expressed as: Here, · denotes the modulus of the vector. The position of the centroid of the evolved spill, c m , depends on time t, while the centroid of the ground value slick c * g does not, because it is an observation at a final time, t * . The centroid of a finite set of N points {x k } k∈N ∈ R n is defined as: where   to a 20 × 20 mesh-grid and decreases the radius of initial blobs to 3.5 · 10 −5 ( • ). Figure 3b) makes visible an underlying structure which is directly related to the invariant dynamical structures, as we will shown later. In the limit r → 0, expression (3) is rewritten as: where x(t) is a trajectory of a fluid parcel, a solution to the system (1), with initial position x 0 at each cell in the grid. Figure 4a) illustrates the results of this calculation in a very fine grid.
One of the goals of this paper is to establish connections between the uncertainty quantifier (5) and invariant dynamical objects that control transport in vector fields. These are geometrical objects that organize particle trajectories schematically into regions corresponding to qualitatively distinct dynamical behaviors. In the context of fluid dynamics these objects are referred to as Lagrangian Coherent Structures (LCS). An essential ingredient of the LCS are hyperbolic trajectories characterized by high contraction and expansion rates. Directions of contraction and expansion define, respectively, stable and unstable directions, which are, respectively, related to the stable and unstable manifolds. In the context of the ocean model and the event described above, Garcia-Sánchez et al. have shown in [2] that there exists a hyperbolic trajectory located on the coastline very close to the accident point, in a detachment configuration, which is related to the phenomena of flow separation. Under this configuration, the stable manifold of the hyperbolic trajectory is aligned with the coast, and the unstable manifold is transversal to it. Figure 5 illustrates, respectively, in blue and red the stable and unstable manifolds of the separation trajectory for this particular case that we studied. Any blob placed in the neighbourhood of the separation trajectory, eventually evolves to become aligned with the unstable manifold, which is an attracting material curve. The configuration of the unstable manifold in Figure 5, is consistent with the observed spill marked in red in Figure 2. Indeed, the observed spill has evolved to be completely aligned with the unstable manifold. The invariant dynamical structures displayed in Figure 5 have been obtained by means of the Lagrangian Descriptor [4,5,6], which measures arc-length of trajectories: In [7], for some generalizations of the expression (6), is proved that stable and unstable manifolds are aligned with singular features of this function. The forward integration highlights the stable manifold, while the backwards integration highlights the unstable manifold. In this work the notion of singular feature is related to an undefined directional derivative in a direction transverse to the manifold curve. In Figure 5, red and blue features are placed on the singular features of M .
The metric given in (5) posses similarities with the forward definition of the function M : The analogy between structures obtained from expressions (5) and (7), is confirmed from figure 4.
Panel a) displays expression (5) results, while panels b) and c) display Eq. (7) results for τ = 1.5 and 2 days respectively. Results presented next in this section are displayed following this type of representation.
Our analysis in this section follows the spirit of the work by [6,7,8,9]. We assume the definition of singular features given there, by considering that these are features of L U Q on which the transversal derivative is not defined. We will prove, for simple selected examples, that stable manifolds are aligned with those singular features of L U Q .
Finally, before beginning our discussion the definition (5) is generalized as follows: Here, x 0 is the initial condition of the trajectory (x 1 (t), x 2 (t), .., x n (t)). In the particular case described in the previous section n = 2, which corresponds to the ocean surface, and p = 2. The coordinates of the target c * g are (x * 1 , x * 2 ) An alternative for this expression that we will use is the following:

The autonomous saddle point
The first example that we analyze is the vector field that corresponds to the Hamiltonian linear saddle case for which x ∈ R 2 , satisfying x = (x, y) and the equations of motions are: where λ > 0. For any initial observation (x 0 , y 0 ), we consider the unique solution of this system passing through the condition (x 0 , y 0 ), which is: For this example, the origin (0, 0) is a hyperbolic fixed point with stable and unstable manifolds: For simplicity we assume, without loss of generality, that t 0 = 0 (this is possible for autonomous systems). We consider the target in the position (x * , y * ) and apply (11) to (9) to obtain: Regrouping terms, we get In this expression, ω = e −λt , which always satisfies ω > 0. We explore separately the first and second terms. For the factor |1 − aω| p there exists a t L such that if t > t L , then aω << 1 and (1−aω) > 0. This is always the case for a < 0 and is a plausible assumption for a > 0, if a << ω −1 .
We recall that a = x * /x 0 and that therefore such t L exists if x 0 = 0. In this case positiveness is guaranteed for sufficiently large t, i.e, a Taylor series around ω = 0, attained if t 1 and t > t L , is performed for the binomial: We recall that p ≤ 1 and ω = e −λt . This yields, We analyse next the second term in Eq. (14). The sign of (ω − b) depends crucially on the sign of b given that for t 1 the term ω can be as small as we like. Let us consider t L such that if t > t L and b < 0, then (e −λt − b) > 0. The following is valid for t above this lower value: Therefore, when t 1, b < 0 and e −λt 1, Now, we consider b > 0 and t > t L satisfying (ω − b) < 0. Therefore: Therefore, when t 1, b > 0 and e −λt 1, Finally, we conclude that Therefore, since b = y * /y 0 Thus, we can approximate L U Q as where the dominant term is |x 0 | p e λtp . Hence, to leading order, the stable manifold at x = 0 is aligned with a singular feature of L U Q for 'sufficiently large' t. This statement however must be considered in the sense that expression expresses a good approximation for L U Q as far as x 0 = 0, and that for any |x 0 | > 0 is valid for a sufficient large t satisfying, t > t L .
We have shown that Eq. (9) is able to highlight the stable manifolds for the autonomous saddle point. Next, we analyse what happens when we consider Eq. (8), and in which cases it provides information about the stable manifold of the autonomous saddle point. When we apply (11) to (8), it yields We remark that our following calculations will be restricted to p > 1 and integer. Rewriting the Recalling that p > 1, ω > 0 and that (1 − aω) > 0 for sufficiently small ω, i.e. sufficiently large t, and the case b < 0 and (ω − b) > 0, a Taylor series around ω = 0, for p integer number, for the second factor: Alternatively, considering the case b > 0 and (ω − b) < 0, a Taylor series around ω = 0 for the second factor is as follows: Therefore in general: We notice that since p > 1, the terms ω p−1 , ω p have positive exponents and if ω 1 in principle are much smaller than the first term ω −1 and therefore: Hence, the stable manifold at x = 0 is aligned with a singular feature of this approximation to L U Q for 'sufficiently large' t. This statement though, must be considered with care, because in Eq. (22) the terms in ω p−1 , ω p , are multiplied by |x 0 | 1−p which has a negative exponent, and therefore have a singularity at x 0 = 0. For this reason neglecting these terms versus the first one if we are very close to x 0 = 0 would require checking that the products, such as, |x 0 | 1−p ω p−1 are really small.
This implies that Eq. (23) is correctly approximating L U Q for ω 1 as far as we are sufficiently away from 0 in x 0 , i.e. |x 0 | 1−p ω p−1 1 =⇒ ω p−1 |x 0 | p−1 . In practice this condition is satisfied for any grid (x 0 , y 0 ) used in later figures that exclude x 0 = 0. Also we can state that for sufficiently large t, L U Q is very close in almost all the domain to the function given in Eq. (23), which possess a "singular feature" aligned with the stable manifold. Fig.7 a)

The autonomous rotated saddle point
This second case that we explore is the vector field of the rotated linear saddle for which the equations of motion are: The general solution to this system is: where a and b depend on the initial conditions x 0 and y 0 as follows: Here, a = 0 corresponds to the stable manifold of the hyperbolic fixed point placed at the origin and b = 0 corresponds to its unstable manifold.
For any initial observation (x 0 , y 0 ) and final target observation (x * , y * ), we introduce the solution (25) into Eq.(9) obtaining: We expand next the first term in Eq.(26). In particular, we consider the case in which a > 0 for which it is always possible to find a t > t L in which (ae λt + be −λt − x * ) > 0: Here we have used ω = e −λt . The Taylor expansion of the binomial when ω 1, i.e., when t 1, is Thus, The stable manifold at a = 0 is aligned with a singular feature also at a = 0.
As for Eq.(8) considering p > 1, p = 1/δ and ω = e −λt leads: As before the leading term in ω leads to: and similar considerations apply in neglected terms to the ones made for Eq. (22), but now regarding to singularities at a = 0. Fig.8 a)

Discrete maps
These findings in the previous two subsections can be easily extended to discrete time dynamical systems, which are also useful in applications. Discrete time dynamical systems are defined as maps.

The autonomous saddle point
Consider the following linear, area-preserving autonomous map: For an initial condition (x 0 , y 0 ), the unique solution of this system is: As for the continuous time case, the origin (0, 0) is a hyperbolic fixed point with stable and unstable manifolds: We apply (37) to (9) to obtain: Regrouping terms, we get Considering that the transformation e λt → λ n can be directly applied to all the calculations performed in the continuous time case, we recover from Eq. (14): When p < 1, the dominant term is |x 0 | p λ np . Hence, the stable manifold at x = 0 is aligned with singular features of L U Q for 'sufficiently large' iteration n.
The same is applicable to L U Q (8). It yields, , when, n 1 and pδ < 1.
We note that the stable manifold at x = 0 is aligned with singular features of L U Q for a 'sufficiently large' iteration n. The same issues as before regarding singularities on the manifold position apply. for p = 2. In b), we illustrate the same representation but for Eq.(9) when p = 0.1. It can be appreciated how the stable manifold is aligned with the singular feature.

The autonomous rotated saddle point
We consider the following discrete dynamical system: It is easy to see that the stable and the unstable manifolds are given by the vectors (1, −1) and (1, 1) respectively. The solution of this system yields to,    x n = aλ n + bλ −n y n = aλ n − bλ −n , Again considering the transformation e λt → λ n and the use of previous results for the continuous time case, we recover from Eq.(33): Therefore, the dominant term is L U Q (x, n, p) ∼ |a| p |λ| np .
Since a = x 0 +y 0 2 , there is a singular feature at x = −y, i.e., in the subspace generated by (1, −1). Hence, the stable manifold is aligned with a singular feature of L U Q .
The same is applicable to L U Q (8). It yields, Hence, the stable manifold is aligned with a singular feature of L U Q . Fig.10 a)

The Duffing equation
Results in the previous sections are generalized to the autonomous nonlinear case by means of the Moser's theorem [10]. This theorem applies to analytic two-dimensional symplectic maps having a hyperbolic fixed point or, similarly, to two-dimensional time-periodic Hamiltonian vector fields having a hyperbolic periodic orbit (which can be reduced to the former case considering a Poincaré map). The case of a Hamiltonian nonlinear autonomous system is a one-parameter family of symplectic maps, and therefore Moser's theorem applies. Following proofs sketched by [8,7] results may be extended to the case of non-autonomous nonlinear dynamical systems by utilizing results like the Hartman-Groβman theorem.
This section discusses further results on the Lagrangian Uncertainty Quantifier by considering the evaluation of (8) over a vector field obtained from the nonlinear periodically forced Duffing Prior to discuss outputs of (8) into Eq. (46), we discuss the structure of invariant manifolds of hyperbolic trajectories in Eq. (46) for the case = 0.1 and the persistence versus this time dependent perturbation of tori present in the unforced version of Eq. (46), i. e, = 0. In the perturbed case, the hyperbolic fixed point placed at the origin, becomes a hyperbolic periodic trajectory [4], and their stable and unstable manifold can be highlighted by the Equation 6 appeared in Section 2 [6]. Additionally, as discussed in [7] a scale factor 1/(2τ ) applied to (6) converts the expression to an average, which in compact Hamiltonian systems like (46) converges for τ → ∞ and when this convergence is observed, level curves correspond to invariant structures of the dynamical system. Convergence of means are computationally verifiable on tori, however on hyperbolic sets as discussed in [9], rounding computational errors practically prevent convergence. Figure 11   values of L U Q are found on the stable manifold that are optimal pathways towards the unstable manifolds. In panel b) we show the results for the target x * = (1, 0) which is within the right tori like structure highlighted in Fig. 11b). Accordingly, uncertainty values are very low for initial observations (x 0 , y 0 ) in the corresponding tori region, but very high for the tori like region at the left side. Indeed initial observations in this region will never go near a final observation in the right tori like structure, and therefore this model is structurally uncertain for those observations, i.e. the model is inadequate to represent those. Panel c) shows the computation of the uncertainty for a target x * = (0, 1) outside the geometry of the unstable invariant manifold displayed in Fig. 12a).
It is remarkable the persistence in all these examples of an structure on the uncertainty field with singular features linked to the stable manifold independently of the target value x * . Figure 14 shows the evaluation of Eq.(9) with p = 0.1, and different targets or integration periods. Panels a) and b) have the same integration periods, but different targets. In them the black arrow points out different spurious features that seem singular but do not correspond to any invariant manifold. Figure 14 c) shows the same than panel a) for an integration period of 50. The black arrow marks a region, which from the analysis of Fig. 11, is known to be covered by tori, while the structure attained from Eq.(9) does not highlight this. Indeed, the ergodic partition theory discussed in [11,12] and implemented in [7] for Lagrangian Descriptors such as that in Eq. (6), shows that this ability to highlight tori, requires averaging along trajectories, while expressions Eq.

Discussion
There is much interest in uncertainty quantification involving trajectories in ocean data sets. Recent efforts in this direction are [13,14,15,16]. Vieira et al. [16] have developed a clustering method to partition the space of trajectory data sets into distinct flow regions. This method contains free parameters and they use uncertainty quantification to assess the parameter dependence on the partitions they obtain. The work [15] discusses approaches for uncertainty quantification from a geometrical perspective, and uncertainty quantifiers based on distances are proposed, but no connections are proposed with invariant dynamical structures. Similarly in [14] uncertainty quantifiers based on a descriptive statistics of distances between modelled trajectories and observations are used, but no relations are presented between these and the invariant dynamical structures. Results appeared in [17] also discuss along these lines, but in the context of general models, not related to ocean data. In this work authors have established links between uncertainty quantification and invariant manifolds under explicit mathematical assumptions for the "true" model. The results in our work do not require such assumptions. Equations similar to (8) or (9) have been used in the literature to highlight Lagrangian structures in oceanic flows. For instance [18,19] have done so. Prants in [18] proposed to use the arc-length D to study the displacement of particles in the coasts of Japan.
Here, D represents the relative displacement of a particle from its initial position (x 0 , y 0 ) to certain final one (x f , y f ). This expression is the analogue to Eq.(8) with p = 2 and target adjusted to each initial condition. However, in this work no connections are established between Eq. (47) and uncertainty quantification.
This article explores the implications of a definition for Uncertainty Quantification recently proposed in oceanic contexts [2]. It is found that for this definition, which is associated to forward Uncertainty Quantification, stable invariant manifolds of hyperbolic trajectories of the underlying flow, provide a structure for it. That is, we found that the proposed Uncertainty Quantifier is a function that contains a very rich structure, which is related to these well known structures from dynamical systems theory. dynamical structures have been used to provide a framework for discussion on structural uncertainty, which is related to inadequate models. This vision enriches traditional descriptions of UQ on which structure is discussed just in terms of means and/or statistical moments of distributions.
The findings described in this article are particularly interesting because they have important environmental applications. Nowadays, multiple ocean data sources are available and in this context our results allow a quantitative comparison of the transport properties associated to them. In this way, discriminating the level of performance of different data sources will help to gain precision in the description of dispersion of contaminants, determination of waste and plastic sources, etc.