Lagrangian Data-Driven Reduced Order Modeling of Finite Time Lyapunov Exponents

This paper proposes a Lagrangian data-driven reduced order model (ROM) for an efficient and relatively accurate numerical simulation of the finite time Lyapunov exponent (FTLE) field. To generate the ba- sis, the new Lagrangian ROM explicitly uses Lagrangian data (i.e., the FTLE field) in addition to the Eulerian data (i.e., the velocity field); the Eulerian ROM, on the other hand, uses only Eulerian data. The Lagrangian ROM and the Eulerian ROM are compared in the numerical simulation of the quasi-geostrophic equations. It is shown that the new Lagrangian ROM outperforms the Eulerian ROM with respect to both Eulerian (velocity) and Lagrangian (FTLE) accuracy criteria. Furthermore, the online CPU time of the Lagrangian ROM is dramatically lower than the CPU time of the corresponding fine resolution numerical simulation. Our numerical investigation suggests that Lagrangian ROMs could be used as an efficient and relative ac- curate alternative to fine resolution numerical simulations to generate the FTLE field.


Introduction
Material transport by ocean flows is a ubiquitous process that has profound environmental and societal impact. Oceanic transport is important in climate and regulating the Earth's thermal budget [27]. Oceanic mixing and stirring has significant ecological implications [43]. The accurate prediction of ocean transport is critical for effective decision-making during potentially catastrophic pollution release events [57], such as oil spills [71,68] or radioactive waste [12].
Ocean flows comprise spatially complex and temporally transient features, and thus understanding how flow transport is organized and predicting where things go remains a formidable scientific challenge for a variety of reasons. Unsteady transport processes, due to nonlinear advection, require precise analysis so that material transport is accurately defined and predicted. Moreover, accurate numerical modeling of the flow field is highly complex due to nonlinearity and uncertainty.
However, in the past two decades, the fluid dynamics and dynamical systems communities have made progress in developing novel Lagrangian (trajectory-based) approaches to understanding transport in complex timedependent flows. These methods fall into four main categories: geometric, probabilistic, topological, and complexity. Under the umbrella of geometric methods are the techniques of invariant manifolds (of fixed points or larger invariant sets), and their extension to time-varying analogues, i.e., the most influential material surfaces, which have come to be known as Lagrangian coherent structures (LCSs). There are an increasing number of studies that apply various concepts of LCSs, based on the classic right Cauchy-Green tensor, to describe and predict the time evolution of Lagrangian features in geophysical systems [9,8,42,58,68,71,85,92,93]. The most widely used diagnostic for providing candidate LCSs is the finite-time Lyapunov exponents (FTLE) field, calculated from advected virtual particles based on the time-varying velocity field. The method of FTLE and LCS has proven to be particularly useful in studying transport in time-dependent systems and has found a variety of applications including aircraft performance [89], experimental fluid mechanics [74], heating, ventilation, and cooling efficiency [7,84], and medicine [44,88].
LCS methods are well suited to address challenges of quantifying and predicting oceanic transport by providing unequivocal definitions of transport structures, enabling rigorous understanding and new discoveries (e.g., [6]).
Such methods also provide new metrics for comparing predictions of different models. Moreover, with efficient analytics and effective visualizations, their data reduction capabilities support understanding and forecasting of material transport, thereby aiding decision making strategies.
One of the major hurdles faced by the development of the FTLE and LCS methods is their high computational cost. We schematically illustrate the computation of the FTLE field in Algorihm 1.
The computational cost of the FTLE field generated in Step (2) of Algorithm 1 is generally high. Indeed, a large number of particle trajectories must be integrated to construct a particle flow map [11,79]. Furthermore, long time integration might be needed to compute the FTLE field. Finally, it is often necessary to compute a sequence of FTLE fields in time to visualize unsteady flows [11].
Different approaches have been proposed to reduce the computational cost of the FTLE field. Most of these approaches assume that the velocity field in Step (1) of Algorithm 1 is given and only attempt to reduce the computational cost of the FTLE field generated in Step (2) of Algorithm 1 [11]. We argue, however, that reducing the cost of the velocity field used in the FTLE field computation is also needed. This is the case when, e.g., the FTLE field is used in realistic ocean and atmospheric settings that generally pose extremely heavy burdens on the available computational resources. Thus, the following natural question arises: "Can more efficient computational models be used to generate the velocity field in Step (1) of Algorithm 1?" In this paper, we investigate whether reduced order models (ROMs) can be used to significantly decrease the computational cost of a direct numerical simulation (DNS) of the velocity field in Step (1) of Algorithm 1, without compromising its numerical accuracy. Of course, reduced order modeling of fluid flow has become a well established research field (see, e.g., [45,47,67,73,94]). ROMs have also been used in ocean and atmospheric modeling (see, e.g., [21,87]). Thus, using ROMs to decrease the computational cost of the velocity field in Step (1) of Algorithm 1 is a natural strategy. This, in turn, will result in a significant decrease of the overall cost of the computation of the FTLE field in Algorithm 1. We note, however, that to ensure the success of this natural strategy, the following question needs to be addressed: "Is the ROM velocity field in Step (1) of Algorithm 1 accurate enough to ensure an accurate computation of the FTLE field in Step (2) of Algorithm 1?" Indeed, reducing the computational cost of Step (1) of Algorithm 1 is certainly desirable; however, if this yields an inaccurate velocity field approximation, then the FTLE field computed in Step (2) of Algorithm 1 would also be inaccurate [39,69].
In this paper, we propose a Lagrangian data-driven ROM for the computation of FTLE field. First, we describe in Section 3 an Eulerian ROM, i.e., a ROM that uses only Eulerian (velocity) data to generate the ROM basis, which is a standard approach in data-driven reduced order modeling [10,36,45,47,54,55,61,67,72,73]. The ROM that we propose in Section 4 employs a fundamentally different strategy: In addition to the Eulerian (velocity) data, this ROM also uses Lagrangian (FTLE field) data to generate the ROM basis. The resulting Lagrangian data-driven ROM that we put forth in Section 4 is, to our knowledge, new.
The rationale for constructing the Lagrangian data-driven ROM in Section 4 is straightforward: Standard ROMs (such as that in Section 3) are generally built from Eulerian data, e.g., streamfunction, velocity, or vorticity fields. Thus, standard ROMs can efficiently approximate the corresponding Eulerian fields. However, we conjecture that, to approximate Lagrangian data, e.g., the FTLE field, the ROM might benefit from including additional Lagrangian data. This is precisely how the ROM that we put forth in Section 4 is built: In addition to the vorticity field (which is Eulerian data), it also uses the FTLE field (which is Lagrangian data). In Section 5, we investigate whether this new Lagrangian ROM yields better approximations than the standard Eulerian ROM (described in Section 3).
Previous Relevant Work To our knowledge, there is only little work on reduced order modeling for the FTLE calculation [2,56,90].
The Lagrangian ROM that we propose is different from the ROMs used in [2]. Indeed, the ROM used in [2] for the FTLE computation is an Eulerian ROM, which utilizes the newly developed optimally time-dependent (OTD) modes [3,83] as ROM basis. In contrast, the ROM that we propose in this paper is a Lagrangian ROM, since it employs Lagrangian information (i.e., FTLE field) to build the ROM basis.
The Lagrangian ROM that we propose in this paper is also different from the model used in [56]. In fact, in Section 2.2 of [56], the authors state that they are not using an actual ROM to compute the FTLE field. Instead, they are using an "empirical method" to compute an approximation of the velocity field for a flow past a circular cylinder. This empirical method utilizes heuristic approximations of the ROM coefficients without building and integrating in time the actual ROM. In contrast, in this paper we propose ROMs that are built around standard Galerkin methods, in which a lowdimensional set of ordinary differential equations (ODEs) are first discretized in time and then run over the time interval of interest.
We also note that the new Lagrangian ROM that we propose is different from the ROM used in [90]. The two ROMs are built on the same fundamental principle: Lagrangian data can be used to generate ROMs for the computation of Lagrangian fields. We note, however, that the strategy used in [90] to incorporate Lagrangian information in the ROMs is different from the strategy that we use in this paper. Indeed, the FTLE field is used in [90] only to determine the number of ROM basis functions. In contrast, we explicitly use the FTLE field in the actual construction of the ROM basis.
Finally, we note that the new Lagrangian ROM that we propose is related to the ROM data weighting approaches [4,5,15,18,23,37,48,49,50,77] (which, in turn, belong to the wider class of goal oriented ROMs [13,15]), in which the input data (i.e., the snapshots) are weighted according to a userchosen criterion. The new Lagrangian ROM is also related to the structurepreserving ROM approaches [16,17,35]. However, we emphasize that, to our knowledge, none of these techniques use the FTLE field, as we do in this paper.
The rest of the paper is organized as follows: In Section 2, we describe the FTLE field and outline some of the algorithms currently used for its computation. In Section 3, we outline the derivation of Eulerian ROMs. In Section 4, we propose a novel Lagrangian ROM. In Section 5, we investigate whether the new Lagrangian ROM and the Eulerian ROM can be used for an accurate computation of the FTLE field. Finally, in Section 6, we present conclusions and outline future research directions. 5

