2.1. Eulerian and Lagrangian perspectives

Transport of a passive scalar in an underlying flow field $\boldsymbol{u}(\boldsymbol{x},t)$ , here assumed to be known deterministically, is described from an Eulerian perspective by the ADE,

(2.1a) \begin{align} \frac {\partial c}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u} c) &= \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \boldsymbol{D} \boldsymbol{\nabla }c \right ), \\[-9pt] \nonumber \end{align}

(2.1b) \begin{align} c(\boldsymbol{x},0) &=c_0(\boldsymbol{x}) , \\[9pt] \nonumber \end{align}

where $c(\boldsymbol{x},t)$ is the concentration at position $\boldsymbol{x}$ (in $d$ spatial dimensions) and $\boldsymbol{D}(\boldsymbol{x},t)$ is the scalar diffusion tensor. For a given initial condition $c_0(\boldsymbol{x})$ , the ADE can be related to the FPE for the PDF $f_{\boldsymbol{{X}}}(\boldsymbol{x};t) = c(\boldsymbol{x},t)/\int _{\mathbb{R}^d} \text{d}\boldsymbol{x} \, c(\boldsymbol{x},t)$ of the locations of the solute particles composing the plume (Risken Reference Risken1996; Noetinger et al. Reference Noetinger, Roubinet, Russian, Le Borgne, Delay, Dentz, De Dreuzy and Gouze2016):

(2.2a) \begin{align} \frac {\partial f_{\boldsymbol{{X}}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\mu }{\kern-1pt}f_{\boldsymbol{{X}}}) &= \boldsymbol{\nabla }\boldsymbol{\nabla} ^\top : \boldsymbol{D} f_{\boldsymbol{{X}}} , \\[-9pt] \nonumber \end{align}

(2.2b) \begin{align} f_{\boldsymbol{{X}}}(\boldsymbol{x};0) &= f^0_{\boldsymbol{{X}}}(\boldsymbol{x}). \\[9pt] \nonumber \end{align}

Here the drift is $\boldsymbol{\mu }(\boldsymbol{x},t)=\boldsymbol{u}(\boldsymbol{x},t)+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{D}(\boldsymbol{x},t)$ and : denotes the inner tensor product.

From a Lagrangian perspective, the equivalence of the FPE and Langevin dynamics (Risken Reference Risken1996) is used to describe the evolution of particles of solute (Pope Reference Pope1985). In this sense, an equivalent formulation to (2.2) in Lagrangian form is given by the stochastic Langevin (or Itō) equation for particle trajectories $\boldsymbol{{X}}(t)$ , i.e.

(2.3a) \begin{align} \text{d}\boldsymbol{{X}}(t) &= \boldsymbol{\mu }(\boldsymbol{{X}}(t),t)\text{d}t + \boldsymbol{\varSigma }(\boldsymbol{{X}}(t),t) \text{d} \boldsymbol{W}_t, \\[-9pt] \nonumber \end{align}

(2.3b) \begin{align} \boldsymbol{{X}}(0) &= \boldsymbol{X}_0, \\[9pt] \nonumber \end{align}

with standard Wiener increments $\text{d}\boldsymbol{W}_t$ , $2\boldsymbol{D}=\boldsymbol{\varSigma }\boldsymbol{\varSigma }^\top$ , and the PDF of $\boldsymbol{X}_0$ given by the Eulerian initial condition $f^0_{\boldsymbol{{X}}}(\boldsymbol{x})$ . The PDF of particle positions, which corresponds to the concentration field normalised by the total (conserved) mass, is obtained through a sampling process that can be formally expressed in the continuum limit as

(2.4) \begin{align} f_{\boldsymbol{{X}}}(\boldsymbol{x};t) = \overline {{\delta [\boldsymbol{x}-\boldsymbol{{X}}(t)]}}. \end{align}

Here, $\overline {(\boldsymbol{\cdot })}$ indicates the average over all trajectories computed through (2.3), i.e. averaging over the (diffusive) noise and over the variability in the initial condition, for each position $\boldsymbol{x}$ at a given time $t$ .

2.2. Augmented phase space

The framework proposed here is built upon the construction of a trajectory-based formulation for which the PDF $f_{\boldsymbol{X}}$ of the trajectories $\boldsymbol{X}(t)$ solves the FPE (2.2), but where any stochasticity is time independent. To this end, we write

(2.5a) \begin{align} \frac {\text{d}\boldsymbol{X}(t)}{\text{d}t} &= \boldsymbol{u}(\boldsymbol{X}(t),t) + \boldsymbol{s}(\boldsymbol{X}(t),\boldsymbol{\varXi },t), \\[-9pt] \nonumber \end{align}

(2.5b) \begin{align} \boldsymbol{X}(0) &= \boldsymbol{X}_0, \\[9pt] \nonumber \end{align}

with $\boldsymbol{\varXi }$ a vector of $d$ random coefficients with PDF $f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }):\mathbb{R}^d \rightarrow \mathbb{R}^+$ and $ \boldsymbol{s} (\boldsymbol{x},\boldsymbol{\xi },t): \mathbb{R}^d\times \mathbb{R}^d \times \mathbb{R}^+ \rightarrow \mathbb{R}^d$ a deterministic function of the random positions, random coefficients and time. Here, $d$ is the spatial dimensionality. The sole sources of uncertainty in (2.5a ) are the initial condition $\boldsymbol{X}_0$ and the random coefficients $\boldsymbol{\varXi }$ , and (2.5) is therefore a random differential equation (Soong Reference Soong1973), as opposed to a stochastic differential equation like (2.3): while (2.5) is amenable to conventional calculus, (2.3) requires stochastic calculus. The trajectories of the present framework do not represent physical trajectories and only equality of the one-point PDF (as opposed to multi-point statistics, such as temporal correlations along trajectories) will be enforced, i.e. $f_{\boldsymbol{X}} = f_{\boldsymbol{{X}}}$ should solve the FPE (2.2). We show in § 3 how to chose the source term $\boldsymbol{s}$ and the PDF $f_{\boldsymbol{\varXi }}$ of the random coefficients so as to achieve this.

First, it is convenient to study (2.5) for arbitrary $\boldsymbol{s}$ and $\boldsymbol{\varXi }$ . The time evolution of the joint PDF of the particle locations and random coefficients $f_{\boldsymbol{X}\boldsymbol{\varXi }}(\boldsymbol{x},\boldsymbol{\xi };t):\mathbb{R}^d\times \mathbb{R}^d \times \mathbb{R}^+ \rightarrow \mathbb{R}^+$ is governed by a deterministic hyperbolic Liouville equation, which in conservative form reads as

(2.6a) \begin{align} \frac {\partial f_{\boldsymbol{X}\boldsymbol{\varXi }}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot } ([ \boldsymbol{u}(\boldsymbol{x},t) + \boldsymbol{s}(\boldsymbol{x},\boldsymbol{\xi },t) ] f_{\boldsymbol{X}\boldsymbol{\varXi }} ) &= 0, \\[-9pt] \nonumber \end{align}

(2.6b) \begin{align} f_{\boldsymbol{X}\boldsymbol{\varXi }}(\boldsymbol{x},\boldsymbol{\xi };0) = f^0_{\boldsymbol{X}\boldsymbol{\varXi }}(\boldsymbol{x},\boldsymbol{\xi }) &= f^0_{\boldsymbol{X}}(\boldsymbol{x})f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }), \\[9pt] \nonumber \end{align}

where the coefficients play the role of parameters and at the initial time, the coefficients are independent and uncorrelated with the initial locations. The derivation of (2.6a ) from (2.5a ) is a direct application of the Liouville theorem (Moyal Reference Moyal1949). Alternatively, in Appendix A we provide a derivation using a distributional approach (Venturi & Karniadakis Reference Venturi and Karniadakis2012; Venturi et al. Reference Venturi, Tartakovsky, Tartakovsky and Karniadakis2013; Tartakovsky & Gremaud Reference Tartakovsky and Gremaud2015). After marginalising over the coefficients, the normalised concentration field is recovered:

(2.7) \begin{align} f_{\boldsymbol{X}}(\boldsymbol{x};t)=\int _{\mathbb{R}^d} \text{d}\boldsymbol{\xi } \, f_{\boldsymbol{X}\boldsymbol{\varXi }}(\boldsymbol{x},\boldsymbol{\xi };t). \end{align}

A strong advantage of this formulation relates to the hyperbolicity of the Liouville equation, which permits the use of the method of characteristics (Domínguez-Vázquez et al. Reference Domínguez-Vázquez, Castiblanco-Ballesteros, Jacobs and Tartakovsky2024a ). Points $(\boldsymbol{x},\boldsymbol{\xi })$ can be thought of elements of the augmented phase space spanned by the particle locations and random coefficient values. Let $\hat {\boldsymbol{x}}(t;\boldsymbol{\xi },\boldsymbol{x}_0,t_0)$ be a trajectory along a characteristic in this augmented phase space, parameterised by the values $\boldsymbol{\xi }$ of the random coefficients and the initial position $\boldsymbol{x}_0$ at time $t_0$ . We allow here for arbitrary $t_0$ because it is useful in discussing backward trajectories below. We omit $t_0$ when $t_0=0$ and we sometimes further omit $\boldsymbol{\xi }$ and $\boldsymbol{x}_0$ for brevity. We also introduce $\hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}(t; \boldsymbol{\xi },\boldsymbol{x}_0,t_0)=f_{\boldsymbol{X}\boldsymbol{\varXi }}(\hat {\boldsymbol{x}}(t;\boldsymbol{\xi },\boldsymbol{x}_0,t_0),\boldsymbol{\xi };t)$ , the joint PDF along the characteristics. Characteristic lines $\hat {\boldsymbol{x}}$ describe the path of what we call support particles, analogously to the more standard use of the term ‘particles’ in relation to advective–diffusive (Langevin) trajectories. These support particles correspond to tracking discrete points of the support of the joint PDF in the augmented phase space, carrying their joint PDF value $\hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}$ along their corresponding characteristic lines. We find the following system of characteristic lines to solve (2.6) in this picture:

(2.8a) \begin{align} \frac {\text{d}\hat {\boldsymbol{x}}(t)}{\text{d}t} &= \boldsymbol{u} \big(\hat {\boldsymbol{x}}(t),t \big) + \boldsymbol{s}\big(\hat {\boldsymbol{x}}(t),\boldsymbol{\xi },t \big), \\[-9pt] \nonumber \end{align}

(2.8b) \begin{align} \frac {\text{d}\hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}(t)}{\text{d}t} &= -\hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}(t) \boldsymbol{\nabla }\boldsymbol{\cdot } \big[ \boldsymbol{u}\big(\hat {\boldsymbol{x}}(t),t \big) +\boldsymbol{s}\big(\hat {\boldsymbol{x}}(t),\boldsymbol{\xi },t \big) \big], \\[-9pt] \nonumber \end{align}

(2.8c) \begin{align} \hat {\boldsymbol{x}}(t_0)&=\boldsymbol{x}_0, \quad \hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}(t_0) = \hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}^0 = f^0_{\boldsymbol{X}}(\boldsymbol{x}_0)f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }). \\[12pt] \nonumber \end{align}

This system of equations combines properties of the Fokker–Planck and Langevin equations, in the sense that we solve directly for the (joint) PDF, but in Lagrangian form along characteristics. Equations (2.6) or (2.8) define a crypto-deterministic (Whittaker Reference Whittaker1943) description of the problem of passive scalar transport, i.e. all the uncertainty is contained in the initial conditions in the augmented phase space.

