# Chemical front propagation in periodic flows:

FKPP vs G

###### Abstract

We investigate the influence of steady periodic flows on the propagation of chemical fronts in an infinite channel domain. We focus on the sharp front arising in Fisher–Kolmogorov–Petrovskii–Piskunov (FKPP) type models in the limit of small molecular diffusivity and fast reaction (large Péclet and Damköhler numbers, and ) and on its heuristic approximation by the G equation. We introduce a variational formulation that expresses the two front speeds in terms of periodic trajectories minimizing the time of travel across the period of the flow, under a constraint that differs between the FKPP and G equations. This formulation shows that the FKPP front speed is greater than or equal to the G equation front speed. We study the two front speeds for a class of cellular vortex flows used in experiments. Using a numerical implementation of the variational formulation, we show that the differences between the two front speeds are modest for a broad range of parameters. However, large differences appear when a strong mean flow opposes front propagation; in particular, we identify a range of parameters for which FKPP fronts can propagate against the flow while G fronts cannot. We verify our computations against closed-form expressions derived for and for .

ront propagation, large deviations, WKB, cellular flows, Hamilton–Jacobi, homogenisation, variational principles

76V05, 76R99, 35K57

## 1 Introduction

A classical model for the concentration of spreading reacting chemicals is the FKPP, or FK for short, equation named after the classical work by Fisher [13] and Kolmogorov, Petrovskii and Piskunov [20] based on logistic growth and diffusion. Numerous environmental and engineering applications, from the dynamics of ocean plankton to combustion [35, 27], motivate its extension to include the effect of an incompressible background steady flow . The FK equation considered here then takes the non-dimensional form

(FK) |

The reaction term or, more generally, any function that satisfies with for , for and . The non-dimensional parameters are the Péclet and Damköhler numbers

(1) |

where and are the characteristic speed and lengthscale of the flow, the molecular diffusivity, and the reaction time. Motivated by experiments, we focus on two-dimensional channel domains with parallel, impenetrable walls where and take the front-like initial and boundary conditions

where denotes the indicator function. In the absence of advection, (FK) admits front solutions that propagate from the left to the right of the channel at the ‘bare’ speed

(2) |

When the flow is spatially periodic, front solutions persist as pulsating fronts [39, 40, 3], changing periodically in time as they travel at a speed , so that

(3) |

where is the spatial period of the flow.

When advection and reaction dominate over molecular diffusivity, i.e. when

(4) |

the front interface is sharp and can be approximated by a single curve (in 2D as assumed here) where all the reaction takes place. In these conditions, a heuristic model is often used in place of (FK). In this model, the front is the zero-level curve , say, where satisfies the Hamilton–Jacobi equation

(G) |

termed G equation [38] (see also [33, 18]). This model is popular in the combustion science literature (e.g. [28] and references therein). For , the front speed predicted by (G) is obviously , matching the speed predicted by (FK). For spatially periodic , (G) predicts pulsating front solutions propagating with a speed that in general differs from [41, 4]. The relation between the two speeds and is the subject of this paper.

Majda and Souganidis [24] showed that in the limit (4) the leading-order can be deduced from the long-time solution of a certain Hamilton–Jacobi equation. This long-time solution is obtained by applying the asymptotic procedure of homogenisation [21, 11] which exploits spatial scale separation to express in terms of the eigenvalue of a nonlinear cell problem posed over a single period of the flow. A similar procedure can be applied to (G), leading to a different nonlinear eigenvalue cell problem for . The two nonlinear cell problems are significantly simplified for the special case of shear flows [8]. For more general flows explicit analytical expressions are not available and the two cell problems need to be solved numerically. However, these computations can be rather challenging (see e.g. [19] for the nonlinear cell problem related to ).

In this paper, we take an alternative approach that relies on the variational representation of the two front speeds and . For (FK), this approach was introduced by Freidlin and collaborators (see [16, Ch. 10], [14, Ch. 6] and [15]) to establish an expression for in terms of a single trajectory that minimises an action functional. This was subsequently exploited in [36] to obtain explicit results for cellular flows by carrying out a minimisation over periodic trajectories. For (G), Fermat’s principle in a moving medium determines . The variational formulations enable us to express and in closely related terms, and hence to understand and compute their difference for a large class of steady periodic .