FTLE Computation
Constructing efficient algorithms is critical for the development of the FTLE methods (see, e.g., Section 7 of the survey in [42]). In this section, we first briefly describe the FTLE (Section 2.1) and then we outline several of the recent algorithms for the FTLE field computation (Section 2.2).

FTLE Background
The FTLE field has become a useful Lagrangian diagnostic to identify LCSs in fluid flows, providing useful information about large-scale flow patterns and transport and mixing phenomena in flow domains [42,89]. Next, we briefly describe the FTLE field.
Consider a Euclidean domain Ω ⊆ R n , where n is the dimension of the spatial domain. A typical dynamical system is given in the form of a velocity field v(x, t) defined on Ω × R. In this case, the trajectories are given by the solutions of the ODE systemẋ = v(x, t). Each trajectory is a function of time but it also depends on the initial position x 0 and the initial time t 0 . How the trajectory changes as the initial time and the initial position change is a central interest in this paper. For this reason, the trajectory which is solution of the initial value problem is denoted as x(t; t 0 , x 0 ), which emphasizes the explicit dependence on the initial condition. It is also convenient to define the flow map to further highlight the dependence on the initial position x 0 . For a given initial time t 0 and a given final time t, the flow map is the function Consider two particles, simultaneously released at time t 0 ; one at location x, the other at location x + δx. Under the effect of the flow map, the small displacement vector between two particles, δx, changes. As illustrated in Fig. 1(a), after an elapsed time T , the new vector between the two particles is The closer the initial displacement vector δx is to the ξ 2 direction, the more it will be stretched to that maximum perturbation.
where φ t 0 +T t 0 is the flow map for the vector field (1) from time t 0 to t 0 + T , is the Jacobian of the flow map, and · is the usual Euclidean norm. Consider the right Cauchy-Green strain tensor, Since the flow domains for both numerical experiments in Section 5 are twodimensional, we consider incompressible two-dimensional flows (i.e., n = 2). The eigenvalues λ i and normalized eigenvectors ξ i of C satisfy [42], where the (x, t 0 , T ) dependence of C, λ i , and ξ i is understood. As illustrated in Fig. 1(b), the two eigenvectors, ξ 1 and ξ 2 , are carried along by the flow to the two vectors r 1 and r 2 , respectively, where The lengths of r i are scaled by a factor √ λ i compared with the normalized eigenvectors ξ i . The maximum possible separation between the released particles after a time interval T , assuming a sufficiently small initial distance where λ max = λ 2 . The FTLE, with t 0 and T fixed, is considered a scalar field of the Lyapunov exponent as a function of initial position, x, The LCS are ridges in the FTLE field. They are co-dimension 1 manifolds that maximize the FTLE in the transverse direction. This notion is more formally defined in [59,89]. The objective is to generalize the concept of stable and unstable invariant manifolds (of hyperbolic invariant objects) to time-dependent systems. Whether in the ocean, in the atmosphere, or in any other dynamical system, LCS are the instantaneous most important material surfaces which serve as a time-varying scaffold for Lagrangian transport.

