Finite-time barriers to front propagation in two-dimensional fluid flows

Finite-time barriers to front propagation in two-dimensional fluid flows

John R. Mahoney    Kevin A. Mitchell School of Natural Sciences, University of California, Merced, California, 95343
July 6, 2019

Recent theoretical and experimental investigations have demonstrated the role of certain invariant manifolds, termed burning invariant manifolds (BIMs), as one-way dynamical barriers to reaction fronts propagating within a flowing fluid. These barriers form one-dimensional curves in a two-dimensional fluid flow. In prior studies, the fluid velocity field was required to be either time-independent or time-periodic. In the present study, we develop an approach to identify prominent one-way barriers based only on fluid velocity data over a finite time interval, which may have arbitrary time-dependence. We call such a barrier a burning Lagrangian coherent structure (bLCS) in analogy to Lagrangian coherent structures (LCSs) commonly used in passive advection. Our approach is based on the variational formulation of LCSs using curves of stationary “Lagrangian shear”, introduced by Farazmand, Blazevski, and Haller [Physica D 278-279, 44 (2014)] in the context of passive advection. We numerically validate our technique by demonstrating that the bLCS closely tracks the BIM for a time-independent, double-vortex channel flow with an opposing “wind”.

Looking into a pocket near the edge of a lazy stream, one may see froth and small plant particulates twisting and evolving in sinuous strands. While streams are always subject to change—a branch falls in, a rock rolls, a fish swims nearby—the observed pattern changes little due to these differences in the stream’s current. An explanation for this is that the pattern organizes around a robust set of curves that form a “skeleton” for particulate advection.
While it is clear to the eye that such curves exist, it is another thing to define them precisely. A variety of techniques have been developed for this purpose, and while the original techniques were limited to flows that were either constant or periodic in time, the most general modern techniques can analyze flows with arbitrary time-dependence.
Considering now not a stream, but the entire ocean, one can find a teaming algae population within large eddies. The “living pattern” created by the swirling collection of replicating organisms is different from the passive particulates because the algae cluster is active—replication increases the size of the algae patch. There is a corresponding new set of “active” curves underlying the pattern formed by this growing front.
In this paper, we develop a technique to extract the most attracting and repelling curves for such active material within an arbitrary flow defined over a finite time interval. We imagine eventual application to improving management of oceanic algae blooms or to informing the design of combustion systems.

I Introduction

A growing interest has recently developed in fluid systems that combine the phenomenon of (chaotic) advection with that of front propagation. This interest is due in part to geophysical and technological applications, such as the growth of plankton blooms in oceanic flows Neufeld and Hernandez-Garcia (2009); Scotti and Pineda (2007); Sandulescu et al. (2008), the growth of the ozone hole in the atmosphere, and the spreading of chemical reactions in microfluidic devices. Another important driver of this field has been a steady accumulation of experimental observations and measurements of (Belousov-Zhabotinsky) chemical reaction fronts spreading in driven, laboratory-scale flows Paoletti and Solomon (2005, 2005); Paoletti et al. (2006); Schwartz and Solomon (2008); Boehmer and Solomon (2008); O’Malley et al. (2011); Pocheau and Harambat (2006, 2008); von Kameke et al. (2013), as well as early theoretical predictions on such advection-reaction-diffusion systems Abel et al. (2001, 2002); Cencini et al. (2003).

A recent advance in understanding advection-reaction-diffusion (ARD) systems is the realization that invariant manifold theory, familiar in the analysis of passive advective mixing MacKay et al. (1984); Ottino (1989); Rom-Kedar and Wiggins (1990); Wiggins (1992), is also applicable to the analysis of front propagation in fluids Mahoney et al. (2012); Mitchell and Mahoney (2012); Bargteil and Solomon (2012); Mahoney and Mitchell (2013); Mahoney et al. (2015); Gowen and Solomon (2014). Such manifolds have been called burning invariant manifolds, or BIMs, to distinguish them from their analogous, but distinct, advective counterparts. The most important property of BIMs is that they form time-invariant, or time-periodic, one-way barriers to front propagation in fluid flows that are themselves either time-invariant or time-periodic, respectively. The unstable BIMs are attracting structures, in that an initial stimulation near such a BIM ultimately creates a reaction front in the fluid that converges upon the BIM. In the most striking cases, this front persists for arbitrarily long times, forming a frozen or pinned front, whose profile follows the BIM Mahoney et al. (2015). Alternatively, stable BIMs are repelling structures. Initial stimulation points on either side have qualitatively distinct future behavior. Various other aspects of passive invariant manifolds extend to BIMs. For example, BIMs lead to a theory of turnstiles for front propagation in time-periodic fluid flows Mahoney and Mitchell (2013).

