T L E B O R G N E

E V I L L E R M A U X

The lamellar description of mixing in porous media

We develop a general framework for modeling mixing in porous media flows, in which the scalar mixture is represented as an ensemble of lamellae evolving through stretching, diffusion and coalescence. Detailed numerical simulations in Darcy scale heterogeneous permeability fields are used to analyze the lamella deformation process, which controls the local concentration gradients and thus the evolution of the concentration mixture through stretching enhanced diffusion. The corresponding Lagrangian deformation process is shown to be well modeled by a Langevin equation with multiplicative noise, which can be coupled with diffusion to predict the temporal evolution of the concentration probability density function (PDF). At late times, lamella interaction is enforced by confinement of the mixture within the dispersion area. This process is shown to be well represented by a random aggregation model, which quantifies the frequency of lamella coalescence and allows predicting the temporal evolution of the concentration PDF in this regime. The proposed theoretical framework provides an accurate prediction of the concentration PDFs at all investigated times, heterogeneity levels and Péclet numbers. In particular, it relates the temporal behavior of mixing, as quantified by concentration moments, scalar dissipation rate or spatial increments of concentration, to the degree of structural heterogeneity.

Introduction

The prediction of the concentration distribution of scalar mixtures in porous media is a key element for the quantification of mixing and reactive transport processes (e.g. [START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF], and the assessment of uncertainty and risk for contaminant transport in the subsurface (e.g. [START_REF] Tartakovsky | Assessment and management of risk in subsurface hydrology: A review and perspective[END_REF]). The development of predictive models for these purposes requires upscaling the impact of flow heterogeneity on effective transport. To this end, one would ideally wish to rely on solid guidelines to describe the entire concentration field in a given medium, for given injection conditions. In the context of porous media, the discussion has historically been mostly confined around a particular aspect of the overall problem, namely the dispersive properties of the medium. By 'dispersive', we mean the description of the extent of the spatial spreading of the mixture, or of its residence time in the vicinity of a given point in the medium (see, e.g, [START_REF] Danckwerts | Continuous flow systems[END_REF][START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] De Gennes | Hydrodynamic dispersion in unsaturated porous media[END_REF][START_REF] Gelhar | Three-dimensional stochastic analysis of macrodispersion in aquifers[END_REF][START_REF] Bouchaud | A simple model for hydrodynamic dispersion[END_REF][START_REF] Koch | The effect of order on dispersion in porous media[END_REF][START_REF] Cushman | Nonequilibrium statistical mechanics of preasymptotic dispersion[END_REF][START_REF] Dagan | Theory of solute transport by groundwater[END_REF][START_REF] Attinger | Temporal behavior of a solute cloud in a chemically heterogeneous porous medium[END_REF][START_REF] Dentz | Temporal behavior of a solute cloud in a heterogeneous porous medium, 1, point-like injection[END_REF][START_REF] Delgado | Longitudinal and transverse dispersion in porous media[END_REF]. Field measurements and laboratory experiments have thus prompted research to understand the scale dependence of global dispersion coefficients and its relation to the medium heterogeneity, the influence of the possible presence of traps and dead zones for the transported substances, as well as the validity of the advection-dispersion equation to describe large scale solute transport (e.g. [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Neuman | Perspective on theories of anomalous transport in heterogeneous media[END_REF]. However, the spatial extent of a mixture does not tell much about its content, but rather its average concentration, and the knowledge of the spreading properties alone is insufficient to quantify mixing (e.g., [START_REF] Kitanidis | The concept of the dilution index[END_REF][START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF][START_REF] Le Borgne | Persistence of incomplete mixing: A key to anomalous transport[END_REF], that is the full concentration distribution, and the local spatial structure of the concentration field in the medium.

Mixing processes in porous media are particularly important for controlling fluid-fluid as well as fluid-solid chemical reactions (e.g. [START_REF] Cirpka | Two-dimenstional concentration distribution for mixing-controlled bioreactive transport in steady-state[END_REF][START_REF] Cirpka | Concentration statistics for mixing-controlled reactive transport in random heterogeneous media[END_REF][START_REF] Rolle | Enhancement of dilution and transverse reactive mixing in porous media: experiments and modelbased interpretation[END_REF][START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF]de Anna et al. 2014). Dissolution or precipitation processes, for instance, are often triggered by the mixing of different waters respectively at chemical equilibrium (e.g. [START_REF] Tartakovsky | Mixing-induced precipitation: Experimental study and multiscale numerical analysis[END_REF]. Different mixing measures have been studied to investigate these processes in porous media [START_REF] Kitanidis | The concept of the dilution index[END_REF][START_REF] Fiori | Concentration fluctuations in aquifer transport: a rigorous firstorder solution and applications[END_REF][START_REF] Fiori | The lagrangian concentration approach for determining dilution in aquifer transport: Theoretical analysis and comparison with field experiments[END_REF]Bolster et al. submitted;[START_REF] Le Borgne | Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF][START_REF] Chiogna | Mixing, entropy and reactive solute transport[END_REF][START_REF] De Barros | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF]. The dynamics of mixing can be fully quantified by the temporal evolution of the concentration probability density function (PDF), as the mixture evolves towards uniformity under the combined action of diffusion and advection. A common approach used is to assume an ad-hoc distribution, which is parametrized from analytical derivation of concentration moments obtained under different approximations (e.g., [START_REF] Fiorotto | Concentration statistics in heterogeneous aquifers for finite peclet values[END_REF][START_REF] Bellin | Probability density function of non-reactive solute concentration in heterogeneous porous formations[END_REF][START_REF] Tonina | Effects of pore-scale dispersion, degree of heterogeneity, sampling size, and source volume on the concentration moments of conservative solutes in heterogeneous formations[END_REF][START_REF] Cirpka | Concentration statistics for mixing-controlled reactive transport in random heterogeneous media[END_REF]. Other approaches include mapping closures to transfer the stochasticity of the flow field to the concentration distribution (e.g., [START_REF] Chen | Probability distribution of a stochastically advected scalar field[END_REF][START_REF] Dentz | Probability density functions for passive scalars dispersed in random velocity fields[END_REF][START_REF] Dentz | Concentration statistics for transport in heterogeneous media due to stochastic fluctuations of the center of mass velocity[END_REF]) and phenomenological approaches that employ stochastic mixing models [START_REF] Valiño | A binomial Langevin model for turbulent mixing[END_REF][START_REF] Fedotov | Stochastic mixing model with power law decay of variance[END_REF] to simulate the stochastic evolution of concentration using surrogate mechanisms, see also [START_REF] Pope | Turbulent Flows[END_REF] and [START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF] for an overview.

The study of scalar mixing in a simple stirring protocol by [START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF] has nevertheless demonstrated that these mechanisms can be quantified from first principles using the powerful concept of the 'diffusive strip', or 'lamellae' representation of the scalar field [START_REF] Ranz | Application of a stretch model to mixing, diffusion and reaction in laminar and turbulent flows[END_REF]. Scalar lamellae are naturally formed by the repeated action of advective motions. In this representation the coupling of diffusion to the stretching experienced by the lamellae is straightforward. The complete PDF of concentration can thus be computed from the diffusion equation in a moving frame, and an additional principle of lamellae random coalescence initially proposed by [START_REF] Villermaux | Mixing as an aggregation process[END_REF] allows for the description of the mixture content up to its final uniformity. In the same line of thought, we have recently highlighted the role of stretching, diffusion and coalescence mechanisms in governing the evolution of concentration mixtures in Darcy scale porous media [START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Le Borgne | Stretching, coalescence and mixing in porous media[END_REF]. While the porous media structure does not deform in general, the 'frozen' flow heterogeneity at different scales generates similar folding and stretching mechanisms as in chaotic or turbulent flows [START_REF] Metcalfe | A partially open porous media flow with chaotic advection: towards a model of coupled fields[END_REF][START_REF] De Barros | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF][START_REF] Mays | Plume spreading in groundwater by stretching and folding[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF]. In this framework, the understanding of the basic composition rule of the dispersing mixture provides an analytical description of the overall concentration distribution.

We explore here the potential of this theory for predicting mixing in Darcy scale heterogeneous porous medium over a range of Péclet numbers and structural heterogeneities. The main challenge when dealing with these highly heterogeneous flow systems is that dispersion and mixing occur simultaneously. This leads to a range of mixing behaviors at the ends of which lie two turbulent flow equivalents: the ever dispersing mixture scenario [START_REF] Duplat | A nonsequential turbulent mixing process[END_REF] for large Péclet numbers, and the confined mixture scenario [START_REF] Duplat | Mixing by random stirring in confined mixtures[END_REF] for small Péclet numbers. In the first case, the concentration content can be determined from the stretching history of individual lamellae alone, while for the second case, the concentration content is given by the frequency of coalescence of neighboring lamellae. The large flow heterogeneity induced by the underlying permeability structures implies here the existence of an extended intermediate regime, characterized by partial lamellae coalescence and the persistence of concentration lacunarities in the mixture. This regime is also characterized by pre-asymptotic super-diffusive dispersion, which is known to persist over a large range of temporal and spatial scales in heterogeneous porous media (e.g. [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Neuman | Perspective on theories of anomalous transport in heterogeneous media[END_REF]. As discussed in this study, super-diffusive spreading has a direct impact on the stretching behavior and on the evolution of concentrations within the scalar mixture.

We investigate these processes in the light of a lamellar mixing model using numerical simulations of the full Darcy scale flow and transport problem, which allows for the detailed measurement of the lamellae deformation and concentration contents of the mixture. By tagging a line of passive scalar perpendicular to the mean flow in a well characterized heterogeneous permeability field, we analyze in details the Lagrangian stretching distribution and correlations (sections 2 to 5), and the evolution of concentration PDFs in relation to the geometry of material lines and their coalescence frequency in sections 6 and 7. Furthermore, we discuss the impact of the medium heterogeneity on the mixing behavior, as quantified by different mixing measures in sections 8 and 9, as well as on the spatial structure of the concentration fields through the analysis of the PDF of concentration increments in section 10.

Flow and transport in heterogeneous permeability fields

In this section, we introduce the basic flow and transport equations and their boundary and initial conditions, as well as their numerical solution. We consider two dimensional heterogeneous conductivity fields, where the conductivity K(x) is represented as a lognormally distributed random field. The joint distribution of the point values of log-conductivity f (x) = ln K(x) is multivariate Gaussian and characterized by a mean f variance σ 2 f and correlation length λ. The correlation function of the fluctuations of f (x) is taken as Gaussian C f (x) = exp[-x 2 /(2λ 2 )]. This type of random conductivity is generic as it allows representing a range of spatial variability and correlation properties observed in natural permeability fields (e.g. [START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF]. As such, it represents a reference field for theories of flow and transport in heterogeneous porous media (e.g. [START_REF] Fiori | Concentration fluctuations in transport by groundwater: Comparison between theory and field experiments[END_REF][START_REF] Cirpka | Two-dimenstional concentration distribution for mixing-controlled bioreactive transport in steady-state[END_REF][START_REF] Le Borgne | Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF]. Since porosity is known to vary over a much narrower range than permeability (e.g. [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF], the impact of porosity heterogeneity on solute transport is generally considered subdominant compared to spatial variability in hydraulic conductivity. Thus, we consider here a constant porosity in order to restrain the space of parameters to be explored. The methodology may be however extended to consider porosity variations as well.

Spatial variability of the conductivity K(x) is mapped onto the divergence-free flow field u(x) via the Darcy equation u(x) = -K(x)∇h(x), with h(x) the hydraulic head. We assume that fluid and porous medium are both incompressible so that continuity is expressed by ∇ • u(x) = 0. The flow is driven by a uniform in the mean head gradient due to constant head boundary conditions at the lateral domain boundaries. At the horizontal boundaries no-flow conditions are applied. The mean head gradient and the geometric mean conductivity are set to one such that the average flow velocity is unity. Note that all spatial dimensions are normalized to the permeability correlation length and all temporal scales to the advection time so that results may be applicable to corresponding field cases upon appropriate rescaling. We consider here non-reactive transport through the divergence-free Darcy-flow field u(x) described above. The transport of concentration c(x, t) is described by the Fokker-Planck equation

∂c(x, t) ∂t + u(x) • ∇c(x, t) -D∇ 2 c(x, t) = 0. (2.1)
In general, dispersion is a tensorial quantity and depends on the Darcy velocity u(x).

For the sake of simplicity, and to to focus on the impact of advective heterogeneity, it dispersion is represented here the constant dispersion coefficient D. Time is rescaled with constant porosity φ such that t = t /φ, where t is the simulation time. As boundary conditions for concentration, we impose n v • D∇c(x, t) = 0 at the vertical boundaries, and zero normal flux, that is n h • [u(x) -D∇]c(x, t) = 0 at the horizontal boundaries, where n h and n h are the outwards pointing unit normal vector at the vertical and horizontal boundaries. The initial condition is given by an instantaneous uniform line injection at

x 1 = λ such that c(x, 0) = L -1 δ(x 1 -λ)
with L the extension of the domain perpendicular to the mean flow direction. Notice that we choose an initial plume that integrates a representative part of the heterogeneity spectrum in order to systematically quantify the impact of stochastic flow fluctuations on the mixing behavior. In the opposite case of point-like solute injection, the initial mixing behavior is expected to be different [START_REF] De Barros | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF].

The characteristic advection time scale is defined by τ a = λ/u, where u is the mean transport velocity. The characteristic diffusion time is τ D = λ 2 /D. The Péclet number

P e = λu D (2.2)
compares the diffusive and advective time scales. The transport equation (2.1) is solved by random walk particle tracking using N p = 10 7 particles based on the equivalent Langevin equation

dx(t|a) dt = u[x(t|a)] + √ 2Dξ(t), (2.3)
with ξ(t) a Gaussian white noise of zero mean and unit variance, and a is the initial tracer injection location, x(t = 0|a) = a. The flow equation ∆h(x)+∇f (x)•∇h(x) = 0 is solved numerically with a finite difference scheme. The flow domain is the two-dimensional rectangle Ω = {x|x 1 ∈ [0, 819.2λ], x 2 ∈ [0, 102.4λ]}. Both flow and transport are solved numerically using the numerical modeling platform H2OLAB (http://h2olab.inria.fr).

The probability density function (PDF) p c (c, t) of concentration point values is obtained from spatially sampling all concentration values in the domain that are above a concentration threshold . This defines the characteristic volume V c (t) that is occupied by the solute

V c (t) = c(x,t)> dx.
(2.4) Figure 1a displays an example of the simulated velocity field for a moderately heterogeneous permeability field, with σ 2 lnK = 1. The corresponding concentration fields at t 1 = 12.5τ a and t 2 = 112.5τ a are shown in figure 1b and figure 1c, for an initial injection with a uniform concentration c 0 at longitudinal position x 0 = λ, for different Péclet numbers. The scalar line, initially located perpendicular to the main flow direction, distorts into a brush with a typical streamwise width σ, increasing in time. The brush is made of a collection of strips, aligned on average with the direction of the mean flow. Similar behaviors can be observed for solute plumes originating from different extended initial distributions. The deformation action of the flow field on the fluid support leads eventually to the creation of laminary structures such as the ones observed here for a uniform initial concentration along a line.

For large Péclet numbers, the scalar concentration is close to uniform along the stretching direction, and varies rapidly along the compressive one, thus forming a lamella-like topology [START_REF] Ottino | The Kinematics of Mixing: Stretching, Chaos, and Transport[END_REF]. Hence, the concentration field can be seen as a set of lamellae, whose concentration depends on the interaction between local stretching and diffusion. Compression perpendicular to the stretching direction enhances diffusive smoothing by steepening the scalar gradient. For small Péclet numbers, neighbouring lamellae overlap by diffusive coalescence. The concentration field is then composed of lamella aggregates and concentration lacunarities in areas with small lamella density. The permeability field heterogeneity has a strong impact on the concentration distribution as it determines the distribution of stretching rates along the mixing front and the frequency of lamella coalescence, see Figure 2. In the following, we formalize the role to these different processes in building up the overall concentration distribution.

Local interaction between stretching and diffusion

We recall here the basic interaction between stretching and diffusion at the scale of a lamella, leading to an expression for the concentration PDF across a strip that depends on the stretching history [START_REF] Ranz | Application of a stretch model to mixing, diffusion and reaction in laminar and turbulent flows[END_REF][START_REF] Duplat | A nonsequential turbulent mixing process[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF]. Considering a single lamella, the temporal evolution of the concentration c(ζ, t) at a position ζ along the direction perpendicular to the lamella is

∂c ∂t + ζ s ds dt ∂c ∂ζ = D ∂ 2 c ∂ζ 2 . (3.1)
where D is the diffusion coefficient and s(t) is the distance between two material particles in the direction perpendicular to the lamella. By mass conservation under pure advection, s(t) is related to the lamella length (t) by s(t) (t) = s 0 0 , where s 0 and 0 are the initial lamella dimensions. Therefore s -1 ds/dt = --1 d /dt. For the initial condition c(ζ, t = 0) = c 0 exp(-ζ 2 ), the solution is given by the Gaussian profile,

c(ζ, t) = c 0 exp -ζ 2 ρ(t) 2 s 2 0 (1+4τ (t)) 1 + 4τ (t) , (3.2)
where we defined the dimensionless time

τ (t) = D s 2 0 t 0 dt ρ(t ) 2 (3.3) from the relative strip elongation ρ(t) = (t)/ 0 . The maximum lamella concentration θ in ζ = 0 is θ(t) = c 0 1 + 4τ (t) (3.4)
The concentration PDF in the strip corresponding to the Gaussian profile (3.2) is given by

p(c|θ) = 1 2c ln(θ/ ) ln(θ/c) (3.5)
with θ c. The concentration range is [ , θ], with and θ the minimum and maximum concentrations in the strip. Equation (3.5) provides the relationship between the concentration PDF of a lamella and its elongation history, through the temporal evolution of the maximum concentration (3.4) in the strip.

Statistics of Lagrangian elongations

In this section we analyze the stretching of lamellae in the heterogeneous flow field using the explicit numerical simulations. The stretching process is then represented by a multiplicative stochastic process that models the series of Lagrangian deformation events, and thus the stretching history.

The history of relative lamella elongation ρ(t) can be quantified by computing the pair separation (t) = |x(t|a + ∆a) -x(t|a)| of neighbouring particles initially located at a distance 0 = |∆a| on the injection line. Figure 3 displays the location of 5 × 10 5 advective particles at different times, initially separated by a distance 0 = 1.25 × 10 -4 λ. An example of Lagrangian elongation history for a particle pair is given in Figure 4. The stretching rate γ(t) in the following is defined by

γ(t) = 1 ρ(t) dρ(t) dt . (4.1)
Strong stretching events are typically followed by compression events, with a small effect on the net elongation (Figure 4). Significant variations of the elongation occur during localized events characterized by a succession of positive elongation rates. In between these events, the total elongation rate remains approximately constant. The dimensionless time τ (t), which controls the maximum lamella concentration (3.3), also changes significantly during the positive stretching events. In between these events, it evolves linearly with time as the elongation is approximately constant. Successive stretching rates appear to be correlated over relatively short times. The correlation function of the fluctuations of γ(t), namely γ (t) = γ(t) -γ(t) is defined by

C γ (τ ) = 1 σ 2 γ (t) γ (t)γ (t + τ ) , σ 2 γ (t) = γ (t) 2 , (4.2)
where the angular brackets denote here the average over all particle pairs. As illustrated in Figure 5 Figure 6. (a) PDF of elongations at times t1 = 0.25τa , t2 = 12.5τa, t3 = 50τa and t4 = 112.5τa, (b) PDF of elongations at time t4 = 112.5τa for σ 2 f = 0.25 (green diamonds), σ 2 f = 1 (red triangles), σ 2 f = 4 (blue disks). The lognormal distributions predicted by (4.9) are shown in dashed lines with {µ = 0.55, ν = 0.1} for σ 2 f = 0.25, {µ = 0.7, ν = 0.2} for σ 2 f = 1 and {µ = 1., ν = 0.35} for σ 2 f = 4.

Figure 6a

, p ρ (ρ, t) has a significant probability of values less than one, which correspond to compression processes. Furthermore, its width increases with increasing time. As expected, increasing the degree of heterogeneity broadens the distribution of stretching rates shown in Figure 6b.

Non-linear stretching

The average elongation is found, from our simulations, to be well fitted by a power law in time as The profound nature of this behavior and the origin of its connection with the structure of the permeability field is out of the scope of the present work. We take equation (4.3) as an empirical fact, and note that the exponent α depends on the permeability field heterogeneity σ 2 f , as illustrated in Figure 7. Note that the power law behavior does not hold at early times, t < τ a , since the elongation ρ should converge to one as time goes to zero. This nonlinear stretching law has important implications for the geometry of the line. During its transport, the advected line is confined in a domain whose longitudinal extent grows as the advective spreading length σ(t) and whose lateral extent is constant, see also Figure 3.

The advective spreading length σ(t) is defined in terms of the longitudinal position of purely advectively transported particles x 1 (t), this means D = 0 in (2.3), as

σ 2 (t) = [x 1 (t) -x 1 (t) ] 2 .
(4.4)

The angular brackets denote the average over all particles in the initial line. As shown in Figure 8, the advective spatial variance increases non-linearly over the time scales of simulations, which is characteristic of non-Fickian spreading processes in heterogeneous media [START_REF] Le Borgne | Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF]. The anomalous scaling of advective spreading is in general persistent in time as the advective sampling of the velocity heterogeneity is not fostered by diffusion. Hence, the convergence to Fickian advective dispersion is expected to occur at a characteristic advective time corresponding to the lowest velocity in the domain. This time may be thus very large for strongly heterogeneous media. Advective spreading in porous media has been shown to be characterized by strong intermittent properties, which can be related to the flow field heterogeneity by correlated continuous time random walk approaches (de Anna et al. 2013). The non-linear evolution of the spatial variance σ(t) in the non-Fickian regime is often approximated as a power law,

σ(t) ∝ t β/2 , (4.5) 1 0 0 1 0 1 1 0 2 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3 σ 2 ~ t 1 . 0 5 σ 2 ~ t 1 . 2 σ 2 / λ 2 t / τ a σ 2 ~ t 1 . 3 5 Figure 8
. Temporal evolution of the advective spatial variance σ 2 (t), which characterizes advective spreading for σ 2 f = 0.25 (green curve), σ 2 f = 1 (blue curve) and σ 2 f = 4 (red curve).

with 1 β 2 (Figure 8). The power law evolution of the mean elongation ρ(t) and the advective spreading length σ(t) implies that the advected line support is fractal, as shown as follows. Considering a fractal line support, the number of boxes of size r needed to cover the line at a line extension σ(t) is

N (r, t) ∝ r σ(t) -d f , (4.6)
Thus, with increasing time, the number of boxes of unit size needed to cover the line increases as σ(t) d f . At the same time N (r, t) increases according to ρ(t) so that

ρ(t) ∝ σ(t) d f . (4.7)
Thus, since both ρ(t) and σ(t) evolve as a power law of time, the line should be characterized by a stable fractal dimension d f such as,

d f = 2α β . (4.8)
The fractal dimension of the advected line computed by box counting, is found to be well defined over two orders of magnitude in space, as seen in figure 9 showing the number of boxes N (r) ∼ r -d f of size r needed to cover the line support [START_REF] Villermaux | On the geometry of turbulent mixing[END_REF]. It is independent of time and ranges from d f = 1.3 for the moderate heterogeneity to d f = 1.8 for the strong heterogeneity, which is consistent with (4.8). Equation (4.8) thus provides a key relationship between stretching, spreading and the geometrical properties of the scalar field.

r/λ N (r, t)/N (1, t) r -1.3 r -1.6 r -1.8 Figure 9
. Above : Number of boxes N (r) ∼ r -d f of size r needed to cover the line support (figure 3) at times t1 = 12.5τa, t2 = 50τa and t3 = 112.5τa. Below : Same as above at time t3 = 112.5τa for different permeability field heterogeneities: σ 2 f = 0.25 (green curve), σ 2 f = 1 (blue curve) and σ 2 f = 4 (red curve).

Stochastic model of Lagrangian elongations

As discussed by [START_REF] Duplat | A nonsequential turbulent mixing process[END_REF], different types of stochastic models can be considered for modeling stretching processes. The statistics of Lagrangian elongation suggest that the considered elongation process is characterized by relatively short range correlated elongation rates γ(t), (4.1), as shown by the exponential decrease of the correlation of Lagrangian stretching rates (figure 5). Thus, we model the fluctuations of γ (t) as γ (t) = σ γ (t)ξ(t) with ξ(t) a Gaussian white noise, which renders the evolution of ρ(t) a multiplicative stochastic process through the definition of (4.1). Thus, a simple stochastic model of Lagrangian elongation, which reflects the multiplicative nature of the stretching process, and satisfies the temporal power law behavior of the mean elongation (4.3), established for t > τ a , is the following Langevin equation [START_REF] Le Borgne | Stretching, coalescence and mixing in porous media[END_REF],

1 ρ dρ dt = µ t + 2ν t • ξ(t), (4.9)
where ξ(t) is taken as a Gaussian white noise, with ξ(t) = 0 and ξ(t)ξ(t ) = δ(t -t ), which results from the observed short-range temporal correlation of computed Lagrangian stretching rates γ(t) seen in Figure 4. Considering this equation in the Stratonovich interpretation [START_REF] Risken | The Fokker-Planck Equation[END_REF], we obtain by integration of (4.9),

ρ(t) = t τ a µ exp   t τa dt 2ν t ξ(t )   (4.10)
with the initial condition ρ(t = τ a ) = 1. The mean and mean squared elongations are given by

ρ(t) = t τ a µ+ν , ρ(t) 2 = t τ a 2µ+4ν . (4.11)
Thus the exponent α in (4.3) can be expressed by α = µ + ν. The parameters µ and ν are obtained from the numerically computed Lagrangian elongations as {µ = 0.6, ν = 0.05} for σ 2 f = 0.25, {µ = 0.7, ν = 0.1} for σ 2 f = 1, and {µ = 1.1, ν = 0.15} for σ 2 f = 4. The multiplicative noise structure of (4.9) implies that z = ln(ρ) is normally distributed and therefore p ρ (ρ, t) is given by the lognormal distribution

p ρ (ρ, t) = 1 ρ exp -[ln(ρ)-µ ln(t/τa)] 2 4ν ln(t/τa) 4πν ln(t/τ a ) . (4.12)
This agrees well with the elongation PDFs obtained numerically, see Figure 6.

We can now determine the non-dimensional time τ (t), (3.3), in terms of ρ(t). To this end, we note that the Langevin equation (4.9) may be written as

d ln(ρ) du = µ + √ 2ν • ξ(u), d ln(t) du = 1, (4.13)
which renders ρ a function of u = ln(t). Thus, we can rewrite equation (3.3) for τ (t) as

τ (t) = D s 2 0 u 0 du exp(u )ρ(u ) 2 = D s 2 0 u 0 du exp[u + 2 ln(ρ)]. (4.14)
Using (4.13) to express ln(ρ) on the right side, we obtain

τ (t) = D s 2 0 u 0 du exp [u + (2µ + 4ν)u ] f (u ), (4.15) where f (u ) = exp 2 √ 2ν u 0 dyξ(y) -4νu
. By a shift of the integration variable, the integral expression of τ (t) can be written as

τ (t) = D s 2 0 exp[u(1 + 2µ + 4ν)] u 0 du exp[-u (1 + 2µ + 4ν)]f (u -u ). (4.16)
The main contribution of the integrand comes from u = 0 because the random function f (u) > 0 is in average equal to 1, its standard deviation goes as exp(4νu). Thus, we obtain approximately

τ (t) ≈ D s 2 0 exp[u(1 + 2µ + 4ν)]f (u) 1 + 2µ + 4ν = D s 2 0 ρ(t) 2 t 1 + 2µ + 4ν . (4.17)
Using this approximation, the maximum lamella concentration θ(t) given by (3.4) can be related to the current lamella elongation as

θ(t) ≈ ρ(t) -1 (At/τ a ) -1/2 , A = Dτ a s 2 0 c 2 0 4 1 + 2µ + 4ν , (4.18)
for ρ(t) 1. This scaling trivially expresses mass conservation in two dimensions since at the lamella scale, the product θ √ Dtρ remains constant in time. From (4.9), we obtain for θ(t) the evolution equation

1 θ dθ dt = - µ t - 1 2t + 2ν t • ξ(t), (4.19) with θ(t = τ a ) = A -1/2 ∼ c 0 . The PDF of θ(t) is lognormal, p θ (θ, t) = exp -[ln(θ/c0)-mz] 2 2σ 2 z θ 2πσ 2 z . (4.20)
where

m z (t) = -µ ln(t/τ a ) - 1 2 ln(t/τ a ), σ 2 z (t) = 2ν ln(t/τ a ) (4.21)
are the mean and variance of z(t) = ln[θ(t)]. The moments of the maximum lamella concentration across the stretched lamellae are given by

θ(t) q = c q 0 exp(qm z + q 2 σ 2 z /2) = c q 0 t τ a -qµ-q 2 +q 2 ν . (4.22)
A similar construction can be found in [START_REF] Kalda | Simple model of intermittent passive scalar turbulence[END_REF] and [START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF].

Mixing scale, mixing time

The derived model of Lagrangian elongations can be used to quantify the temporal evolution of the local mixing scale, whose evolution is governed by the competition between diffusive broadening and substrate compression perpendicular to the elongation direction [START_REF] Le Borgne | Persistence of incomplete mixing: A key to anomalous transport[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF]). These two effects equilibrate at the mixing time, when the transverse lamellae thickness has reached the Batchelor scale [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in a turbulent fluid. part 1. general discussion and the case of small conductivity[END_REF]). This competition is expressed by the advection-diffusion equation (3.1). The local scale m (t) is defined from the variance of the concentration profile across a lamellae [START_REF] Villermaux | Mixing by porous media[END_REF]) as a measure of the lamellae 'thickness',

2 m (t) = dζζ 2 c(ζ, t)/ dζc(ζ, t). Using expression (3.2) for c(ζ, t), we have, 2 m (t) = s 2 0 2ρ(t) 2 + 2s 2 0 τ (t)
ρ(t) 2 . Further using approximation (4.17) for τ (t) yields

2 m (t) = s 2 0 2ρ(t) 2 + 2Dt 1 + 2µ + 4ν . (5.1)
Note that the elongation ρ(t) is lognormally distributed and given by (4.12). We obtain for the evolution of the average mixing scale s m (t) = m (t) from (5.1) the approximate expression

s m (t) ≈ s 2 0 2 t τ a -2µ+2ν + 2Dt 1 + 2µ + 4ν
.

(5.2)

Figure 10 shows the temporal evolution of the average mixing scale for different Péclet numbers and permeability field variances. Since the compression rate is stronger than diffusive growth at early time, the mixing scale s m decays until it reaches the Bachelor scale at which the diffusive growth balances compressive effects due to stretching. This occurs at the mixing time τ m when,

ds m (τ m ) dt = 0, τ m τ a ∝ P e 1 1+2µ-2ν 0 ,
(5.3)

where we defined the strip Péclet number P e 0 = s 2 0 /τ a D. The ratio between the mixing time τ m and the diffusion time τ D can be expressed as τ m /τ D = (s 0 /λ) 2/(1+2µ-2ν) P e -(2µ-2ν)/(1+2µ-2ν) . Since the exponent -(2µ -2ν)/(1 + 2µ -2ν) is negative, the difference between τ m and τ D is expected to increase with Péclet number. For large Péclet numbers P e 1, the mixing time occurs much before the diffusion time. On the other hand, for small Péclet numbers P e 1, the diffusion time may occur before the mixing time. In this case, we expect stretching to have a minor effect on the concentration field, as diffusion would be the dominating process.

The mixing scale at time τ m s B = s m (τ m ) is the Batchelor scale, given by

s B s 0 ∝ P e 0 2ν-µ 1+2µ-2ν , (5.4)
Before the mixing time, the decay of the decay of the average mixing scale is controlled by the elongation process, s m ∼ ρ -1 ∼ t -µ-ν , which depends on the permeability field heterogeneity (Figure 10b). After the mixing time, the average mixing scale grows as

s m = 2Dt
1+2µ+4ν for all degrees of heterogeneity. The explicit relation (5.1) between the mixing scale m (t) and the elongation ρ(t) allows deriving the PDF of m (t) by mapping from the PDF (4.12) of ρ(t). Thus, we for P e = 10 4 at times t = 2τa, 4τa, 8τa, 16τa, 32τa, 64τa, 128τa.

obtain the shifted lognormal PDF

p m ( m , t) = m exp - [ln( 2 m -2Dt 1+2µ+4ν )-ln(s 2 0 /2)+2µ ln(t/τa)] 2 16ν ln(t/τa) [ 2 m -2Dt/(1 + 2µ + 4ν)] 4πν ln(t/τ a ) , (5.5)
for m 2Dt 1+2µ+4ν and zero else. This lower limit corresponds to the most stretched fraction of lamellae which have attained their mixing time and are growing diffusively. The temporal evolution of the mixing scale PDF is illustrated in figure 11 for P e = 10 4 . At early times, the distribution of mixing scales is broad reflecting the broad distribution of stretching rates and line elongation. With increasing time, the peak of the PDF moves towards smaller values of m as a consequence of lamella compression. For times larger than the mixing time τ m diffusion wins over compression in the evolution of the mixing scale, m = 2Dt 1+2µ+4ν . The random nature of the stretching process has only a marginal influence on the mixing scale, which evolves now deterministically. Therefore, the PDF p m ( m ) converges asymptotically to a Dirac-delta distribution

p m ( m ) = δ m - 2Dt 1 + 2µ + 4ν (5.6)
expressing that all lamellae ultimately experience the same history as they have sampled the space of possible stretching values.

Concentration PDF of an ensemble of non-interacting lamellae

For large Péclet numbers, a well developed regime exists in which lamellae do not overlap as illustrated in Figure 1b. We describe here a corresponding mixing model which maps the PDF of lamella elongations into the PDF of concentrations. The predictions are compared to the concentration data of the full flow and transport simulations.

In this regime, the concentration PDF p c (c, t) is obtained by sampling concentration over individual lamellae, whose PDF (3.5) is conditional to the respective maximum concentration θ(t). The maximum concentration describes a stochastic process due to the stochasticity of the elongation ρ(t). The PDF of θ(t) is denoted by p θ (θ, t). The concentration PDF is thus written as

p c (c, t) = ∞ c dθp c (c|θ)p θ (θ, t). (6.1)
Using (4.20) in (6.1) and rescaling the integration variable, we can write the concentration PDF as

p c (c, t) = ∞ ln(c) dz exp -(z-mz) 2 2σ 2 z 2c 2πσ 2 z [z -ln( )][z -ln(c)] (6.2)
where we used (3.5) for p c (c|θ). For concentrations ln(c) < m z , the PDF can be approximated by

p c (c, t) ≈ 1 2c [m z -ln( )][m z -ln(c)] . (6.3)
Thus, it shows the 1/c-decay characteristic for the diffusive profile (3.5) across the lamella.

For concentrations ln(c) > m z , we obtain the approximation

p c (c, t) ≈ exp -[ln(c)-mz] 2 2σ 2 z c 8 ln(c/ )[ln(c) -m z ] . (6.4)
For large concentrations, the PDF decreases essentially as the lognormal distribution. The concentration PDF (6.2) is parameterized only by the parameters µ and ν in (4.9), which are obtained from the evolution of ρ(t). Its predictions are in good agreement with the numerical simulations for all cases, in which the stretched lamellae that constitute the line evolve on their own and do not overlap, i.e. at early times and high Péclet numbers. This is illustrated in Figure 12). The model predicts correctly the temporal evolution of the shape of the concentration PDF, as it evolves towards uniformity.

The behavior of p c (c, t) for small concentrations is given by p c ∝ 1/[c ln(1/c)], which is characteristic of the Gaussian concentration profiles across each lamella as given in (3.2). The PDF of large concentrations is controlled by the lognormal PDF of maximum lamella concentrations p m (θ), as expressed by (6.4). Comparison with the prediction of a purely diffusive model (inset of figure12b) shows that flow heterogeneity accelerates considerably the temporal decay of concentration levels through stretching enhanced diffusive mixing. Furthermore, the model predicts accurately the impact of heterogeneity on mixing and represents well the decrease of the maximum concentration as the degree of heterogeneity increases. This is illustrated in Figure 12b.

A corresponding similar example of broad concentration statistics reflecting solely the distributed histories among non interacting lamellae is provided by mixing in the dispersing turbulent plume studied by [START_REF] Duplat | A nonsequential turbulent mixing process[END_REF].

Concentration PDF of an ensemble of interacting lamellae

Interactions between lamellae are enforced at later times and for smaller Péclet numbers because the average elongation ρ ∝ t α grows faster than the longitudinal size of the dispersion area σ ∝ t α/d f in which the line is confined. Recall that d f 1 as given in (4.8). When the lamella width is larger than the distance between lamellae, the concentration fields around each lamella start coalescing to form lamella aggregates as illustrated in Figure 13. Due to continued diffusive broadening these lamella aggregates may overlap to form bundles. We develop here a mixing theory in which lamella interaction is modeled as a random aggregation process with partial or complete coalescence.

The coalescence processes and model predictions are tested by comparison to detailed numerical simulations.

Concentration PDF of lamella aggregates

This process can be described by the random aggregation process proposed in [START_REF] Villermaux | Mixing as an aggregation process[END_REF]. We consider an aggregate of n lamellae, characterized by the maximum concentration c m . The value of c m results from the addition of concentration levels of the individual lamella. Assuming that the lamellae that overlap have uncorrelated stretching histories, these concentration levels can be considered as uncorrelated random variables, and the coalescence can be described as a random aggregation process (e.g., [START_REF] Villermaux | Mixing as an aggregation process[END_REF]. The concentration PDF p m [c m , t|n] then is given by a Gamma distribution Figure 14. PDF of local maximum concentrations p(cm, t) Péclet number P e = 80. The PDF is plotted for times t1 = 12.5τa (red triangles), t2 = 37.5τa (green diamonds), t3 = 62.5τa (blue squares) and t4 = 112.5τa (black circles). The dashed lines represent the prediction of the Gamma distribution (7.1) for an average number of overlaps obtained from equation (7.9).

of order n independently of the initial distribution,

p m (c m , t|n) = c n-1 m θ(t) n Γ[n] exp - c m θ(t) , (7.1)
where θ(t) is the average concentration of an elementary lamella at time t given by (4.22) for q = 1, which derives from the stretching enhanced diffusive regime described previously.

We approximate the concentration profile around each aggregate by a Gaussian spatial profile. Thus, the concentration PDF p c (c|c m ) conditioned on the maximum concentration c m is given by (3.5). The concentration PDF across lamella aggregates containing n(t) lamellae is then given by for c . To test the validity of equation ( 7.1), we evaluate numerically the PDF of local concentration maxima. To do so, we compute the PDF p m (c m , x 1 , t) of concentration maxima at positions x 1 along the mean flow by sampling concentration values along x 2 , the direction perpendicular to the mean flow. The PDF of concentration maxima p m (c m , t) then is obtained by averaging p m (c m , x 1 , t) over x 1 . Figure 14 shows the resulting concentration PDF p m (c m , t). The distribution of local maxima narrows down with time, following the characteristic shape of the Gamma distribution. The prediction of equation ( 7.1), is in good agreement with the numerical simulations. The Gamma distribution of local concentration maxima thus supports the random addition nature of the lamella aggregation process. In order to obtain an estimate for the distribution of local maxima from (7.1), we require the number of lamellae across bundles n. We set the number of lamellae equal to the average number of lamellae n across bundles. In general, the number of lamellae is dependent on the position x 1 at which the cross-section is taken. This is discussed in the following Section.

p c (c, t|n) = c0 c dc m c n-1 m θ(t) n Γ[n]

Analysis of the distribution of lamella overlaps

The number of lamellae n[x(t|b), r, t] within a radius r around the position x(t|b) in the purely advectively transported line is obtained by counting all particles at locations x(t|a) in the line that are within a radius r and subsequent normalization by r, as illustrated in figure 15. Note that x(t|a) is the trajectory of a purely advectively transported particle, i.e., D = 0 in (2.3). The initial particle distribution is a uniform line at x 1 (t = 0|a) = λ and x 2 (t = 0|a) = a. Note also that the number of particles within a radius r of the position x(t|b) includes the particles inside the same lamella, whose number increases according to r. The normalization by r corrects for this. Figure 16 shows maps of n[x(t|b), r, t] along the stretched and folded line, as well as the projection of the number of lamellae n(x 1 , r, t) at a position x 1 along the mean flow direction at time t = 112.5τ a and for r = λ/4, λ and 4λ. The number of overlaps is distributed along the line with large overlap numbers in regions of high lamellae density, where the line is folded several times over itself. For small r (see Figure 16a), the number of overlaps n(x 1 , r, t) is found to be approximately constant along the longitudinal direction. When r reaches the characteristic heterogeneity length λ, the number of overlaps increases drastically at the center of the plume as illustrated in Figures 16b and17c This phenomenon can be seen as a percolation process [START_REF] Villermaux | Mixing by porous media[END_REF]. While for small r lamella aggregates are isolated from each other, they become connected with increasing r to form lamella bundles. This process is most active a the center of the plume where the density of lamellae tends to be larger. The transition to percolation occurs approximately at r = λ. Note that lamellae within a radius r coalesce only if the lamella width (5.1) is equal to r. This means for the present regime that r ≈ √ Dt. Thus, percolation occurs at the characteristic diffusion time τ D = λ 2 /D. For r > λ, i.e., at times t > τ D , n(x 1 , r, t) follows a distribution which is peaked at the center of the plume and is sharply cut-off with distance from the peak position, see Figure 17a. In this regime, the projected density of purely advectively transported particles c a (x 1 , t) is proportional to the density n(x 1 , r, t) of lamellae at x 1 , c a (x 1 , t) ≈ n(x 1 , r, t)/n(r, t) with n(r, t) the average number of lamellae within a radius r. Thus, we obtain for n(x 1 , r, t) the approximation n(x 1 , r, t) ≈ n(r, t)c a (x 1 , t), (7.3) Figure 17b compares n(x 1 , r, t)/n(r, t) and c a (x 1 , t) at time t = 112, 5τ a , and confirms (7.3). The PDF p n (n, t, r) of the number of lamellae in the time regime t > τ D can be estimated as,

p n (n, t, r) = dx 1 (c a ) dc a ca=n/n(r,t) . (7.4)
The longitudinal concentration distribution c a (x 1 , t) is often approximated by a Gaussian distribution in the Fickian dispersion framework. However, for the purely advective case, the upstream and downstream tails are determined respectively by the minimum and maximum velocities. Thus, the longitudinal concentration distribution does not present infinite tails, as does the Gaussian distribution. In fact, according to Figure 17b, we approximate c a (x, t) by a linear distribution with |dx(c a )/dc a | = const. Thus for t > τ D , p n (n, t, r) is given by the uniform distribution p n (n, r, t) = 1 2n(r, t) for 0 n < 2n(r, t). For times t < τ D , the spatial profile of n(x 1 , r, t) is uniform and thus the PDF p n (n, t, r) can be approximated by the delta density

p n (n, t, r) = δ[n -n(r, t)].
(7.6)

Temporal evolution of the average number of lamella overlaps

The average number n(r, t) of lamellae with concentration θ(t) having coalesced in a radius r is such that it restores, by addition, the average mixture concentration c(x, t) so that c(x, t) = n(r, t) θ(t) .

(7.7)

For times below the percolation threshold t < τ D , the number of lamellae within a radius r is conditioned by the geometry of the folded and stretched line support. In order to determine the number of lamellae n(r, t), we consider the number n p (r, t) of particle pairs in the purely advectively transported line that are separated by a distance smaller than r. Thus, the number of lamellae within a radius r is given by n(r, t) = n p (r, t)/r, where the normalization by r is done to correct for the particles which belong to the same lamella as in the previous section. The number of particle pairs n p (r, t) over a fractal support scales as n p (r, t) ∼ r d f where d f is the fractal dimension of the line support [START_REF] Grassberger | Characterization of strange attractors[END_REF]. Therefore the average of lamella overlaps behaves as n(r, t) ∼ r d f -1 .

(7.8)

Hence, the average number of lamella overlaps is determined by the fractal dimension of the line support, which is set by the stretching and spreading dynamics through equation (4.8). This spatial scaling is consistent with numerical estimates (figure 18). Compression perpendicular to the direction of elongation implies that lamellae may approach each other as close as the Batchelor scale in equation (5.4). Hence, the first lamella overlap is expected to occur at the mixing time, when diffusive growth overcomes compression. Thus, at t = τ m and r = s B (equations (5.3) and (5.4)), the number of overlap is expected to be on average equal to one. After the mixing time, the width of the diffusive profile of a single lamella (5.1) increases like √ Dt, which sets r. Thus, we . Average number of overlaps n(r, t) for times t1 = 12.5τa (red triangles), t2 = 37.5τa (green diamonds), t3 = 62.5τa (blue squares) and t4 = 112.5τa (black disks). The dashed line represents the prediction of equation (7.9).

obtain for the temporal evolution of the mean number of coalescing lamellae,

n c (t) = Dt/s 2 B (d f -1)/2 , for τ m < t < τ D . (7.9)
The mixing time τ m (equation (5.3)) sets the characteristic time for the transition from stretching enhanced mixing to lamellae coalescence. Note that the time range τ m < t < τ D spans two orders of magnitude for P e = 8 and four orders of magnitude for P e = 8 × 10 4 . Hence, over this range of time scales, we expect the concentration field to be composed of aggregates of coalesced lamellae in areas where the line is densely folded over itself, and concentration lacunarities in low flow zones which have not been filled by diffusion.

For times above the percolation threshold t > τ D , most of the concentration lacunarities have been filled at the center of the plume. Hence, the average concentration in the plume is estimated from the longitudinal plume extension σ as c ∝ 1/σ. Thus, in this regime, the average number of colaescences that restores the average concentration is obtained from equation (7.7) as

n c (t) ∝ 1 σ(t) θ(t)
, for t > τ D (7.10)

Concentration PDF

The overall concentration PDF depends on the PDF of the number of lamellae in a given distance r (7.4). The density of lamella coalescences at a given time t is obtained by setting r = √ 2Dt and substituting n(r, t) by n c (t) in (7.5) and (7.6). Partial Coalescence Regime For τ m < t < τ D , the number of overlaps is approximately constant along the line and the density of lamellae is such that it allows an essentially diffusive concentration profile within the aggregate. Thus the concentration PDF is obtained by setting n = n c (t) in (7.2) as 7.11) with n c given by (7.9). For concentration c < n c θ , we can approximate the PDF (7.11) by

p c (c, t) = c0 c dc m c nc-1 m θ nc Γ(n c ) exp(-c m / θ ) 2c ln(c m / ) ln (c m /c) , (
p c (c, t) ≈ 1 2c ln[(n c -1) θ / ] ln[(n c -1) θ /c] (7.12)
Thus, for small concentrations, the PDF is dominated by the diffusive concentration profile that persists in the gaps between the lamella aggregates. For larger concentrations c > n c θ , the PDF (7.11) behaves as

p c (c, t) ≈ c nc-1 2 θ nc Γ(n c ) exp(-c/ θ ) ln(c/ ) √ π c/ θ -(n c -1) , (7.13)
with a large concentrations exponential decay characteristic of Gamma distribution.

The prediction of equation ( 7.11) for the temporal evolution of the concentration PDF is found to be in good agreement with the numerical simulations as illustrated in Figure 19. In particular, the large concentration behavior of the PDF shows the exponential decay described in (7.13) (see inset of figure 19a). This property, which is not observed in the first regime when the line does not interact with itself (see inset of figure 12a), is a characteristics of the coalescence process described by (7.1). As outlined above, the low concentration behavior (7.12) results from the diffusive profiles around the lamella aggregates. The comparison with the prediction of (6.2) (see inset of figure 19b), which does not account for lamella interactions, shows that the coalescence process tends to increase the intermediate concentration probabilities around the average at the expense of the probability of high and low concentration levels. Percolating Coalescence Regime For t > τ D , the concentration field is composed of bundles of aggregated lamellae, whose number varies along the mean flow direction according to (7.5). Aggregates with a number of lamellae that is smaller than or equal to n c (τ D ) = (λ/s B ) d f -1 , see (7.9), are characterized by a diffusive profile across the aggregate and thus the concentration PDF is given by p c (c, t|n), see (7.2). For n > n c (τ D ), this means above the percolation threshold, the lamella density in the aggregate is so large that there is no space for a diffusive profile to develop between lamellae and therefore, the concentration PDF is given by the PDF (7.1) of concentration maxima p m (c, t|n) with n = n c (t) given by (7.10)

p c (c, t) = c nc-1 θ nc Γ(n c ) exp[-c/ θ ]. (7.14)
We expect this distribution of concentration values around the center of the plume because there the lamella density is highest, see Figure 16. This prediction is verified in figure 20, where the concentration PDF is computed at the center of the plume for a Péclet number P e = 8 at different times. In order to analyze the concentration PDF in a region where full percolation occurs, we evaluate the concentration statistics at the central longitudinal positions x, defined by

x 1 (x 2 )-σ x1 (x 2 )/2 < x 1 < x 1 (x 2 )-σ x1 (x 2 )/2
where x 1 (x 2 ) and σ x1 (x 2 ) are respectively the local center of mass and spatial standard deviation of the concentration field at transverse position x 2 . The prediction of equation (7.14) are found to be in good agreement with the numerical results, with the average coalescence number given by equation (7.10).

The global concentration PDF in this regime integrates percolating coalescence processes in the center of the plume and diffusion in open areas at the back and fore fronts of the plume. Thus, it is given by the weighted sum of the concentration PDFs (7.11) and (7.14) in the partial coalescence and in the percolating coalescence regimes as described in Appendix A. The prediction of equation (A 2) for the temporal evolution of the concentration PDF is found to be in good agreement with the numerical simula- tions (figure 21). The two coalescence models are compared in the inset figure 21a. Both models predict a 1/c scaling at small concentrations, which is characteristic of diffusive profiles, and an exponential decay at large concentrations, which results from the coalescence process. The percolating coalescence model (t > τ D , equation (A 2)) shows a smaller contribution from the 1/c diffusive process, which is acting only on the plume edges. Furthermore, in this model, intermediate concentrations have a close to uniform probability, which results from the uniform distribution of coalescence numbers within the plume.

A corresponding similar example of a mixture built by a random aggregation process in a turbulent channel was studied by [START_REF] Duplat | Mixing by random stirring in confined mixtures[END_REF].

Concentration moments

The knowledge of the concentration PDF at any times allows predicting the temporal evolution of the concentration moments. In this section, we compare these theoretical predictions to the result of numerical simulations in all regimes and we develop explicit expressions for the asymptotic scaling behaviors of the concentration moments.

The moments of concentration are defined by

c(x, t) q = 1 V c (t) c(x,t)> dxc(x, t) q , (8.1)
where V c (t) is the characteristic volume (2.4) occupied by the solute. This expression can be written in terms of the concentration PDF p c (c, t) as

c(x, t) q = ∞ 0 dcc q p c (c, t). (8.2)
The average scalar dissipation rate is defined by [START_REF] Zeldovich | The asymptotic law of heat transfer at small velocities in the finite domain problem[END_REF])

χ(t) = 1 V c (t) c(x,t)> dx2D [∇c(x, t)] 2 . (8.3)
This expression for the scalar dissipation rate can be approximated in terms of the mean and mean squared concentration by noting that

c(x, t) ≈ V c (t) -1 , c(x,t)> dx2D [∇c(x, t)] 2 ≈ - d dt c(x,t)> dxc(x, t) 2 , (8.4)
where the latter relation follows directly from the Fokker-Planck equation (2.1). Using these relations, the scalar dissipation rate can be approximated as (see also [START_REF] Villermaux | Mixing by porous media[END_REF])

χ(t) ≈ c(x, t) 2 d ln c(x, t) dt - d c(x, t) 2 dt (8.5)
thus making a link between χ, c and c 2 . Figure 22 compares the simulated temporal evolutions of the concentration moments and scalar dissipation rates to the predictions obtained from expressions (6.2), (7.11), and (A 2) for the concentration PDF in the different mixing regimes corresponding to different Péclet numbers in the case of the intermediate heterogeneity, σ 2 ln K = 1. The predicted and observed temporal evolutions are in very good agreement, which confirms the validity of the theoretical mixing approaches developed for the different regimes. In the following, we derive explicit expressions for the asymptotic scaling behavior of concentration moments and scalar dissipation rate.

Non-interacting lamellae

The concentration PDF in the non-interacting regime is given by (6.1) in terms of the diffusive profiles of individual lamellae weighted by the probability of their maximum concentration. For concentrations ln(c) < m z and ln(c) > m z with m z given by (4.21), we derived the approximations (6.3) and (6.4), respectively. Thus, we write the q-th concentration moment as

c(x, t) q = exp(mz) 0 dcc q p c (c, t) + ∞ exp(mz) dcc q p c (c, t) (8.6)
Inserting the approximations (6.3) and (6.4) and evaluating the remaining integrations gives the approximate expression c(x, t) q = π q

(t/τ a ) -qµ-q 2 -µ ln(t/τ a ) -ln( ) + π 2 (t/τ a ) -qµ-q 2 +q 2 ν erfc(q/ √ 2) 8q[-µ ln(t/τ a ) -ln( )] . (8.7)
Note that the minimum concentration (t/τ a ) -µ and thus -µ ln(t/τ a ) -ln( ) > 0. It scales in leading order as c(x, t) q ∼ θ(t) q ∼ t -qµ-q 2 +q 2 ν , see (4.22). Thus, we obtain for the mean and mean squared concentrations the approximate scaling

c(x, t) ≈ t τ a -µ+ν-1/2 , c(x, t) 2 ≈ t τ a -2µ+4ν-1 . (8.8)
The scalar dissipation rate scales as

χ(t) ≈ t τ a -2µ+4ν-2
.

(8.9)

Partial Coalescence Regime

In the partial coalescence regime, τ m < t < τ D , the temporal scaling of the concentration moments is affected by the interaction between lamellae. The concentration PDF in this regime is given by (7.11) and can be approximated by (7.12) for c < n c θ and by (7.13) for c > n c θ . As above, we write the q-th concentration moment as

c q (x, t) = (nc-1) θ 0 dcc q p c (c, t) + ∞ (nc-1) θ dcc q p c (c, t) (8.10)
Inserting the approximations (7.12) and (7.13), we obtain the following approximate expression for the concentration moments

c(x, t) q ≈ (n c -1) q θ(t) q √ πerf q ln[n c -1) θ(t) / ] 2 q ln[n c -1) θ(t) / ] + √ π 2 q ln[(n c -1 + q) θ(t) / ] θ(t) q Γ(n c + q, n c -1) Γ(n c ) . (8.11)
At large times, we obtain at leading order the asymptotic scaling behavior

c(x, t) q ∼ n q c θ(t) q ∼ t q d f -1 2 -qµ-q 2 +qν , (8.12)
where we used q = 1 in (4.22) to obtain θ(t) ∼ (t/τ a ) -µ-1/2+ν . This gives for the mean and mean squared concentrations at leading order at large times,

c(x, t) ∼ t -µ+ν+d f /2-1 , c(x, t) 2 ∼ t -2µ+2ν+d f -2 . (8.13)
For the scalar dissipation rate, we obtain accordingly the asymptotic scaling:

χ(t) ∼ t -2µ+2ν+d f -3 . (8.14)

Percolating Coalescence Regime

In the percolation coalescence regime for t > τ D the concentration PDF at the center of the plume, which determine the concentration moments in this regime, is given by the Gamma-PDF (7.14). Thus, the concentration moments are given by the moments of the Gamma-distribution as

t/τ a < c > < c > < c 2 > < c 2 > < χ > < χ > < c > < χ > a) < c 2 > c) b)
c(x, t) q = θ(t) q Γ(n c + q) Γ(n c ) , (8.15)
where the number n c of lamella aggregates is given by (7.10). Thus, we obtain for the mean c(x, t) = θ(t) n c and mean squared c(x, t) 2 = θ(t) 2 n c (n c + 1) concentrations