FTLE Algorithms
In this section, we summarize some of the approaches used to reduce the computational cost of the FTLE field generated in Step (2) of Algorithm 1. We start by presenting standard methods for the FTLE field computation. Then, we outline some of the algorithms that aim at increasing the efficiency of the FTLE computation.
Standard FTLE Algorithms One of the most popular FTLE algorithms centers around the finite difference approximation of the Jacobian of the flow map, Dφ t 0 +T t 0 , which is used to define the Cauchy-Green strain tensor C in (4) (see, e.g., Section 3.3 in [42] and Section II in [58]). Next, we briefly describe this algorithm.
Let us consider a Cartesian mesh with grid points (x ij , y ij ) ∈ R 2 . We compute the final position, In practice, the final position is often obtained by numerically integrating the velocity field (1). The matrix representation of the Jacobian of the flow map based on a finite difference approximation is given by where δ = max δ x , δ y is the maximum grid size. To compute the FTLE (8), one can then find the largest singular value of the matrix M [26], which approximates Dφ t 0 +T t 0 for a sufficiently small grid spacing δ. Alternatively, as suggested by [41], one can search for the largest eigenvalue of M ⊤ M. The algorithm above is therefore identical to the algorithm proposed as the twodimensional direct Lyapunov exponent (DLE) in [41], the two-dimensional algorithm of [89] or the n-dimensional algorithm of [59].
The major advantage of this algorithm is its robustness. Regardless of the grid size δ, an existing LCS passing between two grid points creates a visible shade in the approximated FTLE field (see, e.g., Fig. 1 in [58]). Thus, even on coarse meshes, one is able to identify a rough estimate of the LCS structure. We note, however, that the FTLE can be underestimated if the grid size is too large.
We also note that the FTLE is typically computed at grid resolutions much finer than the flow resolution [51,89]. Indeed, the FTLE indicates long term particle separation and the LCS have a fine, foliated structure, even for flows defined on a coarse mesh. In time-periodic systems, the LCS are known to form tangles and arbitrary small lobes [75]. For arbitrary time dependence, the LCS can also form meanders of arbitrary small size, provided that the integration time is long enough [28,76]. Thus, one needs to use a fine mesh resolution to capture the small spatial features of the FTLE and LCS, even when the flow mesh resolution is relatively coarse.
Other FTLE algorithms are mentioned in, e.g., Section 3.3 in [42]. Rather than integrating particle trajectories (a Lagrangian approach), an Eulerian approach to computing the flow map, and subsequently the FTLE, was proposed in [60]. In two dimensions, the method involves solving two level set (Liouville) equations, which are PDEs. An advantage over Lagrangian-type methods is that no interpolation of the velocity field (to non-mesh points) is needed. Another approach to computing FTLE involves a set-oriented, transfer operator approach [91]. The finite-time flow map has associated with it an operator which evolves densities, known as the Perron-Frobenius operator. The FTLE value for a finite-size set can be defined in terms of the covariance of an initially uniform density and its image under the Perron-Frobenius operator. This approach provides new interpretations of the FTLE and a link between geometric and probabilistic approaches.
Efficient FTLE Algorithms Adaptive mesh refinement was used in [34] to significantly reduce the number of particle integrations needed for the computation of FTLE fields in two and three dimensions. Mesh adaptivity was utilized in [58,64] to significantly reduce the number of particle integrations needed for the computation of FTLE fields in two and three dimensions. An unstructured high-order hp/spectral-element method and a corresponding high-order particle tracking method for FTLE field computation were developed in [80]. An efficient method for computing FTLE fields in unsteady flows was proposed in [11], where the particle flow map was approximated with bidirectional and unidirectional intermediate flow maps, eliminating redundant particle integrations in neighboring flow map calculations. This flow map composition method was made more accurate in [62] by approximating the intermediate short-time flow maps using polynomial basis functions (e.g., Legendre polynomials) which are then composed to form a long-time flow map. In [19], the computation of FTLE fields for bluff body flows was accelerated using GPUs and APUs along with the method of intermediate flow maps proposed in [11].
FTLE Software General purpose software for the FTLE computation has also been developed. ManGen [57,59] calculates the FTLE field and advects material curves in two-dimensional velocity fields. An updated version, Newman [28], is a parallel code that computes the FTLE fields and supports analytic and dataset velocity definitions. Newman is the software used for the results in this paper and is available at https://github.com/RossDynamics/Newmanv3.1. FlowTK [1] is a CUDA parallel computing platform that calculates the FTLE in two and three dimensions on Cartesian and unstructured grids. LCS Tool [70] is a computational platform for LCS and FTLE fields. Finally, the Particle Tracking and Analysis TOolbox (PaTATO) [30] is a Matlab code that aim at making Lagrangian particle tracking and analysis methods widely accessible.