To recover normalised concentrations, marginalisation of the joint PDF is needed with (2.7). Due to the relationship (2.8a ) between coefficients $\boldsymbol{\xi }$ and positions $\hat {\boldsymbol{x}}$ along characteristics, this marginalisation may also be performed directly in Lagrangian form. Indeed, we can write the following Lagrangian relation (Cramér Reference Cramér1946):

(2.9) \begin{align} &\hat {f}_{\boldsymbol{X}\boldsymbol{\varXi }}(t;\boldsymbol{\xi },\boldsymbol{x}_0) = J(t;\boldsymbol{\xi },\boldsymbol{x}_0)^{-1} f^0_{\boldsymbol{X}}(\boldsymbol{x}_0)f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }), &J(t;\boldsymbol{\xi },\boldsymbol{x}_0)=\left | \frac {\partial (\hat {\boldsymbol{x}},\boldsymbol{\xi })}{\partial (\boldsymbol{x}_0,\boldsymbol{\xi })} \right | = | \boldsymbol{\nabla} _{\boldsymbol{x}_0}\hat {\boldsymbol{x}}|. \end{align}

Here $J(t)$ is the determinant of the transformation Jacobian from the initial to final time. Note that, according to (2.8b ), $J(t)$ is unity if and only if $\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}+\boldsymbol{s})=0$ . Relation (2.9) offers an alternative to solving (2.8b ) based on a transformation of random variables between the initial and later time. Using this, the marginalisation (2.7) becomes

(2.10) \begin{align} f_{\boldsymbol{X}}(\boldsymbol{x};t)=\int _{\mathbb{R}^d} \text{d}\boldsymbol{\xi } \, J \big[t;\boldsymbol{\xi },\hat {\boldsymbol{x}}_{0}(\boldsymbol{\xi },\boldsymbol{x},t) \big]^{-1}f^0_{\boldsymbol{X}} \big[\hat {\boldsymbol{x}}_{0}(\boldsymbol{\xi },\boldsymbol{x},t) \big]f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }), \end{align}

where we have introduced the backward map $\hat {\boldsymbol{x}}_{0}(\boldsymbol{\xi },\boldsymbol{x},t)$ for $t_0=0$ . For arbitrary $t_0$ , the backward map $\hat {\boldsymbol{x}}_{0}(t_0;\boldsymbol{\xi },\boldsymbol{x},t)$ is the solution of the characteristics (2.8a ) given a specified final condition, i.e. for an initial time argument of $\hat {\boldsymbol{x}}$ (here $t$ ) larger than the target time along the characteristic (here $t_0$ ). Thus, we have the general equivalence $\hat {\boldsymbol{x}}_{0}(t_0;\boldsymbol{\xi },\boldsymbol{x},t) \equiv \hat {\boldsymbol{x}}(t_0;\boldsymbol{\xi },\boldsymbol{x},t)$ , but we introduce the backward notation for intuitive clarity and to allow us to omit $t_0=0$ for notational simplicity. In practice, the backward map can be computed in one of two ways: (i) by integrating (2.8a ) backward as just described; or (ii) given the forward solution, by inverting $\hat {\boldsymbol{x}}[t;\boldsymbol{\xi },\hat {\boldsymbol{x}}_0(t_0;\boldsymbol{\xi },\boldsymbol{x},t),t_0]=\boldsymbol{x}$ .

Since the sole sources of stochasticity in the description are the initial condition and the random coefficients, moments may be directly obtained from the characteristics by averaging over their PDFs (see Appendix B). Moreover, due to linearity, quantities associated with any initial condition can be constructed by summing over the solutions for deterministic (i.e. Dirac-delta) initial conditions, after integrating over the random coefficients for each initial point. We make explicit use of this fact in § 4.

3.1. Closed two-moment description

We are now ready to determine the source term $\boldsymbol{s}$ such that the Liouville equation (2.5) marginalised along the coefficients leads to the FPE (2.2). To this end, we derive conditions for $\boldsymbol{s}$ by comparing moments of the FPE with those of the system (2.5). Taking moments of the FPE will in general lead to a closure problem. Instead, we locally linearise the flow in order to obtain a closed system. Details on the derivation of a closed two-moment description of $\boldsymbol{X}(t)$ based on (2.5) are given in Appendix B. We focus here on the case of a spatially independent dispersion tensor, and we seek a simple form of $\boldsymbol{s}$ capable of matching the first two moments of the FPE (2.2). In this case, it turns out to be sufficient to take a noise term that is independent of position and linear in the random coefficient values, i.e.

(3.1) \begin{align} \boldsymbol{s}(\boldsymbol{\xi },t) = \boldsymbol{\varphi }(t)\boldsymbol{\xi }, \end{align}

where, by definition, $\boldsymbol{\varphi }(t):\mathbb{R}^+ \rightarrow \mathbb{R}^{d}\times \mathbb{R}^d$ is the matrix containing the first derivatives of $\boldsymbol{s}$ with respect to the coefficients $\boldsymbol{\xi }$ . We also assume for simplicity that the underlying flow field is steady (i.e. the Eulerian field $\boldsymbol{u}$ is independent of time $t$ ), two-dimensional ( $d=2$ ) and incompressible ( $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0$ ). Then, the results in Appendix B read

(3.2a) \begin{align} \frac {\text{d}\overline {\boldsymbol{X}}(t)}{\text{d}t} &= \boldsymbol{u} \big[\overline {\boldsymbol{X}}(t) \big] + \boldsymbol{\varphi }(t)\overline {\boldsymbol{\varXi }} , \\[-9pt] \nonumber \end{align}

(3.2b) \begin{align} \frac {\text{d}\boldsymbol{\kappa }(t)}{\text{d}t} &= \boldsymbol{\epsilon }(t)\boldsymbol{\kappa }(t) \ + \boldsymbol{\kappa }(t) \boldsymbol{\epsilon }(t)^\top + \boldsymbol{\varphi }(t)\boldsymbol{\zeta }(t) + [\boldsymbol{\varphi }(t)\boldsymbol{\zeta }(t)]^\top , \\[-9pt] \nonumber \end{align}

(3.2c) \begin{align} \frac {\text{d}\boldsymbol{\zeta }(t)}{\text{d}t} &= \boldsymbol{\zeta }(t) \boldsymbol{\epsilon }(t)^\top + \boldsymbol{\theta }\boldsymbol{\varphi }(t)^\top , \\[-9pt] \nonumber \end{align}

(3.2d) \begin{align} \overline {\boldsymbol{X}}(0) &= \boldsymbol{} \overline {\boldsymbol{X}}^0 , \quad \boldsymbol{\kappa }(0) = \boldsymbol{\kappa }^0 , \quad \boldsymbol{\zeta }(0) = \boldsymbol{0}. \\[9pt] \nonumber \end{align}

For each given initial position, this system involves the velocity gradient tensor at the average particle location,

(3.3) \begin{align} \boldsymbol{\epsilon }(t) = (\boldsymbol{\nabla }\boldsymbol{u} [\overline {\boldsymbol{X}}(t)])^\top , \end{align}

along with the first two moments of the phase space variables, i.e. the averages $\overline {\boldsymbol{X}}(t)$ and $\overline {\boldsymbol{\varXi }}$ and the correlations $\kappa _{\textit{ij}}(t)=\overline {X_i^\prime X_{\kern-1.5pt j}^\prime }(t)$ , $\zeta _{\textit{ij}}(t)=\overline {\varXi _i^\prime X_{\kern-1.5pt j}^\prime }(t)$ and $\theta _{\textit{ij}}=\overline {\varXi _i^\prime \varXi _j^\prime }$ . Inspired by previous results for Brownian motion and the Ornstein-Uhlenbeck process, see Domínguez-Vázquez et al. (Reference Domínguez-Vázquez, Jacobs and Tartakovsky2024b ) and Appendix C, we take the coefficients to follow independent standard normal distributions, so that $\overline {\boldsymbol{\varXi }}=\boldsymbol{0}$ and $\boldsymbol{\theta }$ is the identity matrix, and

(3.4) \begin{equation} f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }) = \frac {\exp \left (-\frac {1}{2} \sum _{i=1}^d \xi _i^2\right )}{(2\pi )^{d/2}}. \end{equation}

This turns out to be sufficient to capture the more general, linearised process examined here.

We now equate the results for the first two moments from (3.2a ) and (3.2b ) to the equivalent results for the FPE or the ADE; see, for example, De Barros et al. (Reference De Barros, Dentz, Koch and Nowak2012) or Aquino & Bolster (Reference Aquino and Bolster2017) where the moment equations of the ADE are analysed. This leads to the following two conditions:

(3.5a) \begin{align} \boldsymbol{\varphi }(t)\overline {\boldsymbol{\varXi }}&= \boldsymbol{0}, \\[-9pt] \nonumber \end{align}

(3.5b) \begin{align} \boldsymbol{\varphi }(t)\boldsymbol{\zeta }(t) + [\boldsymbol{\varphi }(t) \boldsymbol{\zeta }(t)]^\top &= 2\boldsymbol{D}. \\[9pt] \nonumber \end{align}

The first condition states that the contribution of the noise to the average drift is zero and it is fulfilled by $\overline {\boldsymbol{\varXi }}=\boldsymbol{0}$ . The second sets the contribution of the noise to the second moments. To do so, it is convenient to employ a streamline coordinate system, i.e. a locally orthogonal reference frame where one of the basis vectors is locally aligned with the velocity vector, as shown below. In what follows we will further simplify the description by considering spatially constant and isotropic dispersion, so that $\boldsymbol{D}=D\delta _{\textit{ij}}$ , where $D$ is a dispersion (or diffusion) coefficient.

3.2. Streamline coordinate frame

To determine the matrix $\boldsymbol{\varphi }(t)$ from (3.5b ) explicitly, we move to the streamline coordinate frame, where it is upper triangular. To see this, consider first the distance vector between two nearby streamlines in an orthogonal frame that follows one of them and has the $x$ -direction oriented along it (Winter Reference Winter1982). In this frame, we locally linearise the flow field following Dentz et al. (Reference Dentz, Lester, Le Borgne and De Barros2016). We then restore dispersion to this picture by looking at the evolution of the position $\boldsymbol{z}(t)$ of a diffusive particle in this frame, whose evolution is governed by the following FPE (see Dentz et al. Reference Dentz, de Barros, Le Borgne and Lester2018 and Appendix D):