Despite the successful application of BIMs to ARD systems, their applicability has been limited to flows that are either time-independent or time-periodic. The present paper addresses this limitation by developing a theory of prominent, one-way barriers to front propagation in general fluid flows specified over a finite interval of time, with no assumption on the time-dependence over this interval. These barriers, called burning Lagrangian coherent structures (bLCSs), are either locally most repelling or most attracting. We validate this theory by numerically demonstrating that the bLCS reproduces the BIMs for time-independent flows considered over a sufficiently long time interval. A future publication will assess this theory in the general case of unsteady flows.

Our work is motivated by recent progress in the use of finite-time Lyapunov exponent (FTLE) fields Pierrehumbert and Yang (1993); von Hardenberg et al. (2000); Doerner et al. (1999); Haller (2001, 2002) and the associated Lagrangian coherent structures (LCSs) for passive advection in unsteady flows Haller (2000); Haller and Yuan (2000); Haller (2001, 2002, 2011); Shadden et al. (2005, 2006). A common approach in advection studies is to assume that ridges of the (forward) FTLE field are (repelling) LCSs that approximately separate trajectories with different future behavior. The FTLE ridge approach has its deficiencies, however, and can produce both false negatives and false positives when identifying LCSs Haller (2011), so it is best used, perhaps, as an intuitive diagnostic tool. That said, if one were to implement the FTLE-ridge approach for front-propagation dynamics, one would immediately encounter a challenge of dimensions. This is because the front-element dynamics in two spatial dimensions is naturally represented as a three-dimensional dynamical system, with and the two spatial coordinates and the orientation of the front element. (See Sect. II.) Note that the objects of our interest, the bLCSs, are still one-dimensional curves, just like the BIMs they seek to generalize; that is, a bLCS is not a codimension-one object in our study. Thus, a possible implementation of the FTLE-ridge approach to the front-element dynamics would not seek to follow 1D ridges in a 2D space, or even 2D ridges in a 3D space, but rather to follow 1D ridges in a 3D space. We see no obvious reason why such curves in the 3D phase space would correspond to fronts in the -position space, since a curve in the 3D space must satisfy the front-compatibility criterion (see Sect. II) for it to be the lift of a front in -space. It remains an open question, then, whether the 1D FTLE ridges would even yield fronts in the -position space, and assuming they did, these would undoubtedly suffer the same deficiencies noted in Ref. Haller (2011) for passive advection.

Instead of seeking FTLE ridges, we base our approach on the recent work of Farazmand, Blazevski, and Haller (FBH) Farazmand et al. (2014), which uses a variational principle to define the LCS. The advantage of a variational formulation is that it is based from the start on a search for curves in -space. These curves could be material lines, as in FBH Farazmand et al. (2014), or they could be fronts, as in the current study. The fact that the evolution of each front element depends on its orientation is naturally accommodated within the variational approach. Following FBH Farazmand et al. (2014), we base the variational integral on the Lagrangian shear, generalized to the front propagation dynamics. The output of the variational approach, when combined with the front compatibility criterion, is a flow on a 2D surface, which we call the shearless surface, embedded within the 3D phase space. This flow foliates the shearless surface into shearless fronts, which are candidates for a bLCS. This reduction to a 2D surface tames the challenge of dimensions discussed above. One can, for example, plot burning FTLE (bFTLE) fields on the 2D surface. However, the dimensional challenges have not been completely eliminated, since the shearless surface can have a complicated geometry within the 3D phase space, with folds and intersections that are not easily represented in a 2D projection. This geometry raises new considerations in the search for bLCSs that do not arise when searching for LCSs in passive advection.

This article has the following structure. Section II reviews our prior dynamical systems approach to front propagation in fluid flows and introduces the concept of BIMs as one-way barriers. Section III contains the core derivations. The main result is Theorem 2 of Sect. III.7, which identifies perfect shearless fronts as the candidates for a bLCS. Sect. III.8 describes how to compute these shearless fronts, and Sect. III.9 formulates the selection of the (repelling) bLCS from the set of shearless fronts as finding that front which maximizes the average normal repulsion. Section IV applies the bLCS approach to a time-independent double-vortex in a channel with wind. We study both short (Sect. IV.2) and long (Sect. IV.3) evolution times to show that the bLCS closely reproduces the BIM of the time-independent flow.

Ii Preliminaries

Following previous studies Mahoney et al. (2012); Mitchell and Mahoney (2012); Bargteil and Solomon (2012); Mahoney and Mitchell (2013); Mahoney et al. (2015); Gowen and Solomon (2014), we give a brief account of front propagation from the dynamical systems perspective. We make the sharp front assumption, in which the front marks a clear delineation between a reacted or “burned” region and an “unburned” region. Then, a front in a two-dimensional fluid flow can be represented as a curve in the three-dimensional -space, where is an arbitrary smooth parameterization. Here denotes the orientation of the normal and tangent to the front, where . (Fig. 1.) The normal vector points in the propagation, or burning, direction of the front.

Figure 1: (color online) Each front element evolves according to Eqs. (2).

Not all curves in -space represent fronts. There must be pointwise agreement between the tangent vector of the curve and the orientation given by . We call this the front-compatibility criterion


