Tomás Aquino

Camille Bouchez

Tanguy Le Borgne

Dilution of reactive plumes: evolution of concentration statistics under diffusion and nonlinear reaction

Introduction

Quantifying the variability in solute concentrations in the subsurface is of primary importance to a large range of applications, including reactive transport modeling (e.g., [START_REF] Carl | Reactive transport modeling: An essential tool and a new research approach for the earth sciences[END_REF]), characterizing the natural chemical signature of groundwater systems (e.g., [START_REF] Guadagnini | Probabilistic assessment of spatial heterogeneity of natural background concentrations in large-scale groundwater bodies through functional geostatistics[END_REF]), and assessing human health risks associated with contaminant transport (e.g., [START_REF] Erica | Evaluating effective reaction rates of kinetically driven solutes in large-scale, statistically anisotropic media: Human health risk implications[END_REF]). In particular, when reaction rates depend nonlinearly on available concentrations, effective reaction rates can differ by orders of magnitude from reaction rates determined under well-mixed conditions [START_REF] Battiato | On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media[END_REF][START_REF] Battiato | Applicability regimes for macroscopic models of reactive transport in porous media[END_REF][START_REF] Hubert | Enhanced and non-monotonic effective kinetics of solute pulses under michaelis-menten reactions[END_REF][START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF]. In such cases, knowledge of mean concentrations is not sufficient 1 to predict overall reaction rates, which requires additional knowledge about the spatial variability of concentration values and local reaction rates.

Concentration probability density function (PDF) methods aim to quantify transport and reaction processes in terms of the spatial statistics (i.e., frequency of occurrence of given values) of transported scalars [START_REF] Ee O'brien | The probability density function (pdf) approach to reacting turbulent flows[END_REF][START_REF] Stephen B Pope | Pdf methods for turbulent reactive flows[END_REF][START_REF] Stephen B Pope | Lagrangian pdf methods for turbulent flows[END_REF][START_REF] Pr Van Slooten | Advances in pdf modeling for inhomogeneous turbulent flows[END_REF]. These methods were originally developed in the context of turbulent flows, where they have been extensively employed to model both conservative transport [START_REF] Steven R Hanna | The exponential probability density function and concentration fluctuations in smoke plumes[END_REF][START_REF] Yee | The shape of the probability density function of short-term concentration fluctuations of plumes in the atmospheric boundary layer[END_REF][START_REF] Yee | A simple model for the probability density function of concentration fluctuations in atmospheric plumes[END_REF][START_REF] Alisse | Experimental probability density functions of small-scale fluctuations in the stably stratified atmosphere[END_REF][START_REF] Munro | The high concentration tails of the probability density function of a dispersing scalar in the atmosphere[END_REF][START_REF] Bakosi | Probability density function modeling of scalar mixing from concentrated sources in turbulent channel flow[END_REF] and reactive transport with an emphasis on combustion [START_REF] Dopazo | An approach to the autoignition of a turbulent mixture[END_REF][START_REF] Givi | Probability density function calculations in turbulent chemically reacting round jets, mixing layers and one-dimensional reactors[END_REF][START_REF] Ms | Calculations of premixed turbulent flames by pdf methods[END_REF][START_REF] Rodney | Computation of turbulent reactive flows: first-principles macro/micromixing models using probability density function methods[END_REF][START_REF] Connell | Progress in probability density function methods for turbulent reacting flows[END_REF][START_REF] Broyda | Probability density functions for advective-reactive transport in radial flow[END_REF]. More recently, these methods have gained popularity in the context of porous media to study conservative [START_REF] Guadagnini | Probabilistic assessment of spatial heterogeneity of natural background concentrations in large-scale groundwater bodies through functional geostatistics[END_REF][START_REF] Kapoor | Concentration fluctuations and dilution in aquifers[END_REF][START_REF] Fiori | Concentration fluctuations in transport by groundwater: Comparison between theory and field experiments[END_REF][START_REF] Fiori | Concentration fluctuations in aquifer transport: A rigorous firstorder solution and applications[END_REF][START_REF] Shvidler | Probability density functions for solute transport in random field[END_REF][START_REF] Dentz | Concentration statistics for transport in random media[END_REF][START_REF] Dentz | Probability density functions for passive scalars dispersed in random velocity fields[END_REF][START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Villermaux | Mixing versus stirring[END_REF][START_REF] Icardi | Probability density function (pdf) models for particle transport in porous media[END_REF][START_REF] Hamada | Diffusion limited mixing in confined media[END_REF] and reactive [START_REF] Cirpka | Concentration statistics for mixing-controlled reactive transport in random heterogeneous media[END_REF][START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF][START_REF] De | Anomalous kinetics in diffusion limited reactions linked to non-Gaussian concentration probability distribution function[END_REF][START_REF] Bellin | On the local concentration probability density function of solutes reacting upon mixing[END_REF][START_REF] Cirpka | Probability density function of steady state concentration in two-dimensional heterogeneous porous media[END_REF][START_REF] Chiogna | Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume[END_REF][START_REF] Suciu | A fokker-planck approach for probability distributions of species concentrations transported in heterogeneous media[END_REF] mixing. The statistics of concentrations have been studied in relation to different types of fluctuations. The first approach is to study concentration fluctuations with respect to an average transport model, such as the advectiondispersion equation [START_REF] Kapoor | Concentration fluctuations and dilution in aquifers[END_REF][START_REF] Fiori | Concentration fluctuations in transport by groundwater: Comparison between theory and field experiments[END_REF][START_REF] Fiori | Concentration fluctuations in aquifer transport: A rigorous firstorder solution and applications[END_REF], in terms of a stochastic transport model such as the continuous time random walk [START_REF] Dentz | Concentration statistics for transport in random media[END_REF][START_REF] Sund | Effective models for transport in complex heterogeneous hydrologic systems[END_REF]. Regarding, for instance, average concentrations at a given spatial position, fluctuations typically arise from variability in particle transit times due to heterogeneity. A second approach is to characterize uncertainty in concentrations due to uncertainty in model parameters, such as flow velocities and reaction rates [START_REF] Dentz | Probability density functions for passive scalars dispersed in random velocity fields[END_REF][START_REF] Peter | Stochastic analysis of effective rate constant for heterogeneous reactions[END_REF][START_REF] Tartakovsky | Probability density functions for advectivereactive transport with uncertain reaction rates[END_REF][START_REF] Daniel | Pdf equations for advective-reactive transport in heterogeneous porous media with uncertain properties[END_REF][START_REF] Venturi | Exact pdf equations and closure approximations for advective-reactive transport[END_REF]. These methods quantify the variability in concentrations due to uncertainty in physical and chemical parameters. A third point of view is to consider the variability in concentrations over the spatial extension of a solute plume which arises from known transport and reaction processes. Even with homogeneous physical and chemical parameters, mean concentrations and associated fluctuations evolve under the combined effect of transport and reaction.

In this context, the concentration PDF provides a useful tool to quantify the spatial frequency of occurrence of different solute concentration values, and to provide information about the mixing state of solute plumes without explicitly resolving the corresponding spatial profiles [START_REF] Villermaux | Mixing versus stirring[END_REF].

Here, we follow this last approach to study how the dilution of a reactive plume by diffusive mixing affects concentration statistics (Fig. 1). In this context, evolution equations have been derived for the conservative concentration PDF by quantifying the formation and interactions of lamellar structures in concentration fields subject to advection and diffusion [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Villermaux | Mixing versus stirring[END_REF][START_REF] Le Borgne | Scalar gradients in stirred mixtures and the deconstruction of random fields[END_REF]. Recent developments have extended these theories to mixing-limited bimolecular reactions [START_REF] Guilbert | Chemical reaction for mixing studies[END_REF]. Here, we consider a different reactive transport problem in which mixing does not act to bring segregated reactants into contact but rather to dilute a reactive plume. This system is relevant to a broad range of reactive transport problems [START_REF] Hubert | Enhanced and non-monotonic effective kinetics of solute pulses under michaelis-menten reactions[END_REF][START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], including situations where transported solute plumes react nonlinearly with a solid phase as they get diluted by mixing. This leads to new dynamics that fundamentally differ from what is known regarding mixing-limited reactions across diffusing or dispersing fronts [START_REF] De | The filamentary structure of mixing fronts and its control on reaction kinetics in porous media flows[END_REF][START_REF] Le Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF][START_REF] Jiménez-Martínez | Pore-scale mechanisms for the enhancement of mixing in unsaturated porous media and implications for chemical reactions[END_REF][START_REF] Bandopadhyay | Enhanced reaction kinetics and reactive mixing scale dynamics in mixing fronts under shear flow for arbitrary damköhler numbers[END_REF].

We investigate the coupling of dilution by diffusive mixing with the nonlinear degradation of a single species locally undergoing diffusion in one spatial dimensional and decaying at a rate proportional to an arbitrary power of concentration (Fig. 1). We generalize the results of [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] for the evolution of total mass and mean concentration to the evolution of the full concentration PDF, and we obtain a dynamical equation for the latter in terms of the scalar dissipation rate [START_REF] Yb Zeldovich | The asymptotic law of heat transfer at small velocities in the finite domain problem[END_REF][START_REF] Julio | The kinematics of mixing: stretching, chaos, and transport[END_REF][START_REF] Le Borgne | Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF][START_REF] Barros | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF][START_REF] Nicholas B Engdahl | Scalar dissipation rates in nonconservative transport systems[END_REF], expressed as a function of concentration. In addition to the mean concentration, the concentration PDF formulation allows us to quantify higher-order concentration moments, such as the concentration variance, which provides information about the variability in concentration values observed throughout a spatial domain. Building on the weak-coupling approximation developed in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], we relate the evolution of the concentration PDF and its moments to the evolution of the maximum concentration, and quantify the resulting anomalous kinetics.

In the interest of developing and illustrating the new method proposed here, the class of reactiondiffusion problems we analyze introduces significant simplifications. In particular, we consider a onedimensional problem and neglect heterogeneous advection effects. In addition, the weak-coupling approach involves a Gaussian approximation for the solute plumes at late times. Although reaction under diffusion or dispersion in the absence of advective heterogeneity has been studied in its own right in the context of porous media [START_REF] Battiato | On breakdown of macroscopic models of mixing-controlled heterogeneous reactions in porous media[END_REF][START_REF] Deepashree | Experimental study of bimolecular reaction kinetics in porous media[END_REF][START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF][START_REF] Valdes-Parada | On the effective diffusivity under chemical reaction in porous media[END_REF][START_REF] Fj Valdés-Parada | On diffusion, dispersion and reaction in porous media[END_REF][START_REF] Pereira | Multi-scale modeling of diffusion and reaction-diffusion phenomena in catalytic porous layers: comparison with the 1d approach[END_REF], these conditions may be seen as fundamental building blocks in the theory of mixing in heterogeneous flows. In the presence of heterogeneous flow fields characteristic of porous media, solute plumes are deformed into a collection of elongated structures termed lamellae (Fig. 1) [START_REF] Villermaux | Mixing versus stirring[END_REF]. Using the Ranz transform [START_REF] William E Ranz | Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows[END_REF], the resulting compression-diffusionreaction equation can be mapped onto a one-dimensional diffusion-reaction equation, which motivates m X k R J k q H W L X 5 1 e T E x E h S Y c K 9 X 2 3 E T 7 K Z a a E U 4 n h Y 5 R N M F k i P u 0 b a n A E V V + O r t 6 g k 6 s 0 k N h L G 0 J j W b q 7 4 k U R 0 q N o 8 B 2 R l g P 1 K I 3 F f / z 2 k a H l 3 7 K R G I 0 F W S + K D Q c 6 R h N I 0 A 9 J i n R f G w J J p L Z W x E Z Y I m J t k E V b A j e 4 s v L p F E p e 2 f l 8 7 t K q X q V x Z G H I z i G U / D g A q p w C z W o A w E J z / A K b 8 6 j 8 + K 8 O x / z 1 p y T z R z C H z i f P 5 I 1 k p A = < / l a t e x i t > Figure 1: Illustration of the reactive transport problem under consideration (adapted from a simulation of transport at the Darcy scale in a heterogeneous permeability field [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF]). We consider a plume of reactive solute (green) that undergoes nonlinear decay either by reaction with the solid, with another solute in excess everywhere, or with another solute that is well mixed inside the plume. The plume is diluted by mixing with the non-reactive resident fluid (black) while it reacts, leading to non-trivial effective reaction rates [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF]. The latter are governed by the coupling of nonlinear reaction with the diffusive flux along concentration gradients. Here, we study the one-dimensional diffusion-reaction problem, which represents the first step towards solving the full advection-diffusion-reaction system.