(3.6) \begin{align} & \frac {\partial f_{\boldsymbol{Z}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\epsilon } \boldsymbol{z} f_{\boldsymbol{Z}} = D {\nabla} ^2 f_{\boldsymbol{Z}}, \qquad f_{\boldsymbol{Z}}(\boldsymbol{z};t_0) = f^0_{\boldsymbol{Z}}(\boldsymbol{z}). \end{align}

Here $\boldsymbol{\epsilon }(t;\boldsymbol{{x}}_0,t_0)$ is the velocity gradient tensor in the streamline coordinate system,

(3.7) \begin{align} \boldsymbol{\epsilon }(t;\boldsymbol{{x}}_0,t_0) = \begin{bmatrix} \alpha (t;\boldsymbol{{x}}_0,t_0) & \gamma (t;\boldsymbol{{x}}_0,t_0) \\ 0 & -\alpha (t;\boldsymbol{{x}}_0,t_0) \end{bmatrix} \!, \end{align}

whose components are the shear rate $\gamma (t;\boldsymbol{{x}}_0,t_0)$ and the elongation rate $\alpha (t;\boldsymbol{{x}}_0,t_0)$ along the streamline; see (D8). This description can be extended to three-dimensional flows by the use of the so-called Protean coordinate frame (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018).

Using (3.5b ) and (3.2c ) and exploiting the upper-triangular form of the velocity gradient tensor in this frame, we obtain the components of $\boldsymbol{\varphi }(t;\boldsymbol{x}_0,t_0)$ as

(3.8) \begin{equation} \begin{aligned} &\varphi _{11}(t) = \frac {D}{\zeta _{11}(t)} \left (\frac {\zeta _{21}(t)^2}{\zeta _{22}(t)^2}+1\right ), &&\varphi _{12}(t) = -\frac {D\zeta _{21}(t)}{\zeta _{22}(t)^2},\\ &\varphi _{21}(t) = 0, &&\varphi _{22}(t) = \frac {D}{\zeta _{22}(t)}, \end{aligned} \end{equation}

where the components of the correlation matrix $\boldsymbol{\zeta }(t;\boldsymbol{x}_0,t_0)$ obey the Bernoulli ((3.9a ) and (3.9d )) and linear ((3.9b ) and (3.9c )) differential equations

(3.9a) \begin{align} &\frac {\text{d}\zeta _{11}(t)}{\text{d}t} -\alpha (t)\zeta _{11}(t) = \frac {D}{\zeta _{11}(t)} \left (\frac {\zeta _{21}(t)^2}{\zeta _{22}(t)^2}+1\right )\! , &\zeta _{11}(t_0) &=0, \\[-9pt] \nonumber\end{align}

(3.9b) \begin{align} &\frac {\text{d}\zeta _{12}(t)}{\text{d}t} + \alpha (t)\zeta _{12}(t) = 0, &\zeta _{12}(t_0) &= 0, \\[-9pt] \nonumber \end{align}

(3.9c) \begin{align} &\frac {\text{d}\zeta _{21}(t)}{\text{d}t} + \left ( \frac {D}{\zeta _{22}(t)^2} - \alpha (t) \right ) \zeta _{21}(t) = \gamma (t)\zeta _{22}(t), &\zeta _{21}(t_0) &=0, \\[-9pt] \nonumber \end{align}

(3.9d) \begin{align} &\frac {\text{d}\zeta _{22}(t)}{\text{d}t} +\alpha (t)\zeta _{22}(t) = \frac {D}{\zeta _{22}(t)}, &\zeta _{22}(t_0) &=0, \\[9pt] \nonumber \end{align}

with solutions

(3.10) \begin{equation} \begin{aligned} &\zeta _{11}(t) = \left [ 2D p(t)\!\int _{t_0}^{t} \frac {\text{d}t^\prime }{p(t^\prime )} \left ( \frac {\zeta _{21}(t^\prime )^2}{\zeta _{22}(t^\prime )^2 } + 1\right ) \right ]^{1/2} , &\zeta _{12}(t) &= 0,\\ &\zeta _{21}(t) = \frac {1}{q(t)}\int _{t_0}^{t} \text{d}t^\prime q(t^\prime ) \gamma (t^\prime ) \zeta _{22}(t^\prime ), &\zeta _{22}(t) &= \left ( \frac {2D}{ p(t)} \int _{t_0}^t \text{d}t^\prime p(t^\prime )\right )^{1/2}, \end{aligned} \end{equation}

where

(3.11) \begin{align} p(t) = \textrm{e}^{2 \int _{t_0}^{t} \text{d}t \alpha (t)}, \quad q(t) = \exp {\left [ \int _{t_0}^{t} \text{d}t \left ( \frac {D}{\zeta _{22}(t)^2} -\alpha (t) \right ) \right ]}. \end{align}

With the source term $\boldsymbol{\varphi }(t)$ defined by (3.8), the hyperbolic description is fully characterised. Using (2.6) and (3.1), the Liouville master equation in conservative form is

(3.12) \begin{align} \frac {\partial f_{\boldsymbol{Z}\boldsymbol{\varXi }}}{\partial t} &+ \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \boldsymbol{\epsilon }\boldsymbol{z} + \boldsymbol{\varphi }\boldsymbol{\xi }\right ) f_{\boldsymbol{Z}\boldsymbol{\varXi }} = 0 , \qquad f_{\boldsymbol{Z}\boldsymbol{\varXi }}(\boldsymbol{z},\boldsymbol{\xi };t_0) = f^0_{\boldsymbol{Z}}(\boldsymbol{z})f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }), \end{align}

and the corresponding system of random differential equations(2.5) reads

(3.13a) \begin{align} \frac {\text{d}Z_1(t)}{\text{d}t} &= \alpha (t)Z_1 + \gamma (t) Z_2 + \varXi _1 \varphi _{11}(t) + \varXi _2 \varphi _{12}(t), &Z_1(t_0) &= Z_1^0, \\[-9pt] \nonumber\end{align}

(3.13b) \begin{align} \frac {\text{d}Z_2(t)}{\text{d}t} &= -\alpha (t)Z_2 + \varXi _2 \varphi _{22}(t), &Z_2(t_0) &= Z_2^0. \\[9pt] \nonumber \end{align}

For a given choice of distribution for the initial condition $\boldsymbol{Z}^0$ and random variables $\boldsymbol{\varXi }$ , the system (3.13) can be solved numerically using Monte Carlo method. Since for any given value of these random variables, (3.13) becomes an ODE system, we can say more about this system analytically using the method of characteristics as before. The characteristic lines $\hat {\boldsymbol{z}}(t;\boldsymbol{\xi },\boldsymbol{x}_0,t_0)$ of (3.12) solve

(3.14a) \begin{align} \frac {\text{d}\hat {\boldsymbol{z}}(t)}{\text{d}t} &= \boldsymbol{\epsilon }(t)\hat {\boldsymbol{z}}(t)+ \boldsymbol{\varphi }(t)\boldsymbol{\xi }, &\hat {\boldsymbol{z}}(t_0)& = \boldsymbol{z}_0, \\[-9pt] \nonumber \end{align}

(3.14b) \begin{align} \frac {\text{d}\hat {f}_{\boldsymbol{Z}\boldsymbol{\varXi }}(t)}{\text{d}t} &= 0, &\hat {f}_{\boldsymbol{Z}\boldsymbol{\varXi }}(t_0) &= f^0_{\boldsymbol{Z}}( \boldsymbol{z}_0 )f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }). \\[9pt] \nonumber \end{align}

With the previous results, the analytical solution for these characteristic lines can be expressed as the sum of a homogeneous and a particular solution,

(3.15) \begin{align} \hat {\boldsymbol{z}}(t) = \boldsymbol{B}(t) \boldsymbol{z}_0 + \boldsymbol{\zeta }(t)^\top \boldsymbol{\xi }, \end{align}

such that the homogeneous part relates to the deformations without considering diffusion, i.e. to the kinematics of the flow according to (D6). An analytical solution for this part is given in Dentz et al. (Reference Dentz, Lester, Le Borgne and De Barros2016, Reference Dentz, de Barros, Le Borgne and Lester2018), from which

(3.16) \begin{equation} \begin{aligned} B_{11}(t) &= \frac {v(t)}{v(t_0)} , &B_{12}(t) &= v(t_0)v(t) \int _{t_0}^t\text{d}t^\prime \frac {\gamma (t^\prime )}{v(t^\prime )^{2}} , \\ B_{21}(t) &= 0 , &B_{22}(t) &= \frac {v(t_0)}{v(t)}. \end{aligned} \end{equation}

In terms of the moment description (3.2), it corresponds to the average path of the blob of solute. The particular solution, on the other hand, encodes the effect of diffusion in terms of the random coefficients. It produces deviations from the average path, such that characteristic lines with large values of the coefficients $\boldsymbol{\xi }$ deviate more from the average solution encoded in the homogeneous part.

Then, the complete solution to the characteristic system (3.14), marginalised along the coefficients to obtain the concentration field in the streamline coordinate system, is

(3.17) \begin{align} f_{\boldsymbol{Z}}(\boldsymbol{z};t) &= \int _{\mathbb{R}^d}\text{d}\boldsymbol{\xi } \, f^0_{\boldsymbol{Z}} \big[\hat {\boldsymbol{z}}_0(\boldsymbol{\xi },\boldsymbol{z},t) \big]f_{\boldsymbol{\varXi }}(\boldsymbol{\xi }), \end{align}

where from (3.15) the backward map $\hat {\boldsymbol{z}}_0(\boldsymbol{\xi },\boldsymbol{z},t)$ is simply

(3.18) \begin{align} \hat {\boldsymbol{z}}_0(\boldsymbol{\xi },\boldsymbol{z},t) = \boldsymbol{B}(t)^{-1} \big ( {\boldsymbol{z}} - \boldsymbol{\zeta }(t)^\top \boldsymbol{\xi } \big ). \end{align}

In the particular case of a deterministic (point) initial condition, where $f_{\boldsymbol{Z}}^0(\boldsymbol{z})=\delta (\boldsymbol{z})$ , the properties of the Dirac delta can be used to simplify the integral in (3.17), leading to the Gaussian distribution

(3.19) \begin{align} f_{\boldsymbol{Z}}(\boldsymbol{z};t) = \frac {1}{2\pi \det [\boldsymbol{\zeta }(t)]}\exp \left (-\frac {1}{2}\boldsymbol{z}^\top \big [ \boldsymbol{\zeta }(t)^\top \boldsymbol{\zeta }(t) \big ]^{-1} \boldsymbol{z} \right )\!, \end{align}

with centred second moments $\boldsymbol{\kappa }(t) = \boldsymbol{\zeta }(t)^\top \boldsymbol{\zeta }(t)$ and $\boldsymbol{\kappa }(0)=\boldsymbol{\kappa }^0=\boldsymbol{0}$ . Note that the choice of independent, normally distributed random coefficients ensures this result is exact for a linear flow with a constant velocity gradient tensor (Van Kampen Reference Van Kampen1992; Risken Reference Risken1996; Bolster, Dentz & Le Borgne Reference Bolster, Dentz and Le Borgne2011; Dentz et al. Reference Dentz, de Barros, Le Borgne and Lester2018). More generally, for other forms of $f_{\boldsymbol{Z}}^0$ of interest, (3.17) can be easily computed numerically.

5.1. Constant strain-rate tensor

We first consider two canonical cases in which the velocity gradient tensor $\boldsymbol{\epsilon }$ is constant in space.

5.1.1. Simple shear flow

In the case of a pure linear shear flow, the velocity field is given by

(5.1) \begin{align} \boldsymbol{u}(\boldsymbol{x})=\gamma x_2 \boldsymbol{e}_1. \end{align}

For this case, (3.8) and (3.10) give

(5.2) \begin{align} \begin{split} \varphi _{11}(t) &= \frac { D \left (1+\dfrac {1}{4}\gamma ^2 t^2\right )}{\sqrt { 2 D t +\dfrac {1}{6}D\gamma ^2 t^3}} , \quad \varphi _{12}(t) = -\frac {\gamma \sqrt {D t}}{2\sqrt {2}}, \\ \varphi _{21}(t) &= 0 , \qquad \qquad \qquad \qquad \varphi _{22}(t) = \sqrt {\frac {D}{2t}}, \end{split} \end{align}

and

(5.3) \begin{equation} \begin{aligned} \zeta _{11}(t) &= \sqrt {2 D t + \frac {1}{6}D\gamma ^2t^3} , &\zeta _{12}(t) &= 0 , \\ \zeta _{21}(t) &= \frac {1}{2}{\gamma }t\sqrt {2 D t}, &\zeta _{22}(t) &= \sqrt {2Dt}. \end{aligned} \end{equation}

These results can be used to obtain the characteristics from (3.15):

(5.4) \begin{align} \hat {\boldsymbol{z}}(t;\boldsymbol{\xi },\boldsymbol{z}_0) &= \begin{bmatrix} z_{0,1} + \gamma t z_{0,2} + \xi _1\sqrt {2 D t+\dfrac {1}{6}D\gamma ^2 t^3} + \dfrac {1}{2} \xi _2 \gamma t \sqrt {2 D t} \\[9pt] z_{0,2} + \xi _2 \sqrt {2 D t} \end{bmatrix} \! . \end{align}