The front propagation speed in the local fluid frame, which we call the burning speed, is assumed to be isotropic and homogeneous throughout the fluid. Furthermore, we assume that is independent of the local curvature of the front or any other front property, i.e. is entirely constant in our analysis. Then each front element evolves independently of its neighboring front elements according to (see Fig. 1)


where is the fluid velocity field, dot denotes the time derivative, comma denotes differentiation with respect to , and repeated indices are summed. The fixed points of this dynamical system are called burning fixed points. For time-independent or time-periodic fluid velocity fields , burning fixed points may have any combination of stable (S) and unstable (U) directions: SSS, SSU, SUU, UUU. The one-dimensional stable and unstable manifolds of SUU and SSU burning fixed points, respectively, play a special role and are called burning invariant manifolds (BIMs). BIMs satisfy the front-compatibility criterion and may thus be viewed as (virtual) fronts within the fluid. As such, they form one-way barriers to front propagation, owing to the fact that no front can catch up to and surpass another front moving in the same direction (what we call the no-passing lemma Mitchell and Mahoney (2012).) BIMs are thus natural, dynamically defined, geometric objects that restrict the propagation of fronts in an advecting fluid flow that is either time-independent or time-periodic.

Iii Coherent structures as solutions to an optimization problem

iii.1 Lagrangian shear for fronts

Figure 2: (color online) The definition of Lagrangian shear and normal repulsion . An initial normal displacement of a front element (left) evolves under into a displacement whose tangential projection (blue) is times the Lagrangian shear . The normal projection (blue) is times the normal repulsion.

We adapt the definition of Lagrangian shear in FBH Farazmand et al. (2014) to the front propagation scenario. We first assume that is given between an initial time and a final time , with no requirement that the flow be time-independent or periodic. Then the evolution of a front element over this time interval is obtained by solving Eqs. (2). We denote the flow map acting on the 3D phase space by . Consider now an initial front at time , and suppose the element at a given position is displaced in the normal direction an infinitesimal distance . (The orientation of the front element is not perturbed.) Under the time evolution , the point maps forward to , and the relative displacement between the original and perturbed front elements maps forward to a vector that will in general have a nonzero tangential projection onto the time-evolved front. (See Fig. 2.) We define the Lagrangian shear over the time interval to be this projected length divided by .

To compute , we first note that the 3D normal displacement vector at evolves forward to at , where is the gradient of , with components . Since , shown in Fig. 2, is the -component of , we have




is the restricted projection operator from the -tangent space into the -tangent plane. Note . We then define


to be the restriction of to the -tangent plane. Note also that Eq. (1) implies




Then Eq. (3) becomes


We next need an expression for the tangent vector to the front at . The 3D tangent vector to the curve at time is of the form . This vector maps forward under to a vector whose -component is proportional to . Note that this fact is actually true regardless of the value of , since depends on the curvature of and not its tangent. We may thus set to obtain


where we again used the front-compatibility criterion Eq. (1). We are now in a position to compute via


The minus sign on conforms to the convention of FBH Farazmand et al. (2014). To summarize, we may express the Lagrangian shear as




Note that and are both symmetric matrices. We call the projected Cauchy-Green tensor. Note that even though is a matrix acting on the -tangent plane, it still depends on the orientation coordinate . In the limit that the burning speed vanishes, the projected Cauchy-Green tensor reduces to the usual Cauchy-Green tensor defined for passive advection. Finally, note that we have written either as a function of or , since does not depend on the magnitude of and since is the angle of , as given by Eq. (1).

iii.2 Normal repulsion for fronts

The bLCS will be defined to locally maximize the average normal repulsion . Adapting the definition of normal repulsion in Ref. Haller (2011) to the front propagation scenario, we define to be the normal component of the time-evolved displacement in Fig. 2. That is, the normal repulsion is a local measure of the degree to which a nearby front is pushed away, with indicating repulsion and attraction. To compute , we shall need the normal at the final time , which we claim is given by


To prove this, we need only show that . From Eqs. (10) and (16), is proportional to


Using Eq. (8), the normal repulsion is


In summary, the normal repulsion is a function on the 3D phase space given by


iii.3 Variational problem

The variational problem is now stated as follows. We seek a curve , parameterized as , that makes the average Lagrangian shear


stationary under infinitesimal variations, where the average is taken over the interval , choosing the initial point to have parameter . The endpoints of are assumed fixed under the variations. Note that does not depend on how the curve is parameterized, since the numerator and denominator in Eq. (13) are both second order in . However, the average in Eq. (20) does depend on the parameterization through the line element . Different choices of the parameterization will change the weighting of the average along . The stationarity condition of should be understood not only in terms of variations in but also in terms of variations of the parameterization along , with the interval of integration held constant. (In contrast, one could use the euclidean line element for the average, in which case would be independent of the parameterization of .) We shall find it useful to change the parameterization of the stationary curves, so we next summarize the reparameterization of trajectories in the Hamiltonian and Lagrangian formalisms.