Eulerian ROMs
In this section, we outline the derivation of Eulerian ROMs for FTLE computation. Since both numerical experiments in Section 5 use the quasigeostrophic equations (QGE) as mathematical model, in what follows we will also use the QGE to illustrate the ROM development. We emphasize, however, that the new Lagrangian ROM that we propose in Section 4 could be applied to other mathematical settings (e.g., the Navier-Stokes equations or the Boussinesq equations).

Quasi-Geostrophic Equations (QGE)
In what follows, we use the non-dimensional QGE as mathematical model. The QGE, which have been a computational paradigm for FTLE/LCS testing [20,63,92], have the following form: where ω is the vorticity, ψ is the streamfunction, Re is the Reynolds number, and Ro is the Rossby number. The velocity can be computed from the streamfunction according to the following formula, e.g., used in (1): For details regarding the parameters and non-dimensionalization of the QGE (10), we refer the reader to, e.g., [29,53,65,81]. Indeed, the QGE (10) are equations (4) and (10) in [81] and equation (1) in [65].

Proper Orthogonal Decomposition (POD)
One of the most popular reduced order modeling techniques is the POD. In this section, we briefly outline the POD for the QGE (10). For more details, we refer the reader to, e.g., [45,47,67,73]. The POD starts by collecting the snapshots {ω 1 h , . . . , ω M h }, which are, e.g., finite element (FE) approximations of the vorticity in the QGE (10) at M different time instances. The POD seeks a low-dimensional basis that approximates the snapshots optimally with respect to a certain norm. Probably the most popular inner product is the L 2 inner product : The solution of the resulting minimization problem is equivalent to the solution of the eigenvalue problem where ϕ j and λ j denote the vector of the FE coefficients of the POD basis functions and the POD eigenvalues, respectively, Y denotes the snapshot matrix, whose columns correspond to the FE coefficients of the snapshots, M h denotes the FE mass matrix, and N is the dimension of the FE space. The eigenvalues are real and non-negative, so they can be ordered as follows:λ 1 ≥ λ 2 ≥ . . . ≥ λ R ≥ λ R+1 = . . . = λ N = 0. The POD vorticity basis consists of the normalized functions {ϕ j } r j=1 , which correspond to the first r ≤ N largest eigenvalues. Thus, the ROM vorticity space is defined as X r := span{ϕ 1 , . . . , ϕ r }. We follow [81] and define the POD streamfunction basis as the normalized functions {φ j } r j=1 , which are chosen such that −∆φ j = ϕ j , j = 1, . . . , r .