L p d H f 9 T G r Q P a f F 3 W t M F F 1 q U = " > A A A C K X i c b V B N S w M x E M 3 6 b f 2 q e v Q S L I I X y 6 4 o e i z q w a O C t Y V u K d n s b D u Y T Z Y k W y 3 F v + P F v + J F Q V G v / h H T j 4 N a H w Q e 7 8 1 M Z l 6 U C T T W 9 z + 8 q e m Z 2 b n 5 h c X C 0 v L K 6 l p x f e P a q F x z q H I l l K 5 H z I B A C V W L V k A 9 0 8 D S S E A t u j k d + L U u a I N K X t l e B s 2 U t S U m y J l 1 U q t Y C a V C G Y O 0 9 A x F P h B p 1 K M p 3 q F s 0 1 u 0 H S q V 3 A v D g p v K L X a B a j A 4 b E h E j n G r W P L L / h B 0 k g R j U i J j X L S K L 2 G s e J 6 6 C V w w Y x q B n 9 l m n 2 m L X M B 9 I c w N Z I z f s D Y 0 H J U s B d P s D y + 9 p z t O i W m i t H t u g 6 H 6 s 6 P P U m N 6 a e Q q U 2 Y 7 5 q 8 3 E P / z G r l N j p t 9 l F l u Q f L R R 0 k u q F V 0 E B u N U Q O 3 o u c I 4 x r d r p R 3 m H a R u H A L L o T g 7 8 m T 5 H q / H B y U D y / 3 S 5 W T c R w L Z I t s k 1 0 S k C N S I e f k g l Q J J w / k i b y S N + / R e / b e v c 9 R 6 Z Q 3 7 t k k v + B 9 f Q M 0 G q c 2 < / l a t e x i t >
K W 8 t E E d o i k k v V D b C m n A n a M s x w 2 o 0 V x V H A a S e Y 3 G R + 5 5 E q z a R 4 M N O Y + h E e C R Y y g o 2 V B u X K n R R Z F i t k Y 2 Q h V t 2 a O w d a J V 5 O q p C j O S h / 9 Y e S J B E V h n C s d c 9 z Y + O n W B l G O J 2 V + o m m M S Y T P K I 9 S w W O q P b T + e 0 z d G q V I Q q l s k 8 Y N F d / J 1 I c a T 2 N A j s Z Y T P W y 1 4 m / u f 1 E h N e + S k T c W K o I I t F Y c K R k S g r A g 2 Z o s T w q S W Y K G Z v R W S M l e 3 A 1 l W y J X j L X 1
the present study.

The paper is structured as follows. Section 2 introduces a general definition for the concentration PDF associated with a given spatial profile and explores its key properties. In Section 3, in light of these concepts, we derive a governing equation for the dynamics of the concentration PDF of a solute undergoing diffusion-induced dilution in one spatial dimension and nonlinear reaction described locally by a power-law decay of concentration. In this context, we also discuss the relationship between this simplified one-dimensional problem and reaction-diffusion in the presence of heterogeneous flows typically found in porous media, based on the Ranz transform. We then study the one-dimensional problem and the associated effective kinetics in terms of the weak-coupling approximation in Section 4.

Conclusions and an overall discussion of the results are presented in Section 5. Some additional technical derivations are provided in the Appendices.

The concentration PDF

This section is concerned with providing a general definition of the concentration PDF associated with arbitrary spatial concentration profiles. We first present the general multispecies formulation, and then discuss the single-species case. We provide a detailed discussion of the properties of the concentration PDF, as well as of its numerical computation. The approach developed here, which forms the basis for the subsequent sections, builds on that introduced in [START_REF] Aquino | The diffusing-velocity random walk: a spatial-Markov formulation of heterogeneous advection and diffusion[END_REF] regarding Eulerian velocity PDFs. The main concepts behind the concentration PDF, which quantifies the spatial frequency of occurrence of different concentration values, are illustrated in Fig. 2 for a one-dimensional, single-species concentration profile.

Before proceeding, we introduce some notation that will be employed throughout. The concentration C(x; t) at a given location x and time t, with C = (C 1 , . . . , C ns ), is defined so that C(x; t) dx is the mass in an infinitesimal volume dx centered at x = (x 1 , . . . , x d ) of each chemical species 1, . . . , n s .