We begin with the simple case of shear flows before examining in detail a two-parameter family of periodic cellular flows, given by with streamfunction

(5) |

This is used as a testbed in numerous experimental studies of advection–diffusion–reaction (e.g., [30, 32, 2, 26, 22]). The classic cellular flow introduced in [31] corresponds to a zero mean velocity and to . When confined between walls at and , this flow consists of a one-dimensional infinite array of periodic cells composed of two vortices of opposite circulation. These vortices are bounded by the separatrix streamline that connects a network of hyperbolic stagnation points (see Fig. 1(a)). All streamlines remain closed when and but the symmetry is broken. For , the number of hyperbolic stagnation points doubles and the periodic cell consists of four vortices rotating in alternatively clockwise and anticlockwise directions (see Fig. 1(b)). The topology of the streamlines changes drastically for a non-zero mean velocity : an open channel, bounded by the separatrices and , traverses the domain, splitting apart the row of closed vortices. As the value of increases, the width of the open channel increases (see Fig. 1(c) for and Fig. 1(d) for ). For large enough, the hyperbolic stagnation points and closed streamlines disappear.

Our aim is to determine the effect of flow structures on the value of the two front speeds and and on their difference. To achieve this, we develop and implement a highly accurate numerical method that is based on the efficient discretisation of a pair of variational principles that we obtain. Computations of the two front speeds are complemented by a set of explicit expressions derived in the limit of small and large values of and various values of and .

The paper is organised as follows. In section 2, we provide a brief derivation of the two nonlinear cell problems that determine and . In section 3, we introduce an alternative characterisation in the form of a pair of variational principles from which we deduce that . The two principles greatly simplify for shear flows in which case . Section 4 is devoted to flows with streamfunction (5). The numerical scheme employed for the computations is described in the Appendix. The paper ends with a discussion in section 5.

## 2 Front speed

### 2.1 Equation (Fk)

Gärtner and Freidlin [17] showed that for initial conditions sufficiently close to a step function, the speed of the front associated with (FK) can be deduced by the long-time behaviour of the solution near the front’s leading edge. There and so that (FK) becomes

(6) |

For and , the solution can be sought in the WKBJ (Wentzel–Kramers–Brillouin–Jeffreys) or geometric-optics form

(7) |

Collecting the terms with the same powers in , we find that at leading order satisfies the Hamilton–Jacobi equation

(8) |

the Hamiltonian. The step-function initial conditions correspond to for and for , and the boundary conditions to at , . The front is then identified as the location where (7) neither grows nor decays exponentially with time. It is therefore the level curve

(9) |

In the long-time limit, the solution to (8) converges to that of the homogenised Hamilton–Jacobi equation

(10) |

The effective Hamiltonian, , may be derived from a nonlinear eigenvalue problem, obtained by writing the solution to (8) as the multiscale expansion

(11) |

Here is the slow variable describing the speed of a moving frame of reference and is the fast variable.
We emphasise
the particular form of (11), with a leading-order term
that is independent of
and involves
that depends on only.^{1}^{1}1Note that
may be interpreted as
the Freidlin–Wentzell (small-noise, large-) large-deviation rate function
for the position of fluid particles that have been displaced
by advection and diffusion to a distance in a time
(see [16],[14, Ch. 6] and [15] for rigorous treatments).
The next order involves where
while the boundary conditions at , imply that there, .
Substituting (11) into (8) and equating powers of yields at leading order
the nonlinear eigenvalue problem

(12) |

with the prime denoting derivative with respect to the first argument, can be treated as a parameter and

(13) |

is the eigenvalue. It can be shown that is unique, non-negative, real and convex in (see [21, 9] for proofs) and therefore and are related via a Legendre transform

(14) |

### 2.2 Equation (G)