c(x, t) = 1 σ(t) , c(x, t) 2 = 1 σ 2 (t) + θ(t) σ(t) (8.16)
where we used (7.10). Inserting the scaling (4.4) of σ(t) in these expression, we obtain in leading order the asymptotic scaling, c(x, t) ∼ t -β/2 , c(x, t) 2 ∼ t -β . (8.17)

Thus, the scalar dissipation rate scales asymptotically as

χ(t) ∼ t -β-1 (8.18)

Impact of heterogeneity on mixing

The expressions derived for the scaling exponents in the different regimes quantify the impact of the average stretching rate (exponent µ), the variability of stretching rates (exponent ν) and the line geometry (fractal dimension d f ) on the temporal behavior of the concentration moments and of the scalar dissipation rate. The asymptotic scaling exponents are in good agreements with the results of the numerical simulations for the different structural heterogeneities under consideration, as shown in Figure 23 for the second moment of concentration. These expressions allow understanding the effect of the permeability field heterogeneity on the effective mixing rates. The structural heterogeneity clearly leads to enhanced mixing as the stretching exponent increases with the permeability field variance. On the other hand, the coalescence process tends to decrease the absolute value of the scaling exponent characterizing the temporal evolution of the second moment of concentration (as seen by comparing equations (8.8) and (8.13)). Thus, the coalescence process tends to slow down the concentration variance decay. In the first regime, diffusion mixes the solute contained in the lamella with the surrounding zero concentration fluid. In the coalescence regime, the action of diffusion is partly to merge lamellae of different concentrations with a reduced net effect on mixing compared to the stretching regime.

