Dependence of advection-diffusion-reaction on flow coherent structures

A study on an advection-diffusion-reaction system is presented. Variability of the reaction process in such a system triggered by a highly localized source is quantified. It is found, for geophysically motivated parameter regimes, that the difference in bulk concentration subject to realizations of different source locations is highly correlated with the local flow topology of the source. Such flow topologies can be highlighted by Lagrangian coherent structures. Reaction is relatively enhanced in regions of strong stretching, and relatively suppressed in regions where vortices are present. In any case, the presence of a divergence-free background flow helps speed up the reaction process, especially when the flow is time-dependent. Probability density of various quantities characterizing the reaction processes is also obtained. This reveals the inherent complexity of the reaction-diffusion process subject to nonlinear background stirring.


I. INTRODUCTION
The studies on passive scalar mixing in advection-diffusion systems have a long history and naturally bifurcate into two prongs.1][12][13] Indeed, the marriage between the two prongs of studies on passive scalar mixing is on the finite-time behavior, where the decay of scalar variance is strongly influenced by local Lagrangian stretching.Recently, it is found that finite-time scalar variance is highly dependent on flow topology, with anomalous scalings of decay rates in different zones demarcated by LCS. 14 When the problem of interest is only passive scalar diffusion, and with finite diffusivity, as the time scale surpasses the Lagrangian decorrelation time, molecular diffusion will inevitably dominate the scalar behavior, and homogenization theory applies. 2,15 n shorter time, scalar behavior exhibits strong intermittency, and such behaviors are shown to be highly dependent on LCS. 14 This indicates that for real physical processes (with finite diffusion) that are shorter than the Lagrangian decorrelation time, LCS could be a valuable tool, not only describing the general mixing and transport patterns, but also as an indicator for distinct dynamical behaviors.One example of the role that LCS plays in identifying the change of dynamical behaviors in a biogeochemical problem is the primary production of phytoplankton near an island wake. 16In this study, two distinct cases have been generated with different biological inputs.It is seen that LCS, as characterized by finite-size Lyapunov exponents (FSLE), 17 plays an important role in incubating a phytoplankton bloom.a) Wenbo.Tang@asu.eduAlthough the mechanism for the two distinct scenarios of phytoplankton bloom was explained in the aforementioned study, the range of behaviors has not been explored.The fine boundaries of the emergence of one behavior from the other can only be understood with a broader range of realizations, and such behaviors are likely tied to the LCS.Motivated by this open question, we study the variability on scalar reaction in nonlinear aperiodic flows.In fact, the role that LCS plays may differ by reaction types.As such, we start from the simplest type of reaction -the diffusionreaction model governed by the Fisher-Kolmogorov-Petrovsky-Piskunov (FKPP) equation, coupled with background stirring.The reaction process is a very simple Verhulst model, and the long-term behavior is the saturation of the reactive scalar at its carrying capacity.However, at intermediate times, there is still considerable variability among the realizations.This variability is not yet important as far as the long-term behavior of the scalar concentration in the FKPP system is concerned.But the finite-time variability may play an important role in the overall reactive system when the FKPP model is coupled with other reaction processes.Henceforth, our goal in this study is to (1) generate a class of realizations characterizing the range of behaviors in the reaction process for the FKPP model at finite-time; (2) study the dependence of such processes with LCS; (3) study the dependence of such variabilities on the reaction and diffusion time scales; (4) generalize a Lagrangian approach for scalar variability for studies on other reaction processes.
Chemical reaction in advection-diffusion-reaction (ADR) systems has been studied from various angles previously.In a decaying reaction process, the decay rates of reactive scalars have been tied to flow stretching statistics, hence it establishes a connection of reaction with the bulk flow behavior. 18he mean and variance of scalar concentration in FKPP system subject to spatially varying carrying capacity have been studied in a random shear flow and in two-dimensional turbulence. 19,20 n terms of low order nonlinear dynamics, the propagation of reaction fronts in FKPP systems has been studied via burning invariant manifolds. 21,22 he propagation of reaction fronts in turbulence has also been studied via direct numerical simulations (DNS). 23In terms of chemical reaction in a simple vortex, variability of reaction was parameterized by the geometrical properties of the vortex. 24This study intends to complement the above aspects by linking the variability in reaction process with more complex flow geometry.More generally, advances in the dynamical systems approach to chemical and biological activities in fluid flows have been summarized in recent texts. 25,26 8][29] In advection-dispersion and advectionreaction systems, the probability density function for one-point statistics of scalar concentration field can be mapped to distribution of the random velocity field (as compared to LCS in this study), so the tail distributions can be quantified. 30,31 e rest of the paper is outlined as follows.In Sec.II we provide the background for the governing KFPP equation and the analytical tool.In Sec.III we present the results and discuss the various aspects of the dependence on flow topology and governing parameters with explanations.In Sec.IV we further discuss the implications of our results and draw conclusions.

