Dispersion and reaction in random flows: single realisation vs ensemble average

TL;DR

This paper investigates the dispersion of passive scalars in random, spatially periodic flows, contrasting single realisation and ensemble-averaged descriptions, and explores implications for reactive scalars modeled by FKPP equations.

Contribution

It demonstrates that while effective diffusivities are the same for single realisations and ensemble averages, their large-deviation rate functions differ, affecting tail predictions and reactive front behavior.

Findings

Single-realisation and ensemble-average diffusivities are identical.

Rate functions differ, with ensemble overestimating tail concentrations.

Reactive scalar fronts are well-described by single-realisation rate functions.

Abstract

We examine the dispersion of a passive scalar released in an incompressible fluid flow in an unbounded domain. The flow is assumed to be spatially periodic, with zero spatial average, and random in time, in the manner of the random-phase alternating sine flow which we use as an exemplar. In the long-time limit, the scalar concentration takes the same, predictable form for almost all realisations of the flow, with a Gaussian core characterised by an effective diffusivity, and large-deviation tails characterised by a rate function (which can be evaluated by computing the largest Lyapunov exponent of a family of random-in-time partial differential equations). We contrast this single-realisation description with that which applies to the average of the concentration over an ensemble of flow realisations. We show that the single-realisation and ensemble-average effective diffusivities are identical but that the corresponding rate functions are not, and that the ensemble-averaged description overestimates the concentration in the tails compared with that obtained for single-flow realisations. This difference has a marked impact for scalars reacting according to the Fisher--Kolmogorov--Petrovskii--Piskunov (FKPP) model. Such scalars form an expanding front whose shape is approximately independent of the flow realisation and can be deduced from the single-realisation large-deviation rate function. We test our predictions against numerical simulations of the alternating sine flow.

Equations

$${\partial }_{t}C+\mathbf{u}\cdot \nabla C=\kappa {\nabla }^{2}C,$$

$$\mathrm{d}\mathbf{X}=\mathbf{u}\left(\mathbf{X}(t),t\right)\mathrm{d}t+\sqrt{2\kappa }\mathrm{d}\mathbf{W},$$

$$C\left(\mathbf{x},t\right)=\mathbb{E}\left[\delta (\mathbf{X}\left(t\right)-\mathbf{x}\right)],$$

$$⟨C\left(\mathbf{x},t\right)⟩=⟨\mathbb{E}\left[\delta (\mathbf{X}\left(t\right)-\mathbf{x}\right]⟩,$$

$${\partial }_{t}C=\{\begin{array}{ll}-4a\sin \left(y+{\phi }_n\right){\partial }_{x}C& \text{for }t\in \left[n,n+1/4\right)\\ -4a\sin \left(x+{\psi }_n\right){\partial }_{y}C& \text{for }t\in \left[n+1/4,n+1/2\right)\\ 2\kappa {\nabla }^{2}C& \text{for }t\in \left[n+1/2,n+1\right)\end{array}$$

$$C_{n+\frac{1}{4}}\left(x,y\right)$$

$$C_{n+\frac{1}{2}}\left(x,y\right)$$

$$C_{n+1}$$

$${\partial }_{t}C+\mathbf{v}\cdot \nabla C=\nabla \cdot \left(\mathrm{K}\cdot \nabla C\right),$$

$$\mathbf{v}$$

$$\mathrm{K}$$

$$\overline{\mathbf{v}}$$

$$\overline{\mathrm{K}}$$

$$\mathrm{K}=\overline{\mathrm{K}}+\underset{t\to \infty }{\lim }\frac{⟨\mathbb{E}\left[\mathbf{X}\left(t\right)\right]⟩\otimes ⟨\mathbb{E}\left[\mathbf{X}\left(t\right)\right]⟩-⟨\mathbb{E}\left[\mathbf{X}\left(t\right)\right]\otimes \mathbb{E}\left[\mathbf{X}\left(t\right)\right]⟩}{2t}.$$

$$\mathbf{v}=\overline{\mathbf{v}}=0$$

$$\mathrm{K}=\overline{\mathrm{K}}-\underset{t\to \infty }{\lim }\frac{⟨\mathbb{E}\left[\mathbf{X}\left(t\right)\right]\otimes \mathbb{E}\left[\mathbf{X}\left(t\right)\right]⟩}{2t}$$

$$\mathbf{\chi }\left(\mathbf{x},t\right)=\mathbb{E}\left[\mathbf{X}\left(t\right)|\mathbf{X}\left(0\right)=\mathbf{x}\right]-\mathbf{x},$$

$${\partial }_t(\mathbf{\chi }+\mathbf{x})-\left(\mathbf{u}\cdot \nabla \right)(\mathbf{\chi }+\mathbf{x})=\kappa \Delta (\mathbf{\chi }+\mathbf{x}),$$

$${\partial }_t\mathbf{\chi }-\left(\mathbf{u}\cdot \nabla \right)\mathbf{\chi }=\kappa \Delta \mathbf{\chi }+\mathbf{u},$$

$$\sigma \left(t\right)=⟨\mathcal{A}\left[|\mathbf{\chi }{|}^2\left(\mathbf{x},t\right)\right]⟩,$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}t}=-2\kappa ⟨\mathcal{A}\left[\Vert \Vert \nabla \mathbf{\chi }\Vert {\Vert }^2\right]⟩+2⟨\mathcal{A}\left[\mathbf{u}\cdot \mathbf{\chi }\right]⟩$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d}t}\le -2\kappa \sigma +2U\sqrt{\sigma },$$

