Mixing as a correlated aggregation process

Abstract

Mixing describes the process by which solutes evolve from an initial heterogeneous state to uniformity under the stirring action of a fluid flow. Fluid stretching forms thin scalar lamellae that coalesce due to molecular diffusion. Owing to the linearity of the advection–diffusion equation, coalescence can be envisioned as an aggregation process. Here, we demonstrate that in smooth two-dimensional chaotic flows, mixing obeys a correlated aggregation process, where the spatial distribution of the number of lamellae in aggregates is highly correlated with their elongation, and is set by the fractal properties of the advected material lines. We show that the presence of correlations makes mixing less efficient than a completely random aggregation process because lamellae with similar elongations and scalar levels tend to remain isolated from each other. We show that correlated aggregation is uniquely determined by a single exponent that quantifies the effective number of random aggregation events. These findings expand aggregation theories to a larger class of systems, which have relevance to various fundamental and applied mixing problems.

1. Introduction

Mixing of solutes by the stirring action of a fluid flow is ubiquitous to many natural and industrial processes (Ottino 1990). The evolution of the solute concentration $c$ in an incompressible velocity field $\boldsymbol {v}$ is governed by the conservation equation

(1.1)\begin{equation} \partial_t c + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}c = \kappa\,\nabla^2 c, \end{equation}

with $\kappa$ the molecular diffusivity. While advected by the flow, the solute also mixes with its surrounding due to the irreversible effect of molecular diffusion, and concentration tends to homogenize.

Although fully linear, (1.1) is characterized by rich spatio-temporal behaviours that have garnered considerable attention over the years (Rothstein, Henry & Gollub 1999; Warhaft 2000; Falkovich, Gaweȩdzki & Vergassola 2001). As illustrated in figure 1(a), an initial blob of passive solute repeatedly stirred in a chaotic flow produces scalar filaments that stretch exponentially in time (Aref 1984). Their thickness is blocked at the so-called Batchelor scale, at which molecular diffusion compensates the effect of stretching on scalar gradients (Batchelor 1959). Owing to flow incompressibility, these filaments are also exponentially compressed onto each other, leading to their diffusive coalescence (figure 1b). These three processes – stretching, diffusion and coalescence – are key to quantifying and predicting the statistics of $c$ in the homogenizing mixture.

Figure 1. (a) Mixing of a diffusive scalar by a random stirring protocol (time sequence top to bottom), evidencing the apparition of stretched scalar filaments (adapted from Villermaux 2012). (b) Blow-up on the coalescence of neighbouring filaments under the action of compression (adapted from Duplat & Villermaux 2008). (c) Concentration profile of a scalar field showing the coexistence of solitary filaments and bundles of filaments. The scalar concentration $c$ is obtained by the superposition of individual lamellae in a bundle of size $s_a$. All lamellae have a Gaussian shape with decaying maximum concentration $\theta _i$ (see (2.14)) and width tending to $s_B$.

Filament stretching and diffusion can be quantified effectively in a Lagrangian frame aligned with the local directions of elongation and compression of the flow field (Ranz 1979; Balkovsky & Fouxon 1999). Such an approach, referred to as Lagrangian stretching or lamellar theory, provides a theoretical link between the stretching statistics of the flow (its advective properties) and the evolution of concentration in isolated lamellae. It yields good estimates of the decay of the scalar concentration variance (Haynes & Vanneste 2005; Tsang, Antonsen & Ott 2005) and of the shape of the scalar probability density function (p.d.f.) in a range of chaotic flows (Meunier & Villermaux 2010). However, it is a priori limited to the early stage of mixing, when scalar filaments evolve far from each other, well before coalescence (Fereday et al. 2002; Villermaux & Duplat 2003).

Interestingly, the Lagrangian stretching framework has been shown (Haynes & Vanneste 2005; Tsang et al. 2005) to provide correct mixing rate predictions even at asymptotic time, for smooth flow in the so-called Batchelor regime, that is, when scalar gradients persist only at small scales. Indeed, in periodic domains, asymptotic mixing rates have been shown to be controlled either ‘locally’ by the stretching of material lines, or ‘globally’ by macro-dispersive flow properties, depending on the domain size with respect to the velocity correlation length. In the ‘local’ case, it is still unknown why the Lagrangian stretching framework would remain accurate after coalescence time, when stretched filaments are not isolated from each other. In addition, it is unknown how stretching statistics control the shape of the asymptotic scalar p.d.f. Some theoretical results exist on the limiting shape of the scalar p.d.f. (Sinai & Yakhot 1989), but owing to closures that are not trivial to relate to physical properties of the flow (Sukhatme 2004). In turn, Pierrehumbert (1994) and followers (Rothstein et al. 1999; Sukhatme & Pierrehumbert 2002) have documented the self-similarity of asymptotic scalar p.d.f.s through the emergence of a ‘strange eigenmode’, but a physically based modelling of the shape of the scalar p.d.f. is still an open question.

Recognizing the additive nature of scalar coalescence in bounded domains, Villermaux & Duplat (2003) have proposed describing the late-time evolution of concentration p.d.f.s by a random aggregation mechanism, whereby the concentration field results from a sum of independent filaments that individually obey Lagrangian stretching dynamics, and that overlap by diffusion. This random aggregation scenario accurately captures the shape of asymptotic scalar p.d.f.s (Duplat & Villermaux 2008; Meunier & Villermaux 2010) if one uses the rate of filament aggregation as a fitting parameter. However, no independent measure of filament aggregation in chaotic flows has been obtained so far to validate the random aggregation hypothesis. In addition, it is still unclear if aggregation processes tend to accelerate or decelerate mixing compared to stretching.

Recently, Heyman, Lester & Le Borgne (2021) have proposed that scalar aggregation could obey a correlated process, where poorly stretched regions of the flow containing weakly mixed solutes remain weakly aggregated in proportion. Such intimate correlation between stretching and aggregation has the advantage of explaining the persistent role of Lagrangian stretching statistics on asymptotic scalar mixing that is largely observed in the Batchelor regime (Tsang et al. 2005).

Here, we develop a general aggregation theory for scalar mixing in smooth, two-dimensional and fully chaotic flows in the Batchelor regime. The theory is based on the description of the spatial clustering of advected material lines, and its correlation with local elongation rates. The theory allows us to describe the asymptotic evolution of the scalar p.d.f. and its moments at any Péclet number based on the sole knowledge of Lagrangian stretching statistics and fractal dimensions of the advective material lines. It may be generalized to a wider range of flows and boundary conditions.

The paper is organized as follows. In § 2, we recall important concepts and variables pertaining to Lagrangian stretching and aggregation frameworks. In § 3, we present the two-dimensional incompressible chaotic flows used to illustrate and validate the general aggregation theory, and their numerical resolution. In § 4, we focus on the geometry of material lines advected by a chaotic flow, and we derive the evolution of the p.d.f. of the number of advective filaments in aggregates, as well as the statistics of elongation in aggregates. In § 5, we describe how the scalar concentration p.d.f. and its moment can be obtained from the knowledge of both the distribution of aggregate size and filaments statistics in aggregates.

2. Lagrangian stretching and aggregation

2.1. Stretching and dilution

We focus on smooth, two-dimensional, incompressible and time-varying flows, for which there exist two Lyapunov exponents of opposite signs. Advection thus creates elongated one-dimensional structures, called lamellae or filaments, that are the backbone of scalar mixing.

The temporal evolution of single scalar filaments is trackable in a Lagrangian frame with a coordinate system $(x,y)$ advected with the flow and aligned with the directions of compression ($x$) and elongation ($y$) (Ranz 1979; Villermaux 2019). Because of their elongated shape, the concentration of lamellae is approximately constant in the $y$ direction. Thus $\partial _y c$ is negligible compared to $\partial _x c$, and the two-dimensional advection–diffusion problem (1.1) simplifies to a one-dimensional advection–diffusion equation

(2.1)\begin{equation} \partial_t c + u(x)\,\partial_x c = \kappa\,\partial^2_x c, \end{equation}

with $u=-x\,\lambda (t)$ the velocity at which solute particles are compressed in the $x$ direction, and $\lambda (t) \geq 0$ the stretching rate. In two-dimensional incompressible flows, the stretching rate $\lambda (t)$ in the $y$ direction leads to a compression rate $-\lambda (t)$ in the $x$ direction. In chaotic flows, the mean stretching rate $\mu _\lambda$ is positive and equal to the Lyapunov exponent. Note that the notation

(2.2)\begin{equation} \mu_{{\bullet}}=\langle \,\bullet \rangle \end{equation}

is used throughout the paper for a mean quantity. We define the total lamellar elongation as