Galerkin ROM (G-ROM)
The ROM approximations of the vorticity and streamfunction are where {a j (t)} r j=1 are the sought time-varying ROM coefficients. We emphasize that, with the choices in (14)- (16), once the coefficients a j are determined from (10a), equation (10b) is automatically satisfied. Replacing the vorticity ω by ω r in the QGE (10a) and then using a Galerkin projection onto X r , we obtain the Galerkin ROM (G-ROM) for the QGE: ∀ i = 1, . . . , r, The G-ROM (17) yields the following autonomous dynamical system for the vector of time coefficients, a(t): where b, A, and B correspond to the constant, linear, and quadratic terms in the numerical discretization of the QGE (10), respectively. The finite dimensional system (18) can be written componentwise as follows: For all i = 1, . . . , r, A im a m (t) + r m=1 r n=1 B imn a m (t) a n (t), where

Lagrangian ROMs
In this section, we propose new ROMs for the FTLE computation. In order to construct the ROM modes, the new ROMs use not only Eulerian data (vorticity field) as the standard G-ROM (17)

FTLE-ROMs
The main tools that we use to construct the new ROMs for FTLE computation are two novel inner products, which are fundamentally different from the standard L 2 inner product (12) used to develop the Eulerian G-ROM (17). These new inner products are FTLE weighted inner products (·, ·) F T LE , which aim at including both Eulerian data (i.e., the vorticity field) and Lagrangian data (i.e., the FTLE field) in the ROM basis generation. These two inner products generate ROM basis functions that are different from the standard POD modes (built with the standard L 2 inner product). These two new bases will yield two new ROMs, which we present in Sections 4.1.1 and 4.1.2.

α-ROM
The first FTL weighted inner product that we consider is where λ 1 and λ 2 are FTLE fields. The parameter α in (23) is a weighting parameter that determines the weight associated with the FTLE field's contribution to the inner product: When α = 0, the FTLE field does not play any role, so the inner product (23) is the standard L 2 inner product (12). When α > 0, the FTLE field plays a significant role: The higher the α value, the more important the FTLE contribution to the inner product (23). We use the new FTLE weighted inner product (23) to generate the ROM basis for a new FTLE-ROM. First, we collect snapshots that consist of both vorticity and FTLE FE approximations. (Note that this is different from the standard ROM basis generation in Section 3, where only vorticity snapshots were collected.) Then, we construct the new ROM basis that approximates the snapshots optimally with respect to the norm (Again, we note that this is different from the approach in Section 3, which uses the norm ω .) Finally, from the resulting ROM basis functions, we only use their vorticity components in the ROM (17).
The new ROM for the FTLE computation is the G-ROM (17) in which the ROM basis is generated by using the new FTLE weighted inner product (23) instead of the standard L 2 inner product (12). In what follows, we will denote with α-ROM the resulting new ROM.