We use an underscore to denote quantities whose components relate to concentration species (of which there are a number n s ), and boldface to denote spatial vectors (which have components according to vector a denotes its Euclidean norm. We consider a spatial domain Ω(t) where the concentrations are defined for each time t, and denote by X(t) a random variable with a uniform distribution over Ω(t).

The average with respect to X(t), i.e., the spatial average, is denoted by • . We also introduce, for each time t, a random variable C(t) taking concentration values according to their spatial probability of occurrence, and denote the average with respect to it (i.e., over concentrations) by • . As will be discussed in more detail below, the two averages, while conceptually different, yield equivalent results.

That is, for an arbitrary function f (c, t) of concentrations c and times t,

f (C[X(t), t]) = 1 |Ω(t)| Ω(t) dx f (C[X(t), t]) = ∞ 0 dc f (c, t)p(c; t) = f [C(t), t] , (1)
where, at each time t, p(•; t) is the concentration PDF, to be formally defined momentarily.

General formulation for a multispecies system

At each time t, the concentration PDF p(•; t) is defined so that p(c; t) dc is the probability, in the sense of spatial frequency of occurrence, of finding a concentration in the infinitesimal vicinity dc = ns j=1 dc j of c in the domain Ω(t). According to this definition, we have

p(c; t) = δ[c -C[X(t); t]] = 1 |Ω(t)| Ω(t) dx δ[c -C(x; t)], (2)
where δ(•) is the Dirac delta. Intuitively, this expression encodes the fact that the contributions to the probability of having a value c of concentration at time t correspond to points x where C(x; t) = c, which are "counted" by the Dirac delta in accordance with the continuous nature of concentrations (for the general theory of PDFs as averages of Dirac deltas, see, e.g., [START_REF] Feller | An introduction to probability theory and its applications[END_REF]). This corresponds to the continuum limit described in Fig. 2.

The concentration PDF contains information about the full structure of point statistical variability of concentrations. According to the definition, Eq. ( 2), moments of concentration, such as the component-wise mean concentrations C(t) and variances σ c 2 = C 2 (t) -C(t) 2 (representing not a spatial variance but variability with respect to the mean concentration), may be computed given knowledge of the PDF, and agree with those obtained from the spatial profile. In particular,

C(t) = ∞ 0 dc cp(c; t) = M (t) |Ω(t)| = C[X(t); t], (3)
where the total component masses are given by

M (t) = Ω(t) dx C(x; t). (4)
It is important to note that Ω(t) may be time-dependent, but its size |Ω(t)| must be finite at all times, because a uniformly-distributed random variable X(t) does not exist in an infinite domain.

Intuitively, this is related to the fact that the concentration PDF is defined in terms of the spatial frequency of occurrence of concentration values. Thus, if the total component masses are finite, most regions in an infinite domain Ω(t) must have a vanishing concentration, and p(c; t) = 0 for all c = 0, so that we must have trivially p(c; t) = δ(c) by normalization. In practice, as will be discussed in more detail below, this means that in order to define a non-trivial concentration PDF we must consider either a finite, fixed reference domain Ω, or a minimum concentration threshold c m such that Ω(t) is the region where concentrations are above this threshold. Finally, we note that, in an infinite periodic system, the concentration PDF may be naturally defined over a single period. The previous argument can be circumvented in this case because the total mass in such an infinite system is infinite, although the mass per period is finite.

Single-species formulation

Writing the concentration PDF in a more useful analytical form requires changing variables in the Dirac delta in Eq. ( 2), in order to identify the spatial points x that contribute to the integral for a given value of concentration c. In other words, we must invert the relations C(x; t) = c, which form a set of n s constraints in d dimensions. These constraints are not necessarily independent, and this is a difficult problem in general from an analytical perspective. We now derive a simple analytical form for the single-species case. We write C = C 1 and c = c 1 for the single species component. A numerical computation procedure, also valid for the general multispecies PDF, is discussed in Appendix A.

Assuming C(•; t) to be smooth, the Dirac delta can be expressed as a simple-layer integral [START_REF] Hörmander | The analysis of linear partial differential operators I: Distribution theory and Fourier analysis[END_REF],

δ[c -C(x; t)] = Λ(c;t) dS(y) δ(x -y) |∇C(y; t)| , (5)
where dS(x) is the (d -1)-area element at point x on the (d -1)-surface consisting of the points in the domain at time t at which the concentration has value c, and which do not constitute a local extremum (or zero-gradient inflection point). Formally, this surface is defined as Λ(c; t) = {x ∈ Ω(t) | C(x; t) = c ∧ |∇C(x; t)| = 0}. Note that, for the one-dimensional example in Fig. 2, Λ(c) = 2 for all concentrations strictly between the minimum and maximum values. This is discussed in more detail in Section 2.4. Substituting Eq. ( 5) in Eq. ( 2), we find

p(c; t) = Λ(c;t) dS(x) |Ω(t)||∇C(x; t)| . (6)
This form corresponds to the intuitive notion of a change of variables mentioned above: it expresses the contribution to the probability density of a given value c of concentration as a spatial integral over the spatial region where the concentration takes the required value (see also Fig. 2). This spatial region has dimensionality one less than the full space due to the constraint C(x; t) = c. This result holds for smooth, non-constant C(•; t). The contributions of concentrations found near spatial concentration extrema and in extended regions of constant concentration are discussed below.

Equation [START_REF] Hubert | Enhanced and non-monotonic effective kinetics of solute pulses under michaelis-menten reactions[END_REF] can also be expressed so as to highlight the role of concentration gradients associated with a given concentration in determining the value of the PDF. To this end, consider the harmonic average g h (c; t) of the concentration gradient magnitudes in Λ(c; t). By definition, it is given by

g h (c; t) =    1 |Λ(c; t)| Λ(c;t) dS(x) |∇C(x; t)|    -1 . (7)
Substituting this definition in Eq. ( 6) immediately yields the simple form

p(c; t) = |Λ(c; t)| |Ω(t)|g h (c; t) . (8)
The proportionality to the inverse of a concentration gradient magnitude is a consequence of the fact that p(c; t) is a density, i.e., a probability per unit concentration. Spatial regions where the concentration gradient is low correspond to a higher density around the corresponding concentration value due to the slow variation. To understand this, consider a surface Λ(c; t) over which the gradient magnitude g h (c; t) = |∇C| is homogeneous and therefore equal to its harmonic average. A small variation dc of concentration occurs over a distance dz = dc/g h (c; t) along the concentration gradient.

The gradient is locally perpendicular to the surface of constant concentration, so that the spatial This is illustrated for a one-dimensional profile in Fig. 2. More generally, if the gradient is not constant over the surface, the same argument can be applied locally, using the area element dS introduced above and the volume element dSdz. The inverse proportionality to the harmonic average then arises due to the inverse proportionality of the local contributions to the inverse gradient magnitude, which must be integrated over. This provides an intuitive explanation for Eq. ( 6) and recovers Eq. (8).

Extrema and constant-concentration regions

We now study the impact on the concentration PDF of local extrema of the spatial concentration profile, at which |∇C(x; t)| = 0. Near an extremum c 0 at position x 0 and time t, we have C(x 0 + ∆x; t) ≈ c 0 + ∆x • H(x 0 ; t)∆x/2 and ∇C(x; t) ≈ H(x 0 ; t)∆x, where H(x 0 ; t) is the Hessian matrix of concentration at position x 0 and time t, with components H ij (x 0 ; t) = ∇ i ∇ j C(x; t)| x=x0 . We find that in d = 1 the PDF of concentration behaves in a vicinity ∆c of c 0 as

p(c 0 + ∆c; t) ≈ |Ω(t)| ∆c ∂ 2 C(x; t) ∂x 2 x=x0 -1 , (9)
to first order in |∆c| c 0 , for ∆c < 0 near a maximum and ∆c > 0 near a minimum. This type of divergent behavior due to a spatial maximum can be observed in Fig. 2. We note that this is an integrable divergence, that is, the integral of p(•; t) converges near c 0 as necessary for a PDF. This result is a direct consequence of the fact that for a smooth function near an extremum the behavior of the function is quadratic, and the behavior of its derivative is linear. We have assumed the extremum is non-degenerate, so that the second derivative is nonzero; otherwise, higher order terms are important and the behavior is different.

In arbitrary dimension, we provide an argument based on dimensional analysis. First, note that

p(c 0 + ∆c; t) ∝ |∆c| d/2-1 |Ω(t)||H| 1/2 . ( 10
)
Thus, there is no divergence in d 2. In d = 2, extrema of the spatial concentration correspond to extrema of the concentration PDF, and in d = 3 they correspond to zeros. As before, we have assumed that the minimum is non-degenerate, i.e., |H| = 0. As a simple example of degeneracy, consider a concentration field which is constant along the second and third dimensions, and exhibits a regular extremum along the first dimension. As expected, the divergence as the extremum is neared then depends only on the distance ∆x along the first dimension, and the square-root divergence of the concentration PDF characteristic of d = 1 is recovered.

If C(•; t) is constant in space (i.e., the concentration is homogeneous), or if we relax the smoothness assumption so that |∇C(•; t)| may be zero in open sets by letting C(•; t) be non-differentiable at their boundary (i.e, piecewise-smooth), these sets may contribute additional terms to the concentration PDF.

Consider first that the concentration is piecewise-constant. Each set Ω i (t) where the concentration takes a given position-independent value c i (t) at time t is associated with a probability

|Ω i (t)|/|Ω(t)|.
Note that there must be a finite number N (t) of such constant-concentration sets at each time, and thus a finite number of corresponding concentration values, because open sets have nonzero measure and sets corresponding to different concentrations are disjoint, so that

|Ω(t)| = N (t) i=1 |Ω i (t)|. Therefore, p(c; t) = N (t) i=1 |Ω i (t)| |Ω(t)| δ[c -c i (t)]. (11)
The Dirac deltas express the fact that concentration takes a single specific value within spatiallyhomogeneous regions. If the concentration is piecewise-smooth and non-constant over certain parts of the domain, Eq. ( 8) may be applied separately to each such region, and the contributions of all constant and non-constant regions summed over to obtain the full concentration PDF.

One-dimensional Gaussian profile

As a relevant example which will play a central role in Section 3, we consider a single-species Gaussian concentration profile in d = 1 spatial dimension (Fig. 1), characterized by spatial variance σ 2 (t) and total mass M (t):

C(x; t) = M (t) 2πσ 2 (t) exp - x 2 2σ 2 (t) . ( 12
)
The corresponding regions of constant concentration c are zero-dimensional, Λ(c; t) = {x|x = ±x c (c; t)}.

Inverting C(x; t) = c for x, x c (c; t) = 2σ 2 (t) ln[c M (t)/c], (13)
where

c M (t) = M (t) 2πσ 2 (t) (14)
is the maximum concentration at time t. The zero-dimensional area element on Λ(c; t) is a sum of point measures concentrated at the two points ±x c (c; t) where the concentration takes a given value,

dS(x) = (δ[x + x c (c; t)] + δ[x -x c (c; t)]) dx. For all c < c M (t), |Λ(c; t)| = 2 because it consists of
two points. The gradient of concentration magnitudes at ±x c (c; t) have the same value, so that their harmonic average is equal to this value and given by

g h (c; t) = |∇C(x; t)| x=±xc(c;t) = 2c 2 ln[c M (t)/c] σ 2 (t) . ( 15
)
As discussed in Section 2, we must consider a finite domain in order for the concentration PDF to be well-defined. This can be conveniently achieved through one of two possible approaches. The first considers a fixed (i.e., time-independent) reference domain Ω centered at x = 0. In this case,

a time-dependent minimum concentration c m (t) = c M (t) exp{-|Ω| 2 /[8σ 2 (t)]
} is found at the fixed domain boundaries. Alternatively, we fix a (time-independent) minimum concentration c m , and consider Ω(t) = {x|C(x; t) > c m }, the time-dependent spatial domain where concentration are above the minimum value. In this case, c m can be thought of as a lower detection limit, and the concentration PDF as being computed over the region of space where non-zero concentration is detected. In either case, using Eq. ( 8) leads to the concentration PDF

p(c; t) = H[c M (t) -c]H[c -c m (t)] 2c ln[c M (t)/c] ln[c M (t)/c m (t)] , (16)
where under the second approach c m (t) = c m is constant. This PDF is illustrated in Fig. 2, and it agrees with the result reported, e.g., in [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF]. By noting that |Ω(t)| = 2x c [c m (t); t] and using Eqs. ( 13) and ( 14), we find

|Ω(t)| = 2M (t) c M (t) ln[c M (t)/c m (t)] π , (17)
and the concentration PDF can thus also be expressed as

p(c; t) = M (t)H[c M (t) -c]H[c -c m (t)] |Ω(t)|cc M (t) π ln[c M (t)/c] . (18)
Note that, for fixed t, setting c = c M (t) -∆c and expanding to first order in ∆c recovers the expected ∆c -1/2 divergence due to the presence of a spatial maximum discussed in Section 2.3. For small concentrations compared to the peak value, this PDF scales as 1/c with logarithmic corrections.

Using these results and Eq. ( 12), it can be verified by direct computation that the mean concentration is given by

C(t) = C(X; t) ≈ M (t) |Ω(t)| = c M (t) 2 π ln[c M (t)/c m (t)] , (19)
and for the second moment of concentration we have

C 2 (t) = C 2 (X; t) ≈ M (t)c M (t) √ 2|Ω(t)| = c 2 M (t) 2 π 2 ln[c M (t)/c m (t)] , (20)
where the approximation corresponds to approximating the total mass M (t) in the infinite domain by the mass within Ω(t), over which concentration is larger than c m (t). Note that, in particular, this approximation requires c M (t) c m (t), that is, the lower detection limit must be small compared to the peak concentration. Under the same approximation, the variance of concentration is thus

σ 2 c (t) ≈ πc 2 M (t) 4 ln[c M (t)/c m (t)] 2 ln[c M (t)/c m (t)] π -1 , (21)
or equivalently

σ 2 c (t) C(t) 2 ≈ 2 ln[c M (t)/c m (t)] π -1. ( 22
)
Note that the mean and variance of concentration represent the mean value and the associated variability in the frequency of occurrence of concentration values found throughout the domain. They differ from the more common metrics of mean plume position, which here is zero, and spatial or plume variance, which here is σ 2 (t).

Concentration PDF evolution under nonlinear reaction and diffusion

In this section, we apply the concepts developed above to the single-species, one-dimensional reactiondiffusion problem

∂C(x; t) ∂t = D ∂ 2 C(x; t) ∂x 2 -r[C(x; t)], (23)
where the reaction rate as a function of concentration has the form

r(c) = κc β , (24)
with β > 0 the order of the reaction and κ the reaction rate constant in units of [C 1-β T -1 ]. Understanding the role and dynamics of concentration statistics under one-dimensional reaction-diffusion represents a fundamental first step towards treating more complex reaction problems in heterogeneous flow fields. In this connection, we discuss the Ranz transform, which was developed to quantify mixing and reaction in the context of turbulence and combustion problems [START_REF] William E Ranz | Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows[END_REF][START_REF] Frank | The coherent flame model for turbulent chemical reactions[END_REF] and has more recently been employed to study these phenomena in porous media [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Le Borgne | Stretching, coalescence, and mixing in porous media[END_REF]. The goal of the Ranz transform is to exploit the fact that material elements in an incompressible flow tend to orient themselves according to the principal directions of strain-induced stretching (Fig. 1), forming thin lamellar structures (lines in 2D or sheets in 3D) [START_REF] De | The filamentary structure of mixing fronts and its control on reaction kinetics in porous media flows[END_REF][START_REF] Keith | The effect of homogeneous turbulence on material lines and surfaces[END_REF][START_REF] Betchov | An inequality concerning the production of vorticity in isotropic turbulence[END_REF][START_REF] Wj Cocke | Turbulent hydrodynamic line stretching: consequences of isotropy[END_REF][START_REF] Wm T Ashurst | Alignment of vorticity and scalar gradient with strain rate in simulated navier-stokes turbulence[END_REF][START_REF] Ss Girimaji | Material-element deformation in isotropic turbulence[END_REF][START_REF] Martínez-Ruiz | The diffusive sheet method for scalar mixing[END_REF][START_REF] Martínez-Ruiz | The diffusive sheet method for scalar mixing[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF][START_REF] Heyman | Stretching and folding sustain microscale chemical gradients in porous media[END_REF].

Due to incompressibility, the thickness of such a material element along the principal axis of compression is given by s(t) = s(0)A(0)/A(t), where A(t) is its area (or length in 2D). Assuming that local compression is well described by a first-order Taylor expansion of the flow, and neglecting concentration gradients on the surface in comparison to normal gradients along the compression direction ζ, the advection-diffusion-reaction equation becomes [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Martínez-Ruiz | The diffusive sheet method for scalar mixing[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF]]

∂C ∂t = - ζ s ds dt ∂C ∂ζ + D ∂ 2 C ∂ζ 2 -r(C). (25)
These approximations provide a very good description of transport at high Péclet number [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Martínez-Ruiz | The diffusive sheet method for scalar mixing[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF].

Applying the Ranz transform [START_REF] William E Ranz | Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows[END_REF],

τ = t 0 dt s 2 (0) s 2 (t ) , ξ = s(0) s(t) ζ, (26)
reduces Eq. ( 25) to a diffusion-reaction equation in the new time τ and rescaled transverse coordinate ξ, which takes into account the history of advective deformation:

∂C ∂t = D ∂ 2 C ∂ξ 2 - s 2 s 2 0 r(C). (27)
Thus, the full advection-diffusion-reaction problem reduces to a collection of one-dimensional diffusionreaction problems over lamellae. Solving this problem remains highly nontrivial, as it requires determining the statistics of thicknesses s(t) and those of lamella coalescence (overlapping) [START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Le Borgne | Scalar gradients in stirred mixtures and the deconstruction of random fields[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Le Borgne | Stretching, coalescence, and mixing in porous media[END_REF][START_REF] Duplat | Mixing by random stirring in confined mixtures[END_REF][START_REF] Daniel R Lester | Chaotic mixing in three-dimensional porous media[END_REF], and the solution of the one-dimensional reaction-diffusion problem [START_REF] Shvidler | Probability density functions for solute transport in random field[END_REF] for given s(τ ), where the reaction rate becomes time-dependent according to a rescaling by s 2 (τ )/s 2 0 . Treating the full problem is beyond the scope of this work. Here, we restrict ourselves to the purely-diffusive case, Eq. ( 23), which corresponds to setting s 2 /s 0 = 1 in Eq. ( 27), with the power-law reaction rate [START_REF] Kapoor | Concentration fluctuations and dilution in aquifers[END_REF].

The dynamics remain nontrivial even under this simplified scenario. When the reaction is nonlinear (order β nonzero and different from unity), the effective kinetics, representing the overall reaction rate for a given amount of available mass, depend on the spatial variability of concentration. In other words, the global reaction rate dM (t)/dt for the total mass is not given by κM β (t), except when β = 1 or when the system is well mixed, i.e., uniformly characterized by a single concentration value.

Similar considerations apply to the reaction rate in terms of the mean concentration instead of the total mass. Nonetheless, knowledge of the point statistics of concentration as encoded in the concentration PDF is sufficient to fully characterize the global reaction rate, which can be obtained by integrating Eq. ( 23) over the domain Ω(t). Using integration by parts, the diffusive term can be seen to contribute 2D|∂C/∂x| evaluated at the domain boundaries, which we neglect under the assumption of a small lower detection limit c m as above, we obtain

dM (t) dt = - ∞ 0 dc r(c)p(c; t) = -κ ∞ 0 dc c β p(c; t), (28)
where in the last equality we have expressed the integral of the reactive contribution in terms of the concentration PDF using its definition, Eq. ( 2).

In what follows, we will first derive a dynamical equation for the evolution of the concentration PDF, in terms of a concentration-dependent scalar dissipation rate. We will then discuss and employ the Gaussian spatial profile arising from the weak-coupling approximation of [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] to directly compute the evolution of the concentration PDF using the results of Section 2.4, and compare it to the PDF obtained from numerical solutions of Eq. ( 23).

Dynamical equation for the concentration PDF

The concentration PDF evolves under the combined action of transport and reaction processes. The structure of the reaction-diffusion equation [START_REF] Broyda | Probability density functions for advective-reactive transport in radial flow[END_REF], where transport and reaction terms locally contribute independently to the evolution of the concentration field, ensures that the same holds for the evolution of the concentration PDF for each value of concentration. We thus seek to express its evolution as

∂p(c; t) ∂t = ∆p D (c; t) + ∆p R (c; t), (29)
where the two contributions respectively represent the effect of diffusion and reaction. Note that, while this independence holds locally, the global reaction dynamics (i.e., the evolution of the total mass and other global quantities such as the concentration mean and variance) are affected by the dynamical evolution of the plume through the combined effect of transport and reaction. These dynamics will be discussed in detail in Section 4 using a weak-coupling approximation. We focus in what follows on the case of a fixed minimum detection limit c m and a time-varying domain Ω(t) where concentrations are above c m . This choice avoids the introduction of an arbitrary domain devoid of physical meaning, and corresponds more closely to real-world situations, where a lower detection limit always exists.

Nonetheless, note that the introduction of a fixed reference domain may be appropriate if one is in fact interested in a specific spatial region, and the formalism can be easily adapted to that case.

As shown in Appendix B, the reactive contribution in Eq. ( 29) is given by It should be clear that, in general, knowledge of the frequencies of occurrence of concentration values, as expressed by p(c; t), is not sufficient to determine their evolution under transport. This is because concentration fluxes due to advection and diffusion are sensitive to the spatial gradients in the concentration profiles. Thus, without additional information, an equation for the concentration PDF, which encodes only the point statistics of concentration, cannot be closed. However, it is possible to express the concentration PDF's evolution in terms of an additional mixing metric, which corresponds to the scalar dissipation rate (see, e.g., [START_REF] Julio | The kinematics of mixing: stretching, chaos, and transport[END_REF][START_REF] Le Borgne | Non-Fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media[END_REF]),

χ x (x; t) = (∇C) • (D∇C), (31)
expressed as a function of concentration. We make the simplifying assumption that the initial concentration profile C(x, 0) is symmetric about the origin and monotonically decreasing away from it (as happens, e.g., for a Gaussian profile or a square pulse). In that case, these properties hold for all times, so that a given concentration value is associated with a single concentration gradient magnitude, simplifying the description. Under these assumptions, the concentration-dependent scalar dissipation rate is given by

χ(c; t) = χ x [±x c (c; t); t] = Dg 2 h (c; t), (32)
where the harmonic average gradient magnitude g h is equal to the gradient magnitude at both points ±x c where the concentration is equal to c. As shown in Appendix B, we have for the transport contribution:

∆p D (c; t) = p(c; t) 1 4 ∂ ln χ(c; t) ∂c ∂χ(c; t) ∂c - 1 2
∂ 2 χ(c; t) ∂c 2 - p(c m ; t) 2 ∂χ(c; t) ∂c c=cm , (33)
where the first and second terms arise from the impact of diffusion on the concentration gradients associated with a given concentration value, and the third term is due to the change in domain size due to concentrations dropping below c m by diffusion.

Substituting Eqs. ( 30) and [START_REF] Hamada | Diffusion limited mixing in confined media[END_REF] in Eq. ( 29), we conclude that the combined effect of transport and reaction leads to

∂p(c; t) ∂t = ∂r(c)p(c; t) ∂c + 1 4 ∂ ln χ(c; t) ∂c ∂χ(c; t) ∂c - 1 2
∂ 2 χ(c; t) ∂c 2 + r(c m ) - 1 2 ∂χ(c; t) ∂c c=cm p(c m ; t) p(c; t). (34)
This equation shows that the dynamical evolution of the concentration PDF for this problem is fully As discussed in the introduction, the scalar dissipation rate has been studied in connection with mixing in porous media. By clarifying its role in the evolution of the concentration PDF, this approach provides a promising departure point for upscaling of reaction dynamics under transport.

Concentration dynamics under the weak-coupling approximation

The weak-coupling approximation developed in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] simplifies the full reaction-diffusion dynamics of Eq. ( 23) by assuming that, after a time on the order of the diffusion time, the concentration profile remains Gaussian as for the conservative problem, but with a spatial variance and mass that are affected by reaction. We briefly review this approach before applying the theory and concepts developed here to analyze the evolution of concentration statistics.

Problem setup and nondimensionalization

We will focus here, as in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], on a rectangular pulse initial condition, corresponding to constant concentration c 0 = M 0 /s 0 over a domain Ω(t) centered at the origin and of initial width s 0 = |Ω(0)|. This setup correspond to a Dirac delta initial condition for the PDF, p(c; 0) = δ(c-c 0 ). The initial condition allows for meaningful direct comparison with the well-mixed batch problem in the domain Ω(0), which corresponds to the purely-reactive scenario in the absence of diffusion. In that case, D = 0, so that the scalar dissipation rate χ ≡ 0, and we denote the associated uniform solution for concentration as a function of time as c B (t).

In addition, in order to meaningfully compare the behavior of the system across different parameter values, it is convenient to nondimensionalize concentration according to the initial concentration c 0 , position according to the initial width s 0 , and time according to the characteristic reaction time [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] τ R = κ -1 c 1-β 0 . Thus, denoting nondimensionalized quantities by an asterisk, nondimensional time is given by t * /τ R , positions by x * = x/s 0 , and concentrations by C * = C/c 0 , seen as a function of nondimensional times and positions. Functions of these quantities, such as the concentration PDF, are nondimensionalized accordingly; further details may be found in Appendix D. In what follows, we work exclusively in nondimensional units and drop the asterisks for notational brevity. The nondimensional transport equation is given by

∂C(x; t) ∂t = 1 2 Da ∂ 2 C(x; t) ∂x 2 -r[C(x, t)], (35)
where we have introduced the Damköhler number

Da = τ D τ R , τ D = s 2 0 2D . (36)
A large Damköhler number corresponds to fast reaction compared to the diffusion time τ D associated with the nondimensionalization lengthscale s 0 , and a low Damköhler number to slow reaction. Note that Eq. ( 34) for the evolution of the concentration PDF remains unchanged when all intervening quantities are nondimensionalized.

Early times -Batch dynamics

Under the weak-coupling approximation, up to times on the order of the diffusion time, the concentration distribution is approximated as not having deformed substantially with respect to the initial condition, here a square pulse of unit width in nondimensional coordinates. Thus, up to nondimensional time t = Da, the system is taken to evolve according to reaction in a batch reactor. Since variations in the spatial distribution up to this time are considered negligible and not resolved, under this approximation the concentration PDF remains a Dirac delta for t Da, p(c; t

) ≈ δ[c -c B (t)].
We have (see also Appendix E):

dc B dt = -c β B . (37)
The special case β = 1 corresponds to linear reactions. Decay is then exponential, c B (t) = exp(-t).

In this classical case, reaction and transport dynamics are fully decoupled: the exact solution of the full reactive transport problem is simply the conservative transport solution, multiplied by exp(-t) to account for mass loss due to reaction. For nonlinear reactions, β = 1, integrating Eq. ( 37) yields

c B (t) = [1 -(1 -β)t] 1 1-β . (38)
In nondimensional units, the total mass and concentration coincide for batch dynamics, because the width of the pulse is approximated as constant and equal to the initial unit width in this time window. Thus, for t < Da we approximate M (t) ≈ c B (t). Note that for β < 1 the batch reaction is fully depleted in a finite time

t B f = 1 1 -β . ( 39
)
If t B f < Da, the full dynamics are thus simply treated as a well-mixed batch reactor.

Late times -Gaussian profile

For times t < Da, [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] have verified that the dynamics are reasonably well approximated by a well-mixed batch reactor, disregarding transport and spatial variability in concentration. We thus focus now on times t Da, for which spatial concentration variability cannot be neglected. Starting at t = Da, when diffusion has had sufficient time to deform the concentration pulse, the weak-coupling approximation consists in approximating the spatial profile as a Gaussian, whose properties are governed by the combined effect of reaction and diffusion as detailed below. Consider first Eq. ( 34) for the dynamical evolution of the concentration PDF. Under the Gaussian approximation, the scalar dissipation rate can be computed analytically. From Eqs. ( 12) and ( 13) for the Gaussian profile and (77b) for the scalar dissipation rate, we find

χ(c; t) = 2πc 2 M (t)c 2 ln[c M (t)/c] Da M 2 (t) H[c M (t) -c]H[c -c m ], (40)
where we have used Eq. ( 14) to express the spatial variance in terms of the mass and peak concentration.

Note that, expanding the logarithm in Taylor series for small ∆c near the peak value, c = c M (t) -∆c, we find that the scalar dissipation rate approaches zero linearly with ∆c. At low concentrations compared to the peak value, it scales as c 2 with logarithmic corrections. Using Eq. ( 40), Eq. ( 34) for the evolution of the concentration PDF becomes

∂p(c; t) ∂t = ∂c β p(c; t) ∂c + 1 + 2 ln[c M (t)/c] 2 ln[c M (t)/c] + 1 -2 ln c M (t) c m + c β-1 m c m p(c m ; t) πc 2 M (t)p(c; t) Da M 2 (t) . ( 41
)
Consider first the conservative case, corresponding to omitting the reactive terms involving c β and c β m . As explained above, we consider a Gaussian initial condition at t = Da, and take the detection limit to be sufficiently low for the mass corresponding to concentrations above c m to be approximately equal to the total mass for all times of interest. Then, M (t) = 1, and the peak concentration is given by Eq. ( 14) with a spatial variance σ 2 (t) = σ 2 (0) + (t -Da)/ Da. Using Eq. ( 16) for the Gaussian concentration PDF, it can be verified by direct calculation that Eq. ( 34) holds identically for all c and t. As a direct verification of the combined effect of reaction and diffusion on the concentration PDF, consider the particular case of a linear decay reaction, r(c) = c, for which analytical solutions are available. In this case, the concentration profile remains exactly Gaussian as before, but with a total mass that evolves in time as M (t) = exp(-t). The spatial variance remains unchanged with respect to the conservative case, and Eq. ( 16) remains valid with c M (t) given by Eq. ( 14) as before. Once again, direct computation of both sides of Eq. [START_REF] Sund | Effective models for transport in complex heterogeneous hydrologic systems[END_REF] shows that they are identical for all values of c and t.

Equation (41) quantifies the full evolution of the concentration PDF and sheds light on the fundamental governing processes. However, solving it directly is challenging. We will thus focus instead on extending the weak-coupling approach to directly compute the concentration PDF and associated quantities based on the evolution of the total mass M (t) and the peak concentration c M (t). This will allow us to obtain additional analytical insights into the dynamics. [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] focused on the evolution of the total mass and spatial variance. For the purpose of determining the concentration PDF, and in particular the mean concentration above a detection threshold and the associated variability, it is more convenient to quantify the dynamics in terms of the evolution of total mass and peak concentration, as the latter plays a central role (see Eqs. ( 12), [START_REF] Givi | Probability density function calculations in turbulent chemically reacting round jets, mixing layers and one-dimensional reactors[END_REF], and ( 21)). At time t = Da, when we begin approximating the concentration profile as Gaussian, the mass and peak are approximated as

M (Da) = c M (Da) = c B (Da) = [1 -(1 -β) Da] 1 1-β , (42)
where c B (Da) is computed according to the batch solution [START_REF] Cirpka | Probability density function of steady state concentration in two-dimensional heterogeneous porous media[END_REF]. We note that in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], instead of matching the peak value to the batch dynamics at t = Da, the spatial variance was set to match that of the initial square pulse, σ 2 (Da) = 1/12. This leads to qualitatively similar results and makes no difference at late times. The present choice enforces continuity of c M (t) at t = Da, while the previous choice enforces continuity of σ 2 (t) instead. Due to Eq. ( 14), the approximation of a sharp transition between batch dynamics and a Gaussian profile cannot enforce the continuity of both, and for the current choice we have σ 2 (Da) = 1/(2π), roughly double the initial spatial variance.

The total mass dynamics are then obtained from the local reaction kinetics. That is, integrating Eq. ( 35) in space for fixed time and using Eq. ( 12) for the Gaussian profile, neglecting mass present at concentrations below the detection limit c m , we obtain the mass evolution for β > 0 as

dM (t) dt = - M β (t) √ β 2πσ 2 (t) 1-β 2 . ( 43
)
In order to obtain an equation for the peak concentration, we employ Eq. ( 35) at x = 0. This gives

dc M (t) dt = - c M (t) 2 Da σ 2 (t) -c β M (t). ( 44
)
Using Eq. ( 14) to express the spatial variance in terms of the mass and peak concentration, we obtain

dM (t) dt = - M (t) √ β c β-1 M (t), (45a)
dc M (t) dt = - πc 3 M (t) Da M 2 (t) -c β M (t). (45b)
Note that Eq. (45a) for the overall reaction rate can also be obtained directly from the concentration PDF by using Eqs. ( 18) and [START_REF] Dentz | Concentration statistics for transport in random media[END_REF]. The result follows from setting c m = 0 in the lower integration limit, which, for β > 0, is equivalent to neglecting the mass present at concentrations lower than c m , as was done for the spatial profile in Eq. ( 43). This confirms that, as expected, knowledge of the concentration PDF directly determines overall reaction rates and mass evolution, without need for additional detailed knowledge about the spatial profile.

Equations [START_REF] Venturi | Exact pdf equations and closure approximations for advective-reactive transport[END_REF], with initial conditions at time t = Da according to Eq. ( 42), fully determine the evolution of total mass and peak concentration in this regime. Given knowledge of the mass and peak concentration, the concentration PDF for t Da, based on the Gaussian assumption, is given by Eq. ( 16). In the examples considered here, we fix the minimum concentration, corresponding to an arbitrary lower detection limit, to c m = 10 -6 . Recall that, in nondimensional units, this refers to a detection limit six orders of magnitude below the initial concentration characterizing the rectangular pulse injection. The results are qualitatively similar for different c m , so long as the peak concentration c M (t) remains large compared to c m for the times of interest. If c M (t) approaches c m , an appreciable amount of mass loss will start occurring as it drops below the detection limit. The system will then appear completely depleted, with all detectable mass and concentrations reaching zero.

In [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], good agreement was shown to hold between the total mass and spatial variance predicted by the weak-coupling approximation and numerical simulations of Eq. [START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Here, we verify the predictions of this approach, cast in terms of the total mass and peak concentration, for the full concentration PDF, as well as the time-evolution of the associated mean concentration and concentration variance. To this end, we first conduct resolved numerical simulations of the coupled diffusion and reaction problem, Eq. [START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. To compute the quantities associated with model predictions, we employ Eqs. [START_REF] Venturi | Exact pdf equations and closure approximations for advective-reactive transport[END_REF], which constitute ordinary rather than partial differential equations and can be efficiently solved numerically.

The numerical methods used for both these computations are detailed in Appendix C. Furthermore, we obtain analytical predictions for late-time behaviors under the weak-coupling approximation.

In order to test the applicability of the Gaussian assumption for predicting concentration variability at times t Da, we note that the concentration PDF associated with a Gaussian profile, Eq. ( 16) can be expressed in terms of concentrations rescaled by their (time-dependent) peak value, c = c/c M (t), as

p (c ; t) = c M (t)p[c M (t)c ; t] = f p (c ) H(1 -c )H[c -c m (t)] ln[1/c m (t)] , (46)
where c m (t) = c m /c M (t), and

f p (c ) = 1 2c ln(1/c ) . ( 47
)
The concentration PDF p (c ; t) of rescaled concentration is independent of time except through the lower cutoff c m (t) (on which it depends in terms of both its range and normalization). The functional dependency f p (c ) on rescaled concentration under the Gaussian approximation is both timeindependent and independent of the system parameters β and Da. In order to compare this to the resolved numerical results, below we compute the functional dependency according to Eq. ( 46), i.e., as

ln[1/c m (t)]c M (t)p[c M (t)c ; t].
Similarly, to assess the quality of the Gaussian approximation regarding the scalar dissipation rate, we consider the rescaled form (see Eq. ( 40))

χ (c ; t) = Da M 2 (t) 2πc 4 M (t) χ[c M (t)c , t] = f χ (c )H(1 -c )H[c -c m (t)], (48)
where the functional dependency on rescaled concentration is given by

f χ (c ) = c 2 ln(1/c ). (49)
Below, we employ Eq. ( 48) to compute the functional dependency from numerical simulations, similarly to the case above for the concentration PDF. Comparison of the concentration PDF functional forms provides a metric for the performance of the Gaussian approximation regarding the relative frequency of occurrence of different concentration values, and comparison of the scalar dissipation rate functional forms provides complementary information about the approximation's performance regarding the gradient values associated with each concentration value.

As they correspond to qualitatively different regimes, we consider the cases β < 1, 1 < β < 3, and β > 3 separately below. In what follows, we also discuss analytical approximations for late-time behaviors. Before proceeding, we note that the starting point for these results is the fact that, for t Da, the concentration profile may be approximated as Gaussian with a spatial variance that increases diffusively as σ 2 (t) ≈ t/ Da [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF]. Substituting into Eq. ( 43) and solving the resulting equation, we obtain an approximation for the mass evolution for t Da,

M (t) ≈ M (Da) 1-β - 2 √ β 1 -β 3 -β 2π Da 1-β 2 t 3-β 2 -Da 3-β 2 1 1-β , (50)
valid for β = 1, 3. Recall that the case β = 1, corresponding to a first-order decay reaction, can be fully solved analytically. The case β = 3, which we do not analyze in detail here for brevity, behaves qualitatively similarly to β > 3, although extremely slow (logarithmic) decay of mass occurs at late times [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF].

Reaction order β < 1 -Full-depletion regime

For β < 1, reaction rates increase with decreasing concentration, leading to full depletion in finite time. Under fully-mixed conditions, the depletion time t B f is given by Eq. [START_REF] Chiogna | Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume[END_REF]. Recall that for sufficiently large Da full depletion occurs before transport must be taken into account. We thus focus on Damköhler numbers such that Da < t B f .

Under the combined effect of reaction and diffusion, mass is still fully depleted in a finite time t f , different from t B f . Setting M (t f ) = 0 in Eq. ( 50) and solving for t f yields [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] t f ≈ Da

3-β 2 + √ β 2 3 -β 1 -β Da 2π 1-β 2 [1 -(1 -β) Da] 2 3-β . ( 51
)
For sufficiently small detection limit c m , we expect this prediction to provide a good estimate of the time at which the peak, mean concentration, and concentration variance reach zero. Interestingly, the depletion time decreases with decreasing Da. This decrease is due to the fact that, because in this regime reaction is faster at lower concentrations, stronger diffusion (smaller Da) leads to faster depletion when compared to a system with the same reaction rate. For this reason, both mass and peak concentration decay faster than the batch solution. In particular, for small Da 1, we find a simple power-law scaling with Da,

t f ≈ √ β 2 3 -β 1 -β 2 3-β Da 2π 1-β 3-β . ( 52
)
Note however that t f / Da increases with decreasing Da, meaning that with increasing diffusion the depletion time occurs increasingly later than the time t = Da at which the Gaussian approximation kicks in. 9) and ( 47)).

To illustrate these results, we take β = 1/2 as a representative example. First, we compare the functional dependency of the concentration PDF (see Eqs. ( 46) and ( 47)) under the weak-coupling approximation to that obtained from resolved numerical simulations of Eq. ( 35) (see Appendix A and Appendix C). We focus on times t such that Da t < t f (see Eq. ( 51 40) and ( 49)).

of concentration. This discrepancy at low concentrations also leads to a corresponding overprediction of higher concentrations, due to the normalization condition for the PDF. Although the bulk of the mass distribution is generally well captured, as discussed below, high Da leads to significant depletion at early times and results in a worse approximation.

Similar conclusions are reached from consideration of the functional dependency of the scalar dissipation rate (see Eqs. ( 48) and ( 49)), shown in Fig. 4. The Gaussian approximation again provides good estimates, but higher scalar dissipation rates, corresponding to higher gradients, are observed at low concentration values. The discrepancy is more pronounced in the same cases and for the same reasons as discussed regarding the concentration PDF, although in this case normalization does not play a role and thus no associated discrepancy at higher concentration values is observed.

Next, we turn to the evolution of the mass and peak concentration, shown in Fig. 5. As predicted by the weak-coupling approximation, both quantities reach zero in finite time as complete depletion occurs. However, due to the nonlinear nature of the reaction and the presence of diffusion, the peak dynamics differ from the total mass dynamics, as captured by Eqs. [START_REF] Venturi | Exact pdf equations and closure approximations for advective-reactive transport[END_REF]. Note the onset of a diffusive scaling of peak decay for low Damköhler, showing that the dynamics are controlled by diffusion at times t ∼ Da (corresponding to dimensional times on the order of the diffusion time). For late times, as concentration becomes smaller, reaction becomes more efficient. The dynamics depart from the diffusive scaling and full depletion eventually occurs. The slight underprediction of the depletion time t f by the analytical approximation compared to the numerically-computed weak-coupling approximation results is due mostly to the breakdown of the linear variance growth approximation just before complete depletion, as observed in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF]. When comparing the predictions of the weak-coupling approximation to the resolved numerical simulations, we find overall good agreement. Regarding peak concentrations, the approximation of a sharp transition between batch and Gaussian-profile dynamics at t = Da leads to some discrepancy at times t ∼ Da, which becomes more pronounced with increasing Da. As the Damköhler number increases and the depletion time approaches Da from above (t f / Da ∼ 1, see Eq. ( 51)), stronger discrepancies are observed for both the mean concentration and the total mass, due to significant depletion before an approximately Gaussian profile is achieved. Nonetheless, the predictions for t f still provide a reasonable approximation. Recall that, for higher Da > t B f = 2 (Eq. ( 39)), diffusion does not play a role before depletion the dynamics become well described by simple batch kinetics. 45)), and dash-dotted lines show the purely-diffusive behavior. The vertical dashed lines show approximate analytical results for the depletion time (Eq. ( 51)).