Using the moment description (3.2), the evolution of the covariance matrix $\boldsymbol{\kappa }(t)$ can also be obtained analytically, i.e.

(5.5a) \begin{align} \kappa _{11}(t) &= \kappa _{11}^0 + 2\gamma t\kappa _{12}^0 + \gamma ^2 t^2\kappa _{22}^0 + 2 D t + \frac {2}{3}D \gamma ^2 t^3 , \\[-9pt] \nonumber \end{align}

(5.5b) \begin{align} \kappa _{12}(t) &= \kappa _{12}^0 + \gamma t \kappa _{22}^0 + D\gamma t^2 , \\[-9pt] \nonumber \end{align}

(5.5c) \begin{align} \kappa _{22}(t) &= \kappa _{22}^0+ 2 D t , \\[9pt] \nonumber \end{align}

in agreement with the results in Bolster et al. (Reference Bolster, Dentz and Le Borgne2011) and Dentz et al. (Reference Dentz, de Barros, Le Borgne and Lester2018). We note that the asymptotic decay of the maximum concentration $c_m(t) \propto 1/\sqrt {\det {[\boldsymbol{\kappa }}(t)]}$ at long times in a simple shear flow is $c_m(t) \propto t^{-2}$ , which is not captured by the lamellae model because diffusion along the lamella orientation is neglected (Souzy et al. Reference Souzy, Zaier, Lhuissier, Le Borgne and Metzger2018).

Relations (5.2) and (5.3) can be used to describe the evolution of any initial concentration profile. For example, the case of a point source is given analytically by (3.19), which leads to the classical result (Risken Reference Risken1996)

(5.6) \begin{align} f_{\boldsymbol{Z}}(\boldsymbol{z};t)=\frac {\exp \left ( -\dfrac {1}{2}\boldsymbol{z}^\top \boldsymbol{\kappa }^{-1}(t)\boldsymbol{z} \right )}{2 \pi \sqrt {\det [\boldsymbol{\kappa }(t)]}}, \end{align}

with

(5.7) \begin{align} {\kappa }_{11}(t) = 2 D t + \frac {2}{3}D \gamma ^2 t^3 , \qquad {\kappa }_{12}(t) = D\gamma t^2 , \qquad {\kappa }_{22}(t) =2 D t . \end{align}

If we consider an initial condition composed by $N_I$ discrete point injection locations $\tilde {\boldsymbol{{x}}}_k$ along the horizontal axis,

(5.8) \begin{align} f_{\boldsymbol{X}}^0(\boldsymbol{x}) = \frac {1}{N_I}\sum _{k=1}^{N_I} \delta (\boldsymbol{x}-\tilde {\boldsymbol{x}}_k), \end{align}

the PDF of solute particle locations is given analytically by

(5.9) \begin{align} f_{\boldsymbol{X}}(\boldsymbol{x};t) = \frac {1}{N_I} \sum _{k=1}^{N_I}\frac {\exp \left [ -\dfrac {1}{2}(\boldsymbol{x}-\tilde {\boldsymbol{x}}_k)^\top \boldsymbol{\kappa }(t)^{-1}(\boldsymbol{x}-\tilde {\boldsymbol{x}}_k) \right ]}{2 \pi \sqrt {\det [\boldsymbol{\kappa }(t)]}}. \end{align}

An example with $P=2$ is shown in figure 2(a–d), evaluated at discrete support particle locations $\hat {\boldsymbol{x}}_k = \tilde {\boldsymbol{{x}}}_k + \hat {\boldsymbol{z}}$ , using (5.4) for $\hat {\boldsymbol{z}}$ with $\boldsymbol{z}_0=\boldsymbol{0}$ , see also (4.1), where by virtue of (5.1) one has $\boldsymbol{A}(t)=\boldsymbol{\mathcal{I}}$ , the identity matrix. The blob evolution of arbitrary initial conditions can also be computed with (4.4). In figure 2(e–h) we show the evolution of a multimodal distribution based on skew-normal and Gaussian distributions used as an initial condition, by integrating (4.4) numerically with the trapezoidal rule.

Figure 2. Evolution of the PDF of solute particle locations $f_{\boldsymbol{X}}$ in a shear flow from a deterministic initial condition composed of two injection points according to (5.8) with $(\tilde {x}_{1,1},\tilde {x}_{2,1})=(-2,0)$ and $(\tilde {x}_{1,2},\tilde {x}_{2,2})=(2,0)$ , with $D=0.2 \ \text{m}^2 \text{s}^{-1}$ at (a) $t=0 \ \text{s}$ , (b) $t=1.67 \ \text{s}$ , (c) $t=3.33 \ \text{s}$ and (d) $t=5 \ \text{s}$ ; and a multimodal initial condition with $D=2\times 10 ^{-3} \ \text{m}^2 \text{s}^{-1}$ at (e) $t=0 \ \text{s}$ , (f) $t=1.33 \ \text{s}$ , (g) $t=2.67 \ \text{s}$ and (h) $t=4 \ \text{s}$ . For both cases, $\gamma =1 \ \text{s}^{-1}$ .

5.1.2. Straining flow with shear

Consider now a velocity field given by

(5.10) \begin{align} \boldsymbol{u}(\boldsymbol{x}) = (\alpha x_1 + \gamma x_2) \boldsymbol{e}_1 - \alpha x_2 \boldsymbol{e}_2. \end{align}

This corresponds to the most general form of gradient tensor (3.7) needed to treat the incompressible, time-dependent case in the streamline frame. This example is of special importance because at the initial time $\boldsymbol{\zeta }(0)=\boldsymbol{0}$ in general, and the source term $\boldsymbol{\varphi }(t)$ given by (3.8) becomes unbounded, causing numerical issues in the computation of characteristic lines. This occurs because the characteristics deviate from the local streamline infinitely fast at arbitrarily small times, and it is consistent with the small-time limit of the FPE, which leads to infinitely fast initial velocities (Zaslavsky Reference Zaslavsky2002). This singularity can be avoided by assuming a constant strain-rate tensor during a small initial time step, over which the characteristics can be obtained analytically with the following results, before continuing the computation numerically.

In this case, (3.8) and (3.10) give

(5.11a) \begin{align} \varphi _{11}(t) &= \frac { \sqrt {D} \big ( \textrm{e}^{\alpha t} \big(4\alpha ^2+\gamma ^2 \big)\sinh ^2[\alpha t] + \alpha \gamma ^2 \textrm{e}^{\alpha t } \big(1-\textrm{e}^{2\alpha t}\big)t + \alpha ^2 \gamma ^2 \textrm{e}^{3\alpha t}t^2 \big )}{ 2 \sqrt {\alpha \textrm{e}^{3 \alpha t} \sinh ^3[\alpha t] \left ( -2 \alpha ^2 \left (\gamma ^2 t^2+2\right )+\left (4 \alpha ^2+\gamma ^2\right ) \cosh [2 \alpha t] + \gamma ^2\right )} } , \\[-9pt] \nonumber \end{align}

(5.11b) \begin{align} \varphi _{12}(t) &= -\frac { \sqrt {D} \gamma \textrm{e}^{\alpha t} \left (\textrm{e}^{2 \alpha t} (2 \alpha t-1)+1\right )}{2 \sqrt {\alpha } \left (\textrm{e}^{2 \alpha t}-1\right )^{3/2}} , \\[-9pt] \nonumber \end{align}

(5.11c) \begin{align} \varphi _{21}(t) &= 0 , \\[-9pt] \nonumber \end{align}

(5.11d) \begin{align} \varphi _{22}(t) &= \sqrt {\frac {\alpha D}{1-\textrm{e}^{-2 \alpha t}}} , \\[9pt] \nonumber \end{align}

and

(5.12a) \begin{align} \zeta _{11}(t) &= \frac {\sqrt {D \left [ \left (4 \alpha ^2+\gamma ^2\right ) \cosh (2 \alpha t) - 2 \alpha ^2 \left (\gamma ^2 t^2+2\right ) -\gamma ^2\right ]}}{ \alpha ^{3/2} \sqrt {2(1-\textrm{e}^{-2\alpha t})}} , \\[-9pt] \nonumber \end{align}

(5.12b) \begin{align} \qquad \zeta _{12}(t) &= 0 , \\[-9pt] \nonumber \end{align}

(5.12c) \begin{align} \zeta _{21}(t) &= \frac { \gamma \sqrt {D} \left ( 2 \alpha t -1+ \textrm{e}^{-2\alpha t}\right )}{2 \alpha ^{3/2}\sqrt {1 - \textrm{e}^{-2 \alpha t}}} , \\[-9pt] \nonumber \end{align}

(5.12d) \begin{align} \zeta _{22}(t) &= \sqrt {\frac {D}{\alpha }\left (1-\textrm{e}^{-2 \alpha t} \right )}. \\[9pt] \nonumber \end{align}

With these source terms and the moment description (3.2) for an arbitrary initial condition, we obtain the following result:

(5.13a) \begin{align} \begin{split} \kappa _{11}(t) &= \kappa _{11}^0\textrm{e}^{2\alpha t} + \kappa _{12}^0 \frac {\gamma }{\alpha }\big(\textrm{e}^{2\alpha t}-1 \big) +\kappa _{22}^0\frac {\gamma ^2}{\alpha ^2}\sinh ^2 (\alpha t) + \frac {D}{\alpha }\big(\textrm{e}^{2\alpha t} -1 \big) \\ &+ \frac {D \gamma ^2}{2\alpha ^3}\sinh (2\alpha t) - \frac {D \gamma ^2}{\alpha ^2}t , \end{split} \\[-9pt] \nonumber \end{align}

(5.13b) \begin{align} \kappa _{12}(t) &= \kappa _{12}^0 + \kappa _{22}^0 \frac {\gamma }{2\alpha }\big(1-\textrm{e}^{-2\alpha t}\big)+\frac {D \gamma }{2\alpha ^2}\big ( 2\alpha t +\textrm{e}^{-2\alpha t}-1 \big ) , \\[-9pt] \nonumber \end{align}

(5.13c) \begin{align} \kappa _{22}(t) &= \kappa _{22}^0\textrm{e}^{-2\alpha t} + \frac {D}{\alpha }\big (1-\textrm{e}^{-2\alpha t} \big ) . \\[9pt] \nonumber \end{align}

These results reduce to the previous example in the limit $\alpha \to 0$ . For $|\alpha | \gt 0$ , the long-term scaling of the maximum concentration is exponential, $c_m(t) \propto 1/\sqrt {\det {\kappa (t)}} \propto \exp (-|\alpha | t)$ . This result remains true for pure straining ( $\gamma =0$ ) and agrees with the lamella theory prediction (Villermaux Reference Villermaux2019).

5.2. Vortex flow

Here we consider the idealised vortex flow discussed in Meunier & Villermaux (Reference Meunier and Villermaux2003), given in polar coordinates by

(5.14) \begin{align} \boldsymbol{u}(r,\theta ) = \frac {\varGamma }{2\pi r}\boldsymbol{e}_{\theta }. \end{align}

The only non-zero component of the strain-rate tensor (3.7) in the streamline coordinate system is $\epsilon _{12} = \gamma (r) \equiv \gamma (\boldsymbol{{x}}_0)$ , where the radius of a streamline is also only dependent on the initial condition by virtue of (5.14). One thus has the relations

(5.15) \begin{align} \gamma (r) = \frac {\varGamma }{\pi r^2} \equiv \gamma (\boldsymbol{{x}}_0) = \frac {\varGamma }{\pi |\boldsymbol{{x}}_0|^2}, \end{align}