Mixing measures

The knowledge of the concentration PDF allows predicting the temporal evolution of different metrics classically used for quantifying mixing. These relie on quantities directly computed from the full concentration content of the mixture, of which they are often used as Ersatz. We examine two classical measures :

Intensity of segregation

The intensity of segregation originally imagined by [START_REF] Danckwerts | The definition and measurement of some characteristics of mixtures[END_REF] is

I(t) = c(x, t) 2 -c(x, t) 2 c(x, t) (c 0 -c(x, t) ) (9.1)
It is a convenient index going from 1 at initial segregation, when the scale is not mixed at all and the PDF of concentration is bimodal, presenting two levels of concentration only (c = 1 for the injection concentration,and c = 0 for the diluting medium), to 0 at complete homogeneity, when the concentration PDF presents a single concentration level, at the average concentration.

The first and second moment of the full concentration PDF are known, and therefore the intensity of segregation.

Dilution index

Another surrogate of the concentration PDF is the dilution index E(t) [START_REF] Kitanidis | The concept of the dilution index[END_REF])

E(t) = exp [H(t)] , (9.2)
where H(t) is the entropy of the normalized concentration distribution p(x, t)

H(t) = -dxp(x, t) ln[p(x, t)]Θ[c(x, t) -], p(x, t) = c(x, t) c(x, t) V c (t) . (9.3)
This index, which increases as the mixture progresses towards uniformity, can be written in terms of the concentration c(x, t) as