A. Lagrangian coherent structures and the kinematic flow model
In the absence of diffusion, passive scalars conserve their density as they transport across the fluid domain.At the same time, they experience stretching and folding induced by the background flow.Characteristics of the emulsion process has been qualitatively described several decades ago. 8he theory on chaotic mixing depends on ideas from nonlinear dynamical systems, where the pattern formation is explained via fixed points and invariant manifolds in steady and time-periodic flows.In recent years, the definition and extraction of finite-time invariant manifolds have become an important tool in the study of chaotic mixing in real flows, because steady and periodic flows rarely exist in the nature.Roughly speaking, the finite-time invariant manifolds serve as locally strongest repelling or attracting material lines/surfaces, and nearby fluid trajectories diverge from or converge to them over the finite-time duration of their advection process.Therefore, these invariant manifolds serve as boundaries for different regions of coherent motion, hence boundaries for LCS.However, since the definition is over finite-time, these invariant manifolds are only valid over the time period where the trajectories are observed.As time progresses, these invariant manifolds may weaken and disappear in their role as organizing centers, and new material lines/surfaces will emerge and strengthen as the organizing centers for chaotic mixing at different times.
There are many tools developed to extract the finite-time invariant manifolds based on the Lagrangian trajectories. 10,11,13,17,32,33 A siple diagnostic tool to identify LCS boundaries is the finite-time Lyapunov exponents (FTLE). 346][37][38] Zones of FTLE have also been used previously to relate Eulerian and Lagrangian chaos in Kelvin-Helmholtz instability. 39e briefly outline the computation for FTLE.Denoting the current position of a trajectory x(t; x 0 , t 0 ) started from location x 0 at time t 0 , the deformation gradient tensor , is then defined as the scalar field that associates with each initial position x 0 the maximal rate of stretching: with λ max (J) denoting the maximum singular value from the singular value decomposition of J.
Because the FTLE field measures the largest separation among nearby trajectories, in forward (resp.backward) time, highlighters of FTLE indicate largest separation (resp.convergence) over the finite time interval of trajectory integration.
In recent studies, it has been emphasized that FTLE ridges can give false positives and negatives since they highlight both stretching and shear in the flow, and a closely related quantity, the "strainlines," is a more appropriate feature that highlights exponential separation away from the material lines. 33,40 t is also worth noting that FTLE alone does not outline the transport barrier that encloses elliptic structures (eddies).If these flow features are of interest, other tools should be used. 13,32,33 Wuse the kinematic model of double gyre flow to stir our fluid domain -its Lagrangian dynamics and properties of the transport barriers have been well studied 12,41 and it intimately relates to the chain vortex flow, the paradigm for both passive scalar mixing 42,43 and ADR systems. 21Here we use a previously considered velocity field, 44 defined as where g(t) = 1 sin (4π t) + 2 sin (2t) controls the characteristics of the background flow.When 1 = 2 = 0, the velocity field is the steady cell flow considered in classical studies.The velocity is periodic if one of the parameters is 0. Periodic velocities have been considered in previous studies on burning invariant manifolds. 21,22 hen both parameters are nonzero the velocity is quasiperiodic.Particle motion in both cases is typically chaotic, but for the first case, integrable regions such as KAM tori typically exist in the chaotic sea, and for the second case, these structures are typically destroyed.In our study, for most of the cases with quasiperiodic velocity, 1 = 0.3 and 2 = 0.1.For the one periodic velocity case of contrast, we set 2 = 0. Note that for a simple (as compared to turbulence) flow geometry such as the double gyre system, the strainlines match very well with the FTLE ridges. 41For this reason, we take direct comparison between the reaction behaviors in our system and the FTLE field.Additionally, we point out that it is not only the reaction behavior near the LCS boundaries, but also the behaviors away from them that are of our interest.

B. The governing equation and methodology
Consider a chaotic flow with length scale L and velocity scale U. A chemical species c is embedded in this background flow, whose density scale is C 0 , diffusivity is κ, and has an autocatalytic reaction with an intrinsic reaction rate k.The governing equation for the ADR system is the classical FKPP equation, 45,46 given in non-dimensional form as where u is the background velocity field (the double gyre flow as specified earlier), Pe = UL/κ is the Peclét number that characterizes diffusion time scale over advection and η = kC 0 L/U characterizes the advection time scale over reaction.Motivated by geophysical applications, consider a velocity scale of 0.1 ∼ 1 m/s, length scale of about 100 km, and an effective horizontal diffusivity at the ocean surface as 10 ∼ 10 3 m 2 /s, the estimated Pe = 10 ∼ 10 4 .The advective time scale is about several hours to a day.In terms of reaction, consider biogeochemical processes on the order of days, a reasonable choice for η is around η = 1.Accordingly, we chose our baseline case with Pe = 1000, η = 1.Variations based on parameter spaces and flow types are also considered.
The long term behavior of this ADR system is trivial -as long as c starts nontrivially, over long time, it will saturate to the carrying capacity (1 in nondimensional units).However, the time history of c approaching this limit may very well be non-trivial, and may lead to variabilities of behaviors when coupled with other processes.The statistics of the time history leads to a range of behaviors on the bulk reaction rates, and we would like to relate it to the LCS, if possible.Strategically, we study the reaction process of a highly compact scalar randomly appearing in our flow domain, but the background stirring is deterministic.This compact scalar is assumed to have equal probability to appear anywhere in the flow domain.Since we are not just interested in the probability density functions and the associated moments, but their spatial structures, we had to sample the fluid domain uniformly to fully understand this variability.
Note that in an ADR system, the dynamics of fronts are governed by the burning invariant manifolds, generally different from the fluid invariant manifolds.However, in our particular case, where reaction time scale is comparable to advection time, and diffusion time scale is much longer, we are safely in a regime where the nondimensional frontal speed v 0 = 2 √ η/Pe 1 (except for Pe around 10), where the burning invariant manifolds converge to the fluid invariant manifolds. 21enceforth, at least for the high Pe cases, it is meaningful to just consider the flow LCS.
There are several levels of information that we extract.At the most fundamental level, scalars are undergoing the ADR process, and the stretching and folding history of the scalar fields will be discussed.To visualize the spatial structures of the variability in reaction, at the secondary level, we construct a scalar field, where we assign to each location x 0 the value of c at time t when starting at t = 0 with a localized excited patch at x 0 .Here, each realization of the ADR process is represented as the average concentration c , where is the domain average of scalar q.This synthesized scalar field is a function of initial source locations and time of simulation.This is the primary result we relate to LCS.Statistics of the time history of c subject to initial conditions, flow types and governing parameters are also studied to explain the overall process.We use DNS for our problem. 47,48 he fluid domain is considered doubly periodic, but some special care must be taken into consideration.To avoid the velocity jump at the boundaries when only two gyres are considered, we extend the domain to simulate four vortex cells in 2 × 2 arrangements.Since the flow is symmetric, we only have to consider initial release in two cells.The chemical concentration is initially released as a Gaussian patch of peak value 1, and a full width at half maximum of 0.053.This corresponds to a dimensional length scale of 5 km.They vary just by spatial locations as specified later.For each realization the scalar density is numerically integrated towards an asymptotic state.For η = 1 this takes about 15 time units for the baseline case.We then study the full range of behaviors by varying Pe and η, and consider a few extreme cases in other flow types.256 grid points are used in both directions of this two-dimensional flow, and we uniformly sample the domain with ∼2000 realizations for each case of fixed Pe and η.