with

(5.16) \begin{align} r = | \boldsymbol{{x}}_0 | . \end{align}

Thus, along a given streamline, characterised by constant radius $r$ , the solute experiences constant shear only. The trajectories corresponding to initial locations $\boldsymbol{{x}}_0$ at time $t_0=0$ are given by

(5.17) \begin{align} \boldsymbol{{x}}(t;\boldsymbol{{x}}_0) = r \cos [ \theta (t;\boldsymbol{{x}}_0) ] \boldsymbol{e}_{1} + r \sin [ \theta (t;\boldsymbol{{x}}_0) ] \boldsymbol{e}_{2} , \end{align}

with

(5.18) \begin{align} \theta (t;\boldsymbol{{x}}_0) = \theta _0(\boldsymbol{{x}}_0) + \frac { \varGamma }{2 \pi r^2} t \ , \qquad \theta _0(\boldsymbol{{x}}_0) = \text{atan}(\mathrm{x}_{0,2} / \mathrm{x}_{0,1}). \end{align}

According to (D5) the entries of $\boldsymbol{A}(t;\boldsymbol{{x}}_0)$ are defined by

(5.19) \begin{align} \boldsymbol{v}(t;\boldsymbol{{x}}_0) = - \frac { \varGamma }{2 \pi r} \sin [\theta (t;\boldsymbol{{x}}_0)] \boldsymbol{e}_{1} + \frac { \varGamma }{2 \pi r} \cos [\theta (t;\boldsymbol{{x}}_0)] \boldsymbol{e}_{2}, \end{align}

leading to

(5.20) \begin{align} \boldsymbol{A}(t;\boldsymbol{{x}}_0) = \begin{bmatrix} -\sin \big[\theta _0(\boldsymbol{{x}}_0)+ \varGamma t /\big(2\pi r^2 \big)\big] & -\cos \big[\theta _0(\boldsymbol{{x}}_0)+\varGamma t /\big(2\pi r^2 \big)\big] \\[6pt] \cos \big[\theta _0(\boldsymbol{{x}}_0)+\varGamma t /\big(2\pi r^2 \big)\big] & -\sin \big[\theta _0(\boldsymbol{{x}}_0)+\varGamma t /\big(2\pi r^2 \big)\big] \end{bmatrix}. \end{align}

Since for each streamline we have constant shear, the entries of $\boldsymbol{\varphi }(t;\boldsymbol{{x}}_0)$ and $\boldsymbol{\zeta }(t;\boldsymbol{{x}}_0)$ are the same as those for shear flow, so that they are given by (5.2) and (5.3) with $\gamma$ given by (5.15). This yields

(5.21) \begin{align} \begin{split} \varphi _{11}(t;\boldsymbol{{x}}_0) &= \frac { \sqrt {3} D \left ( \varGamma ^2 t^2 +4\pi r^4 \right ) }{2\pi r^2 \sqrt {2D \varGamma ^2 t^3 +24D t \pi r^4}} , \quad \varphi _{12}(t;\boldsymbol{{x}}_0) = -\frac {\varGamma \sqrt {D t}}{2 \pi r^2 \sqrt {2}}, \\ \varphi _{21}(t;\boldsymbol{{x}}_0) &= 0 , \qquad \qquad \qquad \qquad \qquad \quad \varphi _{22}(t;\boldsymbol{{x}}_0) = \sqrt {\frac {D}{2t}}, \end{split} \end{align}

and

(5.22) \begin{align} \begin{split} \zeta _{11}(t;\boldsymbol{{x}}_0) &= \sqrt {2 D t + \frac {D \varGamma ^2 t^3}{6\pi r^4}} , \qquad \zeta _{12}(t;\boldsymbol{{x}}_0) = 0 , \\ \zeta _{21}(t;\boldsymbol{{x}}_0) &= \frac {\varGamma t \sqrt {2 D t}}{2 \pi r^2}, \qquad \qquad \quad \zeta _{22}(t;\boldsymbol{{x}}_0) = \sqrt {2Dt}; \end{split} \end{align}

recall (5.16).

The characteristic lines in the streamline coordinate frame are given according to (3.15) as

(5.23) \begin{align} \hat {\boldsymbol{z}}(t;\boldsymbol{\xi },\boldsymbol{z}_0) &= \begin{bmatrix} z_{0,1} +\dfrac {\varGamma z_{0,2}}{\pi r^2}t + \xi _{1}\sqrt {2 D t + \dfrac {D\varGamma ^2 t^3}{6 \pi r^4 }} + \xi _2\dfrac {\varGamma t \sqrt {2 D t}}{2 \pi r^2} \\[9pt] z_{0,2} +\xi _2 \sqrt {2 D t} \end{bmatrix} \! . \end{align}

The forward map (4.1) of support particles can be obtained explicitly by combining (5.17), (5.20) and (5.22) (not shown for brevity). The backward map $\hat {\boldsymbol{x}}_0$ for an arbitrary point can be similarly obtained from (4.6) given the matrix $\boldsymbol{B}(t;\boldsymbol{{x}}_0)$ , which according to (3.16) is

(5.24) \begin{align} \boldsymbol{B}(t;\boldsymbol{{x}}_0) = \begin{bmatrix} 1 & \dfrac {\varGamma }{\pi r^2} t \\[9pt] 0 & 1 \end{bmatrix}\!. \end{align}

For a given initial condition, using the analytical expressions (5.14)–(5.24), the normalised concentration field can be explicitly computed for a Lagrangian unstructured grid using (4.4) or an imposed grid using (4.7).

Numerical algorithms based on the concept of lamellae have been constructed to consider initial concentration profiles that are symmetric along the transversal direction of the lamellae and constant along the longitudinal direction (Meunier & Villermaux Reference Meunier and Villermaux2010; Martínez-Ruiz et al. Reference Martínez-Ruiz, Meunier, Favier, Duchemin and Villermaux2018; Sen et al. Reference Sen, Singh, Heyman, Le Borgne and Bandopadhyay2020; Meunier & Villermaux Reference Meunier and Villermaux2022; Pushenko et al. Reference Pushenko, Scollo, Meunier, Villermaux and Schumacher2025). Additionally, lamella models neglect diffusion along lamellae, a mechanism that can be important at long times or low Péclet numbers. These approximations are not needed in the present framework. To illustrate this, we introduce uniform, normal and skew-normal distributions:

(5.25a) \begin{align} f_{\mathcal{U}}(x;a,b) &= \frac {\mathcal{H}(x-a)\mathcal{H}(b-x)}{b-a} , \\[-9pt] \nonumber \end{align}

(5.25b) \begin{align} f_{\mathcal{N} }(x;\mu ,\sigma ) &= \frac {1}{\sqrt {2\pi }\sigma } \exp \left [ -\frac {(x-\mu )^2}{2 \sigma ^2} \right ]\! , \\[-9pt] \nonumber \end{align}

(5.25c) \begin{align} f_{\mathcal{S} }(x;\mu ,\sigma ,\alpha ) &= \omega ^{-1} f_{\mathcal{N} }((x-\eta )/\omega ;0,1) \big( 1 + \textrm {erf} \big[ \alpha (x-\eta )/(\omega \sqrt {2}) \big] \big) . \\[12pt] \nonumber \end{align}

Here $\mathcal{H}(\boldsymbol{\cdot })$ denotes the Heaviside function and

(5.26) \begin{align} \omega = \frac {\sigma }{\sqrt {1-2\delta ^2/\pi }}, \qquad \eta = \mu - \omega \delta \sqrt {2/\pi }, \qquad \delta &= \frac {\alpha }{\sqrt {1+\alpha ^2}}. \end{align}

Then, we combine these shapes along the two Cartesian directions to form three example initial conditions:

(5.27a) \begin{align} f_{\boldsymbol{X}}^{0 \mathcal{U}}(\boldsymbol{x}) &= f_{\mathcal{U}}(x_1;a_1,b_1)f_{\mathcal{U}}(x_2;a_2,b_2), \\[-9pt] \nonumber \end{align}

(5.27b) \begin{align} f_{\boldsymbol{X}}^{0 \mathcal{UN}}(\boldsymbol{x}) &= f_{\mathcal{U}}(x_1;a_1,b_1)f_{\mathcal{N}}(x_2;\mu _2,\sigma _2), \\[-9pt] \nonumber \end{align}

(5.27c) \begin{align} f_{\boldsymbol{X}}^{0 \mathcal{S}}(\boldsymbol{x}) &= f_{\mathcal{S}}(x_1;\mu _1,\sigma _1,4)f_{\mathcal{S}}(x_2;\mu _2,\sigma _2,4). \\[9pt] \nonumber \end{align}

For comparison purposes, we choose the first two moments of all initial conditions along both Cartesian directions to be the same by imposing $\mu _1=(a_1+b_1)/2$ , $\mu _2=(a_2+b_2)/2$ , $\sigma _{1} = (b_1-a_1)/\sqrt {12}$ and $\sigma _{2} = (b_2-a_2)/\sqrt {12}$ , and we define the width as $s_0=b_2-a_2$ and the length as $l_0=b_1-a_1$ based on the uniform initial condition (5.27a ). Note that the initial condition (5.27a ) corresponds to the case in Meunier & Villermaux (Reference Meunier and Villermaux2003) and (5.27b ) to Meunier & Villermaux (Reference Meunier and Villermaux2010) and Sen et al. (Reference Sen, Singh, Heyman, Le Borgne and Bandopadhyay2020). Concentration profiles for these three initial conditions are compared in figure 3. To highlight the fact that the hyperbolic formulation introduced here allows for choosing the spatial locations where one wishes to compute the concentration profile at a given time, we also depict the evolution of the skewed concentration profile (5.27c ) for different grid choices in figure 4. In particular, Lagrangian grids composed by the positions of support particles (as in figure 3) are shown in figure 4(a–c) at different resolutions; uniform grids are shown in panels (d–f) and three arbitrary unstructured grids are shown in panels (g–i).

Figure 3. Temporal evolution of solute plumes in a vortex flow defined by (a) a uniform initial profile (5.27a ), (b) a uniform-normal initial profile (5.27b ) and (c) a skewed initial profile (5.27c ) at $t=0 \ \text{s}$ , $t=0.35 \ \text{s}$ and $t=0.8 \ \text{s}$ . Based on the set-up of Meunier & Villermaux (Reference Meunier and Villermaux2003), we set $a_1=0.6 \ \text{cm}$ , $b_1=1.8 \ \text{cm}$ , $a_2=-0.11 \ \text{cm}$ , $b_2=0.11 \ \text{cm}$ , $\varGamma =14.2 \ \text{cm}^2\, \text{s}^{-1}$ and we take $D=10^{-3} \ \text{cm}^2 \text{s}^{-1}$ as an intermediate value from those used in Meunier & Villermaux (Reference Meunier and Villermaux2010), which corresponds to $ \textit{Pe}_{\tilde {\gamma }} \approx 6 \times 10^2$ . The Lagrangian grids correspond to support particle locations with a uniform discretisation of $\boldsymbol{{x}}_0$ and $\boldsymbol{\xi }$ .

Figure 4. Temporal evolution of the skewed concentration profile (5.27c ) for different grid types. Top row: concentration profiles at $t=0 \ \text{s}$ , $t=0.35 \ \text{s}$ and $t=0.8 \ \text{s}$ for an initially uniform Lagrangian grid composed of support particles with three levels of refinement: (a) coarse, (b) medium and (c) fine. Middle row: concentration profiles for three uniform grids with different levels of refinement: (d) coarse at $t=0 \ \text{s}$ , (e) medium at $t=0.35 \ \text{s}$ and (f) fine at $t=0.8 \ \text{s}$ . Bottom row: arbitrary unstructured grids are employed to depict the solution at (g) $t=2 \ \text{s}$ , (h) $t=5\ \text{s}$ and (i) $t=10 \ \text{s}$ . The rest of the parameters are given in the caption of figure 3.