Note that defining a mean concentration in terms of a fixed, reference spatial region, as done in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], leads to a mean concentration scaling identical to that of the total mass, since in that case

C(t) = M (t)/|Ω| with constant domain size |Ω|.
Interestingly, when the mean concentration is instead defined, as here, in terms of the temporally-variable region where concentration is above a detection limit, the concentration mean and variance are instead fully controlled by the peak dynamics, see Eqs. ( 19) and ( 21). This leads also to a similar quality of the weak-coupling approximation when compared to the resolved simulations, as shown in Fig. 6. In particular, for low Damköhler numbers, our approach accurately predicts a transition from diffusion-to reaction-controlled dynamics, governed by the peak behavior described above, for both the mean concentration and the concentration variance.

Reaction order 1 < β < 3 -Power-law-decay regime

In this regime, the interplay between reaction and diffusion leads to power-law mass decay at late times. Indeed, for t Da, Eq. ( 50) reduces to [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF]

M (t) ≈ 2π Da √ β 2 3 -β β -1 1 β-1 t -3-β 2(β-1) . (53)
Substituting in Eq. (45b), and noting that, for 1 < β < 3, as the peak concentration decreases the reactive term eventually dominates over the diffusive term, we have for late times

dc M (t) dt ≈ -c β M (t). ( 54
)
Solving this equation, we obtain the late time peak decay