iii.4 Trajectory reparameterization for Hamiltonian systems

Suppose is a particular trajectory for an -independent Hamiltonian and let be the value of the (conserved) Hamiltonian along the trajectory. Here plays the role of time. Define a new trajectory parameterization by


for some phase space function , and denote the reparameterized trajectory by . Then is a trajectory of the Hamiltonian


where the (conserved) value of along the trajectory is .

To prove this fact, we need only demonstrate , where is an arbitrary phase space quantity, denotes the Poisson bracket, and the dot denotes differentiation with respect to . To this end,


where the second to last equality utilizes


and the final equality uses Eq. (21).

iii.5 Trajectory reparameterization for Lagrangian systems

Suppose is a particular trajectory for the -independent Lagrangian . We again let equal the value of the corresponding Hamiltonian along the trajectory. We also define a new parameterization by Eq. (21), except that is expressed as a function of rather than . Then the trajectory , referred to the new parameterization, is a trajectory of the Lagrangian


where and are and re-expressed as functions of


To prove this, we recall the Legendre transforms


where the momentum is


Substituting Eqs. (29) and (30) into Eq. (22) yields Eq. (27).

iii.6 Reparameterizing stationary curves of the average Lagrangian shear

Defining the Lagrangian , Eq. (13) yields


Re-expressing this in terms of the transformed derivative , Eq. (28), yields


which has the same functional form as Eq. (32), since the derivative occurs at the same order (two) in both the numerator and denominator. We now choose


so that Eq. (27) becomes


Equation (31) implies , which is just the derivative of in the direction of increasing radius in the -space. Investigating Eq. (32), however, we see that the dependence of on is solely through the angle of and not . Thus, . Furthermore, since , as well. Thus, Eqs. (29) and (30) imply


Since , the value of along the stationary trajectory also vanishes. These results were obtained by FBH Farazmand et al. (2014) through an alternative derivation.

Theorem 1

Every stationary curve of the average Lagrangian shear , Eq. (20), is also a stationary curve of


where is the (conserved) value of the Lagrangian shear along . The value of along is also conserved and equal to 0. Note that the parameterization of that makes stationary is typically different from the parameterization of that makes stationary.

iii.7 Perfect shearless fronts

Equations (32) and (36) show that is equivalent to . Furthermore, one can show that any curve satisfying is a stationary curve of Eq. (38). Thus, we have the following.

Theorem 2

A curve satisfying , or equivalently , is a stationary curve of


We call such curves perfect shearless fronts (analogous to the perfect shearless barriers of FBH Farazmand et al. (2014)), or just shearless fronts for short. They are the focus of the remainder of this article.

Note that since is symmetric and traceless [see Eq. (15)], it is like a pseudo-Riemannian metric tensor; “pseudo” because it has signature . However, the resulting metric fails to be given purely by a bilinear form, since itself depends on . Such a metric is called a (pseudo-)Finsler metric. Perfect shearless fronts are thus light-like, or null, geodesics of this pseudo-Finsler metric.

iii.8 Computing perfect shearless fronts

In Sects. III.3, III.6, and III.7, fronts were viewed as curves in the -position space for the purpose of the variational problem. Referring now to the full -phase space, a shearless front satisfies both the front compatibility criterion Eq. (1) and . The latter implies that , and hence that is an eigenvector of . One possible approach to finding the shearless fronts, then, is to integrate an eigenvector field of . Let denote a choice of unit eigenvector over -space. Taking proportional to , the front compatibility criterion Eq. (1) becomes




This quantity is zero at if and only if . If the choice of were smooth everywhere, then Eq. (40) would define a smooth two-dimensional constraint manifold in -space. In general, however, there is no continuous choice of available over the entire phase space. The continuity of can break down at points in -space where has degenerate eigenvalues, and hence where there is not a unique eigendirection. Since degeneracy of real symmetric two-by-two matrices is a codimension-two criterion, the degeneracy points generically occur along curves in -space. Even if we remove these degenerate curves from consideration, can not in general be chosen continuously on the remainder of -space. This is due to the fact that may obtain an overall rotation by when transported along a loop encircling a degenerate curve.

Rather than finding the eigenvectors explicitly, an improved approach is to use the constraint equation




This quantity is zero at if and only if is any eigenvector of , not just a particular choice as above. The constraint surface Eq. (42) can thus be understood as the union of the two surfaces defined by Eq. (40) for the two choices of eigenvector . The constraint surface Eq. (42) is a smooth surface without boundary (though potentially with self-intersections), which we call the shearless surface.

A perfect shearless front, lifted into -space via the front compatibility condition, must lie within the shearless surface. A perfect shearless front is thus an integral curve of a vector field tangent to the shearless surface. This vector field, normalized to unity, is




An integral curve of (45) will lie within the shearless surface so long as its initial condition does. Note that these vector fields do not require us to find the eigenvectors of .

We shall also need the vector field that is everywhere orthogonal to Eq. (45), but still tangent to the shearless surface. This field is obtained by replacing by in Eqs. (45) and (46), i.e.,