H(t) = - 1 c(x, t) V c (t) dxc(x, t) ln[c(x, t)]Θ[c(x, t) -] -ln[ c(x, t) V c (t)]. (9.4)
The latter can be written in terms of the concentration PDF as

H(t) = - 1 c(x, t) ∞ 0 dcc ln(c)p c (c, t) -ln[ c(x, t) V c (t)].
(9.5)

Figure 24 displays the temporal evolution of the dilution index and the segregation intensity and their predictions by the mixing models corresponding to the different regimes. The predictions are in good agreement with the numerical simulations.

Spatial structure of the concentration fields

While concentration moments quantify the global mixing behavior, they do not give information on the spatial structure of the concentration field. The latter can be quantified by studying the PDF of concentration increments and their structure functions [START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF]). In this section we present explicit expressions for the in-crement PDFs and structure functions based on the mixing models developed in previous sections. Model predictions are then tested against numerical simulations.

The concentration increments are quantified for different spatial increments ∆x as ∆c(x, t) = c(x + ∆x, t) -c(x, t). Figure 25 and 28 show the PDF of concentration increments p ∆x (∆c) computed for different spatial increments ∆x = |∆x|, for respectively large and small Péclet numbers. The dependence of the concentration increment PDFs on the spatial scale ∆x can be quantified by analyzing the spatial scaling of the structure functions, defined as, |∆c(∆x)| q = d∆cp ∆x (∆c)|∆c| q .