λ-ROM
The second FTLE weighted inner product that we consider is where λ is the time average of the FTLE field, λ. When we use the FTLE weighted inner product (25) to generate the ROM basis, these basis functions approximate the snapshots optimally with respect to the norm Note that, by definition, the FTLE field (8) is always positive. Thus, the FTLE inner product (25) and the associated norm (26) are well defined. The new ROM for the FTLE computation is the G-ROM (17) in which the ROM basis is generated by using the new FTLE weighted inner product (25) instead of the standard L 2 inner product (12). In what follows, we will denote with λ-ROM the resulting new ROM. Table 1 summarizes the Eulerian ROM and Lagrangian ROMs, together with the corresponding inner products.

Eulerian ROMs or Lagrangian ROMs?
In reduced order modeling, the basis is generally built from Eulerian data, e.g., velocity, vorticity, or streamfunction. Thus, ROMs generally are Eulerian ROMs. Of course, these Eulerian ROMs can efficiently approximate the corresponding Eulerian fields. However, Lagrangian fields (e.g., the FTLE field) might be better approximated by Lagrangian ROMs, i.e., ROMs whose bases are built from both Eulerian and Lagrangian data. The new α-ROM and λ-ROM are Lagrangian ROMs, since the FTLE weighted inner products (23) and (25) use both Eulerian data (i.e., the vorticity ω) and Lagrangian data (i.e., the FTLE λ). The input Eulerian data (i.e., the vorticity field) help the resulting FTLE-ROMs yield an accurate approximation of the output Eulerian data (i.e., the vorticity, streamfunction, and velocity fields computed by the FTLE-ROM). On the other hand, the input Lagrangian data (i.e., the FTLE field) "nudge" the FTLE-ROMs toward an accurate approximation of the output Lagrangian data (i.e., the FTLE field obtained from the FTLE-ROM velocity field). Thus, we expect the new α-ROM and λ-ROM to yield more accurate FTLE approximations than the standard G-ROM (17). In Section 5, we investigate the new α-ROM and λ-ROM, as well as the standard G-ROM, in the FTLE computation of a test problem that uses the QGE as mathematical model.
When comparing the pros and cons of Lagrangian ROMs and Eulerian ROMs, one natural question is whether the Lagrangian snapshots used to generate the Lagrangian ROMs are actually redundant. Indeed, the FTLE field is generated from the velocity field, which is generated from the streamfunction field (according to equation (11)), which in turn is generated from the vorticity field (according to equation (10b)). Thus, at first glance it would seem that using both vorticity and FTLE data to generate the ROM basis is a redundant process. Next, we show that this is not the case.
Adding FTLE information to the snapshot matrix generally does not significantly change the eigenvalues and eigenvectors of the eigenproblem (13), from which the ROM basis is generated. We emphasize, however, that adding FTLE information generally changes the relative ordering of these eigenvalues and eigenvectors. This, in turn, yields significantly different truncated ROM bases. This is clearly illustrated in Fig. 2, which, for the QGE test problem used in the numerical investigation in Section 5, displays the ROM basis functions ψ 10 , ψ 20 , and ψ 30 generated with the standard L 2 inner product (12) (i.e., the G-ROM basis functions), the new FTLE weighted inner product (23) (i.e., the new α-ROM basis functions), and the new FTLE weighted inner product (25) (i.e., the new λ-ROM basis functions). The α-ROM basis functions are completely different from the other two sets of ROM basis functions. The λ-ROM basis functions are different from the G-ROM basis functions, although the differences occur mainly for ψ 30 .
We also note that adding FTLE information to the snapshot matrix is similar to the use of snapshot difference quotients in general ROMs. Indeed, although the snapshot difference quotients are just linear combinations of the standard snapshots, including them or not in the set of snapshots can have a significant effect on the resulting ROMs [48].