(2.3)\begin{equation} \rho(t)=\exp\left(\int_0^t \lambda(t')\,\text{d}t'\right). \end{equation}

Ranz (1979) showed that (2.1) transforms to the simple diffusion equation

(2.4)\begin{equation} \partial_\tau c = \partial^2_\xi c \end{equation}

if time is rescaled by

(2.5)\begin{equation} \tau(t)=\frac{\kappa}{s_0^2} \int_0^t \rho(t')^2\,\text{d}t', \end{equation}

and space is rescaled by $\xi =x \rho /s_0$, with $s_0$ the initial lamellar width. For a Gaussian initial condition $c(\xi,0)=\theta _0\exp (-\xi ^2)$, the solution of (2.3) is

(2.6)\begin{equation} c(\xi,\tau)=\frac{\theta_0}{\sqrt{1+4 \tau}} \exp(-(\xi/\sqrt{1+4 \tau})^2). \end{equation}

Equivalently, in the original Lagrangian coordinate system $(x,y)$,

(2.7)\begin{equation} \varTheta(x,t)=\theta(t)\exp(-(x/s(t))^2), \end{equation}

where $\theta$ is the maximum concentration of the lamella, following

(2.8)\begin{equation} \theta(t)=\frac{\theta_0}{\sqrt{1+4 \tau(t)}}, \end{equation}

and $s$ is the unit lamellar thickness, following

(2.9)\begin{equation} s(t)=\frac{s_0 \sqrt{1+4 \tau(t)} }{\rho(t)}. \end{equation}

Note that the scalar mass per unit length under a lamella simplifies to

(2.10)\begin{equation} m(t) = \sqrt{\rm \pi}\,\theta\,s(t) = \sqrt{\rm \pi}\,s_0 \theta_0\,\rho^{-1}(t). \end{equation}

Since in random flows the stretching rate is a random variable of time, the elongation of lamellae and the rescaled time are also randomly distributed. An approximation of the statistics of $\tau$ was proposed (Meunier & Villermaux 2010; Lester, Dentz & Le Borgne 2016) as

(2.11)\begin{equation} \tau \approx \frac{\kappa }{ 2 s_0^2}\,\frac{ t }{ \log \rho}\,(\rho^2-1), \end{equation}

recognizing the fact that the last stretching events have a predominant weight in the stochastic integral (2.4). At large times, $\rho \gg 1$ and $\log \rho /t \to \mu _\lambda$. Thus

(2.12)\begin{equation} \tau \to \frac{\kappa}{ 2 \mu_\lambda s_0^2}\,\rho^{2} = \frac{1}{4} \left( \frac{s_B}{s_0}\,\rho\right)^2, \end{equation}

with the so-called Batchelor scale

(2.13)\begin{equation} s_B=\sqrt{2 \kappa / \mu_\lambda}. \end{equation}

Note that in his seminal paper, Batchelor (1959) uses the mean turbulent strain rate instead of the Lyapunov exponent $\mu _\lambda$ to quantify fluid stretching. Thus filaments dilute at large times, with a peak concentration following

(2.14)\begin{equation} \theta (t) \to \frac{\theta_0 s_0}{ s_B}\,\rho^{-1}(t), \end{equation}

while their thickness tends to the Batchelor scale $s \to s_B$. In contrast, the length of filaments grows as

(2.15)\begin{equation} L(t) = \ell_0\,\mu_\rho (t), \end{equation}

where $\ell _0$ is the initial filament length. Equation (2.14) shows that the concentration of single diffusing filaments is inversely proportional to their elongation. Thus describing the advective stretching statistics is sufficient to describe the mixing of passive scalars. Indeed, the spatial variance of an elongated material line in a flow domain of area $\mathcal {A}$ can be obtained by integration of the lamellar concentration (2.14) squared along the filament path (Meunier & Villermaux 2010), giving

(2.16)\begin{equation} \sigma_c^2 \sim \frac{1}{\mathcal{A}}\,L(t)\,\mu_\rho^{-2} = \frac{1}{\mathcal{A}}\,\mu_\rho^{-1}. \end{equation}

Assuming a normal distribution of stretching with mean $\mu _\lambda t$ and variance $\sigma ^2_\lambda t$, the scalar variance decays can be estimated from the stretching statistics. For flows with moderate stretching variability, where $\sigma ^2_\lambda \leq \mu _\lambda$, we have

(2.17)\begin{equation} \sigma_c^2 \sim \exp(-(\mu_\lambda-\sigma_\lambda^2/2)t). \end{equation}

In contrast, for flows with larger stretching variability where $\sigma ^2_\lambda \geq \mu _\lambda$, the value of $\mu _\rho ^{-1}$ is entirely controlled by small elongations with $\rho \leq 1$ (Balkovsky & Fouxon 1999). In that case, $\mu _\rho ^{-1}\sim \exp (-\mu _\lambda ^2/(2\sigma _\lambda ^2)t)$.

It is important to note that the Lagrangian stretching approximation is correct in the limit of one-dimensional filaments, that is, when the curvature radius is much larger than the diffusing filament width $s_B$. In bounded flow domains, there exist regions where high curvatures develop (Tang & Boozer 1996), which are known (Thiffeault 2004) to be associated with weak stretching rates (low elongations) and thus low mixing (high concentrations). This questions the ability of the Lagrangian stretching framework to capture the tails of the scalar concentration p.d.f. However, at small $s_B$ (small $\kappa$), the impact of these regions on global mixing rates should remain limited. This is confirmed by the good agreement of the Lagrangian stretching theory with observed mixing rates (Haynes & Vanneste 2005).

2.2. Aggregation

The production of curvature tends to produce folds in the material line (figure 1b). The distance between separated parts of the elongated material line is thus reduced exponentially via fluid compression. This mechanism creates a highly foliated structure at late time (figure 1a). Individual filaments are thus no longer isolated, and start to coalesce on length scales comparable to their width $s_B$. This so-called aggregation regime (Villermaux & Duplat 2003) has two properties.

First, filaments tend to accumulate locally due to exponential flow compression. Bundles of aggregated lamellae are thus formed by individual filaments sharing the same diffusive neighbourhood of size ${\sim }s_B$. The time at which aggregation initiates is when the total area of a lamella of length $L(t)$ and width $s_B$ is greater than the available area of the flow domain $\mathcal {A}$ (Garrett 1983), e.g.

(2.18)\begin{equation} L(t)\,s_B \gtrsim \mathcal{A}. \end{equation}

The mean number of lamellae $n$ in each bundle is thus

(2.19)\begin{equation} \mu_{n}(t) \sim \frac{L(t)\,s_B }{ \mathcal{A}}. \end{equation}

Assuming a constant stretching rate $\lambda$, the length of the lamella is $L(t) = \ell _0 \exp (\lambda t)$, and the coalescence time $t_c$ at which (2.18) is first fulfilled is

(2.20)\begin{equation} t_c \sim \frac{1}{\lambda} \log\left(\frac{\mathcal{A}}{\ell_0 s_B}\right). \end{equation}

Second, the linearity of the advection–diffusion equation implies that scalar concentration fields can be decomposed into a sum over the concentration profiles of solitary lamellae (Le Borgne et al. 2017). The scalar concentration $c$ at a given position can be constructed as the sum of concentrations $\theta _i$ of individual lamellae $i$ (and defined by (2.14)) present in a small neighbourhood around this position, of size $s_a$ comparable to the asymptotical size of the lamellae, $s_B$, that is,

(2.21)\begin{equation} c(t) \sim \sum_{i=1}^{n(t)} \theta_{i}(t)\sim \sum_{i=1}^{n(t)} \rho_{i}^{-1}(t), \end{equation}

where the asymptotic relation (2.14) between $\theta$ and $\rho$ was used. Through (2.21), the aggregation process offers an appealing way to model the statistics of the scalar concentration field $c$ via knowledge of the aggregated number of filaments $n$ at the aggregation scale $s_a$, their elongations $\rho _i$, and the possible correlations between these two variables.

Assuming no correlations between $n$ and $\rho _i$, and among $\rho _i$, leads to a random aggregation scenario. The scalar concentration $c$ in a bundle is thus formed by the sum of independent and identically distributed random variables, following the solitary filament concentration p.d.f. The scalar concentration p.d.f. $P_c(c,t)$ thus results from the $n$-convolution of the isolated lamella concentration p.d.f. $P_{1/\rho }(\rho ^{-1},t)$, with the mean number of aggregations $\mu _n$ given by (2.19) (Duplat & Villermaux 2008). If $P_{1/\rho }$ is exponential or gamma distributed, then $P_c$ is a gamma distribution. Note that if $P_{1/\rho }$ is log-normal, then the $n$-convolution p.d.f. is not explicit, although its moments are given by the central limit theorem (Schwartz & Yeh 1982). This random conjecture was shown (Duplat & Villermaux 2008) to correctly describe mixing in turbulent flows, for which the velocity cascade over a large range of scales favours the decorrelation of filament stretching histories.

Below the characteristic velocity length scale (the Batchelor regime), however, filament aggregation seems not to obey completely random dynamics. To see that, it is useful to consider the decay of scalar variance. In the random hypothesis, by the central limit theorem, scalar variance decays as the inverse of the aggregation number, that is, the inverse of the material length (2.19):

(2.22)\begin{equation} \sigma^2_c = \frac{\mu_c^2}{\mu_n} \sim 1/ L(t). \end{equation}

Assuming a normal distribution of stretching rates, $L \sim \exp ((\mu _\lambda +\sigma _\lambda ^2/2)t)$ and thus the asymptotic scalar variance decay exponent is $\mu _\lambda +\sigma _\lambda ^2/2$, larger than the decay exponent of solitary strips ($\mu _\lambda -\sigma _\lambda ^2/2$, (2.17)). This result contradicts observations suggesting the same decay exponent before and after aggregation time in the Batchelor regime (Fereday et al. 2002; Tsang et al. 2005).

We attribute this failure to the existence of correlations in the aggregation process below the characteristic velocity length scale. Indeed, in incompressible flows, lamellar elongation is always balanced with transverse compression (figure 1b), which attracts neighbouring lamellae onto each other. Thus highly stretched lamellae tend to be highly aggregated. Conversely, weakly stretched lamellae have also experienced little compression, thus remaining isolated from the bulk, and their evolution well described by the isolated lamellar theory. Since these lamellae bear high concentration levels, they must dominate the statistics of scalar fluctuations at late time. Such correlations may explain why the Lagrangian stretching prediction for the scalar dissipation remains accurate long after aggregation time.

Evidence of correlated aggregation was also observed experimentally during the chaotic mixing of two dye blobs (Duplat, Jouary & Villermaux 2010). If injected in a concentric manner, the blobs locally experience similar stretching rates and aggregate in a correlated manner. In contrast, when placed at a certain distance from each other (larger than the characteristic length scale of the velocity field), their concentrations obey random aggregation rules. This observation suggests that stretching and aggregation are strongly correlated below the characteristic velocity length scale, while they become uncorrelated above this scale.

The remainder of the paper aims to establish the laws governing correlated aggregation for scalar mixing by smooth chaotic flows, and to deduce their impact on the evolution of scalar concentration p.d.f.

3. Numerical simulation of synthetic chaotic flows

We focus on scalar mixing in the so-called Batchelor regime (Haynes & Vanneste 2005) where the velocity field is smooth, and present a single length scale, much greater than the typical scalar fluctuation scale. The domain size and the velocity length scale are of comparable magnitude, $\mathcal {L}\sim \mathcal {L}_v \sim 1$ – the so-called ‘local’ mixing regime described in Haynes & Vanneste (2005). Thus no scalar gradients can persist at scales larger than the velocity scale.

To illustrate the geometrical features of aggregation, we use two synthetic chaotic transformations that comply with these conditions, namely the baker map and the random sine flow (figure 2). These are simple sequential advective maps that have been widely used to investigate the properties of chaos (Finn & Ott 1988; Ott & Antonsen 1989; Giona, Cerbelli & Adrover 2001; Tsang et al. 2005; Meunier & Villermaux 2010, 2022). We recall their definitions below.

Figure 2. (a) Transformations operated by the incompressible baker map with parameter $a$. First, the domain is cut horizontally at $y=a$, where $a$ is a constant between $0$ and $0.5$. Uniform fluid compression operates on the domain parts $y< a$ and $y>a$ with $a$ and $1-a$, respectively. Then vertical stretching occurs with a factor $1-a$ and $a$ in these two regions, preserving the total area. (b) Transformations operated by the sine flow with amplitude $A$. The flow is an alternation of horizontal and vertical sinusoidal velocity waves with amplitude $A$ and period $2{\rm \pi}$. Random phases are chosen at each time period so that the flow is fully chaotic.

3.1. Incompressible baker map

The incompressible baker map (Ott & Antonsen 1989; Wonhas & Vassilicos 2002) is a discontinuous transformation that operates on a two-dimensional periodic domain $[0,1]\times [0,1]$. The transformation is

(3.1)\begin{equation} \left.\begin{gathered} x_{t+1} = \left\lbrace \begin{array}{@{}ll} a x_t & \text{if } y_t < a, \\ 1-(1-a) x_t & \text{if } y_t > a, \end{array}\right.\\ y_{t+1} = \left\lbrace \begin{array}{@{}ll} y_t / a & \text{if } y_t < a, \\ (1-y_t)/(1-a) & \text{if } y_t > a, \\ \end{array} \right. \end{gathered}\!\right\} \end{equation}

where $a \in [0,0.5]$ is a parameter controlling the heterogeneity of the map. A visual sketch of the map operation is shown in figure 2a.

An advantage of the baker map is that purely vertical scalar patterns (for which $c(x,y)=f(x)$) remain one-dimensional after application of the map, thus simplifying the problem to a single dimension. This simplicity allows for the analytical derivation of many features of the map, as we will show later. Another advantage is that it is possible to explore a wide range of stretching heterogeneity by varying $a$ between 0 and 0.5. Indeed, the first two moments of stretching rate in the baker map are

(3.2)$$\begin{gather} \mu_\lambda/ t =-a\log(a) - (1-a) \log(1-a) , \end{gather}$$

(3.3)$$\begin{gather}\sigma^2_\lambda / t = a(1-a)(\log(1-a)-\log(a))^2 . \end{gather}$$

Thus for $a=0.01$, $\sigma _\lambda ^2/\mu _\lambda = 3.7$, while for $a=0.49$, $\sigma _\lambda ^2/\mu _\lambda =5.7 \times 10^{-4}$. It is important to note that this map involves discontinuous transformations, or ‘cuts’, that are absent in continuous flows such as turbulence but are common in flows through porous media (Lester, Metcalfe & Trefry 2013). This map also prevents the formation of folds and cusps, where the one-dimensional hypothesis of the Lagrangian stretching framework is violated. The correspondence between (1.1) and (2.1) is thus exact in this map.

3.2. Random sine flow

The random sine flow (Pierrehumbert 1994; Haynes & Vanneste 2005) is a continuous transformation operating on a periodic domain $[0,1]\times [0,1]$ (figure 2b). The flow is periodic in time and space, and reads for a given time period $t$ as

(3.4)\begin{equation} \left.\begin{gathered} y_{t'+\delta t} = y_{t'} +A\,\delta t \left\lbrace \begin{array}{@{}ll} \sin( 2{\rm \pi} x_t' + \phi_t) & \text{ for } t< t'< t+1/2, \\ 0 & \text{ for } t+1/2< t'< t+1, \end{array}\right.\\ x_{t'+\delta t} = x_{t'} + A\,\delta t \left\lbrace \begin{array}{@{}ll} 0 & \text{ for } t< t'< t+1/2, \\ \sin( 2{\rm \pi} y_t' + \psi_t) & \text{ for } t+1/2< t'< t+1, \\ \end{array} \right. \end{gathered}\!\right\} \end{equation}

where the amplitude $A$ is a positive constant, $\phi _t,\psi _t$ are random phases that change at each time period $t$, and $\delta t$ is the time step. The flow velocity having a single component, incompressibility is enforced. Scalar transport is continuous and considered on a periodic domain $[0,1]\times [0,1]$, ensuring a local control of mixing rates (Haynes & Vanneste 2005).

As most random chaotic flows, the elongation of material lines in sine flows approximately follows a log-normal distribution with parameters $\mu _\lambda t$ and $\sigma ^2_\lambda t$ that depend on the amplitude $A$. The stretching heterogeneity is much less variable than in the baker map, with ratio $\sigma ^2_\lambda /\mu _\lambda$ ranging from $1$ when $A\to 0$ to $\sigma ^2_\lambda /\mu _\lambda \approx 0.6$ for $A=1.8$. The stretching statistics of random sine flows can be found in Meunier & Villermaux (2022). In Appendix A, we recall useful results concerning the distribution of filament elongation and its moments in the sine flow and the baker map.

3.3. Numerical methods

The hypothesis used throughout the paper is that the concentration field obeying (1.1) is well approximated by a local summation over elementary lamellar concentrations, each of them individually following (2.1). This implies that the statistics of concentrations and their temporal evolution can in principle be inferred from the statistics of lamellar elongation and aggregation. Thus instead of directly solving the two-dimensional advection–diffusion equation (1.1), our approach consists in computing elongation of a deforming Lagrangian support, advected by the flow (figure 3).

Figure 3. Advection of a material filament in the sine flow with $A=0.8$, for (a) $t=3$, (b) $t=5$ and (c) $t=10$.

Numerically, we advect a material line of initial length $\ell _0$ in the velocity field. In practice, the material line is defined by a series of particles $\boldsymbol {x}_i$ linked by segments. In the alternated random sine wave, the positions of particles are tracked through time $t$ via an explicit Euler scheme (figure 3). The particles have constant velocities during a time step $\textrm {d}t =1/2$, so this scheme is exact. The local elongation of segments is estimated as

(3.5)\begin{equation} \rho_i(t) = \rho_i(t-\text{d}t)\,\frac{\|\boldsymbol{x}_i(t+\text{d}t)- \boldsymbol{x}_{i+1}(t+\text{d}t)\|}{\|\boldsymbol{x}_i(t)-\boldsymbol{x}_{i+1}(t)\|}. \end{equation}

The material line is refined after each time step to maintain a maximum local distance between points of $\textrm {d}\kern0.06em x \leq 10^{-3}$.

In the baker map, we start with a single filament aligned with the $y$ direction, so that we need to track only its coordinate $x$. Each operation of the map doubles the number of filaments, with the first half being elongated by a factor $1/a$, and the second half by a factor $1/(1-a)$. We thus keep track of the growing number of filament positions $x_i$ and their elongations $\rho _i$.

The elongated and folded filament is tracked up to the advection time where the total length is $L=10^7\ell _0$, a limit corresponding to our computer memory. Local elongation of the advected material line can then be used to compute its local concentration, via (2.14). Aggregated scalar fields can then be estimated by a local summation of individual lamellar concentrations, via (2.21).

3.4. Reconstruction of scalar fields

As shown by Meunier & Villermaux (2010), it is possible to reconstruct the aggregated scalar field via the superposition of elementary lamellae, represented as Gaussian ellipsoids of short axis $s$, that swipe along the advected material line. This reconstruction is exact in the baker map since no cusps are forming. In the sine flow, it is only approximate in the regions of high curvature, where diffusion is genuinely two-dimensional.

Here, we take an alternative reconstruction approach, which has the advantage of allowing analytical treatment. At late time, lamella widths tend to $s_B$ (2.13), which is the minimum scale of fluctuations of the scalar field (Batchelor 1959). Thus we construct the aggregated scalar field by binning lamellar concentrations on a regular grid of size $s_a \sim s_B$. Consider a box of surface $s_a^2$ centred in position $\boldsymbol {x}$. The aggregated concentration level in this box can be constructed from the sum of the masses per unit length $m_i$ of the $n(\boldsymbol {x})$ individual lamellae of length $\ell _i$ present in this box:

(3.6)\begin{equation} c(\boldsymbol{x}) \approx \frac{1}{s_a^2} \sum_{i=1}^{ n(\boldsymbol{x}) } s_a m_i = \frac{ \sqrt{\rm \pi}\,\theta_0 s_0}{s_a} \sum_{i=1}^{n(\boldsymbol{x})} \rho^{-1}_i, \end{equation}

where we used (2.9) for the evolution of the solute mass per unit length carried by an individual lamella at a given location. Equation (3.6) forms the base of the statistical description of aggregation intended in this paper.

In practice, we choose

(3.7)\begin{equation} s_a=\sqrt{2{\rm \pi}}\,s_B \end{equation}

for the variance of the regular grid reconstruction to match the variance of the true field (i.e. the one obtained by aggregating Gaussian lamellae). The factor $\sqrt {2{\rm \pi} }$ is obtained so that both reconstruction methods provide the same mean-squared concentration in the case of a single lamella.

Note that the aggregated concentration field can be deduced at various Péclet numbers by varying the aggregation scale $s_a$ (3.7), without recomputing the advective filament position. We define the Péclet number as the ratio of diffusive to stretching time scales, e.g.

(3.8)\begin{equation} Pe= \mu_\lambda \mathcal{L}_v^2/ \kappa. \end{equation}

Through (2.13), $s_B= \mathcal {L}_v\sqrt {2/Pe}$ and thus $s_a= 2\sqrt {{\rm \pi} }\,\mathcal {L}_v/Pe$.

This regular grid approximation of the aggregated scalar field disregards the exact Gaussian shape of lamellar concentrations and the variations of their width $s$ (2.5). However, classical box counting methods and fractal dimensions can be used to describe the geometrical and statistical features of the aggregated scalar field. Our goal is to describe the joint statistics of the variables $\rho _i$ and $n$ arising in (3.6), to infer the statistics of $c$.

The approximation (3.6) is compared to direct numerical simulations (DNS) of the advection–diffusion equation for different Péclet numbers in figure 4. The DNS are obtained with the spectral method described in Meunier & Villermaux (2022), with a $2048^2$ grid, and time step 0.1. There is a good qualitative match between the aggregated field and the DNS, a match confirmed by comparing the distributions of concentrations. Thus we conclude that the construction of the aggregated scalar field via (3.6) is able to capture the essential statistical features of the two-dimensional advection–diffusion problem, and offers a convenient way to explore its statistics.

Figure 4. Comparison of log-concentration fields obtained at time $t=10$ (fully aggregated regime): (i) with the aggregation framework (3.6), (ii) with direct numerical simulations (DNS) of (1.1), and (ii) with their p.d.f. The comparison is made for the sine flow ($A=0.8$) for two Péclet numbers, corresponding to (a) $s_a=1/150$ and (b) $s_a=1/50$.

4. Clustering properties of advected material lines

In incompressible flows, the stretching of material elements by velocity gradients is compensated by transverse compression, which causes distances between lamellar elements to decrease exponentially over time. Smaller and smaller scales are thus produced continuously by flow compression. Furthermore, in smooth chaotic flows, the typical scale of variation of velocity gradients is fixed and produces a heterogeneous stretching field for material lines. Dense (black) or diluted (white) regions of material lines are thus created at large scale in the chaotic flow (figure 5). Such heterogeneous structures then cascade to smaller scales under the action of net compression, thus creating a fractal set of one-dimensional objects (lines) clustered around their transverse direction. In two-dimensional incompressible flows, the Hausdorff dimension of this fractal set is necessarily $D_0=2$, as per the Kaplan–Yorke result (Farmer, Ott & Yorke 1983). Higher-order dimensions can be smaller than 2 if stretching is heterogeneous. To illustrate this, let us define a normalized measure $p_k$ with $k=1,\ldots, N$ defining a regular grid of bin size $\epsilon =\mathcal {L}/N$, where $\mathcal {L}$ is the domain size. For instance, $p_k$ may be defined as the local density of lamellae in the bin, e.g. $p_k=n_k/n$, where $n$ is the total number of lamellae.

Figure 5. Fractal geometry of material lines in (a–c) the sine flow ($A=0.8$) and (d–f) the baker map ($a=0.1$), observed at different scales. The red square indicates the area selected for zooming.

The fractal dimension of order $q$ of the measure $p$ is obtained with (Grassberger 1983)

(4.1a,b)\begin{equation} D_q -1 = \lim_{\epsilon \to 0} \frac{1}{q-1}\,\frac{\log I_{q}(\epsilon)}{\log \epsilon},\quad I_q(\epsilon) \equiv \sum_k^{N=\mathcal{L}/\epsilon} p_k^q, \end{equation}

where the subtraction of 1 on the left-hand side accounts for the counting of one-dimensional structures (lamellae) in a two-dimensional domain. This definition implies the following spatial scaling of the integral of the measure:

(4.2)\begin{equation} I_q(\epsilon) \sim \epsilon^{(q-1)(D_q-1)}.\end{equation}

In simple flows such as the baker map, $D_q$ can be obtained (Finn & Ott 1988) by observing the similarity properties of the map, which transfer at small scales the heterogeneity of the measure produced at large scales by a single operation of the map. Characterizing the result of one elementary operation of the map on the measure thus also informs us about the spectrum of fractal dimensions.

We consider the measure of the local number of lamellae in bin $k$, $p_k=n_k/n$. As shown in figure 2(a), an operation of the baker map doubles the total number of these lamellae, while maintaining the same local distribution of lamellae on smaller bins of sizes $a\epsilon$ for $x< a$, and $(1-a)\epsilon$ for $x>a$. This similarity allows one to compute the first fractal dimension explicitly (see Appendix B and Finn & Ott 1988):

(4.3a,b)\begin{equation} D_0=2,\quad D_1 =1+ \frac{2 \log 2}{\log (a^{-1} +(1-a)^{-1})}. \end{equation}

For random flows such as the sine flow, Ott & Antonsen (1989) argue that there exists a general relationship between stretching rate statistics and fractal dimensions as $D_q = f_q(\sigma _\lambda ^2,\mu _\lambda )$, although a closed-form solution is not always as trivial as for the baker map. We show in figure 6 that the ratio $\sigma _\lambda ^2/\mu _\lambda$ is directly related to the fractal dimension $D_1$. This suggests that the fractal clustering of material lines is closely linked to the large-scale heterogeneity of stretching rates. Since the flow is smooth, the heterogeneity created at large scales cascades to smaller scales, conserving its geometrical structure and creating a fractal geometry.

Figure 6. Relation between stretching rate mean $\mu _\lambda$ and variance $\sigma _\lambda ^2$, and fractal dimension $D_1$, in the baker map and sine flow with varying parameters $a$ and $A$.

As suggested by figure 6, the function $f_1$ is different for the baker map and the sine flow. Indeed, the ratio $\sigma _\lambda ^2/\mu _\lambda$ tends to a positive constant in the sine flow when $A\to \infty$, while $\sigma _\lambda ^2/\mu _\lambda \to 0$ in the baker map when $a\to 0.5$. This finite limit comes from the fact that the sine flow is a continuous transformation with no cutting, thus does not tend to a uniform stretching rate. On the contrary, when $A\to 0$, $\sigma _\lambda ^2/\mu _\lambda \to 1$, which is a maximum bound for the ratio in the sine flow (Meunier & Villermaux 2022), thus limiting the possible range of fractal dimensions produced by continuous chaotic flows, compared to discontinuous maps.

4.1. Spatial distribution of $n$

The spatial distribution of the number $n$ of lamella per bin size $\mathcal {A}$ (figure 7) can be obtained as follows. Comparing the mean area occupied by a filament of length $L(t)$ and aggregation scale $s_a$ to the domain surface $\mathcal {A}$, we get an estimate of the mean number of lamellae $\mu _n$ in bundles:

(4.4)\begin{equation} \mu_n \sim L(t)\,s_a / \mathcal{A}. \end{equation}

Higher moments can be obtained from a study of the fractal structure of material lines. To this end, we consider the spatial measure corresponding to the local number of lamellae in each bundle:

(4.5)\begin{equation} p_k \equiv \frac{n_k}{\displaystyle\sum_k n_k} \end{equation}

The Renyi definition (Grassberger 1983) of the fractal dimension of order 2 of this measure is

(4.6)\begin{equation} D_2 -1 \approx \frac{\displaystyle \log \sum_k p_k ^2}{ \log s_a}, \end{equation}

assuming $s_a\to 0$. Replacing (4.5) in the last expression provides

(4.7) \begin{equation} \sum_k \begin{pmatrix}\displaystyle\frac{ n_k}{\displaystyle\sum_k n_k}\end{pmatrix}^2= s_a^{D_2-1}. \end{equation}

Since $\sum _k n_k= N \mu _n$, we have

(4.8) \begin{equation} \sum_i \begin{pmatrix}\displaystyle\frac{n_k}{\displaystyle\sum_k n_k}\end{pmatrix}^2 = \frac{1}{N}\,\langle ({n}/{\mu_n})^2\rangle, \end{equation}

with $N \approx \sqrt {\mathcal {A}}/s_a$ the number of bundles in the flow domain. Since

(4.9)\begin{equation} \sigma^2_{n/\mu_n} = \langle (n/\mu_n)^2\rangle- \langle n/\mu_n\rangle^2, \end{equation}

we have

(4.10)\begin{equation} \sigma^2_{n/\mu_n} = \sqrt{\mathcal{A}}\,s_a^{D_2-2} - 1. \end{equation}

Thus the variance of $n/\mu _n$ reaches a constant at asymptotic times, which is given by the fractal dimension of order 2. The spatial variance of $n$ is then

(4.11)\begin{equation} \sigma^2_{n} = \mu_n^2 \left(\sqrt{\mathcal{A}}\,s_a^{D_2-2} - 1\right). \end{equation}

The predictions of (4.4)–(4.11) are plotted against time in figures 8 and 9 showing good agreement with simulations for a large range of aggregation scales $s_a$ and flow heterogeneity, characterized by the parameters $a$ for the baker map, and $A$ for the sine flow. Higher moments of $n$ can be obtained with similar scaling arguments and linked to fractal dimension of higher order. This is, however, out of the scope of the present study.

Figure 7. Spatial distribution of the number $n$ of lamellae in bundles defined by a regular grid of size (a) $s_a=1/200$ and (b) $s_a=1/50$, in the sine flow with parameter $A=0.5$.

Figure 8. (a) Scaling of the spatial variance of $\log n$ as a function of $a$ in the baker map and theoretical prediction, (4.10). (b) First two moments of $P_n$ through time compared to theoretical predictions, (4.4) and (4.11) in the baker map ($\mathcal {A}=1$).

Figure 9. First two moments of $P_n$ through time compared to theoretical predictions, (4.4) and (4.11) in the sine flow ($\mathcal {A}=1$).

Instead, we observed that the p.d.f. of lamellar aggregation number $P_n (n)$ is well fitted by a gamma distribution (figure 10)

(4.12)\begin{equation} P_n (n) = \frac{1}{\varGamma(k_n)\,\theta^k_n}\,n^{k_n-1} \exp(-n/\theta_n), \end{equation}

Figure 10. Plots of $P(n,t)$ for the baker map and sine flow for $s_a=1/100$. Solid lines stand for the gamma p.d.f. with theoretical moments given by (4.14), and symbols stand for numerical simulations. (a) Baker map $a=0.3$ and variable $t$. (b) Baker map for $t=20$ and variable $a$. (c) Sine flow for $s_a=1/100$ and $A=0.4$.

with $n \geq 0$, and $k_n,\theta _n$ defined by the moments of the distribution of $n$:

(4.13)$$\begin{gather} k_n = \left(\sqrt{\mathcal{A}}\,s_a ^{D_2-2} -1\right)^{-1}, \end{gather}$$

(4.14)$$\begin{gather}\theta_n = \mu_n(t)\left(\sqrt{\mathcal{A}}\,s_a ^{D_2-2} -1\right), \end{gather}$$

where $\mu _n=L(t)\,\mathcal {A}/s_a$. Note that the gamma distribution yields a power-law distribution at small $n$, with exponent $k_n-1$ going from $-1$ to infinity with increasing $D_2$. For small $D_2$ (strong stretching heterogeneity), a significant part of the probability is concentrated at small $n$, thus in non-aggregated regions of the flow. We will show later that this can affects the value of negative moments of $n$.

4.2. Local correlations between $n$ and $\rho$

Having given a precis of the statistics of the bundle size $n$, we now explore the local statistics of lamellar elongations in these bundles. To this end, we define the conditional averaging operator acting in lamellae located in the local neighbourhood of size $s_a$ by

(4.15)\begin{equation} \mu _{X\,|\, n} = \frac{1}{n} \sum_{i=1}^{n} X_i, \end{equation}

where $X$ is a Lagrangian variable transported by lamellae and $n$ is the number of lamellae aggregated in the bundle. The spatial variability of averaged variables is shown in figure 11 for simulations in the baker map and sine flow. The remainder of this section is dedicated to uncovering the behaviour of conditional moments of elongation conditioned to $n$ (e.g. $\mu _{\rho ^{-q}\,|\, n}$). Section 5 will then be dedicated to deriving unconditional probabilities by averaging on the distribution of $n$.

Figure 11. Simulation of aggregation statistics in the baker map ($a=0.3$) and sine flow ($A=0.5$) for $s_a=1/50$: (i) number of lamellae $n$ in bundles, (ii) mean of log-elongation in bundles, and (iii) sum of lamellar concentrations in bundles.

We plot in figure 12 the joint probability $P(n, \mu _{X \,|\, n})$ obtained in the baker map and the sine flow for $X=\rho ^{-1}$ and $X=\log \rho$, the inverse of elongation and the log-elongation of lamella, respectively. Figure 12 suggests that the following scaling holds in both flows:

(4.16)\begin{equation} \log n \sim- \log\mu_{\rho^{-1}\,|\, n}, \end{equation}

which confirms the strong correlation between the number of lamellae in aggregates and their elongation. Indeed, for large times, $c$ must tend to the conserved average scalar concentration $c \to \mu _c$. Thus, owing to (3.6), we must have

(4.17)\begin{equation} n \sim 1/\mu_{\rho^{-1} \,|\, n }, \end{equation}

a scaling that we confirm numerically (figure 12).

Figure 12. Joint p.d.f. (grey scale) of the number of lamellae in a bundle of size $s_a=1/200$, and (i) their mean inverse elongation, (ii) their mean log-elongation, for (a) baker map ($a=0.1$, $t=24$, $D_1=1.57$) and (b) sine flow ($A=0.8$, $t=10$, $D_1=1.74$). The theoretical scalings of the measures (i) and (ii), given by (4.16) and (4.18), respectively, are plotted as solid red lines with the slope indicated in the legend. Dashed red lines are guides for the eye.

Figure 12 also suggests that

(4.18)\begin{equation} \log n \sim (D_1-1)\,\mu_{\log \rho \,|\, n }, \end{equation}

where $D_1$ is the information dimension (Ott & Antonsen 1989) of the measure $n$, given by (4.3a,b) for the baker map. Equation (4.18) can be derived analytically in the case of the baker map. Indeed, by the action of the map, the total number of lamellae increases as $\log n = t \log 2$, while the mean log-elongation of these lamellae is $\mu _{\log \rho } =t ( -\log a - \log (1-a) )/ 2$, leading to a constant ratio

(4.19)\begin{equation} \frac{\log n}{\mu_{\log \rho }}= \frac{2\log 2}{\log a + \log(1-a)}, \end{equation}

which is exactly the value of $D_1-1$ (4.3a,b). Assuming that the partition between $\log n$ and $\mu _{\log \rho }$ is preserved at small scales in each bundle, we have

(4.20)\begin{equation} \log n = ({D_1 - 1})( \mu_{\log {\rho} \,|\, n} - \mu_{\log \rho_c}), \end{equation}

with $\mu _{\log \rho _c}$ a constant standing for the mean elongation at coalescence time.

The constant $\mu _{\log \rho _c}$ can be estimated by comparing the surface covered by the filament at the aggregation scale, $\mathcal {S}=s_a \rho _c L_0$, with the domain area $\mathcal {A}$. The first aggregation event occurs when $\mathcal {S} \approx \mathcal {A}$, that is, when the elongation is $\rho _c \approx \mathcal {A}/(s_a \ell _0)$. Equation (4.20) is indeed verified for the baker map and sine flow with various parameters $a$ and $\mathcal {A}$, with $\mu _{\log \rho _c} = - \log s_a$.

4.3. Distribution of $\log \rho$ in a bundle of size $n$

The two scaling laws $n \sim \mu _{\rho ^{-1}\,|\, n}$ and $\log n\sim (D_1-1)\,\mu _{\log \rho \,|\, n}$ provide key information about the heterogeneity of lamellar elongations inside bundles. Since the ensemble distribution of elongation $P_\rho (\rho )$ has a log-normal shape, we assume that the distribution of elongations inside bundles, denoted $P_{\rho \,|\,n}$, is also log-normally distributed. This implies that $\log \rho$ is normally distributed in bundles, with mean

(4.21)\begin{equation} \mu_{\log \rho\,|\,n} \sim (D_1-1)^{-1} \log n. \end{equation}

Since $\log \mu _{\rho ^{-1} \,|\, n}= -\mu _{\log \rho \,|\,n }+\sigma ^2_{\log \rho \,|\,n}/2 \sim -\log n(\boldsymbol {x})$, the variance of log-elongation in bundles at large $n$ must be

(4.22)\begin{equation} \sigma^2_{\log \rho\,|\,n} \sim \frac{2(2-D_1)}{D_1-1} \log n. \end{equation}

We report in figure 13 the simulated scaling $\mu _{\log \rho \,|\,n}/\log n$ and $\sigma ^2_{\log \rho \,|\,n}/\log n$ obtained asymptotically at large times. When $D_1\to 2$, $\mu _{\log \rho \,|\,n}/\log n \to 1$ while $\sigma ^2_{\log \rho \,|\,n}/\log n \to 0$, meaning that bundles are formed by lamellae of identical elongations. In contrast, when $D_1 \to 1$, both $\mu _{\log \rho \,|\,n}/\log n$ and $\sigma ^2_{B}/\log n$ become infinite, while their ratio is $\sigma ^2_{\log \rho \,|\,n}/\mu _{\log \rho \,|\,n} = 2(2-D_1) \to 0.5$. This limit suggests that the aggregation of lamellae remains correlated to their average elongation, although a fixed amount of stretching variability arises in bundles. Good agreement is found between theoretical prediction (4.21)–(4.22) and numerical simulations of aggregation in the baker map (figure 13). In contrast, the theory captures only qualitatively the behaviour of the random sine flow. This may be due to the continuity of the sine flow, which produces curved lamellar structures that are not exactly one-dimensional.

Figure 13. Scaling of the (a) mean and (b) variance of log-elongation in bundles as a function of the information dimension $D_1$. Circles stands for numerical simulations in baker maps (open circles) and sine flow (filled circles). Solid lines stand for theoretical prediction of the mean (4.21) and variance (4.22). Dashed and dotted lines are plotted to compare the mean with fractal dimensions of other orders.

The stretching variability in bundles is directly linked to the heterogeneity of the chaotic flow, because of the intimate relationship existing between the fractal geometry of the material line and the stretching statistics of fluid elements (Ott & Antonsen 1989). As such, it is impossible to have a single stretching rate per bundle as soon as the chaotic flow is heterogeneous and exhibits a distribution of stretching rates. The absence of stretching variability in bundles ($\sigma ^2_{\log \rho \,|\,n}=0$) implies the uniformity of stretching at large scale ($\sigma ^2_{\rho }=0$). This uniform case is reached when $D_1\to D_0=2$, for instance, in the baker map when $a \to 0.5$. In continuous flow maps such as the sine flow, regions of high and low stretching always coexist, and $\sigma ^2_{\log \rho \,|\,n}>0$.

4.4. Moments of $1/\rho$ in a bundle of size $n$

Having described the first two moments of the distribution of lamella elongation in bundles (4.21)–(4.22), we now assume that the distribution is of log-normal shape. This choice is justified by the fact that elongation is a multiplicative process, thus usually leading to log-normal distributions (Le Borgne, Dentz & Villermaux 2015; Souzy et al. 2020). This allows us to compute the scaling of the $q$ moments of lamella concentrations in bundles, $\theta \,|\,n$. Owing to (2.14), we have

(4.23)\begin{align} \mu_{\theta^q \,|\, n} &\sim \mu_{ \rho^{-q} \,|\, n } \nonumber\\ &= \int_1^\infty \rho^{-q}\,P_{\rho\,|\,n}(\rho)\,\text{d}\rho \nonumber\\ &\approx \int_1^\infty \exp({-(\log\rho-\mu_{\log \rho\,|\,n})^2/ (2\sigma_{\log \rho\,|\,n}^2) - q \log \rho})\,\text{d}\rho. \end{align}

The minimum bound for the integral is taken at $\rho =1$ and not 0, taking into account the fact that lamellar structures cannot be compressed in their longitudinal direction. As a consequence, $P_{\rho \,|\,n}(\rho )$ is truncated for $\rho <1$. Denoting $\varLambda =\log \rho / \log n$, $\tilde {\mu }=\mu _{\log \rho \,|\,n}/\log n$ and $\tilde {\sigma }^2=\sigma _{\log \rho \,|\,n}^2/\log n$, this expression becomes

(4.24)\begin{equation} \mu_{\rho^{-q}\,|\,n} \approx \int_1^\infty \exp({H(\varLambda)\log n})\,\text{d}\rho, \end{equation}

with $H(\varLambda )=-(\varLambda -\tilde {\mu })^2/(2\tilde {\sigma }^2)-q\varLambda$. For large $n$, the value of this integral tends to $\exp ({H(\varLambda ^*)\log n})$, where $\varLambda ^*$ is the value where $H$ takes a maximum, that is, either at $\varLambda ^*=\tilde {\mu }-q\tilde {\sigma }^2$ if $\tilde {\mu }-q\tilde {\sigma }^2>0$, or at $\varLambda ^*=0$ otherwise. Thus

(4.25)\begin{equation} \log \mu_{\rho^{-q} \,|\, n} \approx-({\gamma}_{q,\rho^{-1}\,|\,n} \log n +{\omega}_{q,\rho^{-1}\,|\,n}), \end{equation}

with

(4.26) \begin{equation} \left.\begin{array}{@{}ll} {\gamma}_{q,\rho^{-1}\,|\,n}= q\tilde{\mu} - q^2 \tilde{\sigma}^2/2 & \text{if } \tilde{\mu} > q \tilde{\sigma}^2, \\[2pt] {\gamma}_{q,\rho^{-1}\,|\, n}= \tilde{\mu}^2/(2\tilde{\sigma}^2) & \text{if } \tilde{\mu} \leq q \tilde{\sigma}^2, \end{array}\right\} \end{equation}

and ${\omega }_{q,\rho ^{-1}\,|\,n}=q \log (\mathcal {A}/(s_a \ell _0 ))$ a constant. In particular, we are interested in the exponent $q=2$, which is useful to describe fluctuations around the mean. We have

(4.27)\begin{equation} \tilde{\gamma} \equiv {\gamma}_{2,\rho^{-1}\,|\,n} = \left\lbrace \begin{array}{@{}ll} 2 \tilde\mu_\lambda- 2 \tilde\sigma^2_\lambda & \text{if } \tilde{\mu} > 2 \tilde{\sigma}^2, \\[3pt] \tilde\mu^2/(2\tilde\sigma^2) & \text{if } \tilde{\mu} \leq 2 \tilde{\sigma}^2. \end{array} \right. \end{equation}

The predicted dependence of $\tilde {\gamma }$ upon $D_1$ is reproduced in figure 14. Here, $\tilde {\gamma }$ is bounded between 2 (for $D_1 \to 2$) and 1 for $D_1\approx 1.5$. The prediction agrees reasonably well with numerical simulations of the baker map and sine flow (figure 14). The small discrepancies can be attributed to deviations from log-normally distributed elongation in bundles. We also verify numerically that $\tilde \omega \equiv \omega _{2,\rho ^{-1}\,|\,n}$ is independent of $D_1$ (figure 14b). In figure 15, we verify that $\tilde \gamma$ is independent of the aggregation scale $s_a$. In contrast,

(4.28)\begin{equation} \tilde\omega \approx 2 \log (\mathcal{A}/(l_0 s_a)) \end{equation}

is a function of the aggregation scale only, independent of time and fractal dimension (figure 14b). Having determined both the elongation statistics inside aggregates of size $n$ and the spatial distribution of $n$, we will deduce in the following section the statistics of aggregated scalar levels $c$.

Figure 14. (a) Scaling exponents of the $q$ lamellar concentration moments in bundles (4.23). Numerical estimates are plotted with symbols: red diamonds for $q=2$, and black circles for $q=1$. Unfilled and filled symbols represent simulations in the baker map and sine flow, respectively. Theoretical predictions (4.27) are represented by lines. (b) Intercept $\tilde \omega$ of the scaling exponent of the $q$ lamellar concentration moments in bundles (4.23) in the baker map (empty squares) and sine flow (filled squares), and theoretical prediction (line, (4.28)).

Figure 15. Dependence of $\tilde \gamma$ and $\tilde \omega$ with the aggregation scale in simulations of the baker map (empty symbols, $a=0.2$) and the sine flow (filled symbols, $A=1.2$), and comparison to theoretical prediction (4.28).

5. Aggregated scalar concentrations

In the preceding section, we have shown that the moments of lamellar compression inside a bundle of size $n$ follow

(5.1)$$\begin{gather} \mu_{\rho^{-1} \,|\, n} = \frac{s_a \ell_0}{\mathcal{A}}\,n^{-1}, \end{gather}$$

(5.2)$$\begin{gather}\mu_{\rho^{-2} \,|\, n} = \frac{(s_a \ell_0)^2}{\mathcal{A}^2}\,n^{-\tilde{\gamma}}, \end{gather}$$

with $\tilde {\gamma }\equiv \gamma _{2,\rho ^{-1}\,|\,n}$ a flow-dependent exponent varying with the fractal dimension $D_1$ and taking value between 1 and 2 (figure 14a). The variance of lamellar concentrations inside bundles thus follows:

(5.3)\begin{equation} \sigma^2_{\rho^{-1}\,|\,n} = \frac{(s_a \ell_0)^2}{\mathcal{A}^2}\, (n^{-\tilde{\gamma}}- n^{-2}).\end{equation}

We first assume that the aggregated scalar concentration of a bundle $c\,|\,n$ can be obtained by a sum of independent and identically distributed random variables (2.21), whose statistics have been described previously. In doing so, we make two hypotheses. First, we assume that the elongation statistics of lamellae inside a bundle of size $n$ or among all bundles of size $n$ are comparable, which is true for large $n$ (see Appendix C). Second, we assume that variability in bundle aggregated concentrations arises from independent realizations of the sum. In other words, we assume that the chaotic flow is sufficiently random to shuffle lamellar elongations between each bundles. We will show later that independence is not ensured in the deterministic baker map.

With these hypotheses, the mean aggregated concentration is

(5.4)\begin{equation} \mu_{c\,|\,n} = n \times \frac{\sqrt{\rm \pi}\,\theta_0 s_0}{s_a}\,\mu_{\rho^{-1}\,|\,n} = \frac{\sqrt{\rm \pi}\,\theta_0 \ell_0 s_0}{\mathcal{A}}, \end{equation}

which is also the mean concentration $\mu _c$. The variance of aggregated concentrations over bundles of similar size is

(5.5)\begin{equation} \sigma^2_{c\,|\,n} \approx n \left(\frac{\sqrt{\rm \pi}\,\theta_0 s_0}{s_a}\right)^2 \sigma^2_{\rho^{-1}\,|\,n} =\frac{(\sqrt{\rm \pi}\,\theta_0 \ell_0 s_0)^2}{\mathcal{A}^2}\, (n^{1-\tilde{\gamma}} - n^{-1}). \end{equation}

When $n$ is large, this expression further simplifies to

(5.6)\begin{equation} \sigma^2_{c\,|\,n} \sim n^{1-\tilde{\gamma}}, \end{equation}

with $\tilde {\gamma } \in [1,2]$ given by (4.27).

The scaling obtained with the independent realization hypothesis compares well with numerical simulations of aggregation in random sine flows of various heterogeneity (figure 16). However, it largely underestimates the scaling observed in the deterministic baker map. Indeed, the simplicity and regularity of the deterministic baker map make bundles of similar size not statistically independent. While bundle concentration still results from the addition of variable lamellar concentrations, independent realizations of the sum are not achieved due to the deterministic nature of the map, the same combinations being repeated identically in most bundles, as in the case of a unique realization. In that case, the variance of the sum is the variance of the random variable, and (5.5) transforms into

(5.7)\begin{equation} \sigma^2_{c\,|\,n} \sim \sigma^2_{\rho^{-1},n} \sim n^{-\tilde{\gamma}}. \end{equation}

This scaling indeed fits more accurately the deterministic baker map simulations (figure 16).

Figure 16. Scaling exponent $\xi$ (5.8) of the variance of bundle concentrations knowing $n$ estimated from simulations (dots) and theoretical predictions with the independent realization hypothesis for the sine flow (solid lines, $\xi =\tilde \gamma - 1$, (5.6)) and baker map (dashed lines, $\xi =\tilde \gamma$, (5.7)).

To summarize, the addition of lamellar concentration levels in a bundle yields a concentration whose deviation from the mean decays algebraically with the number of lamellae in the bundle:

(5.8)\begin{equation} \sigma^2_{c \,|\, n} \approx \frac{(\theta_0 \ell_0 s_0)^2}{\mathcal{A}^2}\,n^{-\xi}, \end{equation}

with $\xi =\tilde {\gamma }$ for purely deterministic flows (baker map), and $\xi =\tilde {\gamma }-1$ for random flows (sine flow). We call $\xi$ the correlation exponent, which can take values between 0 and 2 depending on the flow heterogeneity and randomness.

5.1. Distribution of $c$

We now focus on deriving the expression of the p.d.f. of aggregated scalar concentration $P_c$ over the flow domain. The later can be obtained from the conditional p.d.f. of bundle aggregated concentrations knowing their size $P_{c\,|\,n}$ and the p.d.f. of bundle size $P_n$ by the weighted summation

(5.9)\begin{equation} P_c(c) = \int_n \text{d}n\,P_{c\,|\,n}(c)\,P_n(n). \end{equation}

Here, $P_{n}$ was found to be well approximated by a Gamma distribution, with moments given by the fractal characteristics of advected material lines (4.14). In turn, we have determined the scaling of the first two moments of bundle aggregated concentration $c\,|\,n$ (5.8), without specifying the precise shape of the p.d.f. A natural choice for $P_{c\,|\,n}$ is the log-normal distribution, since the sum of log-normally distributed random variables is known (Schwartz & Yeh 1982) to be well approximated by log-normal distributions. Adopting the log-normal shape, the parameters of the distribution are

(5.10)$$\begin{gather} \mu_{\log c\,|\,n} = \log\mu_{c\,|\,n} - \log(\sigma^2_{c\,|\,n}/\mu_{c\,|\,n}^2+1)/2, \end{gather}$$

(5.11)$$\begin{gather}\sigma^2_{\log c\,|\,n} = \log(\sigma^2_{c\,|\,n}/\mu_{c\,|\,n}^2+1). \end{gather}$$

In figure 17, we plot the simulated distribution of aggregated concentration levels compared to the prediction (5.9) for the baker map and sine flow. The agreement is fair in the region near $\mu _c$, but deviates for large $c$. Indeed, this corresponds to lamellae with weak aggregation for which $n\approx 1$. In this region, the solitary strip p.d.f. $P_{\rho ^{-1}}$ describes well the tail of $P_c$ because such high concentration excursions are essentially supported by isolated lamellae, while the correlated aggregation model assumes $n\gg 1$. The presence of these weakly aggregated, high concentration levels is particularly evident at high Péclet number (figure 17b). The scalar concentration p.d.f. is thus the combination of an aggregated core around the mean following (5.9) and tails following the isolated strip concentration p.d.f.

Figure 17. Distributions of aggregated scalar concentrations in the sine flow depending on (a) the sine wave amplitude $A$ ($s_a=1/50$), and (b) the aggregation scale $s_a$ ($A=0.9$). Symbols stand for numerical simulations, solid lines are the aggregation model (5.9), and dashed lines are the isolated strip prediction (2.6). Simulations are all taken at the time when the total filament length reaches $L=10^7 \ell _0$.

In figure 18, we compare the correlated aggregation model with the random aggregation model where $\mu _n \sim L(t)$ (2.19). The random aggregation assumption yields gamma p.d.f.s that are narrowing much faster than the simulated p.d.f.s in the sine flow. In contrast, the correlated model better captures the p.d.f.

Figure 18. Plots of $P_c$ in sine flows at several times ($A=0.8$, $s_a=1/50$) (symbols) compared with the random aggregation (dashed lines) and correlated aggregation (solid lines, (5.9)) models.

From the p.d.f. of aggregated scalar concentration, we now derive its moments. They are directly related to the p.d.f. of $n$, since

(5.12)\begin{equation} \mu_c = \int \text{d}c \sum_n c\,P(c\,|\,n)\,P(n) = \sum_n \mu_{c\,|\,n}\,P(n) = \frac{\theta_0 \ell_0 s_0}{\mathcal{A}} \end{equation}

and

(5.13)\begin{equation} \mu_{c^2} = \int_c \text{d}c \int_n \text{d}n\,c^2\,P(c\,|\,n)\,P(n) = \int_n \text{d}n\,\mu_{c^2\,|\,n}\,P(n) = \frac{(\theta_0 \ell_0 s_0)^2}{\mathcal{A}^2}\, (\mu_{n^{-\xi}} +1). \end{equation}

Thus the scalar variance is

(5.14)\begin{equation} \sigma^2_c = \frac{(\theta_0 \ell_0 s_0)^2}{\mathcal{A}^2}\,\mu_{n^{-\xi}}. \end{equation}

Note that because we chose a gamma distribution for $n$ with parameters $\theta _n$ and $k_n$ defined in (4.14), $\mu _{n^{-\xi }}$ is not defined for all $D_2$. Indeed, the gamma distribution admits a power law tail at small $n$ with exponent $k_n-1$, which affects the value of negative moments: when $k_n(D_2)<\xi (D_2)$, $\mu _n^{-\xi }$ does not exist. However, because fully chaotic flows are space-filling ($D_0=2$), we must have $n\geq 1$ at late time. Thus we restrict the p.d.f. of $n$ to $n\geq 1$ and get

(5.15)\begin{equation} \mu_{n^{-\xi}} \sim (\theta_n)^{-\min(k_n,\xi)}. \end{equation}

An intuitive understanding of this equation can be formulated as follows. If the spatial heterogeneity of $n$ is moderate ($\xi < k_n$), then the average of $n^{-\xi }$ is affected by all values of $n$ in the distribution. In contrast, if the heterogeneity is stronger ($\xi >k_n$), then the probability of having low aggregation regions ($n\approx 1$) is high and controls the value of $\mu _{n^{-\xi }}$. In that case, the average no longer scales with $\xi$, but rather scales with the parameter $k_n$, explaining the minimum in the exponent. The value of the exponent depending on Péclet number and fractal dimension is plotted in figure 19 for the sine flow. In such random flow, the exponent takes values between 0 and 1. For high Péclet numbers and low fractal dimension, $\mu _{n^{-\xi }}$ is governed by weakly aggregated regions and $\min (k_n,\xi )=k_n$. In contrast, for low Péclet number and large fractal dimension, the whole distribution of $n$ plays a role in the determination of $\mu _{n^{-\xi }}$, and $\min (k_n,\xi )=\xi$.

Figure 19. Evolution of the exponent $\min (k_n,\xi )$ (contours) in the sine flow for various Péclet numbers (aggregation scale $s_A$) and sine amplitude $A$ (fractal dimension $D_2$). Here, $k_n$ is obtained with (4.13) and $\xi =\tilde {\gamma }-1$ with (4.27).

Combining (5.15) and (4.14) provides the asymptotic scalar variance decay as a function of the domain area $\mathcal {A}$, the elongation of material lines $L(t)$, and the aggregation scale $s_a=\sqrt {2{\rm \pi} }\,s_B$:

(5.16)\begin{equation} \sigma^2_c(t) = \left(\frac{L(t)\,s_a (\sqrt{\mathcal{A}}\, s_a^{D_2-2}-1)}{\mathcal{A}}\right)^{-\min(k_n,\xi)} \sim L(t)^{-\min(k_n,\xi)}.\end{equation}

Thus in a correlated aggregation scenario, the decay exponent of scalar variance is found to be a fraction of the growth exponent of material lines, respectively $\log 2$ and $\mu _\lambda +\sigma ^2_\lambda /2$ in the baker map and sine flow.

Equation (5.16) suggests that the existence of correlations in the aggregation process may be envisioned as a retardation of the purely random aggregation scenario (2.22). From this perspective, the effective number of random aggregation events $L(t)^{\min (k_n,\xi )}$ is smaller than the total number of aggregation $L(t)$, because the correlation exponent $\min (k_n,\xi )$ is generally smaller than 1 (see figure 19).

In contrast, the random aggregation model (2.22) clearly overestimates the variance decay rate in the sine flow. This is explained by the correlated nature of aggregation, which is less efficient at homogenizing concentration levels than a completely random addition. In other words, small concentration levels have a higher probability of coalescing with other small concentrations than with high concentrations, retarding the homogenization of the mixture. Indeed, the correlated aggregation model proposed herein (5.16) shows decay rates closer to the numerical observations. The solitary lamella model still outperforms the correlated aggregation framework in predicting scalar decay rate, because the latter is dominated by weakly aggregated regions of the flow. These regions are poorly described by statistics obtained in the large-$n$ limit (e.g. (5.5)).

A different picture is observed for scalar decay rates obtained in the baker map (figure 20b), because of the deterministic nature of the map. The solitary lamella model is accurate for fractal dimensions $D_1<1.6$, where stretching heterogeneity is important and weakly aggregated regions of the flow dominate scalar fluctuations. For $D_1>1.6$, the isolated lamella model underpredicts scalar decay due to aggregation effects, as explained in the following. The random aggregation scenario predicts an invariant scalar decay rate ($\gamma _{2,c} = \log 2$) equal to the growth rate of material lines. It thus overpredicts the observed decay rates for $D_1<1.8$, and underpredicts them at larger values of $D_1$. Interestingly, when the flow tends to the uniform case $a\to 0.5$ ($D_1\to 2$), simulations yield scalar decay rate $2\log 2$, larger than the decay rate of isolated lamellae. This acceleration of mixing by aggregation is a pure consequence of the determinism nature of the map, and is well captured by the correlated aggregation scenario (5.7). However, for all baker maps with $D_1<1.9$ ($a<0.2$), $k_n<\xi$ and scalar fluctuations are mainly governed by the regions where $n\sim 1$ for which the correlated aggregation theory is not expected to be accurate. This explains the limited range of validity of the correlated aggregation theory for describing scalar decay rates in the deterministic baker map.

Figure 20. Decay exponent $\gamma _2$ of the variance of aggregated scalar levels with time, as a function of fractal dimension $D_1$ for (a) the baker map and (b) the random sine flow. Dots stand for numerical simulations, and lines from theoretical predictions for isolated lamellae ((2.17) and exponents of $\mu _{\rho ^{-1}}$ in tables 1 and 2), random aggregation ((2.22) and exponents of $1/\mu _{\rho }$ in tables 1 and 2) and correlated aggregation (5.14).

Table 1. Moments of log-normally distributed stretching sampled over infinitesimal fluid elements ($\mu_{\bullet}$), material line ($\mu_{\bullet,L}$).

Table 2. Moments of binomial distributed stretching sampled over infinitesimal fluid elements ($\mu_{\bullet}$), material line ($\mu_{\bullet,L}$).

Finally, note that the simulations do not show the super-exponential decay of scalar fluctuations classically observed for the uniform stretching rate at $a= 0.5$. In fact, the reconstruction of the scalar field by a summation of lamellar maximum concentrations on a fixed grid (3.6) impedes the apparition of the super-exponential mode. As $a\to 0.5$, all lamellae are subjected to similar stretching rates around $\log 2$, thus yielding a scalar variance decaying as $2\log 2$.

6. Conclusions

Scalar mixing by incompressible smooth two-dimensional chaotic flows leads to the emergence of two coupled phenomena that are the elongation of interfaces by stretching, and their aggregation and coalescence at the Batchelor scale by compression. In turbulent scalar mixing, aggregation was well captured by random addition rules (Duplat & Villermaux 2008). Here, we show that for smooth chaotic flows in the Batchelor range, aggregation is highly correlated. This correlated aggregation process significantly reduces the flow mixing efficiency compared to a random hypothesis, maintaining it close to the mixing efficiency obtained for solitary lamellae, and explaining the observed monotonic exponential decay of scalar variance before and after coalescence time (Fereday et al. 2002).

Using two-dimensional chaotic flows as a reference, we measured the aggregation rate of exponentially stretched material lines across a broad range of chaotic flow regimes. We showed that the most elongated lamellae are also the most aggregated ones, due to the fact that larger compression rates attract larger flow regions. The link between elongation and compression, induced by incompressibility, hence generates a direct correlation between elongation and aggregation. The heterogeneity in stretching rates therefore controls the heterogeneity of the number of lamellae in bundles.

We showed that the statistics of aggregated lamella numbers can be predicted exactly from the fractal dimensions of the elongated material line. We then derived a general theoretical framework that captures the effect of correlated aggregation, where lamellae of similar stretching aggregate preferentially, and predict the p.d.f.s of aggregated scalar levels. We find that correlated aggregation can be uniquely characterized by single correlation exponent $\xi \in [0,1]$, which provides a measure of the effective number of aggregation events, compared to the total number. In that sense, correlated aggregation delays the route to uniformity compared to a fully random hypothesis, although it does not alter the fundamental nature of the aggregation process (Villermaux & Duplat 2003).

Our results apply for two-dimensional fully chaotic flows in the Batchelor regime, that is, for smooth velocity fields below the integral scale. These flow fields are representative of a large class of flows, including notably porous media flows (Heyman et al. 2020; Souzy et al. 2020), geophysical flows, and turbulent flows at high Schmidt number. We have also considered small periodic flow domains where no large-scale scalar gradients can appear. When scalar length scales can develop beyond the integral velocity scale in smooth flows, mixing is controlled globally by the slowest dispersing modes (Haynes & Vanneste 2005; Tsang et al. 2005). The aggregation rules determined herein are not supposed to change above the velocity scale. Thus it is possible in principle to obtain the statistics of larger-scale scalar fields by the summation of individual lamellae. To do so, one may consider the mean density of aggregation $\mu _n(t)$ (4.4) to vary spatially through macro-dispersion, while microscale fluctuations would remain governed by (4.11). Note that similar ideas were used to predict the evolution of a dispersing scalar plume in two-dimensional porous media (Le Borgne et al. 2015).

It should also be possible to extend the correlated aggregation theory to three-dimensional flows that produce two-dimensional sheets (Ngan & Vanneste 2011; Martínez-Ruiz et al. 2018; Meunier & Villermaux 2022) that aggregate at late times due to folding. These flows have a unique negative Lyapunov exponent, the mean compression rate normal to the sheets, and thus also simplify to the one-dimensional Lagrangian stretching framework on which aggregation models are framed.