A. Variability on reaction
We summarize the different behaviors in our baseline case with Pe = 1000, η = 1 in Fig. 1. Figure 1(a) shows the time history of the average concentration c over time, for two bounding cases.All realizations have been exhausted towards the asymptotic state where c is almost 1 everywhere so the domain average c is 1.However, it is clearly seen that there is a band of behaviors subject to the initial peak locations, and it is not clear if the probability density of this band is going to take any simple structure.Since the largest variation appears around the middle of the simulation, we show the probability density of c when T = 7.5 in Fig. 1(b) as the dashed curve, and the probability density of T that reaches c = 0.5 in Fig. 1(c) as the solid curve.The vertical dashed line and horizontal solid line in Fig. 1(a) show the slices taken to generate the pdf for Figs.1(b) and 1(c).It is not surprising that the two pdfs are somewhat reciprocal.There is a higher than normal probability that c reaches 0.5 at T = 7.5.Also somewhat shocking is the difference in bulk concentration for the extreme cases at this time, which is about 20% of the average behavior.The implication of this could be profound -if there is a chemical process that takes a comparable time scale of 7.5 time units and coupled with the FKPP dynamics, the difference in extreme cases can be substantial.
To better understand this variability, we plot the LCS as revealed by forward-time FTLE and the spatial structure of c in Fig. 2. Figure 2(a) shows c at T = 7.5 when the variability is the strongest.Each grid cell indicates the peak location for each realization, and the scalar value at this grid is c(T = 7.5) for this realization.Figure 2(a) has a grid resolution of 64 × 32, hence a total of 2048 sources are uniformly sampled in the ensemble, with grid locations g i /32, g j /32 where g i = 1, 2, . . ., 64; g j = 1, 2, . . ., 32 are the integer indices in the x, y directions, respectively.The FTLE field is obtained by integrating deterministic Lagrangian trajectories inside the domain over one time unit so it characterizes the early time Lagrangian dynamics (for Fig. 2(b), integration time is T = 0 → 1).The auxiliary trajectory approach 41 has been employed to improve the accuracy of the FTLE field.We use the forward-time FTLE because they highlight locations of significant initial stretch of the scalar in the fluid domain, so that a larger area covered by the scalar will see growth quicker, contributing to the faster growth of the average concentration.As comparison, backward-time FTLE extracted at this time highlights dynamics between T = 0 → −1 and does not directly relate to the initial scalar stretch behavior.However, as shown later, for scalar concentration field c we will have alignment of scalars to backward-time FTLE, since scalar stretching has already taken place during the time when the backward-time FTLE is computed.As seen, in Fig. 2, the strainlines, conveniently revealed by the FTLE, maps precisely to the variability plot in bulk concentration -at regions where repelling motions are the strongest, the bulk concentration is larger, due to more effective spreading, and hence homogenization of the scalar; at regions where repelling motion is weak, the scalar is less spread, resulting in a poorer reaction process.This indicates that the FTLE field would be a useful tool to represent the range of behaviors in reaction in the ADR system considered at hand.Since the FTLE computation is a lot cheaper, it significantly helps the parameterization of the variability in reaction.
Before discussing any possible modeling of the variability, we analyze the time history of the reaction behaviors in more detail.First, we analyze the range of behaviors again, for the chaotic stirring case and a few other cases.These cases are revealed in Fig. 3.In Fig. 3(a) we plot again the time history of the realizations in black.All realizations have been plotted, resulting in the black band between T = 0 and 15.We consider two extreme cases as comparison.First, we consider the case where c is uniform and equal to the average concentration c of a compact patch release.Effectively, this represents a homogeneous concentration case whose dynamics is simply the ordinary differential equation The analytic solution of this differential equation can be easily obtained.For our case, the role that stirring plays is to help homogenize the scalar field, hence an initially homogenized scalar field can be considered as well-stirred, and will see the fastest growth.The time history of the homogeneous release case is given in the red curve (which is indistinguishable from the black band but can be seen in the inset plot of Fig. 3(a)).It remains as the upper bound (fastest reaction) for all possible realizations/incompressible flow types with the given η.In contrast, we consider the case of a compact patch release in no flow.Since no stirring is present, the scalar spread is the slowest, and the resulting time history, plotted in the blue curve, does show this slow reaction.In fact, it takes about twice the time to reach the asymptotic state.
In addition, we consider two more realizations in steady cell flows.Their behaviors should be between the no-flow case and chaotic stirring.One realization is with the peak concentration released near a hyperbolic fixed point, at the corner of a cell.The scalar peak is helped to redistribute in the cell boundaries quickly, and gradually invades the cell center.The other realization is in the center of the cell, where the fluid undergoes simple circular motion, and the scalar is not helped to redistribute in early time.The time history of these two realizations are shown as the green dashed curves in Fig. 3(a).It is seen that for large c , the realization with faster reaction behaves similar to the slowest reaction case in the black band, whereas the realization with slower reaction lags behind (but still much faster than the no-flow case).
We closely examine the early time behaviors of these curves in the inset plot of Fig. 3(a).The first two time units are shown.It can be seen that the red and blue curves bound the range of reaction behaviors.In terms of the black band, all realizations start close to the blue curve, but the realizations close to the strainlines quickly adapt to approach the red curve.In comparison, on the other end of the spectrum, the realizations furthest away from the strainlines follow the blue curve for significantly longer time, and deviate from it around T = 1.8.With the steady cell case (the two green dashed lines), the early time behavior is quite similar to the chaotic flow case -the realization stretched the most in the hyperbolic region quickly adjusts to follow the red curve, while the realization in the cell center behaves a lot like the no flow case for a long time.However, since the flow is steady, scalars in the realization with the fastest growth will not penetrate into the cell center until late, so its late time behavior lags behind the chaotic stirring ones.
We explain the dynamics of c by first averaging Eq. ( 2) and applying the periodic boundary conditions.Separating the scalar concentration into mean and perturbations (c = c + c p ) for c 2 and considering c p = 0, we obtain At the initial time c(0) = c h (0), where c h (t) denotes the solution to Eq. ( 3) from a homogeneous release.Any inhomogeneity in c (nonzero c p ) will result in a slower speed in c (due to concavity in Eq. ( 4) 26 ).Since the growth is slower, c for any inhomogeneous release lags behind the homogeneous release case, hence the homogeneous case serves as the upper bound (fastest reaction) for the range of behaviors for any flow given a fixed intrinsic reaction rate η.As seen, the larger the inhomogeneity, the slower the reaction for c .Because the no-flow case has no stirring to help on homogenizing the scalar field, it serves as the lower bound (slowest reaction) for the range of behaviors.
Since the classical tool to analyze the simple dynamical system Eq.( 3) is the phase portrait of c against dc/dt, we use a similar idea to study the variability in reaction rates, by analyzing d c /dt.However, it's behavior is very close to the right hand side of Eq. ( 3), we thus highlight the difference by plotting d( c − c h )/dt against c .These curves are plotted in Fig. 3(b).As seen, all cases lag behind the homogeneous case (red line) because the difference in bulk reaction speed has a dip around c = 0.25.The smaller the dip, the faster the bulk reaction.However, as c passes about 0.7 these reaction behaviors quickly converge to a single behavior, indicating an almost uniform late-time dynamics.Henceforth, it is crucial to understand the early-time dynamics, when the scalars are still undergoing stretching, and diffusion is not yet smoothing out the small-scale heterogeneities in c, as they dictate the variability in c .We also show in the inset plot of Fig. 3(b) the small c behavior in log scale.As seen, the faster reacting realizations quickly adjusts to 0 (representing the homogeneous case), and the slower reacting realizations follow the blue curve longer.The small vertical line marks the start of the visible departure of the black band from the blue curve, at c ∼ 0.0035.This location has also been marked in the inset plot of Fig. 3(a) as the horizontal line.The slowest realization reach this value at T ∼ 1.8.
As the c field appears to correlate with the FTLE very well in Fig. 2, we study the time history of the correlation between the two maps, given as where std(q) denote the standard deviation of function q and c denotes the spatial average of the synthesized field c .As seen, this correlation remain high and asymptote to a constant value as T becomes large.(Note that even though at large times the c field approaches a constant, there is still very small scale variation.The correlation function highlights the relation between the structures, hence it does not decay to 0.) In the inset plot of Fig. 3(c), we show the probability density of c and FTLE, where the FTLE is associated with the black curve and the black coordinates, and c is associated with the blue curve and the blue coordinates.These two pdfs, at least for the parts where the probabilities are large, do seem quite similar.As a reference, we use the black horizontal line to mark the FTLE values between 2 and 4, and the blue horizontal line to mark the c values between 0.52 and 0.55.We then locate those regions in Fig. 2. The regions correspond to the band of structures near the FTLE strainlines and those in the high reaction regions.Their pdf, as part of the pdf for the entire domain, also appear to be similar, indicating that it is promising to match the reaction behaviors with the FTLE, structure by structure.