Numerical Results
In this section, we investigate whether the new Lagrangian ROMs (i.e., the α-ROM and the λ-ROM) can be used for an accurate computation of the FTLE field. To this end, we test the α-ROM and the λ-ROM in the numerical simulation of the QGE (10) without an analytical solution. For comparison purposes, we also test the Eulerian G-ROM. As a benchmark, we use the full order model (FOM) (see Algorithm 2), which also serves to generate the snapshots for the ROMs.

Criteria
To compare the three ROMs (i.e., α-ROM, λ-ROM, and G-ROM), we compare the results from Algorithm 3 with those from Algorithm 2. To this end, we use two fundamentally types of criteria: An Eulerian criterion, i.e., the L 2 norm of the errors in the time-averaged streamfunction and a Lagrangian criterion, i.e., the L 2 norm of the errors in the timeaveraged FTLE field computation

Test Problem Setup
In this section, we consider the QGE (10) with a symmetric double-gyre wind forcing given by F = sin π (y −1) , which yields a four-gyre circulation in the time mean. This test problem has been used in, e.g., [22,38,46,66,81,82] as a simplified model for more realistic ocean dynamics. As shown in [38], although a double-gyre wind forcing is used, the long term time-average yields a four-gyre pattern. Indeed, the time-averaged streamfunction contour plot in Fig. 3 clearly displays the four-gyre pattern, 20 which cannot be inferred from a visual inspection of the instantaneous contour plots in the same figure. In [81,82], it was shown that standard coarse resolution numerical simulations or standard ROMs cannot recover the correct four-gyre pattern. Thus, in addition to the Eulerian criterion in Section 5.1, we also use the recovery of the four-gyre pattern to assess the new α-ROM and λ-ROM. In the QGE (10), we use the same parameters as those used in [46,81,82]

Snapshot Generation
For the FOM (see Algorithm 2), we use a spectral method with a 257 × 513 spatial resolution and the adaptive RKF45 time discretization with 0.01 time step. We follow [81,82] and run the FOM up to the dimensionless time T max = 80. In Fig. 4, we plot the time evolution of the spatially averaged kinetic energy, E(t). Figure 4 (see also Fig. 1 in [81]) shows that the flow converges to a statistically steady state, after a short transient interval that ends around t = 10. In Fig. 3, we plotted the instantaneous contour plot for the streamfunction field at t = 40 and t = 60. We emphasize that, although t = 40 and t = 60 are well within the statistically steady state regime, the flow displays a high degree of variability. To generate the ROM basis, we collect 701 snapshots in the time interval [10,80] (on which the statistically steady state regime is attained) at equidistant time intervals. In Fig. 5, we plot the eigenvalues of the correlation matrix used to determine the ROM basis for the G-ROM, λ-ROM, and α-ROM. We note that, for the G-ROM, these eigenvalues represent the percentage of the system's kinetic energy captured by the first ROM modes. For λ-ROM and α-ROM, however, these eigenvalues do not have a straightforward physical interpretation, since the inner products used to assemble the corresponding correlation matrices are not the standard L 2 inner product.

Lagrangian Investigation
In Table 3, we list the L 2 norm of the errors in the time-averaged FTLE (28) for G-ROM, λ-ROM, and α-ROM for different r values. The results in Table 3 show that the α-ROM with high α values (i.e., α = 10 3 and α = 10 4 ) consistently outperform the λ-ROM and the G-ROM for all r values, but especially for the small r values. We also note that the relatively high magnitudes of the errors in Table 3 are due to the errors in the ROM velocity field approximations. Decreasing the magnitude of the errors in the ROM velocity field approximations (e.g., by increasing the ROM dimension, r) would probably decrease the magnitude of the errors in the FTLE field approximation in Table 3. Table 3: L 2 norm of the errors in the time-averaged FTLE (28) for G-ROM (second column) λ-ROM (third column), and α-ROM for α = 1 (fourth column), α = 10 (fifth column), α = 10 2 (sixth column), α = 10 3 (seventh column), and α = 10 4 (eighth column).