$$\sigma \left(t\right)\le \frac{U^2}{{\kappa }^2}{\left(1-{\mathrm{e}}^{-\kappa t}\right)}^2\le \frac{U^2}{{\kappa }^2}.$$

$${\mathbf{v}}_n\equiv \frac{\mathbb{E}\left[{\mathbf{X}}_n\right]}{n}\text{,}{\overline{\mathbf{v}}}_n\equiv \frac{⟨\mathbb{E}\left[{\mathbf{X}}_n\right]⟩}{n},$$

$${\mathrm{K}}_n\equiv$$

$${\overline{\mathrm{K}}}_n\equiv$$

$${\mathbf{X}}_{n+\frac{1}{4}}$$

$${\mathbf{X}}_{n+\frac{1}{2}}$$

$${\mathbf{X}}_{n+1}$$

$$C_n\left(\mathbf{x}\right)=\mathbb{E}\left[\delta \left({\mathbf{X}}_n-\mathbf{x}\right)\right],$$

$$⟨C_n⟩=\frac{1}{{\left(2\pi \right)}^{2n}}\underset{\left(\mathbf{\phi },\mathbf{\psi }\right)\in {\left[0,2\pi \right]}^n}{\text{idotsint}}\mathrm{d}\mathbf{\phi }\mathrm{d}\mathbf{\psi }{C}_n,$$

$${⟨C_n}^{{}^{(1)}}\left(x,y\right)⟩=⟨C_n\left(x-a\sin {\phi }_n,y\right)⟩,$$

$${⟨C_n}^{{}^{(2)}}\left(x,y\right)⟩={⟨C_n}^{{}^{(1)}}\left(x,y-a\sin {\psi }_n\right)⟩.$$

$$⟨C_n\left(\mathbf{x}\right)⟩=⟨\mathbb{E}\left[\delta \left({\tilde{\mathbf{X}}}_n-\mathbf{x}\right)\right]⟩$$

Full text

††preprint: APS/123-QED

I Introduction

We are interested in the dispersion of a passive scalar field under the combined action of advection by a prescribed flow and molecular diffusion. The scalar concentration $C\left(\bm{x},t\right)$ is governed by the advection–diffusion equation

[TABLE]

where $\kappa$ is the molecular diffusivity and $\bm{u}\left(\bm{x},t\right)$ a velocity field satisfying $\nabla\cdot\bm{u}=0$ as appropriate for incompressible flows. In the Lagrangian interpretation of (1), the positions ${\bm{X}}(t)$ of scalar particles obey the stochastic differential equation

[TABLE]

where ${\bm{W}}$ is a multidimensional Brownian motion. When suitably normalised, the concentration in (1) is identified with the probability density function

[TABLE]

of the particle positions, where $\mathbb{E}$ denotes expectation over realisations of the Brownian motion conditioned on ${\bm{X}(0)}={\bm{0}}$ so that $C({\bm{x}},0)=\delta({\bm{x}})$.

We focus on random flows, taking ${\bm{u}}({\bm{x}},t)$ to be a random field with given stationary statistics, as is standard when modelling turbulent dispersion G. I. Taylor (1922); Kraichnan (1966); A. J. Majda and P. R. Kramer (1999); Falkovich et al. (2001). This introduces an ensemble of flow realisations, additional to the ensemble of Brownian motion, with an associated average which we denote by $\mathbin{\left\langle\cdot\right\rangle}$. The averaged scalar concentration, defined as

[TABLE]

is often the main object of investigation (e.g. A. J. Majda and P. R. Kramer, 1999). However, for practical purposes, when a single flow realisation is experienced, it is the single-realisation dynamics of $C$ that matters. For large times and/or on large spatial scales, ergodicity arguments can be invoked to identify $C$ with $\mathbin{\langle C\rangle}$. The relationship between the single-realisation $C$ and the ensemble-averaged $\mathbin{\langle C\rangle}$ is in fact rather subtle; it is the topic of this paper. In statistical-physics parlance, the fixed flow realisation associated with $C$ corresponds to a quenched disorder, while the sampling over realisations associated with $\mathbin{\langle C\rangle}$ corresponds to an annealed disorder.

We investigate the relationship between $C$ and $\mathbin{\langle C\rangle}$ in the context of the initial-value problem for (1) in an unbounded domain, corresponding to the instantaneous release of a localised patch of scalar. In the absence of advection (${\bm{u}}=0$), the concentration behaves asymptotically as the Gaussian $(4\pi\kappa t)^{-d/2}\mathbin{\mathrm{e}}^{-|\bm{x}|^{2}/\left(4\kappa t\right)}$, with $d$ the number of spatial dimensions. For a spatially periodic, time-independent ${\bm{u}}\not=0$, at a coarse-grained level, the core of the concentration is also Gaussian: the dynamics of $C$ is approximately diffusive, with a (tensorial) effective diffusivity which is computed through the solution of the cell problem of homogenisation theory (see, e.g., A. J. Majda and P. R. Kramer (1999); G. A. Pavliotis and A. M. Stuart (2007) and references therein). The tails, however, are non-Gaussian: they take a large-deviation form, characterised by a rate function which is deduced from the solution of an eigenvalue problem (generalising the cell problem of homogenisation) and encodes all the cumulants of the concentration P. H. Haynes and J. Vanneste (2014). These conclusions extend naturally to flows that are also periodic in time.