The long-time solution to equation (G) can be treated similarly. It satisfies the homogenised Hamilton–Jacobi equation

(17) |

with an effective Hamiltonian found as eigenvalue of the nonlinear cell problem

(18) |

where

(19) |

Note that the nonlinearity in is replaced here by . Nevertheless, is unique and convex (details and proofs can be found in [41, 5]). The solution of (17) is then where and and are related via a Legendre transform analogous to (14). Since the front corresponds to , in the long-time limit, the speed of right-propagating (G) fronts is found as the positive solution of or, equivalently, as

(20) |

### 2.3 Comparison

A comparison between expressions (16) and (20) is not possible until we specify and solve the corresponding nonlinear cell problems. These problems are significantly simplified for the special case of shear flows [8]. For more general flows, their solution are pursued using numerical methods. A number of approaches have been proposed (see e.g. [19] which focusses on (16)). However, the computational effort can become high. We now obtain an alternative formulation that sheds light on the difference between the two speeds, is amenable to straightforward numerical computations, and yields explicit expressions in asymptotic limits.

## 3 Variational principles

### 3.1 Equation (Fk)

It is well known (see e.g. [10]) that the solution to (8) may be written as a variational principle involving an action functional associated with the Lagrangian

(21) |

that is dual to the Hamiltonian in (8). For the solution is given by

(22) | ||||

(23) |

where represents a family of smooth trajectories with . From (11) we have

(24) |

where the dependence on the specific value of drops out (e.g. [29]). Together with (22) this determines the function .

Expression (24) can be simplified using the spatial periodicity of the background velocity [36]. Assuming that the minimising trajectory inherits the same periodic structure, we take with and to reduce (24) to

(25) | ||||

Expression (25) provides a direct way to compute the minimising trajectory and, from (15), the corresponding front speed , both numerically and in asymptotic limits. Such computations were carried out in [36] for the specific case of the cellular flow with closed streamlines that we consider further in Section 4. These computations were validated against the numerical evaluation of for finite Péclet and Damköhler numbers obtained from an advection–diffusion eigenvalue problem and direct numerical simulations of (FK) with .

We now obtain an alternative variational characterisation of . Since satisfies , it can be written as extremum of the function

(26) |

for arbitrary variations of the Lagrange multiplier . Here we use that is convex in , so that a single satisfies the constraint enforced by . Using (25) and redefining to absorb a factor , we can rewrite this as

(27) | ||||

(28) |

This can be interpreted as the maximisation of under a constraint enforced by the Lagrange multiplier . Therefore, the front speed predicted by (FK) for , is given as

(29) | ||||

This variational characterisation expresses as the maximum mean velocity achievable by periodic trajectories that are constrained to depart from passive-particle trajectories in a prescribed way.

### 3.2 Equation (G)

An analogous variational characterisation describes the front speed associated with (G). Taking the same initial conditions as for (FK), the front propagates from its initial location at along trajectories that obey Fermat’s principle in a moving medium (e.g. [6]). Thus the front reaches location after a travel time

(30) | ||||

where again we assume that represents a family of smooth trajectories with . In the long-time limit, is large and the front moves at a constant speed given by

(31) |

where once more the dependence on drops out. This characterisation is significantly simplified if we apply the same strategy as before and assume that the minimising trajectory is periodic. Taking with , we obtain that

(32) | ||||

This characterisation of the front speed for (G) closely parallels the characterisation (29) of the front speed for (FK).

For practical computations, it is convenient to rewrite (32) taking as the independent variable, using

(33) |

where denotes the time it takes to reach the point . The minimal travel time over a spatial period is then expressed as

(34) | ||||

and , are taken to be smooth.

### 3.3 Comparison

We now compare the two variational characterisations (29) and (32) for the (FK) and (G) equations. In both the front speeds are expressed in terms of the travel times and which are determined by the periodic trajectories that traverse a spatial period of the flow in the least time. The only difference is that the pointwise constraint on the relative velocity in (32) is replaced by a slacker, time-averaged constraint in (29). An immediate consequence is that