where parameterizes the curve and


iii.9 Burning Lagrangian coherent structures (bLCSs) as shearless fronts of maximal normal repulsion

The normal repulsion of a front element is given by Eq. (19). On the shearless surface, where is an eigenvector of , this simplifies to


where is the eigenvalue of normal to the tangent. We now define the burning finite-time Lyapunov exponent (bFTLE) field as


Note that the specification of the eigenvalue is based on the orientation of its eigenvector (normal to the front) rather than its magnitude; thus need not be the largest eigenvalue.

Since the bFTLE is defined on the two-dimensional shearless surface, we have overcome part of the challenge of dimensions highlighted in Sect. I. In particular, we can plot and visually interpret images of the bFTLE (Fig. 4a). In some cases a bLCS follows a ridge of the bFTLE field. However, we see no reason why a ridge of the bFTLE field need satisfy the front-compatibility criterion, and so we opt not to compute bFTLE ridges. Again, we see the advantage of using the shearless fronts, for which the problem of defining the bLCS becomes the problem of selecting the “best” curve from the one-parameter family of shearless fronts. We choose the repelling bLCS to be that shearless front that (locally) maximizes the average normal repulsion


where is the -euclidean line element and is the euclidean length of the curve. Similarly, the attracting bLCS is that shearless front that (locally) minimizes the average normal repulsion.

There is still a question of what interval length to choose. In keeping with the intuitive idea of following an bFTLE ridge, we shall adopt the approach of Refs. Haller (2011); Farazmand and Haller (2012) and restrict the length of a shearless front to those regions of the shearless surface where the second derivative of the bFTLE field in the direction normal to the front, given by Eqs. (47), is negative. In practice, the exact length used should make little difference, so we use the curvature criterion to tell us approximately where to cut the shearless curves.

iii.10 Cusps in the shearless fronts

The shearless fronts, and hence a bLCS, may exhibit cusps when plotted in -space. This is one consequence of the 3D phase-space geometry that does not occur for traditional LCSs in passive advection. Cusps have already been recognized as important features of BIMs for time-independent and time-periodic flows (see Fig. 3a); cusps change the one-way barrier direction of the BIMs. Here, we determine when a cusp occurs along a shearless front. First, we consider the fold of the shearless surface under the projection from the -space to the -space; a fold is where two branches of this projection meet. Since a shearless front is smooth on the shearless surface, if it remains on one branch of this surface, then its projection to -space is also smooth. Thus, a cusp of a shearless front can only occur at a fold of the shearless surface. What is less obvious, but nevertheless true, is that every time a shearless front passes around a fold, its -projection forms a cusp. To prove this, consider a shearless front as it passes around a fold. Locally, can be expressed as a function of , which by Eqs. (45) satisfies


At the fold of the shearless surface, the 3D normal vector to the surface, equal to , has no component in the direction. That is, . This implies that at a fold, by Eq. (46), and hence , by Eq. (52). More precisely, one sees that, at a fold, transitions from pointing in one direction to pointing in the opposite direction. This is exactly the criterion for forming a cusp—the curve moves forward and then backs up. In summary, then, a cusp occurs at exactly those points where a shearless front passes transversely through a fold in the shearless surface.

Iv A test case

iv.1 Double-vortex in a windy channel

Figure 3: (color online) Double vortex channel flow with vortex strength and burning speed . a) no wind, ; b) . The streamlines of the flow are illustrated by grey curves. The grey regions are “slow zones”, where the fluid speed is less than . The blue and red curves are stable and unstable BIMs respectively. The BIMs begin at burning fixed points on the channel boundary. For plot a, the BIMs have cusps where they intersect the central slow zones. For plot b, the BIMs span the channel without forming a cusp. The arrows normal to the BIMs illustrate the burning direction.

As an example, consider an infinitely long channel flow containing two side-by-side, counter-rotating vortices plus a spatially uniform “wind” flowing from right to left. The streamfunction is


where is the vortex strength and is the wind speed (Fig. 3). Throughout our analysis we use , , and burning speed . The channel width spans to . We choose this flow because it models what can be realized in the laboratory of Tom Solomon Paoletti and Solomon (2005, 2005); Paoletti et al. (2006); Schwartz and Solomon (2008); Boehmer and Solomon (2008); O’Malley et al. (2011). Furthermore, the BIM spans the channel without a cusp, which simplifies our analysis. Though in this paper we keep the wind speed constant, a future study will consider time-varying wind. Importantly, such time-varying wind can be realized in the Solomon lab Gowen and Solomon (2014).