In this paper, we examine the case of spatially-periodic, random-in-time flows. This class includes the widely-studied random alternating sine flow introduced by Pierrehumbert (1994) which we use as our testbed. For these flows, the coarse-grained picture of a Gaussian core and large-deviation tails continues to hold, with a single effective diffusivity and a single rate function describing the concentration $C$ in almost all realisations of the flow. On the other hand, the ensemble-averaged concentration $\mathbin{\langle C\rangle}$ also has a Gaussian core and large deviations tails. We show that, under the assumption of a vanishing spatial average of ${\bm{u}}$, the Gaussian cores of $C$ and $\mathbin{\langle C\rangle}$ coincide; in other words, the single-realisation and ensemble-average effective diffusivities are identical. The large-deviation tails of $C$ and $\mathbin{\langle C\rangle}$ are generally different, however, and we show that $\mathbin{\langle C\rangle}$ overestimates the concentration in the tails. We illustrate these results using numerical simulations of the alternating sine flow.

The differences in tail behaviour take a crucial importance for reacting scalars. When a logistic term is added to (1), leading to an FKPP model, the scalar dynamics is characterised by reaction fronts whose speed can be determined from the large-deviation rate function characterising the dispersion of the non-reacting scalar J. Gärtner and M. I. Freidlin (1979); Friedlin (1985). We demonstrate the value of this observation in the context of flows that are random in time by predicting the location of a chemical front propagating in the alternating sine flow and checking the prediction against a numerical simulation.

The paper is organised as follows. Section II motivates the problem considered by introducing the random alternating sine flow and presenting numerical results which illustrate its action on a passive scalar. Section III establishes the equivalence of $C$ and $\mathbin{\langle C\rangle}$ in the diffusive approximation that applies to the core of the scalar concentration when the spatial average of ${\bm{u}}$ vanishes. We check this prediction in the case of the alternating sine flow for which the effective diffusivity can be computed analytically by exploiting its ensemble-average interpretation. We consider the large-deviation aspects in section IV. We show there how the single-realisation rate function can be obtained by computing the largest Lyapunov exponent of a family of partial differential equations with coefficients which are random-in-time and periodic in space. We carry out this computation for the alternating sine flow and show that the single-realisation rate function so obtained differs from the ensemble-average one (which is given in closed form up to a Legendre transform). We highlight the importance of our findings for reacting scalars by examining front propagation in the FKPP model with the random alternating sine flow in section V. We conclude in section VI.

II A simple model of random dispersion: the alternating sine model

To make the problem concrete, we start by introducing a simple model of random dispersion. We consider the 2D alternating sine flow with random phases popularised by Pierrehumbert (1994) and still widely studied (see e.g. Haynes and Vanneste (2005); Meunier and Villermaux (2010), among others). This is one of the simplest random flows that efficiently stretches and folds material lines, making it uniformly mixing. This feature is crucial to enhance dispersion over that achieved by the sole action of molecular diffusion. The model simplifies both analytical and numerical computations by treating advection and diffusion in alternate steps in a manner akin to standard time-splitting numerical methods A. J. Valocchi and M. Malmstead (1992).

The velocity field of this alternating sine model is spatially periodic and random in time. The scalar field obeys the evolution equation

[TABLE]

where $n\in\mathbb{N}$ labels the iterations, $a$ is the maximum excursion distance in $x$ and $y$ over an iteration, $\left(\varphi_{n},\psi_{n}\right)$ are random phases drawn independently from a uniform distribution in $\left[0,2\pi\right]$, and $\kappa$ is the molecular diffusivity. In (5), time and space have been non-dimensionalised to set the time of each iteration – in effect the correlation time of the velocity – to 1 and the spatial period to $2\pi$.

The time-splitting enables a direct stepwise integration, turning the dynamics into that of a map. Let us introduce $C_{n}\left(\bm{x}\right)=C\left(\bm{x},t=n\right)$, with $\bm{x}=(x,y)$. The solutions of the first two steps in (5) are obtained as $C(x,y,t)=C_{n}(x-4a\sin(y+\varphi_{n})(t-n),y)$ and $C(x,y,t)=C_{n+\frac{1}{4}}(x,y-4a\sin(x+\psi_{n})(t-n-1/4))$ by integrating the advection equations. The solution of the diffusion equation in the third step is simply a convolution of $C_{n+\frac{1}{2}}(\bm{x})$ with a Gaussian diffusion kernel. This gives the 3-step map

[TABLE]

where $\mathcal{K}(\bm{x})=(4\pi\kappa)^{-1}\mathbin{\mathrm{e}}^{-|\bm{x}|^{2}/(4\kappa)}$ and $*$ denotes the convolution product. Numerically, the iterations (6) are carried out by discretising $C$ over a uniform square grid, with the two advection steps (6a)–(6b) performed using a bi-linear interpolation and the convolution of the diffusion step (6c) performed in Fourier space.

Figure 1 illustrates the dynamics (6) by showing the concentration field $C_{n}$ and the intermediate fields $C_{n+\frac{1}{4}}$ and $C_{n+\frac{1}{2}}$ for $n=0,1,2$, starting from a diamond-shaped scalar patch. It can be contrasted with figure 2 which shows the same fields but averaged over $1000$ realisations of the velocity field. The single-realisations and ensemble-averaged concentration are quite different in these early iterates of the alternating sine flow, but our focus is on the long-time (large-$n$) behaviour when the scalar has spread over many flow periods and is well represented by a coarse-grained description. In this limit, $C_{n}$ and $\mathbin{\langle C_{n}\rangle}$ are related in a manner that we elucidate in the rest of the paper.