c M (t) ≈ [(β -1)t] -1 β-1 . ( 55
)
This behavior is characteristic of reaction-dominated conditions, and identical to the late-time batch dynamics, see Eq. [START_REF] Cirpka | Probability density function of steady state concentration in two-dimensional heterogeneous porous media[END_REF]. However, note that the total mass, Eq. ( 53), exhibits a different temporal scaling, as it results from the reactive decay of Gaussian-distributed concentration values undergoing diffusive spreading. This results, in contrast to the β < 1 regime, in slower asymptotic mass decay than predicted by batch dynamics, because in this case reaction is slower at lower concentrations. 45)), and dash-dotted lines show the purely-diffusive behavior. The vertical dashed lines show approximate analytical results for the depletion time (Eq. ( 51)).

We now take β = 2 as a representative example. Again, we start by assessing the performance of the Gaussian assumption for the functional dependency of the concentration PDF (Fig. 7). The approximation provides very good estimates of the frequency of occurrence of different concentrations.

A slight underestimation of low concentrations, and corresponding overestimation of intermediate concentrations, occurs at late times for large Damköhler number. When considering the scalar dissipation rate (Fig. 8), we again observe very good agreement, with a corresponding slight but more discernible overestimation of concentration gradients at low concentrations. In this case, reaction is less efficient at low concentrations, so these results are somewhat counter-intuitive. Indeed, reaction is stronger and therefore has a more pronounced effect compared to diffusion at high concentrations. However, when the corresponding spatial profile is compared to a Gaussian of the same variance and maximum concentration, the latter underestimates the intermediate values but overestimates the tails. This is illustrated in Fig. 9 for the case of Da = 10 2 at time t = 10 Da. This effect is slight, but it is more pronounced at the level of concentration gradients.