For a steady wind, there are exactly two hyperbolic advective fixed points of the flow, assuming the wind is not too large. In Figs. 3a and 3b, these fixed points lie within the gray “semicircles” attached to the top and bottom channel wall, near the channel midpoint. The gray regions are “slow zones”, regions where the fluid speed is less than the burning speed . The bottom fixed point in these figures generates an unstable manifold transverse to the channel wall, while the top point generates a stable manifold. The burning dynamics induces the bottom advective fixed point to split into two SSU burning fixed points, one at each intersection point between the boundary of the slow zone and the channel wall. (Only the rightmost burning fixed point is shown in Fig. 3.) There are also two SUU points near the bottom of the channel that do not concern us here. The right-burning SSU fixed point generates an unstable BIM, shown in red in Fig. 3. Similarly, the advective fixed point at the top of the channel splits into two SUU burning fixed points on the upper channel boundary. The leftmost of these is shown in Fig. 3, along with the attached stable BIM in blue.

For weak or no wind, each BIM forms a cusp before crossing the channel (Fig. 3a.) As the wind speed increases above , each BIM crosses the channel without forming a cusp (Fig. 3b.) The process by which the BIM attaches to the boundary is described in detail in Ref. Mahoney et al. (2015). For , each BIM divides the channel in two. If the entire region left of the unstable BIM in Fig. 3b is burned, then it will remain burned for all time, forming a frozen front along the BIM Mahoney et al. (2015); Gowen and Solomon (2014). In contrast, for , any initially burned region will eventually propagate arbitrarily far down the channel to both the left and the right. Thus, the existence of an unstable BIM that traverses the channel with no cusp has a very clear experimental signature.

A stable BIM traversing the channel without a cusp (Fig. 3b) forms a kind of basin boundary. An initial point stimulation left of the stable BIM will eventually grow to become a frozen front. However, a stimulation right of the stable BIM will eventually fill up the entire channel. In the remainder, we shall focus on the stable BIM.

iv.2 The bLCS for a steady wind—short integration time

Figure 4: (color online) a) The two-dimensional shearless () surface embedded in the three-dimensional phase space. The coloring represents the two-dimensional bFTLE field on this surface. The Cauchy-Green integration time is . The black curve is the stable BIM attached to the rightward propagating burning fixed point at the top of the channel. The computation of the isosurface is restricted to the computation box shown. b) A view of the bFTLE field on the shearless surface projected onto the -plane. The two white patches on the bottom left and right are where the shearless surface has exited the computation box. c) The sign of the bFTLE curvature normal to the shearless fronts (see Eq. (47). Blue is negative and red is positive.

We now apply the method of Sect. III to the steady wind case () with the relatively short Cauchy-Green integration time . We should find a bLCS that follows the initial length of the time-independent BIM. Before finding the bLCS, however, it is helpful to visualize the shearless surface. Fig. 4a shows the surface in 3D. For comparison, the stable BIM is shown in black. Near the burning fixed point, the BIM lies near (though not precisely within) the shearless surface. This is an important first check on the validity of our technique; if the BIM were far from the shearless surface, we would have no chance of recovering it.

At the bottom center of Fig. 4a is a “pitchfork” intersection between two branches of the surface. The differentiability of the field implies that the intersection curve between these two branches must be a shearless front, i.e. a trajectory of Eqs. (45). Thus, no shearless fronts may pass through this intersection.

The coloring in Figs. 4a and 4b represents the bFTLE field Eq. (50). Notice that the BIM follows a bFTLE ridge at the top of Fig. 4b, another critical check on our technique. Curiously, there is a second bFTLE ridge on the left in Fig. 4b. This ridge lies on a separate branch of the pitchfork structure, as seen in Fig. 4a.

We next identify those regions where the second derivative of the bFTLE field normal to the shearless curves is negative. Figure 4c illustrates the sign of this curvature, with blue representing negative and red representing positive. There are two regions of negative curvature, one for each of the bFTLE ridges. We focus on the rightmost region associated with the BIM.

We next compute shearless fronts, targeting the bFTLE ridge on the right. We select initial conditions using the intersection of the shearless surface with a plane normal to the BIM, yielding the bold line of initial points in Fig. 5Note1 (). We integrate Eqs. (45) upward toward the top of the channel and downward toward the bottom. Figures 5a and 5b show these shearless fronts in 3D and 2D, respectively.

We then compute the average normal repulsion, Eq. (51), along the segments of the shearless curves above the dashed line in Fig. 5b. This line is chosen to lie close to the bottom of the negative curvature region in Fig. 4c. The resulting average normal repulsion is shown in the inset of Fig. 5c. It exhibits a clear local maximum, which selects the red shearless curve isolated in Fig. 5c. This curve is the bLCS, and it faithfully follows the initial length of the BIM, as anticipated, before eventually deviating markedly from it.

Figure 5: (color online) a) Shearless fronts computed for Cauchy-Green integration time . The initial points (bold) of the trajectories were computed by intersecting a plane transverse to the BIM (at the red dot) with the shearless surface. The BIM is the thin black curve. b) The shearless fronts are shown projected onto the -plane. c-inset) Plot of the average normal repulsion, Eq. (51), as a function of the -coordinate of the initial point for each trajectory in parts a and b. A smooth maximum is marked in red and determines which shearless front is the bLCS. c) The bLCS is plotted in red next to the BIM in black. The open circle marks the initial condition for generating the bLCS.