We emphasise that we consider scalar dynamics in an unbounded domain. In a bounded domain, achieved e.g. by imposing periodic boundary conditions for $C_{n}$, the evolution is radically different, with a scalar relaxing exponentially fast towards a homogeneous state and adopting the form of a strange eigenmode Pierrehumbert (1994); Tsang et al. (2005); P. H. Haynes and J. Vanneste (2005).

III Diffusive approximation

In the long-time limit, when the passive scalar has sampled the velocity field well, the advection–diffusion equation (1) is often approximated by a diffusion equation

[TABLE]

where ${\bm{v}}$ is a uniform effective velocity and $\mathrm{K}$ a uniform effective-diffusivity tensor (A. J. Majda and P. R. Kramer, 1999; G. A. Pavliotis and A. M. Stuart, 2007). Correspondingly, the concentration field takes a Gaussian form. This is a manifestation of the central limit theorem which applies at distances that are $O(\sqrt{t})$ away from the centre of mass. The diffusive approximation implies the asymptotic linear growth with time of the mean and covariance of the position of particles obeying (2). The effective velocity and diffusivity

[TABLE]

where $\otimes$ denotes the tensor product, then characterise the dispersion for a specific realisation of the velocity field. We expect $\bm{v}$ and $\mathrm{K}$ to be identical for almost all realisations of the velocity field: this is the property of self-averaging which arises provided that the ensemble of velocity field is well sampled in the long-time limit. When if holds, it is clear that $\bm{v}=\mathbin{\langle\bm{v}\rangle}$ and $\mathrm{K}=\mathbin{\langle\mathrm{K}\rangle}$.

The property of self-averaging has been explored, mostly for media that are random in space; though delicate to establish rigorously, it appears to be generic in the absence of anomalous diffusion (e.g. J.-P. and Georges, 1990). In section IV we give an argument for flows that are random in time which indicates that self-averaging holds generically not only for $\bm{v}$ and $\mathrm{K}$ but more generally for all the cumulants of $C$ divided by $t$. Numerical results demonstrate this to apply to the alternating sine model.

The ensemble averaged concentration $\mathbin{\left\langle C\right\rangle}$ is also approximately Gaussian, with corresponding effective velocity and diffusivity defined as

[TABLE]

With the self averaging property, we readily obtain that $\bm{v}=\overline{\bm{v}}$ and

[TABLE]

It is not obvious whether $\mathrm{K}=\overline{\mathrm{K}}$ or not. The difference depends on the size of the fluctuations in the position of the centre of mass $\mathbb{E}\left[\bm{X}\left(t\right)\right]$ between realisations of the flow (see Le Doussal and Machta (1989); J.-P. and Georges (1990) for a discussion in the context of spatially random media]). In the next subsections, we establish that the equality holds for spatially periodic flows provided that the spatial average of $\bm{u}$ vanishes, and we verify this conclusion for the alternating sine model. This equivalence between single-realisation and ensemble-averaged diagnostics of dispersion does not generalise to high-order cumulants as we discuss in section IV.

III.1 Equivalence for spatially periodic random flows with vanishing mean

We restrict our attention to flows $\bm{u}$ that are random in time, with stationary statistics, zero mean, and finite correlation time, and $2\pi$-periodic in space with vanishing spatial average. We show that flows in this class satisfy $\mathrm{K}=\overline{\rm K}$, hence that the single-realisation and ensemble-averaged scalar fields $C$ and $\mathbin{\left\langle C\right\rangle}$ are equivalent in the diffusive approximation. We start by averaging (2) with respect to both the Brownian motion and the velocity realisation. This gives $\mathbin{\langle\mathbb{E}\left[\bm{\bm{X}\left(t\right)}\right]\rangle}=\bm{0}$, since $\mathbin{\langle\bm{u}(\cdot,t)\rangle}=0$, hence

[TABLE]

and

[TABLE]

using (10). We now show that $\left\|\mathbin{\langle\mathbb{E}\left[\bm{X}\left(t\right)\right]\otimes\mathbb{E}\left[\bm{X}\left(t\right)\right]\rangle}\right\|=o\left(t\right)$.

In order to estimate $\mathbb{E}\left[\bm{X}\left(t\right)\right]$, we introduce the expected displacement at time $t$ of a particle initially located at $\bm{x}$,

[TABLE]