To further quantify the influence of radial diffusion on mixing in a vortex, we consider the maximum concentration $c_m(r,t)$ profile as a function of radial position and time both considering diffusion in all directions and neglecting radial diffusion. As we will see, local linearisation of the flow affects the results at late times. We begin by discussing the linearised behaviour and consider its role below. We focus on initial conditions (5.27a ) and (5.27b ), for which the calculations can be carried out analytically. Because the Jacobian of the transformation (2.8a ) is unity, $J(t)=1$ , and the domain is unbounded, concentration gradients cannot grow over time. This holds for incompressible flows with a constant diffusion coefficient by virtue of (2.8b ). Thus, the maximum radial concentration profile for a symmetric initial condition with respect to the $x_1$ axis will travel along those streamlines that at the initial time start along the line $x_2=0$ .

To describe the maximum concentration profile at time $t$ , consider a position $\boldsymbol{x}_m(r,t)$ such that $|\boldsymbol{x}_m(r,t)|=r$ . Since the flow has no radial component, this position corresponds to a maximum if and only if the backward flow map associates it with the position $\hat {\boldsymbol{{x}}}_0[\boldsymbol{x}_m(r,t),t] = (r, 0)^\top$ at the initial time; see (4.5). We first obtain the backward map in the streamline coordinate system for a given point $\boldsymbol{x}$ that at time $t$ coincides with the streamline, i.e. taking $\boldsymbol{z}=\boldsymbol{0}$ in (3.18). We find that

(5.28) \begin{align} \hat {\boldsymbol{z}}_0[\boldsymbol{\xi },\boldsymbol{x}_m(r,t),t] = \begin{bmatrix} \xi _2\dfrac {\varGamma t \sqrt {2 D t}}{2 \pi r^2} - \xi _{1}\sqrt {2 D t + \dfrac {D \varGamma ^2 t^3}{6 \pi ^2 r^4 }} \\[9pt] - \xi _2 \sqrt {2 D t} \end{bmatrix} \!, \end{align}

where $r=|\boldsymbol{x}|=|\hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t)|$ . The full backward map in the Cartesian frame is (see Appendix E)

(5.29) \begin{align} \hat {\boldsymbol{x}}_0(\boldsymbol{\xi },\boldsymbol{x},t) = \hat {\boldsymbol{{x}}}_0 (\boldsymbol{x},t) + \boldsymbol{A} \big[0;\hat {\boldsymbol{{x}}}_0 (\boldsymbol{x},t) \big] \hat {\boldsymbol{z}}_0(\boldsymbol{\xi },\boldsymbol{x},t), \end{align}

which combined with (5.28) gives

(5.30) \begin{align} \hat {\boldsymbol{x}}_0[\boldsymbol{\xi },\boldsymbol{x}_m(r,t),t] = \left ( r + \xi _2 \sqrt {2 D t} \right ) \boldsymbol{e}_1 + \left ( \xi _2\frac {\varGamma t \sqrt {2 D t}}{2 \pi r^2} - \xi _{1}\sqrt {2 D t + \frac {D \varGamma ^2 t^3}{6 \pi ^2 r^4 }} \right ) \boldsymbol{e}_2. \end{align}

Consider now what happens if we neglect radial diffusion. The second component of (5.28) represents deviations in the initial condition along the radial direction, because it is locally perpendicular to the flow direction. If no radial diffusion is present, no radial deviations can occur when tracing backward. Specifically, this corresponds to neglecting the noise contribution along the second component of the characteristic lines in the streamline coordinate system. That is, taking $[\boldsymbol{\zeta }^\top (t) \boldsymbol{\xi }]_2=0$ (and $\boldsymbol{z}=\boldsymbol{0}$ as before) in (3.18) leads to

(5.31) \begin{align} \hat {\boldsymbol{z}}_0[\boldsymbol{\xi },\boldsymbol{x}_m(r,t),t] = \begin{bmatrix} \xi _2\dfrac {\varGamma t \sqrt {2 D t}}{2 \pi r^2} - \xi _{1}\sqrt {2 D t + \dfrac {D \varGamma ^2 t^3}{6 \pi ^2 r^4 }} \\[9pt] 0 \end{bmatrix} \!, \end{align}

from which

(5.32) \begin{align} \hat {\boldsymbol{x}}_0[\boldsymbol{\xi },\boldsymbol{x}_m(r,t),t] = r \boldsymbol{e}_1 + \left ( \xi _2\frac {\varGamma t \sqrt {2 D t}}{2 \pi r^2} - \xi _{1}\sqrt {2 D t + \frac {D \varGamma ^2 t^3}{6 \pi ^2 r^4 }} \right ) \boldsymbol{e}_2. \end{align}

Considering the uniform initial condition (5.27a ), relation (4.7) together with (5.30) yields the maximum concentration profile

(5.33) \begin{align} \frac {c_m(r,t)}{c_0} &= \frac {f_{\boldsymbol{X}}[\boldsymbol{x}_m(r,t);t]}{\text{max}_{\boldsymbol{x}}[f_{\boldsymbol{X}}(\boldsymbol{x};0)]} = \mathcal{B}\left (k_2,\frac {-a}{\sqrt {b^2+1}};\frac {-b}{\sqrt {b^2+1}}\right ) - \mathcal{B}\left (k_1,\frac {-a}{\sqrt {b^2+1}};\frac {-b}{\sqrt {b^2+1}}\right ) \nonumber \\[3pt]& \quad +\mathcal{B}\left (k_2,\frac {a}{\sqrt {b^2+1}};\frac {-b}{\sqrt {b^2+1}}\right ) - \mathcal{B}\left (k_1,\frac {a}{\sqrt {b^2+1}};\frac {-b}{\sqrt {b^2+1}}\right ), \end{align}

in terms of the bivariate normal cumulative distribution according to Owen (Reference Owen1980):

(5.34) \begin{align} \mathcal{B}(\mathrm{k}_1,\mathrm{k}_2;\rho ) = \frac {1}{2\pi \sqrt {1-\rho ^2}}\int _{-\infty }^{\mathrm{k}_2}\text{d}y\,\int _{-\infty }^{\mathrm{k}_1} \text{d}x\, \exp \left [-\frac {x^2-2\rho x y +y^2}{2(1-\rho ^2)} \right ] \! , \end{align}

with

(5.35) \begin{equation} \begin{aligned} a &= \frac {\pi s_0 }{2 \sqrt {D \left (\frac {\varGamma ^2 t^3}{6 r^4}+2 \pi ^2 t\right )}}, & b &= \frac {\sqrt {3} \varGamma t}{\sqrt {\varGamma ^2 t^2+12 \pi ^2 r^4}}, \\ k_1 &= \frac {r-a_1}{\sqrt {2Dt}}, & k_2 &= \frac {b_1-r}{\sqrt {2Dt}}, \end{aligned} \end{equation}

where $c_0$ is the initial maximum concentration anywhere in the domain. If we neglect radial diffusion, combining (4.7) and (5.32) gives

(5.36) \begin{align} \frac {c_m(r,t)}{c_0} = \textrm {erf}\left [ \frac {\sqrt {3}\pi s_0 r^2 }{4\sqrt {D t \left ( 3\pi ^2 r^4 + \varGamma ^2 t^2 \right )}} \right ] \mathcal{H}[r-a_1]\mathcal{H}[b_1-r]. \end{align}

This coincides with the result in Meunier & Villermaux (Reference Meunier and Villermaux2003) obtained under this approximation. The difference between (5.33) and (5.36) can be clearly seen in figure 5(c). To understand when significant differences occur, consider the strip Péclet number based on the second moment $\sigma _1^2$ of the initial condition along the radial direction, $ \textit{Pe}_{\tilde {\gamma }} = \tilde {\gamma } \sigma _1^2/D$ , where $\tilde {\gamma }$ is the average shear over the initial condition. We have $\tilde {\gamma } \approx \varGamma (a_1^{-1}-b_1^{-1})/\pi$ provided that $s_0 \ll l_0$ , obtained by averaging (5.15) over $a_1\leqslant r \leqslant b_1$ . The strip Péclet number for the results in figure 5(a–c) has values $ \textit{Pe}_{\tilde {\gamma }} = 10^5$ , $6\times 10^3$ and $60$ . Diffusion along the radial direction becomes important for times of the order of the time scale $\tau _D = \sigma _1^2/D$ (Dentz et al. Reference Dentz, de Barros, Le Borgne and Lester2018). For the final time in figure 5, $t=20 \, \text{s}$ , we find that $t/t_D = 0.42$ (a), $0.83$ (b) and $1.67$ (c).

Proceeding the same way for the uniform-normal initial condition (5.27b ), the maximum concentration profile considering diffusion in all directions is

(5.37a) \begin{align} \frac {c_m(r,t)}{c_0} &= \frac {\pi s_0 r^2}{2\sqrt {8D \varGamma ^2 t^3 + \pi ^2 r^4 \big( s_0^2 +24 D t \big)}} \left ( \textrm {erf}\left [\frac {b_1-r}{\mathrm{q}(r,t)}\right ] - \textrm {erf}\left [ \frac {a_1-r}{\mathrm{q}(r,t)} \right ] \right ) \! , \\[-9pt] \nonumber \end{align}

(5.37b) \begin{align} \mathrm{q}(r,t) &= 2\sqrt { Dt - \frac {6 D \varGamma ^2 t^4}{ 8 \varGamma ^2 t^3 + \pi ^2 r^4 \big(s_0^2/D + 24t \big)}}, \\[9pt] \nonumber \end{align}

and neglecting radial diffusion, we find that

(5.38) \begin{align} \frac {c_m(r,t)}{c_0} = \frac {\pi s_0 r^2}{\sqrt {8D \varGamma ^2 t^3 + \pi ^2 r^4 \big( s_0^2 +24 D t \big)}}\mathcal{H}[r-a_1]\mathcal{H}[b_1-r]. \end{align}

These results are depicted in figure 5(d–f) for different values of the diffusion coefficient at three different times. The values of the strip Péclet number and times are the same as for figure 5(a–c) and, for the same reason as before, the difference between these solutions is mainly visible in figure 5(f). The solutions for the uniform initial condition considered above are also shown, highlighting the importance of the initial condition even at times $t\gt \tau _D$ .

Figure 5. Maximum concentration along the radial direction for the uniform initial condition (5.27a ) in (a–c) and the uniform-normal initial condition (5.27b ) in (d–f) at times $t=5 \ \textrm {s}$ (circles), $t=10 \ \textrm {s}$ (squares), $t=20 \ \textrm {s}$ (triangles) with $D=5\times 10^{-6}\ \textrm {cm}^2\, \textrm {s}^{-1}$ ( $ \textit{Pe}_{\tilde {\gamma }} \approx 10^5$ ) in (a) and (d); $D=10^{-4}\ \textrm {cm}^2\,\textrm {s}^{-1}$ ( $ \textit{Pe}_{\tilde {\gamma }} \approx 6\times 10^3$ ) in (b) and (e); and $D=10^{-2}\ \textrm {cm}^2\,\textrm {s}^{-1}$ ( $ \textit{Pe}_{\tilde {\gamma }} \approx 60$ ) in (c) and (f). The dash-dotted red lines in (a–c) are (5.36) and in (d–f) they are (5.38); whereas the dashed green lines in (a–c) are (5.33) and in (d–f) they are (5.37a ). The continuous black lines in (d–f) are (5.36). The radius is normalised by the initial vortex radius $a_0=0.3 \ \textrm {cm}$ (Meunier & Villermaux Reference Meunier and Villermaux2003). The remaining parameters are as in figure 3. The green curves consider radial diffusion, whereas the red and black curves do not.