(35) |

While (29) and (32) are useful for comparisons of this type, for numerical computations we found it convenient to use (25) and (34) instead. Eq. (25) is useful for (FK) when, as is the case in section 4, we are interested in computing for a range of values of : the simple dependence of on means that the condition gives an explicit variational formula for as a function of with the endpoint condition as sole constraint.

The variational characterisation (32) is also useful to establish a necessary condition for the existence of right-propagating front solutions for the (G) equation. It is easy to see from the constraint in (32) that

(36) |

For smaller , there are no right-propagating (G) fronts. From (35) we then expect that, for a range of , there exist right-propagating fronts for (FK) but not for (G). We provide explicit examples confirming this in section 4.3.

#### Shear flows

It is easy to show that for shear flows with velocity , . For (FK), the Euler–Lagrange equations associated with the functional in (25) can be written as

(37) |

where and are two constants. The minimum of the functional is then achieved when , where is a constant to be determined. It follows that as imposed by the endpoint condition. The functional then reduces to . Its minimum is non zero for , the maximum velocity in the channel, and given by with such that . Thus,

(38) |

and solving (15) gives the front speed .

On the other hand, the pointwise constraint (32) of the velocity may be parameterised so that

(39) |

where has the same period as . The minimum value of is obtained by maximising . This is achieved for , and , i.e. for trajectories that follow the (straight) streamline associated with maximal flow velocity. We deduce that

(40) |

We therefore conclude that (FK) and (G) are equivalent in describing the long-time speed of propagation. This was previously argued to be the case in [1] and can be inferred from the analysis in [8]; the variational principles (29) and (32) show that this is true in a straightforward way. It is clear that a right-propagating front is obtained for both (FK) and (G) provided that , and that the front is stationary for .

## 4 Front speeds for periodic flows

For more general flows, closed-form formulas are not available. We use the variational problems (25) and (34) whose solutions are easy to approximate numerically. We obtain numerical approximations by discretising trajectories, action functional and constraints and determining the optimal solutions by minimisation. The numerical procedure is detailed in Appendix A. We use this procedure to compute the front speeds for (FK) and (G) and a range of two-dimensional periodic flows. We now describe the results.

### 4.1 Cellular flow

We first compute the solutions for the closed cellular flow with streamfunction (5) and . Figure 2 shows characteristic examples of minimising trajectories obtained for three different values of . For large values of , the periodic trajectories for (FK) and (G) are close to the straight line . In this case, the two trajectories are practically indistinguishable. A larger difference is obtained for small values of , in which case both trajectories follow closely a streamline near the separatrix . In all cases it is clear that the trajectories are invariant under the transformations and .

Figure 3 shows the behaviour of the front speeds for (FK) and (G) as a function of . Clearly, there is a difference between and which is more marked for smaller values of . However, this difference is small: (G) only slightly underpredicts the front speed of (FK). The behaviour of and and their difference can be captured by explicit expressions obtained in two asymptotic limits.

#### 4.1.1 Small- asymptotics

The first asymptotic limit corresponds to . We find an approximation to in this limit by approximating in (25) for . We previously found [36] that the minimising periodic trajectory in (25) may be divided into two regions that we now describe. In region I, and therefore we may seek a regular expansion in powers of of the form

(41) |

where, without loss of generality, we take . In region II, and so we take

(42) |

where with . We then exploit the symmetries that characterises the streamfunction to extend the trajectory over the whole time period .

Substituting (41) and (42) into (25) gives a sequence of integrals corresponding to successive powers of . Minimising each yields

(43a) | ||||

(43b) |

Thus at in Region II, the minimising trajectory follows exactly the streamlines. The two solutions can be matched in their common region of validity, given by (and corresponding to ), to obtain

(44a) | ||||

(44b) |

At this order, the only non-zero contribution to the integral in (25) comes from the behaviour in Region I. We use (44a) to obtain that and thus

(45) |

since . Solving finally gives the approximation

(46) |

Here, is the principal branch of the Lambert W function [7]. The above results were previously derived in [36] and included here for completeness.