Next, we consider the evolution of total mass and peak concentration, see Fig. 10. The analytical approximation captures the asymptotic scalings correctly, although it slightly underpredicts the mass and overpredicts the peak concentration when compared to the numerically-computed weak-coupling approximation results. These discrepancies result in this case from small corrections to the σ 2 (t) ≈ t/ Da diffusive spatial variance growth approximation of [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF], implicit in Eq. ( 50) for the approximate late-time mass evolution. Again, the onset of a purely-diffusive scaling of peak decay is visible for low Damköhler before the asymptotic regime is reached. The weak-coupling approximation provides very good estimates of the resolved dynamics, even for Da ∼ 1. As expected, the largest discrepancies occur for times t ∼ Da, around which the approximation of a sharp transition between batch and Gaussianprofile dynamics has the most impact. In addition, a sharp depletion of mass can be observed at late times for the high-Da cases due to the detection limit. As expected, the transition from the dynamics predicted by the weak-coupling approximation to this sharp depletion occurs when the peak concentration becomes comparable to c m .

As for the β < 1 case, the concentration mean and variance (Fig. 11) are controlled by the peak dynamics, see Eqs. ( 19) and [START_REF] Rodney | Computation of turbulent reactive flows: first-principles macro/micromixing models using probability density function methods[END_REF]. However, the changes in the size of the domain where concentration is larger than the detection limit lead to logarithmic corrections to the purely-reactive power-law scaling observed for the peak concentration. As for the mass and peak, the weak-coupling approximation 9) and ( 47)).

provides very good predictions of the full dynamics, with the quality of the approximation improving for late times t 10 Da. The discrepancies in the analytical vs numerical results for the weak-coupling approximation result from the discrepancies in the computation of the peak concentration discussed above.

Reaction order β > 3 -Incomplete-depletion regime

In this case, the substantial reaction slowdown associated with decreasing concentrations, enhanced by diffusion-induced dilution, prevents the complete depletion of mass by reaction, even for arbitrarily large times. Taking the limit t → ∞ in Eq. ( 50) yields the leftover mass [7]

M ∞ ≈ 1 + 1 + 2 (β -3) (2π) β-1 β (β -1) Da -1 β-1 . (56)
Substituting in Eq. (45b) as before, we now conclude that the late-time dominant contribution is diffusive, so that

dc M (t) dt ≈ - πc 3 M (t) Da M 2 ∞ . (57)
Solving this equation leads, for large times, to

c M (t) ≈ M ∞ Da 2πt . (58)
As the mass becomes constant at late times, the decay of peak concentration exhibits the signature of diffusive plume spreading. The mass in a well-mixed batch decreases asymptotically to zero, meaning the diffusing system leads to less efficient reaction, as expected since the reaction slows down strongly with decreasing concentration in this case. In contrast, however, the diffusion-governed peak decay is asymptotically faster than the corresponding batch prediction, Eq. [START_REF] Cirpka | Probability density function of steady state concentration in two-dimensional heterogeneous porous media[END_REF].

To illustrate this regime, we take β = 4 as a representative example. As predicted by the weakcoupling approximation, reaction is sufficiently inefficient at low concentrations for diffusive spreading 40) and ( 49)).

to completely arrest mass decay. This means that the system dynamics are essentially dominated by diffusion. Correspondingly, both the concentration PDF and the scalar dissipation rate show excellent agreement with the Gaussian prediction, as shown in Figs. 12 and 13, respectively.

Next, we consider the evolution of total mass and peak concentration, see Fig. 14. In this case, the analytical asymptotic predictions provide very accurate estimates of the weak-coupling dynamics. As expected from the analytical discussion, the results become indistinguishable from a purely diffusive system for sufficiently low Da. However, when comparing to the resolved numerical simulations, we find that the weak-coupling approximation underpredicts the asymptotic mass at higher Da. The main reason for the discrepancy is that, because the reaction is particularly inefficient at low concentrations, the system transitions to the Gaussian regime earlier in this case. Indeed, as shown by the dash-dotted lines in Fig. 14, assuming the transition occurs at t = Da /10 rather than at t = Da leads to improved estimates of the asymptotic mass. To obtain this prediction, we proceed as for Eq. ( 56) but set the initial condition of the weak-coupling Gaussian regime at t = Da /10 (see Eq. ( 42)), leading to

M ∞ ≈ 1 + 1 + 2 • 10 β-1 2 (β -3) (2π) β-1 β (β -1) Da 10 -1 β-1 . ( 59
)
Note that the scaling of asymptotic mass with Da, and the relationship between asymptotic peak dynamics and asymptotic mass, Eq. ( 58), remain unchanged.

The concentration mean and variance for this case are shown in Fig. 15. The considerations pertaining to the 1 < β < 3 case regarding the scaling of concentration mean and variance in terms of the peak concentration remain valid here, but in this case we observe logarithmically-corrected diffusive rather than reactive scalings. Once again, the low-Da case is indistinguishable from purely-diffusive dynamics. In addition, due to the low efficiency of reaction at low concentrations, the scaling behavior of both mean concentration and concentration variance remains purely diffusive, even for large Da.

Overall, as for the total mass and peak concentration, the analytical approximation provides very good asymptotic estimates of the weak-coupling approximation, but better predictions of the fully-resolved dynamics are obtained by assuming the transition to the Gaussian regime occurs at time t = Da /10 rather than at t = Da. 35), and the dashed colored line shows a Gaussian profile with the same variance σ 2 (t) and maximum concentration c M (t). The horizontal dashed line shows the lower detection limit.

Discussion and conclusions

We have presented a detailed discussion of the theoretical and numerical properties of the concentration PDF of a reactive plume evolving under the coupled action of dilution and nonlinear reaction. The concentration PDF encodes the full point statistical variability of concentration values found throughout a spatial domain, and is therefore of central interest to both fundamental and applied reactive transport and risk assessment problems.

As a first step towards solving the full advection-diffusion-reaction system, we have conducted a detailed analysis of nonlinear decay of a single solute undergoing one-dimensional diffusion. We have shown that the weak-coupling approximation developed in [START_REF] Le Traon | Effective kinetics driven by dynamic concentration gradients under coupled transport and reaction[END_REF] to quantify the evolution of total mass and mean concentration in a fixed reference domain can be extended to quantify the evolution of mean concentrations above a lower detection limit, along with the associated variability. In particular, we have found that the late-time Gaussian profile approximation leads to good predictions of the concentration PDF. Even for this simple one-dimensional diffusive system, the effective reaction kinetics governing these quantities exhibit rich dynamics. In contrast to well-mixed batch kinetics, the decay of concentrations at different Damköhler numbers can be dominated by diffusive or reactive effects at different times. Late-time concentrations can exhibit complete extinction, power-law decay, or complete arrest of reaction depending on the reaction order. The Gaussian approximation underlying the weak-coupling approach allows us to predict the associated concentration variance. Remarkably, the dynamics of the maximum concentration value, which provide the key control on the mean and variance of concentrations above a given detection threshold, differ qualitatively from the dynamics of total mass. In particular, for low Damköhler numbers, both the mean concentration and associated variability exhibit diffusion-controlled decay for times on the order of the diffusion time, for any reaction order. This diffusive control persists across all times and Damköhler numbers for strongly superlinear decay reactions (of order β > 3), whereas reaction is the fundamental mechanism controlling late-time and/or large-Da dynamics for lower reaction orders.

We have also derived a dynamical equation for the evolution of the concentration PDF in the nonlinear reaction-diffusion problem, showing that the impact of transport is encoded in the scalar dissipation rate seen as a function of concentration. It is important to note that the generalization of this equation to multiple dimensions, variable advection and/or diffusion tensor, and arbitrary reaction rates, for which concentration is no longer necessarily monotonically decreasing, is not straightforward. 53) and ( 55)), and dash-dotted lines show the purely-diffusive behavior.

Nonetheless, we believe the approach developed here, highlighting the role of the scalar dissipation rate, has the potential to form the basis for upscaling procedures in more complex systems. In particular, as discussed in Section 3, the dynamics of reaction under advection-diffusion in heterogeneous flows, in both two-and three-dimensional porous media, can be reduced to one-dimensional reaction-diffusion equations along the principal direction of compression on a collection of lamellar structures. With this in mind, further work will focus on the generalization of the approach developed here to the dynamics of Gaussian profiles subject to stretching-enhanced diffusion over an ensemble of lamellae. 19) and ( 21), using Eq. ( 55) for the peak concentration), and dash-dotted lines show the purelydiffusive behavior. estimated as

p(c; t; N g ) = 1 N g Ng i=1 δ[c -c i (t)], (60)
where N g = |Ω|/V g is the number of grid points and c i is the average concentration in cell i. Note that the c i associated with different cells are not necessarily all different.