B. Reaction and Lagrangian stretching
To understand the cause of the range of c behaviors we study four representing realizations in the chaotic flow case.These cases are chosen to be initiated near a hyperbolic core (denoted as a black dot), the interior (white dot), a piece of repelling structure not near an attracting one (red dot), a piece of the attracting structure not near a repelling one (magenta dot).The colored dots are shown in Fig. 2. In Fig. 4, the scalar fields c are shown at T = 1.This is the time where c is around 0.002 for all cases, but the individual scalar structures vary significantly.In Fig. 4, the panels (a)-(d) correspond to the releases at black, white, red, and magenta dots, respectively.For reference the forward-and backward-time FTLE fields are indicated as the black and white iso-contours (integration started at T = 1 for one time unit).The titles of each panel indicate the average concentration c at this time.Note that due to the north-south symmetry of flow we only show two cells.The scalar field is not necessarily symmetric, but their behaviors are well described in the two cells shown.At T = 1, c that corresponding to the uniform density case is 0.008606 and that corresponding to the no-flow case is 0.001768.
Later in the simulations, at T = 7.5, we plot in Fig. 5 the scalar density field c for the four cases again.At this time, the three cases with the initial peak not in the interior (Figs.5( are relatively homogeneous compared to the case started in the interior (Fig. 5(b)), resulting in the huge difference in c for the four cases (note the scales in colorbar).The forward-and backwardtime FTLE corresponding to this time are also plotted to show the alignment of scalar gradient with the attracting LCS (integration started at T = 7.5 for one time unit).As a reference, at T = 7.5, c that corresponding to the uniform density case is 0.5902 and that corresponding to the no-flow case is 0.1960.Note that at the three times where LCS is shown (T = 0, 1, 7.5), the structures happen to look similar.This is due to our choice of g(t) has a predominate frequency 4π associated with the larger 1 , leaving the coherent structures to slosh back and forth with a period of about 0.5 time units.The LCS vary quite a bit within each period.
To better understand why the variability is intimately related to LCS, we observe that at early time the dominant scalar behavior is the decay from the peak at the source.After some time, reaction overcomes the decay and c increases everywhere in the entire domain.Based on Eq. ( 4  the evolution of c is critically determined by the variance c 2 .In classical approaches to scalar decay in turbulent flows without reaction, scalar variance is only based on the scalar perturbation field because the scalar is conserved over the domain of interest.Here since the mean field is also changing with time due to reaction, we analyze both of the two components c 2 and c 2 p .Multiplying Eq. ( 2) by c and c p and taking domain average we obtain We show the evolution of c 2 and c 2 p along with the contributions from individual terms on the right-hand side of Eqs. ( 6) and ( 7) in Fig. 6.
Figure 6(a) shows the relevant terms in Eq. ( 6) for the four cases.Shown in black are c 2 .Since they are plotted against c , they overlap and appear as a straight line in log scale.The blue curves are the first two terms on the right-hand side of Eq. ( 6), capturing change in c 2 due to reaction in the mean.As they are simple functions of c they again overlap.The green curves are 2η c c 2 p , where the behavior of the four cases are now distinct against c .Among the green curves, the thick solid curve corresponds to the realization released at the hyperbolic core, the thick dashed curve at the interior, the dashed-dotted curve at the repeller and the dotted curve at the attractor, respectively.Note that the difference among these curves do not alter the shape of c 2 against c , they only affect the speed at which a realization moves along the phase space described in Fig. 6(a).The larger the perturbation variance, the slower the scalar homogenizes, hence the slower the realization moves along the curves (from left to right).The red vertical line marks c = 0.002, corresponding to T ∼ 1 for all cases.local minimum, due to reaction surpassing dissipation, leading to the growth of c 2 p .Then, at large time, as the scalar c saturates, c 2 p decays again.Figure 6(c) shows the evolution of c 2 , which is the total contribution from Eqs. ( 6) and (7).The line styles are the same as Fig. 6(a) for the individual cases.The black curves are c 2 , the blue curves are contribution from reaction 2η c 2 − 2η c 3 , and the green curves are scalar dissipation 2 |∇c p | 2 /Pe.It is clear that for c 2 , which critically determines the growth of c in Eq. ( 4), the early contribution is due to scalar dissipation.As the mean dynamics surpasses scalar dissipation, and as c p decays in large time, the late contribution is mainly due to the mean contributions.Henceforth, the evolution of c at large c can roughly be considered uniform.In Fig. 6(c), the red solid line corresponds to c = 0.002, indicating where T ∼ 1.The entire reaction process can roughly be separated into four phases, as indicated by the three red dashed lines.In the first phase, when c < 0.001, c 2 undergoes decay, mainly due to contribution from scalar dissipation.In the second phase, scalar patch has reached its Batchelor scale width, where the peak generally decays to be much below 1 (except for the realization released in the interior), so c 2 grows due to the lengthening of the patch as well as (approximately) linear reaction.This phase roughly stops at c = 0.01.In the third phase, lengthening has created close packing of thin filaments, which now starts to homogenize (with the exception of some transport barriers which are hard to penetrate).The homogenization process in this insufficiently mixed stage contributes to the individual realizations approaching the point where c 2 = 0.01 at c = 0.1 (from different pathways and at different times), following an approximate power-law growth.Finally, in the fourth phase when c > 0.1, stirring does not help differentiate reaction rates considerably anymore.Note the different scaling behaviors in the black curves at different phases separated by the three red dashed lines.Also worth nothing is that for the realization released in the interior (thick black dashed line), there is no strong distinction between phase two and three, mainly since in the interior the strain is not large, so the scalar patch decays to a length scale comparable to the domain size after the first phase, and already starts to behave like phase three from then on.
Figure 7 illustrates the scalar fields for two realizations at the transition of the phases.The top row in Fig. 7 is for the realization released at the hyperbolic core and the bottom row in Fig. 7 is for the realization released in the interior.The left column shows the transition from the first phase to the second, where the scalar patches decay toward the Batchelor scale width.The center column shows the transition from the second phase to the third, where nearby filaments get very close and start homogenizing.In another interpretation the scalar filaments have filled the zones where they were initially released and start homogenizing in that zone.(Note that the behavior in the bottom row behaves differently as the scalar patch stays relatively confined, but it does fill the interior of low dispersion region and so the behavior is homogenization within this region.)The right column shows the transition from the third phase to the fourth, where homogenization within one zone is about to complete and the high concentration region finally invade the unoccupied regions (either from chaotic zone to interior or from interior to chaotic zone).Now we consider properties of the flow that control the early-and late-time dynamics of the scalars.At early time in the first phase, since the scalar is still highly concentrated, we consider the reduced problem where the scalar is undergoing reaction-diffusion in a local linearized flow.When the local flow is linear, assuming that the average dimensionless stretching rate is λ, after time T the length of the scalar patch starting from size L 0 is roughly L 0 e λT .Without diffusion the scalar width should be L 0 e −λT .Considering diffusion, when strain and diffusion balance the Batchelor scale is achieved at L B = 1/ √ λPe.For all four realizations, at the end of phase one, the scalar length roughly corresponds to L 0 e λT with λ given by the local FTLE, where L 0 = 0.053 is the original patch size.The scalar width also reached the Batchelor scale as given by λ.
Consider a flow of the form û = −λ x, v = λ ŷ.Here x and ŷ are local coordinates along which the principle axes of the strain align.At the very beginning of time, scalar peak is of order 1 and scalar gradient is concentrated near the peak so the major contribution to the integral |∇c p | 2 is limited to a small area where c ∼ 1.We hence approximate the scalar dynamics as pure advection-diffusion governed by with natural boundary conditions at x = ∞.The inclusion of the linear strain flow helps us understand the variability in the decay at this time.We must include finite diffusion here since without it, the energy dissipation rate inside a linear strain flow remains constant, contrary to what is observed in Fig. 6.At early time, we assume that the kernel (Green's function) takes the form which is the composition of the effect of exponential separation of trajectories and the decay in scalar concentration.This expression is valid up to first order in t and the total density is conserved.
Convoluting the kernel with the initial Gaussian distribution we have where σ 2 is the variance of the initial profile.Note that we are able to integrate to infinity because we assume that the scalar is well within the size of the domain, so effectively the contribution outside of the domain is negligible.The scalar dissipation rate based on this expression is (the averaging is based on the periodic domain size as we have one patch in every 2 × 2 square) The early time behavior (within the first phase) of Eq. ( 11) is the decrease of dissipation rate, with the case of larger λ decreasing slower.Note closely at the left of Figs.6(b) and 6(c), for the dissipation rate, the thick solid curve (hyperbolic core) is indeed slightly larger than the thick dashed curve (interior), consistent with the above analyses.After some time, in the second phase, an initially circular patch of scalar will be lengthened in the ŷ direction, along which the gradient will decrease.Because of the increase in patch length and maintenance of scalar width, scalar peak falls well below 1 and we can omit the ŷ derivatives and include a linear approximation to reaction (since (1 − c) ∼ 1; for realizations in the interior, peak concentration has not decayed enough, so the effective η for realizations starting from the interior will likely be smaller in the equation below, we hence use η eff ) where η e f f = η(1 − c) and c is roughly the average concentration inside the patch.Note that Fourier transform of the above equation is still a partial differential equation.We proceed by separation of variables and analyze the eigenvalue spectrum.Placing the reaction term η eff c together with c t , the eigenvalue function X ( x) satisfies with natural boundary conditions at x = ∞; and the time dependence T(t) for mode δ satisfies where δ ≥ 0. Equation ( 13) has closed form solution when δ is integer multiple of λ.We focus on the slowest decaying eigenfunction with δ = 0, where X = c 1 exp(−λPe x2 /2). 26Note that the peak c 1 decays as e −λt due to the stretching along ŷ and conservation of scalar in the absence of reaction.Therefore, the scalar dissipation rate as a whole behaves as O(e −2λt e λt e η ef f t ) = O(e (η ef f −λ)t ).Since the variation in λ over the domain (between 0 and 5 in this flow) is much larger than that of η eff (between 0 and 1), this relation roughly explains the behavior in phase two -the larger λ is, the smaller its associated dissipation rate will be.Note that within the second phase, the exponential rate of stretching is not purely given by the local FTLE -as the scalar lengthens it tries to fill up the coherent structure where it started, which generally has a range of FTLE values.Henceforth, the filament width varies within a coherent structure and there is a limit in length that a filament can grow to.This makes direct extrapolation of Eq. ( 12) to represent the entire second phase difficult as it is based on a constant local stretching rate.
After yet some other time, in the third phase, there is sufficient homogenization within each zone, scalar density within each zone can roughly be treated as constant, hence the dependence of c 2 with respect to c is simply quadratic, if the angled bracket denoting spatial average is within each zone.Since scalar value is different in the zones, average over the entire domain will modify the functional relationship to be different from quadratic, but as the density fields approach closer to carrying capacity, this dependence will become quadratic.Henceforth, we observe that the scaling law dependence of c 2 on c to be smoothly transitioning from where the realizations end in phase two towards a quadratic dependence at the start of phase four.
The observation that the evolution of c 2 becomes uniform after c reaches 0.1 in phase four is not too surprising.Consider the extreme case where the scalar is most inhomogeneous.The scalar would have value c = 1 inside the patch and c = 0 outside the patch.In order to achieve an average of c = 0.1, the patch size has to be roughly 0.45 in two cells.This is already comparable to the size of the cell.Motivated by this argument, we explain the large c uniformity via scalar length scales.
First, we obtain the integral length scale of the scalar field, as a function of time.This length scale is given by The results are plotted in Fig. 8(a).The black curves (line styles the same as those in Fig. 6) show the growth in scalar integral scale as c increases.At c = 0.1, the four cases have L c between 0.1 and 0.3, indicating that the scalar size has increased considerably and local flow structures alone no longer dictate scalar evolution.The scalar samples different flow structures as sources to perturbations, instead of major conveyers, as when the scalar size is small.Second, we contrast the fractal dimension between the ADR system and pure chaotic stirring.It was shown that, for chaotic stirring, in the limit where Pe → ∞, the Hausdorff dimension for scalar gradient is less than the spatial dimension of the problem, i.e., the measure is fractal. 49,50 he more the scalar gradient is fractal, the less homogenized the scalar field (in the sense that the structures are dominated by filaments, instead of patches).With finite diffusion, the Hausdorff dimension is the same as the physical dimension of the problem.Henceforth, we use information dimension to characterize the non-uniformity of the scalar field.Specifically, we compute the information dimension, defined as where i loops through the entire domain and P i denotes a normalized scalar density at box i such that N i=1 P i = 1.We expect that the information dimension of our scalar should approach the spatial dimension at large c when the scalar field becomes more uniform.A similar scale-dependent measure was developed to quantify variability of scalars at the sea surface. 51Here, since our smallest scale can be taken as the grid size, we compute the information dimension based on the finest resolution of each realization.The dependence of D 1 on c is shown in Fig. 8(b).Clearly, at c = 0.1, all cases have reached a dimension of 2, equal to the capacity dimension (indicating that the scalar has already been homogenized considerably).

C. Parametric dependence on LCS
In this section we discuss the modeling of the variability observed in the baseline case using the FTLE field.This model is a mix of analytical results at early times and numerical observations at late times.Motivated by discussions in Subsection III B and the traditional approach to autocatalytic reaction, we model the range of average behavior with two degrees of freedom.(Note also that the uniform behavior in phase four starts at c = 0.1 for c 2 , where as d c /dt is not uniform until at a much larger c , cf.Fig. 3(b), making it desirable to include c 2 in the model.) We split f into three phases.In the first phase, the process is dominated by the decay of the scalar peak.Since the initial scalar patch has finite width, FTLE varies even within this small patch.To make the computation of c 2 analytically tractable, we use a weighted FTLE, defined as , as the average stretching rate λ for the local flow.
In the first phase, since the contribution of reaction is small, we assume that the scalar concentration changes its structure (but not the amplitude) the same way as pure diffusion.We use Eq.(10)  to estimate c 2 and subsequently evolve c based on Eq. ( 16).We then adjust the amplitude of the concentration to match the computed c as we advance, as a way to follow the contribution from reaction.Based on Eq. (10) in this phase, where c a (t) is the adjustable amplitude due to contribution from reaction which starts at 1, and σ 2 = 1/200π 2 is the variance of the initial patch.This phase is stopped at T 1 (λ) when c = 0.001.As explained in Subsection III B, direct extrapolation of exponential growth of filament based on local FTLE is not ideal for the modeling of the second phase.Since global chaotic mixing already starts in this phase (in the sense that the scalar patch has grown long enough to sample different FTLE within a coherent structure), we model the second and the third phase together, as smooth transition from a ( c , c 2 ) pair to the point (0.1, 0.01) in phase space, as a straight line in log scale (simple power law function of c ).The resulting c 2 at every c is used to advance c to the next value.Given that c 2 = γ at the end of the first phase, which is used in Eq. ( 16) to advance forward.Finally, in the last phase, we numerically solve for the problem since the functional dependence of c 2 on c now is the simple quadratic relation.Figure 9 shows the results of reaction rates based on the FTLE model.In Figs.9(a) and 9(b), the fastest and slowest reaction processes from the model are compared with the DNS results.For both panels, the black solid curve is for the release at hyperbolic core based on DNS.The red solid curve is comparison from model, Eq. ( 16).The black dashed curve is for the release at the interior based on DNS.The red dashed curve is from the model.Figure 9(a) shows evolution of c with time, where the red curves closely approximate the black curves from DNS. Figure 9(b) shows the functional dependence of c 2 on c .It is apparent that the approximation is better for the interior release case as compared to that for the hyperbolic release case.This is not too surprising as the process in the interior release case is the filling up of the interior, and since the patch scale is already large, homogenization (third phase) already starts at the end of the first phase.It is also clear that the first and last phases are well approximated with the current approach.A better approximation for the chaotic zone in the second phase will lead to better modeling of this variability.Figure 9(c) shows c at T = 7.5 for all realizations based on the model.It is not surprising that it will bear close resemblance with FTLE and the simulated c , as it is purely based on F T L E w .However, the important point is that by this modeling approach, we essentially reconstruct the probabilistic structure and range of behaviors for the reaction rates at any given time, not entirely by data fitting, but by a blend of analyses and observations.This also suggests a systematic approach to the studies on uncertainties in finite-time diffusion-reaction processes with the presence of strong coherent structures -based on behaviors of limiting cases, and careful analyses of the phases of diffusionreaction, one can obtain the generic features within each phase, hence develop models that grasp the essential dynamics within each phase.

D. Dependence of variability on governing parameters
We use F T L E w to compute correlation between LCS and c .Because the scalar patches start with finite width, and the unweighted FTLE strainlines (roughly approximated as ridges in this flow) are infinitesimally thin, the weighted average helps giving FTLE highlighters a sense of width.This effectively increases the correlation between LCS and c .We discuss the relation between c and LCS dependent on variations in Pe and η based on this correlation.In addition, we measure the range of variability as the maximum gap between the fastest and slowest realizations in a given pair of (Pe, η) at any given time.This gap is defined as Figure 10 shows the variation of the correlation as functions of time.Figure 10(a) shows such a correlation subject to changes in Pe while keeping η = 1.Since the range of behaviors of interest are between Pe = 10 ∼ 10 000, we generate four cases, with a uniform increase by a factor of 10 between each case.Specially, for Pe = 10 000, to resolve the fine structures in the concentration field, we decrease the grid size by a factor of 2. In this figure, the black curve shows the case with Pe = 10 000, the blue curve shows Pe = 1000, the green curve shows Pe = 100, and the red curve shows Pe = 10.Apparently, except for Pe = 10, the correlation between LCS and c remains quite high.Note that at large times when c almost saturates, this high correlation is merely an artifact that the variance in c is small.Since the highest variation in c is seen at half of the time when most realizations reach c = 0.5, we mark this time by scattered dots in respective colors.For Fig. 10(a), this time is T = 7.5.This clearly indicates that Pe = 10 (red curve) is an outlier.In the lower inset plot, we show the c field at this time for Pe = 10.Strong diffusion at this Pe renders many of the realizations in the filamentary structures to be indistinguishable from their neighbors, and only the realizations starting near the hyperbolic cores have higher reaction rates from the rest of the domain.Note that Pe = 10 is a case when the flow LCS do not correspond to the burning fronts dictated by burning invariant manifolds.Whether the structures seen in this plot have any relation with the burning invariant manifolds is a subject of future study.In the upper inset plot, we show the maximum gap in c as a function of Pe for η = 1 to measure the variability of reaction for this (Pe, η) pair.At low Pe, because large diffusivity homogenizes the scalar concentration faster, we see more resemblance of a homogenized release case, hence the variability (gap) is very small.As Pe increases, this variability also increases.The fact that this gap starts to approach a constant limit is probably due to the finite width of the initial patch -early-time scalar behavior is somewhat similar if the patch is released close enough to a LCS barrier.We expect to see more distinct initial scalar behaviors with smaller patch sizes.Figure 10(b) shows the variation of correlation subject to changes in η while keeping Pe = 1000.Due to the linear relation between saturation time and 1/η, hence the long time required to compute cases with small η, we only sample the η space around 1. The black curve shows the case with η = 4, the green curve shows η = 2, the blue curve shows η = 1, the red curve shows η = 1/2, and the magenta curve shows η = 1/4.We only run the simulations up to when the scalar fields saturate, hence the simulation is the longest for η = 1/4 and the shortest for η = 4. Since for all cases, the correlation at the end of the respective simulations start to approach their respective asymptotic values, we extrapolate the data to T = 60 (when the η = 1/4 case finishes) using the last value of correlation at the end of individual cases.These extrapolated data are shown as dashed lines.The scattered dots again mark when the variability in each case is the strongest.We show the case with the smallest asymptotic correlation (black curve with η = 4) at T = 1.875 (time of maximum gap) in the lower inset plot of Fig. 10(b).Compared to that in Fig. 10(a), this figure is a lot more similar to the LCS.We also plot the gap as a function of η at Pe = 1000.Interestingly, a maximum is attained at η = 2 for the cases considered.This could be explained as follows.At small η, reaction is slow.It takes longer for the mean concentration to reach c = 0.5 and hence diffusion better homogenizes the field, hence shortening the gap.At large η, reaction dominates over advection and diffusion, hence the behavior of the concentration field is more like a reaction front from a circular patch release.Such a case is not affected much by release locations (inherently by the flow topology), hence the gap is also smaller.
To summarize, in order to observe larger variability among realizations due to local flow topology, one wants Pe large and η to be order 1.

IV. DISCUSSIONS AND CONCLUSIONS
We studied the variability in the reaction processes associated with the Fisher-Kolmogorov-Petrovsky-Piskunov equation, subject to chaotic background stirring from a double gyre model.We find that the variability subject to randomness in the initial locations of chemical patches are strongly tied to the repelling Lagrangian coherent structures as highlighted by the FTLE field.This variability is mainly arising from the early time stretching of the scalar when it is still highly localized, hence it intimately relates to the LCS (since the LCS highlights early-time local stretching).A parameterization of this variability based on FTLE is formulated, and its dependence on the controlling parameters Pe and η are also discussed.These findings suggest that LCS can be used to better understand the uncertainties associated with chemical processes due to flow advectiondiffusion, and it can be used as an efficient tool in parameterizing such uncertainties.
We note that the FKPP system is elementary, and serves as the entry point for our understandings of the role that LCS plays in advection-diffusion-reaction systems.Other chemical processes may lead to more complicated behaviors and subtle roles that LCS plays, and should be carefully studied on a case by case basis for the particular reaction terms involved.Studies on more complex chemical interactions are currently underway.We also note that the double gyre system is elementary in that it provides simple flow geometry where strainlines correspond well with the FTLE ridges.In more complex flow geometries, such as those given by turbulent flows, this simple picture may be violated, and the role that LCS plays and the way to parameterize variability in chemical reactions with LCS may require more careful considerations.Such studies are also in our plan.Finally, we state that the current problem is studied in two dimensions, where it is still possible for us to exhaust the realizations to characterize the range of behaviors.How does our results extrapolate to three dimensions and with more complex reactions and flows is also being studied.

FIG. 1 .
FIG. 1. Range of behaviors of the reaction process.(a) Variation of the average concentration c in the periodic domain.The two black curves are the bounding cases of the fastest and slowest reaction behaviors among the realizations.(b) Probability density of c at T = 7.5.The range of behavior lies on the vertical dashed line in (a).(c) Probability density of the reaction time when c = 0.5.The range of behavior lies on the horizontal solid line in (a).

FIG. 2 .
FIG. 2. Comparison between the variability in c and LCS.Each point in (a) marks the center of a concentration blob initially released in the flow (a realization).The color value is the associated c of this realization at T = 7.5.(b) Repelling LCS marked by the forward-time FTLE field.The four colored dots mark the initial conditions of four representative realizations discussed later.Note the range of c at this time and the strong correlation between the two panels.

FIG. 3 .
FIG. 3. Range of behaviors of the reaction process.(a) Placement of the reaction processes in the chaotic flow with other flow types.The inset plot shows the early-time behavior of these cases.(b) Difference in reaction rates between the actual realizations and the nominal logistic equation case red line).The inset plot shows the early time behavior in log scale.The vertical line is a reference line at c = 0.0035.This line is also shown in the inset plot in (a).The blue line is the realization with no flow.(c) Correlation between the structures in Fig. 2 as a function of time.The inset plot shows the comparison between probability density of FTLE (black) and c (blue).The black horizontal bar marks the range of FTLE between 2 and 4. The blue horizontal bar marks the range of c between 0.52 and 0.55.
FIG. 4. Chemical concentration at T = 1.The four realizations are initially released at (a) black dot, (b) white dot, (c) red dot, and (d) magenta dot in Fig. 2. The black (white) curves are the iso-contours for forward (backward)-time FTLE indicating repelling (attracting) structures.The titles indicate the average chemical content c .
FIG. 5. Same as Fig. 4 but at T = 7.5.Note the major difference between (b) and the other cases (a, c, and d).

FIG. 6 .
FIG.6.Evolution of individual contributions from Eqs. (6) and(7) to the scalar variance c 2 , c 2 p , and c 2 , as functions of c .(a) Terms relevant to Eq. (6).The black curve is c 2 , which is a straight line in log scale as a function of c .The blue curve is 2η c 2 − 2η c 3 .The green curves are 2η c c 2 p for the four cases, plotted in thick solid (hyperbolic core), thick dashed (interior), dashed-dotted (repeller), and dotted (attractor) lines.The red vertical line references c = 0.002, where all cases are around T = 1.(b) Terms relevant to Eq. (7).The black curves are c 2 p .The blue curves are 2η c 2 p − 4η c c 2 p − 2η c 3 p and the green curves are 2 |∇c p | 2 /Pe.Line styles are the same as those in panel (a).The red vertical line again marks c = 0.002.(c) Terms relevant to evolution of c 2 .The black curves are c 2 , the blue curves are 2η c 2 − 2η c 3 and the green curves are 2 |∇c p | 2 /Pe.Line styles are the same as those in panel (a).The red solid vertical line marks c = 0.002, the three red dashed vertical lines mark c = 0.001, 0.01, and 0.1, respectively.All realizations reach c = 0.1 between T = 5 and 5.5.

Figure 6 (FIG. 7 .
FIG. 7. Phases of the entire process for two realizations.Top row: release at hyperbolic core.Bottom row: release at interior.Left column: end of the first phase dominant by scalar decay.Center column: end of the second phase dominant by filament stretching.Right column: end of the third phase dominant by homogenization among neighboring filaments.End of the fourth phase is uniform concentration and is not shown.

FIG. 8 .
FIG. 8. (a) Integral scale the scalar.line styles are the same as those in Fig. 6.The red vertical line indicates c = 0.1.(b) Information dimension of the scalar.Line styles and the vertical line are the same as (a).

FIG. 9 .
FIG. 9. Comparison between simulations and models based on FTLE.Line styles for panels (a) and (b): solid curves: fastest reacting case; dashed curves: slowest reacting case; black curves: DNS results; red curves: modeled from FTLE.(a) c against time.(b) c 2 against c .(c) Modeled c for all realizations at T = 7.5.

FIG. 10 .
FIG. 10.(a) Correlation between FT L E w and c as a function of time T for different Pe.Black: Pe = 10 000; blue: Pe = 1000: green: Pe = 100; and red: Pe = 10.The scattered points are at T = 7.5, marking the time the largest variation appears in the c field.The top inset plot shows the largest gap of c between extreme realizations at any given time as a function of Pe.The bottom inset plot shows c for Pe = 10 at T = 7.5.(b) Correlation between FT L E w and c as a function of time T for different η.Black: η = 4; green: η = 2; blue: η = 1; red: η = 1/2: and magenta: η = 1/4.The scattered points are at the respective times when the variability within each c field is the largest.This time scale is inversely proportional to the reaction rate η.The simulations for different η are only up to when c gets close to 1 everywhere, hence the correlation beyond this time is extrapolated as the last value of each case (where all curves have started to show asymptotic behavior).The top inset plot shows the largest gap of c between extreme realizations at any given time as a function of η.The bottom inset plot shows c for η = 4 at T = 1.875.