ROM Computational Efficiency
To generate the FTLE fields, the following computational environments were used: One computing cluster composed of 5 nodes, each node comprised of dual, quad core, hyperthreaded 2.4GHz Intel Xeon E5620 CPUs (16 processor threads), 24GB RAM, and a 40Gbps InfiniBand host card and cable. Five nodes at 12 threads per node were utilized, for a total of 60 threads, and each thread was allocated 4749mb of memory. To generate the ROM velocity fields, the following computational environment was used: One Dell laptop with a single 2.70 GHZ CPU, running on a 64-bit linux system.
The computational cost of the DNS and ROMs has two components: (i) the computational cost of generating the underlying velocity field; and (ii) the computational cost of generating the FTLE field. To generate the velocity field, the online CPU times of both the Lagrangian ROMs (i.e., α-ROM and λ-ROM) and Eulerian ROM (i.e., G-ROM) are O(10 3 ) lower than the CPU time of the corresponding DNS. To generate the FTLE field, the ROM computational cost is higher than the DNS computational cost, especially when relative low r values are used. We note, however, that this computational cost increase is much lower for the Lagrangian ROM (i.e., α-ROM with α = 10 4 ) than for the Eulerian ROM (i.e., G-ROM): For example, for r = 30, the FTLE field computation with the Eulerian ROM is about 35% slower than the corresponding DNS computation; on the other hand, the FTLE field computation with the Lagrangian ROM is only 12% slower than the DNS computation.
To conclude, the Lagrangian ROM and the Eulerian ROM are significantly more efficient than the DNS in the velocity field computation. In the FTLE field computation, the Lagrangian ROM is only slightly slower than the DNS; the Eulerian ROM, however, is significantly slower than the DNS. We believe that this increase in CPU time for ROMs is because the ROM velocity field is less accurate than the DNS velocity field; we plan to investigate this in a future study.

Conclusions and Outlook
In this paper, we proposed a Lagrangian data-driven ROM for an efficient and relatively accurate numerical simulation of the FTLE field. We explicitly used Lagrangian data (FTLE field) in the inner product utilized to construct the ROM basis of the new Lagrangian ROM. We compared the Lagrangian ROM with the Eulerian ROM in the numerical simulation of the QGE.
The new Lagrangian ROM outperformed the Eulerian ROM with respect to both Eulerian (velocity field) and Lagrangian (FTLE field) accuracy criteria. The online CPU time of the Lagrangian ROM was dramatically lower than the CPU time of the corresponding fine resolution numerical simulation. Our numerical investigation showed that Lagrangian ROMs could be used as an efficient and relatively accurate alternative to fine resolution numerical simulations to generate the FTLE field.
The Lagrangian ROMs investigated in this paper (see also [2,56,90]) showed promise for speeding up the computation of Lagrangian fields, such as the FTLE. There are, however, numerous research directions that could provide improvements both in the efficiency and the accuracy of these Lagrangian ROMs. Next, we outline some of these research directions.
The Lagrangian α-ROM proposed in this paper clearly outperformed the Eulerian ROM in the numerical simulation of the QGE, especially for higher α values. However, finding the optimal α value is still an open question. One could use mathematical approaches (e.g., error estimates that yield optimal parameter scalings), physical approaches (e.g., physical scaling arguments), or a combination of those to determine the optimal α value.
To construct the Lagrangian λ-ROM, we used the FTLE inner product (25). This can be interpreted as spatial weighting, since the POD basis functions are constructed by weighting more the spatial regions where the temporal average FTLE field has higher values. We could, however, use instead temporal weighting, i.e., use the standard inner product (12), but include multiple copies of snapshots with a high FTLE spatial average. We note that temporal weighting for ROMs of fluid flows (without mixing) was used in [18]. Extending temporal weighting to ROMs for FTLE field computation could also be investigated.
The new Lagrangian ROM and the corresponding inner product (23) could also be adapted to and used for speeding up the computation of other structures that characterize transport and mixing. For example, instead of geometric approaches (such as the FTLE/LCS fields), one could approximate probabilistic measures, such as the almost invariant sets [25,24,31,33,32,40,91].
In this paper, we used the POD basis to generate the new Lagrangian ROM. We note, however, that we could use other types of bases that are appropriate for structure dominated systems, such as the dynamic mode decomposition (DMD) [78,86].
The new Lagrangian ROMs could be used to speed up the computation of the FTLE field for other test problems, such as the flow past a cylinder [14,52,56]. For example, Lagrangian data such as the FTLE field can shed new light on mixing for the flow past a cylinder problem, which is relevant to engineering and geophysical applications. Since Eulerian ROMs have been used for efficient and relatively accurate approximations of the velocity field for flow past a cylinder, Lagrangian ROMs could also be used for speeding up the computation of the corresponding FTLE field. 27