(10.1) 10.1. Non-interacting lamellae

For large Péclet numbers, the concentration increments PDFs are sharply peaked at c = 0, and present non Gaussian tails (figure 25). When increasing the spatial increment ∆x, the increment PDFs widen as larger concentration differences are encountered. However, over a characteristic scale ∆x > η, the increment PDFs are found to stabilize and become independent of ∆x (η = 0.4λ in figure 25). This implies that the concentration values are uncorrelated for distances larger than η [START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF]). For uncorrelated concentration values, the PDF of concentration increments can be obtained from the convolution of the global concentration PDF,

p ∆x>η (∆c) = dc p c (c , t)p c (|∆c| -c , t). (10.2)
where p c is given by equation (6.2). For distances |∆x| smaller than η, the concentration PDF may be estimated by a simple rescaling of the concentration scale by a factor |∆x|/η,

p ∆x<η (∆c) = η ∆x p c η ∆x ∆c (10.3)
which corresponds to an approximation of the local Gaussian distribution in lamellae by a triangular function. In the large Péclet regime, this approximation is found to give good predictions for all spatial increments (figure 25). The shape of the concentration increment PDF is in fact independent of the spatial increment ∆x when represented as a function of ∆c/ |∆c| . This is consistent with the scaling form of the increment PDF (10.3) for |∆x| smaller than η, and its independence on the spatial increment for |∆x| larger than η. Equation (10.3) implies that the structure functions (figure 26) scale as |∆c| q ∼ |∆x| q (10.4) for ∆x < η and are constant for ∆x > η.