Further discretizing concentration into bins leads to a probability mass function with values for the probability of concentration in each bin k+1) [, k 0:

B k = [c (k) , c (
p k (t; N g ) = 1 N g Ng i=1 H[c i (t) -c (k) ]H[c (k+1) -c i (t)], (61)
where H(•) is the Heaviside step function. This means that the probability of finding a concentration value in bin k is the fraction of cells where the concentration falls within bin k. If we take concentration bin widths ∆c to be constant, c (k) = k∆c for k 0, and approximate the PDF of concentration by dividing probabilities by ∆c (as in Fig. 2), we obtain

p(c; t; N g , ∆c) = 1 N g ∆c Ng i=1 H[c i (t) -k∆c]H[(k + 1)∆c -c i (t)]. (62)
These discretization procedures generalize directly to a multispecies system. In that case, concentration bins refer to the simultaneous attainment of concentration values of each species. The associated probabilities are computed as above by counting spatial cells where these values occur simultaneously.

The procedure also generalizes directly to non-uniform and/or time-dependent spatial cell sizes and/or concentration bins. More involved techniques employing kernel reconstructions of the concentration field may be formalized in a similar manner [START_REF] Vlad I Morariu | Automatic online tuning for fast gaussian summation[END_REF][START_REF] Fernàndez | Optimal reconstruction of concentrations, gradients and reaction rates from particle distributions[END_REF][START_REF] Sole-Mari | Particle density estimation with grid-projected and boundary-corrected adaptive kernels[END_REF].

The maximum resolution of the discretized concentration PDF described above is given by 1/(N g ∆c),

and the maximum PDF value is 1/∆c. We now discuss the impact of local spatial extrema, which can be associated with divergences of the continuous PDF as discussed in Section 2.3, on the discretized computation. It is also important to note that sources of error typically come into play in the On the other hand, smooth extrema correspond to divergences only for spatial dimension d = 1, for which they lead to an inverse-square-root divergence, as shown in Section 2.3. In this case, averaging the concentration over a range ∆c near the spatial extreme value c 0 yields p(c 0 ; t; ∆c) ∝ 1/ √ ∆c, so that p(c 0 ; t; a∆c) ≈ √ ap(c 0 ; t; ∆c). Thus, observing these behaviors in a numerical computation is a signature of the presence of a spatial extremum, and the scaling behavior with concentration discretization refinement indicates its type.

B Dynamical equation for the concentration PDF

In this appendix, we provide a detailed derivation of the dynamical equation for the concentration PDF discussed in the main text, Eq. ( 34), by explicitly determining the diffusive transport and reaction contributions in Eq. ( 29). We first consider the effect of diffusion, ∆p D (c; t). For one-dimensional diffusion and nonlinear concentration decay, assuming a symmetric initial condition about the origin and monotonically decreasing with distance from the latter, the spatial concentration profile retains these properties for all times. Therefore, we have |Λ(c; t)| = 2 for all concentrations within the range observed at time t, corresponding to the two points x = ±x c (c; t) where C(x; t) = c. Furthermore, the concentration gradient magnitude is the same at ±x c , so that its harmonic average is simply g h (c; t) = ∇C[-x c (c; t), t] = -∇C[x c (c; t), t]. Using Eq. ( 8) and the fact that |Λ(c; t)| is time-independent within 40) and ( 49)).

the concentration range observed, we can write the change in the concentration PDF as

∂p(c; t) ∂t = -p(c; t) ∂ ln |Ω(t)| ∂t + ∂ ln g h (c; t) ∂t . (63)
If a fixed reference volume is considered, the term corresponding to the change of |Ω(t)| in time is zero. As discussed in the main text, we focus here on the case of a minimum detection limit c m and a time-varying domain Ω(t) where c > c m .

In order to compute the terms in square brackets, we consider the time evolution of quantities on a given concentration surface (in one dimension, at the points ±x c (c; t)). By definition, the change in time of concentration over such a surface is zero, so that

and we find

∂ ln g h (c; t) ∂t = 1 2 ∂ 2 χ(c; t) ∂c 2 - 1 4 ∂ ln χ(c; t) ∂c ∂χ(c; t) ∂c . (67)
The change in domain volume Ω(t) = 2x c (c m ; t) due to dilution of concentration below the detection limit c m obeys 56) and ( 58)). Dash-dotted lines show analytical results assuming the transition to the Gaussian regime occurs at time t = Da /10 rather than at t = Da (Eq. ( 59)).

∂|Ω(t)| ∂t = 2 ∂x c (c m ; t) ∂t , (68)
Using Eq. ( 65), this leads to 19) and ( 21), using Eqs. ( 56) and ( 58) for the peak concentration). Dash-dotted lines show analytical results assuming the transition to the Gaussian regime occurs at time t = Da /10 rather than at t = Da (Eqs. ( 19) and ( 21), using Eqs. ( 58) and ( 59) for the peak concentration).

∂|Ω(t)| ∂t = 1 g h (c m ; t) ∂χ(c; t) ∂c c=cm , (69)
Note that integration of the right-hand side from c = c m to c = ∞ yields zero, which ensures the reactive contribution conserves probability for arbitrary r(c). Substituting the effects of transport and reaction, Eqs. ( 33) and [START_REF] Hamada | Diffusion limited mixing in confined media[END_REF], in Eq. ( 29) leads to the dynamical Eq. ( 34) for the evolution of the concentration PDF under one-dimensional diffusion and nonlinear decay.

C Numerical methods

In this appendix, we provide details on the numerical methods used to integrate the weak-coupling equations ( 45) and the reaction-diffusion equation [START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Regarding Eqs. [START_REF] Venturi | Exact pdf equations and closure approximations for advective-reactive transport[END_REF], which are ordinary differential equations, we implemented a standard fourth-order Runge-Kutta method in the C++ language.

This method was chosen for its simplicity of implementation and high accuracy, and also because, as an explicit method, it provides a convenient approach to integrate these nonlinear equations without requiring numerical root-finding methods. We employed a time step ∆t = 10 -2 min{Da, 1/ Da} for the temporal discretization, which we verified led to consistently converged results.

For the fully-coupled reaction-diffusion problem, Eq. ( 35), we employed the py-pde open-source Python package for solving partial differential equations [START_REF] Zwicker | py-pde: A python package for solving partial differential equations[END_REF]. We used a regular finite difference discretization of a one-dimensional domain of half-width L and second-order centered differences for the spatial derivative approximations. For the time integration, we employed an explicit Forward Euler scheme. We set reflecting boundary conditions at the edges of the computational domain, but we verified that the latter was sufficiently large that no appreciable mass reached the edges, rendering the choice of boundary conditions irrelevant. Since the late-time variance growth is approximately diffusive, L = 10 t m / Da, where t m is the maximum simulation time, may be used as a simple estimate of necessary domain size. However, because of the lower detection limit c m = 10 -6 used in the computation of the quantities of interest, we found that in practice it was never necessary to use L > 1500 for the simulations conducted here. We note that, for β < 1, where complete depletion of concentrations can happen in finite time, the increase of reaction rates with decreasing concentration values can lead to numerical issues, because very low concentrations can drop below zero within a time step. We avoid this issue by setting negative concentrations to zero before computing reaction rates.

We chose the spatial and temporal discretizations so as to ensure good accuracy while maintaining reasonable simulation times. The discretization parameters for different system parameters β and Da are summarized in Table 1 for the mass and concentration peak, mean, and variance calculations, and in Table 2 for the concentration PDF and scalar dissipation rate.

The concentration PDF was obtained by counting discretized spatial locations where the concentration value fell within prescribed bins (see also Section A). In order to accurately resolve both low and high concentrations, we employed n logarithmically-spaced concentration bins for concentrations between the lower detection limit c m and 2c M (t)/3, where the time-dependent peak value c M (t) was determined from the numerical concentration profiles, and n h linearly-space bins for the remaining concentrations between 2c M (t)/3 and c M (t). For the lower, intermediate, and higher time examined in each case, we employed (n , n h ) = (20, 10), [START_REF] Alisse | Experimental probability density functions of small-scale fluctuations in the stably stratified atmosphere[END_REF][START_REF] Ee O'brien | The probability density function (pdf) approach to reacting turbulent flows[END_REF], and [START_REF] Stephen B Pope | Lagrangian pdf methods for turbulent flows[END_REF][START_REF] Hubert | Enhanced and non-monotonic effective kinetics of solute pulses under michaelis-menten reactions[END_REF], respectively. The scalar dissipation rate was calculated according to Eq. (77b) by numerically computing the spatial derivative at each discretized spatial location (in the rising limb of the symmetric concentration profile), using second-order central differences. The corresponding concentration values at each spatial location were recorded and used to obtain the scalar dissipation rate as a function of concentration.

  Di↵usive flux < l a t e x i t s h a 1 _ b a s e 6 4 = " f O s W u y 5 o z A h / C j 8 h u C R W 4 1 c s 5 P 0 = " > A A A B 9 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e w G R Y 9 B P X i M Y B 6 Q r G F 2 M p s M m Z 1 d 5 h E N S / 7 D i w d F v P o v 3 v w b J 8 k e N L G g o a j q p r s r S D h T 2 n W / n d z K 6 t r 6 R n 6 z s L W 9 s 7 t X 3 D 9 o q N h I Q u s k 5 r F s B V h R z g S t a 6 Y 5 b S W S 4 i j g t B k M r 6 d + c 0 S l Y r G 4 1 + O E + h H u C x Y y g r W V H m 5 Y G B r F R h S F 3 D x 1 i y W 3 7 M 6 A l o

  t e x i t s h a 1 _ b a s e 6 4 = " h W L m a 7 P q g S E z N P B x e G b w / y z f 1 / Q = " > A A A C B n i c b V D J S g N B E O 1 x j X G L e h S h M Q h e D D M S 0 Y s Q 9 O I x g l k g E 0 N N T 0 / S p G e h u 0 Y I Q 0 5 e / B U v H h T x 6 j d 4 8 2 / s L A d N f F D w e K + K q n p e I o V G 2 / 6 2 F h a X l l d W c 2 v 5 9 Y 3 N r e 3 C z m 5 d x 6 l i v M Z i G a u m B 5 p L E f E a C p S 8 m S g O o S d 5 w + t f j / z G A 1 d a x N E d D h L e D q E b i U A w Q C N 1 C g d u o I B l P h t m P g 4 v T 9 w + J A l Q d u 9 6 H K F T K N o l e w w 6 T 5 w p K Z I p q p 3C l + v H L A 1 5 h E y C 1 i 3 H T r C d g U L B J B / m 3 V T z B F g f u r x l aA Q h 1 + 1 s / M a Q H h n F p 0 G s T E V I x + r v i Q x C r Q e h Z zp D w J 6 e 9 U b i f 1 4 r x e C i n Y k o S Z F H b L I o S C X F m I 4 y o b 5 Q n K E c G A J M C X M r Z T 0 w u a B J L m 9 C c G Z f n i f 1 0 5 J T L p 3 d l o u V q 2 k c O b J P D s k x c c g 5 q Z A b U i U 1 w s g j e S a v 5 M 1 6 s l 6 s d + t j 0 r p g T W f 2 y B 9 Y n z 9 T R p k I < / l a t e x i t > Dilution by mixing with nonreactive resident fluid < l a t e x i t s h a 1 _ b a s e 6 4 = " s T u n

  t e x i t s h a 1 _ b a s e 6 4 = " / o Y E + g F w y P 7 V N P W R 5 N j m 0 L l u B e s = " > A A A B + 3 i c b V D L S g M x F L 1 T X 7 W + a l 2 6 C R b B V Z k p i i 6 L b l x J B f u A d i i Z N N O G Z p I h y Y h l 6 K + 4 c a G I W 3 / E n X 9 j p p 2 F t h 4 I H M 6 5 h 3 t z g p g z b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c N

  4 l 7 X r N O 6 9 d 3 N e r j e u 8 j i I c w w m c g Q e X 0 I B b a E I L C D z B M 7 z C m z N z X p x 3 5 2 M x W n D y z B H 8 g f P 5 A x / 3 l I A = < / l a t e x i t >

Figure 2 :

 2 Figure 2: Illustration of the key concepts underlying the concentration PDF. (a): Example spatial concentration profile, a Gaussian with unit variance and unit height, restricted to the spatial domain Ω = [-3, 3] of size |Ω| = 6. (b): Corresponding concentration PDF. Given a spatial profile, the concentration PDF quantifies the spatial frequency of occurrence of different concentration values c = C(x) occurring at some position(s) x. When concentration values are discretized into bins of width ∆c, the PDF value p(c) is approximated by the fraction of the domain of the domain where the concentration c falls in this range, divided by ∆c to obtain a probability density (bars in panel b). The continuous PDF (solid line in panel b) is obtained in the limit of infinitesimal bin widths dc. Each infinitesimal region dx around a point where C(x) = c contributes a probability density (dx/dc)/|Ω| = 1/(|Ω||dC(x)/dx|). This causes a divergence of the PDF at the maximum concentration c M = 1, where the gradient is zero, and contributes to a large value at the minimum concentration c m = C(3) = C(-3), where the gradient is small. These concepts are quantified and explored in detail throughout the text.

  the qualitative behavior near an extremum is determined by the Hessian determinant |H| and the distance |∆x| to the maximum; thus, the gradient near an extremum obeys |∇C| ∝ |H| 1/d |∆x|, and the concentration changes as |∆c| ∝ |∇c||∆x| ∝ |H| 1/d ∆x 2 . The second equation gives |∆x| ∝ |H| -1/(2d) |∆c| 1/2 , and substituting in the first leads to p(c 0 + ∆c) ≈ |Λ(c 0 + ∆c; t)|/[|Ω(t)||∇C|] ∝ |Λ(c 0 + ∆c; t)|/[|Ω(t)||H| 1/(2d) |∆c| 1/2 ]. Substituting |Λ(c 0 + ∆c; t)| ∝ |∆x| d-1 for the (d -1)-area of the surface of concentration c 0 + ∆c yields

  ∆p R (c; t) = ∂r(c)p(c; t) ∂c + p(c; t)p(c m ; t)r(c m ). (30) Note that it is completely determined given knowledge of the concentration PDF and local reaction rate as a function of concentration. The first term represents the net change in the probability density of concentration c due to reaction, and the second term guarantees the PDF remains normalized by accounting for reactive decay of concentrations below the threshold value c m .

  determined given knowledge of the local reaction rates r(c) and the concentration-dependent scalar dissipation rate χ(c; t). The first (reactive) term on the right-hand side has the form of an advective contribution: the reaction rate r(c) functions as an advective "velocity" (units of inverse concentration per time) generating a net flux of probability towards lower concentrations. The remaining contributions have the form of concentration-dependent sink/source terms. The first two terms in square brackets, which are governed by the scalar dissipation rate, encode the impact on the PDF of changes in concentration values and the associated gradients due to diffusive transport. The remaining terms concern the reactive and diffusive fluxes of probability towards concentration values below the detection limit. These affect the concentration PDF through the normalization across the remaining concentration values, since c M (t) cm dc p(c; t) = 1 by definition.

Figure 3 :

 3 Figure 3: Functional dependency of the concentration PDF in terms of concentration rescaled by peak value, for β = 1/2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 10 -1 , and (c): fast reaction, Da = 1. Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Guilbert | Chemical reaction for mixing studies[END_REF] under the weak-coupling model. The dashed lines show a pure 1/c scaling for reference at low concentrations and the scaling of the inverse-square-root divergence near the maximum concentration value (see Eqs. (9) and (47)).

Figure 4 :

 4 Figure 4: Functional dependency of the scalar dissipation rate in terms of concentration rescaled by peak value, for β = 1/2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 10 -1 , and (c): fast reaction, Da = 1. Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Le Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF] under the weak-coupling model. The dashed lines show a pure c 2 scaling for reference at low concentrations and the linear approach to zero near the peak concentration (see Eqs. (40) and (49)).