The shape of the initial condition does not influence the asymptotic scaling of the decay of the maximum concentration, although the value of the maximum concentration at a given time depends on the initial condition. Neglecting radial diffusion, on the other hand, does affect the asymptotic scaling. To quantify this, denote by $r_m(t)$ the radius where the maximum concentration occurs for a given time, for which an implicit relation may be found using (5.33) and (5.37a ) for the initial conditions (5.27a ) and (5.27b ), respectively. For sufficiently long times, the maximum concentration anywhere in the domain evolves as $c_m(t) = c_m[r_m(t),t] \propto t^{-2}$ . If radial diffusion is neglected, the maximum concentration occurs at the maximum radius, i.e. one simply has $r_m(t)=b_1$ for all times from (5.36) and (5.38). Considering the uniform-normal initial condition (5.27b ), and neglecting the influence of radial diffusion, the maximum concentration profile at long times according to (5.38) evolves as $c_m(t) \propto t^{-3/2}$ instead of $t^{-2}$ . The same is true for the uniform initial condition (5.27a ) according to (5.36), making use of $\textrm {erf}(x)\approx 2x/\sqrt {\pi }$ for $x\ll 1$ .

The nonlinearity of the flow (5.14) restricts the validity of the scalings of the maximum concentration derived above to times shorter than the diffusion time, regardless of the Péclet number and the consideration of radial diffusion. After the combination of stretching and diffusion has homogenised the plume in the radial support $a_1\lesssim r \lesssim b_1$ , the solute forms an axisymmetric annular-shaped profile. This is accompanied by an abrupt change in the scaling of the maximum concentration. Eventually, a disk-shaped plume is formed, leading to a pure diffusion regime along the radial direction, characterised by a scaling $c_m(t) \propto t^{-1}$ .

Figure 6(a–e) shows the time evolution of the solute plume computed with the relinearisation procedure proposed in § 4. The behaviour of the maximum concentration with and without relinearisation is shown in figure 6(f) and compared with the solution obtained from a standard finite-difference PDE solver, using second-order central discretisation in space and a total variation diminishing, third-order Runge–Kutta method in time (Gottlieb & Shu Reference Gottlieb and Shu1998). The late-time dynamics of the concentration field is governed by the nonlinearities in the flow, which are only correctly captured if relinearisation is applied. It should be noted that (5.14) is a simplification of the Lamb–Oseen vortex model that is accurate far from the vortex core (Meunier & Villermaux Reference Meunier and Villermaux2003). We restrict the time frame in figure 6(f) such that the maximum of the solute plume is far from the centre, so that the artificial singularity at $r=0$ does not affect the results significantly. Correspondingly, we limit the velocities in the finite-difference solver to the value at $r=0.05 \ \text{cm}$ . For this reason, the $c_m(t) \propto t^{-1}$ scaling expected for a radially diffusing, disk-shaped plume cannot be observed cleanly under this model.

Figure 6. Time evolution of the uniform initial condition (5.25a ) for $D=10^{-2} \ \text{cm}^2\,\text{s}^{-1}$ , which corresponds to $ \textit{Pe}_{\tilde {\gamma }} \approx 60$ at (a) $t= 0 \ \text{s}$ , (b) $t= 1 \ \text{s}$ , (c) $t= 4 \ \text{s}$ , (d) $t= 6 \ \text{s}$ and (e) $t=10 \ \text{s}$ computed with the relinearisation procedure with a linear interpolant with nodes distributed uniformly in polar coordinates (with $301\times 601$ points in the radial and azimuthal directions, respectively) for a radial domain of $r \in [0, 3] \ \text{cm}$ and the full circumference. Relinearisation is applied at intervals of $\Delta t = 0.02 \ \text{s}$ . Panel (f) shows the maximum concentration assuming linearity of the drift and neglecting radial diffusion using (5.36) (dotted red line), assuming linearity of the drift but considering radial diffusion with (5.33) (dashed green line), using the relinearisation procedure while accounting for radial diffusion (continuous black line with squared symbols) and using the PDE solver in a $[-3,3]\times [-3,3] \ \text{cm}^2$ domain discretised in $401\times 401$ points with a time step of $\Delta t =6\times 10^{-5} \ \text{s}$ (dot-dashed blue line with triangular symbols).

5.3. Bidimensional heterogeneous Darcy flow

Here we consider the evolution of the concentration profile of a solute in a heterogeneous porous medium. In particular, we compare the hyperbolic formulation with the particle tracking simulations based on the Langevin system (2.3) with the concentration profile reconstructed according to a discretised version of (2.4) using square bins.

We consider a Darcy (continuum) scale, steady-state velocity field $\boldsymbol{u}(\boldsymbol{x})$ described by the Darcy equation (Bear Reference Bear1972),

(5.39) \begin{align} \boldsymbol{u}(\boldsymbol{x}) = -K(\boldsymbol{x})\boldsymbol{\nabla }\mathrm{h} (\boldsymbol{x}), \end{align}

where $K(\boldsymbol{x})$ is the (isotropic) hydraulic conductivity and $\mathrm{h}(\boldsymbol{x})$ the hydraulic head. We take the flow field to be incompressible, such that the continuity equation reads $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$ , leading to an explicit equation for the hydraulic head field,

(5.40) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }K(\boldsymbol{x})\boldsymbol{\nabla }\mathrm{h} (\boldsymbol{x}) = 0. \end{align}

We take $K(\boldsymbol{x})$ as a single realisation of a random field such that the log-conductivity field, $ f (\boldsymbol{x}) = \ln [ K ( \boldsymbol{x} )/K_0 ]$ , with $K_0$ an arbitrary reference permeability, has a log-Gaussian point distribution, ensemble mean $\mu _f$ , ensemble variance $\sigma _{f}^2$ and Gaussian spatial correlation structure with isotropic correlation length $l_c$ . Since the correlation length $l_c$ provides a characteristic spatial length scale, it is convenient to introduce a Péclet number $ \textit{Pe}=\overline {u}l_c/D$ , where $\overline {u}$ is the magnitude of the average of the Darcy velocity field (5.39) and $D$ is the dispersion coefficient. We have $ \textit{Pe}=\tau _D/\tau _A$ , where $\tau _D=\ell _c^2/D$ and $\tau _A=\ell _c/\overline {u}$ are the diffusion and advection time scales associated with this length scale. We normalise lengths by $\ell _c$ , times by $\tau _A$ and velocities by $\overline {u}$ , so that, due to linearity of the flow, the numerical values of these quantities do not play a role and we work with non-dimensional magnitudes $t/\tau _A$ , $\boldsymbol{x}/l_c$ and $\boldsymbol{u}/\overline {u}$ . Then, the flow field is obtained for a permeability field realisation with $\mu _f=1$ and $\sigma _{f}^2=1/4$ , and an imposed mean Darcy flow velocity aligned with the horizontal direction. To generate the permeability field and compute the flow field, we make use of the freely available solver developed in Zhou & Hansen (Reference Zhou and Hansen2022, Reference Zhou and Hansen2024). The permeability field is computed using a spectral approach and the flow is solved in a uniform grid using a finite volume method in terms of interfacial Darcy fluxes by combining (5.39) and (5.40) in integral form. We consider a bidimensional domain of $25\times 100$ correlation lengths along the horizontal and vertical directions, respectively, with periodic boundary conditions and discretised in $20$ cells per correlation length along each direction. A third-order interpolant is then used to compute the resulting flow field at arbitrary locations. For transport, the Péclet number is set by varying the value of the diffusion coefficient $D$ . We consider times up to $t_{\textit{end}} = 60 \tau _A$ and Péclet numbers $ \textit{Pe}=10^3$ and $10^4$ . These times are small compared with $\tau _D$ , respectively $t_{\textit{end}}/\tau _D = 6\times 10^{-2}$ and $6\times 10^{-3}$ , and we do not apply the relinearisation procedure proposed in § 4 here.

Because they are commonly employed and play an important role in theoretical and practical applications, as well as to illustrate certain advantages of the formulation introduced here, we focus on singular initial conditions, characterised by the presence of infinite spatial gradients. In particular, we consider the point and line initial injections

(5.41a) \begin{align} f_{\boldsymbol{X}}^p(\boldsymbol{x}; \boldsymbol{a}) &= \delta (\boldsymbol{x} - \boldsymbol{a}), \\[-9pt] \nonumber\end{align}

(5.41b) \begin{align} f_{\boldsymbol{X}}^l(\boldsymbol{x};\boldsymbol{b}) &= \frac {\mathcal{H}(x_2-b_1)\mathcal{H}(b_2-x_2)}{b_2-b_1} \delta (x_1-b_3), \\[9pt] \nonumber \end{align}

at the location $\boldsymbol{x} = \boldsymbol{a}$ and along the line $b_1 \leqslant x_{2} \leqslant b_2$ and ${x}_{1} = b_3$ . respectively. In the numerical examples below, we take $\boldsymbol{a}/l_c = (2,12.5)^\top$ for (5.41a ) and $\boldsymbol{b}/l_c = (8.5,16.5,2)^\top$ for (5.41b ). Note that the strip Péclet number along degenerate (zero-width) directions of the initial condition has zero width. These types of initial conditions present numerical difficulties for Eulerian treatments of the FPE (2.2a ), because regularisation of the Dirac delta and Heaviside functions is required (Suarez, Jacobs & Don Reference Suarez, Jacobs and Don2014; Suarez & Jacobs Reference Suarez and Jacobs2017; Wissink et al. Reference Wissink, Jacobs, Ryan, Don and van der Weide2018; Domínguez-Vázquez et al. Reference Domínguez-Vázquez, Jacobs and Tartakovsky2021). Particle tracking simulations, on the other hand, are easily applicable for these initial conditions. However, the reconstruction of the PDF $f_{\boldsymbol{X}}$ is sensitive to the coexistence of spatial scales with large disparity on the solute plume, which complicates the search for an optimal bandwidth in case of the use of kernel density estimators, or bin size in the case of histograms. As a consequence, practical implementations of the reconstruction (2.4) introduce errors that are hard to estimate without the knowledge of the true (unknown) concentration profile. It is not the intention of this paper to use state-of-the-art kernel density estimators for reconstructing the particle PDF $f_{\boldsymbol{X}}$ (Fernàndez-Garcia & Sánchez-Vila Reference Fernàndez-Garcia and Sánchez-Vila2011; Pedretti & Fernàndez-Garcia Reference Pedretti and Fernàndez-Garcia2013; Sole-Mari et al. Reference Sole-Mari, Bolster, Fernàndez-Garcia and Sanchez-Vila2019; Benson et al. Reference Benson, Bolster, Pankavich and Schmidt2021). Instead, we use histograms to reconstruct $f_{\boldsymbol{X}}$ from Monte Carlo computations of (2.3) using the Euler-Maruyama method, in order to compare and validate the proposed formulation. This way, we eliminate additional intricacies related to PDF reconstruction and focus on the illustration of some of the advantages of the hyperbolic formulation, such as its ability to handle singular profiles, targeting locations of interest and reducing the computational cost. All particle tracking results in this section employ $N=10^6$ particles.

Now we proceed to describe the evolution of the PDF of solute particle locations with relations (4.6) and (4.7). We look at an arbitrary point $\boldsymbol{x}$ at time $t$ , such that the backward flow map is $ \hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t)=\boldsymbol{{x}}(0;\boldsymbol{x},t)$ . Then, we define the matrix