Coalescence regime

For smaller Péclet numbers, the correlation length η is much larger than in the noninteracting regime (figure 30) due to the large number of coalescences experienced by lamellae in this regime. For ∆x > η, the concentration increments are differences of independent values and the PDF of concentration increment is derived by a simple convolution as in the previous regime (equation (10.2)). Below the correlation length, ∆x < η, concentration fluctuations result from the random aggregation process and are thus highly correlated. Based on the elementary construction process of concentration fields, we derive here a new framework for quantifying concentration increments under well developed coalescence processes.

As outlined in Section 7.1, in the percolating coalescence regime, the concentration

∆x/λ |∆c(∆x)| q / |∆c(η)| q η q = 6 q = 1 Figure 26.
Structure functions |∆C| q computed as a function of ∆x, for q ranging from 1 to 6, for a permeability field variance σ 2 ln K = 1 and a Péclet number P e = 8 × 10 3 at time t = 112.5τa. The structure functions are normalized by their asymptotic value for large ∆x, |∆c(η)| q . c(x, t) at a point x results from the random addition of concentration levels of n c (t) lamellae within a radius r = √ Dt. Hence, two disks with radius r = √ Dt, located at x and x + ∆x with ∆x < η have the lamellae in common which are located in their intersection (figure 27). Thus, c(x, t) and c(x + ∆x, t) are each composed by the random addition of lamellae from the disjoint set of the two disk, and random addition of lamellae in the intersection. The latter contribution cancels out in the concentration increment, which now can be written as the difference of the two independent concentrations c (x + ∆x, t) and c (x, t), which are obtained by random addition of n ∆x lamellae in the respective disjoint disk sets, ∆c(x, t|∆x) = c (x + ∆x, t) -c (x, t), (10.5)

