Path-integrated Lagrangian measures from the velocity gradient tensor

Nonlin. Processes Geophys., 20, 987–991, 2013 www.nonlin-processes-geophys.net/20/987/2013/ doi:10.5194/npg-20-987-2013 © Author(s) 2013. CC Attribution 3.0 License. Nonlinear Processes in Geophysics Open Access Path-integrated Lagrangian measures from the velocity gradient tensor V. Pérez-Muñuzuri1and F. Huhn2 1Group of Nonlinear Physics, Faculty of Physics, University of Santiago de Compostela, 15782 Santiago de Compostela, Spain 2Institute for Mechanical Systems, ETH Zürich, Tannenstrasse 3, 8092 Zürich, Switzerland Correspondence to: V. Pérez-Muñuzuri ([email protected]) Received: 2 May 2013 – Revised: 18 September 2013 – Accepted: 25 October 2013 – Published: 15 November 2013 Abstract. Spatial maps of the finite-time Lyapunov exponent (FTLE) have been used extensively to study LCS in two-dimensional dynamical systems, in particular with application to transport in unsteady fluid flows. We use the time-periodic double-gyre model to compare spatial fields of FTLE and the path-integrated Eulerian Okubo–Weiss parameter (OW). Both fields correlate strongly, and by solving the dynamics of the deformation gradient tensor, a theoretical relationship between both magnitudes has been obtained. While for long integration times more and more FTLE ridges appear that do not seem to coincide with the stable manifold, ridges in the field of path-integrated OW represent fewer additional structures. 1 Introduction The importance of Lagrangian analysis to understanding complex transport problems in fluids has been established during the last decade (see Griffa et al., 2007 and Neufeld and Hernández-García, 2009 and references therein for a review). Lagrangian analysis is directly linked to the dynamical systems approach to transport that analyzes the phase space of the dynamical system driving the trajectories in a flow. This dynamical system can be a fluid flow given by velocity data from different sources as models or measured data, or it can be another dynamical system not related to fluids (Tanaka and Ross, 2009; Kuehl and Chelidze, 2010). Among other techniques, finite-time Lyapunov exponents (FTLE) are used extensively to quantify mixing, and especially to extract persisting transport patterns in the flow, so-called Lagrangian coherent structures (LCS) (Peacock and Dabiri, 2010). The FTLE at a given location xmeasures the maximum stretching rate over the interval τ=t−t0of trajectories starting near the point xat time t0(Shadden et al., 2005). Ridges in the FTLE field are used to estimate finite-time invariant manifolds in the flow that separate dynamically different regions. Repelling (attracting) LCS for τ > 0 (τ < 0) are time-dependent generalizations of the stable (unstable) manifolds of hyperbolic fixed points of the system. These structures govern the stretching and folding mechanisms that control flow mixing. A variety of definitions of LCS in aperiodic flows have been developed and discussed during the last years (Haller, 2000, 2011; Shadden et al., 2005; Branicki and Wiggins, 2010). However, the identification of hyperbolic LCS with the FTLE method exhibits some drawbacks. Especially for long integration times, the FTLE field is noisy, with neighbor initial conditions having substantially different FTLE values, and a variety of structures appear that may be difficult to interpret. The FTLE method is also known to produce false positive ridges in high shear regions (Haller, 2002) that are not generated by hyperbolic structures. This encourages one to study other Lagrangian descriptors and compare them to the FTLE approach. In contrast to the Lagrangian approach, a standard Eulerian technique is the Okubo–Weiss (OW) parameter (Okubo, 1970) that separates strain-dominated regions from vorticitydominated regions in instantaneous velocity fields (see, e.g., Koh and Legras, 2002). It is defined in terms of the eigenvalues of the velocity gradient tensor A=∇v(x). For an incompressible flow they can be written as 2λ±= ±s2 n+s2 s−ω21/2, where sn,ssand ωare the normal and the shear component of the strain, and the relative vorticity PublishedbyCopernicus Publicationsonbehalfofthe EuropeanGeosciencesUnion& theAmerican GeophysicalUnion. 988 V. Pérez-Muñuzuri and F. Huhn: Path-integrated Lagrangian measures of the flow. The sign of the Okubo–Weiss parameter Q=s2 n+s2 s−ω2(1) determineswhether the eigenvalues are real orimaginary. For Q > 0 we have real eigenvalues corresponding to exponential stretching in a hyperbolic region of the flow. For Q < 0 the eigenvalues are imaginary, corresponding to a rotational movement in an elliptic region of the flow. The evolving nature of the time-dependent flow is neglected in this approach, and only instantaneous flow structures can be identified that are not generally relevant to finite-time transport. Although only a vague indicator for structures in time-dependent flows, the Okubo–Weiss (OW) parameter is still widely used, in particular for the detection of mesoscale vortices in the ocean (e.g., Chelton et al., 2007). Here, we aim to study a Lagrangian version of the scalar Eulerian OW parameter by averaging OW along trajectories. QP(τ) =ZQ1/2(`)d`. (2) Several related methods using properties of trajectories have been reported. The mesohyperbolicity criterion of Mezíc et al. (2010) is a Lagrangian measure based on the eigenvalues of the gradient of the flow map. It can be used to divide the flow into mesoelliptic regions (complex eigenvalues) and mesohyperbolic regions (real eigenvalues), and converges to the OW criterion for vanishing integration times. Rypina et al. (2011) propose a Lagrangian measure based on the complexity of trajectories that is related to the fractal dimension of the trajectories. They obtain a scalar field in which structures similar to FTLE-based LCS can be observed. In a similar approach, Jiménez and Mancho (2009) propose that the geometry of transport can be observed in spatial plots of the arclength of trajectories. Similar to the FTLE method, the mentioned methods rely on the availability of a complete spatio-temporal data set of velocity fields. Inthis paper, weperform a simplenumerical experimentin the double-gyre flow and study Eulerian quantities, the OW parameter and its two components, total strain and vorticity, which are integrated along finite time trajectories. We compare the spatial fields of these quantities, FTLE fields and the numerically obtained unstable manifold, and find that the field of the path-integrated OW parameter exhibits ridges that mark the unstable manifold. We propose using the pathintegrated OW for an LCS estimation. 2 Model The periodically varying double-gyre flow was discussed in Shadden et al. (2005). It is used as a standard test case for LCS and can be considered a local view of a gulf stream ocean front. In this case, the flow is described by the stream function, ψ(x,y,t) =Asin(πf (x,t))sin(πy), (3) where f (x,t) =a(t)x2+b(t)x, a(t) =sin ωt, b(t) =1−2sin ωt, over the domain D=[0,2]×[0,1]. For =0 the system can be thought of as a time-independent 2-D Hamiltonian system. For this case there is a heteroclinic connection of the unstable manifold of the fixed point (1,1)with the stable manifold of the fixed point (1,0). For 6=0 the gyres conversely expand and contract periodically in the xdirection such that the rectangle enclosing the gyres remains invariant. The periodic perturbation leads to mixing between the two gyres. The respective advective velocity v(x)is obtained as v=ez×∇ψ. For the numerical experiments, initially a cluster of 600× 300 tracer particles was regularly spaced in the domain D. Then, the Lagrangian trajectories of these particles are calculated by integrating the equations above using a 4th-order Runge–Kutta scheme and a fixed time step of 1t =10−3. In order to characterize the coherent Lagrangian transport of this flow, the finite-time Lyapunov exponent σand the Okubo–Weiss parameter Qwere calculated as indicated in Appendix A. Lagrangian coherent structures (LCS) were extracted following the algorithms described in Appendix B. 3 Results We compare the transport structures obtained by different Lagrangian measures in Fig. 1. The FTLE and QPfields for the double-gyre flow are shown in Fig. 1a, b after integration time τ=3T. Similar Lagrangian coherent structures appear in both patterns. As in previous studies, the observed ridge location of σcorresponds to the stable manifold, which also appears for the QPfield at the same location. The most noticeable observation by visual inspection is that many more ridges are present in the σfield than in the QPfield. Comparing the maps of total strain and vorticity Fig. 1c, d with those of σand QP, a strong relationship has been observed. This correlation is stronger for the strain than for the vorticity. In general, there are very low values of σand QP(similarly for the strain) inside the areas with larger values of vorticity. The initial areas of vorticity located around the position of the gyres are advected by the flow and remain encircled by the coherent structures in places where mixing is very small. The observed Lagrangian coherent structures are stretched and folded as time goes on for both σand QPmaps. Figure 2a, b compare σand QPLCS with the directly computed stable manifold WSin Fig. 2d for the double-gyre model for τ=3T. We use the stable manifold in Fig. 2d as a ground Nonlin. Processes Geophys., 20, 987–991, 2013 www.nonlin-processes-geophys.net/20/987/2013/ V. Pérez-Muñuzuri and F. Huhn: Path-integrated Lagrangian measures 989 (a) FTLE × 104(b) QP (c) Total Strain (d) Vorticity 0.5 1 1.5 2 0.5 1 1.5 1 2 3 1 2 3 Fig. 1. Finite-time Lyapunov exponent σand Okubo–Weiss QP fields (upper left and right images, respectively). The total strain s2=s2 n+s2 sand vorticity ω2fields are shown below (left and right images, respectively). Except for the FTLE field, the rest of the magnitudes were integrated along the Lagrangian path followed by tracer particles during the integration time τ=3T. Set of parameters: A=0.1, =0.2, T=10 and ω=2π/T . truth standard for the LCS in this example. Note that QP ridges approximate the manifold WSat least as well as the σridges. Even if the parameters of the ridge extraction are varied, we observe more disconnected spurious ridges in the σfield than in the QPfield. This is congruent with the observation of more ridges in the σfield in Fig. 1. Figure 2c shows the ridges extracted from the total strain field. This quantity is objective in the sense of frame invariance. The extracted ridges are similar to the ridges of QPand show few spurious ridges as well. Differences between the total strain ridges and the manifold WScan be seen in the center of the left gyre at ∼(0.7,0.7), where a part of the manifold is not represented as a ridge. The main observation here is that many spurious ridges in the σfield do not appear in the QPand total strain fields. The statistical relationship between σand QPis quantified in Fig. 3 in terms of the probability distributions of each magnitude at different values of τ. Note that as time integration increases, σdistribution becomes broader, confirming the wide range of σvalues shown in Fig. 1, and its maximum shifts to higher FTLE values. For τ=4T,σdistribution saturates at larger values as spurious ridges develop. However, QPdistributions remain narrower and the peak (and mean) remains approximately at the same values as τincreases (except for τ=T). QPvalues different from zero concentrate mostly in the coherent structures for any value of τ, which Fig. 2. Ridge extraction from (a) FTLE σ,(b) Okubo–Weiss QP, and (c) total strain fields. The ridge detection algorithm used values larger than 50% of maximum for FTLE and 20% for OW and total strain. (d) The stable manifold WS(x0,y0)of the moving instantaneous saddle point around (x0,y0)=(1,0)after the integration time τ=3T(see Appendix B). Set of parameters as in Fig. 1. are closer to the stable manifold than σones, as was shown in Fig. 2. The pdfs of total strain s2behave similarly to the ones for QP. Figure 4 shows the Pearson correlation parameter rbetween σand QPfields calculated for various values of the integration time τand . For comparison, we show the correlation values between σand the total strain s2=s2 n+s2 sand vorticity ω2,r(σ,s2)and r(σ,ω2), respectively. Note that r(σ,s2) > r(σ,ω2)and r(σ,s2)≈r(σ,QP)for any value of τ. As the integration time τincreases, r(σ,QP)diminishes. Results not shown here demonstrate that correlation values near zero are obtained when the Okubo–Weiss parameter is calculated as an Eulerian magnitude at any instant of time. To shed some light on this different behaviour of r(σ,s2) and r(σ,ω2)with τ, we focus on a theoretical relationship between these magnitudes. Following the notation described in Appendix A, the dynamics of the deformation gradient tensor is determined by dF(t)/dt=A(t)F(t). Starting at some initial time from F(t0)=I, the general form of F can be written formally (Chevillard and Meneveau, 2006; Falkovich et al., 2001) in terms of the time-ordered or pathordered exponential function, F(t) =expP  t Z t0 dtA(t) .(4) The index Pmeans that in a Taylor expansion of the exponential, all matrices are ordered such that later times appear www.nonlin-processes-geophys.net/20/987/2013/ Nonlin. Processes Geophys., 20, 987–991, 2013 990 V. Pérez-Muñuzuri and F. Huhn: Path-integrated Lagrangian measures 0 0.5 1 0 0.05 0.1 (a) σ,QP PDF 0 0.5 1 0 0.05 0.1 (b) σ,QP PDF 0 0.5 1 0 0.05 0.1 (c) σ,QP PDF 0 0.5 1 0 0.05 0.1 (d) σ,QP PDF Fig. 3. Probability distribution functions of σ(blue solid line), QP (red dashed line), and s2(black dotted line) for τ=T ,2T , 3Tand 4T(a–d), respectively. All quantities were normalized to one for comparison and PDFs were normalized to have the same unitary area. Set of parameters as in Fig. 1. on the left. Accordingly, for the Cauchy–Green tensor, C(t) =FT(t)F(t) (5) =expP  t Z t0 dtAT(t) expP  t Z t0 dtA(t) . Then, if the gradient tensor Aremains constant during a time scale ξ≤t−t0, it follows that C(t) =eξATeξA=e2ξS,(6) where S=(AT+A)/2 is the strain-rate tensor. Finally, using the similarity transformations logC=V(logC’)V−1, S=W(S’)W−1, logC’ =diag(logλi),S’ =diag(χi)(λi,χi eigenvalues of Cand S, respectively; V,Wmatrix of eigenvectors of logC’ and S’, respectively) (Horn and Johnson, 1991; Davies and Higham, 2003), we conclude that eig(S)∝log(eig(C)). (7) In other words, FTLE values should be proportional to the eigenvalues of the strain-rate tensor. We see this in our numerical experiment, since the FTLE shows a clearly larger correlation with the total strain than with the vorticity. In the theory above, the crucial step is represented by the Kolmogorov timescale ξ(Chevillard and Meneveau, 2006), during which the gradient tensor Aremains constant. Note in Fig. 4 that as the integration time τ→0, r(σ,s2)→1, as expected according to Eq. (7). The extreme case occurs for the stationary double-gyre flow (i.e. =0 in Eq. 3). For this model, both the vorticity and the total strain fields attain similar distribution values; however, r(σ,ω2)→0 and r(σ,s2)≈1. 0123 0 0.5 1(a) τ/T Correl. Coeff. 0 1 2 3 0 0.5 1(b) τ/T Correl. Coeff. Fig. 4. Correlation coefficients between σand the QPparameter (circles), σand the total strain s2=s2 n+s2 s(squares), and σand the flow vorticity −ω2(rhombs), as a function of the integration time τ.(a) =0.2 and (b) =0.5. Rest of parameters as in Fig. 1. 4 Conclusions Path-integrated Lagrangian measures have been compared to FTLE fields. The Okubo–Weiss parameter integrated along the path followed by particle trajectories QPgives rise to Lagrangian coherent structures similar to those observed in FTLE fields. FTLE distributions are broader than QPdistributions as time integration increases. Comparing these LCS with the directly numerically computed stable manifold, we observe less spurious ridges in QPfields. As a consequence of that, the correlation between FTLE and QPpatterns decreases with increasing time. We observed a better correlation with the path-integrated strain than with the vorticity. A theoretical relationship between FTLE and the strain seems to agree with that observation, mainly for small values of τ. Although integration along the path can lead to clear LCS patterns, these new variables do not necessarily account for the same information as their Eulerian counterparts. Scalar variables are calculated instantaneously and summed up with time, losing the directional information that was present in the original tensor (Lapeyre et al., 1999). Thus, for example, the integrated strain field does not necessarily account for real flow stretching. It should be interesting to assess the effect of this other part on the flow description. Such work is in progress. Appendix A Lagrangian methods In order to define the Lagrangian measures used in this paper, the notation is as follows: a initial regular grid of tracers at initial positions x0(t0)is advected by the velocity field v(x,t) to a final position x(t0+τ) after a finite-time τ,x(t0+τ) = φt0+τ t0(x0(t0)). The gradient of the flow map (the deformation gradient tensor) is F(x0)=∇φt0+τ t0(x0)=∇(x(t0+τ)). FTLE values are then computed from the trajectories of Lagrangian tracers in the flow as σ(x0,t0,τ) =1 τlogpλmax(C(x0)), (A1) Nonlin. Processes Geophys., 20, 987–991, 2013 www.nonlin-processes-geophys.net/20/987/2013/ V. Pérez-Muñuzuri and F. Huhn: Path-integrated Lagrangian measures 991 where λmax is the largest eigenvalue of the right Cauchy– Green deformation tensor C=FTF. By definition Cis a symmetric tensor and has real eigenvalues. √λmax denotes the ratio of stretching between two initially close tracers in the direction of the largest stretching. For the computation of FTLE fields, the integration time τmust be predefined. Basically, the time τhas to be long enough to allow trajectories to explore the Lagrangian coherent structures present in the flow. Appendix B LCS and stable manifolds Following Haller (2002) and Shadden et al. (2005), Lagrangian coherent structures (LCS) can be estimated as local maximizing surfaces of the FTLE. Here we compute LCS in terms of second-derivative ridges of the FTLE field σor the Okubo–Weiss QP, at each time t(Sadlo and Peikert, 2007). In order to filter spurious ridges or weak LCS, a minimum threshold for both fields was used. This threshold was chosen as the value at which the probability distribution function of σand QPhas a maximum. The stable manifold WS(x0,y0)of the moving saddle point around (x0,y0)=(1,0)for the double-gyre model is defined as the set of points that converge to (x0,y0)under forward iteration of Eq. (3) (Hobson, 1993; Mancho et al., 2003). Numerically, the stable manifold is grown by adding one point at each step of the integration near the saddle point and integrating backwards Eq. (3). Acknowledgements. This work was supported in part by the Ministerio de Educación y Ciencia under research grant no. FIS2010-21023 and ICMAT Severo Ochoa project SEV-2011- 0087. We gratefully acknowledge the helpful comments of two anonymous reviewers. Edited by: E. Hernández-García Reviewed by: two anonymous referees References Branicki, M. and Wiggins, S.: Finite-time Lagrangian transport analysis: stable and unstable manifolds of hyperbolic trajectories and finite-time Lyapunov exponents, Nonlin. Processes Geophys., 17, 1–36, doi:10.5194/npg-17-1-2010, 2010. Chelton, D. B., Schlax, M. G., Samelson, R. M., and de Szoeke, M. A.: Global observations of large oceanic eddies, Geohpys. Res. Lett., 34, L15606, doi:10.1029/2007GL030812, 2007. Chevillard, L. and Meneveau, C.: Lagrangian dynamics and statistical geometric structure of turbulence, Phys. Rev. Lett., 97, 174501, doi:10.1103/PhysRevLett.97.174501, 2006. Davies, P. I. and Higham, N. J.: A Schur-Parlett algorithm for computing matrix functions, SIAM J. Matrix Anal. Appl., 25, 464– 485, 2003. Falkovich, G., Gawedzki, K., and Vergassola, M.: Particles and fields in fluid turbulence, Rev. Mod. Phys., 73, 913–975, 2001. Griffa, A., Kirwan, A. D., Mariano, A. J., Ozgokmen, T., and Rossby, T. (Eds.): Lagrangian analysis and prediction of coastal and ocean dynamics, Cambridge Univ. Press, 2007. Haller, G.: Finding finite-time invariant manifolds in twodimensional velocity fields, Chaos, 10, 99–108, 2000. Haller, G.: Lagrangian coherent structures from approximate velocity data, Phys. Fluids, 14, 1851–1861, 2002. Haller, G.: A variational theory of hyperbolic Lagrangian coherent structures, Physica D, 240, 574–598, 2011. Hobson, D.: An efficient method for computing invariant manifolds for planar maps, J. Comput. Phys., 104, 14–22, 1993. Horn, R. A. and Johnson, C. R.: Topics in matrix analysis, Cambridge University Press, Cambridge, UK, 1991. Koh, T. Y. and Legras, B.: Hyperbolic lines and the stratospheric polar vortex, Chaos, 12, 382–394, 2002. Kuehl, J. and Chelidze, D.: Identifying invariant manifold using phase space warping and stochastic interrogation, Int. J. Nonlin. Mech., 45, 42–55, 2010. Lapeyre, G., Klein, P., and Hua, B. L.: Does the tracer gradient vector align with the strain eigenvectors in 2D turbulence?, Phys. Fluids, 11, 3729–3737, 1999. Jiménez Madrid, J. A. and Mancho, A. M.: Distinguished trajectories in time dependent vector fields, Chaos, 19, 013111, doi:10.1063/1.3056050, 2009. Mancho, A. M., Small, D., Wiggins, S., and Ide, K.: Computation of stable and unstable manifolds of hyperbolic trajectories in twodimensional, aperiodically time-dependent vector fields, Physica D, 182, 188–222, 2003. Mezíc, I., Loire, S., Fonoberov, V. A., and Hogan, P.: A new mixing diognostic and gulf oil spill movement, Science, 330, 486–489, 2010. Neufeld, Z. and Hernández-García, E.: Chemical and biological processes in fluid flows, Imperial College Press, 2009. Okubo, A.: Horizontal dispersion of flotable particles in the vicinity of velocity singularities such as convergences, Deep-Sea Res., 17, 445–454, 1970. Peacock, T. and Dabiri, J.: Introduction to Focus Issue: Lagrangian CoherentStructures,Chaos,20,017501,doi:10.1063/1.3278173, 2010. Rypina, I. I., Scott, S. E., Pratt, L. J., and Brown, M. G.: Investigating the connection between complexity of isolated trajectories and Lagrangian coherent structures, Nonlin. Processes Geophys., 18, 977–987, doi:10.5194/npg-18-977-2011, 2011. Sadlo, F. and Peikert, R.: Efficient visualization of Lagrangian coherent structures by filtered AMR ridge extraction, IEEE Trans. Vis. Comput. Gr., 13, 1456–1463, 2007. Shadden, S. C., Lekien, F., and Marsden, J.: Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows, Physica D, 212, 271–304, 2005. Tanaka, M. L. and Ross, S. D.: Separatrices and basins of stability from time series data: an application to biodynamics, Nonlinear Dyn., 58, 1–21, 2009. www.nonlin-processes-geophys.net/20/987/2013/ Nonlin. Processes Geophys., 20, 987–991, 2013