iv.3 The bLCS for a steady wind—long integration time

We extract a bLCS for the same dynamics as in Sect. IV.2 but for a longer Cauchy-Green integration time of . The result is summarized in Fig. 6. Clearly the bLCS (red) closely tracks the BIM (black), considerably longer than for the case. In fact, the agreement is quite good all the way to the bottom of the channel. This is exactly what we expect to see as the integration time is increased.

Figure 6: (color online) The bLCS (red) and the BIM (black) for Cauchy-Green integration time . The open circle on the bLCS marks the initial point used to integrate the shearless front. (The neighboring shearless fronts are not shown.) The inset shows the average repulsion on a shearless front as a function of the -coordinate of its initial point. The average was computed along the shearless fronts above the dashed line at .

V Conclusion

We have developed a technique to extract maximally repelling (or attracting) fronts, called burning Lagrangian coherent structures (bLCSs) for the finite-time dynamics of fronts evolving within a flowing fluid. We verified that this technique returns the burning invariant manifolds (BIMs) for the case of a time-independent flow. Clearly more work is needed to establish the relevance of these bLCSs to time-periodic and, most importantly, to time-aperiodic flows. This will be pursued in a separate publication.

This work was supported by the US National Science Foundation under grant CMMI-1201236. The authors gratefully acknowledge discussion with Tom Solomon, George Haller, and Mohammad Farazmand.


  • Neufeld and Hernandez-Garcia (2009) Z. Neufeld and E. Hernandez-Garcia, Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach (Imperial College Press, 2009).
  • Scotti and Pineda (2007) A. Scotti and J. Pineda, “Plankton accumulation and transport in propagating nonlinear internal fronts,” J. Marine Res., 65, 117–145 (2007).
  • Sandulescu et al. (2008) M. Sandulescu, C. Lopez, E. Hernandez-Garcia,  and U. Feudel, “Biological activity in the wake of an island close to a coastal upwelling,” Ecological Complexity, 5, 228–237 (2008).
  • Paoletti and Solomon (2005) M. S. Paoletti and T. H. Solomon, “Experimental studies of front propagation and mode-locking in an advection-reaction-diffusion system,” Euro. Phys. Lett., 69, 819 (2005a).
  • Paoletti and Solomon (2005) M. S. Paoletti and T. H. Solomon, “Front propagation and mode-locking in an advection-reaction-diffusion system,” Phys. Rev. E, 72, 046204 (2005b).
  • Paoletti et al. (2006) M. S. Paoletti, C. R. Nugent,  and T. H. Solomon, ‘‘Synchronization of oscillating reactions in an extended fluid system,” Phys. Rev. Lett., 96, 124101 (2006).
  • Schwartz and Solomon (2008) M. E. Schwartz and T. H. Solomon, “Chemical reaction fronts in ordered and disordered cellular flows with opposing winds,” Phys. Rev. Lett., 100, 028302 (2008).
  • Boehmer and Solomon (2008) J. R. Boehmer and T. H. Solomon, “Fronts and trigger wave patterns in an array of oscillating vortices,” Euro. Phys. Lett., 83, 58002 (2008).
  • O’Malley et al. (2011) G. M. O’Malley, M. S. Paoletti, M. E. Schwartz,  and T. H. Solomon, “Pinning and mode-locking of reaction fronts by vortices,” Communications in Nonlinear Science and Numerical Simulation, 16, 4558 – 4563 (2011), ISSN 1007-5704.
  • Pocheau and Harambat (2006) A. Pocheau and F. Harambat, “Effective front propagation in steady cellular flows: A least time criterion,” Phys. Rev. E, 73, 065304 (2006).
  • Pocheau and Harambat (2008) A. Pocheau and F. Harambat, “Front propagation in a laminar cellular flow: Shapes, velocities, and least time criterion,” Phys. Rev. E, 77, 036304 (2008).
  • von Kameke et al. (2013) A. von Kameke, F. Huhn, A. Muñuzuri,  and V. Pérez-Muñuzuri, “Measurement of large spiral and target waves in chemical reaction-diffusion-advection systems: Turbulent diffusion enhances pattern formation,” Phys. Rev. Lett., 110, 088302 (2013).
  • Abel et al. (2001) M. Abel, A. Celani, D. Vergni,  and A. Vulpiani, “Front propagation in laminar flows,” Phys. Rev. E, 64, 046307 (2001).
  • Abel et al. (2002) M. Abel, M. Cencini, D. Vergni,  and A. Vulpiani, “Front speed enhancement in cellular flows,” Chaos, 12, 481–488 (2002), ISSN 10541500.
  • Cencini et al. (2003) M. Cencini, A. Torcini, D. Vergni,  and A. Vulpiani, “Thin front propagation in steady and unsteady cellular flows,” Phys. Fluids, 15, 679–688 (2003), ISSN 10706631.
  • MacKay et al. (1984) R. S. MacKay, J. D. Meiss,  and I. C. Percival, “Transport in hamiltonian systems,” Physica D, 13, 55 (1984).
  • Ottino (1989) J. M. Ottino, The kinematics of mixing: stretching, chaos, and transport (Cambridge University Press, United Kingdom, 1989).
  • Rom-Kedar and Wiggins (1990) V. Rom-Kedar and S. Wiggins, “Transport in two-dimensional maps,” Archive for Rational Mechanics and Analysis, 109, 239 (1990).
  • Wiggins (1992) S. Wiggins, Chaotic Transport in Dynamical Systems (Springer-Verlag, New York, 1992).
  • Mahoney et al. (2012) J. Mahoney, D. Bargteil, M. Kingsbury, K. Mitchell,  and T. Solomon, “Invariant barriers to reactive front propagation in fluid flows,” Euro. Phys. Lett., 98, 44005 (2012).
  • Mitchell and Mahoney (2012) K. A. Mitchell and J. Mahoney, “Invariant manifolds and the geometry of front propagation in fluid flows,” Chaos, 22, 037104 (2012).
  • Bargteil and Solomon (2012) D. Bargteil and T. Solomon, “Barriers to front propagation in ordered and disordered vortex flows,” Chaos: An Interdisciplinary Journal of Nonlinear Science, 22, 037103 (2012).
  • Mahoney and Mitchell (2013) J. R. Mahoney and K. A. Mitchell, “A turnstile mechanism for fronts propagating in fluid flows,” Chaos: An Interdisciplinary Journal of Nonlinear Science, 23, 043106 (2013).
  • Mahoney et al. (2015) J. R. Mahoney, J. Li, K. A. Mitchell, C. Boyer,  and T. Solomon, “Frozen reaction fronts in steady flows: a burning-invariant-manifold perspective,” submitted (2015).
  • Gowen and Solomon (2014) S. Gowen and T. Solomon, “Experimental studies of coherent structures in an advection-reaction-diffusion system,” submitted to Chaos (2014).
  • Pierrehumbert and Yang (1993) Pierrehumbert and Yang, “Global chaotic mixing on isentropic surfaces,” J. Atmos. Sci., 50 (1993).
  • von Hardenberg et al. (2000) von Hardenberg, Lunkeit,  and Provenzale, “Transient chaotic mixing during a baroclinic life cycle,” Chaos, 10 (2000).
  • Doerner et al. (1999) Doerner, Hubinger, Martienssen, Grossmann,  and Thomae, “Stable manifolds and predictability of dynamical systems,” Chaos Solitons Fractals, 10 (1999).
  • Haller (2001) G. Haller, “Distinguished material surfaces and coherent structures in three-dimensional fluid flows,” Physica D, 149, 248 – 277 (2001), ISSN 0167-2789.
  • Haller (2002) G. Haller, “Lagrangian coherent structures from approximate velocity data,” Phys. Fluids, 14, 1851–1861 (2002), ISSN 10706631.
  • Haller (2000) G. Haller, “Finding finite-time invariant manifolds in two-dimensional velocity fields,” Chaos, 10, 99–108 (2000), ISSN 10541500.
  • Haller and Yuan (2000) G. Haller and G. Yuan, “Lagrangian coherent structures and mixing in two-dimensional turbulence,” Physica D, 147, 352 – 370 (2000), ISSN 0167-2789.
  • Haller (2011) G. Haller, “A variational theory of hyperbolic lagrangian coherent structures,” Physica D: Nonlinear Phenomena, 240, 574 – 598 (2011), ISSN 0167-2789.
  • Shadden et al. (2005) S. C. Shadden, F. Lekien,  and J. E. Marsden, “Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimensional aperiodic flows,” Physica D, 212, 271 – 304 (2005), ISSN 0167-2789.
  • Shadden et al. (2006) S. C. Shadden, J. O. Dabiri,  and J. E. Marsden, “Lagrangian analysis of fluid transport in empirical vortex ring flows,” Physics of Fluids, 18, 047105 (2006).
  • Farazmand et al. (2014) M. Farazmand, D. Blazevski,  and G. Haller, “Shearless transport barriers in unsteady two-dimensional flows and maps,” Physica D: Nonlinear Phenomena, 278-279, 44 – 57 (2014), ISSN 0167-2789.
  • Farazmand and Haller (2012) M. Farazmand and G. Haller, “Computing lagrangian coherent structures from their variational theory,” Chaos: An Interdisciplinary Journal of Nonlinear Science, 22, 013128 (2012).
  • (38) To generate initial conditions for the shearless fronts, we first choose a point along the BIM. Here, we transect the BIM with a rectangle whose normal is . For the case, the rectangle has dimensions 0.1 in and in . We compute on this rectangle in a grid. To find a contour of points where , we first use MATLAB’s contourc function. Because of the large curvature found in , we then improve the accuracy of these points with Newton’s method (fsolve).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description