The typical size of the non overlapping regions of two lamellae aggregates whose centers are separated by the distance ∆x is approximately ∆x (figure 27). Thus, the number of independent lamellae n ∆x participating in the construction of the concentration increment is expected to be, from equation (7.8),

n ∆x = ∆x s B d f -1 (10.6)
The concentration levels c (x, t) and c (x + ∆x, t) are statistically independent because they are composed of independent lamellae. Thus, the PDF of the concentration increment ∆c(x, t) can now be be determined by the convolution of the PDFs of the concentrations around the independent aggregates, p ∆x<η (∆c) = dc p c (c , t|n ∆x )p c (∆c -c , t|n ∆x ). (10.7)

In the percolating coalescence regime the concentration PDF is characterized by the Gamma PDF of maximum concentrations given by (7.1) as

p c (c|n) = c n-1 θ n Γ(n) exp - c θ (10.8)
Inserting (10.8) into (10.7) gives the increment PDF for ∆x < η as

p ∆x<η (∆c|n) = 1 √ π θ 2n∆x Γ(n ∆x ) |∆c| θ 2 n∆x-1/2 K n∆x-1 2 |∆c| θ , ( 10
.9)

Similarly, using the latter in expression (10.1), we have for the moments

|∆c| q = 2 q θ q √ πΓ(n ∆x ) Γ q + 2n ∆x 2 Γ 1 + q 2 (10.10)
For large n ∆x 1, the structure functions can be approximated by

|∆c| q ≈ 2 q θ q n q 2 ∆x √ π Γ 1 + q 2 (10.11)
From equation (10.6), we obtain, for ∆x < η,

|∆c| q ≈ ∆x s B (d f -1) q 2 2 q θ q √ π Γ 1 + q 2 (10.12)
These predictions are in good agreement with the results of the numerical simulation for P e = 8, in the center of the plume where full percolation is established in this regime (figure 28 and 29). Hence, expressions (10.3) and (10.9) provide analytical predictions for the concentration increment PDFs in the two end-member regimes of non-interacting lamellae and percolating coalescence. The latter provide a quantitative description of the spatial organization of the concentration fields and the distribution of concentration gradients. As a corollary, we obtain analytical predictions for the corresponding structure functions (equations (10.4) and (10.12)).

The partial coalescence regime, which combines both lamellae aggregates and concentration lacunarities, is the subject of current investigations and will be discussed in a subsequent study. Note that, in the stretching and percolating coalescence regimes, the concentration increments PDFs can be derived directly from the concentration PDFs since neighboring concentration levels are respectively fully correlated and uncorrelated. This is not the case for the partial coalescence regime which is characterized by strong spatial intermittency and thus requires an explicit quantification of the spatial correlation properties of the concentration field over all scales.

Coarsening scale

The characteristic correlation scale of the concentration field η may be estimated from the analysis of the scalar variance of the coarse grained concentration fields obtained by averaging the fully resolved concentration field over a coarse grained scale r c [START_REF] Villermaux | Coarse grained scale of turbulent mixtures[END_REF]. For r c < η, the scalar variance is not significantly altered by averaging, while for r c > η, the scalar variance strongly decays with r c . Thus, [START_REF] Villermaux | Coarse grained scale of turbulent mixtures[END_REF] refer to η as the coarse grained scale. It is defined here as the averaging scale for which the variance of the coarse grained concentration is half of the finescale scalar variance. The coarse grained scale η defined this way is shown in figure 30 for σ 2 ln k = 1 and for different Péclet numbers. The estimates of η obtained by this method are consistent with the scale above which the concentration increment PDFs are invariant (figure 25 and 28). For large Péclet numbers η is set by the transverse size of a single lamella. As discussed in section 5, the lamella width evolves diffusively for t > τ m and therefore η ∼ √ t. For small Péclet numbers, the coarsening scale is given by the typical transverse size of lamella bundles [START_REF] Villermaux | Mixing by porous media[END_REF], which are composed of n c (t) diffusively growing lamella, where n c (t) is given by (7.9). Therefore η evolves as η ∼ n c (t)t 1/2 ∼ t d f /2 (10.13) With d f = 1.6, we expect η ∼ t 0.8 , explaining the superdiffusive growth of the coarsening scale observed in figure 30.

Conclusions

A scalar mixture can be viewed as a set of stretched lamellae, possibly interacting through a random aggregation process. We have shown here the relevance of this view- point, which helps deciphering the distinct roles played by advection, disorder, and molecular diffusion in the structure of a mixture progressing in a porous medium. The numerical analysis presented here brings new insights into the precise quantification of the role of stretching and coalescence in the composition of the scalar mixture, through direct measurement of line deformation and statistics of lamella overlaps. To this respect, the use of numerical simulations is unique, since it allows, as opposed to experiments, to track individual material lines even in regions where the concentration field is blurred by molecular diffusion, and therefore to establish without ambiguity the mechanism of random overlap of the lamellae. The scalar mixture is represented as a set of stretched lamellae whose rates of diffusive smoothing build up the overall mixture concentration distribution. Compression perpendicular to the stretching direction enhances diffusive smoothing by steepening the scalar gradient. Hence, the concentration distribution is controlled by the distribution and history of elongations along the line. A detailed analysis of the Lagrangian elongation rate distribution and correlation shows that the latter is characterized by a multistep stretching process with decaying stretching rates. This implies that local line elongations is lognormally distributed, with a mean and variance which evolve with a power law in time. The advected line is not entirely 'space filling' in the sense that it does not cover the dispersion plane homogeneously. The geometry of the advected line is characterized by a constant fractal dimension 1 < d f < 2, which is related to the temporal scaling of the characteristic dispersion length and the average line elongation.

At late times, interaction between lamellae is enforced by confinement within the dispersion area, since the length of the line increases faster than the dispersion length. Lamellae aggregates are formed by diffusive coalescence. The decorrelation of the stretching histories of the lamellae, which form the aggregates, leads to a random aggregation process, and the characteristic concentration of lamellae aggregates is shown to follow a Gamma distribution. As time increases, lamellae aggregates become connected through a process analogous to percolation. This leads to a drastic increase of coalescence numbers in the center of the plume and a wide distribution of overlap numbers over the plume.

This scenario provides analytical expressions for the concentration distributions (both in one point, and for the concentration differences at a given spatial distance) in the different regimes experienced by the scalar mixture as it progresses towards uniformity. The predictions, which are entirely parametrized by the elongation distribution, line geometry and diffusion coefficient, are in good agreement with numerical simulations for all investigated permeability field variances, Péclet numbers and times. These expressions characterize the impact of the permeability field heterogeneity on the temporal evolution of mixing, as quantified by the concentration moments, scalar dissipation rate or entropy. The effective mixing behavior is thus directly related to the mean elongation exponent µ, the elongation variance exponent ν and the line fractal dimension d f . The proposed framework is expected to be particularly relevant for upscaling chemical reactivity in porous media flows since the distribution of concentration gradients resulting from the interaction of stretching and coalescence governs the mixing of chemical elements liable to react. Furthermore, this modeling framework is directly applicable for risk assessment studies, in which uncertainty in transported contaminant concentration is given by the concentration PDFs. While the present study focuses on Darcy scale heterogeneity, the developed theory is expected to be applicable at pore scale, as very similar stretching and coalescence processes are induced by the pore scale flow heterogeneity de Anna et al. (2014). Extension of the theory to 3D porous media is currently under investigation. We anticipate that the additional dimension can be treated by considering stretched sheets instead of lamellae as elementary objects. Furthermore, the stretching behavior is likely to be impacted by the additional degree of freedom.

Figure 1 .

 1 Figure 1. (a) Computed velocity field for a realization of a permeability field with σ 2 f = 1. Corresponding concentration fields c(x, t) for (b) P e = 8 × 10 4 and (c) P e = 8 × 10 2 at times t1 = 12.5τa (red) and t2 = 112.5τa (blue), for an initial injection with a uniform concentration c0 at x1 = λ. The plume is transported from left to right with a mean velocity u and a dispersion length σ(t) increasing with time.