We obtain an approximation for in a similar way. The periodic trajectory associated with the variational principle (32) are divided into the same two regions as above. The regular expansions are this time more naturally expressed in powers of so that in region I where , we take

(47) |

where . In region II, and so we take

(48) |

where and once more extend the behaviour over the whole using symmetry.

The periodic trajectory is now obtained by substituting (47) and (48) inside the pointwise constraint in (32) from where we obtain equations for each power of . This leads to two sets of equations

(49a) | ||||

(49b) |

where and arise when parameterising the constraint (32) in polar coordinates. The minimum value of , denoted by , is obtained by maximising , and . This gives and leads to

(50a) | ||||

(50b) |

since , where is a constant to be determined. Matching between the solutions at in their common region of validity, given by (the same cell corner as above), yields an expression for . Using (32), we deduce that

(51) |

and . The order of the error is estimated by matching the solutions at (calculations not shown).

Figure 3 shows that expressions (46) and (51) are in excellent agreement with our numerical solutions; the same is true for expressions (44) and (50) describing the trajectories (not shown). We may use as to further approximate (46) as . This approximation highlights the leading-order difference between (46) and (51). However, this is only a rough approximation which cannot, for instance, capture the non-monotonic behaviour of that arises for small values (not shown). Note that both derivations of (46) and (51) tacitly assume that . This is easily shown to be the case once the behaviour of the trajectory over the whole (rather than a quarter) spatial period of the flow is taken into account.

#### 4.1.2 Large- asymptotics

A second asymptotic limit corresponds to . We extend the approach in [36] and take the minimising trajectory associated with the functional in (25) to be at leading order a straight line with higher order corrections given by a regular expansion in :

(52) |

where and . Here, is a constant and and are -periodic functions (with zero mean). Substituting (52) into (25) gives a sequence of integrals corresponding to successive powers of , obtained using a symbolic algebra package. These are in turn minimised up to with respect to , , , and (contributions from , , and cancel) yielding

(53) |

Introducing (53) into (25) we obtain

(54) |

after a few manipulations. This leads to the asymptotics of the speed

(55) |

with the first two terms previously derived in [36].

In a similar manner, the minimising trajectory associated with the variational principle (32) for (G) is at leading order a straight line. Using the alternative variational characterisation (34), we write the trajectory in terms of and take a regular expansion in powers of :

(56a) | ||||

(56b) |

where . The ’s are -periodic functions satisfying while for all . We substitute these inside the pointwise constraint in (34) from where we obtain equations for each power of . This leads to expressions for which are in turn used to minimise . Up to and after a few manipulations carried out with a symbolic algebra package we obtain that

(57a) | ||||

(57b) |

where and are -periodic and therefore do not contribute to the value of . Note that the difference between the two trajectories obtained in (53) and (57) only appears at . We finally use (34) to deduce that

(58) |

Comparing expressions (55) and (58) confirms that the difference between the front speeds for the (FK) and (G) equation is very small: equation (G) only slightly underpredicts the front speed. This is confirmed in Figure 3 which focuses on verifying (55) and (58). It is clear that the two approximations (55) and (58) are in excellent agreement with the numerical results; however, they are too close apart to distinguish.

### 4.2 Perturbed cellular flow

We now investigate the effect of perturbing the basic cellular flow by taking for in the streamfunction (5), keeping . The perturbation breaks a symmetry of the streamfunction. Characteristic examples of trajectories associated with (FK) and (G) are shown in Figure 4 (top row) for two values of corresponding to distinctly different flow topologies. The trajectories remain symmetric for the transformation . Qualitatively, they are similar to those obtained for , following closely the straight line when is large and the separatrix when is small. Despite the more complex flow structure, the difference between the (FK) and (G) trajectories remains small.

Figure 5 (top) shows the behaviour of as a function of . For , the value of does not greatly differ from the corresponding value obtained for . A significant difference is obtained for . For large , increases quadratically with . This can be shown by generalising the asymptotic result (55) to find, after a lengthy computation, that for