Figure 5 :

 5 Figure 5: Temporal evolution of (a): total mass and (b): peak concentration for β = 1/2. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weakcoupling model (Eqs. (45)), and dash-dotted lines show the purely-diffusive behavior. The vertical dashed lines show approximate analytical results for the depletion time (Eq. (51)).

Figure 6 :

 6 Figure 6: Temporal evolution of (a): mean concentration and (b): concentration variance for β = 1/2. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weak-coupling model (Eqs. (45)), and dash-dotted lines show the purely-diffusive behavior. The vertical dashed lines show approximate analytical results for the depletion time (Eq. (51)).

Figure 7 :

 7 Figure 7: Functional dependency of the concentration PDF in terms of concentration rescaled by peak value, for β = 2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 1, and (c): fast reaction, Da = 10 2 . Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Guilbert | Chemical reaction for mixing studies[END_REF] under the weak-coupling model. The dashed lines show a pure 1/c scaling for reference at low concentrations and the scaling of the inverse-square-root divergence near the maximum concentration value (see Eqs. (9) and (47)).

Figure 8 :

 8 Figure 8: Functional dependency of the scalar dissipation rate in terms of concentration rescaled by peak value, for β = 2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 1, and (c): fast reaction, Da = 10 2 . Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Le Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF] under the weak-coupling model. The dashed lines show a pure c 2 scaling for reference at low concentrations and the linear approach to zero near the peak concentration (see Eqs. (40) and (49)).

Figure 9 :

 9 Figure 9: Concentration profile for β = 2 and Da = 10 2 , at time t = 10 Da. The solid line shows the profile computed from resolved numerical simulations of Eq. (35), and the dashed colored line shows a Gaussian profile with the same variance σ 2 (t) and maximum concentration c M (t). The horizontal dashed line shows the lower detection limit.

Figure 10 :

 10 Figure 10: Temporal evolution of (a): total mass and (b): peak concentration for β = 2. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weakcoupling model (Eqs. (45)), dashed lines show approximate asymptotic analytical results (Eqs. (53) and (55)), and dash-dotted lines show the purely-diffusive behavior.

Figure 11 :

 11 Figure 11: Temporal evolution of (a): mean concentration and (b): concentration variance for β = 2. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weak-coupling model (Eqs. (45)), dashed lines show approximate asymptotic analytical results (Eqs. (19) and (21), using Eq. (55) for the peak concentration), and dash-dotted lines show the purelydiffusive behavior.

Figure 12 :

 12 Figure 12: Functional dependency of the concentration PDF in terms of concentration rescaled by peak value, for β = 2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 1, and (c): fast reaction, Da = 10 2 . Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Guilbert | Chemical reaction for mixing studies[END_REF] under the weak-coupling model. The dashed lines show a pure 1/c scaling for reference at low concentrations and the scaling of the inverse-square-root divergence near the maximum concentration value (see Eqs. (9) and (47)).

Figure 13 :

 13 Figure 13: Functional dependency of the scalar dissipation rate in terms of concentration rescaled by peak value, for β = 2 and (a): slow reaction, Da = 10 -2 , (b): intermediate reaction, Da = 1, and (c): fast reaction, Da = 10 2 . Colors stand for different times, with markers showing the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF]. Solid lines represent the analytical solution[START_REF] Le Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF] under the weak-coupling model. The dashed lines show a pure c 2 scaling for reference at low concentrations and the linear approach to zero near the peak concentration (see Eqs. (40) and (49)).

Figure 14 :

 14 Figure 14: Temporal evolution of (a): total mass and (b): peak concentration for β = 4. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weakcoupling model (Eqs. (45)), and dashed lines show approximate asymptotic analytical results (Eqs. (56) and (58)). Dash-dotted lines show analytical results assuming the transition to the Gaussian regime occurs at time t = Da /10 rather than at t = Da (Eq. (59)).

2 ∂χ 4 ∂ 1 2 ∂. ( 71 )( 73 )Figure 15 :

 2412717315 Figure 15: Temporal evolution of (a): mean concentration and (b): concentration variance for β = 4. Colors stand for different Damköhler numbers, markers show the results of resolved numerical simulations of the fully-coupled reaction-diffusion problem (35), solid lines are numerical results based on the weak-coupling model (Eqs. (45)), and dashed lines show approximate asymptotic analytical results (Eqs. (19) and (21), using Eqs. (56) and (58) for the peak concentration). Dash-dotted lines show analytical results assuming the transition to the Gaussian regime occurs at time t = Da /10 rather than at t = Da (Eqs. (19) and (21), using Eqs. (58) and (59) for the peak concentration).

  concentration into account. Using the same techniques as before, we obtain for the change in the domain size due to reaction:∂ ln |Ω(t)| ∂t = -p(c m ; t)r(c m ).(74)Thus, the complete effect of reaction is ∆p R (c; t) = ∂r(c)p(c; t) ∂c + p(c; t)p(c m ; t)r(c m ).

Table 1 :

 1 Discretization parameters used in computing solutions of Eq.[START_REF] Sánchez-Vila | Conditional probability density functions of concentrations for mixing-controlled reactive transport in heterogeneous aquifers[END_REF] to determine the mass and concentration peak, mean, and variance. For these computations, we employed a temporal discretization ∆t = a min{Da, 1/ Da} and a spatial discretization ∆x = b ∆t/ Da. The corresponding values of (a, b) for different reaction orders β and Damköhler numbers Da are given in the table.

	β \ Da	10 -2	10 -1	10 0	10 1	10 2
	1/2 2	(10 -4 , 1) (10 -2 , 1)	(10 -4 , 1) (10 -2 , 1)	(5 • 10 -7 , 1) (10 -2 , 1)	-(10 -1 , 2)	-(10, 4)
	4	(10 -3 , 1)	(10 -2 , 1)	(10 -2 , 1)	(10 -1 , 1)	(10 -1 , 1)

Table 2 :

 2 Discretization parameters used in computing solutions of Eq. (35) to determine the concentration PDF and scalar dissipation rate. For these computations, we employed a temporal discretization ∆t = a min{Da, 1/ Da} and a spatial discretization ∆x = ∆t/ Da. The corresponding values of a for different reaction orders β and Damköhler numbers Da are given in the table.

	β \ Da	10 -2	10 -1	10 0	10 1	10 2
	1/2 2 4	10 -4 10 -3 10 -4	10 -4 --	5 • 10 -6 10 -4 10 -4	---	-10 -2 10 -1

Availability of data and material: This work did not make use of or generate any datasets.

Code availability Numerical simulations used the py-pde open-source Python package for solving partial differential equations [START_REF] Zwicker | py-pde: A python package for solving partial differential equations[END_REF]. Please contact the authors for any additional information or requests.

Declarations

Funding: TA was supported by a Marie Sk lodowska Curie Individual Fellowship, funded by the European Union's Horizon 2020 research and innovation programme under the project ChemicalWalks 838426. CB acknowledges funding from Region Bretagne and Rennes Metropole. TLB gratefully acknowledges funding by the ERC under the project ReactiveFronts 648377.

Conflicts of interest/Competing interests:

The authors declare that they have no conflict of interest.

A Discretization

When computing the concentration PDF numerically, discretizations are typically employed both spatially and for the concentration values. As a simple example that highlights the central concepts, consider a regular spatial discretization into a grid with constant cell volume V g , with each grid cell associated with the average concentration within it. According to the argument for constant concentration regions developed in the previous section, for each fixed time t, the concentration PDF is then

D Problem setup and nondimensionalization

This appendix provides additional details on the nondimensionalization used in Section 4. Denoting nondimensionalized quantities by an asterisk, we have

where the characteristic reaction time τ R = κ -1 c 1-β 0

. The minimum and maximum concentrations c M (t) and c m (t) are normalized in the same manner. Similarly, the nondimensional spatial variance

The concentration PDF and scalar dissipation are then nondimensionalized accordingly as

with the Damköhler number Da = τ D /τ R .

Note that, in nondimensional units, the pulse initial condition is given by a unit-width rectangle centered at the origin,

which implies p * (c * ; 0) = δ(c * -1). The batch concentration c B (t) for the well-mixed problem is nondimensionalized as above, and c B * (0) = 1.

E Batch dynamics

Here, we provide some details on the equations governing the well-mixed batch problem discussed in

Multiplying through by c and integrating over c (using integration by parts on the right hand side), we recover the standard well-mixed rate law for the batch concentration as a function of time,

Substituting Eq. ( 76) for the rate yields Eq. [START_REF] Bellin | On the local concentration probability density function of solutes reacting upon mixing[END_REF].

Once c B (t) drops below c m , at some time t m , the domain Ω(t) where concentrations are above this detection limit becomes empty, and the concentration PDF becomes ill-defined. By convention, we can set c B (t > t m ) = 0 and p(c; t > t m ) = δ(c), which conveys the meaning that concentrations are zero everywhere (below the detection limit).