Figure 2 .

 2 Figure 2. Concentration fields c(x, t) at times t1 = 12.5τa (red) and t2 = 112.5τa (blue) for P e = 8 × 10 2 and for different permeability field variances : (a) σ 2 f = 0.25, (b) σ 2 f = 1 and (c) σ 2 f = 4.

Figure 3 .Figure 4 .

 34 Figure3. Advected line at times t1 = 2.5τa, t2 = 37.5τa, t3 = 112.5τa for σ 2 f = 1. The initial position of the line is x1 = λ. The red color intensity is proportional to ρ for lamella elongations ρ 10. The blue color intensity is proportional to 1/ρ for ρ < 10. Hence, the red zones represent the parts of the line that are most elongated while the blue zones represent the portions that have been less stretched.

Figure 7 .

 7 Figure 7. Temporal evolution of the mean elongation ρ(t) for different permeability field variances : σ 2 ln K = 0.25 (green diamonds), σ 2 ln K = 1 (red triangles), σ 2 ln K = 4 (blue disks).

Figure 10 .

 10 Figure 10. Temporal evolution of the average mixing scale (equation (5.2)) a) for different Péclet numbers P e = 1, 10, 10 2 , 10 3 , 10 4 and σ 2 ln K = 1, b) for different permeability field variances σ 2 ln K = 0.25, 1, 4 and P e = 10 4 .

Figure 11 .

 11 Figure11. PDF of mixing scales p m ( m) for P e = 10 4 at times t = 2τa, 4τa, 8τa, 16τa, 32τa, 64τa, 128τa.

Figure 12 .

 12 Figure12. (a) Concentration PDFs at times t = 12.5τa,37.5τa, 62.5τa, 87.5τa, 112.5τa for a Péclet number P e = 8 × 10 4 and a permeability field variance σ 2 f = 1. The simulations are represented by blue continuous lines and the predictions obtained from the independent lamella model (6.2) are represented by black dash dotted lines. Inset: same plot in semi-log representation, emphasizing the non-exponential concentration decay. (b) Concentration PDFs at time 112.5τa for a Péclet number P e = 8 × 10 4 and permeability field variances of σ 2 f = 0.25 (green line), σ 2 f = 1 (blue line) and σ 2 f = 4 (red line). The predictions obtained from the independent lamella model (6.2) are represented by black dash-dotted lines. The inset shows the prediction of the purely diffusive model (3.5), with θ = c0/ √ 4πDt).

Figure 13 .

 13 Figure 13. Concentration field and advected line at time t = 112.5τa and for a Péclet number P e = 800. The zoom in inset illustrates the concentration build up in areas with a large density of lamellae.

Figure 15 .

 15 Figure 15. Illustration of the method used to estimate the number of overlap along the line n[x(t|b), r, t] for a diffusion scale r = 4λ at time t = 112, 5τa.

  .

Figure 16 .

 16 Figure 16. Top row: Map of the number of overlaps n[x(t|b), r, t] along the line for a) r = λ/4, b) r = λ and c) r = 4λ at time t = 112, 5τa. Bottom row: projection of the number of overlap as a function of the longitudinal position n(x1, r, t) for each case.

  Figure 17. a) Projection of the number of overlap as a function of the longitudinal position n(x1, r, t) for r = λ/4 to r = 32λ. b) Comparison of n(x1, r, t)/n(r, t) and ca(x, t) for r = 8λ.

  Figure18. Average number of overlaps n(r, t) for times t1 = 12.5τa (red triangles), t2 = 37.5τa (green diamonds), t3 = 62.5τa (blue squares) and t4 = 112.5τa (black disks). The dashed line represents the prediction of equation (7.9).

Figure 19 .

 19 Figure 19. (a) Concentration PDFs at times t = 12.5τa, 37.5τa, 62.5τa, 87.5τa, 112.5τa for a Péclet number P e = 8 × 10 2 and a permeability field variance σ 2 ln K = 1 (blue continuous lines). The predictions of the coalescence model (equation (7.11)) are shown in black dash dotted lines. The inset displays the same graph in semilog representation, showing the exponential decay of the concentration PDFs in this regime. (b) Concentration PDFs for different permeability variances, σ 2 ln K = 0.25 (green line), σ 2 ln K = 1 (blue line) and σ 2 ln K = 4 (red line), at time 112.5τa for a Péclet number P e = 8 × 10 2 . The predictions of the coalescence model (equation (7.11)) are shown in black dash dotted lines. The inset shows the predictions of the model of non interacting lamellae (equation (6.2)) as black dashed lines. In this regime the coalescence of lamellae tends to increase the probability of occurrence of intermediate concentrations.

Figure 21 .

 21 Figure 20.(a) Concentration PDFs at the center of the plume, x(y) -σx(y)/2 < x < x(y) -σx(y)/2, at times t = 12.5τa, 37.5τa, 62.5τa, 87.5τa, 112.5τa for a Péclet number P e = 8 and a permeability field variance σ 2 ln K = 1 (blue continuous lines). The predictions of the percolating coalescence model (equation (A 2)) are shown in black dash dotted lines. Inset: temporal evolution of the average coalescence number nc(t) corresponding to each

Figure 22 .

 22 Figure 22. Concentration moments c (black lines), ˙c2 ¸(blue lines), scalar dissipation rate χ (red lines) for (a) P e = 8 × 10 3 (b) P e = 80 (c) P e = 8. The results of numerical simulations are shown as continuous lines, while the theoretical predictions are displayed as dashed line. In figure a) the theoretical predictions are obtained from the independent lamellae model (equation (6.2)). In figure b) the theoretical predictions are obtained from the partial coalescence model (equation (7.11)). In figure c) the theoretical predictions are obtained from the percolating coalescence model (equation (A 2)).

Figure 23 .

 23 Figure 23. Temporal evolution of the second moment of concentration ˙c2 ¸for the different permeability field variances, σ 2 ln K = 0.25 (green line), σ 2 ln K = 1 (blue line) and σ 2 ln K = 4 (red line), and for (a) P e = 8× 10 3 , (b) P e = 80 and (c) P e = 8. The dashed lines show the predicted scalings for (a) the independent lamella regime, (8.8), (b) the partial coalescence regime, (8.13), and (c) the percolating coalescence regime, (8.17).

Figure 24 .

 24 Figure 24. Dilution index E(t) (green lines) and segregation intensity I(t) for (a) P e = 8 × 10 3 (b) P e = 80 (c) P e = 8. The results of numerical simulations are shown as continuous lines, while the theoretical predictions are displayed as dashed line. In figure a) the theoretical predictions are obtained from the independent lamellae model (equation (6.2)). In figure b) the theoretical predictions are obtained from the partial coalescence model (equation (7.11)). In figure c) the theoretical predictions are obtained from the percolating coalescence model (equation (A 2)).

Figure 25 .

 25 Figure25. PDF of concentration increments p(∆c) at time t = 112.5τa for a permeability field variance σ 2 ln K = 1 and a Péclet number P e = 8 × 10 3 . The concentration increment PDFs, computed for ∆x = λ/8, λ/4, λ/2, λ, 2λ and 4λ,are rescaled by an arbitrary value for clarity. The predictions obtained from equation (10.3), with η = 0.4λ, show a good agreement for all increments |∆x|.

Figure 27 .

 27 Figure 27. Concentration field and underlying advectively deformed line at time t = 87.5τa for a permeability field variance σ 2 ln K = 1 and a Péclet number P e = 8. The inset illustrates the common and independent lamellae constructing the concentration of two points separated by a distance ∆x. The lamellae contributing to the concentration of a given point are located in a disk of radius r = √ Dt around this point (in red). The common lamellae, which cancel out in the calculation of the concentration increments, are located at the intersection of the two disks of radius r centered on the two points of interest.

Figure 28 .

 28 Figure28. PDF of concentration increments p(∆c) at time t = 87.5τa for a permeability field variance σ 2 ln K = 1 and a Péclet number P e = 8, computed at the center of the plume defined by x(y) -σx(y)/2 < x < x(y) -σx(y)/2. The concentration increment PDFs, computed for ∆x = λ/8, λ/4, λ/2, λ, 2λ, 4λ and 8λ, are rescaled by an arbitrary value for clarity. The predictions obtained from equation (10.9) are shown in dashed lines.

Figure 29 .

 29 Figure29. Structure functions |∆c| q computed as a function of ∆x, for q ranging from 1 to 6, for the increment PDF shown in figure28. The structure functions are normalized by their asymptotic value for large ∆x, |∆c(η)| q . The predictions of equation (10.11) are displayed in black dashed lines.

Figure 30 .

 30 Figure30. Temporal evolution of η for Péclet numbers equal to P e = 80 (triangles), P e = 800 (disks) and P e = 8000 (squares). η is estimated as the coarse graining scale for which the concentration variance is decreased by half.

Acknowledgments

Numerical simulations were performed with the numerical platform H20LAB (http://h2olab.inria.fr). TLB acknowledges the funding of the Marie Curie European Reintegration Grant (ERG) ReactiveFlows ( 230947). MD acknowledges the funding of the European Research Council (ERC) through the project MHetScale (617511).

Appendix A. Global Concentration Statistics in the Percolating Coalescence Regime

A.1. Concentration PDF

In the percolating coalescence regime, t > τ D , the number of coalescences is uniformly distributed according to (7.5). For n > n c (τ D ), this means above the percolation threshold, the lamella density in the aggregate is so large that there is no space for a diffusive profile to develop between lamellae and therefore, the concentration PDF is given by the PDF (7.1) of concentration maxima p m (c, t|n). For n n c (τ D ) the concentration PDF is given by (7.11) as in the partial coalescence regime because diffusive profiles can develop between the lamella aggregates. Hence, the global concentration PDF in the percolation regime may be estimated as the

where n c (t) in this regime is given by (7.10). Expression (A 1) can be written in a more compact form by using exp(-x) n k=0

x k k! = Γ(n+1,x) Γ(n+1) with Γ(n) the Gamma function, and Γ(n, x) the incomplete Gamma-function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF],

A.2. Concentration Moments Thus, we obtain for the concentration moments

where we used (8.12) and (8.15). Similarly, we obtain for the mean squared concentration

(A 4) Using (7.10) for n c (t), we obtain for the mean and mean square concentration c(x, t) ≈ 1 σ(t) ∼ t -β/2 , c(x, t) 2 ≈ 4 3σ(t) 2 ∼ t -β .

(A 5)

The scalar dissipation rate scales accordingly as 6)