such that \mathbb{E}\left[\bm{X}\left(t\right)\right]={\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\left(0,t\right). The expectation in (13) satisfies the backward Kolmogorov equation (e.g. Øksendal, 1998, Ch. 8)

[TABLE]

hence $\textstyle\bm{\chi}$ satisfies the corresponding forced equation

[TABLE]

with the initial condition {\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\left(\bm{x},0\right)=\bm{0}. Since $\bm{u}$ and the initial condition {\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\left(\bm{x},0\right) are both spatially periodic, the solution $\textstyle\bm{\chi}$ of Eq. (15) is spatially periodic. The cell-integrated variance

[TABLE]

where $\mathcal{A}$ denotes spatial averaging over a periodic cell, is found from (15) to satisfy

[TABLE]

using integration by parts and the incompressibility condition $\nabla\cdot\bm{u}=0$ and where $\|\|.\|\|$ denotes the Frobenius norm. We bound the right-hand side of (17) in terms of $\sigma$. Averaging Eq. (15) and using that $\mathcal{A}\left[\bm{u}\right]=0$ gives \mathcal{A}\left[{\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\right]=0. This enables us to use the Poincaré inequality for the term \mathcal{A}\left[\|\|\nabla{\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\|\|^{2}\right] and the Cauchy–Schwarz inequality for the term \mathcal{A}\left[\bm{u}\cdot{\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}\right]. This yields the differential inequality

[TABLE]

where the domain is taken as $[0,2\pi]^{d}$ so that the Poincaré constant is $1$ and $U=\mathbin{\langle\mathcal{A}[|{\bm{u}}|^{2}]\rangle}^{1/2}$ is the root-mean-square velocity. Integrating in time, we obtain

[TABLE]

The Laplacian operator in Eq. (15) ensures that the solution $\textstyle\bm{\chi}$ is smooth, so the boundedness of $\sigma$ guaranteed by (19) implies that \mathbin{\left\langle|{\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}|^{2}\left(0,t\right)\right\rangle} is also bounded. Finally, observing that \left\|\mathbin{\langle\mathbb{E}\left[\bm{X}\left(t\right)\right]\otimes\mathbb{E}\left[\bm{X}\left(t\right)\right]\rangle}\right\|=O\left(\mathbin{\left\langle|{\mathchoice{\raisebox{0.0pt}{\displaystyle\bm{\chi}}}{\raisebox{0.0pt}{\textstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptstyle\bm{\chi}}}{\raisebox{0.0pt}{\scriptscriptstyle\bm{\chi}}}}|^{2}\left(0,t\right)\right\rangle}\right)=O(1)=o\left(t\right) we conclude that $\mathrm{K}=\overline{\mathrm{K}}$.

We note that the upper bound $U^{2}/\kappa^{2}$ found in (19) is far from optimal. We also note that the argument above relies crucially on the assumption $\mathcal{A}\left[\bm{u}\right]={\bf 0}$. Without this premise, it is possible that $\mathrm{K}\neq\overline{\mathrm{K}}$. As an example, consider the case of a spatially uniform 2D flow with constant magnitude $U$ whose direction changes randomly every unit time. Such a flow has a non-zero spatial average but its ensemble average vanishes. Computing the two diffusivity tensors is straightforward and gives $\mathrm{K}=\kappa\mathbb{I}_{d}$ and $\overline{\mathrm{K}}=\left(\kappa+U/4\right)\mathbb{I}_{d}$, where $\mathbb{I}_{d}$ is the $d\times d$ identity matrix.

III.2 Application to the alternating sine model

In this subsection, we compare the diffusive approximation of $C_{n}$ and $\mathbin{\left\langle C_{n}\right\rangle}$ for the alternating sine model introduced in section II. We introduce the finite-time velocities

[TABLE]

and the finite-time diffusivity tensors

[TABLE]

which converge to $\bm{v}$, $\overline{\bm{v}}$, $\mathrm{K}$ and $\overline{\rm K}$ in the large-$n$ limit.

The discrete-time equivalent of (2) associated with the map (6) is the random walk for the position $\bm{X}_{n}=\left(X_{n},Y_{n}\right)$ of a single scalar particle,

[TABLE]

where the $\bm{W}_{n}$ are independent Gaussian variables with zero mean and unit variance. Starting from $\bm{X}_{0}=\bm{0}$, the scalar field $C_{n}$ satisfying (6), once suitably normalised, is the probability density function

[TABLE]

of the particle position, with $\mathbb{E}$ denoting the expectation with respect to the Gaussian variables $\bm{W}_{n}$ associated to the molecular diffusion. Taking ${\bm{X}}_{0}={\bm{0}}$ corresponds to a Dirac initial condition $C_{0}\left(\bm{x}\right)=\delta\left(\bm{x}\right)$, which we employ in what follows.

For the alternating sine flow, each realisation is defined by a sequence of phases $(\varphi_{n},\psi_{n})$, so the ensemble average is

[TABLE]

where $\bm{\varphi}=\left(\varphi_{0},\dots,\varphi_{n-1}\right)$ and $\bm{\psi}=\left(\psi_{0},\dots,\psi_{n-1}\right)$. We note that

[TABLE]

which is readily established by the substitution $\phi_{n}\mapsto\phi_{n}-y$ in the integral (24) defining the ensemble average. Similarly,

[TABLE]

As a result, $\mathbin{\left\langle C_{n}\right\rangle}$ can be obtained as the ensemble-averaged probability density function

[TABLE]

associated with the modified random walk

[TABLE]

This has the advantage over (22) that the position after $n$ steps can be written explicitly as

[TABLE]

where $\bm{W}$ is a single Gaussian variable with zero mean and unit variance.

Replacing ${\bm{X}}_{n}$ by $\tilde{\bm{X}}_{n}$ leaves the ensemble-averaged statistics unchanged, as Eq. (27) shows, while making their computation straightforward. We can in particular compute $\overline{\bm{v}}_{n}$ and $\overline{\rm K}_{n}$ by substituting $\tilde{\bm{X}}_{n}$ for ${\bm{X}}_{n}$ in (20) and (21) to obtain

[TABLE]

An argument similar to that in section III.1 but adapted the discrete-time setup can be used to show that ${\bm{v}}=\overline{\bm{v}}$ and $\mathrm{K}=\overline{\mathrm{K}}$ so that $C_{n}$ and $\mathbin{\left\langle C_{n}\right\rangle}$ shares the same diffusive approximation. We do not reproduce this argument here. Instead, we confirm its conclusion numerically by sampling the random walk (22) numerically, estimating $\bm{v}_{n}$ and $\mathrm{K}_{n}$, and verifying their convergence towards $\overline{\bm{v}}$ and $\overline{\mathrm{K}}$ in (30) as shown in Figure 3. This illustrates the equivalence between single-realisation and ensemble-averaged predictions in the diffusive approximation when only the first two cumulants matter. We next investigate the behaviour of higher-order cumulants by considering the large-deviation approximation of the concentration.

IV Large-deviation approximation

IV.1 General formulation

The effective velocity ${\bm{v}}$ and effective diffusivity $\mathrm{K}$ characterise only the core of the scalar distribution, occupying a $O(\sqrt{t})$ range around the centre of mass. At larger distances, the concentration field is generally not Gaussian; instead, it takes the large-deviation form

[TABLE]

valid for $|\bm{x}-\bm{v}t|=O(t)$. Here $g$ is the rate function (or Cramér function) and the symbol $\asymp$ denotes the asymptotic equivalence of the logarithms of the expressions on both sides (this ignores a constant multiple of $t^{-d/2}$ which can be determined by mass conservation). The approximation (31) holds for almost all realisations of the flow ${\bm{u}}$ for a broad class of flows. The Gaussian approximation is recovered from (31) by Taylor expanding $g$ around its minimum, $\bm{\xi}_{*}$ say, to find

[TABLE]

Eq. (31) goes beyond this and captures all the cumulants of the scalar function. Indeed, it can be shown that the scaled (single-realisation) cumulant generating function,

[TABLE]

is the Legendre transform of the rate function P. H. Haynes and J. Vanneste (2014):

[TABLE]

The Legendre duality implies that $g$ is the Legendre transform of $f$. This provides a means to evaluate $g$ since $f$ can be computed as the Lyapunov exponent of a random-in-time partial differential equation. To see this, we note that the function

[TABLE]

satisfies the backward Kolmogorov equation (e.g. Øksendal, 1998, Ch. 8)

[TABLE]

In view of the initial condition $w(\bm{x},\bm{q},0)=\exp(\bm{q}\cdot\bm{x})$ and the spatial periodicity of ${\bm{u}}$, the solution of (36) can be sought in form

[TABLE]

where $\phi(\bm{x},\bm{q},t)$ has the same spatial periodicity as $\bm{u}$. Introducing this into Eq. (36), $\phi$ is found to satisfy

[TABLE]

In the limit $t\to\infty$, $\phi$ grows exponentially with $t$ at a rate – the Lyapunov exponent of (38) – that can be identified from Eqs. (33), (35) and (37) as $f(\bm{q})$. From the theory of Lyapunov exponents in finite Oseledec (1968) and infinite Ruelle (1982) dimensions, $f(\bm{q})$ can be expected to be independent of the realisation of the flow under some assumptions of ergodicity. Thus we can obtain the rate function $g$ governing the large-deviation statistics of the scalar concentration in (almost all) single realisations of the velocity field by solving (38) for $\phi$ with periodic boundary conditions and an arbitrary initial condition for a range of values of $\bm{q}$, identifying $f(\bm{q})$ as the growth rate of $\phi$, i.e., as

[TABLE]

for some norm $\|\cdot\|$, then carrying out a numerical Legendre transform to deduce $g$. We note that the independence of $f(\bm{q})$ on flow realisation implies that the scaled cumulant generating function and hence all the cumulants of $C$ themselves, when divided by $t$, are independent of the flow realisation; in other words, that they are self-averaging.

The large-deviation characterisation of the concentration $C(\bm{x},t)$ in single realisations of the velocity has a counterpart for the ensemble-averaged concentration $\mathbin{\langle C(\bm{x},t)\rangle}$. Specifically, there is an ensemble-average large-deviation rate function $\overline{g}$ such that

[TABLE]

Its Legendre dual is the ensemble-average cumulant generating function

[TABLE]

Note that its single-realisation counterpart can be rewritten as

[TABLE]

on ensemble averaging Eq. (33) since (33) applies to almost all realisations of ${\bm{u}}$. Thus the difference between the two generating functions stems from the lack of commutation between ensemble average and logarithm. Note that the definitions (41)–(42), the concavity of the log and Jensen’s inequality imply that $f\leq\overline{f}$, which guarantees that

[TABLE]

using the order-reversing property of the Legendre transform. Thus, the ensemble-averaged scalar field always overestimates the tail concentration in single realisations of the flow at large time.

A focus of this paper is the relationship between $f$ and $\overline{f}$ and, correspondingly, between $g$ and $\overline{g}$ The equalities $\bm{v}=\overline{\bm{v}}$ and $\mathrm{K}=\overline{\mathrm{K}}$ established in the previous sections imply that Taylor expansions at $\bm{\xi}=\bm{\xi}_{*}$ of $g$ and $\bar{g}$ are identical up to and including quadratic terms. The same is true for $f$ and $\overline{f}$ by virtue of the Legendre duality with $g$ and $\overline{g}$. The equalities do not however hold for higher order terms: in general, $g\not=\overline{g}$ and $f\not=\overline{f}$. Thus the single-realisation and ensemble-averaged concentration fields, which are identical in the diffusive approximation, differ, leading to differences in the cumulants of order greater than 2. Physically, this means that the tails of the concentration in single realisations, unlike the core, cannot be predicted on the basis of ensemble-averaged statistics. We next demonstrate this explicitly for the alternating sine flow.

IV.2 Application to the alternating sine model

For the alternating sine flow, we use the discrete version of the large-deviation rate functions $g$ and $\overline{g}$, such that

[TABLE]

and

[TABLE]

The ensemble-averaged $\overline{f}$ can be computed explicitly, replacing the positions $\bm{X}_{n}$ by $\tilde{\bm{X}}_{n}$ in Eq. (29), exploiting their statistical equivalence (27). This gives

[TABLE]

where ${\bm{q}}=(q_{x},q_{y})$, on using the integral definition of the modified Bessel function of the first kind $I_{0}$ (DLMF, , Eq. 10.32.1). Note that a Taylor expansion using that $I_{0}(aq_{x})\sim 1-a^{2}q_{x}^{2}/4$ as $q_{x}\to 0$ and similarly for $I_{0}(aq_{y})$ recovers the diffusive approximation (30).

Obtaining the single-realisation $f$ is not as straightforward: this requires estimating numerically the Lyapunov exponent of the discrete-in-time equivalent of (38), namely

[TABLE]

This can be integrated in time to obtain the 3-step discrete map

[TABLE]

where $\phi_{n}(x,y)=\phi(x,y,n)$, $\widetilde{\mathcal{K}}(\bm{x})=(4\pi\kappa)^{-1}\mathbin{\mathrm{e}}^{\bm{q}\cdot\bm{x}-|\bm{x}|^{2}/(4\kappa)}$ and $*$ denotes the convolution product.

For a given $\bm{q}$, we estimate $f(\bm{q})$ by iterating (48) numerically for $n$ up to $10^{5}$, starting from $\phi_{0}(x,y)=1$. The function $\phi_{n}(x,y)$ is discretised on a $128^{2}$ uniform spatial grid. The first two steps, (48a)–(48b), are performed in a semi-Lagrangian fashion, as a circular permutation of the spatial grid indices; the third step, (48c), is performed spectrally, in Fourier space. The growth rate $f(\bm{q})$ is approximated as $n^{-1}\log\|\phi_{n}\|$ for $n$ large enough that the approximation has converged. The convergence is rather slow, with $n=10^{5}$ iterations required here. We carry out this numerical estimation of $f(\bm{q})$ for $100^{2}$ values of $(q_{x},q_{y})$ linearly spaced in $[0,1]^{2}$, using symmetries to obtain $f({\bm{q}})$ in $[-1,1]^{2}$.

We show the results for the rate functions $g$ and $\overline{g}$, obtained from $f$ and $\overline{f}$ by numerical Legendre transform, in Fig.  4. The two functions are markedly different away from the centre ${\bm{\xi}}={\bm{\xi}}_{*}={\bm{0}}$ where the diffusion approximation applies and around which they are approximately equal. The function $g$ grows substantially faster than $\overline{g}$ (cf. Eq. (43)). Consequently, the ensemble-averaged $\mathbin{\langle C_{n}\rangle}$ considerably overestimates the tail concentration in single realisations of the flow. The departure from axis-symmetry, present for both $g$ and $\overline{g}$, is the result of the flow geometry: fast dispersion stems from particles experiencing maximum advection at each step, which can lead to the maximum displacements $(\pm a,\pm a)$ at each time step (ignoring diffusion). This indicates that for sufficiently large ${\bm{\xi}}$ (and sufficiently small $\kappa$), the contours of $g$ and $\overline{g}$ are approximately square as would be the case for advective displacements randomly chosen from $(\pm a,\pm a)$ at each iteration (cf. A. Tzella and J. Vanneste, 2016). This effect is much more marked for the ensemble-averaged concentration than for single realisations.

The broad conclusion is that the ensemble-averaged concentration is a poor predictor of the behaviour of dispersing passive scalars when it comes to their low-concentration tails. While the tail behaviour might seem of marginal interest, it plays a crucial role in the presence of chemical reactions, in particular by determining the dynamics of some reaction fronts. In this instance, the differences between single-realisation and ensemble-averaged dynamics can be dramatic as we discuss next.

V Application to reactive front propagation

V.1 FKPP model

We now consider the dynamics of a reacting scalar governed by the classical FKPP model of Fisher Fisher (1937) and Kolmogorov et al. Kolmogorov et al. (1937) which, in the presence of advection, reads

[TABLE]

where $\gamma>0$ is the reaction rate and the concentration $C$ has been normalised so that its saturation value be $1$. According to this model, in an infinite domain, an initially localised patch of scalar spreads behind a reaction front.

For ${\bm{u}}=0$ and for $t\gg 1$, the front is circular in two dimensions (spherical in three) with a radius growing at the constant speed $s=2\sqrt{\kappa\gamma}$. For a spatially periodic flow, the front can be approximated by a self-similar convex curve (or surface) defined by

[TABLE]

where $g$ is the large-deviation rate function introduced in (31) for the passive scalar problem (1) (i.e. $\gamma=0$ in (49)) (J. Gärtner and M. I. Freidlin, 1979; Friedlin, 1985; Xin, 2000; A. Tzella and J. Vanneste, 2014; A. Tzella and J. Vanneste, 2015). The remarkable connection between the large-deviation tail of a passive scalar and the FKPP reaction front embodied by (50) can be understood heuristically as follows. The front’s propagation is controlled by its leading edge, where $C\ll 1$ so that (49) can be linearised about $C=0$; the solution is then that of the passive-scalar problem (1), approximately $\exp\left(-tg({\bm{x}}/t)\right)$, multiplied by the factor $\exp(\gamma t)$ associated with the (linearised) reaction, hence

[TABLE]

Eq. (50) follows by observing that the front marks the transition between exponential growth and exponential decay.

For the random-in-time flows considered in this paper, $g$ is clearly the single-realisation rate function, so the shape of the front can be determined from (50) using that $g$ is the Legendre transform of the Lyapunov $f$ exponent of the random partial differential equation (38). The ensemble-averaged $\overline{g}$, while easier to evaluate, is only relevant to the front problem in the limit of small reaction rate $\gamma\ll 1$, when the diffusive approximation applies, $g\approx\overline{g}$, and the front is the ellipse (or ellipsoid)

[TABLE]

as follows from Taylor expanding (50).

V.2 Reacting alternating sine model

We illustrate how (50) predicts the shape of a reacting front by considering a reacting version of the alternating sine model (6). This adds the reaction step

[TABLE]

where $C_{n+\frac{3}{4}}=\mathcal{K}*C_{n+\frac{1}{2}}$, to the advective and diffusive steps in (6), corresponding to the time-$1$ map associated with the FKPP reaction $\partial_{t}C=\gamma C(1-C)$. The front location is determined by the discrete-time equivalent of (50), that is,

[TABLE]

Note that this equation makes it possible to interpret the contours of $g$ in Fig. 4 as giving the shape of fronts for different values of $\gamma$.

Figure 5 shows the concentration field obtained by direct numerical integration of the reacting alternating sine model starting from a localised patch with $a=\pi$, $\kappa=0.1$ and $\gamma=2$. The result is compared with the prediction (54) for the front location. The prediction which substitutes the ensemble averaged $\overline{g}$ (computed in subsection IV.2) is also shown. The figure confirms the validity of (54) and illustrates the inadequacy of the prediction based on $\overline{g}$.

VI Summary and conclusion

This paper investigates the large-time statistics of passive scalars dispersing under the combined action of molecular diffusion and advection by single realisations of flows taken to be spatially periodic and random-in-time. The focus is on the instantaneous release of a localised patch of scalar in an unbounded domain. After some time, the scalar patch extends over many flow periods and has sampled flow correlation time such that it is well represented using a coarse-grained description, and the velocity $\bm{v}$ and diffusivity $\mathrm{K}$ are identical in almost all flow realisations (that is, self-averaging applies). This happens for times large compared to both the correlation time of the flow, and the typical time for dispersion across a periodic cell which depends on both molecular diffusivity value and flow strength.

The single-realisation statistics are compared to their ensemble-averaged counterparts which are more easily accessible to analysis and prominently studied in the literature. Under the assumption of a vanishing spatial average of the flow, we show that the Gaussian cores of the ensemble-average and single-realisation scalar distributions are identical. We introduce a large-deviation description which captures the tails of the concentration, hence higher-order cumulants, and show that these are in general overestimated by the ensemble-average prediction. In the single-realisation problem, we demonstrate that the large-deviation rate function can be obtained by computing the largest Lyapunov exponent for the spatially periodic solution of a family of partial differential equations with coefficients that are random-in-time. We illustrate these results by computing the single-realisation rate function for the widely-used alternating sine flow and comparing it with its ensemble-average counterpart which is obtained analytically.

The tail behaviour of the concentration of a dispersing passive scalar has a zeroth-order manifestation in the presence of a chemical reaction of FKPP type. This is because the large-time dynamics of the reaction front that is then formed is directly related to the large-deviation rate function characterising the dispersion. We consider this problem which further highlights the importance of the distinction between ensemble-average and single-realisation predictions for the scalar-concentration tails of the scalar distributions. For the random alternating sine flow, we successfully predict the location of the reaction front using the single-realisation rate function computed for the non-reacting scalar. The corresponding prediction based on the ensemble-average rate function fails.

In conclusion, we emphasise the existence of a hierarchy of behaviours for scalars dispersing in random flows: at the most basic level, the self-averaging of the cumulants scaled by $t^{-1}$ ($\bm{v}$, $\mathrm{K}$ and higher order) can be expected to hold for typical flows. Next, the equivalence between single-realisation and ensemble-average statistics in a diffusive approximation (that is, the equalities $\bm{v}=\overline{\bm{v}}$ and $\mathrm{K}=\overline{\mathrm{K}}$) is more delicate but applies to periodic flows with vanishing spatial average. We expect it to also hold for spatially random flows with finite correlation length and vanishing spatial average. Finally, the analogous equivalence for higher-order cumulants, which characterise the large-deviation approximation, does not hold, even for such a well-behaved model as the alternating sine flow. We also note that the ensemble average overestimates the value of the scaled cumulants of order 3 and higher; in physical terms, it overestimates the concentration in the tails of the scalar distribution.

We finally observe that the distinction between single-realisation and ensemble-average predictions of scalar concentrations is an issue of broad interest, arising whenever stochastic models are used to represent complex velocity fields. While ensemble-averaged statistics are usually more amenable to analysis, single-realisation statistics are often more meaningful for practical applications such as risk assessment. In view of this broad interest, it would be desirable to extend the analysis of the present paper to other classes of flows, and in particular to flows that are random in space Isichenko (1992). We also note that the problem posed here in unbounded domains has a counterpart in bounded domains, in which case the focus will be on the differences between the rate of decay of the scalar concentration in single realisation and in ensemble averages P. H. Haynes and J. Vanneste (2014).

Acknowledgements.

This research is funded by grant EP/S023291/1 of the UK Engineering and Physical Sciences Research Council. We thank Alexandra Tzella and the two anonymous referees for useful comments.