(5.42) \begin{align} \boldsymbol{P}[t;\boldsymbol{x}] = \boldsymbol{A} \big[0;\hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) \big]\boldsymbol{B} \big[t;\hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) \big]^{-1}\boldsymbol{\zeta } \big[t;\hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) \big]^\top , \end{align}

and from (4.6) the backward map is simply

(5.43) \begin{align} \hat {\boldsymbol{x}}_0(\boldsymbol{\xi },\boldsymbol{x},t) = \hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) - \boldsymbol{P}[t;\boldsymbol{x}] \boldsymbol{\xi }. \end{align}

Finally, combining (5.43) with (4.7) and making use of the Dirac-delta properties, one can obtain for the initial condition (5.41a ), the PDF of solute particle locations:

(5.44) \begin{align} f_{\boldsymbol{X}}(\boldsymbol{x};t) = \frac {\exp \left [- \frac {1}{2}\left ( \hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) - \boldsymbol{a} \right )^\top \left ( \boldsymbol{P}[t;\boldsymbol{x}]^\top \boldsymbol{P}[t;\boldsymbol{x}] \right )^{-1} \left ( \hat {\boldsymbol{{x}}}_0(\boldsymbol{x},t) - \boldsymbol{a} \right ) \right ]}{2\pi \det (\boldsymbol{P}[t;\boldsymbol{x}])}. \end{align}

This is a semi-analytic result, in the sense that once the matrix $\boldsymbol{P}$ and the initial position $\hat {\boldsymbol{{x}}}_0$ according to the backward flow map are known, the concentration field follows analytically. Because of the spatial dependence of the velocity field, these two quantities are to be calculated numerically using standard ODE solvers along streamlines. If we repeat the same procedure for (5.41b ), using the properties of both the Dirac delta and Heaviside functions, we find that

(5.45a) \begin{align} f_{\boldsymbol{X}}(\boldsymbol{x};t) &= \frac {\exp \left (-\mathrm{C}^2[b_3]\right )}{\sqrt {8\pi \big(P_{11}^2+P_{12}^2 \big)}}\frac {\textrm {erf}\left (\mathrm{D}[b_3] - \mathrm{E}[b_1]\right ) - \textrm {erf}\left ( \mathrm{D}[b_3] - \mathrm{E}[b_2] \right ) }{b_2-b_1}, \\[-9pt] \nonumber \end{align}

(5.45b) \begin{align} \mathrm{C}[b] &= \frac {b - \hat {\mathrm{x}}_{0,1}}{\sqrt {2 \big(P_{11}^2+P_{12}^2 \big)}}, \quad \mathrm{D}[b] = \frac {P_{11} P_{21} + P_{12} P_{22}}{ \det (\boldsymbol{P})}\mathrm{C}[b], \\[-9pt] \nonumber \end{align}

(5.45c) \begin{align} \mathrm{E}[b] &= \frac {b-\hat {\mathrm{x}}_{0,2}}{\det (\boldsymbol{P})}\sqrt {\frac {P_{11}^2 + P_{12}^2}{2}}, \\[9pt] \nonumber \end{align}

where we have omitted the arguments of $\boldsymbol{P}$ , its matrix elements $P_{\textit{ij}}$ and backward flow map components $\hat {\mathrm{x}}_{0,i}$ for brevity.

Figure 7. Evolution of the PDF of solute particle positions $f_{\boldsymbol{X}}$ in a heterogeneous medium for the initial condition (5.41a ) for Péclet number s (a) $ \textit{Pe}=10^4$ and (b) $ \textit{Pe}=10^3$ . The plots show the PDF reconstructed from Monte Carlo particle tracking simulations using histograms with bin side (a) $h=l_c/50$ and (b) $h=l_c/20$ for several times (contour maps), and the evolution of the PDF along the locations of the particle trajectory starting at ${x}_{0,1}=a_1$ , ${x}_{0,2}=a_2$ (coloured continuous line) computed with (5.44).

For the point-based injection (5.41a ), we illustrate the method by tracing the evolution of the PDF $f_{\boldsymbol{X}}$ along the corresponding streamline, $\boldsymbol{{x}}(t;\boldsymbol{a})$ , demonstrating the ability of the hyperbolic formulation to target particular locations of interest. This is shown in figure 7(a) for $ \textit{Pe}=10^4$ and figure 7(b) for $ \textit{Pe}=10^3$ , along with the classical particle tracking solutions for the concentration profile. The background shows the permeability field and a quiver plot of the flow field.

By interpolating the PDF $f_{\boldsymbol{X}}$ at the location of the trajectory for any given time, both approaches can be compared more directly, as shown in figure 8. It should be noted that the optimal bin size (or kernel bandwidth) for reconstructing the particle tracking PDF depends non-trivially on time, space and the Péclet number, complicating the reconstruction of the PDF from particle tracking simulations. Whist it is beyond the scope of the present work to assess state-of-the-art reconstruction techniques, which can improve the particle tracking solution, figure 8 shows how typical features such as sampling noise are avoided when using the present approach. We note in this context that the accuracy of the PDF at a given location under the hyperbolic formulation does not depend on how many points are used to discretise the solution, whereas under classical particle tracking Monte Carlo method, sampling leads to an error that is typically $\sim N^{-1/2}$ . For the spatially targeted examples shown here, the computations with the present approach are faster by orders of magnitude. In this particular case, the hyperbolic formulation only requires integrating a total of two ODEs for the trajectory $\boldsymbol{{x}}(t;\boldsymbol{a})$ and the rest of the magnitudes are given by closed-form relations or quadratures: $\boldsymbol{\zeta }$ from (3.10), $\boldsymbol{B}$ from (3.16) and $\boldsymbol{A}$ from (D5). The solution is then directly computed using (5.44).

Figure 8. Evolution of the maximum of the PDF of solute particle positions $f_{\boldsymbol{X}}$ in a heterogeneous medium for the initial condition (5.41a ), computed with the hyperbolic formulation (5.44) (continuous lines) and reconstructed from Monte Carlo particle tracking simulations with finer (dashed lines) and coarser (dash-dotted lines) histogram bin sizes. Monte Carlo results are also interpolated at the location of the trajectory from figure 7 with a zeroth-order interpolation for Péclet numbers $ \textit{Pe}=10^3$ (squares) and $ \textit{Pe}=10^4$ (circles). The tendencies shown correspond to the analysis in Dentz et al. (Reference Dentz, de Barros, Le Borgne and Lester2018).

At time $t/\tau _A\approx 30$ , a systematic error that we attribute to linearisation of the flow field can be observed in figure 8, especially for $ \textit{Pe}=10^3$ . As already noted, we do not employ relinearisation in this case. The final simulation times are small compared with the diffusion time ( $t_{\textit{end}}/\tau _D = 6\times 10^{-2}$ for $ \textit{Pe}=10^3$ ), but these errors illustrate that this is a global criterion and since the flow is heterogeneous, local variability of the flow field can lead to larger errors depending on the history of the velocity gradient tensor that each particular point of the plume samples over time.

Similarly, we compute from (5.45a ) the evolution of the normalised concentration field along streamlines associated with the line injection (5.41b ), i.e. at positions along an advected material line corresponding to this initial condition. These positions are given by $\boldsymbol{{x}}[t;\boldsymbol{{x}}_0(s^0)]$ , where $\boldsymbol{{x}}_{0}(s^0) = (b_3, \, b_1 + (b_2-b_1)s^0 )^\top$ , with $s^0\in [0, 1]$ the normalised position along the material line initially. At the initial time, we discretise the corresponding material line uniformly by taking values of $s^0$ given by $s^0_i=(i-1)/(N_s-1)$ , $i=1,\ldots , N_s$ . The normalised position $s(t)$ along the line at later times is obtained as follows. First, join the positions $\boldsymbol{{x}}_i=\boldsymbol{{x}}[t;\boldsymbol{{x}}_0(s_i^0)]$ corresponding to adjacent $s_i(t)$ by straight segments of length $l_i(t)=|\boldsymbol{{x}}[t;\boldsymbol{{x}}_0(s_{i+1}^0)]-\boldsymbol{{x}}[t;\boldsymbol{{x}}_0(s_i^0)]|$ . Then, the total length up to $\boldsymbol{{x}}_i$ is $L_i(t) = \sum _{j=1}^{i-1} l_i(t)$ and the corresponding normalised length at $\boldsymbol{{x}}_i$ is $s_i(t) = L_i(t)/L_{N_s}(t)$ . The continuous version $s(t)$ can be constructed by interpolating along the straight line segments. If desired, it could be used to resample the line uniformly as it evolves in time.

To post-process the particle tracking simulations, we employ square bins of constant side $h$ for a given time but adapt it for different times as the plume evolves and changes its characteristic spatial size. For the hyperbolic formulation, we subdivide the material line using $N_s=5 \times 10^3$ points, whereas for the particle tracking simulations, we use $N=10^6$ particles. Results at different times for $ \textit{Pe}=10^3$ obtained using the hyperbolic formulation are shown in figure 9. A visual comparison with classical particle tracking results is shown in figure 10 for two different times and Péclet numbers $ \textit{Pe}=10^3$ and $10^4$ . A more quantitative comparison with the results from (5.45a ) with the reconstructed PDF interpolated at the locations of the material line can be found in figure 11, plotted along the normalised position along the material line, showing good agreement. Similar to the case of the point injection, a total of $2N_s$ ODEs are required to compute the trajectories that define the material line; then, $\boldsymbol{\zeta }$ is computed from (3.10), $\boldsymbol{B}$ from (3.16) and $\boldsymbol{A}$ from (D5); finally, the closed-form relation (5.45a ) provides the normalised concentration profile. As discussed for figure 8, these results illustrate the benefits of the present formulation in avoiding sampling effects associated with reconstructing the concentration field from classical particle tracking, although we note again that we refrain from considering more involved reconstruction methods.

Figure 9. Evolution of the PDF of solute particle positions $f_{\boldsymbol{X}}$ in a heterogeneous medium for the initial condition (5.41b ) and Péclet number $ \textit{Pe}=10^3$ , computed with (5.45a ) along the locations of a material line at times $t/\tau _A=0$ , $2.7$ , $24.9$ and $60$ .

Figure 10. Evolution of the PDF of solute particle positions $f_{\boldsymbol{X}}$ in a heterogeneous medium for the initial condition (5.41b ) and Péclet number $ \textit{Pe}=10^3$ , computed with (5.45a ) for (a) $t=2.7 \tau _A$ and (b) $t=24.9 \tau _A$ ; and with Monte Carlo particle tracking simulations reconstructed with histograms for (c) $t=2.7 \tau _A$ with bin side $h=l_c/50$ and $309 \times 394$ bins, and (d) $t=24.9 \tau _A$ with $h=l_c/20$ and $750\times 245$ bins.

Figure 11. Evolution of the PDF of solute particle positions $f_{\boldsymbol{X}}$ in a heterogeneous medium for the initial condition (5.41b ), computed with (5.45a ) (dotted green line) and Monte Carlo particle tracking simulations (continuous red line) for (a) $t=2.7\tau _A$ and $ \textit{Pe}=10^3$ , (b) $t=24.9\tau _A$ and $ \textit{Pe}=10^3$ , (c) $t=2.7\tau _A$ and $ \textit{Pe}=10^4$ , and (d) $t=24.9\tau _A$ and $ \textit{Pe}=10^4$ . To reconstruct the PDFs from the particle tracking simulations, we use histograms with a bin side of (a) $h=l_c/50$ , (b) $h=l_c/20$ , (c) $h=l_c/100$ and (d) $h=l_c/50$ ; we then interpolate the PDF values at the locations of the advected material line from figure 10(a,b) with a zeroth-order interpolation scheme and plot it along $s(t)$ accordingly.