hagués arribat fins aquest punt.

Agreement No. 617511 (MHetScale).

First, I would like to thank my two advisors. To Philippe, thanks for making me discover the applied sciences, the research world and opening new horizons. I have appreciated every piece of advice and every hint that you have given me during these three years. Cette aventure a commencé grâce à toi, et je te dois en grande partie d'être arrivé jusqu'ici. Merci pour tout.

To Marco, thank you for believing in me since the first instant, offering me a shot and welcoming me like family. I am grateful for these three years, for the work done together and the advances made. Thanks for all the professional (and personal) advices. I have learnt a lot and I think I am much wiser than I was three years ago. I owe it to you. I would like to thank my committee for accepting to assist and to revise the manuscript. Thank you for accepting the invitation so promptly, I look forward to the thesis defense and the discussions that will follow. To Alberto Guadagnini, Philippe Ackerer, Séverin Pistre, Michel Quintard, Insa Neuweiler and Markus Holzner, thank you.

I also want to thank my mid-term thesis committee members, thanks to Pietro De Anna, Bernard Célérier, and Pierre Camps for making sure that everything was going the right way and for the insightful conversations. A special thank you to you Pietro: you have always closely watched my work and given me useful advices.

Merci à l'équipe TMP ainsi que tout le personnel de la MSE, pour m'avoir accueilli et pour tous ces bons moments. Je pense notamment à Gérard, Richard, Olivier, Maria, Francesca, Filip, Anna, Sofia, Delphine, Roger, Massimiliano, Marie, ainsi que Olivier et Valentin de Voxaya. Un merci spécial à Anna et Filip avec qui j'ai passé des moments mémorables qui m'ont particulièrement aidé au début de la thèse et sans qui, le commencement eut été beaucoup plus ardu. I gràcies a la Maria and thanks to Francesca per les pauses and the coffee breaks.

Quiero agradecer a toda la gente del CSIC de Barcelona por haberme acogido como uno de los suyos. Muchas gracias a Silvia, Jordi B., Quim, Alba, Berta, Yoar, Sheila, Ester, Linda, Laura S., Laura M., Carme, M.A., Francesco, Jordi P., Cristina, Rotman, Cristian, Lucía. Y a los que llegaron después os agradezco también todos los momentos pasados, Robert, Max, Lidia, Tybaud, Jing-Jing, Andrea, Sonia, Nacho, Ahmad, Ashkan, Sandra, Elina, Lucía. Y agradezco también a la gente de la UPC y especialmente a Michela, las dos Sònias, Guillem, Mireia, Carme, Oriol y Nuria.

I also want to thank more specifically our group. Thanks to Alessandro, Vivien, Tomás, John, Ran, Laz, Kevin, Maria, Filippo and Juan, for all the discussions, the meetings and the very good times. I have enjoyed belonging to our group.

Un merci spécial à Vivien avec qui j'ai partagé pendant longtemps le bureau, pour de nombreuses conversations et quelques fou-rires. Gracias también a Ale con quien he compartido buenos momentos (lo siento no haber venido más a menudo a merendar!).

I would also like to thank my very good friends here in Barcelona. Thanks to Tomás, Laz, Filippo, Lidia, Michela. To Tomás: we have arrived together in that city and have been very good friends since. Thanks so much for the hand during the last lap of the race. Thanks for the all the good times and the support, appreciate it. I'll just slip a small "tiiinyy Riiiiiiiick" in! A Laz, tal vez la persona que conozco desde hace más tiempo aquí y uno de mis mejores amigos, gracias por todos los momentos en el trabajo y fuera del trabajo. No olvidaré los momentos en Córcega, Génova, Viena, Rotterdam, Washington, Filadelfia y Nueva York. ¡Muchísimas gracias por todo y por el apoyo! ("Here with my homies, just writing the acknowlegments..."). To Filippo, lots of thanks, I didn't think I would lose you to Switzerland! Thank you for all the good moments, I don't think there's a single person that could be unhappy when you are in the area. A warm hug to Valentina and the little Nina! Still can't believe you're a dad now! A la Lídia, moltíssimes gràcies per tot. Per la companyia al despatx, per les converses gracioses i pels moments una mica més durs... Merci pour le soutien durant ces deux dernières semaines qui ont été particulièrement compliquées.

Y gracias por el ambiente más relajado en el despacho en esta época estresante.

A Miki, ¡muchas gracias por tu presencia aunque no nos hemos podido ver tanto por varias razones! ¡Gracias por sacarme a tomar el aire durante este último sprint que ha sido escribir la tesis! ¡¿Ahora a disfrutar, no?! ¡Con toda la C-family! I also want to thank Sebastian who has become a good friend of mine and who offered me his help when I needed it Most. Je veux remercier également mes amis de Montpellier, je pense notamment à Polo, Cynthia, Anais, Abou, John, Yohan, Chloé, Julien.

Je veux aussi remercier Ad, Al, Olf, Miche, Swan et Rita, certains de mes plus proches amis, pour les bons moments passés ensembles. Juste ces quelques mots. À Al « heureusement que j't'ai toi ». À Flo « quelle prochaine paire de shoes on s'achète? ». À Ad « à toutes les nuits blanches » et une petite line : « cette société n'est qu'une enclume, j'ai couru après le fric quitte à y laisser ma plume ». À Swan, j'ai hâte de repasser à MTP et qu'on se voit. À Rita, merci, c'est sans doute grâce à toi que j'ai essayé une voie plus appliquée. Et Michele, tu seras toujours notre italien préféré! Et pour finir : « I don't have friends, I got family ».

Je veux aussi remercier des amis proches, certains que je n'ai pas vu depuis des années, d'autres que je vois plus souvent, mais que je n'oublie pas. Je pense à Margaux, Amy, Adrien, Tomas, Damien, Amandine, Maria, Julie.

Moltíssimes gràcies a l'Antònia, al Josep Lluís, a l'Anna, al Bernat, a la Irene, la Montse, i la Carme. Em vau acollir des del principi com un altre membre de la família i us ho agraeixo molt. Em vau recolzar moltíssim durant tota l'època de tesi. A l'Anna i al Bernat: sense vosaltres aquesta tesi no hagués quedat tan bonica, gràcies! Je remercie également ma famille qui a toujours été là pour moi malgré la distance. À mon père, merci d'avoir été là, merci d'être venu ici à Barcelone et de m'avoir soutenu durant ces trois ans. À ma grand-mère Marguerite, qui malgré la distance est toujours dans mon coeur. À Frère: tu as toujours été mon modèle, un exemple d'acharnement et toujours les émotions à fleur de peau, les bonnes comme les plus tristes. Les mots me manquent mais je peux probablement résumer tout ça avec un « on va donner 24 heures de répit à l'humanité ». À Soeur, merci pour tout, tout le soutien, toujours. Je sais que je peux toujours compter sur toi, pour tout type de situation. Merci pour le soutien malgré la distance, je n'y serais pas arrivé sans toi. Et je sais que toi tu réussiras, ils peuvent te planter dans le dos y'aura pas d'sang, parce que c'est du son qui coule dans tes veines, en BPM. Et quelles que soient les épreuves, on les affronte ensemble parce qu'on est QLF. À ma mère, merci d'avoir toujours été là pour moi. La réalisation et l'écriture de cette thèse n'aurait jamais été possible sans toi. Merci pour le soutien, toujours, merci de toujours être fière de moi, de respecter mes choix, d'être là quand j'en ai besoin. Regarde jusqu'où on est arrivés, comme quoi au final, la persévérance a payé. Merci Mam's, de m'avoir guidé à chacun de mes pas. Si je suis la personne que je suis aujourd'hui, c'est à toi que je le dois. À ma grand-mère Francette, qui a toujours été fière de moi et qui me manque énormément.

Résumé

Les mécanismes responsables du transport hydrodynamique anormal (nonfickéen) peuvent être rattachés à la complexité de la géométrie du milieu à l'échelle des pores. Dans cette thèse, nous étudions la dynamique des vitesses de particules à l'échelle des pores. À l'aide de simulations de suivi de particules effectuées sur un échantillon numérisé de grès de Berea, nous présentons une analyse détaillée de l'évolution des statististiques de vitesse lagrangienne et eulérienne et de leur dépendance aux conditions initiales. Le long de leur ligne de courant, les particules voient leur vitesses évoluer temporellement d'une manière complexe, appelée intermittence, alors que leur série de vitesses spatiales présente, elle, des fluctuations régulières. La distribution spatiale des vitesses des particules converge rapidement vers l'état stationnaire. Ces résultats conduisent à des processus markoviens qui permettent de prédire les fluctuations de vitesse dans le réseau poral. Ces processus, associés à la tortuosité et à la distance de corrélation de vitesse, permettent de paramétrer un modèle de marche aléatoire dans le temps (CTRW) et de réaliser le changement d'échelle pour simuler le transport à l'échelle de Darcy. Le modèle, comme tout modèle issu d'un changement d'échelle, repose sur la définition d'un volume élémentaire représentatif (VER). Nous montrons qu'un VER basé sur les statistiques de vitesse permet de définir un support pertinent pour la modélisation du transport hydrodynamique pré-asymptotique à asymptotique, et ainsi d'éviter les limitations associées à l'équation d'advection-dispersion fickéenne. Enfin, nous étudions l'impact de l'hétérogénéité du réseau poral sur le volume de mélange et la masse du produit d'une réaction bimoléculaire. .5 Spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the fluxweighted initial particle distribution p 0 (v) to the steady s-Lagrangian PDF P(v) at different distances. And, spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the uniform initial particle distribution p 0 (v) to the steady s-Lagrangian PDF P(v) at growing distances. . . . . . . . . . . . . . . . .

3.6

Temporal evolution of the ensemble t-Lagrangian PDF P(v,t) from the flux-weighted initial particle distribution p 0 (v) to the steady t-Lagrangian PDF P(v) at different times. And, temporal evolution of the ensemble t-Lagrangian PDF P(v,t) from the uniform initial particle distribution p 0 (v) to the steady t-Lagrangian PDF P(v) at growing times. . . . . . . . . . . . 3.8 Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for direct particle tracking and Markov model simulations from the uniform initial particle distribution p 0 (v) to the steady s-Lagrangian PDF P(v) at distances s = 0, 8/3, 16/3, 8 p . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Evolution of the mean s-Lagrangian velocity v s (s) for a uniform initial distribution and expression (3.33) with c = 2.5 p . . . . . . . . . . . . . . . 3.10 Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for both direct particle tracking and Bernoulli process simulations from the uniform initial particle distribution p 0 (v) to the steady s-Lagrangian PDF P(v) for distance s = 0, 8/3, 10 p . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for direct particle tracking and Ornstein-Uhlenbeck process simulations from the uniform initial particle distribution p 0 (v) to the steady s-Lagrangian PDF P(v) (red) at distance s = 0, 4/3, 4, 8, 10 p . . . . . . . . . . . . . . . . . . in the first porous medium (symbols), and from the hydrodynamic dispersion coefficient (black solid line). The analytical solution overpredicts the product total mass in the medium due to the assumption of complete mixing between reactants. Note that similar behaviors have been observed by [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF] The prediction and proper understanding of the hydrodynamics of transport in porous media play an important role in numerous environmental and industrial applications. These include groundwater [START_REF] Fu | The use of zero-valent iron for groundwater remediation and wastewater treatment: a review[END_REF] and soil remediation [START_REF] Yoshida | Land-surface contamination by radionuclides from the Fukushima Daiichi nuclear power plant accident[END_REF][START_REF] Rabiet | Positive gadolinium anomalies in wastewater treatment plant effluents and aquatic environment in the Hérault watershed (south France)[END_REF], geothermal energy and petroleum production, nuclear waste disposal [START_REF] Ewing | Nuclear waste disposal-pyrochlore (a 2 b 2 o 7): Nuclear waste form for the immobilization of plutonium and "minor" actinides[END_REF], geologic CO 2 storage [START_REF] Szulczewski | Lifetime of carbon capture and storage as a climate-change mitigation technology[END_REF], and sustainable exploitation of groundwater [START_REF] Harvey | Arsenic mobility and groundwater extraction in bangladesh[END_REF]. Geological and engineered porous media are characterized by spatial heterogeneity which makes the understanding of flow, transport and reaction processes across scales a challenge. This thesis focuses on upscaling these processes from the pore to the Darcy scale.

The transport of mass through porous media has traditionally been modeled using a Fickian transport paradigm that models the impact of pore-scale velocity fluctuations in terms of hydrodynamic dispersion [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. Thus, Darcy scale transport is described by an advection-dispersion equation (ADE) characterized by the average pore velocity and dispersion. This framework predicts that localized solute injections evolve to a Gaussian-shaped distributions whose second and second-centered moments increase linearly with time, breakthrough curves are characterized by inverse Gaussian-shaped distributions, which decay rapidly after the mean solution arrival. However, non-Fickian (also called anomalous) transport behaviors have been observed in disordered media ranging from geological media [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]Neuman and Tartakovsky, 2009b), to biological tissues [START_REF] Guigas | Sampling the cell with anomalous diffusion -the discovery of slowness[END_REF], and amorphous semiconductors [START_REF] Scher | Anomalous transit-time dispersion in amorphous solids[END_REF]. In the case of groundwater, Fickian-breaking features have been observed across the scales: at the field scale [START_REF] Sudicky | A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process[END_REF][START_REF] Garabedian | Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts: 2. analysis of spatial moments for a nonreactive tracer[END_REF][START_REF] Le Borgne | Non-fickian dispersion in porous media: 2. model validation from measurements at different scales[END_REF][START_REF] Boggs | Transport of tritium and four organic compounds during a natural-gradient experiment (made-2)[END_REF] and Darcy scale [START_REF] Koch | Anomalous diffusion in heterogeneous porous media[END_REF][START_REF] Levy | Measurement and analysis of non-Fickian dispersion in heterogeneous porous media[END_REF], but also at the pore scale (Cortis et al., 2004;[START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF]. A Vertically-integrated concentration. The solute plume deviates from Fickian behavior with early breakthrough and late tailing. This is reflected in the non-Gaussian vertically averaged concentration.

non-exhaustive list of the classical signatures of anomalous transport includes non-linear scaling of the mean squared displacement, early breakthrough and late tailing of solute, non-Gaussian features such as tailed plume shapes (see Figure 1.1), and heavy-tailed arrival time distributions [START_REF] Bouchaud | Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]Neuman and Tartakovsky, 2009b). These behaviors, arising from complex mechanisms, render the modeling of transport in heterogeneous media challenging. At the millimeter scale, flow and transport can be fully resolved through direct simulations and the mechanisms involved are fully understood. However, obtaining predictions at larger, relevant scales is often much more challenging. This complexity lies first in the impossibility of obtaining exhaustive knowledge of cross-scale geometrical features, and second in the absence of sufficiently powerful computational resources for performing flow and transport simulations on very large, resolved domains. In addition, even in an hypothetical setting where an accurate micron-scale description of rock architecture for cubic-meters scale would be achievable, such direct simulations would provide little insight into the actual mechanisms that are responsible for the large-scale behavior of solute.

Thus, if we wish to model and understand the transport of solutes and particulates such as contaminants, colloids, or viruses and microbial agents, we must develop a methodology that correctly accounts for spatial heterogeneity. A further issue is that, in most situations, the architecture of the medium is not accurately known. This is especially true for subsurface flows where obtaining precise information on the geometry of the rock through experiments is prohibitively complicated. To tackle these issues, effective parameters such as permeability or porosity are defined as appropriate averages over a volume. This leads directly to the concept of a representative elementary volume (REV). The REV is defined as the minimum volume of a sample over which the measured property of interest becomes constant and independent of the size of the domain (see, for example, [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Yang | Upscaling for adiabatic solid-fluid reactions in porous medium using a volume averaging theory[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF]. It is important to note that the scale at which the REV property is defined should be related to the scale at which the information is used. For instance, permeability defined over a geological formation will differ strongly from that obtained in the laboratory [START_REF] Bear | Dynamics of fluids in porous media[END_REF].

The previous considerations lead naturally to the notion of upscaling. One is often interested in large-scale behavior emerging from small-scale mechanisms; therefore, the concept behind upscaling techniques is to develop methodologies that predict the behavior of the transport at the scale of interest, while taking into account in an efficient manner the subscale complexity. Due on the one hand to the limited information available for the characterization of geological media, and on the other to the goal of describing largescale features efficiently, the upscaling technique may involve the use of probabilistic tools, with features such as medium heterogeneity being described effectively through statistical properties.

The traditional Eulerian description of Darcy scale solute transport is given by the advection-dispersion equation (ADE),

∂ ∂t φ c(x x x,t) + ∇ • φ q q qc(x x x,t) -∇ • φ D∇c(x x x,t) = 0, (1.1)
where c(x x x,t) is the solute concentration at position x x x at time t, where φ is the porosity, q q q is the Darcy velocity, which describes the drift rate of the plume, and D is the dispersion tensor. For isotropic dispersion, which occurs in particular in homogeneous media, the latter is diagonal with nonzero elements equal to the dispersion coefficient D, but it may also account for anisotropic dispersion along the flow and transverse directions [START_REF] Bear | On the tensor form of dispersion in porous media[END_REF].

This equation describes mass balance between REVs, with concentration change (the first term) being impacted by advective flux (second term) and Fickian dispersion (third term).

It is valid at scales much larger than the REV. Since the REV is usually defined for the porosity, a standard assumption is to consider that this choice of REV holds for advection and dispersion properties [START_REF] Whitaker | The Method of Volume Averaging[END_REF].

The validity of Fickian dispersion is linked to the central limit theorem. In a Lagrangian sense, if all solute parcels sample the full velocity heterogeneity of the medium, the resulting spreading converges, under certain conditions, to a normal distribution [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF]. In practice, however, this convergence may be jeopardized at the scale of interest by 1.1. STATE OF THE ART 7 a number of mechanisms which may undermine the validity of the classical REV definition. In the presence of preferential flow paths, such as in fractured media, communication between distant points in space may occur on a short timescale. Conversely, stagnation zones such as dead-end pore spaces or low-permeability inclusions may lead to anomalously long retention times. Upscaled models typically represent these effects through nonlocal features. Local models (such as the ADE) allow for local solute exchange only.

In other words, mass transfer described by differential operators, or equivalently to Eulerian numerical models where only nearby nodes communicate. This holds also for the temporal behavior, in the sense that there is no long-term memory of past plume history.

In contrast, spatially nonlocal models allow for communication between locations separated by large distances, while temporally nonlocal models allow for long-range memory in time (see Figure 1.2 for an example of non-local model). Even in the absence of such features, the validity of Fick's law requires a time τ D = 2 c /D, where D is the dispersion coefficient and c is the characteristic length of the REV. This dispersion timescale represents the time needed for diffusion to homogenize solute within the REV, so that all the corresponding heterogeneity is sampled. At times t < τ D (called pre-asymptotic times), the Fickian paradigm typically breaks down because pore-scale variability is not fully sampled and thus not averaged out. We briefly introduce three of the most commonly used nonlocal models for predicting transport in heterogeneous media: fractional advection dispersion equations (fADEs), multi-rate mass transfer (MRMT) and continuous time random walks (CTRWs). In this thesis, we will focus on continuous time random walk models. fADEs are generalizations of the ADE allowing for non-local features. Rather than the short-range temporal or spatial correlations implicit in the ADE, which arise from the conditions of the classical central limit theorem and lead to Gaussian-shaped plumes and regular derivative operators, fADE models describe long-range correlations through fractional derivatives. These are integral operators which describe long-range memory or spatial correlations through heavy-tailed memory or spatial kernels, which reflect the generalized central limit theorem [START_REF] Meerschaert | Stochastic models for fractional calculus[END_REF]. This translates into solute being transported over anomalously large distances over short times (spatial fADEs), retained for anomalously long times (temporal fADEs), or a combination of both (spatiotemporal fADEs).

The multi-rate mass transfer (MRMT) approach models non-Fickian transport due to mass transfer between mobile and immobile continua [START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and[END_REF][START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF][START_REF] Geiger | A novel multi-rate dual-porosity model for improved simulation of fractured and multiporosity reservoirs[END_REF]. These could describe, for example, mass transfer between low-permeability inclusions and the highly conductive background in groundwater flow and mass transfer between fracture and matrix as well as between mobile and immobile microporosity on the pore scale [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]. Using this framework, complex features of solute retention times may be captured, which reflect themselves in anomalous transport features of the plume and breakthrough curves as discussed above.

Continuous time random walks (CTRWs) are a generalization of the classical random walk, first developed by [START_REF] Montroll | Random walks on lattices[END_REF], generalized by [START_REF] Scher | Stochastic transport in a disordered solid. I. Theory[END_REF], and later introduced in hydrogeological modeling, see [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF] for a review. The classical random walk approach, due to [START_REF] Einstein | On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat[END_REF] and [START_REF] Smoluchowski | Mathematical Theory of the Kinetics of the Coagulation of Colloidal Solutions[END_REF], describes particle trajectories in terms of stochastic recursion relations for the position x x x n and the time t n after n steps,

x x x n+1 = x x x n + η η η n , t n+1 = t n + ∆t. (1.2)
During the nth step, particle positions change by a constant distance along the vector η η η n , and the constant duration of a step is ∆t. Along each spatial dimension, the components of η η η n are equal to ±1 independently with probability 1/2. It is important to note that, while in the original models transported particles were physical entities (such as pollen in water for classical Brownian motion [START_REF] Einstein | On the motion of small particles suspended in liquids at rest required by the molecular-kinetic theory of heat[END_REF]), they are now typically employed as abstract Lagrangian entities representing discretized solute mass [START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF].

CTRW models allow for stochastic displacements and jump durations, and generalize this description as

x x x n+1 = x x x n + ξ ξ ξ n , t n+1 = t n + τ n , (1.3)

STATE OF THE ART

where the ξ ξ ξ n and τ n are independent and identically distributed for each step n and sampled from a joint distribution. The CTRW framework thus aims to describe the mechanisms leading to anomalous transport in terms of the statistics of particle displacements and their (potentially coupled) durations. We note also that fADEs and MRMT models have been shown to be equivalent to CTRW for certain choices of the joint distribution of displacement and duration [START_REF] Schumer | Fractal mobile/immobile solute transport[END_REF][START_REF] Dentz | Transport behavior of a passive solute in continuous time random walks and multirate mass transfer[END_REF][START_REF] Benson | Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes[END_REF]. A one-dimensional illustrative example of a one-dimensional CTRW is given in Figure 1.2. For further details on the CTRW framework, the reader is referred to Section 2.3.

Note that Fickian transport may, however, be recovered at large temporal and spatial scales. This occurs when particles have experienced the full velocity heterogeneity, and corresponds to times that are larger than the memory of the initial conditions [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF]. To better understand the pore-scale mechanisms (see Figure 1.3) leading to this phenomenon, many authors have investigated and modeled the Lagrangian velocities experienced by the solute particles. de Josselin de [START_REF] De Josselin De Jong | Longitudinal and transverse diffusion in granular deposits[END_REF] and [START_REF] Saffman | A theory of dispersion in a porous medium[END_REF] were the first to consider the concept of residence time of particles in the modeling of solute motion.

The models they considered were similar to the CTRW introduced above in the sense that particle displacements over characteristic distances were coupled to random increments related to pore-scale velocity distributions. More recently, [START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF] (at the pore-scale) and [START_REF] Meyer | Particle-based transport model with Markovian velocity processes for tracer dispersion in highly heterogeneous porous media[END_REF] (at the Darcy-scale) studied the temporal Lagrangian velocity series and modeled them with non-linear Langevin equations. However, the complexity of the temporal behavior of the particles velocities requires the introdution of velocity dependent parameters to quantify the evolution of the particle velocities. The need for this level of complexity is tied to the intermittency exhibited by the temporal Lagrangian velocity series: the temporal velocity evolution of a particle exhibits long periods of low velocity punctuated by short peaks of high velocity [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF]. These behaviors are directly linked to the spatial organization of the velocities: particles cross high velocity regions rapidly while spending more time in low velocity areas. This wide range of transit times is at the root of anomalous transport features [START_REF] Carrel | Pore-scale hydrodynamics in a progressively bioclogged three-dimensional porous medium: 3-d particle tracking experiments and stochastic transport modeling[END_REF][START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF]. However, [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] have shown that spatial velocity series do not display such intermittency and rather exhibit velocity changes with a regular frequency. The corresponding spatial velocity persistence, arising from the pore structure of the medium, has lead to the investigation of the spatial evolution of the Lagrangian velocity statistics in the framework of CTRWs with constant spatial increments [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]. CTRW approaches can naturally represent this type of behavior with a waiting time distribution that accounts for the spatial correlation scale. A particle (blue) experiences a preferential channel of low tortuosity, whereas another particle (red) experience a more tortuous pathway. Finally, the last particle particle (green) is trapped in a region of low velocity and will need a large time to escape. These mechanisms, together with the heterogeneity of the velocity field they experience, enhance anomalous transport.

Understanding these phenomena requires an accurate characterization and understanding of the dynamics of Eulerian and Lagrangian velocity statistics, leading to a series of recent studies [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF]Siena et al., 2014;[START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF][START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF]Jin et al., 2016a;[START_REF] Matyka | Power-exponential velocity distributions in disordered porous media[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF].

Together with these concepts, it is important to consider the role of the initial velocity distribution of solute particles. While velocity statistics obtained in numerical simulations and experiments have often been assumed to be stationary in time, the injection mode and volume size have been shown to impact the early-time velocity statistics [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF].

In some cases, the injection velocity distribution largely differs from the steady-state, which can affect the pre-asymptotic dynamics of the solute velocity distribution. This effect has been studied for Darcy scale transport in fractured and porous media [START_REF] Hyman | Influence of injection mode on transport properties in kilometer-scale three-dimensional discrete fracture networks[END_REF][START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF][START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF].

Open Questions

The understanding and modeling of transport in heterogeneous media have recently seen significant progress. However, there remain a large number of open issues. Some of the

OBJECTIVES

key open problems refer to:

1. the relation between temporally-and spatially-sampled pore scale velocities, together with their connection to the Eulerian statistics;

2. the modeling of spatial velocity series experienced by solute particles;

3. the upscaling of purely advective, pre-asymptotic transport through continuous time random walks;

4. the definition of a representative elementary volume allowing for the upscaling of anomalous transport;

5. the impact of pore-scale heterogeneity on the mixing volumes between reactants and the resulting reactive transport dynamics.

Objectives

The two Main objectives of this thesis are:

• To understand and characterize the impact of pore-scale velocity heterogeneity on spreading and mixing of the solute.

• To develop an effective model for the upscaling of transport from the pore to the Darcy scale.

These two main objectives can be broken down into Specific objectives towards answering the Open questions. These objectives are:

1. to mathematically relate the equidistantly and isochronously sampled velocity statistics along a particle streamline and to validate the relation via pore-scale direct simulations on digitized real rock samples;

2. to develop stochastic processes for modeling particle spatial velocity series;

3. to identify the spatial information necessary for the parameterization of CTRW type upscaling models;

4. to derive a CTRW for the upscaling of purely advective preasymptotic transport;

5. to define a representative elementary volume in terms of the Eulerian and Lagrangian velocity statistics for the upscaling of anomalous transport;

6. to investigate the behavior of mixing and reaction in media displaying high pore-scale heterogeneity.

CHAPTER 1. INTRODUCTION

Thesis Structure

This thesis is structured in seven chapters, which are summarized below. All the chapters containing the research performed towards this thesis are presented in the format of scientific publications. Such papers are either published (Chapter 3), submitted to peer reviewed journals (Chapter 4), or to be submitted (Chapters 5 and 6). These chapters can be read separately and may therefore contain repeated concepts. Nevertheless, all are interconnected by the two principal objectives of this thesis, see Section 1.3.

• In Chapter 2, a description of the methodology used in this thesis for flow and transport is provided. First, we detail the image acquisition technique and the successive steps involved in solving the flow and determining the Eulerian velocity field. Then, we present the implementation of the different mechanisms of a particle tracking algorithm: the streamline tracing method, the diffusion process, and the computation of reaction. Finally, we introduce the CTRW framework that will be used for the upscaling of transport in the following chapters.

• Chapter 3 presents a detailed analysis of the evolution of Lagrangian and Eulerian statistics and their dependence on the injection condition, with the aim of deriving a framework for the efficient upscaling of transport from the pore to the Darcy scale.

The study presented here is based on velocity data obtained from computational fluid dynamics simulations of Stokes flow and advective particle tracking in the threedimensional pore structure obtained from high resolution X-ray microtomography of a Berea sandstone sample. While isochronously sampled velocity series show intermittent behavior, equidistant series vary in a short-correlation pattern. Both statistics evolve toward stationary states, which are related to the Eulerian velocity statistics. The equidistantly sampled Lagrangian velocity distribution converges within only a few pore lengths. These findings indicate that the equidistant velocity series can be represented by an ergodic Markov process, and a stochastic Markov model for the equidistant velocity magnitudes captures the evolution of the Lagrangian velocity statistics. The model is parameterized by the Eulerian velocity distribution and a relaxation length scale, which can be related to hydraulic properties and the medium geometry. These findings lay the foundation for a predictive stochastic approach to upscale solute dispersion in complex porous media from the pore to the Darcy scale.

• In Chapter 4, we study the upscaling of purely advective pore-scale dispersion in terms of the Eulerian velocity distribution and advective tortuosity, both flow attributes, and of the average pore length, a medium attribute. The stochastic particle • In Chapter 5, we define REVs in terms of the statistics of the medium velocity process. In order to achieve this, we propose two types of REV definitions that are based on Eulerian and Lagrangian statistics respectively, which we term velocity REVs (v-REVs). According to the first definition, a volume is considered an Eulerian velocity REV if it is sufficiently large for the Eulerian velocity PDF to be constant and not to evolve anymore with the size of the domain. The second definition is the Lagrangian counterpart of the first definition. In order for a volume to be a Lagrangian velocity REV, any initial velocity distribution defined therein must evolve and reach a steady state before trajectories leave the volume. We then show that, when combined, these definitions allow for accurate upscaling via the CTRW framework. To illustrate and test these definitions and their consequences, we use a mm 3 sample of a Berea sandstone as a potential v-REV candidate. The sample, which is an REV for porosity, happens to fulfill the Eulerian v-REV definition since the Eulerian velocity distribution converges to the steady state within 1/8 of the total sample size. We also observe that a randomly selected velocity distribution at the inlet leads to an already stationary velocity distribution after a distance of about two thirds of the sample, guaranteeing that the sample is also a Lagrangian REV. Then, we upscale transport through a one-dimensional CTRW model based on these REV properties. Finally, we discuss the impact of such REV definitions on the upscaling of both purely advective and advective-diffusive transport to the Darcy scale through the CTRW transport description.

• In Chapter 6, we simulate an instantaneous bimolecular reaction, A + B → C, in an heterogeneous two dimensional porous domain. A particles are injected at the inlet of the generated medium, while B particles occupy the rest of the domain. The reaction CHAPTER 1. INTRODUCTION occurs at the interface. We first solve the flow, before performing the reactive transport simulations. These simulations are performed with a reactive particle tracking algorithm that sequentially applies the advective, diffusive, and reactive steps. Simultaneously, we discuss the dispersive lamella methodology of [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF], which is based on the computation of the width of the mixing front for the prediction of the overall product mass. We show that this methodology can be extended to complex pore-scale transport systems in heterogeneous porous media. We also discuss the impact of the porosity experienced by the fluid plume when it is not representative of the average porosity of the full sample. Lastly, we validate the predictions of the dispersive lamella methodology against the reactive particle tracking results.

• Chapter 7 presents the main conclusions of the thesis and discusses future work.

• At the end of the document, three Appendix sections present derivations that are not detailed in the main text and additional material concerning the valorization of this work. Au chapitre 5, nous proposons une définition de volume élémentaire représentatif (VER) en termes des statistiques des vitesses contenues dans le milieu. Pour y parvenir, nous proposons deux types de définitions de VER basés respectivement sur les statistiques eulériennes et lagrangiennes. Nous les appelons VERs de vitesse (VERV). Selon la première définition, un volume est considéré comme un VERV eulérien s'il est suffisamment grand pour que sa distribution eulérienne des vitesses soit constante et ne puisse plus évoluer en fonction de la taille du domaine. La deuxième définition est l'homologue lagrangien de la première définition. Pour qu'un volume soit un VERV lagrangien, toute distribution de vitesse initiale définie dans ce VERV doit évoluer et atteindre un état stationnaire avant que les trajectoires n'atteignent la fin du domaine. Nous montrons ensuite que, lorsqu'elles sont combinées, ces définitions permettent une modélisation à plus grande échelle précise grâce à des modèles de type CTRW. Pour illustrer et expérimenter ces définitions et leurs conséquences, nous utilisons le même échantillon de grès de Berea comme candidat potentiel de VERV. L'échantillon, qui est un VER pour la porosité, satisfait la définition eulérienne du VERV car la distribution de vitesse eulérienne converge vers l'état stationnaire après seulement 1/8 de la taille totale de l'échantillon. Nous observons également qu'une distribution de vitesse sélectionnée aléatoirement à l'entrée du domaine converge vers une distribution de vitesse déjà stationnaire après une distance d'environ deux tiers de l'échantillon.

Résumé Etendu -Français

Ceci garantit que l'échantillon est également un VERV lagrangien. Nous utilisons alors un modèle CTRW unidimensionnel basé sur les propriétés de VERV pour prédire le transport hydrodynamique. Enfin, nous discutons l'impact de ces nouvelles définitions sur le changement d'échelle de problèmes de transport purement advectif et advectif-diffusif à l'échelle de Darcy.

Au chapitre 6, nous simulons une réaction bimoléculaire instantanée, deux réactifs A + B réagissent pour produire C, dans un domaine poreux hétérogène en deux dimensions.

Des particules A sont injectées à l'entrée du support généré, tandis que des particules B occupent le reste du domaine. La réaction se produit à l'interface. Nous résolvons d'abord l'écoulement, avant de réaliser les simulations de transport réactif. Ces simulations sont effectuées avec un algorithme de suivi de particules réactif qui effectue de manière séquentielle les étapes d'advection, de diffusion et de réaction. Parallèlement, nous utilisons une méthodologie appelée dispersion lamellaire qui repose sur le calcul de la largeur du volume de mélange pour la prédiction de la masse globale de réaction chimique produite. Nous montrons que cette méthodologie peut être étendue à des systèmes complexes où le transport s'effectue dans des milieux poreux hétérogènes. Nous étudions également l'impact de la porosité sur les mesures de mélange lorsque le soluté n'échantillonne pas (ou pas encore) toute la porosité de l'échantillon. Enfin, nous testons les prédictions de cette méthode avec les résultats d'une modélisation par suivi de particules réactif.

En conclusion, nous avons fourni une analyse complète de la dynamique stochastique des vitesses de particules de soluté à l'échelle des pores.

CHAPTER 2

Methodology for Flow and Transport

Chapter 2

Methodology for Flow and Transport

As stated in the introduction section, our study goes through the performance of flow and transport simulations at pore-scale on both real and artificial medium samples. Therefore, this section aims at giving a comprehensive overview of image acquisition, flow resolution and tool development for transport simulations. However, it is important to mention that for easy readability, most of the technical details of the methodology are not provided in the current section and are presented in Appendix 8.1. This chapter is organized in the following minor sections. Section 2.1.1 provides details of image acquisition. In Section 2.1.2, the steps performed to solve flow and obtain velocity fields are described. Section 2.2 focuses on the methodology we use for the computation of the transport. The accurate descriptions of advection, diffusion and reaction processes are provided in Sections 2.2.1, 2.2.2 and, 2.2.3 respectively. Finally, in Section 2.3 we introduce the continuous time random walk (CTRW) framework that will be used for the upscaling of the transport in the following chapters.

Flow

In this section, the steps required for calculating the velocity fields of the samples under consideration are presented. First, the 3D and 2D geometries used in this thesis for flow and transport simulations are described, and some particularities regarding their acquisition are given. Then, the most important details of the flow solving are provided.

Geometries

In order to perform the flow and transport simulations described in Chapters 3 and 5, a digitized and segmented Berea sandstone sample extracted from a complete core was selected. The three-dimensional geometry of this sample is displayed in

Solving the Flow

The details of the flow simulations on both 3D and 2D cases are available in Sections 3 and 6, respectively. Therefore, to avoid redundancies, just a brief description of the flow simulations is hereby provided simply to highlight most crucial details.

The 3D digitized Berea sandstone sample and the 2D artificially generated geometry displayed in Figure 2.1 were fine-meshed with cubes and squares respectively to capture accurately the pore space. Then we solved the Stokes and continuity equations (that dictate fluid flows under low Reynolds number) for constant viscosity and density. The system of equations is solved via a finite-volume scheme implemented in OpenFOAM [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. We prescribed pressure boundary conditions at the inlet and at the outlet of the domain and implemented no-slip conditions at the solid-void interfaces and at the remaining boundaries of the domain. The simulations result in flow fields where the velocity values are obtained at each interface of the mesh in the normal direction (see Figure 2.2 for an illustrative example). Once the velocity fields have been obtained, we then need a stream- line tracing method that tracks particles using the previously calculated Eulerian velocities at the interfaces. Such method is described in the following section.

Transport

Advective Particle Tracking

In this section we start by providing a brief overview of the existing particle tracking methods before going through the details of the algorithm we use in this report in Section 2.2.1.2.

Bases

Computing accurately the properties of the transport requires an effective particle tracking algorithm. Indeed, generating particle path lines and travel-time information through particle tracking techniques can be extremely handy for the analyses of complex two-and threedimensional groundwater flow systems. They are usually based on a previously obtained (usually through numerical flow simulations) Eulerian velocity field on the geometry of interest. These particle tracking schemes and methods have been incorporated into formal solute transport models to mimic the advective mechanism of the transport [START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF][START_REF] Most | Evolution and persistence of cross-directional statistical dependence during finite-Péclet transport through a real porous medium[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF][START_REF] Bijeljic | Pore-scale modeling of longitudinal dispersion[END_REF]. They can also be used to generate streamlines and tracking the solute as an ensemble of particles and therefore evaluate the advective characteristics of the transport [START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF].

The simplest particle tracking method probably consists in, first computing every particle's position velocity component and second, moving the particle to a new location by adding the displacement in each direction to the original position of the particle. The displacement in each direction is simply being computed by as ∆x = u∆t where u is the velocity component and ∆t is a previously elected time step. This can be performed for every component, resulting in a new location for the particle after time step ∆t. Unquestionably, while being conceptually easy to understand, these types of methods rely on a velocity field completely known and on the capability of the scheme to know the three velocity components at any of the domain locations. For analytical flow field solution, the velocity at every point of the field is given and therefore, no further calculations are required. However, numerical flow field solutions obtained from, for instance, finite differences, finite volumes, finite elements or even from their continuous and discontinuous Galerkin extension, lack some kind of velocity interpolation inside every mesh element.

The type of flow considered can also drastically impact on the choice of velocity interpolation. Numerous approaches exist for each case. For finite differences and volumes in case of steady flows and regular hexaedron meshing, three of the most used approaches might be the simple linear interpolation, the multilinear interpolation and the step function interpolation (see, [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF], for example). The idea behind the linear interpolation is that each velocity component varies linearly with distance between two interfaces of the mesh in its component direction. The multilinear interpolation acts in a similar way but assumes that each velocity component is a linear function of the three others' coordinate direction. The step function simply assumes that the velocity component is constant within a cell and changes abruptly at the next interface. These methodologies give valid results for steady flows computed on squared and cubic meshes. For unstructured quadrilateral and triangular meshes, extensions have been proposed. [START_REF] Matringe | Streamline tracing on general triangular or quadrilateral grids[END_REF], for example, developed a streamline technique based on the existence of a stream function. For more complex methodologies, such as Galerkin and discontinuous Galerkin, other methods exist (see, for example, [START_REF] Sun | Streamline tracing on unstructured grids[END_REF]. For non-steady flow, researchers have added a linear temporal dependence to account for the rapidly changing velocity field in transient flows, see [START_REF] Suk | Multidimensional finite-element particle tracking method for solving complex transient flow problems[END_REF][START_REF] Maier | An accurate method for transient particle tracking[END_REF].

As discussed in Section 2.1, we elected a finite volume scheme onto a regular hexaedron mesh for the computation of our steady-state flow. We therefore focused on the simple linear interpolation (see, for example, [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF]) that fits our approach and allows for CHAPTER 2. METHODOLOGY FOR FLOW AND TRANSPORT an analytical expression for particle trajectory within the cell. We discuss it in detail in the following section.

Streamline Tracing Algorithm

Pollock algorithm The particle tracking method we use in this thesis is based on a modified version of the Pollock algorithm for pathline reconstruction [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF]. In the current section we summarize the Pollock algorithm and highlight its main limitations. To sum up, from the finite volume scheme used, a continuity equation can be written in every cell. In what follows, the six cell faces are referred to as x 1 , x 2 , y 1 , y 2 , z 1 and z 2 , and the ∆x, ∆y and ∆z represent the length of a cell in the x, y and z directions respectively. Face x 1 is the face perpendicular to the x direction at x = x 1 . Similar definitions hold for the other five faces. The continuity equation in a cell reads as:

u 2 -u 1 ∆x + v 2 -v 1 ∆y + w 2 -w 1 ∆z = 0, (2.1)
where u 1 is the velocity component of the face perpendicular to the x-direction at position x 1 and the u 2 , v 1 , v 2 , w 1 and w 2 velocity components are defined similarly. Such average velocity components are obtained by dividing the volume flow rate across the face by the crossectional area of the face. For simplicity, we assume that there is no rate of created or consumed volume by internal wells or sinks. Then, to compute the streamlines, one needs a method to obtain the velocity at any point within this mesh cell. We use the simple linear interpolation introduced earlier to write:

u = u 2 -u 1 ∆x (x -x 1 ) + u 1 , (2.2) v = v 2 -v 1 ∆y (y -y 1 ) + v 1 , (2.3) w = w 2 -w 1 ∆z (z -z 1 ) + w 1 . (2.4)
This linear interpolation produces a continuous velocity vector field within each individual grid cell that identically satisfies the differential conservation of mass equation everywhere within the cell. The fact that the velocity vector field within each cell satisfies the differential mass balance equation is important because it assures that path lines will distribute water across the flow field in a way that is consistent with the overall movement of water across the system indicated by the solution of the finite-volume flow equations. Now, let us consider the movement of a particle through a three-dimensional finitedifference cell. The rate of change in the particle's x-component of velocity for instance as it moves through the cell is given by Equation (2.2). Interestingly, such rate of change 2.2. TRANSPORT 29 corresponds to the time derivative of the x-direction trajectory within a cell and can be rewritten as:

ẋ(t) = u 2 -u 1 ∆x (x(t) -x 1 ) + u 1 , (2.5)
which is equal to

ẋ(t) = u 2 -u 1 ∆x (x(t) -x 1 + u 1 ∆x u 2 -u 1 ), (2.6)
integrating the latter with respect to time gives

t 0 ẋ(t) x(t) -x 1 + u 1 ∆x u 2 -u 1 dt = t 0 u 2 -u 1 ∆x dt , (2.7) which is identical to ln x(t) -x 1 + u 1 ∆x u 2 -u 1 x 0 -x 1 + u 1 ∆x u 2 -u 1 = u 2 -u 1 ∆x t.
(2.8)

Finally, solving for x(t) gives the x-component trajectory equation:

x(t) = x 1 - u 1 ∆x u 2 -u 1 + u 1 ∆x u 2 -u 1 + (x p -x 1 ) e u 2 -u 1 ∆x t .
(2.9)

Solving similarly for the yand z-components gives the following set of equations

x(t) = x 1 - u 1 ∆x u 2 -u 1 + u 1 ∆x u 2 -u 1 + (x p -x 1 ) e u 2 -u 1 ∆x t ,
(2.10) (2.12) which completely describe the trajectory of a particle within a cell given any initial location (x p , y p , z p ). Now that the trajectory equations are known for any time t, the time t e needed for the particle to exit the cell needs to be computed. This can be done by solving this set of equations for t. However the exit interface for the particle needs first to be identified.

y(t) = y 1 - v 1 ∆y v 2 -v 1 + v 1 ∆y v 2 -v 1 + (y p -y 1 ) e v 2 -v 1 ∆y t , (2.11) z(t) = z 1 - w 1 ∆z w 2 -w 1 + u 1 ∆x w 2 -w 1 + (z p -z 1 ) e w 2 -w 1 ∆z t ,
Indeed, if we want to calculate for instance the exit time on the x-axis, we need to know if the particle will leave the voxel by the left side or by the right side. Using the Pollock algorithm, there are actually four possible situations. We give details for one of the possible situations while leaving the rest of them in the appendix, see Section 8.1.1. We start by considering a situation where velocities are positive on both faces. This is the situation displayed in Figure 2.3. This case is possibly the simplest of all since it is obvious that the particle will get out by the right side, through the plane x 2 . To compute the time needed to reach this point, we only need to replace x(t) by x 2 in Equation (2.10), and solving it for t we obtain

∆t x = ∆x u 2 -u 1 ln u 2 ∆x u 1 ∆x + (u 2 -u 1 )(x p -x 1 ) , (2.13)
where ∆t x is the time needed by the particle to get out of the cell by the plane x 2 . The three other cases are detailed in Section 8.1.1. After having identified the eventual exit interface in every direction and computed the times ∆t x , ∆t y , and ∆t z that the particle needs to exit the cell in the x-, yor z-direction, respectively, the actual time needed to exit can simply be computed as ∆t := min{∆t x , ∆t y , ∆t z }.

(2.14)

The exit coordinates can then be obtained by substituting t by ∆t:

x out = x(∆t) = x 1 - u 1 ∆x u 2 -u 1 + u 1 ∆x u 2 -u 1 + (x p -x 1 ) e u 2 -u 1 ∆x ∆t , (2.15)
y out = y(∆t) = y 1 - v 1 ∆y v 2 -v 1 + v 1 ∆y v 2 -v 1 + (y p -y 1 ) e v 2 -v 1 ∆y ∆t , (2.16
)

z out = z(∆t) = z 1 - w 1 ∆z w 2 -w 1 + u 1 ∆x w 2 -w 1 + (z p -z 1 ) e w 2 -w 1 ∆z ∆t .
(2.17)

With these three equations, we obtain both exit and arrival coordinates in the next cell. We then have to reiterate the process to get the streamlines (and the velocities) in the entire domain. Therefore, the whole process for one particle can be summarized with the following algorithm.

1-Read starting location.

2-Assign particle to a grid cell.

3-Compute cell face velocity components. [START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF] found that for porous media, the linear interpolation of the method does not account for the no-slip boundary conditions at the solid-void interfaces. This leads to high discrepancies in low porosity media such as the materials that we consider. In next section we discuss an extension of this algorithm that allows for a quadratic velocity interpolation close to the solid boundaries. This method allows for an accurate velocity interpolation even in the presence of no-slip boundary conditions.

Extended Pollock Algorithm

This algorithm aims at giving an accurate tracing technique even in the presence of solid boundaries. A classical textbook example to understand the need for a quadratic interpolation is the Hagen-Poiseuille velocity profile in a pore

u(y) = u 0 1 - y r 2 , (2.18)
where u 0 is the velocity in the middle of the channel. A quadratic interpolation is imposed in the direction perpendicular to the boundaries to take into account the no-slip boundary condition at the solid-void interfaces. [START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF] then proposed to extend the equations to all the boundary situations that can be encountered in a three-dimensional mesh. A mesh cell is surrounded by 6 other voxels. Considering that a cell can either be void or solid, there are 2 6 = 64 different possible situations. As an example, a void cell could be surrounded by one solid voxel on its left and 5 void voxels at the remaining locations. In this section we derive the equations for one of these situations. Derivations of the remaining possibilities can be found in the appendix, see Section 8.1.2. For cells that are only surrounded by void cells we use the standard Pollock algorithm introduced in the previous paragraph. When a voxel is surrounded by a solid voxel on its left and 5 void voxels, see Figure 2.6, we impose a quadratic interpolation in the normal direction to the rock interface as:

u(x) = u 1 ∆x 2 (x 2 -x) 2 . (2.19)
The velocity interpolation of the two remaining components can be obtained from the zero divergence criterion

∂ u(x, y) ∂ x + ∂ v(x, y) ∂ y + ∂ w(x, y) ∂ z = 0, (2.20)
see appendix, Section 8.1.2.1. We obtain the following quadratic interpolation for the remaining velocity components:

v = 2v 1 ∆x (x 2 -x) + 2(v 2 -v 1 ) ∆x∆y (x 2 -x)(y -y 1 ), (2.21) w = 2w 1 ∆x (x 2 -x) + 2(w 2 -w 1 ) ∆x∆z (x 2 -x)(z -z 1 ). (2.22)
This set of equations, together with Equation (2.19), fulfills the divergence criterion (2.20).

Then, we integrate the equation with respect to time (derivations are available in the appendix, see Section 8.1.2.1) to obtain the full trajectories:

x(t) = x 2 - ∆x 2 (x 2 -x p ) ∆x 2 + u 1 (x 2 -x p )t , (2.23) y(t) = y 1 - v 1 ∆y v 2 -v 1 + v 1 ∆y + (v 2 -v 1 )(y p -y 1 ) (v 2 -v 1 )∆y 1 + u 1 (x 2 -x p ) ∆x 2 t 2∆x(v 2 -v 1 ) u 1 ∆y - ∆x 2 u 1 (x 2 -x p ) ,
(2.24)

z(t) = z 1 - w 1 ∆z w 2 -w 1 + w 1 ∆z + (w 2 -w 1 )(z p -z 1 ) (w 2 -w 1 )∆z 1 + u 1 (x 2 -x p ) ∆x 2 t 2∆x(w 2 -w 1 ) u 1 ∆z - ∆x 2 u 1 (x 2 -x p )
.

(2.25)

Finally, solving these equations for t gives the time of flight of the particle to reach the exit of the cell as

∆t x = ∆x 2 u 1 1 x 2 -x 1 - 1 x 2 -x p ,
(2.26)

∆t y = ∆x 2 u 1 (x 2 -x p ) ∆y 2 v 2 v 1 ∆y + (v 2 -v 1 )(y p -y 1 ) u 1 ∆y 2∆x(v 2 -v 1 ) (2.27) ∆t z = ∆x 2 u 1 (x 2 -x p ) ∆z 2 w 2 w 1 ∆z + (w 2 -w 1 )(w p -w 1 ) u 1 ∆z 2∆x(w 2 -w 1 )
.

(2.28)

The 5 other similar situations (the solid voxel is either on top, below, on the left, behind or in front of the cell) can be solved by symmetry. The derivation of the 57 other cases (2 solid cells or more) are either derived or explained in the appendix, see Section 8.1.2. In summary, an analytical solution for any situation can be derived, allowing for an accurate computation of the particle's position evolution in time. The algorithm has also the advantage to be fast since, due to its analytical trajectory solution within a cell, it only needs one step to compute the exit location of the particle. This largely differs from the number of steps an Euler scheme (or even a higher order scheme such as RK4) would require for a similar accuracy. The final algorithm can be summarized as follows: This algorithm can be repeated until reaching the end of the domain or the criterion of choice for the end of the simulation. This tool allows us to compute accurately and efficiently the advective motion of any particle in the domain. An illustrative example of the extended Pollock algorithm is displayed in Figure 2.7. In the following, we detail the methodology to add a diffusive motion to these displacements in order to be able to simulate any Péclet number transport situation.

Diffusion

To simulate transport simulation under finite Péclet number conditions, one needs to insert a process to account for the diffusive motion that particles experience. In this section we discuss how to account for this diffusive displacement through a random walk methodology. CHAPTER 2. METHODOLOGY FOR FLOW AND TRANSPORT

Bases

Accounting for diffusion in conservative and reactive particle tracking transport models is often done through random walks [START_REF] Ahlstrom | Multicomponent mass transport model: theory and numerical implementation (discrete-parcel-random-walk version)[END_REF][START_REF] Ackerer | Random-walk method to simulate pollutant transport in alluvial aquifers or fractured rocks[END_REF][START_REF] Tompson | Numerical simulation of chemical migration in physically and chemically heterogeneous porous media[END_REF][START_REF] Maier | Simulation of flow through bead packs using the lattice Boltzmann method[END_REF][START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF]. The Langevin equation for the particle motion is given by [START_REF] Risken | The Fokker-Planck Equation[END_REF])

dx dt = v[x(t)] + √ 2Dξ (t) (2.29)
with D the diffusion coefficient and ξ a Gaussian white noise with zero mean and the correlation

ξ (t)ξ (t ) = δ i j δ (t -t ), (2.30)
where δ i j is the Kronecker delta. For simplicity, we consider d = 1 dimension. Let us integrate now the trajectory over the time of flight t f that is calculated advectively,

x(t + t f ) = x(t) + ∆x f + √ 2D t+t f t ξ (t )dt , (2.31)
where we denote

∆x ξ (t f ) = √ 2D t+t f t ξ (t )dt . (2.32)
It is also Gaussian distributed (because the sum of Gaussian random variables is again a Gaussian random variable, or, more generally, the sum of random variables with finite variance converge towards Gaussian random variable, according to the Central Limit Theorem).

The mean of ∆x ξ (t f ) is zero and its variance is given by

∆x ξ (t f ) 2 = 2Dt f . (2.33)
Thus, the advective-diffusive displacement during the time t f is given by

x(t + t f ) = x(t) + ∆x f + 2Dt f η(t), (2.34)
where η is a Gaussian random variable with 0 mean and unit variance.

Implementation

In order to discuss the results in function of the Péclet number, we implemented diffusion in our particle tracking code. It is implemented as follows. For every time step each particle does two jumps, first an advective step and then a diffusion step. The advective step is performed as described in Section 2.2.1.2, the particle follows its streamline during time ∆t. Knowing this ∆t, the diffusion jump in each direction can then be computed as

∆x di f f = √ 2∆tDξ 1 , (2.35
)

∆y di f f = √ 2∆tDξ 2 , (2.36) ∆z di f f = √ 2∆tDξ 3 , (2.37)
where D is the diffusion coefficient, and ξ 1 , ξ 2 and, ξ 3 are organized according to a normal distribution of mean and variance values equal to 0 and 1, respectively. These diffusion jumps can then be added to the component of the vector position x n as

x n+1 = x n + ∆x adv + ∆x di f f , (2.38) y n+1 = y n + ∆y adv + ∆y di f f , (2.39) z n+1 = z n + ∆z adv + ∆z di f f . (2.40)
This allows to add the diffusive motion into our particle tracking algorithm. Figure 2.8 displays an illustrative example of a particle performing sequentially the advective motion and the diffusive motion. In practice we use an equivalent scheme which we detail in the CHAPTER 2. METHODOLOGY FOR FLOW AND TRANSPORT appendix, see Section 8.1.3. The aim is to avoid the cost of generating normal random numbers. We also discuss there the choice of the time step that now need to be taken small enough for the scheme to be equivalent to the diffusion equation and also discuss the boundary conditions that we use at the solid-void interfaces. The Péclet number for a simulation can be changed by varying either the diffusion coefficient D or by multiplying the velocity field by a scalar.

Reaction

In this section, we aim to add an irreversible bimolecular reaction A+B → k C to our particle tracking algorithm. We first briefly discuss in Section 2.2.3.1 the methodologies that exist before introducing our reactive particle tracking in Section 2.2.3.2.

Bases

Most of the existing reactive particle tracking methods (particle tracking coupled with reaction) assume that a A particle and a B particle may react according to some give probability if the distance that separates them is smaller than some reaction radius r [START_REF] Schmidt | A kernel-based Lagrangian method for imperfectly-mixed chemical reactions[END_REF][START_REF] Edery | Particle tracking model of bimolecular reactive transport in porous media[END_REF]. There exist different manners of defining this reaction radius and this probability of reaction. [START_REF] Edery | Particle tracking model of bimolecular reactive transport in porous media[END_REF] forces particles to react whenever there are within this radius r. However, they rely on physical properties of the medium to make the reaction radius r representative of the degree of mixing between reactants. Therefore, r has to be modified in function of the size and of the properties of the medium. Other studies, based on kernel density estimators [START_REF] Schmidt | A kernel-based Lagrangian method for imperfectly-mixed chemical reactions[END_REF] use a time-varying reaction radius. However, Perez et al. (2019a) state that these models suffer from an overprediction of the reaction product at late time due to an overly large radius r(t). This is because they compensate the lack of particles in their simulations by a large reaction radius r. Other authors (see, for example, [START_REF] Benson | Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes[END_REF][START_REF] Ding | Modeling bimolecular reactions and transport in porous media via particle tracking[END_REF] have used a co-location probability and a conditional probability of reaction to control the number of reactant pairs. However, their methodology connects the reaction rate to the total initial number of particles N. This means that changing the number of particles for a given simulation implies solving a different chemical problem. Perez et al. (2019a) proposed a reactive particle tracking methodology that they prove to be equivalent to the advection diffusion reaction equation (ADRE). We use this method to insert this irreversible bimolecular reaction into our particle tracking algorithm. In next section we discuss this methodology and summarize its implementation.

TRANSPORT 39

Implementation

In our simulations the concentration of each type of agent is expressed in term of number densities. The Eulerian reactive transport formulation is described by the advection diffusion reaction equations:

∂ c A (x,t) ∂t + ∇ • v(x,t)c A (x,t) -D∇ 2 c A = -kc A (x,t)c B (x,t), (2.41) ∂ c B (x,t) ∂t + ∇ • v(x,t)c B (x,t) -D∇ 2 c B = -kc B (x,t)c A (x,t), (2.42) ∂ c C (x,t) ∂t + ∇ • v(x,t)c C (x,t) -D∇ 2 c C = kc A (x,t)c B (x,t), (2.43)
where v(x) is the velocity vector at location x, c A (x,t), c B (x,t), and c C (x,t) are the concentrations at location x and time t of species A, B, and C respectively. The reaction process that we implement is equivalent to this system of equations (see, Perez et al., 2019a). The ADRE reactive transport formulation assumes that its support scale is well-mixed and that the reaction rate k is smaller than the mass transfer rate (the Damköhler number is smaller than 1). The use of a Lagrangian particle based model implies that the volume is well-mixed since the concentration is distributed between all the particles. We chose the reaction radius to be r = √ 24D∆t since on that scale the mixing can be assumed to be uniform (Perez et al., 2019a). The reaction is then performed as follows.

At each time step, the position of the particles is updated according to the advection and the diffusion processes introduced in Sections 2.2.1.2 and 2.2.2.2 respectively. Then the distances between the A particles and the B particles are computed. We focus on the point of view of a B particle (see Figure 2.9a for an illustrative example) and denote by N A (t) the number of A particles that are located in a disk of radius r around the B particle.

Then, the B particle has an equal chance to react with any of the N A (t) particles since we do not limit the mass transfer inside the support volume. Let p r (∆t) = k∆t N 0 ∆V be the probability for a pair of particle to react during the time ∆t with N 0 the total number of particles. The probability for the B particle to survive is then p s (∆t) = 1p r (∆t). Thus, considering that the B particle can react with any of the particles N A (t) of the support volume, its probability of surviving at the time t after a time step ∆t can be expressed as:

P s (t, ∆t) = [1 -p r (∆t)] N A (t) .
(2.44)

Therefore, its reaction probability can simply be computed as t) .

P r (t, ∆t) = 1 -P s (t, ∆t) = 1 -[1 -p r (∆t)] N A (
(2.45) Then, in practice, a Bernoulli trial determines if the reaction occurs. In the eventuality it happens, the B particle and the closest A are removed from the system while a C particle is placed between the two particles, see Figure 2.9b. The C particle then starts following the flow.

Continuous Time Random Walks

The continuous time random walk framework considers particle motion as a random walk in space and time,

x n+1 = x n + ξ n , t n+1 = t n + τ n . (2.46)
It is continuous in time in the sense that the particle time is a continuous random variable as opposed to a discrete time random walk, for which time evolves in discrete steps. Transition length and time are in general coupled, but independent at subsequent steps. Their joint distribution is denoted here by ψ(x,t). The marginal distributions in space and time are denoted by Λ(x) and ψ(t), respectively. The particle distribution p(x,t) in this framework

CONTINUOUS TIME RANDOM WALKS

is given by

p(x,t) = t 0 dt R(x,t ) ∞ t-t dt ψ(t ), (2.47)
where R(x,t) denotes the frequency at which particles arrive at the position x at time t, the second term denotes the probability that a particle's transition takes longer than tt , or in other words that the particle stays at site x until time t after it arrived at time t . The R(x,t)

satisfies the Chapman-Kolmogorov type equation

R(x,t) = ρ(x)δ (t) + dx t 0 dt ψ(x -x ,t -t )R(x ,t ).
(2.48)

It can be read as follows. The probability per time that the particle arrives at (x,t) is

given by the probability that it arrives at (x ,t ) times the probability to make a transition to (x,t). The first term on the right side denotes the initial condition ρ(x) = p(x,t = 0).

Equations (2.47) and (2.48) can be combined into a single governing equation for p(x,t) using Laplace transforms. The Laplace transforms of (2.47) and (2.48) read as

p(x, λ ) = R(x, λ ) 1 -ψ(λ ) 1 -ψ(λ ) , (2.49) R(x, λ ) = ρ(x) + dx ψ(x -x , λ ) R(x, λ ), (2.50)
respectively. Solving (2.49) for R(x, λ ) and inserting the resulting expression into (2.50)

gives λ p(x, λ ) = ρ(x) + dx K (x -x , λ ) p(x , λ ) -p(x, λ ) , (2.51)
where we defined the memory kernel

K (x, λ ) = λ ψ(x, λ ) 1 -ψ(λ ) . (2.52)
Inverse Laplace transform of (2.51) gives the generalized master equation

∂ p(x,t) ∂t = dx t 0 dt K (x -x ,t -t ) p(x ,t ) -p(x,t ) . (2.53)
This equation describes the evolution of the particle distribution p(x,t) under memory.

Changes in the particle density at time t (left side) are related to particle densities at much CHAPTER 2. METHODOLOGY FOR FLOW AND TRANSPORT earlier times t . Transport is in general history dependent. If the space and time increment are independent and the time increment is distributed exponentially, this means

ψ(x,t) = Λ(x)τ -1 0 exp(-t/τ 0 ), (2.54)
the memory kernel reduces to

K (x,t) = τ -1 0 Λ(x)δ (t).
(2.55)

Thus, the generalized master equation reduces to the local in time master equation

∂ p(x,t) ∂t = dx τ -1 0 Λ(x -x p(x ,t ) -p(x,t ) .
(2.56)

The CTRW framework has been employed to model a wide variety of history-dependent processes. The underlying physical mechanisms are comprised in the joint distribution ψ(x,t) of transition lengths and times. In this thesis we extend the CTRW framework for memory-dependent hydrodynamic transport and the upscaling of anomalous dispersion from the pore to the Darcy scale.

CHAPTER 3

Stochastic Dynamics of Lagrangian Pore-Scale Velocities in Three-Dimensional Porous Media

Chapter 3

Stochastic Dynamics of Lagrangian Pore-Scale Velocities in Three-Dimensional Porous Media

Introduction

Understanding the dynamics of pore-scale flow and transport is a central issue for the modeling and upscaling of porous media phenomena and processes from the pore to the Darcy scale such as hydrodynamic dispersion, the filtration of bacteria and colloids and the mixing of dissolved chemicals and reactions between them. The sound upscaling of these processes and their modeling on the Darcy scale contribute to the understanding of the hydrodynamics of porous media and play an important role in environmental and industrial applications such as groundwater and soil remediation, the assessment of geological gas and waste storage, geothermal energy and petroleum production.

Pore-scale flow heterogeneity is the cause of hydrodynamic dispersion but also of preasymptotic non-Fickian transport. These phenomena are directly linked to Lagrangian velocity statistics. In fact, dispersion in heterogeneous flows such as turbulent flow and flow through heterogeneous porous media is quantified in terms of the covariance of Lagrangian velocities [START_REF] Taylor | Diffusion by continuous movements[END_REF][START_REF] Dagan | Theory of solute transport by groundwater[END_REF]. Pre-asymptotic behaviors such as early and late solute arrivals compared to Fickian predictions and non-linear scaling of solute dispersion [START_REF] Berkowitz | The role of probabilistic approaches to transport theory in heterogeneous media[END_REF][START_REF] Scher | The dynamical foundation of fractal stream chemistry: The origin of extremely long retention times[END_REF][START_REF] Levy | Measurement and analysis of non-Fickian dispersion in heterogeneous porous media[END_REF][START_REF] Becker | Interpreting tracer breakthrough tailing from different forced-gradient tracer experiment configurations in fractured bedrock[END_REF][START_REF] Gouze | Non-Fickian dispersion in porous media: 1. multiscale measurements using single-well injection withdrawal tracer tests[END_REF][START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF] can be traced back to particle retention in low velocity zones and fast transport in regions of high velocities, which give rise to broad distributions of solute residence times. The Fickian limit may be approached only at time values that are much larger than the largest residence times [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF]. The concept of residence or transi-tion times in the modeling of pore-scale particle motion was employed in the pioneering studies of de Josselin de [START_REF] De Josselin De Jong | Longitudinal and transverse diffusion in granular deposits[END_REF] and [START_REF] Saffman | A theory of dispersion in a porous medium[END_REF]. The models proposed by these authors are similar to time-domain and continuous time random walk approaches [START_REF] Delay | Simulating solute transport in porous or fractured formations using random walk particle tracking[END_REF][START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum transport models[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF] in that they consider particle motion through transitions over the characteristic pore lengths characterized by random time increments that depend on the distribution of pore-scale velocities. Recent experimental and numerical studies have shown that the occurrence of non-Fickian particle dispersion due to long advective residence times is directly linked to intermittency in the Lagrangian velocity time series [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF][START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF][START_REF] Carrel | Pore-scale hydrodynamics in a progressively bioclogged three-dimensional porous medium: 3-d particle tracking experiments and stochastic transport modeling[END_REF]. Thus, the understanding of these phenomena requires a sound characterization and understanding of the dynamics of Lagrangian and Eulerian pore-scale velocities, which have been the subject of a series of recent studies [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF]Siena et al., 2014;[START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF][START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]Jin et al., 2016a;[START_REF] Matyka | Power-exponential velocity distributions in disordered porous media[END_REF].

It is frequently assumed that the velocity statistics obtained in experiments and numerical simulations are stationary, implying that they do not evolve in time. The experimental particle tracking velocimetry data of [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] have shown that the distribution of initial particle velocities can in fact differ from the stationary velocity distribution depending on the injection volume and injection mode. This means, that the velocity distribution evolves in time, depending on the initial particle placement within the sample. The dependence of the Lagrangian velocity statistics on the initial particle velocity distribution was studied by Le [START_REF] Le Borgne | Characterization of the velocity field organization in heterogeneous media by conditional correlation[END_REF] for Darcy scale flow, by [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF] in a theoretical work, and analyzed by [START_REF] Hyman | Influence of injection mode on transport properties in kilometer-scale three-dimensional discrete fracture networks[END_REF] and [START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF] for particle motion in random fracture networks. Based on experimental particle tracking velocimetry data from three-dimensional bead packs [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] analyzed particle velocities sampled equidistantly along particle trajectories, which removes the intermittency observed for isochronous velocity series. These authors model the mean and displacement variance as well as the velocity increment statistics based on a Markov model for equidistant velocities. The evolution of equidistant velocity series reflects the spatial organization of a steady flow field in that they vary on the characteristic heterogeneity length scales. [START_REF] Shapiro | Stochastic analysis of solute arrival time in heterogeneous porous media[END_REF] and [START_REF] Cvetkovic | An exact solution of solute transport by one-dimensional random velocity fields[END_REF] proposed to analyze equidistant velocities as a basis to systematically quantify flow and travel time statistics in heterogeneous media (see also, Le [START_REF] Le Borgne | Characterization of the velocity field organization in heterogeneous media by conditional correlation[END_REF][START_REF] Gotovac | Flow and travel time statistics in highly heterogeneous porous media[END_REF].

The presented study is based on pore-scale velocity data obtained for flow in a threedimensional Berea sandstone sample, whose structure has been imaged by high resolution CHAPTER 3. STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES X-ray microtomography. We systematically quantify the evolution of particles moving along streamlines both in time (isochronous sampling) and distance (equidistant sampling) and study the effect of the injection conditions. We provide explicit relations between the different statistics, discuss the issues of ergodicity and stationarity of the measured velocity series in time and distance, and the impact of the finiteness of the rock sample under consideration. The insights gained from this analysis lay the basis for the stochastic description of the equidistantly sampled velocity series in terms of an ergodic Markov chain. We consider three different stochastic models and study their capability of predicting the evolution of the Lagrangian velocity statistics. The impact of diffusion on pore-scale particle motion is discussed in section 3.5 (see also, [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF][START_REF] Most | Evolution and persistence of cross-directional statistical dependence during finite-Péclet transport through a real porous medium[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF].

The paper is organized as follows. Section 3.2 presents the methodology underlying this study. It details the flow and particle transport equations, summarizes briefly the acquisition and segmentation of the rock sample, and explains the numerical solution method for the flow and particle tracking problems. Section 3.3 provides a comprehensive analysis of the statistics of the Lagrangian velocity magnitude both in time and distance. Then, we discuss the relations among them as well as their evolution toward stationarity. Section 3.4

investigates the capability of Markov models for predicting the stochastic dynamics of Lagrangian velocities. The conclusions on the stochastic description of the spatially sampled Lagrangian velocity dynamics as an ergodic Markov chain and the implications on the upscaling of pre-asymptotic hydrodynamic transport are given in Section 3.5.

Methodology

In this paper, we analyze the statistical properties of Lagrangian velocities for purely advective transport in pore-scale flows. Pore-scale flow in general is governed by the Navier-Stokes equation. For the pore-scale flow scenarios under consideration here, the Reynolds number Re = v c p /ν, with a characteristic pore velocity v c , a characteristic pore length p and the kinematic viscosity ν, is smaller than 1. Thus the pore-scale flow velocity or Eulerian velocity v(x) can be obtained by solving the Stokes equation

∇ 2 v(x) = 1 ν ∇p(x), (3.1)
where p(x) is the fluid pressure. Conservation of volume is expressed by ∇ • v(x) = 0. We specify constant pressure at the inlet and outlet boundaries and no-slip at the void-solid boundaries and the remaining domain boundaries. The porous rock sample and numerical solution of the pore-scale flow problem are described in sections 3.2.1 and 3.2.2.1 below.

The trajectory x(t, a) of a particle that is initially located at x(t = 0, a) = a is given by the advection equation

dx(t, a) dt = v[x(t, a)]. (3.2)
The Lagrangian velocity in the following is denoted by v(t, a) = v[x(t, a)] and its magnitude by

v t (t, a) = v[x(t, a)] . The initial velocity magnitude is denoted by v(t = 0, a) = v 0 (a).
The Eulerian velocity magnitude is denoted by

v e (x) = v[x]
. The distribution of initial particle positions is denoted by ρ(a). We consider here two different initial distributions at the inlet plane at x 1 = 0. The uniform distribution spreads particles uniformly in the pore space, this means

ρ(a) = I Ω 0 (a) V 0 , (3.3)
where Ω 0 denotes the domain in which particles are injected and V 0 its volume. The in- This injection condition represents a constant finite concentration pulse in the injection plane, such that the number of injected particles is proportional to the local flow velocity.

dicator function I Ω 0 (a) is equal to 1 if a ∈ Ω 0
The numerical particle tracking method is described in Section 3.2.2.2. Before, however, we discuss the methodology of streamwise velocity sampling for the statistical analysis of Lagrangian velocity magnitudes.

Rock Sample

The analysis of the velocity field was performed using a volume of 0.95 mm 3 cropped into a digital representation of a Berea sandstone (Upper Berea Sandstone unit, Ohio, USA) core sample of length 10 mm and diameter 6 mm. The Berea sandstone is a (quarried) sedimentary rock composed of well-sorted quartz grains held together by silica-rich cement displaying intermediate porosity and permeability values as well as intermediate pore-scale structural heterogeneity (tortuosity, pore-size distribution, etc.) of the pore network compared to standard reservoir rocks, while showing remarkable macroscopic homogeneity.

Because of this, it is a "rock standard" which is widely used as a proxy of mildly heterogeneous rock for experimental works by academic and petroleum industry and thus results

CHAPTER 3. STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES

can be easily compared (e.g., [START_REF] Bijeljic | Pore-scale modeling of longitudinal dispersion[END_REF][START_REF] Bijeljic | Signature of non-Fickian solute transport in complex heterogeneous porous media[END_REF][START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]. Furthermore, the characteristics of the pore size distribution compared to both the imagery technique resolution and the image size limitation for Navier-Stokes simulations make this material ideal for investigating a mildly heterogeneous natural material. The image was acquired at the BM5 beamline at the European Synchrotron Radiation Facility (Grenoble, France) using X-Ray micro-tomography. The 3-D volume was reconstructed from 3495 X-rays projections using the single distance phase retrieval algorithm [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF][START_REF] Sanchez | Threedimensional synchrotron virtual paleohistology: a new insight into the world of fossil bone microstructures[END_REF]. Since the Berea sandstone is a monocrystalline rock, we were able to relate the grey scale X-ray absorption directly to the porosity and transform the images into binary images (void and solid) using segmentation processes [START_REF] Smal | An automatic segmentation algorithm for retrieving sub-resolution porosity from x-ray tomography images[END_REF]. The details of the data processing can be found in [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]. The cubic subset of 300 3 voxels was fine-grained (each voxel was divided by 3 in each direction) giving a 900 3 voxels of volume 1.05 µm 3 . The characteristic pore length is p = 1.5 • 10 -4 m.

Numerical Simulations

Flow

Details of the flow simulation can be found in [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]. For completeness, we summarize them in this section. Generating the computational mesh that discretizes the geometry usually encounters two main problems. The first one is to create a mesh that is equivalent to the real digitized images while the second is to make it fine enough to get a high resolution of the flow field. In order to avoid the smoothing and averaging procedure that often takes place in the OpenFOAM mesh creation, we use an algorithm that generates a mesh composed of cubes that fit exactly the voxels of the digitized sample. To obtain a fine resolution, we divide every 3D cell of the mesh in 27 cubes resulting in a size of 1.05 µm 3 for each cell.

Then we compute the single phase pore-scale flow by solving the Stokes and continuity equations for constant viscosity and density. The equations are solved via a finite volume scheme using the SIMPLE algorithm of OpenFOAM. This algorithm, based on a pressurevelocity coupling, solves the Stokes equation iteratively and allows us to obtain steady state pressure and velocity fields. Convergence is reached when the difference between the current and the previous step is smaller than a criterion.

Flow is solved by imposing pressure boundary conditions at the inlet and at the outlet. No-slip conditions are implemented at the interfaces of the solid phase and at the boundaries of the domain. We also add 20 layers at the inlet and at the outlet of the domain to minimize boundary effects. Once convergence is reached, we extract the velocity field. The Right panel: 3D sample volume with 20 particle streamlines. Dark grey and light grey color denotes the pore space and the solid phase, respectively. The color scale of the streamlines denotes the velocity magnitude; from white for low velocity values to blue for high velocity values. The streamline colored from black to red is the one for which the space and time velocity series is reported in Figure 2.

In order to study particle displacements larger than the longitudinal sample size, particles are reinjected at the inlet boundary when they leave the flow domain at the outlet. The reinjection is processed as follow. When a given particle reaches the end of the domain at

x 1 = L, its velocity magnitude v L (a) = v t (t, a)| x 1 (t,a
)=L is computed. Then, the pore space ⌦ v 0 at the inlet plane where the flow velocity magnitude values are

v 0 = ku(x)| x 1 =0 k 2 [v L (a) v, v L (a) + v] is identified and the particle is reinjected randomly in ⌦ v 0 , v ⇡ v L (a)
200 . This procedure guarantees continuity of velocity and velocity statistics at reinjection and makes sure that particle velocities do not decorrelate artificially. As reported in the following the evolution of Lagrangian velocity statistics towards their respective steady state is not a ected by reinjection in the sense that there is no noticeable acceleration due to a potential artificial decorrelation.

Lagrangian and Eulerian Velocity Statistics

In this section we introduce and discuss the Eulerian and Lagrangian velocity statistics used to analyze and understand pore-scale particle motion. We define velocity statistics sampled isochronously and equidistantly along streamlines and the relations between them. We first give a brief account of the literature on Lagrangian velocities and their use. The concept of isochronously sampled particle velocities was used by [START_REF] Taylor | Diffusion by continuous movements[END_REF] to quantify di usion by continuous movements, more specifically by turbulent motion. A detailed statistical characterization of isochronous Lagrangian and Eulerian velocities was introduced by Lumley [1962]. [START_REF] Shapiro | Stochastic analysis of solute arrival time in heterogeneous porous media[END_REF] proposed and analyze the statistics of Lagrangian velocities sampled equidistantly along the mean flow direction. Le [START_REF] Le Borgne | Characterization of the velocity field organization in heterogeneous media by conditional correlation[END_REF]

Particle Tracking

The numerical solution of equation (3.2) for the particle trajectories, or, equivalently, streamlines of the pore-scale flow field, requires the interpolation of the flow velocities, which are defined at the faces of the finite volume voxels. Linear interpolation of each velocity component between opposing faces is volume preserving, this means ∇ • u(x) inside each voxel.

Linear interpolation has been used for particle tracking in Darcy scale heterogeneous flow fields on a routine basis [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF]. However, [START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF] found that the linear interpolation does not respect the no-slip boundary condition at the void-solid interface. Thus, in the void voxels in contact with the solid voxels these authors replaced the linear by a quadratic velocity interpolation, which is the implementation employed here to interpolate velocity values in the void voxels. Particle trajectories are simulated until exiting the physical domain, or reaching a given distance or a given elapsed time (Figure 3.1).

The particle tracking solver probes the velocity statistics using regular sampling in space or time along the streamline as described in section 3.3.1.

CHAPTER 3. STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES

In order to study particle displacements larger than the longitudinal sample size, particles are reinjected at the inlet boundary when they leave the flow domain at the outlet. The reinjection is processed as follows. When a given particle reaches the end of the domain at x 1 = L, its velocity magnitude v L (a) = v t (t, a)| x 1 (t,a)=L is computed. Then, the pore space Ω v 0 at the inlet plane where the flow velocity magnitude values are a) 200 . This procedure guarantees continuity of velocity and velocity statistics at reinjection and makes sure that particle velocities do not decorrelate artificially. As reported in the following the evolution of Lagrangian velocity statistics toward their respective steady state is not affected by reinjection in the sense that there is no noticeable acceleration due to a potential artificial decorrelation.

v 0 = u(x)| x 1 =0 ∈ [v L (a) -∆v, v L (a) + ∆v] is identified and the particle is reinjected randomly in Ω v 0 , ∆v ≈ v L (

Lagrangian and Eulerian Velocity Statistics

In this section we introduce and discuss the Eulerian and Lagrangian velocity statistics used to analyze and understand pore-scale particle motion. We define velocity statistics sampled isochronously and equidistantly along streamlines and the relations between them. We first give a brief account of the literature on Lagrangian velocities and their use. The concept of isochronously sampled particle velocities was used by [START_REF] Taylor | Diffusion by continuous movements[END_REF] to quantify diffusion by continuous movements, more specifically by turbulent motion. A detailed statistical characterization of isochronous Lagrangian and Eulerian velocities was introduced by [START_REF] Lumley | The mathematical nature of the problem of relating Lagrangian and Eulerian statistical functions in turbulence[END_REF]. [START_REF] Shapiro | Stochastic analysis of solute arrival time in heterogeneous porous media[END_REF] 2012) for spatiotemporally varying flow. Recently, the statistics of isochronously sampled Lagrangian velocity were analyzed for pore-scale particle motion [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF]Siena et al., 2014;[START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF][START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF], which are directly related to the particle dispersion [START_REF] Taylor | Diffusion by continuous movements[END_REF][START_REF] Kubo | Statistical Physics II, Non-Equilibrium Statistical Mechanics[END_REF]. [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF] and [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] considered particle velocities sampled equidistantly along trajectories, which reflects the spatial organization of pore-scale flow. In the following we detail different methods to sample velocity statistics and their properties.

LAGRANGIAN AND EULERIAN VELOCITY STATISTICS

Streamwise Velocity Sampling

We consider two sampling methods to characterize particle velocities along the streamlines: isochronous and equidistant. First, the t(ime)-Lagrangian velocity magnitude or speed is defined as v t (t, a) ≡ v[x(t, a)] . The velocity time series {v t (i∆t, a)} ∞ i=0 , with ∆t being a constant time increment, is obtained by isochronous sampling along a particle trajectory. [START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF] modeled the Lagrangian velocity time series as Markov processes in order to quantify particle motion in heterogeneous velocity fields. Yet, isochronous velocity series in steady heterogeneous flow fields have been shown to display intermittency [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF]. Figure 3.2 shows an isochronously sampled velocity series. It is characterized by long periods of low velocity values and short peaks of high velocity values. The origin of this intermittent behavior lies in the spatial organization of the flow.

The steady Eulerian velocity field varies on a length scale of the order of the average pore length p . Thus, significant changes of the flow velocity along a trajectory occur at times p /v. This explains the temporal persistence of low velocity magnitudes and high frequency of change of high flow velocities. In order to account for the spatial organization of the velocity field, we consider particle velocities sampled equidistantly along trajectories [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF]. The travel distance s(t) along a particle trajectory is given by ds(t, a) dt = v t (t, a).

(3.5)

Performing the variable transform t → s in (3.2) gives the following set of equations describing the particle trajectory:

dx(s, a) ds = v[x(s, a)] v[x(s, a)] , dt(s, a) ds = 1 v s (s, a) . (3.6)
The s-Lagrangian velocity magnitude is defined by v s (s

) ≡ v[x(s, a)] . The initial speed is denoted by v s (s = 0, a) = v 0 (a).
The velocity series {v s (i∆s, a)}, with ∆s being a constant space increment, is obtained by equidistant sampling along particle trajectories. Unlike for isochronous sampling, here velocities are sampled independently of their magnitude since this velocity value does not impact the sampling distance. Note that the system of equations (3.6) describes particle motion as a process in which the particle position is incremented by a constant value and the particle time by a variable transition time. In this sense it describes a time-domain random walk [START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF]. The s-Lagrangian velocity series shown in Figure 3.2 does not display intermittent patterns. The signal seems stationary and is characterized by a characteristic correlation scale v . To determine this correlation distance we consider the velocity covariance function for an injection into the flux, which is defined as

C v (s) = 1 L L 0 ρ(a) v s (s + s, a) -µ v s (s , a) -µ ds da, (3.7)
where ρ(a) is given by (3.4) and µ is the mean velocity

µ = 1 L L 0 ρ(a)v s (s , a)ds da. (3.8)
The velocity autocorrelation function

A v (s) = C v (s)/C v (0) is shown in Figure 3.3. The correlation length is defined as v = ∞ 0 A v (s)ds. (3.9)
We find that the velocity correlation length v equals the average pore length, p = 1.5 • 10 -4 m. Note that the sample size is about (6 p ) 3 and that the average streamline length is around 10. can be characterized by isochronous sampling along a single streamline labeled by a. This sampling mode defines the streamwise t-Lagrangian velocity PDF P(v, T, a)

P(v, T, a) = 1 T T 0 δ [v -v t (t, a)]dt, (3.10)
which in general depends on the sampling time T . In the following, the statistics obtained by isochronous sampling are marked by a hat. The statistics of the velocity series v s (s, a)

illustrated in the bottom panel of Figure 3.2 is characterized by equidistant sampling,

P(v, L, a) = 1 L L 0 δ [v -v s (s, a)]ds, (3.11)
where L is the sampling length. The PDF P(v, L, a) is referred to in the following as streamwise s-Lagrangian velocity PDF. The relation between the streamwise s-and t-Lagrangian velocity PDFs defined in (3.11) and (3.10) is obtained by the variable change s → t according to the map (3.5), which gives

P(v, L, a) = vT (L) L P[v, T (L), a], (3.12)
where T (L) is the time that the particle needs to travel the distance L along the streamline and Therefore, under ergodic conditions, the velocity statistics sampled between an ensemble of particles and along a single streamline are equivalent. Ergodicity can only be achieved if first, the sampling distance or sampling time along a streamline is large enough for the particle to experience the full velocity spectrum, and second, if the ensemble of particles is large enough to contain the full velocity statistics. The stationary s-and t-Lagrangian ensemble statistics are defined by (3.15) respectively. In practice, the initial volume V 0 is of course finite. In order to achieve ergodicity it needs to be chosen large enough to contain the significant velocity statistics;

T (L) = L 0 ds v s (s) . ( 3
P(v) = lim V 0 →∞ 1 V 0 Ω 0 v 0 (a) v 0 (a) δ [v -v s (s, a)]da, (3.14) P(v) = lim V 0 →∞ 1 V 0 Ω 0 δ [v -v t (t, a)]da,

LAGRANGIAN AND EULERIAN VELOCITY STATISTICS

57 see also Appendix 3.6.1. Thus ergodicity can be expressed as

P(v) = lim L→∞ P(v, L, a) ≡ P(v), P(v) = lim T →∞ P(v, T, a) ≡ P(v). (3.16)
Under ergodic conditions, the flux-weighting relation (3.12) implies for the stationary ensemble statistics

P(v) = v v e P(v), (3.17)
where v e is the mean Eulerian velocity magnitude. We now consider the relation between the Lagrangian PDFs and the Eulerian velocity PDF sampled over an infinite domain

P e (v) = δ [v -v e (x)]dx. (3.18)
The stationary Lagrangian PDF P(v) = P e (v) because of volume conservation as discussed in section 3.6.3.2. Thus, expression (3.17) implies for the stationary s-Lagrangian PDF

P(v) = v v e P e (v). (3.19)
This is a key relationship for the prediction of particle velocity statistics because the Eulerian velocity PDF can be determined independently from transport.

A detailed discussion on the relations between s-Lagrangian and t-Lagrangian statistics for finite sampling domains is given in Appendix (3.6.1). There it is shown that the Lagrangian statistics for the rock sample under consideration are stationary and ergodic.

Evolution of the Lagrangian Velocity Statistics and Stationarity

In the previous sections, we have seen that the s-Lagrangian and t-Lagrangian velocity statistics evolve asymptotically to different steady state distributions, which are related through flux weighting according to (3.58). In this section, we study in detail the evolution of the respective statistics from uniform and flux-weighted initial conditions. Note that since, the Lagrangian and Eulerian quantities are related through equations (3.48) and

(3.56), studying the Lagrangian statistics evolution includes studying the Eulerian statistics evolution. In the following we only refer to Lagrangian distributions.

Evolution of the s-Lagrangian Velocity Statistics

The s-Lagrangian velocity distribution for an arbitrary initial particle distribution ρ(a) is defined by The results were computed using 10 7 particles.

p(v, s) = δ [v -v s (s,

Evolution of the Lagrangian Velocity Statistics and Stationarity

In the previous sections, we have seen that the s-Lagrangian and t-Lagrangian velocity statistics evolve asymptotically to di erent steady state distributions, which are related through flux weighting according to (30). In this section, we study in detail the evolution of the respective statistics from uniform and flux-weighted initial conditions. Note that since, the Lagrangian and Eulerian quantities are related through Equations ( 20) and ( 28), studying the Lagrangian statistics evolution includes studying the Eulerian statistics evolution. In the following we only refer to Lagrangian distributions.

Evolution of the s-Lagrangian velocity statistics

The s-Lagrangian velocity distribution for an arbitrary initial particle distribution ⇢(a)

is defined by

p(v, s) = Z [v v s (s, a)] ⇢(a)da. ( 32
)
The initial velocity distribution is p 0 (v) = p(v, t = 0). Clearly, fully stationary conditions such that the s-Lagrangian PDF p(v, s) = P(v) at all s, are achieved for the initial condition

p 0 (v) = P(v).
For the flux-weighted initial distribution (4), we recover expression (15).

Figure 8 shows the evolution of p(v, s) for the flux-weighted and uniform initial particle distributions (4) and (3). The PDF evolves from both initial distributions towards its steady state P(v). For the flux-weighted initial particle distribution, p 0 (v) is skewed towards high velocity values compared to p 0 (v) for the uniform injection with a high probability weight at low velocities. For both initial distributions, the steady state P(v) is reached after a distance of s ⇡ 7`c. We note that the high velocity part of the PDF converges faster to the steady state than the low velocity part.

The initial velocity PDF corresponding to the flux-weighted injection is closer to the steady state PDF P(v) than the one corresponding to the uniform injection. This can be understood as follows. First, note that the s-Lagrangian PDF P(v, s) is related to the s-Eulerian

PDF P e (v, s) through flux-weighting, see Eq. ( 20). Let us now assume that the s-Eulerian The initial velocity distribution is p 0 (v) = p(v,t = 0). We consider the uniform and fluxweighted initial particle distributions (3.3) and (3.4). For an ergodic injection domain Ω 0 , p 0 (v) = P e (v) for the uniform injection and p 0 (v) = P(v) is equal to the stationary s-Lagrangian PDF for the flux-weighted injection. While the injection domain here is not large enough to be ergodic, the initial distribution under flux-weighted conditions is close to the stationary s-Lagrangian PDF as shown in Figure 3.5. Figure 3.5 shows the evolution of p(v, s) for the uniform and flux-weighted initial particle distributions (3.3) and (3.4).

velocity PDF P e (v, s) = P e (v) is stationary. This implies that 1 Z [v v e (x)]dx = 1 Z [v v e (x)]dx, (33)
The PDF evolves from both initial distributions toward its steady state P(v). For the fluxweighted initial particle distribution, p 0 (v) is skewed toward high velocity values compared to p 0 (v) for the uniform injection with a high probability weight at low velocities. For both initial distributions, the steady state P(v) is reached after a distance of s ≈ 7 p . We note that the high velocity part of the PDF converges faster to the steady state than the low velocity part.

Evolution of the t-Lagrangian Velocity Statistics

The t-Lagrangian velocity PDF for an arbitrary initial particle distribution is defined by

p(v,t) = δ [v -v t (t, a)]ρ(a)da. (3.21)
The initial velocity distribution is p0 (v) = p(v,t = 0), which is identical to the initial s-Lagrangian velocity PDF p 0 (v). As in the previous section, we consider the uniform and which means that the statistics contained in the domain ⌦(s) are invariant. This condition can be achieved in an ergodic medium if the injection domain is large enough to contain a representative part of the heterogeneity spectrum. Then, also the s-Lagrangian P(v, s) = P(v) is stationary and given by

P(v) = 1 V 0 Z ⌦ 0 v 0 (a) hv 0 (a)i [v v s (s, a)]da. ( 34
)
Thus, if the injection domain ⌦ 0 is ergodic, the flux-weighted injection protocol (4) represents a steady initial condition for p(v, s). For the sample under consideration, flux-weighting the velocities in the injection domain gives an initial velocity PDF p 0 (v), which is (only) relatively close to the steady state PDF P(v) because the injection domain here is not large enough to be ergodic.

Evolution of the t-Lagrangian velocity statistics

The t-Lagrangian velocity PDF for an arbitrary initial particle distribution is defined by

p(v, t) = Z [v v t (t, a)]⇢(a)da. ( 35
)
The initial velocity distribution is p0 (v) = p(v, t = 0), which is equal to the initial s-Lagrangian velocity PDF p 0 (v). For ⇢(a) given by (3), p(v, t) is equal to P(v, t) defined in ( 23). The t-Lagrangian p(v, t) is related to the s-Lagrangian statistics as

p(v, t) = 1 v 1 Z 0 p(v, t, s)ds, (36)
where p(v, t, s) is the joint PDF of velocity v s (s, a) and particle time t(s, a),

p(v, t, s) = Z [v v s (s, a)] [t t(s, a)]⇢(a)da, (37)
see Appendix A.3. As outlined in Section 3.4, under stationary conditions the two statistics are related by flux-weighting, see Eq. ( 30).

Figure 9 shows the evolution of p(v, t) for the uniform and flux-weighted initial con- As expected, p(v,t) evolves toward the steady state distribution P(v) from both initial distributions. The time for convergence toward the steady state is t > 10 4 τ c . Since the average time needed to reach the outlet of the sample is on the order of 6τ c we use the reinjection procedure to keep all particles in the domain. As for the s-Lagrangian statistics, also here, the high velocity part of p(v,t) converges faster than the low velocity part.

Synthesis

In summary, we distinguish between s-Lagrangian statistics, which are sampled equidistantly along particle trajectories, and t-Lagrangian statistics, which are sampled isochronously along particle trajectories. Moreover, we distinguish velocity PDFs that are sampled along single streamlines and velocity PDFs that are obtained by sampling from an ensemble of particles, as well as mixed sampling between particles and along streamlines. We find that the streamwise and ensemble sampled statistics eventually converge after a given streamwise travel distance or streamwise travel time. The convergence of streamwise and ensemble statistics to the same steady state distributions indicates that the underlying velocity Stationary conditions are achieved for the t-Lagrangian velocity statistics in case of a uniform injection into an ergodic subdomain Ω 0 . For the s-Lagrangian statistics this corresponds to a flux-weighted injection. In the case of ergodic conditions, the s-and t-Lagrangian steady state statistics P(v) and P(v) can be obtained by volumetric sampling over an ergodic subdomain because of their relations to the Eulerian velocity PDFs [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]. In the following, we model the s-Lagrangian velocity series as a stationary and ergodic Markov process in order to capture the evolution of the s-Lagrangian velocity statistics and its dependence on the initial conditions.

Markov model

We model the s-Lagrangian velocity series v s (s) as a stationary and ergodic Markov process. This means that the process v s (s) is fully characterized by the velocity transition probability r(v, ss |v ), which denotes the PDF of v s (s) given that v s (s ) = v . Both r(v, s|v ) and p(v, s) satisfy the Chapman-Kolmogorov equation

p(v, s) = ∞ 0 r(v, s -s |v )p(v , s )dv . (3.22)
The steady state distribution is an Eigenfunction of r(v, s|v ),

P(v) = ∞ 0 r(v, s|v )P(v )dv . (3.23)
Furthermore, the transition probability converges to the steady state distribution in the limit of s p , lim s→∞ r(v, s|v ) = P(v).

(3.24)

This implies that lim s→∞ p(v, s) = P(v), independent of the initial condition p 0 (v). Note also that the joint PDF P(v, ss , v ) of v and v under stationary conditions is given by

p(v, s -s , v ) = r(v, s -s |v )P(v ).
(3.25)

Note that [START_REF] Gotovac | Flow and travel time statistics in highly heterogeneous porous media[END_REF] studied the statistical properties of the inverse Lagrangian velocity 1/v s (s), which is termed slowness. Evolution equations for the PDF of slowness can be deduced from the Markov model for v s (s) by variable transformation. In the following, we first construct the transition probability empirically from the direct numerical simulations. This is then used to propagate the s-Lagrangian velocity statistics from uniform initial conditions. Second, we use a Markov model based on a Bernoulli process for the persistence of velocities. Third, we employ an Ornstein-Uhlenbeck (OU) process for v s (s). The simulated results of these Markov models are then compared to the simulation data presented in the previous sections.

MARKOV MODEL

Empirical Transition Probability

In this section, we use empirically computed transition probabilities for the modeling of velocity series in space. This approach is conceptually similar to the work of [START_REF] Benke | Modeling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF] and [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum transport models[END_REF] for fractured rock. In order to determine the velocity transition probability r(v, ss |v ), we discretize the velocity interval

[v , v u ]
between minimum and maximum velocities v and v u sampled in the domain into n bins of 62 CHAPTER 3. STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES

width ∆v i = v i+1 -v i , where v = v 1 and v u = v n such that v j = v 1 + j-1 ∑ i=1 ∆v j . (3.26)
The empirical transition probability is given by mixed streamwise and ensemble sampling as

T i j (∆s) = 1 V 0 Ω 0 1 L L 0 I(v j ≤ v(s + ∆s, a) < v j + ∆v j )| v i ≤v(s ,a)<v i +∆v i ρ(a)
T i j (∆s) = v j +∆v j v j v i +∆v i v i p(v, ∆s, v )dvdv v i +∆v i v i P(v )dv . (3.28)
Note that the empirical determination of the transition probability requires stationary initial conditions. As the initial velocity distribution p 0 (v) here is not stationary, as discussed in the previous section, sampling should start once stationary conditions are achieved at approximately s = 7 p , see Figure 3.12. In practice, mixed sampling along the trajectories of lengths s p guarantees stationary conditions. Figure 3.7 shows the transition matrix T i j (∆s) computed using 10 6 particle trajectories for different lags ∆s, and n = 100 logarithmically spaced velocity bins such that v i+1 = v i exp(1/n). The smaller the ∆s the higher the correlation and the more diagonal is the matrix. In principle, any ∆s would provide a good estimate for the transition matrix T i j (∆s) if the velocity correlation were exponential because in this case, the slope of the correlation function would equal its value everywhere and it could be uniquely characterized by the correlation length p . This is not the case here. Figure 3.3 shows that the correlation function drops for small distances ∆s p faster than for larger distances. This means that for small ∆s an exponential fit simulates a shorter correlation length that the actual full correlation function. This implies that a transition matrix determined at short ∆s sees only this sharp drop and thus underestimates the true correlation. The lag distance ∆s needs to be large enough such that the correlation information can be sampled. This means here ∆s ≥ p . For estimating the evolution of the s-Lagrangian statistics we choose ∆s = 8 p /3. velocity statistics p(v, s), which reaches the steady state distribution after the same distance as the data form the direct numerical simulations. Thus, the evolution of the s-Lagrangian velocities can be well represented as a Markov process. In the following section, we consider a Bernoulli process as a Markov relaxation model for this evolution.

Bernoulli Process

We (3.30)

In the limit ∆s → 0, we obtain the evolution equation [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF])

∂ p(v, s) ∂ s = - 1 c [p(v, s) -P(v)] , (3.31)
whose solution for the initial condition p 0 (v) is

p(v, s) = P(v) + exp(-s/ c ) [p 0 (v) -P(v)] . (3.32)
Thus, we obtain for the mean velocity v s (s), the explicit analytical expression

v s (s) = ( v 0 -v s ) exp(-s/ c ) + v s , (3.33)
where v 0 is the average initial velocity and v s the average s-Lagrangian steady state velocity. We use this expression to estimate the characteristic length scale from the direct numerical simulations as shown in Figure 3.9. We find c = 2.5 p , which is of the order of the pore length. The Bernoulli model converges to the steady state distribution after s ≈ 3 c , but does not reproduce the velocity PDFs at intermediate distances. Note that the Bernoulli model uses the same convergence rate -1 c for all velocities. While this model represents the evolution of the high velocity part of p(v, s) is relatively well, the low velocity part evolves faster than the data from the direct numerical simulations. This indicates that the convergence rate may be velocity dependent. In the next section, we model this behavior with an alternative Markov model for v s (s).

Ornstein-Uhlenbeck Process

In this section, we consider a velocity Markov model for the evolution of v s (s) that is based on the Ornstein-Uhlenbeck (OU) process [START_REF] Gardiner | Stochastic Methods[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF], where ξ (s) is a Gaussian white noise characterized by zero mean, ξ (t) = 0 and covariance ξ (s)ξ (s ) = δ (ss ). The angular brackets denote the noise average over all realizations of ξ (s). The Ornstein-Uhlenbeck process has originally been considered for the modeling of the stochastic (in time) velocity fluctuations of Brownian particles [START_REF] Langevin | Sur la théorie du mouvement brownien[END_REF][START_REF] Uhlenbeck | On the theory of the Brownian motion[END_REF][START_REF] Risken | The Fokker-Planck Equation[END_REF], later also for particle velocities in turbulent flows [START_REF] Pope | Turbulent Flows[END_REF]. From a mathematical point of view, the OU model is a stationary Gaussian Markov process. Its increments are Gaussian random variables. Its distribution φ (w, s) satisfies the Fokker-Planck equation [START_REF] Risken | The Fokker-Planck Equation[END_REF])

dw(s) ds = --1 c w(s) + 2 -1 c ξ (s), (
∂ φ (w, s) ∂ s --1 c ∂ wφ (w, s) ∂ w --1 c ∂ 2 φ (w, s) ∂ w 2 = 0. (3.35)
It relaxes from any initial distribution φ 0 (w) to a Gaussian steady state distribution φ (w), which has zero mean and unit variance for the specific process (3.34). We use this process here to model the stochastic evolution of the particle velocity v s (s) and the relaxation of its statistics from any initial distribution p 0 (v) toward the steady state P(v). This requires to map v s (s) onto w(s) through their steady state PDFs. This is done through the Smirnov

MARKOV MODEL

67 transform [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF]

w(s) = Φ -1 (Π[v s (s)]) ≡ M [v s (s)], v s (s) = Π -1 (Φ[w(s)]), (3.36)
where Π(v) and Φ(w) are the cumulative distributions of v(s) and w(s),

Π(v) = v 0 P(v )dv , Φ(w) = w -∞ φ (w )dw . (3.37)
The latter is given by Φ(w) = [1 + erf(w/ √ 2)]/2. As φ (s) is a unit Gaussian, this map generates w(s) as the normal score of v s (s). This map guarantees that p(v, s) evolves from any initial distribution p 0 (v) toward its steady state P(v) on the relaxation scale c , which is set equal to the one determined for the Bernoulli model in the previous section. Note that all the normal scores w(s) evolve with the same rate -1 c . Note that the transition probability r w (w, ∆s|w ) is given by the Gaussian distribution [START_REF] Risken | The Fokker-Planck Equation[END_REF]

r w (w, ∆s|w ) = exp -[w-w exp(-∆s/ c )] 2 2[1-exp(-2∆s/ c )] 2π [1 -exp(-2∆s/ c )] . (3.38)
The transition probability for the velocity process r(v, ∆s|v ) is given in terms of r w (w, ∆s|w )

according to the map (3.36)

r(v, ∆s|v ) = r w [M (v), ∆s|M (v )] dM (v) dv , (3.39)
which in general leads to velocity-dependent convergence rates for p(v, s).

The process (3.34) is solved for an ensemble of particles. The initial values w(s = 0) = w 0 are obtained from v s (s = 0) = v 0 by the map (3.36) as w 0 = Φ -1 [Π(v 0 )], where the v 0 are distributed according to p 0 (v). Once w 0 is obtained, the process (3.34) is solved numerically using an explicit Euler scheme,

w n+1 = w n --1 c w n ∆s + 2 -1 c ∆sξ n , (3.40)
where w n = w(n∆s) and ξ n is a Gaussian random variable with 0 mean and unit variance.

The value w n+1 is transformed back to the velocity v s (n∆s) via Eq. (3.36) at every step.

The discretization of scheme (3.40) is chosen such that ∆s ≤ c /10.

Synthesis

In summary, based on the stationary and ergodic properties of the s-Lagrangian velocity series, we model their stochastic dynamics as an ergodic Markov chain. We consider three different Markov models. First, a Markov model based on an empirical transition probability, which is obtained from conditional equidistant velocity sampling along streamlines.

This model naturally reproduces the evolution of the s-Lagrangian velocity PDF and confirms the Markovian nature of the velocity transitions. Second, we consider a Bernoulli velocity model, which at each step either persists at the velocity of the previous step or changes to a new velocity, which is randomly sampled from the stationary s-Lagrangian PDF. This model yields an evolution of the s-Lagrangian PDF from an initial to the stationary PDF. However, it uses the same convergence rate for all velocity classes, which does not capture the evolution at small velocities. Third, we consider a velocity transition model that is based on an OU process for the normal scores of the s-Lagrangian velocities.

This process correctly predicts the full evolution of the s-Lagrangian velocity PDF and is parameterized by the stationary Lagrangian velocity PDF and a characteristic relaxation scale c . The former is related to the Eulerian velocity PDF, a flow attribute; the latter is of the order of the characteristic pore-length. Thus, this stochastic velocity model can be parameterized in terms of hydraulic and geometric characteristics of the porous medium.

Conclusions

We have presented a comprehensive analysis of Lagrangian pore-scale velocity series. Even though the study is based on velocity data in the three-dimensional pore-structure obtained from X-ray microtomography of a Berea sandstone sample, the presented methods and results are valid for particle motion in steady pore-scale flows in general. Our analysis has revealed the stochastic dynamics of particle velocities and led to the formulation of a predictive modeling approach for the velocity evolution based on Markov processes for the streamwise Lagrangian velocities. These results are part of the endeavor of setting up an upscaling framework for hydrodynamic flow and transport from the pore to the Darcy scale. The past years have seen a significant increase of experimental and numerical porescale studies along with improved imaging techniques and computational resources. The presented methods for the statistical analysis of pore-scale velocity data provide new tools for the interpretation of such experimental and numerical data and their use in the upscaling of flow and transport.

The evolution of the velocity statistics represents a key feature that needs to be accounted for both in the interpretation of experimental and numerical velocity data and in 3.5. CONCLUSIONS 69 the modeling and upscaling of particle transport. For example in particle tracking and particle imaging velocimetry, the measured velocity distributions may be dependent on the initial preparation, this means on the seeding of the injection volume with particles, and not be representative of the porous sample. Furthermore, data analysis often invokes stationarity of the measured particle velocities, which in general, however, is not the case and depends again on the injection condition. Pore-scale velocity variability is at the origin of hydrodynamic dispersion and other transport phenomena observed on the Darcy scale.

The Kubo formula provides a measure for hydrodynamic dispersion in terms of the time integral of the t-Lagrangian velocity covariance,

D(t) = t 0 v t (t )v t (t) dt , (3.41)
where v t (t) denotes the fluctuation of the t-Lagrangian velocity around its mean. Dispersion in general evolves in time and depends on the initial conditions (stationary or nonstationary) and the time evolution of the velocity statistics. For example, dispersion at times smaller than the advection time τ v is ballistic and given by D(t) = σ 2 0 t where σ 2 0 is the variance of the PDF of initial velocities p 0 (v), which clearly depends on the initial velocity distribution. The asymptotic behavior is determined by the velocity correlation time and velocity variance, which are related to the intermittent temporal velocity signals. The tailing of particle breakthrough curves is determined by the occurrence of low velocities and their spatial persistence. Thus, the retention phenomena also depends on the evolution of the velocity statistics and the initial preparation of the system. For example, a uniform initial particle distribution emphasizes more the low end of the velocity spectrum than a flux-weighted. Thus the corresponding breakthrough curves, or residence time distributions in a sample may be significantly different depending on the initial conditions.

These behaviors as they evolve in time hold a certain complexity, which is reflected in the intermittent features of Lagrangian velocity time series. This complexity can be removed by applying a different sampling protocol, namely by sampling equidistantly along particle trajectories. This streamwise spatial point of view provides a significant simplification of otherwise complex phenomena and thus opens new possibilities for transport modeling and upscaling. The formulation of the s-Lagrangian velocity magnitude as an ergodic Markov chain renders particle motion naturally as a (correlated) continuous time random walk [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]Le Borgne et al., 2008;[START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF] or timedomain random walk [START_REF] Benke | Modeling conservative tracer transport in fracture networks with a hybrid approach based on the Boltzmann transport equation[END_REF][START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum transport models[END_REF], because streamwise particle motion can be modeled in terms of fixed spatial steps ∆s, which take tortuosity [START_REF] Koponen | Tortuous flow in porous media[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF].

The OU process for the normal scores of the s-Lagrangian velocities correctly predicts the full evolution of the s-Lagrangian velocity PDF and is parameterized by the stationary Lagrangian velocity PDF and the characteristic correlation scale. The former is related to the Eulerian velocity PDF, a flow attribute, the latter is of the order of the characteristic pore-length. This stochastic velocity model can be parameterized in terms of hydraulic and geometric characteristics of the porous medium. Thus, it is a predictive model in the sense that it can be based on the characterization of transport independent quantities, which is an important step for flow and transport upscaling from the pore to the Darcy scale. While significant progress has been made [START_REF] Bibliography De Anna | Prediction of velocity distribution from pore structure in simple porous media[END_REF][START_REF] Alim | Local pore size correlations determine flow distributions in porous media[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF], the relation between pore structure and pore velocity distribution still remains an open issue. In addition to the OU velocity model, we consider a Bernoulli process, which reproduces the velocity evolution qualitatively, but lacks the correct convergence rates for low and intermediate velocities. Yet, due to its simplicity, it may serve to obtain fast qualitative estimates of transport features related to the evolution of particle velocity statistics. Finally, note that the methodology used here applies to transport in steady flow through heterogeneous media in general, for which a relaxation of the Lagrangian velocity statistics toward a steady state can be observed such as Darcy scale fractured and porous media [START_REF] Cvetkovic | Analysis of nonlinear effects on tracer migration in heterogeneous aquifers using Lagrangian travel time statistics[END_REF][START_REF] Le Borgne | Characterization of the velocity field organization in heterogeneous media by conditional correlation[END_REF][START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF][START_REF] Kang | Anomalous transport in disordered fracture networks: Spatial Markov model for dispersion with variable injection modes[END_REF].

We consider here purely advective particle motion and do not account for the effect of diffusion on particle motion. Thus, the derived stochastic framework is directly relevant for advection-dominated pore-scale transport. In fact, practically relevant pore-scale Péclet numbers may range from 10 -2 to 10 6 [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF]). In the presented Markov models, velocity transitions occur essentially with a fixed spatial 3.6. APPENDIX 71 frequency which is given by the inverse velocity correlation length. In the presence of diffusion, velocity transitions can also occur due to particle transitions between streamlines, a process that is related to a constant frequency in time, namely the inverse diffusion time over the characteristic velocity length scale. Advective and diffusive velocity transitions may depend on a local Péclet number. Thus, the derived Markov model provides a basis to account for the impact of velocity variability and diffusion on hydrodynamic dispersion.

Furthermore, pore-scale flow variability has an impact on processes such as the filtration of colloidal particles and bacteria [START_REF] Liang | Evaluation of bacterial run and tumble motility parameters through trajectory analysis[END_REF] as well as mixing between dissolved chemicals [START_REF] Kree | Scalar mixtures in porous media[END_REF], while these processes are also affected by other factors such as volume exclusion and interactions with the solid matrix as well as diffusion, for example, the derived stochastic model for Lagrangian particle velocities may serve as a starting point to account systematically for the effect of hydrodynamic variability.

In summary, the fact that the Lagrangian velocity statistics are stationary allows for the stochastic description of the s-Lagrangian velocity dynamics as an ergodic Markov chain. This stochastic framework renders particle motion as a correlated continuous time random walk. The consequences of the stochastic s-Lagrangian velocity dynamics for the prediction of pre-asymptotic spatial and temporal transport characteristics and their systematic upscaling are studied elsewhere.

Appendix

Equidistant and Isochronous Ensemble Statistics

In this section, we discuss the relations between the s-and t-Lagrangian ensemble statistics for finite sampling domains.

Equidistant Sampling: s-Lagrangian Statistics

The ensemble s-Lagrangian velocity PDF is obtained by sampling the velocity magnitude v s (s, a) at a given streamline distance s in the flux-weighted ensemble of particles comprised in the injection domain

Ω 0 P(v, s) = 1 V 0 Ω 0 v 0 (a) v 0 (a) δ [v -v s (s, a)]da, (3.43)
each particle being weighted by its initial velocity. This PDF can be computed for any distance s ≥ 0 and in general evolves with distance. The mixed s-Lagrangian PDF is defined by sampling along streamlines and between particles as

P m (v, L) = 1 L L 0 1 V 0 Ω 0 v(a) v(a) δ [v -v s (s, a)]dads.
(3.44)

CHAPTER 3. STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES

This method samples more statistics than its ensemble and streamwise counterparts since it integrates over all particles labeled by a and distances s traveled. Under ergodic conditions, sampling along streamlines and ensemble sampling are equivalent. As discussed in the main text, this can be achieved for sampling distances L and initial volumes V 0 large enough that a representative part of the velocity variability can be experienced. Under these conditions, the streamwise, ensemble and mixed s-Lagrangian PDFs are identical and independent of L and a (streamwise) and s (ensemble),

P(v) = P(v) = P m (v).
(3.45)

Figure 3.12 shows the streamwise, ensemble and mixed s-Lagrangian PDFs for L ≈ 10 8 p in the streamwise case and s = 7 p for the ensemble. The three different PDFs are in very good agreement, which implies that the reinjection method detailed in Section 3.2.2.2 is ergodic (convergence of the streamwise PDF), and that stationary conditions are already attained within the sample size (convergence of the ensemble PDF). The evolution of the ensemble s-Lagrangian PDF is analyzed in detail in Section 3.3.2. As the different statistics are identical in steady state, we refer to the steady state distribution as P(v).

Equidistant Sampling: s-Eulerian Statistics

The fluid flow map ϕ s , defined as ϕ s : a → x(s, a) (3.46) maps the initial particle position a on the particle position x(s, a) at distance s according to (3.6). The s-Eulerian velocity PDF is obtained by volumetric sampling of the Eulerian velocity magnitude v e (x) in the subdomain Ω(s) = ϕ s (Ω 0 ) ⊂ Ω f , where Ω f is the flow domain,

P e (v, s) = 1 V (s) Ω(s) δ [v -v e (x)]dx. (3.47)
The relation between the s-Eulerian velocity PDF P e (v, s) and the ensemble s-Lagrangian velocity PDF P(v, s) is obtained by using the map (3.46) in order to transform the integration variable a → x in (3.43). Note that the map (3.46) is not volume preserving. Thus V (s) = V (t). In section 3.6.3.1 we show that

P(v, s) = v µ e (
P(v) = v v e P e (v).
(3.50)

10 -2 10 -1 10 -4 10 -3 10 -2 10 -1 10 0 10 1 pdf v/ v Figure 3
.12: Comparison between the streamwise P(v) (orange dashed line), ensemble P(v) (blue circles, computed with 10 7 particles), and mixed P m (v) (solid black line, computed with 10 5 particles) of the s-Lagrangian velocity PDFs. PDF = probability density function.

Isochronous Sampling: t-Lagrangian Statistics

Sampling of the velocity magnitude v t (t, a) between particles gives the ensemble t-Lagrangian

PDF P(v,t) = 1 V 0 Ω 0 δ [v -v t (t, a)]da, (3.51)
where Ω 0 is the fluid domain in which particles are initially placed, and V 0 is its volume. The mixed t-Lagrangian PDF samples velocity magnitudes both between particles and along particle trajectories

Pm (v, T ) = 1 V 0 Ω 0 1 T T 0 δ [v -v t (t, a)]dtda. (3.52)
The mixed method is often used for the empirical determination of velocity statistics from particle tracking velocimetry (PTV) because it yields better statistics than either sampling along a single trajectory or between particles. Under ergodic conditions, the sampling of the velocity magnitude along a streamline for a long enough time T is identical to ensemble sampling for a large enough time t. The sampling time T and initial domain Ω 0 need to be large enough such that the sampled velocity variability is representative. Under these conditions, P(v, T, a) = P(v) is independent of a and T and P(v,t) = P(v) is independent of the sampling time t such that

P(v) = P(v) = Pm (v).
(3.53) Figure 3.13, shows P(v, T, a), P(v,t,V 0 ) and Pm (v, T,V 0 ). The three statistics are in good agreement, which indicates that in the numerical simulations, ergodic conditions are reached. (3.57)

(t) = ϕ t (Ω 0 ), Pe (v,t) = 1 V 0 Ω(t) δ [v -v e (x
In the following, we refer to the steady state distribution as P(v) because the three statistics are identical in the steady state.

Relations Between Isochronous and Equidistant Statistics

We have seen in (3.57) that the statistics of the isochronously sampled Lagrangian velocity P(v) along a trajectory equals the Eulerian velocity statistics Pe (v) under ergodic conditions and that the streamwise s-and t-Lagrangian velocity statistics are related through a flux weighted relation (3.17). Furthermore, under ergodic conditions, we know from (3.45) that the ensemble equals the streamwise s-Lagrangian PDF, P(v) = P(v), and from (3.53) that the ensemble is identical to the streamwise t-Lagrangian PDF, P(v) = P(v). Thus, we obtain from (3.17) the following relation 3.6. APPENDIX 77

P(v) = v v e P(v). ( 3

Relations Between s-and t-Lagrangian Statistics

In the following, we provide some details on the derivation of the relations between s-and t-Lagrangian statistics reported on in the previous section.

3.6.3.1 Relation Between the s-Eulerian Velocity PDF P e (v, s) and the Ensemble s-Lagrangian Velocity PDF P(v, s)

We derive the relation between the s-Eulerian velocity PDF P e (v) and the ensemble s-Lagrangian velocity PDF P(v, s). P(v, s) is defined as:

P(v, s) = 1 V 0 Ω 0 v(a) v(a) δ [v -v s (s, a)]da, (3.60)
which is originally defined on the initial injection domain Ω 0 and is parameterized through flux weighting. V 0 is the volume of Ω 0 , and a is the initial position of particle a. We can derive this definition to any domain Ω(s), which is the domain occupied by the particles after they have traveled the distance s along their streamline, this is achieved using the map a → x(s, a). Doing so we obtain:

P(v, s) = 1 V 0 Ω(s) J(s, a) -1 v(a) v(a) δ [v -v e (x(s, a))]dx, (3.61)
where J(s, a) = dx(s, a)/da is the Jacobian of the transformation. It can be determined as follows. First, we note that its derivative is given by (Batchelor, 2000, p., 75)

d ds J(s, a) = J(s, a)∇ • v[x(s, a)] v s (s, a) , (3.62)
which can be expanded to .

d ds J(s, a) = -J(s, a) v[x(s, a)] • ∇v s (s, a) v s (s, a) 2 , ( 3
(3.66)

Which leads to

P(v, s) = 1 V 0 Ω(s) v(a) v(a) v(x) v(a) δ [v -v e (x(s, a))]dx, (3.67) P(v, s) = 1 V 0 Ω(s) v(x) v(a) δ [v -v e (x(s, a))]dx, (3.68)
and because of the Dirac-delta,

P(v, s) = 1 V 0 v v(a) Ω(s) δ [v -v e (x(s, a))]dx. (3.69) Furthermore, we have 1 v(a) = V 0 v(x) V (s)
, which can be seen as follows. First, we observe that

Ω 0 v(a) v(a) da = V 0 (3.70)
by definition. Using the transformation a → x(s; a), we can obtain

Ω(s) J -1 v(a) v(a) dx = Ω(s) v(x) v(a) v(a) v(a) dx = Ω(s) v(x) v(a) dx = V 0 , (3.71) which implies that v(a) = V (s) V 0 1 V (s) Ω(s) v(x)dx = V (s) v(x) V 0 . (3.72)
Inserting this in Eq. 3.69 gives

P(v, s) = 1 V 0 v v(x) V 0 V (s) Ω(s) δ [v -v s (s, a)]dx = v v(x) P e (v, s). (3.73)
This shows that the ensemble s-Langrangian velocity PDF is related to the s-Eulerian velocity PDF through flux weighting.

3.6.3.2 Relation Between the t-Eulerian Velocity PDF P e (v,t) and the Ensemble t-Lagrangian Velocity PDF P(v,t)

Keeping the same spirit we derive the relation between the Eulerian temporal velocity PDF Pe (v,t) and the ensemble temporal velocity PDF P(v,t). Starting with the definition of P(v,t):

P(v,t) = 1 V 0 Ω 0 δ [v -v t (t, a)]da, (3.74)
3.6. APPENDIX 79 and changing variable according to the map a → x(t; a) we obtain

P(v,t) = 1 V 0 Ω(t) J -1 δ (v -v e [x(t, a)])dx, (3.75) P(v,t) = 1 V 0 Ω(t) δ (v -v e [x(t, a)])dx, (3.76)
because the map is volume conserving and therefore J = 1. And again from volume conservation V 0 = V (t) which gives for Eq. (3.76)

P(v,t) = 1 V (t) Ω(t) δ (v -v e [x(t, a)])dx = Pe (v,t).
(3.77) 3.6.3.3 Relation Between the Ensemble s-Lagrangian Velocity PDF P(v, s) and the Ensemble t-Lagrangian Velocity PDF P(v,t)

Then, we can derive the relation between ensemble s-Lagrangian velocity PDF P(v, s) and the ensemble t-Lagrangian velocity PDF P(v,t), starting with the definition of P(v,t) we have: (3.80)

p(v,t) = Ω 0 δ (v -v e [x(t), a)])ρ(a)da, (3.78)
And then using the fact that

δ [ f (x)] = ∑ i 1 f (x i ) δ (x -x i ), (3.81) we obtain δ [s -s(t, a)] = 1 v e [x(s; a)] δ [t -t(s, a)]. (3.82)
Inserting this in Eq. 3.80 leads to where p(v,t, s) is the joint PDF of velocity v s (s, a) and particle time t(s, a),

p(v,t) = Ω 0 1 v e [x(s; a)] δ [t -t(s, a)]]δ [v -v e [x(s, a)]]ρ(a)dads, (3.83) = 1 v Ω 0 δ [t -t(s, a)]δ [v -v e [x(
p(v,t, s) = δ [v -v s (s, a)]δ [t -t(s, a)]ρ(a)da. (3.86)
Chapter 4

Upscaling of Anomalous Pore-Scale Dispersion

Introduction

Upscaling hydrodynamic transport is a critical step for modeling solute dispersion in porous media. Since the pioneering works of de Josselin de Jong (1958) and [START_REF] Saffman | A theory of dispersion in a porous medium[END_REF], different approaches have been used for deriving dispersion coefficients and advectiondispersion models for asymptotic spreading and mixing in heterogeneous media [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Brenner | Macrotransport Processes[END_REF][START_REF] Salles | Taylor dispersion in porous media. determination of the dispersion tensor[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF]. However, the asymptotic regime is often not reached for space and time scales relevant at the laboratory scale, for environmental or industrial applications [START_REF] Levy | Measurement and analysis of non-Fickian dispersion in heterogeneous porous media[END_REF][START_REF] Le Borgne | Non-fickian dispersion in porous media: 2. model validation from measurements at different scales[END_REF][START_REF] Moroni | Analysis of dispersion in porous media via matched-index particle tracking velocimetry experiments[END_REF], which makes it important to account for preasymptotic transport, which can in general not be characterized by constant hydrodynamic dispersion, and thus may be termed non-Fickian or anomalous. Non-Fickian preasymptotic dispersion is caused by incomplete mixing on the support scale and thus incomplete sampling of the velocity heterogeneity due to spatial heterogeneity, which characterizes natural systems [START_REF] Dentz | Temporal behavior of a solute cloud in a heterogeneous porous medium: 1. point-like injection[END_REF][START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF](Dentz et al., , 2011b;;[START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Nicolaides | Anomalous physical transport in complex networks[END_REF][START_REF] Wood | The role of scaling laws in upscaling[END_REF]Le Borgne et al., 2011). The conditions under which the behavior can asymptotically be described by hydrodynamic dispersion and the transition to such a regime was discussed in [START_REF] Salles | Taylor dispersion in porous media. determination of the dispersion tensor[END_REF], [START_REF] Dentz | Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport[END_REF], and [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF]. For systems characterized by large or infinite Péclet numbers, non-Fickian behavior may be related to a broad distribution of velocity values. Recent pore-scale transport studies [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF][START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF][START_REF] Carrel | Pore-scale hydrodynamics in a progressively bioclogged three-dimensional porous medium: 3-d particle tracking experiments and stochastic transport modeling[END_REF] showed that observed intermittency of temporal velocity series along individual streamlines are closely related to the occurrence of anomalous dispersion.

INTRODUCTION

85

The quantification of pre-asymptotic dispersion and its causes in the medium and flow properties is a critical issue for upscaling hydrodynamic transport from the pore to the Darcy scale. Pre-asymptotic (non-Fickian) dispersion on the pore and Darcy scales have been modeled by a variety of non-local approaches (Neuman and Tartakovsky, 2009a), such as the multirate mass transfer (MRMT) approach [START_REF] Haggerty | Multiple-rate mass transfer for modeling diffusion and[END_REF][START_REF] Carrera | On matrix diffusion: formulations, solution methods and qualitative effects[END_REF], volume averaging and two-equation formulations for transport [START_REF] Cherblanc | Two-domain description of solute transport in heterogeneous porous media: Comparison between theoretical predictions and numerical experiments[END_REF][START_REF] Davit | Equivalence between volume averaging and moments matching techniques for mass transport models in porous media[END_REF][START_REF] Porta | Upscaling solute transport in porous media from the pore scale to dual-and multicontinuum formulations[END_REF], the continuous time and time-domain random walk approaches [START_REF] Berkowitz | On characterization of anomalous dispersion in porous and fractured media[END_REF][START_REF] Dentz | Transport behavior of a passive solute in continuous time random walks and multirate mass transfer[END_REF][START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF][START_REF] Wright | Upscaling mixing in highly heterogeneous porous media via a spatial Markov model[END_REF][START_REF] Sund | Pre-asymptotic transport upscaling in inertial and unsteady flows through porous media[END_REF][START_REF] Bibliography Sund | Upscaling of dilution and mixing using a trajectory based spatial Markov random walk model in a periodic flow domain[END_REF][START_REF] Sherman | Relationship between pore size and velocity probability distributions in stochastically generated porous media[END_REF], see also the recent review by [START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF]. A critical step for implementing these non-local models concerns the relation between the velocity statistics that are controlled by the pore-scale structure, and the macroscopic transport process. et al., 2013;Siena et al., 2014;[START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF][START_REF] Jin | Statistics of highly heterogeneous flow fields confined to three-dimensional random porous media[END_REF][START_REF] Matyka | Power-exponential velocity distributions in disordered porous media[END_REF][START_REF] Bibliography De Anna | Prediction of velocity distribution from pore structure in simple porous media[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF].

A strategy to systematically upscale (advective) transport from the pore to the Darcy scale consists in identifying the stochastic dynamics of particle velocities, formalizing the link between Lagrangian and Eulerian statistics, and relating the flow statistics to statistical pore-scale properties. [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] and Puyguiraud et al. (2019a) linked observed intermittent patterns in the temporal velocity series to the spatial persistence of pore-scale velocities. These authors showed that velocity series sampled equidistantly along streamlines do not exhibit such intermittent patterns and can be represented by a spatial Markov process. Puyguiraud et al. (2019a) showed for the Berea sandstone sample under consideration in this paper that the spatial velocity series can be represented by an ergodic and stationary Markov process at the sample scale. This observation implies that upscaled transport can be understood and modeled in terms of time-domain or continous time random walks.

In this paper we use the representation of equidistant particle velocities as a stationary Markov process to upscale particle motion and solute transport in the framework of timedomain random walks in terms of the pore-scale velocity distribution and characteristic length scale. We employ two velocity Markov models of different complexity. The first CHAPTER 4. UPSCALING OF ANOMALOUS PORE-SCALE DISPERSION is based on a Bernoulli process for the prediction of the velocity series, the second on an Ornstein-Uhlenbeck velocity process for the normal scores of velocity [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF]Puyguiraud et al., 2019a). The resulting time-domain random walk models are used to predict breakthrough curves, displacement mean and variance as well as the full spatial particle distributions or propagators from full three-dimensional flow and particle tracking simulations for a Berea sandstone sample.

The paper is organized as follows. The methodology we use is detailed in Section 4.2. We specify the flow equation and the transport equation that are solved, and summarize the details about the image acquisition, the flow simulation and the particle tracking simulations. Then, we present the stochastic particle model and describe the parameterization of the velocity process. In Section 4.3, we compare the transport data of the three-dimensional direct particle tracking simulations to the predictions of the upscaled models for uniform and flux weighted injection modes. The conclusions are presented in Section 4.4.

Methodology

In this section, we first present the basic equations for the three-dimensional direct numerical simulation (DNS) of flow and particle motion (particle tracking simulations) at pore scale. Then, we detail the upscaling methodology in the framework of a stochastic model and finally, we provide a summary of the numerical methodology.

Flow and Particle Motion

The Navier-Stokes momentum balance equation is classically used to model pore-scale flow v(x) of an incompressible fluid. At low values of the Reynolds number, the inertial forces are negligible in comparison to the viscous forces and the momentum equation reduces to, [START_REF] Leal | Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes[END_REF]: We consider purely advective transport. Thus, the trajectory of a particle originally located at x(t = 0, a) = a is described by

∇ 2 v(x) = 1 ν ∇p(x), ( 4
dx(t, a) dt = v[x(t, a)], (4.2)
where v[x(t, a)] is the Lagrangian velocity. Its magnitude is v t (t, a) = v[x(t, a)] . The travel distance s(t, a) along a particle trajectory until time t and the travel time t(s, a) up to a streamwise distance t(s, a) are given by

ds(t, a) dt = v t (t, a), dt(s, a) ds = 1 v s (s, a) , (4.3a)
where we defined v s (s, a) = v t [t(s), a]. We perform a variable change from time to streamwise distance, which renders time as t(s, a) a dependent variable. The transform from t → s implies setting dt = ds/v s (s, a) in Eq. (4.2). This gives for the particle position as a function of distance s the evolution equation

d x(s, a) ds = ω ω ω(s, a), ω ω ω(s, a) = v[x(s, a)] v s (s, a) , (4.3b)
where ω ω ω(s, a) denotes the unit vector in the flow direction, v[x(s, a)] is denoted the s(pace)-Lagrangian velocity [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]Puyguiraud et al., 2019a) because it is the particle velocity at a given spatial distance s along the particle trajectory, and its magnitude v[x(s, a)] is equal to v s (s, a). We will refer in the following to v s (s) simply as particle velocity. Equation (4.3) describes the motion of a particle along a given streamline as a time-domain random walk [START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum transport models[END_REF][START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF] in that particles perform transitions over a fixed streamwise distance in variable time, which depends on the local velocity. Particle motion can be solved alternatively by integrating Eq. (4.2) in time or by integrating the system of equations (4.3) in streamwise distance. The numerical simulations performed in this paper use the former, the upscaling methodology presented in the next section uses the latter. The distribution of initial particle positions is denoted by ρ(a). We consider here two different initial particle distributions, uniform and flux-weighted, in order to probe the impact of the initial condition on average particle transport. The uniform and flux-weighted initial distribution read as

ρ(a) = 1 V 0 I(a ∈ Ω 0 ), ρ(a) = v e (a) Ω 0 dxv e (x) I(a ∈ Ω 0 ), (4.4)
where Ω 0 is the injection domain, V 0 its volume; I(•) is the indicator function which is 1 if its argument is true and 0 otherwise.

In the following, we study the transport behavior in terms of breakthrough curves and spatial particle distributions and subsequently quantify dispersion by analyzing the first and second displacement moments. The breakthrough curve at a control plane located at position x 1 in the mean flow direction, is defined in terms of the first passage time

τ(x 1 , a) = min[t|x 1 (t, a) ≥ x 1 ]. (4.5)
where x 1 (t, a) denotes the position of particle a after time t in mean flow direction. The breakthrough curve is equal to the PDF of the first passage times,

f (t, x 1 ) = Ω 0 daρ(a)δ [t -τ(x 1 , a)] . (4.6)
The breakthrough curves contain information on the residence times within the rock sample or the volume between the inlet and control plane, and the concentration in the effluent fluid. This information is useful for modeling reactive transport, for instance for applications to design aquifer decontamination or model laboratory dissolution-precipitation experiments. Furthermore, we consider the spatial particle distribution, also called propagator, which is defined by

g(x 1 ,t) = Ω 0 daρ(a)δ [x 1 -x 1 (t, a)] . (4.7)
This quantity gives information on the dispersion of a solute or particle cloud. Likewise this information can be used in the modeling of reactive transport and deployment of a reactant species, as well as for assessment of propagators in NMR imaging of flow and transport in porous media. The mean displacement and displacement variance are defined by

m 1 (t) = Ω 0 daρ(a)x 1 (t, a), (4.8) σ 2 (t) = Ω 0 daρ(a) [x 1 (t, a) -m 1 (t)] 2 . (4.9)
They measure the center of mass position and spatial variance of the particle distribution g(x 1 ,t). The spatial variance is a measure for hydrodynamic dispersion. If its asymptotic evolution is linear, its growth rate is equal to the hydrodynamic dispersion coefficient.

Stochastic Model

We formulate a stochastic model for particle motion in the mean flow direction x 1 of the coordinate system based on the stochastic representation of the s-Lagrangian velocity magnitude v s (s, a) as a stationary and ergodic Markov process v s (s) (Puyguiraud et al., 2019a).

The Markov process for v s (s) is characterized by the PDF r(v, ss |v ) to make a transition from v = v s (s ) at distance s to v = v s (s) at s > s and the steady state distribution p s (v). Both the velocity PDF p(v, s) and the transition PDF r(v, s|v ) satisfy the Chapman-Kolmogorov equation [START_REF] Risken | The Fokker-Planck Equation[END_REF])

p(v, s) = ∞ 0 r(v, s -s |v )p(v , s )dv . (4.10)
The transition probability converges to the steady-state distribution in the limit of s c with c a characteristic velocity correlation scale,

lim s→∞ r(v, s|v ) = p s (v). (4.11)
The characteristic correlation scale c ≈ 2.5 p where p is the characteristic pore length. Eq.

(4.11) implies that the distribution converges to the steady-state PDF P(v) independently of the initial condition p 0 (v). The steady state distribution p s (v) is related to the Eulerian velocity PDF p e (v) through flux weighting [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]Puyguiraud et al., 2019a)

p s (v) = vp e (v) v e , (4.12)
where v e is the mean Eulerian velocity. In this framework, the irregular particle motion described by (4.3) is represented by the stochastic evolution equations

d x1 (s) ds = χ -1 , dt(s) ds = 1 v s (s) , (4.13a)
where x1 (s) indicates the position of the particle in the mean flow direction (denoted by the subscript 1 similarly to Section 4.2.1). Note that the displacement rate ω 1 (s) in 1-direction in general fluctuates with s. We represent it here by its average ω 1 (s) = χ -1 , where χ is the advective tortuosity given by [START_REF] Koponen | Tortuous flow in porous media[END_REF][START_REF] Ghanbarian | Tortuosity in porous media: a critical review[END_REF])

χ = v e (x) v 1 (x) . (4.13b)
The advective tortuosity compares the distance s along the streamline with the average linear distance in the mean flow direction x 1 (s) , see Appendix 4.5.1 for details. It is an

CHAPTER 4. UPSCALING OF ANOMALOUS PORE-SCALE DISPERSION

indicator fort the complexity of the pore space and gives information on advective excursions transverse to the mean flow direction.

The stochastic model (4.13) belongs to the continuous time random walk or timedomain random walk class of models because the time increment varies between random walk steps [START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF]. Different initial particle distributions ρ(a) are in this framework quantified in terms of the corresponding initial velocity distribution p 0 (v), this means in terms of the PDF of velocities v e (x) in the injection domain Ω 0 ,

p 0 (v) = 1 V 0 Ω 0 daρ(a)δ [v -v e (a)] . (4.14)
The breakthrough curve f (t, x 1 ) is in this framework given by

f (t, x 1 ) = δ [t -t(x 1 χ)] , (4.15)
where the angular brackets denote the ensemble average over all particles and t(x 1 χ) denotes the travel time over the streamlines distance s = x 1 χ. The mean displacement and its variance are

m 1 (t) = x1 [s(t)] , (4.16) σ 2 (t) = [ x1 [s(t)] -m 1 (t)] 2 , (4.17)
where s(t) = max [s|t(s) ≤ t]. The spatial particle distribution is accordingly given by

g(x 1 ,t) = δ (x 1 -x1 [s(t)]) . (4.18)
In the following, we briefly review two Markov processes, which model the evolution of p(v, s) from arbitrary initial conditions, a Bernoulli velocity process [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF] and an Ornstein-Uhlenbeck process [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] for the evolution of the normal scores of velocity.

Bernoulli Process

This Markov process for the prediction of the fluid particle velocities is modeled as a Bernoulli process where the velocity changes after a distance ∆s according to a Bernoulli trial. This means that the particle velocity v s (s) does not change with probability p B (∆s) = exp(-∆s/ c ) and changes randomly with probability 1-p B (∆s) to a velocity which is sampled from the steady state PDF p s (v). The transition probability r(v, ∆s|v ) is then expressed by [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF]

r(v, ∆s|v ) = exp(-∆s/ c )δ (v -v ) + [1 -exp(-∆s/ c )]P(v). (4.19)
The Bernoulli process reproduces qualitatively the evolution of the Lagrangian velocity statistics, but underestimates the convergence rate of p(v, s) toward its steady state at low velocities (Puyguiraud et al., 2019a). The Bernoulli process does not account for any velocity dependence of the decorrelation rate and is therefore not capable of capturing the faster decorrelation of low velocities. Nevertheless, we consider the Bernoulli process as a possible evolution model for the particle velocities due to its simplicity. In the following, we refer to this model as Bernoulli model.

Ornstein-Uhlenbeck Process

This velocity Markov process considers the evolution of the normal scores

w(s) = Φ -1 (Π [v s (s)]) ≡ F[v s (s)], (4.20)
where Φ(w) is the cumulative Gaussian distribution and Π(v) the cumulative steady state velocity distribution,

Φ(w) = 1 + erf(w/ √ 2) 2 , Π(v) = v 0 p s (v )dv , (4.21)
where p s (v) is the steady state velocity distribution. The normal scores w(s) follow the Ornstein-Uhlenbeck process [START_REF] Gardiner | Stochastic Methods[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] dw(s) ds = --1 c w(s) + 2 -1 c ξ (s), (4.22) where ξ (s) is a Gaussian white noise characterized by zero mean ξ (t) = 0 and covariance ξ (s)ξ (s ) = δ (ss ). The Ornstein-Uhlenbeck process is a mean reverting process. In the absence of the noise terms, w(s) relaxes exponentially fast towards 0. In the presence of noise, there is a steady state for s c , at which w(s) ∼ ξ (s). The transition PDF for the Ornstein-Uhlenbeck process is given by [START_REF] Gardiner | Stochastic Methods[END_REF]

r w (w, s|w ) = exp -[w-w exp(-s/ c )] 2 2[1-exp(-2s/ c )] 2π [1 -exp(-2s/ c )] . (4.23)
The velocity values v s (s) are obtained from w(s) at any distance s through the Smirnov transform [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF]

v(s) = Π -1 (Φ[w(s)]) ≡ F -1 [w(s)]. (4.24)
The velocity transition PDF r(v, s|v ) is thus given by In the following, we refer to this model as the OU model.

r(v, s|v ) = r w [F(v), F(v )] dF(v) dv . ( 4
The process (4.22) is implemented numerically via an Euler scheme as, (4.26) where w n = w(n∆s) and ξ n is a Gaussian random variable with 0 mean and unit variance. Accurate results are obtained using a discretization ∆s ≤ c /10.

w n+1 = w n --1 c w n ∆s + 2 -1 s ∆sξ n ,

Rock Sample, Flow Field and Velocity Statistics

Here we provide a brief summary of the rock sample, numerical methodology and velocity statistics. Details on the image acquisition and segmentation as well as the flow field computation can be found in [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF]. Details regarding the particle tracking computation can be found in Puyguiraud et al. (2019a).

We use a three-dimensional digitized image (900 3 voxels) of a sample of a Berea sandstone (Upper Berea Sandstone unit, Ohio, USA). Berea sandstone is a sedimentary rock characterized by medium porosity and permeability values as well as medium pore-scale structural heterogeneity compared to common reservoir rocks. Because of these average properties, its simple composition (quasi pure silica) and its remarkable macroscopic homogeneity that allows easy comparisons, Berea sandstones are often used as a reservoir rock standard for experimental/ laboratory works. The image is reconstructed from X-ray microtomography [START_REF] Paganin | Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object[END_REF][START_REF] Sanchez | Threedimensional synchrotron virtual paleohistology: a new insight into the world of fossil bone microstructures[END_REF]. The material density is recorded in a raw grey level image, which is is segmented [START_REF] Smal | An automatic segmentation algorithm for retrieving sub-resolution porosity from x-ray tomography images[END_REF] in order to obtain a binary image mapping the solid and the conne porosity. The porosity of the sample is 0.18. The voxel length is 10 -6 m. The average pore length is 1.5 • 10 -4 m.

The steady-state Navier-Stokes equations are solved using the SIMPLE method implemented in OpenFOAM (simpleFoam) [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF] in order to obtain the velocity components at the center of the voxel surface for the full domain. The mean flow velocity is aligned with the 1-direction of the coordinate system and given by v 1 = 4.9 • 10 -4 m/s. The characteristic time scale is given by τ c = p / v 1 = 3 • 10 -1 s. The streamlines starting at any location in Ω 0 are built from the interpolated velocity using quadratic interpolation at the voxel in contact with the solid and linear interpolation elsewhere [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF][START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF]. The injection domain Ω 0 is a box of an extension of 50 voxels in mean flow direction and 900 voxels in the directions perpendicular. Figure 4.1 illustrates a subdomain of the segmented Berea sandstone image and some example trajectories through the pore-space.

The stochastic models for particle motion described in the previous section require the knowledge of the velocity correlation length c , the steady state velocity PDF p s (v) and the initial velocity PDF p 0 (v). Puyguiraud et al. (2019a) performed a full statistical analysis of the velocity statistics of the rock sample under consideration. There, the velocity correlation length c is found to be 2.5 times the characteristic pore length p . The mean Eulerian velocity magnitude is v e = 8.05 • 10 -4 m/s, which give the advective tortuosity χ = v e / v 1 = 1.64. The mean s-Lagrangian velocity magnitude is v s = 3.4 • 10 -3 . The steady state velocity distribution p s (v), the Eulerian velocity distribution p e (v) and the initial velocity distributions for uniform and flux-weighted injections are shown in Figure 4.2. All velocity distributions show a strong tailing toward low velocities. For the flux-weighted injection, the initial velocity PDF p 0 (v) is close to the steady state PDF p s (v), while for the uniform injection, p 0 (v) is close to the Eulerian velocity PDF p e (v) (Puyguiraud et al., 2019a). The next section studies particle transport through the sample using direct numerical simulations of purely advective particle motion, and its upscaling in terms of the velocity correlation length c and velocity PDF p s (v) in the framework of the velocity Markov models discussed in the previous section. 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10

1 v 0 v -0.8 pdf v/ v e p e (v) p 0 (v) p s (v)
Flux weighted p 0 (v)

Results

We study here the upscaling of the purely advective particle motion in the Berea sandstone sample discussed in the previous section. Hydrodynamic transport in Berea sandstones is known to be non-Fickian at the scale of centimeter sized samples [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF][START_REF] Bijeljic | Signature of non-Fickian solute transport in complex heterogeneous porous media[END_REF]. The direct numerical flow and particle tracking simulations represent the reference data. The large scale behavior is measured in terms of particle breakthrough curves at different control planes, the particle displacement variance or dispersion, and the spatial particle distribution or propagators. These behaviors are then compared to the ones predicted by the stochastic particle models presented in the previous section, which quantify the upscaled particle motion.

Breakthrough Curves

The breakthrough curve denotes the residence time distribution of the solute in the domain. It may be used to infer the likeliness of chemical reactions to occur, and to assess the retention or storage potential of the subdomain, for example. Under uniform and homogeneous flow conditions, the BTC at control plane has an inverse Gaussian shape and

RESULTS

95

decays sharply at long times. Under heterogeneous flow conditions breakthrough curves are characterized by early and late particle arrivals. In the following, we consider BTCs for flux-weighted and uniform injection conditions and compare them to the predictions of the upscaled transport models.

Flux-Weighted Injection

In this section we use a flux-weighted injection at the inlet for the computation of the breakthrough curves. The velocity PDF p 0 (v) at the inlet is close to the stationary PDF p s (v), see Figure 4.2. This implies that the particle velocities are approximately stationary.

We compare the breakthrough curves of the direct simulation described in Section 4.2.3 to the two CTRW models described in Section 4.2.2. We perform the simulations using 10 7 particles in the DNS case, while we used 10 9 for the upscaled models. We compute arrival times at distances x 1 = 6 p which corresponds to the end of the sample, x 1 = 36 p , and x 1 = 200 p . To compute the breakthrough curves at distances larger than the sample size, a particle exiting the sample at the outlet is reinjected at the inlet while conserving the velocity continuity (Puyguiraud et al., 2019a).

Figure 4.3 displays the breakthrough curves from the DNS and the two stochastic models. We observe a strong anomalous behavior characterized by early peak arrivals and long tailing at late times. The late time tails display the power-law t -2 at all distances. The exponent can be predicted from CTRW theory because it is directly linked to the behavior of the low velocity part of the steady Lagrangian velocity PDF, see Appendix 4.5.2. The velocity distribution scales as p s (v) ∝ v β -1 with β = 1 for the small values of v. This implies that f (t, x 1 ) ∝ t -2 . The Bernoulli and OU models perform equally well. The early, intermediate, and late times are well captured even if at the closest control plane the two models do not reproduce the first arrivals perfectly. The two models give similar results because the injection velocity PDF is close to the steady-state PDF and therefore the models only need to be able to preserve this distribution over time, which they are both capable of doing (Puyguiraud et al., 2019a).

Uniform Injection

The results are however different when using a uniform injection as initial condition. Under this condition, the particle velocities are non-stationary.

RESULTS

97

as p 0 (v) ∝ v -0.8 , see Figure 4.2. The initial velocities persist within a distance of about c from the inlet. The long time behavior is dominated by strong particle retention in the vicinity of the inlet and dominated by the transition time over the distance c , (4.27) see also Appendix 4.5.2.

ψ 0 (t) = c t 2 p 0 ( c /t) ∝ t -1.2 ,
Both, the Bernoulli and OU models provide good predictions of the breakthrough curves obtained from the DNS. The peak position and width are well captured. Also the behaviors at intermediate and long times are accurately predicted. Both models give the correct long time tailing, while the Bernoulli model slightly overestimates the tail compared to the DNS. This can be traced back to the observation that the Bernoulli model overestimates the persistence of low velocities (Puyguiraud et al., 2019a).

Particle Distribution

In this section we study the evolution of the spatial particle distributions, or propagators g(x 1 ,t) with focus on the differences in the evolution due to the initial particle distribution. The propagator at a given time maps the spatial heterogeneity of the velocity field which controls the spatial distribution of the mass in the system and, for instance, gives information on the localization of reaction with the solid phase. Thus, together with the breakthrough curve it allows for a spatio-temporal characterization of the solute distribution.

Figure 4.5 shows g(x 1 ,t) for uniform and flux-weighted injection conditions at three different times. In both cases, the particle distributions are asymmetric and characterized by a leading edge and long spatial tail. These behaviors are caused by the broad distribution of particle velocities. For the uniform injection, the proportion of particles in low velocity regions is larger than for the flux-weighted injection. Thus, the tailing at short and intermediate times is stronger in the uniform than in the flux-weighted case. With increasing time, the spatial distributions lose the memory of the initial condition and assume the same shape. Note that this is different from the breakthrough curves, whose long time behavior is dominated by the injection condition.

The OU and Bernoulli models predict the spatial profiles under both injection conditions for times t > τ c . For times t < τ c , the stochastic models do not capture the trailing tail in the case of the flux-weighted injection. At times t < τ c , the tail of the spatial distribution in the direct numerical simulation is determined by the velocity components v 1 (x) < c /t

RESULTS

99

in the mean flow direction, this means by particles that persist in their initial velocity. Thus, if

p 1 (v) = Ω 0 daρ(a)δ [v -v 1 (a)] , (4.28)
is the PDF of the 1-component of the particle velocities in the injection domain, the particle distribution at early times is obtained through the variable transform

x 1 = v 1 t as g(x 1 ,t) = t -1 p 1 (x 1 /t). (4.29)
In the upscaled stochastic particle model, likewise, the tail of the spatial distribution is due to the particle velocities that persist in their initial velocity. The distance traveled at the initial velocity here, however, is x 1 = χv 0 t because v 0 is the initial velocity magnitude.

Thus, the upscaled particle model predicts for the early time distribution

g(x 1 ,t) = χt -1 p 0 (χx 1 /t), (4.30)
where p 0 (v) is the PDF of the velocity magnitude in the injection domain. The PDFs of the 1-component and magnitude are in general different, which explains the difference in the tailing behaviors for small times in the case of flux-weighted injection. For the uniform injection, the distributions of the 1-component and the absolute value of velocity are similar in shape, which explains the good match between the stochastic models and the DNS data.

While the stochastic models correctly capture the memory of the injection condition on the evolution of the spatial distribution, we do not expect them to be valid at short times and distances, for which the behaviors depend on the local details of the velocity fluctuations.

Dispersion

In this section, we consider the displacement mean and variance. The evolution of the displacement mean is an indicator of t-Lagrangian stationarity, while the displacement variance gives information on particle dispersion. We have seen in the previous section that the Bernoulli and OU models perform equally well in the prediction of the spatial profiles.

Thus, here, we compare the DNS data for the displacement mean and variance with the prediction of the stochastic particle model based on the OU model only. is the average over the 1-component of the particle velocities in the injection domain. Thus, the slope for the flux-weighted injection is larger than for the uniform injection. At t > τ c , the initial linear behavior crosses over to the long time behavior, which is independent from the injection conditions. For the steady state velocity distribution p s (v) illustrated in Figure 4.2, for which p s (v) ≈ constant at low velocities, CTRW theory predicts m 1 (t) ∝ t/ ln(t) [START_REF] Comolli | Anomalous dispersion in correlated porous media: a coupled continuous time random walk approach[END_REF]. This is confirmed by the DNS data. Both the OU model and the Bernoulli model (not shown) predict the evolution of the mean velocity, with a slight mismatch at short times for the reasons discussed above.

Figure 4.6 shows the evolution of the displacement variance σ 2 (t) for uniform and flux-weighted initial conditions. At early times, the behaviors are ballistic, this means

σ 2 (t) = σ 2 v 1 t 2 , (4.32)
where σ 2 v 1 is the variance of p 1 (v). As, for the mean, the displacement variance is larger for the flux-weighted than for the uniform injection. At t > τ c , the variance crosses over from the ballistic toward the asymptotic regime. For t τ c , CTRW theory predicts σ 2 (t) ∝ t 2 / ln(t) 3 [START_REF] Comolli | Anomalous dispersion in correlated porous media: a coupled continuous time random walk approach[END_REF]. The behavior is superdiffusive. Both the OU and Bernoulli (not shown) models predict the evolution of the displacement variance with a mismatch in the ballistic early time behaviors because the stochastic models are determined by the statistics of the velocity magnitude.

Conclusions

Using direct three-dimensional pore-scale simulations of flow and transport in a sample of Berea sandstone as a reference case, we have shown that the upscaling of pore-scale dispersion can be accurately performed using a stochastic approach based on velocity Markov models for equidistantly sampled particle velocities. The upscaled model is implemented in the framework of a time-domain or continuous time random walk approach, which describes particle motion in equidistant spatial steps with random transition times. The presented modeling approach is predictive in the sense that it depends on the Eulerian velocity distribution and advective tortuosity, both flow attributes, and the average pore length which is a medium attribute. It is worth noticing that this dependence allows, in turn, inferring information on the velocity statistics and pore length from (experimental) observations of breakthrough curves, spatial particle distributions and/or displacement moments.

Our analysis has shown that the observed transport behaviors are sensitive to the initial distribution of the tracer particles. Breakthrough curve tailing, for example, can depend on the initial velocity distribution. Under this condition, the breakthrough curve tail gives 102 CHAPTER 4. UPSCALING OF ANOMALOUS PORE-SCALE DISPERSION information on the steady-state velocity distribution in the sample only if the injection domain is already representative (large enough) for the initial velocity distribution to be equal to the stationary PDF, otherwise, the breakthrough curve gives information on the local velocities in the injection domain. Similarly, the spatial particle distribution depends at short and intermediate times on the injection condition. At late times, however, the memory of the initial condition diminishes and the shape becomes independent from the injection condition. This is also reflected in the displacement mean and variance. The early time behaviors of the displacement mean and variance give information on the velocity mean and variance in the injection domain. At late times, the displacement mean shows slightly sublinear behavior, the variance being superlinear, which is due to the tail of the steady state velocity distribution toward low velocity values. In this sense, the evolution of the moments can be seen as a scan through the velocity PDF. At short times, it is dominated by the high and intermediate velocity values, which determine the velocity mean and variance, at long times by the low velocities. The upscaled stochastic particle models can be conditioned on the injection condition through the distribution of initial particle velocities and is able to predict the dependence on the initial condition and full evolution of particle dispersion.

We consider two velocity Markov models, the Bernoulli and Ornstein-Uhlenbeck models, which both are parameterized by the velocity correlation length and steady s-Lagrangian velocity PDF. While both models predict the evolution from an initial velocity PDF toward the steady state, they differ in the convergence rates as discussed in (Puyguiraud et al., 2019a). Both processes predict the transport behavior and dependence on the initial distribution, which indicates that here the details of the evolution are secondary compared to the fact that there is an evolution.

The presented analysis and the derived stochastic particle models consider purely advective transport. Thus, they are directly relevant for transport scenarios characterized by high Péclet numbers, such as solute transport at high flow rates and passive particles characterized by low diffusion coefficients. The stochastic model is based on a Markov model for the streamwise particle velocity, this means that velocities are sampled advectively at a constant frequency in space. The breakthrough curve tailing, for example, is due to the persistence of low velocities over a constant length scale, the pore length. For finite Péclet numbers, particle velocities may be decorrelated due to diffusion across streamlines and low advective transition times may be cut off at the characteristic diffusion time. Thus, we expect anomalous behavior to persist in an intermediate regime depending on the Péclet number and to transition towards normal behavior at times larger than the characteristic pore-scale diffusion time.

APPENDIX 103

Hydrodynamic dispersion and other pore-scale phenomena have their origins in porescale velocity fluctuations. Thus, the presented upscaled stochastic model and the associated parameterization can serve as a basis for the systematic quantification of the impact of pore-scale velocity fluctuations on Darcy scale transport phenomena.

Appendix

Tortuosity

We derive here the average of the ω 1 (s, a) along a streamline under ergodic conditions. To this end, we first note that the position x 1 (s, a) can be written by integration of (4.3) as

x 1 (s, a) = s 1 s s 0 ω 1 (s , a)ds . (4.33)
The expression in the square brackets denotes the average of ω 1 (s, a) along a particle trajectory. At the same time, it denotes the ratio of linear to streamwise distance,

ω 1 (s, a) s = lim s→∞ 1 s s 0 ω 1 (s , a)ds = x 1 (s, a) s , (4.34)
where the angular brackets with subscript s denote the streamwise average along a trajectory. The average of ω 1 (s, a) over an ensemble of particles is defined by

ω 1 (s, a) = lim V 0 →∞ 1 V 0 Ω 0 v 1 [x(s, a)] v e [x(s, a)] ρ(a)da. (4.35)
We consider a flux-weighted initial condition, see (4.4). Under ergodic conditions, this initial condition corresponds to the steady state velocity PDF p s (v), which is equal to the flux-weighted Eulerian velocity PDF. This can be seen by using

ρ(a) = 1 V 0 v e (a) v e (x) I(a ∈ Ω 0 ), (4.36)
in the limit V 0 → ∞. Also, [START_REF] Koponen | Tortuous flow in porous media[END_REF] pointed out that it is natural for porous media to consider a flux-weighted average, see also [START_REF] Ghanbarian | Tortuosity in porous media: a critical review[END_REF]. Furthermore, under ergodic conditions, the average over a single particle trajectory is equal to the average over the initial ensemble of particles and so

ω 1 (s, a) s = ω 1 (s, a) = x 1 (s, a) s = χ -1 . (4.37)
Using expression (4.36) in (4.35), we obtain

ω 1 (s, a) = lim V 0 →∞ 1 V 0 Ω 0 v 1 [x(s, a)] v e [x(s, a)] v e (a) v e (x) da. (4.38) CHAPTER 4.

UPSCALING OF ANOMALOUS PORE-SCALE DISPERSION

In order to evaluate this expression, we perform the variable change a → x(s, a),

ω 1 (s, a) = lim V 0 →∞ 1 V 0 Ω(s) v 1 [x(s, a)] v e [x(s, a)] v e (a) v e (x) J(a, s) -1 dx, (4.39)
where J(a, s) is the Jacobian of the transformation. It can be determined by noting that (Batchelor, 2000, p. 75)

dJ(a, s) ds = J(a, s)∇ • v[x(s, a)] v e [x(s, a)] . (4.40)
This differential equation can be integrated by noting that ∇ • v(x) = 0 and

dv e [x(s, a)] ds = ∇v e [x(s, a)] • v[x(s, a)], (4.41)
which follows by using the chain rule and (4.3). Thus, we obtain for the initial condition

J(a, s = 0) = 1 that J(a, s) = v e (a) v e [x(s, a)] . (4.42)
Inserting this expression into (4.38) gives

ω 1 (s, a) = lim V 0 →∞ 1 V 0 Ω(s) v 1 (x) v e (a) dx = v 1 v e . (4.43)
This result is consistent with [START_REF] Koponen | Tortuous flow in porous media[END_REF]. This implies that at s p , we can set

ω 1 (s, a) = χ -1 = v 1 v e .
(4.44)

Continuous Time Random Walk

For transition length of the order of the correlation length c , subsequent particle velocities can be considered independent and thus, the space-time particle motion (4.13a) may be approximated by

x n+1 = x n + c χ , t n+1 = t n + τ n , (4.45)
where x n = x(s n ) with s n = n c . The random transition time τ n is given by

τ n = c v s (s n ) . (4.46)
The time increments for n > 0 is distributed as

ψ(t) = c t 2 p s ( c /t). (4.47) 4.5. APPENDIX 105
For n = 0, the transition time PDF is distributed according to

ψ 0 (t) = c t 2 p 0 ( c /t). (4.48)
Under steady state conditions, this means for p 0 (v) = p s (v) and thus ψ 0 (v) = ψ(v), Equations (4.45) describe a continuous time random walk as discussed in [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF]. Thus, the asymptotic behavior of the breakthrough curves and displacement moments can be predicted based on the scalings of the transition time distribution. For ψ(t) ∝ t -1-β at large times, the breakthrough curves scales as f (t, x 1 ) ∝ t -1-β , the mean displacement scales as m 1 (t) ∝ t and the displacement variance as σ 2 (t) ∝ t 3-β . Note that this scaling for ψ(t) implies that the velocity distribution p s (v) ∝ v β -1 at small velocities.

CHAPTER 5

Is there a Representative Elementary Volume for Anomalous Dispersion? existence of scale separation of spatial medium fluctuations. An REV cannot be defined for continuously hierarchized heterogeneous media such as fractal materials. The REV is typically determined from applying this concept to the material microstructure and specifically to its simplest quantitative notion, which is its porosity. Porosity is easily measurable at laboratory scale and can also be determined at pore-scale by using imaging methods such as computed microtomography which allows characterizing the micro-structures over volumes that are typically larger than the REV. Since the REV can be well-defined for porosity, it is generally assumed that this definition also implies the existence of transport relevant parameters such as the specific discharge and the hydrodynamics dispersion coefficients. The former represents the mean pore velocity, the latter quantifies its fluctuations.

The values of these parameters are considered to be well-defined and constant on the REV scale. Note that assuming the specific discharge constant within the REV implies that the product of permeability and the pressure gradient is constant. This is implied by the Darcy equation, which states

q = - k µ dP(x 1 ) dx 1 , (5.2)
where k is permeability, P(x 1 ) is pressure and µ dynamic viscosity. The assumption that the REVs for the porosity and the specific discharge are the same is not evident. The assumption that the REVs for the porosity and for the hydrodynamics dispersion coefficients 112 CHAPTER 5. IS THERE AN REV FOR ANOMALOUS DISPERSION?

are similar is even less evident because the later encompasses the impact of pore-scale velocity fluctuations. Yet, if this assumption holds, average solute transport can be described by the advection-dispersion equation (ADE) [START_REF] Bear | Dynamics of fluids in porous media[END_REF])

φ ∂ c(x 1 ,t) ∂t + q ∂ c(x 1 ,t) ∂ x 1 -D ∂ 2 c(x 1 ,t) ∂ x 2 1 = 0. (5.3)
This approach described Darcy-scale transport in terms of porosity φ , specific discharge q and the hydrodynamic dispersion coefficient D.

Experimental [START_REF] Moroni | Statistical mechanics with three-dimensional particle tracking velocimetry experiments in the study of anomalous dispersion. II. experiments[END_REF]Cortis and Berkowitz, 2004;[START_REF] Holzner | Intermittent Lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] and numerical [START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF][START_REF] Bijeljic | Signature of non-Fickian solute transport in complex heterogeneous porous media[END_REF][START_REF] Liu | Applicability of the Dual-Domain Model to Nonaggregated Porous Media[END_REF][START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF][START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF]Puyguiraud et al., 2019a;[START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF] porescale studies observed deviations from predictions based on the ADE (5.3). This includes tailing in solute breakthrough curves, non-linear growth of dispersion and non-Gaussian particle distributions and propagators. Such behaviors were modeled based on non-local transport approaches such as multirate mass transfer and continuous time random walks [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Noetinger | Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale[END_REF] as well as fractional dynamics [START_REF] Cushman | Statistical mechanics with three-dimensional particle tracking velocimetry experiments in the study of anomalous dispersion. i. theory[END_REF].

In this paper, we investigate the notion of REV for non-Fickian dispersion. In this study, we scrutinize the assumptions underlying modeling approaches for Fickian and non-Fickian dispersion and the relation with the notion of the REV. This study is organized as follows. In Section 5.2.2.1 we discuss the bases of the ADE framework and we include a critical revision of its limitations. In Section 5.2.2.2 we discuss the framework of continuous time random walk models and investigate what underlying assumptions they are relying on. This leads us to define, in Section 5.2.3 a Eulerian REV for anomalous transport and to illustrated its evaluation from computations performed using a digitized volume of a real rock sample in Section 5.2.5. In Section 5.3.1 we successfully predict the transport on this sample with a CTRW for purely advective transport while in Section 5.3.2 we elaborate on the processes to account for once the diffusion is introduced. In Section 5.3.3 we discuss the convergence to asymptotic dispersion for different Péclet regimes.

Dispersion Upscaling and the Representative Elementary Volume

In this section, we consider the assumptions that underlay descriptions of solute dispersion by advection-dispersion models and continuous time random walks. From these considerations we propose the definition of an REV in terms of the Eulerian flow statistics, and discuss conditions on the Lagrangian velocity statistics.

DISPERSION UPSCALING AND THE REV 113

Porescale Flow and Transport

Porescale flow is described here by the Stokes equation

∇ 2 u(x) = - ∇p(x) µ , (5.4)
together with incompressibility ∇ • u(x) = 0. The mean pressure gradient is aligned with the 1-direction of the coordinate system. We consider here purely advective transport that can be described by the kinematic equation

dx(t, a) dt = v(t, a), (5.5)
where x(t = 0, a) = a and v(t, a) = u[x(t, a)] is the Lagrangian velocity. We disregard here diffusion and focus on particle advective motion along streamlines as the only mechanism by which the velocity field can be sampled. The impact of diffusion on the results presented in the following are discussed below. Equation (5.5) may be transformed into streamlines coordinates t → s, where

ds(t) dt = v(t), (5.6)
where v(t) = |v(t)|. This gives the equivalent system of equations

dx(s, a) ds = v(s, a) v(s, a) , dt(s, a) ds = 1 v(s, a)
.

(5.7)

Dispersion

We consider the conditions under which pore-scale velocity fluctuations can be quantified by the concept of hydrodynamic dispersion and how this relates to the notion of an REV.

Then, we discuss the same questions for continuous time random walk models to upscale anomalous dispersion.

Fickian Dispersion

In order to identify the basic assumption underlying the ADE formulation of Darcy-scale transport, we consider the Langevin equation equivalent to the ADE (5.3), which is given by [START_REF] Gardiner | Stochastic Methods[END_REF])

dx(t) dt = v 1 + v 1 (t), (5.8)
where we decomposed the particle velocity v 1 (t) into its mean v 1 and fluctuation v 1 (t). The mean pore velocity is v 1 = q/φ . The velocity fluctuation v 1 (t) is represented by a stationary Gaussian random process characterized by 0 mean and the covariance function

v 1 (t)v 1 (t ) = 2Dδ (t -t ), (5.9)
where D is the hydrodynamic dispersion coefficient and δ (t) is the Dirac Delta. The angular bracket denotes the average over all noise realizations. With these properties of the velocity fluctuations, Equation (5.8) describes Brownian dynamics.

The representation of the velocity fluctuations as a δ -correlated Gaussian process is based on several conditions. First, the velocity process needs to be stationary and ergodic. This means that its mean and variance depend only on the time lag and not on the absolute time. Second, velocity fluctuations decay exponentially fast on a characteristic correlation time scale τ c . Furthermore, based on the assumption that the velocity distribution has finite variance, the displacement distribution, which is the sum of random velocity increments, converges towards a Gaussian distribution as time increases. This is a consequence of the central limit theorem [START_REF] Gardiner | Stochastic Methods[END_REF] and warrants the modeling of the statistics of v 1 (t)

as Gaussian. The correlation model (5.9) is valid at observation times that are much larger than the correlation scale τ c , which can be related to the characteristic advection time over a characteristic length scale c ,

τ c = c v 1 .
(5.10)

This implies that for time t τ c , particles must have access to the full spectrum of velocity variability. The Langevin equation (5.8), which is valid at t τ c , thus implies that at each random walk step, particles can sample the full spectrum of random velocities. Particles become statistically equal on the time scale τ c . This temporal notion can be related to a spatial REV scale through the length scale c that is assumed to mark the correlation time together with the mean velocity v 1 . Thus, the REV scale is supposed to contain a representative set of flow velocities that particles can sample with equal probability. This is discussed further in Section 5.2.3.

Anomalous Dispersion

As outlined in the previous section, Brownian dynamics describe dispersion at times that are much larger than a typical correlation scale τ c , which is equal to the transition time over an average pore length by the mean flow velocity. In order to scrutinize this condition, let because v(s ) can be considered approximately constant over the length of a pore [START_REF] Saffman | A theory of dispersion in a porous medium[END_REF]. This implies that the persistence time of particles with small velocities may be much larger than suggested by τ c . Indeed, porescale velocity time series have been shown to display intermittent patterns, this means they are characterized by long periods of small velocities and rapid fluctuations of large amplitudes [START_REF] De Anna | Flow intermittency, dispersion, and correlated continuous time random walks in porous media[END_REF][START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF][START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF]Puyguiraud et al., 2019a), see also Figure 5.2. These patterns are indicative of a broad distribution of characteristic time scales. In fact, if the variance of τ is infinite, a sizeable amount of particles exhibits persistence times τ τ c . This means, particles do in general not become statistically equal on τ c , which invalidates the central assumption of the Brownian dynamics approach underlying Fickian dispersion. We will show below, that this property does not invalidate the existence of an REV.

Particle velocities vary on spatial scales imprinted in the medium structure rather than on a fixed times scale [START_REF] Kang | Pore-scale intermittent velocity structure underpinning anomalous transport through 3-d porous media[END_REF]Puyguiraud et al., 2019a). This property is accounted naturally by transport models in terms of continuous and time-domain random walks [START_REF] Berkowitz | Modeling non-Fickian transport in geological formations as a continuous time random walk[END_REF][START_REF] Painter | Upscaling discrete fracture network simulations: An alternative to continuum transport models[END_REF]. In fact, particle motion along CHAPTER 5. IS THERE AN REV FOR ANOMALOUS DISPERSION?

the mean flow direction in these frameworks can be described by

x n+1 = x n + c χ , t n+1 = t n + c v n , (5.12)
where, the particle velocities v n are independent random variables. For n = 0, they are distributed according to an initial velocity distribution p 0 (v). For n ≥ 1, they are independent identically distributed according to the flux-weighted Eulerian velocity distribution [START_REF] Dentz | Continuous time random walks for the evolution of Lagrangian velocities[END_REF])

p s (v) = vp e (v) v e .
(5.13)

This means two things. First, at each random walk step particles can sample the full velocity spectrum meaning that they are statistically equal. Second, the Lagrangian velocity statistics are stationary and ergodic. Particularly, they evolve toward their stationary steady state distribution on the length scale c . For these conditions to hold, it is necessary that particles can sample on the support scale a representative part of the Eulerian velocity distribution. Thus, in the following, we define criteria for the existence of a velocity REV in terms of the convergence of the velocity statistics with increasing support scale. Furthermore, we discuss the issue of ergodicity and stationarity.

Representative Elementary Volume

As discussed in the previous section, the representativeness of the velocity statistics sampled in the support volume and the existence of a stationary velocity distribution are key properties for transport upscaling for both Fickian and non-Fickian dispersion. Thus, for the support scale to be a transport REV, the velocity statistics need to be representative of the (stationary) Eulerian velocity statistics in the medium. To exhibit a representative velocity PDF, a sample must first be an REV for porosity because the Eulerian velocity distribution is linked to the pore size distribution [START_REF] Bibliography De Anna | Prediction of velocity distribution from pore structure in simple porous media[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF] and thus, an evolving porosity would translate into evolving velocity statistics. If the sample is an REV for porosity, then the Eulerian velocity statistics may be representative. We define an REV in terms of the Eulerian velocity PDF in a similar manner as the porosity REV. A sample is considered an Eulerian REV if it is large enough for the Eulerian velocity distribution to become stationary. To quantify the evolution of the Eulerian velocity PDF in function of the support scale, the Eulerian velocity PDF is sampled on growing domains starting from a small volume in the center of the sample to the full sample volume. The

DISPERSION UPSCALING AND THE REV

117 spatially sampled PDF of the magnitude of the flow velocity v e (x) = |v(x)| is defined as

p (v) = 1 φ V Ω dxδ [v -v e (x)] I(x ∈ Ω f ), (5.14)
where Ω is the physical domain on which the PDF is computed, and V and φ are its volume and its porosity, respectively. In order to quantify accurately the convergence of these distributions toward the full sample volume Eulerian PDF p e (v), we define the distance between p (v) and p e (v) based on the Kullback-Leibler (KL) divergence [START_REF] Kullback | On information and sufficiency[END_REF] as:

d KL (p e , p ) = ∞ 0 dvp (v) ln p (v) p e (v) .
(5.15)

The Kullback-Leibler divergence has been used to compare evolving PDFs to a reference distribution, (see for example [START_REF] Bigi | Using kullback-leibler distance for text categorization[END_REF][START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF][START_REF] Lindgren | Universality of probability density distributions in the overlap region in high Reynolds number turbulent boundary layers[END_REF]. When d KL = 0 it means that the distributions are identical (see Appendix, Section 8.3 for more details on the KL divergence). Here we consider a threshold value of ε = 10 -2

as the criterion for when the support volume can be considered an REV.

Lagrangian Ergodicity

As discussed above, convergence of the Eulerian velocity statistics on the support scale is not a sufficient condition for the upscaled random walk models discussed above to hold because these models also assume that particles can sample at each step from the same stationary Lagrangian velocity PDF, independently from the initial velocity distribution. In order to illustrate this, let us consider a porous media model consisting of a distribution of isolated straight capillaries. The support scale may be an REV for the Eulerian velocities.

However, since the flow velocities are constant along streamlines, particles are never able to sample the full velocity spectrum. The issue of Lagrangian ergodicity for porescale flow has been studied in detail in Puyguiraud et al. (2019a) in terms of the evolution of the s-Lagrangian velocity PDF, which is defined by

ps (v, s) = daρ(a)δ [v -v(s, a)], (5.16)
where ρ(a) is the initial particle distribution. We measure convergence of ps (v, s) toward the steady state p s (v) by the KL divergence d KL ( ps , p s ).

Berea Sandstone Sample

Here we illustrate the conditions under which anomalous dispersion can be described by CTRW approaches for a Berea sandstone sample. We first consider the REV definitions in terms of porosity and Eulerian velocity statistics. Then, we discuss the second condition, which refers to ergodicity and stationarity of the Lagrangian velocity series.

Representative Elementary Volume

In this section we study the concept of the REV in the light of the Eulerian velocity magnitude for the Berea sandstone sample illustrated in Figure 5.1 (see also, Puyguiraud et al., 2019a,b).

We first probe if the sample under consideration is a porosity REV. Figure 5.1 displays the evolution of the ratio φ between void and bulk volume computed on a size increasing domain that starts from a single point in the center of the sample to the total sample volume.

The porosity starts from a value of one since the initial volume is situated in the void space, it then quickly evolves towards the medium average porosity φ = 0.182 after about 4.5 p .

Since the sample is an REV for porosity, we now investigate the convergence of the Eulerian velocity statistics with increasing support scale. Figure 5.3a displays the Eulerian velocity PDFs respectively computed on cubes of volume V = 4 • 10 -8 V T , V = 1 • 10 -6 V T , V = 10 -5 V T , V = 0.125V T , and V T , where V T denotes the volume of the full sample. We observe that the distribution evolves toward the full sample distribution as the volume of the cube increases. For V = 4•10 -8 V T we only observe a small range of velocities. V = 1•10 -6 exhibits a distribution that recalls the velocity spectrum sampled in a single pore. The distribution sampled in V = 10 -5 V T corresponds to the average of several pore velocity distributions. Finally, despite being 8 times smaller, the volume V = 0.125V T seems to exhibit the same statistics as V T . To quantify accurately the convergence toward the Eulerian PDF, we use the aforementioned KL distance between the successive p (v) and p e (v). The distance between the full sample velocity PDF and the subsequent growing cubes' velocity PDFs is displayed in Figure 5.3b. For small volumes V the distance to the velocity PDF of full sample is large since the volume only contains a restricted range of the sample velocity spectrum, the distance then decreases quickly as the evolving distribution approaches the reference distribution. A distance d KL < 10 -2 is reached for a volume V ≈ 0.125V T . In other words, the limit distribution is attained. This indicates that the sample contains a Eulerian velocity REV. (Right) The KL distance (orange curve) between the steady-state Lagrangian velocity PDF and the successive Lagrangian velocity PDFs computed at increasing distances (from x = 0 to x = p ) along the mean flow direction.

DISPERSION UPSCALING AND THE REV

CHAPTER 5. IS THERE AN REV FOR ANOMALOUS DISPERSION?

Lagrangian Ergodicity

We have seen that an REV for the Eulerian velocity magnitude exists. This, however, is not sufficient for a CTRW based formulation to be valid for the upscaling of transport as discussed in Section 5.2.2.2. This means that the velocity statistics need to reach the steady state for distances of the order of the velocity correlation length, which sets another scale relevant for the definition of a transport REV. A representative part of the velocity spectrum can be sampled for transitions between REVs. We consider now the convergence from a given initial velocity distribution toward the steady state. To investigate accurately this evolution, we injected at the inlet particles in a set velocity range v ∈ [v l , v u ] and gave the same weight to the whole range (see blue solid curve in Figure 5.4a). The resulting distribution is far from the steady-state distribution. We display the spatial evolution of this distribution in Figure 5.4a. The distribution quickly evolves to the steady-state and is representative of p s (v) after x = 5 p .

We investigate quantitatively this convergence by computing the KL divergence between the steady-state velocity PDF p s (v) and the distributions p s (v, s) computed after different distances s along the particle streamlines (see, also Puyguiraud et al., 2019a).

Figure 5.4b displays the evolution of the KL distance. Since the distance is computed in the direction of the streamlines, we projected it onto the mean flow direction using the full sample tortuosity χ: x = s/χ, where χ is the tortuosity. At small distances the KL divergence between the evolving PDF and the steady-state is large because the injection PDF is very different from the steady state distribution, it then evolves quickly to reach a satisfying threshold d KL [p s (v, s), p s (v)] ≤ ε, where ε = 5 • 10 -2 after a distance x = 4 p in the mean flow direction. Despite the fluctuations that we observe at distance x > 4 p , we consider the convergence to be achieved since the KL distance remains below ε. These small fluctuations are due to the complexity of the geometry that particles encounter.

Since the sample fulfills the stationarity conditions that we require, it satisfies all the necessary criteria for the use of CTRW methods. It is characterized by the steady-state distribution (see, Figure 5.4a, red solid line). Note that stationarity of the particle velocity statistics is the necesssary condition for particle velocity series to be ergodic as discussed in Puyguiraud et al. (2019a).

IMPLICATIONS FOR DARCY SCALE TRANSPORT

Implications for Darcy Scale Transport

Purely Advective Transport

In this section, we compare the results obtained for the upscaling of the transport using an uncorrelated one-dimensional CTRW model with those obtained from the direct particle tracking simulations performed on the 3D binarized digital sample displayed in the inset of Figure 5.1. Details of the particle tracking methodology can be found in Puyguiraud et al. (2019a) and we will just summarize the main features of the calculations: the flow velocity field was computed by solving the Stokes equation using OpenFOAM, the streamline are obtained from a reconstruction algorithm based on local velocity interpolation [START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF], particles are injected at the sample inlet in a small volume displaying a velocity distribution p 0 (v). For the CTRW the initial particle velocity is picked in the initial velocity distribution p 0 (v) and the subsequent ones are picked in the steady-state p s (v), the particle motion and time then follow the set of equations (5.12). From this model definition, one can understand the need for a REV since the model relies on a steady-state distribution that parameterize all the particles. Without this REV definition, this type of CTRW would not be valid. For illustration we consider the dispersion coefficient D(t) in the direction of the mean flow, which is defined in terms of the displacement variance as following:

D(t) = 1 2 d dt x(t) 2 -x(t) 2 .
(5.17)

Figure 5.5 shows the temporal evolution of the dispersion coefficient in the mean flow direction for both uniform and flux weighted boundary conditions. The satisfactory match between the results of the 3D particle tracking simulations and the 1D CTRW model indicates firstly that the CTRW model is a powerful upscaled model and secondly validate the concept of the REV for anomalous dispersion. The model matches the theoretical late time scaling t/ln(t) 3 [START_REF] Comolli | Anomalous dispersion in correlated porous media: a coupled continuous time random walk approach[END_REF]. This model allows for an accurate computation of the transport properties based on the Eulerian velocity distribution, the tortuosity, and the correlation distance of the sample.

For a stationary initial velocity distribution p 0 (v) = p s (v), the CTRW framework gives for the bulk concentration p(x,t) the governing equation

∂ p(x,t) ∂t = t 0 dt K (t -t ) p(x -c /χ,t ) -p(x,t ) , (5.18)
where K (t) is defined through its Laplace transform by

K (x, λ ) = λ ψ(x, λ ) 1 -ψ(λ ) , (5.19)
where ψ(x, λ ) is the Laplace transform of the joint distribution of space and time ψ(x, λ ).

The transition time distribution is given by

ψ(t) = c t 2 p s ( c /t).
(5.20)

First, we notice that p(x,t) = φ c(x,t), where φ is porosity and c(x,t) is the concentration in the fluid phase only. Second Taylor expansion of the right side of (5.18) up to second order in c gives for c(x,t)

φ ∂ c(x,t) ∂t = - t 0 dt K (t -t ) c φ χ ∂ c(x,t ) ∂ x - 2 c φ 2χ 2 ∂ 2 c(x,t ) ∂ x 2 , (5.21)
which describes the evolution of the concentration, see also Appendix 8.2.2.

Advective-Diffusive Transport

The impact of diffusion in particle transitions concerns on one hand the spatial motion, on the other hand the temporal transitions. Here, we focus on the impact on the temporal for the computation of the transition times (see, for example, [START_REF] Russian | Time domain random walks for hydrodynamic transport in heterogeneous media[END_REF]. The definition of the transition times ψ(τ) can then be expressed as

ψ(τ) = P s (v)δ t - τ D τ v τ D + τ v dv.
(5.23)

In Appendix 8.2.1 we derive

ψ(τ) = c t 2 p s c τ D -t τ D t .
(5.24)

The dependence on τ D implies a sharp cut-off of the transition time for t > τ D because p s (v) = 0 for v < 0. The governing equation for the bulk particle density is given by (5.21),

where the kernel now is defined in terms of the transition time distribution (5.24). As ψ(t) has a cut-off for t τ D , all moments exist and Equation (5.21) can be localized in time for t τ D . Thus, we obtain for c(x,t) = p(x,t)/φ the evolution equation (5.25) where the dispersion coefficient D * is given by Dentz et al. ( 2004)

∂ c ∂t (x,t) = -q ∂ c ∂ x (x,t) + D * (τ D ) ∂ 2 c ∂ x 2 (x,t),
D * = 2 c τ τ 2 -τ 2 τ 2 .
(5.26)

Direct Consequences for Dispersion and the ADE Framework

We observed in Section 5.3.1 that in a fully advective scenario, the dispersion keeps growing and exhibits at late times D(t) ≈ t ln 3 (t) scaling. Such observations put in evidence the non-existence of a time scale τ D after which the dispersion stabilizes to an asymptotic value D ∞ . Consequently, this proves the lack of a scale at which the ADE (Eq. (5.3)) becomes valid. Nevertheless, these results were the ones expected due to the lack of diffusion in our simulation, a diffusion that usually brings a cutoff in the transition time distribution. Therefore, in the absence of such cutoff, the dispersion keeps growing. However, natural Péclet numbers at the pore-scale have been shown to vary between 10 -2 and 10 6 (see [START_REF] Bear | Dynamics of fluids in porous media[END_REF][START_REF] Bijeljic | Pore-scale modeling and continuous time random walk analysis of dispersion in porous media[END_REF] evidencing that while there exists a broad range of high Péclet numbers, the diffusion cannot usually be neglected. In the presence of diffusion, a cutoff appears in the transition time distribution. The lower the Péclet number and the earlier this cutoff occurs. This translates directly into a shorter time τ D after which the dispersion converges. To quantify this convergence, we investigate the temporal evolution of the dispersion in function of the Péclet number. To be able to vary the Péclet number of our transport simulation, we need to insert a diffusion process in our particle tracking. Diffusion has often been added to particle tracking algorithms in the form of a random walk, (see for example, [START_REF] Ahlstrom | Multicomponent mass transport model: theory and numerical implementation (discrete-parcel-random-walk version)[END_REF][START_REF] Ackerer | Random-walk method to simulate pollutant transport in alluvial aquifers or fractured rocks[END_REF][START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF]. We use a similar methodology, full details can be found in Section 2.2.1.2. sampled. Then, the dispersion becomes constant. Note that the lower the Péclet number, the earlier the convergence to the asymptotic dispersion D ∞ occurs and the lower its value.

The opposite rises for larger Péclet number cases.

This set of results allows us to investigate the evolution of the dispersion in function of the Péclet number (see Figure 5.7). We observe similar results to the one observed in the literature: the dispersion (undimensionalized as D D 0 by the diffusion coefficient) is constant for small Péclet number (P e < 0.1) since the molecular diffusion is the only mechanism of mixing. The diffusion is limited as the solid grains of the medium act as barriers for particles. This results in a dispersion coefficient lower than 1. Note that at very early time the dispersion goes like D 0 while particle movements have not been impacted yet by the geometry [START_REF] Bijeljic | Pore-scale modeling of longitudinal dispersion[END_REF]. This restricted diffusion is for example described by the reciprocal value of the product of formation factor and porosity (see, for example, [START_REF] Brigham | Experiments on mixing during miscible displacement in porous media[END_REF][START_REF] Sahimi | Flow and transport in porous media and fractured rock: from classical methods to modern approaches[END_REF][START_REF] Bijeljic | Pore-scale modeling of longitudinal dispersion[END_REF] or similarly the relation between D and D 0 can be written as

D = D 0 φ m (5.27)
where φ is the porosity and m is the cementation factor (also called porosity exponent or cementation exponent, (see, for example, [START_REF] Kadhim | A review in correlation between cementation factor and carbonate rock properties[END_REF]). We obtain m ≈ 2.1 (see Figure 5.7) which in the range of value obtained by measurements in the literature. The cementation factor has been reported to vary between 1.3 and 2.5 for most rocks and to be around 2 for sandstone formation [START_REF] Attia | Irreducible water saturation from capillary pressure and electrical resistivity measurements[END_REF] while being around 5 for carbonates.

A set of reported results on a few Berea sandstone cores was 1.79 -1.92 [START_REF] Attia | Effects of petrophysical rock properties on tortuosity factor[END_REF].

The effect of advection on dispersion starts to be visible for P e > 0.1, it then progressively transitions from a slope of D ∞ /D 0 ≈ P 0 e to a slope of D ∞ /D 0 ≈ Pe 2 / ln(Pe) 3 from P e = 0.1 to P e = 20 when advection becomes the main mechanism that impacts dispersion. This asymptotic dispersion coefficient for finite Péclet can be estimated from the evolution of the transient D(t) for infinite Péclet, which evolves here in time as D(t) ≈ t/ ln(t) 3 . The cut-off on the diffusion time stabilizes the dispersion coefficient at time t = τ D (see Figure 5.6), which implies that D ∞ can be approximated by D ∞ ≈ Pe/ ln(Pe) 3 . This means that we expect a dependence of D ∞ /D 0 ≈ Pe 2 / ln(Pe) 3 . The data illustrated in Figure 5.7 confirm this estimate. To understand on what spatial scales the dispersion stabilizes and the transport can become Fickian we investigate the average distance particles need to travel to reach the asymptotic dispersion value. Figure 5.8 displays the average number of pores the particles have to travel before reaching a constant dispersion. At very low Péclet numbers (Pe < 0.1), the asymptotic diffusion coefficient is reached after about half of a pore length. This means that since the diffusion is the leading mechanism, it does not matter which particular throat the particles are experiencing. The small advective motion is not visible in the values of the dispersion. At Pe = 1, the advection starts to have an impact and the particles need to sample a distance of about one pore length. The interesting results happen for intermediate Péclet (10 < Pe < 100), where the particles need to travel between 100 and 5.4. CONCLUSIONS 127 1000 pores for the dispersion to finally become constant. This means that the particles need to sample a large part of the heterogeneity before their spread can become constant in time.

At larger Péclet numbers, the diffusion is very small compared to the advection, this leads to a very large time before a particle be able to sample the full velocity spectrum of the domain, this leads to a very high number of pores traveled before reaching the asymptotic value. Note that the number of pores seems to evolve linearly for Pe > 100 (linear scaling can be see in the black solid line of Figure 5.8). These numbers give insight for the spatial scales at which the dispersion might converge and therefore on what scale the ADE might become valid in function of the Péclet number. Note also that we assume here that the rock geometry exhibits the same pore distribution over distances larger than a thousand pores. This is of course not true in nature and might lead to even larger distances before the full heterogeneity can be sampled.

Conclusions

In this paper we aimed at answering the following question: Is there an REV for anomalous dispersion? We showed that not only this REV definition exists but that it is also a necessary condition for CTRW and TDRW approaches to work. This REV definition needs to be set in terms of the Eulerian velocity statistics in order to allow the upscaling of anomalous dispersion through CTRW type models. Precisely, the sample is considered a Eulerian velocity REV if it is large enough for the Eulerian velocity PDF to be constant when increasing the size of the sample domain. Note that this definition requires that the sample also be a porosity REV since an evolving porosity would prevent the convergence of the Eulerian velocity distribution. The REV is then parameterized by this Eulerian velocity PDF (similarly to a porosity REV being characterized by its mean porosity). However, it is not sufficient, for CTRW and TDRW to work, Lagrangian velocity series are required to be stationary processes that relax toward their stationary distribution within a characteristic length scale. This means that the samples also needs to fulfill the following criteria: first, a steadystate velocity distribution for the sample needs to exist since an evolving distribution would prevent the characterization by a single distribution. This means that observing that a Lagrangian velocity distribution converges to a steady-state PDF constitutes a sufficient condition for the sample to be a Lagrangian v-REV. Also, since the velocity increment process of a CTRW type models is ergodic by construction, the velocity process inside an REV needs to be ergodic as well. We then show that, when combined, these criteria allow for an accurate upscaling via the CTRW framework. The CTRW framework that we described 128 CHAPTER 5. IS THERE AN REV FOR ANOMALOUS DISPERSION? requires the existence of steady velocity conditions. Furthermore, the CTRW requires an ergodicity criterion meaning that an ensemble of particle is able to sample the full velocity spectrum. This specification is guaranteed by the Lagrangian v-REV stationary condition.

To illustrate and experiment these definitions and their consequences, we selected a

Berea sandstone sample as a potential REV candidate. The sample, which is a REV for porosity, happens to fulfill the Eulerian REV definition since the velocity distribution converges within 1/8 of the total sample size. We also observed that a randomly elected particle distribution at the inlet leads to a stationary velocity distribution after a distance of about two thirds of the sample, guaranteeing that the sample is also a Lagrangian v-REV.

After verifying that the ergodicity assumption is fulfilled, we upscaled the hydrodynamic transport process with a one-dimensional CTRW model based on these v-REV properties.

It relies on the steady-state velocity PDF for its parameterization and captures well the temporal evolution of the dispersion (and other transport properties (Puyguiraud et al., 2019a))

in a purely advective case. This confirms that the velocity distribution contains all the necessary information for the prediction of the transport.

We then aimed at investigating the connections this methodology brings to Darcy scale transport through the CTRW framework. By definition, the ADE validity is restricted to asymptotic dispersion since it cannot account for the velocity fluctuations impact under that scale. In contrast, CTRW models do not depend on characteristic scales and can therefore represent reliably scale dependent coefficients such as dispersivity. CTRW models have been used in the past to capture Darcy scale transport evolution with uncorrelated jumps (see Section 5.2.2.2). We showed that this type of models also rely on REVs for their parameterization and that they can predict transport properties such as dispersion. However, while this methodology is a powerful tool for the prediction of anomalous transport, the methodology is still limited by several factors. The first one is that it relies on an accurate knowledge of the velocity distribution inside the considered geometry. While large progresses have been made toward relating the Eulerian velocity distribution to the pore geometry [START_REF] Bibliography De Anna | Prediction of velocity distribution from pore structure in simple porous media[END_REF][START_REF] Alim | Local pore size correlations determine flow distributions in porous media[END_REF][START_REF] Dentz | Mechanisms of dispersion in a porous medium[END_REF], it still remains an open issue. In addition, obtaining an accurate description of the pore space at intermediate to large scales is also challenging. Another issue that rises is to account for the diffusion process that cannot in general be neglected. While the velocity distribution contains all the necessary information for the prediction of the advective transport, diffusion forces mass exchanges between the streamlines at low and intermediate values of the Péclet number. It follows that the characterization of the hydrodynamic transport process is more complex.

The advective velocity transition is governed by the a spatial scale engraved in the pore structure while the diffusive transition happens on a temporal scale. We showed that this can be accounted for by introducing a cutoff in the transition time distribution of the CTRW model. This is similar to the advective-diffusive coupling of time domain random walk models. This allows for the modeling of transport for any Péclet case and for large times a evolution equation for the concentration can be recovered.

Any upscaling model relies on stationary properties that are averaged over an REV. We showed that defining the REV in terms of the velocity statistics allow for the upscaling of anomalous dispersion with CTRW type models. With the absence of REV, any tentative of upscaling through this type of model would lead to erroneous results since the parameter inserted in the model would not be representative of the studied medium.

INTRODUCTION

135 as a consequence the reaction rates can be significantly less than the ones measured under well-mixed conditions [START_REF] De Anna | Mixing and reaction kinetics in porous media: An experimental pore scale quantification[END_REF]Dentz et al., 2011a;[START_REF] Tartakovsky | Effects of incomplete mixing on multicomponent reactive transport[END_REF][START_REF] Willingham | Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments[END_REF][START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF]. Understanding the role of mixing in reactions and determining an appropriate model of mixing is essential for the correct characterization and upscaling of chemical reactions in heterogeneous porous media.

Classical reactive transport models assume complete mixture of reactants at the scale of the average volume. The advection-dispersion-reaction equation (ADRE) defined as

φ ∂ c i (x,t) ∂t = -∇ • [qc i (x,t) -D∇c i (x,t)] -r i , (6.1)
where φ is porosity, c i is the concentration of reactant i, q is the Darcy velocity, D is the dispersion tensor, and r i represents the space-time-dependent rate at which a species i is produced (or eliminated) by the reaction, relies on the assumption that reactants are well-mixed over the representative elementary volume. Chemical models based on this assumption tend to overpredict the reactive behavior or the mixing front in the system [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore scale reactive transport modeling on images of porous media[END_REF][START_REF] Ding | Modeling bimolecular reactions and transport in porous media via particle tracking[END_REF][START_REF] Zhang | Evaluation and linking of effective parameters in particle-based models and continuum models for mixing-limited bimolecular reactions[END_REF]. Such chemical overprediction occurs because the ADRE reaction rate is commonly estimated from batch tests under perfect mixing conditions [START_REF] Berkowitz | Measurements and models of reactive transport in geological media[END_REF]Dentz et al., 2011b). A series of studies has focused on proving the validity of the ADRE to predict mixing-limited chemical reactions by performing averages [START_REF] Porta | Upscaling solute transport in porous media in the presence of an irreversible bimolecular reaction[END_REF], or using fitting parameters calibrated from experimental data [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF]. A series of laboratory experiments [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] De Anna | Mixing and reaction kinetics in porous media: An experimental pore scale quantification[END_REF][START_REF] Willingham | Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments[END_REF][START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF] and field studies [START_REF] Hess | Multispecies reactive tracer test in an aquifer with spatially variable chemical conditions, Cape Cod, Massachusetts: Dispersive transport of bromide and nickel[END_REF][START_REF] Davis | Multispecies reactive tracer test in an aquifer with spatially variable chemical conditions[END_REF] has shown that the ADRE overestimates the quantity of reaction occurring by pointing out that chemical reactions occur locally at the pore scale, where the ADRE ignores incomplete mixing of reactants.

Recent alternative mathematical models have related the effective reactivity to transport at local scales where the incomplete mixing of the solutes can be taken into account, this includes particle-based Lagrangian models [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF][START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore scale reactive transport modeling on images of porous media[END_REF][START_REF] Zhang | Evaluation and linking of effective parameters in particle-based models and continuum models for mixing-limited bimolecular reactions[END_REF][START_REF] Ding | Modeling bimolecular reactions and transport in porous media via particle tracking[END_REF][START_REF] Edery | Modeling bimolecular reactions and transport in porous media[END_REF][START_REF] Benson | Simulation of chemical reaction via particle tracking: Diffusion-limited versus thermodynamic rate-limited regimes[END_REF] and interface deformation models (de Anna et al., 2014a;[START_REF] Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF]. Many of these methods consist in predicting and quantifying the evolution of the interface between reacting species. At early times the interface between the two reactants can be described as an aggregate of complex elongated independent structures called lamellae or diffusive strips. The lamellar approach [START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Ranz | Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows[END_REF] has been the basis of predictions for mixing and chemical reactions in a variety of flow situations ranging from chaotic CHAPTER 6. PREDICTION OF PORE-SCALE REACTION

The size of the cells is similar to the first medium cells. To generate this medium, we allowed for a larger distance between grains which results in a wider range of pore throats.

It therefore differs from the former by its larger heterogeneity that forces the creation of preferential channels. This results in a more heterogeneous flow field (see bottom panel of Figure 6.1). Note that, in order to minimize boundary effects, twenty layers were added at the inlet and outlet of both the geometries.

Flow Simulation

In the following we summarize the methodology to solve the flow field (full details for a 3D case are available in Puyguiraud et al. (2019a)). At the pore scale, under low Reynolds number conditions (Re 1), the flow is governed by the Stokes equation [START_REF] Leal | Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes[END_REF]:

∇ 2 v(x) = 1 µ ∇p(x), (6.2)
where v is the velocity vector, p(x) is the pressure and µ the kinematic viscosity of the fluid. The equation is solved together with the continuity equation:

∇ • v(x) = 0, (6.3)
that guarantees the incompressibility of the fluid. We prescribed pressure boundary conditions at the inlet and outlet, and no-slip conditions at the void-solid interfaces and at the remaining domain boundaries. We then solved the flow with the SIMPLE algorithm [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF] implemented in OpenFOAM. After convergence, we extracted the complete velocity field. This resulted in the mean velocity values being expressed at every interface of the mesh in the normal direction to the face. Figure 6.1 displays the flow fields inside the two domains.

Reactive Transport

In this section we discuss the implementation of the advective displacement, diffusive motion, and reactive step of our particle tracking algorithm. The operators are split in the sense that we perform sequentially the three steps. The position is first updated as (6.4) where x(t) is the vector position of a particle at time t, ∆x a (∆t) and ∆x d (∆t) are the advective and diffusive displacements, respectively. We then proceed to the reactive step.

x(t + ∆t) = x(t) + ∆x a (∆t) + ∆x d (∆t),

Advective Displacement

We use an extension of the Pollock algorithm [START_REF] Pollock | Semianalytical computation of path lines for finite-difference models[END_REF][START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF][START_REF] Puyguiraud | Upscaling of anomalous pore-scale dispersion[END_REF] to account for the no-slip boundary conditions. This allows for an accurate interpolation of pore-scale porous medium velocity. Integrating these velocity equations gives rise to analytical trajectories expressed as functions of time in every cell. Thus, every particle position after any given time can be efficiently and accurately determined. The current location of a particle is then updated as x(t + ∆t) = x(t) + ∆x a (∆t), (6.5)

y(t + ∆t) = y(t) + ∆y a (∆t), (6.6)
where ∆x a (∆t) and ∆y a (∆t) are the x-and y-advective displacements of the particle during time ∆t respectively. During the simulation, the position of every particle is updated accordingly before proceeding to the diffusive motion.

Diffusive Motion

The diffusion process in particle tracking algorithms is often implemented through random walk methods [START_REF] Ahlstrom | Multicomponent mass transport model: theory and numerical implementation (discrete-parcel-random-walk version)[END_REF][START_REF] Ackerer | Random-walk method to simulate pollutant transport in alluvial aquifers or fractured rocks[END_REF][START_REF] Mostaghimi | Simulation of flow and dispersion on pore-space images[END_REF]. The CHAPTER 6. PREDICTION OF PORE-SCALE REACTION equivalence of the Eulerian diffusion equation to the Lagrangian Langevin equation for an ensemble of particles [START_REF] Risken | The Fokker-Planck Equation[END_REF] dx dt = √ 2Dξ (t), (6.7)

with D the diffusion coefficient and ξ a Gaussian white noise, allows for the definition of the diffusive displacements occurring in the xand y-direction during a time step ∆t as

∆x d (∆t) = √ 2D∆tξ 1 , (6.8) ∆y d (∆t) = √ 2D∆tξ 2 , (6.9)
where ξ 1 , and ξ 2 are Gaussian white noises. We follow a similar computation since our diffusive jumps are computed as

∆x d (∆t) = √ 3 √ 2∆tDγ 1 = √ 6∆tDγ 1 , (6.10) ∆y d (∆t) = √ 3 √ 2∆tDγ 2 = √ 6∆tDγ 2 , (6.11)
where γ 1 and γ 2 → U (-1, 1) (uniformly distributed between -1 and 1). The Central Limit Theorem insures that the sum of the jumps converges toward a Gaussian distribution.

Multiplying by √ 3 gives a variance of 1. This scheme allows for a better control on the maximum diffusive movement and avoid the costly numerical generation of Gaussian random numbers. We implemented reflective boundaries at the solid interface and absorbing boundaries at inlet and outlet of the domain.

Reaction Process

Reaction is simulated following the methodology presented in Perez et al. (2019a). At each time step, we record the position of each particle as it migrates through the domain and calculate the distance between a given A particle and a B particle. The probability of reaction P r of the B particle in the time interval [t,t + ∆t], depends on the number N A [x(t)]

of A particles within the interaction well-mixed support volume ∆V centered at the position x(t) of the B particle as 6.12) where the probability of a single reaction event p(∆t) = k∆t/(N 0A ∆V ) depends on the chemistry of the problem characterized by the reaction rate coefficient k, and the initial total number of A particles present in the domain N 0A . The interaction well-mixed support volume ∆V = πr 2 is defined with an effective reaction radius r = √ 24D∆t. The well-mixed conditions mean that all reactant particles within the interaction radius, or support volume,

P r = 1 -exp [-p(∆t)N A [x(t)]] , (

METHODOLOGY

141 have the same probability of reaction in a time interval ∆t. The selection of r relates to the characteristic diffusive particle displacement during time ∆t, which is σ (∆t) = 2dD∆t, where d is the spatial dimension. For r ≤ σ (∆t) the support volume may be considered well mixed if N 0 → ∞. Otherwise, r needs to be larger than σ (∆t) in order to capture the local degree of mixing at which reaction takes place. The limits and criteria for the choice of the reaction radius can be found in Perez et al. (2019a). The reaction occurs if P r is larger than a random number picked uniformly between 0 and 1. After reacting, the A and B particles are removed from the system and a particle C is injected in the middle of the A and B particle locations. The migration of C particles in the domain also follows the transport rules specified in Equation (6.4). The total molar concentration of A can be expressed from the number of particles as (6.13) where m 0 is the initial number of moles that the species A carries, N A [x(t)] is the number of A particles at location x and time t, N 0A is the initial total number of A particles present in the domain and ∆V is the support volume. Then, the total mass of C particles is calculated

m A (t) = m 0 N A (x,t) N 0A ∆V dx,
as m C (t) = m 0 N C (x,t) N 0A ∆V dx, (6.14)
where N C (x,t) is the number of C particles at location x and time t.

Simulation Setup

Initially, the A species is placed uniformly throughout the pore space between the inlet plane x = 0.25mm and x = 4.25mm which corresponds to about two thirds of the domain, and there is no B species in the medium. We use this initial spatial distribution of species A to numerically simulate a medium that is already filled with a species.

Note that placing reactant A everywhere in the domain increases the simulation computational cost because the transport and reaction equations will be applied to particles that are near to the outlet of the domain. Particles close to the domain's exit do not react since they leave the medium after a very short time. To simulate a continuous injection the B particles keep being injected at location x = 0.25mm at all time t > 0. We found that injecting the B particles until the end of the domain slows down the computational simulation while not impacting the reaction since the particles injected in the last part of the domain never react during our simulations. This setup results in the A and B particles interface being located at position 2.5 × 10 -4 m at t = 0, which corresponds to the second column of grains in the geometries. Particles then follow sequentially the advective, diffusive and reactive steps until the end of the simulation. The time step is picked sufficiently CHAPTER 6. PREDICTION OF PORE-SCALE REACTION small to have good accuracy and not to affect the results (see, [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF]. The simulations were performed using number of particles N p = 6.5 × 10 6 with a mean ve-

6.3

The Dispersive Lamella Description of Mixing

Bases

In the Fickian approach the reaction support volume is assumed to be well-mixed. Then, in this framework, the evolution of the concentrations c i can be described by the advectiondispersion-reaction equation (6.1). The global reaction behavior can be characterized by the evolution of the total mass of the reaction product

m C (t) = c C (x,t)dx. (6.15)
For an instantaneous bimolecular reaction in a nearly homogeneous porous medium, [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF] characterized the evolution for the total mass of C from (6.1) as,

m C (t) = c 0 L y φ 4D h t π , (6.16)
where c 0 is a characteristic concentration, and D h is the hydrodynamic dispersion coefficient that describes the spreading of solutes. It is defined as .17) where σ 2 a (t) is the apparent variance which is a measure for the dispersion of the interface and is defined as

D h = lim t→∞ σ 2 a 2t . ( 6
σ 2 a (t) = dx c 0 (x ) dx[x -m x (t)] 2 g(x,t | y ), (6.18) m x (t) = dx c 0 (x )m x (t | y ). (6.19)
The √ t scaling of the evolution of the product C mass in Eq. (6.16) can be described in terms of the reaction rate, which is equal to the diffusive mass flux at the interface between the two reactants. The mass obtained from (6.16) serves as a reference for observed behaviors in spatially variable flows. for Pe = 60 from the RPT simulation in the first porous medium (symbols), and from the hydrodynamic dispersion coefficient (black solid line). The analytical solution overpredicts the product total mass in the medium due to the assumption of complete mixing between reactants. Note that similar behaviors have been observed by [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF]. Figure 6.3 compares the total mass of product calculated from the numerical simulations to the analytical prediction (6.16) for the first medium. The analytical solution overpredicts the product mass with respect to the numerical simulation by ∼ 27%. This finding agrees with the observations made in [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF], where the authors measured the outflow concentrations of the product formed during a mixing-limited reaction. This means that the chemical reaction cannot be well quantified by a the reactive transport description based on the hydrodynamic dispersion. This is because the plume has not sampled yet the full heterogeneity of the velocity field and thus, does not spread linearly with the hydrodynamic dispersion yet. Therefore the amount of mixing experienced by the plume is overestimated when using D h . In the next section we elaborate on the dispersive lamella methodology for the measuring of the plume spreading.

Quantification of Mixing

In this section, we present the concept of the dispersive lamella used to quantify the impact of fluid mixing on chemical reactions. This representation of mixing assumes that solutes tend to organize into structures that are formed by the repeated action of advection as they move through heterogeneous media. These structures, called lamellae, can be seen as a decomposition of the mixing front into point elements that disperse as a result of the diffusive sampling of the vertical velocity contrast. This concept of not independent lamellae, called dispersive lamellae, differs from the stretched lamella model [START_REF] Bandopadhyay | Enhanced reaction kinetics and reactive mixing scale dynamics in mixing fronts under shear flow for arbitrary Damköhler numbers[END_REF][START_REF] Borgne | Impact of fluid deformation on mixinginduced chemical reactions in heterogeneous flows[END_REF][START_REF] Le Borgne | Stretching, coalescence, and mixing in porous media[END_REF][START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Meunier | The diffusive strip method for scalar mixing in two dimensions[END_REF][START_REF] Villermaux | Mixing as an aggregation process[END_REF][START_REF] Ranz | Applications of a stretch model to mixing, diffusion, and reaction in laminar and turbulent flows[END_REF] because of the interacting lamellae. The dispersive lamella approach is based on the concept of effective dispersion and accounts for the action of transverse diffusion in contrast to the stretched lamella model.

Let us start from a pulse line injection at the inlet of the medium. It is composed of simultaneous point injections, with initial conditions distribution c 0 (x). The solute is initial distributed along a line perpendicular to the mean flow direction,

c(x,t = 0) = c 0 (x) = 1 L y δ (x 1 ). (6.20)
The concentration c(x,t) satisfies the advection-diffusion equation The Green function satisfies (6.21) for the initial condition g(x,t = 0|y ) = δ (x)δ (yy ).

∂ c(x,t) ∂t + v(x)∇c(x,t) -D∇ 2 c(x,t) = 0, ( 6
We transform into the coordinate system that moves with the center of mass of the Green

function x = x -m(t | y ), (
( x,t | y ) = exp -( x-x ) 2 2σ 2 y (t) 2πσ 2 x (t) , (6.29) G( ŷ,t | y ) = exp -( ŷ-y ) 2 2σ 2 y (t)
2πσ 2 y (t) , (6.30)

where σ 2 x (t) and σ 2 y (t) are the effective spatial variances. There are defined by Figures 6.6 and 6.7 show the temporal evolution of apparent variance σ 2 a (t) (defined in Eq (6.18)) and effective variance σ 2 e (t) for the two different cases, respectively. The calculated σ 2 a (t) and σ 2 e (t) are similar at early times (t < 0.26τ v and t < 0.35τ v , respectively). In this temporal regime, we find that the behavior in both variances is similar to 2Dt, which suggests a diffusion dominated regime. This observation is reflected in the top-left inset of Figure 6.6 which shows a nearly homogeneous front from the spatial distribution of particles at t < 0.26τ v . The snapshot suggests that the front has not been affected by the advective heterogeneity yet. For later times (0.26τ v < t < 100τ v and 0.35τ v < t < 100τ v , respectively), the apparent variance σ 2 a (t) grows faster than the effective variance σ 2 e (t) because the plume experiences the velocity contrast from the advective field and is therefore deformed. The advective deformation, or spreading, is responsible for the rapid increase of σ 2 a (t) over σ 2 e (t). Note that the impact of advection can be observed in the top-right inset of

σ 2 x (t) = dx dx[x -m x (t | y )] 2 g(x,t | y )c 0 (x ), (6.31) σ 2 y (t) = dx dx[y -m y (t | y )] 2 g(x,t | y )c 0 (x ). ( 6

Reaction Behavior at Pore-Scale

As suggested by our results in Figure 6.3 and the experimental observations in [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF], the reactive transport description based on the hydrodynamic dispersion coefficient does not quantify properly the chemical reaction. For the initial condition c 0 , the concentration of C across the lamella is given by (6.33) and the concentration of C at the interface is This expression accounts for the impact of the interface deformation on the overall reactivity since L y φ σ e (t) gives the area of the mixing zone. Figures 6.8 and 6.9 show the evolution of the C product mass obtained from the RPT model and the estimate (6.36) from the dispersive lamella approach for the first and second medium, respectively. For comparison we also display the evolution of the product mass for an interface that would evolve like the apparent width σ 2 a (t). At early times, (t < 0.26τ v and t < 0.35τ v , respectively) diffusion is the main mechanism driving the reaction, as a result mass predictions using σ 2 a (t) and σ 2 e (t) are similar and agree with the solutions estimated from Equation (6.16) coupled with a constant diffusion coefficient D. At later times (0.26τ v < t < 100τ v and 0.35τ v < t < 100τ v ), advection dominates. We find enhanced-mixing behaviors as the product formation increases rapidly. This occurs due to a greater degree of mixing of reactants locally (illustrated in Figure 6.2). In this regime, the dispersive lamella approach CHAPTER 6. PREDICTION OF PORE-SCALE REACTION coupled with σ 2 e (t) provides an accurate description of the C product mass evolution. The increased reaction behavior that we observed occurs preferentially as the reactant particles sample more of the flow heterogeneity, which increases the width of the interface allowing for more mixing. We find that the mass prediction based on σ 2 a (t) fails to characterize the evolution of the product mass because it overestimates the reactants interface as discussed earlier. The greater degree of heterogeneity of the second medium causes the flow to be more heterogeneous and therefore enhances the mixing. This translates into a larger overall width at late times. Note that it is well known that mixing is a limiting factor for instantaneous reactions [START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF], this behavior is well captured by our model. The particle number for the RPT needs however to be large enough to avoid "numerical incomplete mixing" which prevents particles from reacting because of a too low density. However, results converge with the number of particles. This methodology is a powerful tool since it allows for an accurate computation of the C mass produced while not having to perform the reaction. Moreover, the computation of the effective width does not demand a large amount of particles and therefore requires less computational power. Note that we used the porosity of the total medium for the computation of the mass (see Equation (6.36)), however the results might suffer from discrepancies if the porosity experienced by the plume of particles differs from the actual medium porosity. This is discussed in the next section. Mass predicted using σ 2 a (t) Mass predicted using σ 2 e (t) RPT model Figure 6.9: Evolution of the mass m C (t) for Pe = 40 from the RPT simulation (symbols), from the dispersive lamella parameterized by the apparent variance σ 2 a (t) (blue dashedline) and the effective variance σ 2 e (t) (red line) for the second medium.

θ C ( x,t) = L y c 0 2 | x | 2σ 2 e (t) ,
c C ( x,t) = 1 L y L y 0 dy G( ŷ,t | x )θ C ( x,t) = c 0 2 erfc | x | 2σ 2 e (t

Porosity Experienced by the Particles

In this section we discuss the impact that the porosity experienced by the particles has on the reaction mass produced. Our predictive approach uses the porosity of the medium to correct the mixing volume computed. However, there exist situations where the plume experiences a porosity that differs from the average porosity of the medium. In this eventuality, this method leads to an overprediction of the amount of mixing in the case of the porosity experienced by the plume being smaller that the total porosity or to an underprediction if it is larger.

To investigate this potential effect, we numerically mimic a situation that can happen in the laboratory: we place the injection out of the second medium at a location where the initial plume experiences a porosity φ = 1. We use the setup discussed in Section 6.2.3.4 with the difference that the interface between the A and B particles at t = 0 is originally located at x = 5 × 10 -5 m, which corresponds to the middle of the first empty layers in the geometry. We display in Figure 6.10 the evolving porosity φ (t) experienced by the plume: at early times the interface experiences a very high porosity, then the plume porosity starts decreasing quickly once the first particles reach the first grains, finally the experienced porosity converges to the medium average porosity 0.53. Note that due to the absence of grain at the inlet, the porosity seen by the particles is a little larger than the real porosity of from the first medium porosity to the second medium porosity. Therefore, an estimate of the evolving porosity in time might be handy in these situations.

Summary and Conclusions

In this paper we simulated a bimolecular reaction transport problem in two dimensions: a B reactant injected at the inlet reacts while a A reactant occupies the rest of the complex medium. We performed the simulation using a reactive particle tracking algorithm that performs sequentially the advective, diffusive, and reactive steps. These simulations allow for the computation of the mass of C produced and are used to validate the dispersive lamella methodology proposed by [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF] in heterogeneous media in the presence of very heterogeneous velocity fields. The original dispersive lamella model predicts the amount of C mass produced in function of the measured mixing volume of a line injection and was validated in a tube against reactive random walk simulations [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF]. We validated these predictions against RPT simulations for heterogeneous bead pack media. The predictions captures the reaction produced at all times. Moreover, this methodology can easily be generalized to 3D since the measurement of the mixing front is CHAPTER 6. PREDICTION OF PORE-SCALE REACTION also straightforward in three dimensions. This method gives powerful tools for the estimate of the reaction mass produced in heterogeneous porous media since it allows for its computation based only on the mixing volume and thus, no reaction computation is required. This reaction computation often slows down simulations and therefore limits the maximum computational domain size. This method thus allows for larger computational domains. Therefore, this capability of simulating larger domains can give bases for the modeling of reactive experiments whose geometries are usually larger than the numerically simulated domains. The methodology can also be utilized for the prediction of the chemical mass produced from conservative transport experiment through the observation of the mixing volume. Note, however, that if the porosity experienced by the plume is not representative of the average porosity of the medium the method might lead to inaccuracies. A solution is to consider the porosity experienced by the plume instead of the full medium porosity for the computation of the reaction mass.

Note also that since the prediction of the C mass only relies on an approximation of the effective variance of the front, any upscaling model able to predict this information would allow for a prediction of the reaction. CTRW models [START_REF] Puyguiraud | Upscaling of anomalous pore-scale dispersion[END_REF] and MRMT models [START_REF] Li | Transport upscaling using multirate mass transfer in three-dimensional highly heterogeneous porous media[END_REF] for example, are able to give accurate transport prediction for example in terms of displacement variance. Modifying these models to compute the effective σ e accurately would provide powerful reactive transport upscaling tools. This instantaneous bimolecular reaction A + B → C is a simple reaction and is not representative of most chemical reaction happening in the nature. However, any real complex reaction can be split into many simple instantaneous bimolecular reactions of the type A + B → C that can be treated one by one. Therefore this methodology is applicable to any fluid-fluid chemical reaction. The next step for this methodology might be to address a different type of chemical studies such as fluid-solid reactions.

CHAPTER 7

Conclusions

Chapter 7

Conclusions

General Conclusions

In this thesis, we have provided a comprehensive analysis of the stochastic dynamics of pore-scale velocities experienced by solute particles. We quantified the evolution of Lagrangian and Eulerian statistics and showed that the injection conditions drastically impact the velocity statistics experienced by solute particles. We also derived the relation between the isochronously and equidistantly sampled velocities. We found that, because of their fast convergence and pore-architecture-nested properties, studying spatial velocity statistics might be more handy than investigating their temporal counterparts. These observations lead us to develop processes for the prediction of the spatial velocity evolution. These processes, based on Markovian properties, capture faithfully the full spatial Lagrangian velocity distribution evolution. Knowing that these Lagrangian velocity statistics directly control the displacement of the solute particles, we aimed at relating them to the transport. To do so, we incorporated these velocity processes into continuous time random walk models for the prediction of pre-asymptotic to asymptotic anomalous transport. These types of model fit particularly well our observations since they allow for the mimicking of the spatial velocity persistence of fluid particles. Our effective upscaling approach proves itself to be accurate, computationally fast, and predictive in the sense that it does not rely on transport properties but rather on flow and medium information. It also permits transport predictions for any injection condition since it allows for a conditioning on the initial velocity distribution of the particles. However, any upscaling model relies on support volumes on which parameters are considered constant. Consequently, we defined a representative elementary volume in terms of velocity statistics that allows for the upscaling of anomalous transport. This definition is intrinsically different from the classical REV definitions in the sense that the velocity statistics control the transport at all scales. Thus, continuous time domain random walk models for the prediction of the transport can be defined through this new REV 7.2. OPEN QUESTIONS AND FUTURE WORK 159 definition. Continuous time domain random walk models at the scale of the velocity statistics based-REV capture both pre-asymptotic and asymptotic transport regime observables.

Yet, the impact of diffusion cannot be modeled explicitly from the property of the velocity statistics. Advection is a mechanism that causes a velocity decorrelation on a spatial scale.

Conversely, diffusion provokes velocity changes of solute particles on a temporal scale.

We derived a methodology to account for this diffusive impact and to couple these velocity variations through the time-domain and continuous-time random walk framework. Then, we elaborated on the temporal scales on which Fickian transport might be recovered. This gives powerful tools for fast transport predictions based on flow and medium properties for any Péclet number scenario. After having understood and modeled effectively conservative non-reactive transport at pore-scale in heterogeneous media, we extended the study to reactive transport. We quantified the impact of the pore-scale heterogeneities on chemical reactions and on mixing. We found out that the heterogeneity of the geometry (and therefore of the velocity field) enhances the mixing and the reaction rate which display similar behaviors. Based on the close relationships between mixing and chemical reaction, we validated a predictive method based on estimates of the mixing volume for the reaction mass computation. This methodology captures accurately the changes of rate caused by the heterogeneity of the velocity field. A logical next step would therefore be to develop a methodology for the computation of the mixing volume from these one dimensional upscaling CTRW models. It would then allow for the prediction of transport observables and chemical reactions from flow and medium properties only.

Open Questions and Future Work

In this section, we introduce some of the open questions that arise from this thesis:

• The continuous time random walk models we propose rely on the representative Eulerian distribution of the medium for its transport predictions. While significant progress has been recently made to relate the pore geometry to velocity distributions, it still remains an open issue. Our methodology would strongly benefit from these advances.

• In Chapter 5, we introduced a continuous time random walk that includes a diffusion process. This models needs to be validated versus direct numerical simulations for example.

• Particle image velocimetry (PIV) is a powerful tool to sample the Lagrangian velocities tracer particles experience. Therefore they can be tools to first, provide measured This situation (see Figure 8.3) is a slightly more complicated one, we need to identify the line segment where the velocity sign changes, and judging by the particle's original place (to the left of the segment or to the right of it) the exit will be different. In other words u 2u 1 ∆x (xx 1 ) + u 1 > 0 then x > x 1 -u 1 ∆x u 2u 1 and so x out := x 2 , (8.3) u 2u 1 ∆x (xx 1 ) + u 1 < 0 then x < x 1 -u 1 ∆x u 2u 1 and so x out := x 1 , (8.4) u 2u 1 ∆x (xx 1 ) + u 1 = 0 then x = x 1 -u 1 ∆x u 2u 1 and so ∆t x = +∞. (8.5)

In the first case, ∆t x is then computed the same way as in the case of the two positive velocities; in the second situation it is computed the same way as in the case of the two negative velocities; and finally in the third case, the particle stands at a zero velocity location, which means that there will be no displacement in the x-direction. Integrating with respect to y, we get v(x, y) = 2u 1 ∆x 2 (x 2x)(yy 1 ) + v(x, y 1 ). (8.9)

Since v is a function of (x 2x), we can set v(x, y 1 ) = a(x 2x) with the constant a to be determined. Besides, we can use the following equality derived from the mass equation The 5 other situations where there is only one solid voxel (either on top, below, on the right, behind or in front of our cell) can be solved in a symmetrical way.

u 1 ∆x = v 2 -v

Two of the Neighboring Voxels are Solid

This condition can happen in fifteen cases that can be categorized in two groups: (1) when the two solid voxels are blocking opposing faces (this can happen in three cases, for the x-, y-or z-directions), as can be seen (1) We start with the first case. We did not develop the equations for those cases because they can not happen in our simulations. Indeed, since we made a mesh refinement, those cases have disappeared. This can be seen on Figure 8.7 and 8.8. We however recall that this mesh refinement is not motivated by the simplification of analytic developments but by the gain of accuracy on the Navier-Stokes flow computing, see Chapter 3.2.2.1. And after the mesh refinement, Figure 8.8 : The case of two solid voxels blocking opposing faces has disappeared. Similar situations can happen in the cases where 3, 4, 5 or 6 solid voxels are situated nearby the cell, which will be discussed in Section 8.1.2.3.

(2) The second case can not be avoided by a mesh refinement. We solve it like the others before. Here we must consider the cases where the solid voxels are on the left and CHAPTER 8. APPENDIX below the considered cell. The velocity in the x and y direction are interpolated as u = 2u 2 ∆x 2 ∆y (xx 1 ) 2 (yy 1 ), (8.23)

v = 2v 2
∆y 2 ∆x (xx 1 )(yy 1 ) 2 . (8.24)

We now examine the w equation. In order to obtain it we first used the zero divergence criterion

∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z = 0. (8.25)
We then have

∂ u ∂ x = 4u 2
∆x 2 ∆y (yy 1 )(xx 1 ), (8.26)

∂ v ∂ y = 4v 2
∆x∆y 2 (yy 1 )(xx 1 ), (8.27)

and Then we write w(x, y, z 1 ) = a(xx 1 )(yy 1 ).

∂ w ∂ z = - ∂ u ∂ x - ∂ v ∂ y = 4
(8.30) with a to be determined. However, we know that This situation can happen in twenty different ways depending on the locations of solid voxels and can be categorized in two groups. The first group is when two of the three solid voxels are blocking the same direction (for example both in the x-direction) and the third one is in another direction. The second group describes the situation where three surrounding solid voxels are blocking all three coordinate directions, one solid voxel on the x-, the other one on the yand the last one on the z-direction face. This condition can happen in eight different ways.

w
(1) We start with the first case. As seen before, in the two cells' case, this situation has been avoided thanks to a mesh refinement (see Figures 8.7 and 8.8).

(2) Now we consider the case where the three surrounding solid voxels are blocking all three coordinate directions. We interpolate the velocities as As we have seen before, the cases four, five and six cells are avoided (see Figure 8.7 and 8.8. Besides that, as soon as we reach five solid voxels, there is no flow at all inside the cell.

Extreme Cases

We denote by extreme cases all the situations that can not be handled by the equations below (i.e. when a division by zero appears). In these cases, new equations have to be solved. We do not develop these equations in this report. As stated in the main text, see Section 2.2.2.2, the particle diffusive jump during a time step ∆t in each component direction can be computed as (8.56) where ξ → N (0, 1). However, Gaussian distributed numbers brings two issues that we want to avoid. The first one is the cost of generating Gaussian distributed numbers than can drastically slow down the simulation since many of them have to be generated and that we therefore want to reduce. The second is to have a better control on the norm of the diffusive jump. Indeed while the occurrence of numbers n / ∈ [-4, 4] is very low, a Normal distribution is technically not bounded and a number out of that range can be generated. It can causes particles to jump over solid mesh cells. To avoid these issues we compute the diffusive jump as (8.57) where γ → U (-1, 1) (uniformly distributed between -1 and 1). The Central Limit Theorem insures that the sum of the jumps converges toward a Gaussian distribution as the number of jumps grows. Multiplying by a factor of √ 3 insures that the variance is equal to

∆x di f f = √ 2∆tDξ ,
∆x di f f = √ 3 √ 2∆tDγ = √ 6∆tDγ,

Time Step Choice

The time step is chosen small enough so that the maximum possible jump does not bypass the length of a mesh cell. This choice has for goal to avoid that a particle jump over a solid cell which is nonphysical and to avoid limit errors. In practice, to satisfy the first condition the time step has to be chosen as ∆t < ∆x 2 6D , (8.58)

where ∆x is the length of a cell.

Boundary Conditions

At the void-solid interfaces of the mesh we implemented reflective boundary conditions for the diffusive jumps. In practice a particle that would end up in a solid voxel of the domain is reflected in the normal direction to the jump of the distance it would have traveled in the solid cell.

Depending on the type of results needed we implemented two types of boundary conditions at the inlet and outlet of the domain. The first condition is an absorbing boundary condition that causes the particles to get removed from the system when they reach these locations. The second condition is a reinjection boundary condition that reintroduces the particles at the inlet when they reach the outlet (or at the outlet when they reach the inlet).

Theses conditions are handy when in the need of running very long simulations (computation of moments for example). The reinjection location is chosen such that the velocity difference between the exit and reinjection location is acceptable. We consider the following one-dimensional scheme

x n+1 = x n + c χ , ( 8
.59) (8.60) where the transition times τ n are defined as (8.61) which can be seen as the harmonic mean of the advective and diffusive characteristic time.

t n+1 = t n + τ n ,
τ n = 1 D 2 c + v c = 1 1 τ D + 1 τ v ,
This definition aims to account for both the advective and diffusive impact. It is a function of the Péclet number. To highlight this point, we can be rewrite it as

τ n = τ v P -1 e + 1 , (8.62)
where P e is the Péclet number. This definition follows the time domain random walk methodology for the computation of the transition times. The definition of the transition times ψ(τ) can be expressed as

ψ(τ) = P s (v)δ t - 1 1 τ D + v c dv. (8.63) Setting f (v) = δ (t -1 1 τ D + v c
) and using the properties of the δ -Dirac functions, we can write

f (v) = 1 | f (v)| δ (v -v 0 ) with f (v 0 ) = 0. (8.64)
Then we have

f (v) = 1 c 1 ( 1 τ D + v c ) 2 (8.65)
and since every term is positive:

f (v) = 1 c 1 1 τ D + v

Laplace Transformation

Note that transforming to Laplace space, we have p(x, λ ) = R(x, λ ) 1ψ(λ ) where λ ψ(λ ) 1-ψ(λ ) ≡ M(λ ) is our memory function. [START_REF] Kullback | On information and sufficiency[END_REF], also called divergence, we elaborate on its use in this section. The KL divergence is a measure of how similar (or how different) two probability density functions that are defined on the same event space are. It is classically defined as (using similar notations to the ones used in Chapter 5) 8.83) where it represents the distance between the distributions p 1 and p 2 . It has been widely used to compare different distributions, particularly in information theory (see, for example, [START_REF] Bigi | Using kullback-leibler distance for text categorization[END_REF][START_REF] Robert | Statistical equilibrium states for two-dimensional flows[END_REF][START_REF] Lindgren | Universality of probability density distributions in the overlap region in high Reynolds number turbulent boundary layers[END_REF][START_REF] Carpineto | An information-theoretic approach to automatic query expansion[END_REF]. This distance can be seen as a measure of the relative entropy. As the KL distance decreases (and therefore the relative entropy), the two distributions are more and more similar, and inversely, the larger d KL (p 1 , p 2 ) and the more different the PDFs are. Since d KL (p 1 , p 2 ) = 0 ⇔ p 1 = p 2 by definition of the KL divergence, the two distributions are considered equivalent once the distance reaches 0. We detail in the following how to accurately compute the KL divergence between two distributions and detail thee edge cases that may occur.

d KL (p 1 , p 2 ) = ∞ 0 dvp 2 (v) ln p 2 (v) p 1 (v) , (

Edge Cases

Another details that need to be discussed are the cases where either p 1 (v) or p 2 (v) contribution is zero. In the case that p 2 (v) is zero, the contribution p 2 (v) ln p 2 (v) p 1 (v) is zero, since lim x→0 + x ln(x) = 0. Inversely, if p 1 (v) = 0, d sKL (p 2 , p 1 ) is set to be infinite. This can be an issue for discrete distributions that have a category (or bin) that is empty (such as our velocity distributions that have little statistics at early distances for example) since it would result in an overall distance KL d sKL (p 2 , p 1 ) = ∞. This can be tackled by introducing a back-off probability for p 2 (v i ) is v i does not occur. This type of back-off probability to bypass empty bins and data sparseness issues have been studied and used in statistical language modeling and text categorization (see, for example, [START_REF] Mori | Spoken dialogues with computers[END_REF][START_REF] Bigi | Using kullback-leibler distance for text categorization[END_REF]. The distribution can then be modified as

p 1 (v i ) = α p 1 (v i ) if v i occurs in the PDF, ε otherwise,
where ε is a threshold probability for the empty bins that is chosen to be small. Note that since ∑ n i=0 p 1 (v) has to be equal to 1, the rest of the PDF (non-empty bins) are rescaled by a value α that reduces the weight of the bins while maintaining the shape of the distribution. The value α can easily be estimated from (8.85) where |b i | denotes the number of occurrences of velocities inside the bin b i . Therefore

1 -∑ |b i | =0 α p 1 (v i ) = ∑ |b i |=0 ε,
∑ |b i | =0
and ∑ |b i |=0 represent the sum over the empty and non-empty bins, respectively. Finally, since the values of the KL distance can vary largely from a comparison between two distributions to another comparison between two other distributions, it causes issue when in the need of comparing the speed of two different evolutions (for example in Chapter 5 when we observe the evolution to a limited distribution in the Eulerian sense and in the Lagrangian sense). In this case, a solution is to divide the values given by the KL distance by the KL distance between the target distribution and an empty distribution (i.e. that contains a probability value of ε in every bin). where 0 denotes an empty distribution.

In other cases, another possibility of defining a normalized KL divergence is to divide the KL divergence d sKL (p 1 , p 2 ) by the initial distribution of velocity p 0 1 . Where p 0 1 corresponds to the first computation on a small volume for the Eulerian distribution for example.

1. 1

 1 Illustrative representation of anomalous vs. Fickian transport. . . . . . . . . 1.2 Illustrative example of a non-local, one-dimensional continuous time random walk. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Illustrative representation of pore-scales mechanisms triggering anomalous transport. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Digitized two-dimensional and three-dimensional geometries for flow and transport. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Illustrative example of a resolved velocity field on a squared mesh. . . . . . 2.3 Mesh cell exhibiting positive velocities in the x-direction. . . . . . . . . . . 2.4 Illustrative example of the Pollock algorithm. . . . . . . . . . . . . . . . . 2.5 Parabolic velocity profile in a throat. . . . . . . . . . . . . . . . . . . . . . 2.6 One solid cell is on the right of the voxel. . . . . . . . . . . . . . . . . . . 2.7 Illustrative example of the modified Pollock algorithm. . . . . . . . . . . . 2.8 Sequential advective and diffusive displacement of a particle. . . . . . . . . 2.9 Sequential advective, diffusive and reactive steps of a pair of reactant particles. And the displacement of the particle created by the chemical reaction. 3.1 Numerical cross-section cropped in the Berea sandstone sample used for the flow simulations, and 3-D sample volume with 20 particle streamlines. . 3.2 Time and space velocity magnitude series of a particle traveling through the sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Autocorrelation function of the s-Lagrangian velocity v s (s). . . . . . . . . . 3.4 Streamwise t-Lagrangian PDF P(v), s-Lagrangian PDF P(v) and the flux-weighting relation (3.17). . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES ix 3

3. 7

 7 Spatial velocity transition matrices computed with 10 6 particles and respectively with a spatial lag of ∆s = l p 150 , ∆s = l p 3 , ∆s = l p , and 8l p 3 (upper left to lower right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 12

 12 Comparison between the streamwise P(v), ensemble P(v), and mixed P m (v) of the s-Lagrangian velocity PDFs. . . . . . . . . . . . . . . . . . . 3.13 Comparison of the streamwise P(v), ensemble P(v) and mixed Pm (v) of the t-Lagrangian velocity PDFs. . . . . . . . . . . . . . . . . . . . . . . . 3.14 Ensemble t-Lagrangian PDF P(v) and s-Lagrangian PDF P(v) and the fluxweighting relation (3.58). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Illustration of the subdomain of the Berea sandstone sample under consideration and sample streamlines. . . . . . . . . . . . . . . . . . . . . . . . . x LIST OF FIGURES 4.2 Probability distribution function of the s-Lagrangian velocity p s (v), the initial velocities for uniform and flux-weighted injections p 0 (v), and the Eulerian velocity magnitude p e (v). . . . . . . . . . . . . . . . . . . . . . . . 94 4.3 Breakthrough curves for a flux weighted injection obtained from DNS, Bernoulli CTRW and OU CTRW at planes located at 6 p , 36 p , and 200 p . 96 4.4 Breakthrough curves for uniform injection obtained from the DNS, the Bernoulli CTRW, and the OU CTRW at planes located at 6 p , 36 p , and 200 p . Inset: Comparison between the uniform and the flux weighted BTCs for the direct simulation at the plane at 36 p . . . . . . . . . . . . . . 96 4.5 Evolution of the DNS, OU model, and Bernoulli model mean flow direction propagators for (top panel) uniform injection and (bottom panel) fluxweighted injection at times t = 5 • 10 -2 , 50, 5 • 10 3 τ c . . . . . . . . . . . . . 98 4.6 Evolution of the mean displacement and displacement variance in time for both uniform and flux-weighted injections. The DNS and OU models results are displayed with the t/ ln(t) and t 2 / ln 3 (t) late time scalings. . . . . 100 5.1 Porosity measured in a cube of increasing side length cube (x-axis) centered in the middle of the Berea sandstone sample studied in Section 5.2. The side length of the full sample is 6 p and its porosity is 0.182 (gray dashed line). p ≈ 1, 5 • 10 -4 m denotes the average pore length. . . . . . . 111 5.2 Time series of velocity magnitude experienced by a particle in the threedimensional digitized Berea sandstone sample. . . . . . . . . . . . . . . . 115 5.3 Evolution the Eulerian velocity PDF computed on successive cubes centered in the middle of the sample and evolution of the KL distance between the full sample Eulerian velocity PDF and the successive Eulerian velocity PDFs computed on growing volumes. . . . . . . . . . . . . . . . . . . . . 119 5.4 the Lagrangian s-velocity PDF along an ensemble of particle streamlines at different distances together with the steady-state Lagrangian s-velocity PDF. And evolution of the KL distance between the steady-state Lagrangian velocity PDF and successive Lagrangian velocity PDFs computed at growing distances in the mean flow direction. . . . . . . . . . . . . . . . . . . . 119 5.5 Temporal evolution of the mean flow direction dispersion for the direct particle tracking simulation and the uncorrelated CTRW model for uniform and flux weighted injection conditions together with the theoretical late time scaling t ln 3 (t) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 LIST OF FIGURES xi 5.6 Temporal evolution of the mean flow direction dispersion computed from the direct particle tracking simulation. We display the fully advective case P e = ∞, P e = 700, and P e = 40. . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Evolution of the asymptotic normalized dispersion coefficient (D/D 0 ) in function of the Péclet number together with dispersion scaling for high Péclet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Average number of pores traveled before reaching asymptotic dispersion in function of the Péclet number. . . . . . . . . . . . . . . . . . . . . . . . . 6.1 (Top) Velocity magnitude field (m s -1 ) in the first synthetic heterogeneous porous medium. (Bottom) Velocity magnitude field (m s -1 ) in the second synthetic heterogeneous porous medium. . . . . . . . . . . . . . . . . . . . 6.2 Illustrative zoomed snapshot of the reactive particle tracking simulation in the second medium. White particles represent the A reactant already in the medium, the red particles depict the B reactant entering the domain, and the green dots symbolize the C species being formed at the interface. The fingering induced by the heterogeneity of the flow field enhances the mixing between the two reactants and therefore accelerates the chemical reaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Evolution of the product mass m C (t) for Pe = 60 from the RPT simulation

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Evolution of the concentration distribution evolving from a point injection. . 6.5 (left) Concentration distribution integrated in the x-direction (red diamonds) and predicted Gaussian concentration using σ x and m x (t | x ) (black solid line). (right) Concentration distribution integrated in x (red diamonds) and predicted Gaussian concentration using σ y and m y (t | y ) (black solid line). Both plots correspond to a point injection at y = 4.5 × 10 -4 m at t = 68τ v for Pe = 60 in the first medium, which corresponds to the second column of grains. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Evolution of σ 2 e (t) and σ 2 a (t) from the RPT numerical simulation for the studied case characterized by Pe = 60. . . . . . . . . . . . . . . . . . . . . xii LIST OF FIGURES 6.7 Evolution of σ 2 e (t) (red solid line) and σ 2 a (t) (blue dashed line) from the RPT numerical simulation for the case of interest (second medium, Pe = 40). The lower and upper gray solid lines indicate the 2Dt and 2D h t scaling respectively. D = molecular diffusion. D h = hydrodynamic dispersion. . . . 149 6.8 Evolution of the mass m C (t) for Pe = 60 from the RPT simulation (symbols), from the dispersive lamella parameterized by the apparent variance σ 2 a (t) (blue dashed-line) and the effective variance σ 2 e (t) (red line) for the first medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 6.9 Evolution of the mass m C (t) for Pe = 40 from the RPT simulation (symbols), from the dispersive lamella parameterized by the apparent variance σ 2 a (t) (blue dashed-line) and the effective variance σ 2 e (t) (red line) for the second medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.10 Evolution of the porosity experienced by the plume φ (t) originated from a line injection at x = 5 × 10 -5 in the second medium, which corresponds to the middle of the first empty layers in the geometry. . . . . . . . . . . . . . 152 6.11 Evolution of the mass m C (t) for Pe = 40 from the RPT simulation, from the dispersive lamella parameterized by the apparent variance σ 2 a (t) and the effective variance σ 2 e (t) and their respective corrections using the porosity experienced by the plume in time φ (t) for an initial interface between A and B outside of the domain. . . . . . . . . . . . . . . . . . . . . . . . . . 153 8.1 Negative velocities on both faces . . . . . . . . . . . . . . . . . . . . . . . 164 8.2 Positive velocity on the x 1 -face, negative on the x 2 -one. . . . . . . . . . . . 165 8.3 Negative velocity on the x 1 -face, positive on the x 2 -one. . . . . . . . . . . . 165 8.4 One solid cell is on the right of the voxel. . . . . . . . . . . . . . . . . . . 166 8.5 There are two solids next to the empty one, and they are on the same axis. . 168 8.6 There are two solids next to the empty one, blocking two different directions.168 8.7 Mesh before the refinement. . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.1: Illustrative representation of anomalous transport (blue solid line) vs. Fickian transport (red dashed line) at five distinct times. (Top): Evolution of the plume; (Bottom):Vertically-integrated concentration. The solute plume deviates from Fickian behavior with early breakthrough and late tailing. This is reflected in the non-Gaussian vertically averaged concentration.

Figure 1

 1 Figure 1.2: Illustrative example of a non-local, one-dimensional continuous time random walk.Each row displays the state of the system at a different time. A particle (red) jumps from a position x at time t 1 to arrive at a far location after a short time interval at time t 2 , illustrating spatial non-locality. A particle (blue) jumps from a position x at time t 1 to arrive at a nearby location after a long time interval at time t 4 , illustrating temporal non-locality.

Figure 1

 1 Figure1.3: Illustrative representation of some pore-scales mechanisms triggering anomalous transport. A particle (blue) experiences a preferential channel of low tortuosity, whereas another particle (red) experience a more tortuous pathway. Finally, the last particle particle (green) is trapped in a region of low velocity and will need a large time to escape. These mechanisms, together with the heterogeneity of the velocity field they experience, enhance anomalous transport.

  as a time domain random walk by equidistant spatial steps according to random velocities and thus random transition times. Particle velocities are modeled as stationary Markov processes, which evolve along streamlines on the mean pore length. The streamwise motion is projected onto the mean flow direction according to the tortuosity. The stochastic particle model accurately predicts the non-Fickian transport dynamics obtained from direct numerical simulations of particle transport in a three-dimensional digitized Berea sandstone sample. It captures all features of transport and sheds light on the dependence of the upscaled transport behavior on the flow distribution and on the initial particle distribution, which are both critical for the accurate modeling of dispersion from the pre-asymptotic to the asymptotic regimes.

  Figure 2.1. The details of the image acquisition and the segmentation procedure can be found in Chapter 3 or in Gjetvaj et al. (2015); Puyguiraud et al. (2019a). On the other hand, the flow, conservative transport, and reactive transport simulations studied in Chapter 6 were carried out on a two-dimensional rectangular domain composed of solid disks. In this case, such geometry (displayed in Figure 2.1b), was preferred in order to mimic the porous medium matrix phase. Details of the 2D rectangular geometry generation can be found in Chapter 6.

Figure 2 . 1 :

 21 Figure 2.1: On the left, three dimensional digitized Berea sandstone sample. The gray color represents the solid phase while the black denotes the void pore phase. On the right, a section of the two dimensional artificial geometry. The red disks represents the solid phase while the blue background symbolizes the void pore phase.

Figure 2 . 2 :

 22 Figure 2.2: Illustrative example of a resolved velocity field on a squared mesh. Granite color cells represent the solid phase while the white cells denote the pore space. The arrows represent the velocities obtained from the flow simulation.

CHAPTER 2 .

 2 Figure 2.3: Mesh cell exhibiting positive velocities in the x-direction.

Figure 2

 2 Figure 2.4: Example of the Pollock algorithm in two dimensions: a particle (blue dot) originating from location (x p , y p ) at time t follows its streamline to reach its cell exit location (x e , y e ) at time t e .

CHAPTER 2 .

 2 Figure 2.5: Parabolic velocity profile in a throat.

CHAPTER 2 .

 2 Figure 2.6: One solid cell is on the right of the voxel.

Figure 2 . 7 :

 27 Figure 2.7: Illustrative example of the modified Pollock algorithm. A particle (blue dot) originating from the bottom left of the mesh follows its streamline through different cell types. The streamline colours denote the type of cell (and therefore the type of equation to use) that the particle crosses. The green colour denotes the classical Pollock case while the others represent 3 of the new situations.

Figure 2

 2 Figure 2.8: Illustrative example of a particle starting at location (x(t), y(t)) performing sequentially the advective and diffusive displacements.

Figure 2

 2 Figure 2.9: (left panel) An A particle (blue dot) and a B particle (orange dot) subsequently perform their advective (black solid line) and diffusive (red dashed line) jumps before they collide and react. The orange disk surrounding the B particle indicates its support volume ∆V = πr 2 where r is the reaction radius. (right panel) The A and B particles have reacted with each other to form a C particle injected between the two of them. This particle then starts following the streamlines of the flow.

  and 0 otherwise. This injection condition represents the initial condition of a spatially uniform concentration distribution. The fluxweighted initial distribution distributes particles weighted by their initial velocity as ρ(a) = v 0 (a)I Ω 0 (a)

  Confidential manuscript submitted to Water Resources Research

Figure 1 :

 1 Figure 1: Left panel: numerical cross-section cropped in the Berea sandstone sample used for the flow simulations. White and black color denotes the pore space and the solid phase, respectively.Right panel: 3D sample volume with 20 particle streamlines. Dark grey and light grey color denotes the pore space and the solid phase, respectively. The color scale of the streamlines denotes the velocity magnitude; from white for low velocity values to blue for high velocity values. The streamline colored from black to red is the one for which the space and time velocity series is reported in Figure2.

Figure 3

 3 Figure 3.1: (left panel): Numerical cross-section cropped in the Berea sandstone sample used for the flow simulations. White and black color denotes the pore space and the solid phase, respectively. (right panel): 3-D sample volume with 20 particle streamlines. Dark gray and light gray colors denote the pore space and the solid phase, respectively. The color scale of the streamlines denotes the velocity magnitude; from white for low velocity values to blue for high velocity values. The streamline colored from black to red is the one for which the space and time velocity series is reported in Figure 3.2.

  proposed and analyzed the statistics of Lagrangian velocities sampled equidistantly along the mean flow direction. Le Borgne et al. (2007) considered the evolution of the probability density function (PDF) of such Lagrangian velocities and Gotovac et al. (2009) used them as the basis to analyze flow and travel time statistics in heterogeneous porous media. Cvetkovic et al. (1991) compared isochronously and equidistantly sampled Lagrangian velocities for one-dimensional steady flow, Cvetkovic et al. (

  Figure 1: aaaaa

  Figure 3.3: Autocorrelation function of the s-Lagrangian velocity v s (s).

Figure 3

 3 Figure 3.4: Streamwise t-Lagrangian PDF P(v) (full circles), s-Lagrangian PDF P(v) (open circles) and the flux-weighting relation (3.17) (solid line). PDF = probability density function.

  a)] ρ(a)da. (3.20) Confidential manuscript submitted to Water Resources Research

Figure 8 :

 8 Figure 8: Left panel: Spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the (pink open circles) flux-weighted initial particle distribution p 0 (v) to the (red solid line) steady s-Lagrangian PDF P(v). Right panel: Spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the (navy blue open circles) uniform initial particle distribution p 0 (v) to the (red solid line) steady s-Lagrangian PDF P(v). The results were computed using 10 7 particles.

Figure 3

 3 Figure3.5: (left panel): Spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the flux-weighted initial particle distribution p 0 (v) (pink open circles) to the steady s-Lagrangian PDF P(v) (red solid line) at distances s = 0, 1/5, 2/5, 5 p . (right panel): Spatial evolution of the ensemble s-Lagrangian PDF P(v, s) from the uniform initial particle distribution p 0 (v) (navy blue open circles) to the steady s-Lagrangian PDF P(v) (red solid line) at distances s = 0, 4/3, 4, 8 p . The results were computed using 10 7 particles. PDF = probability density function.

Figure 9 :

 9 Figure 9: Left panel: Temporal evolution of the ensemble t-Lagrangian PDF P(v, t) from the (red open circles) flux-weighted initial particle distribution p 0 (v) to the (navy blue solid line) steady t-Lagrangian PDF P(v). Right panel: Temporal evolution of the ensemble t-Lagrangian PDF P(v, t) from the (light blue open circles) uniform initial particle distribution p 0 (v) to the (navy blue solid line) steady t-Lagrangian PDF P(v). The results were computed using 10 6 particles.

Figure 3

 3 Figure 3.6: (left panel): Temporal evolution of the ensemble t-Lagrangian PDF P(v,t) from the flux-weighted initial particle distribution p 0 (v) (red open circles) to the steady t-Lagrangian PDF P(v) (navy blue solid line) at times t = 0, 5, 50, 500, 10 4 τ c . (right panel): Temporal evolution of the ensemble t-Lagrangian PDF P(v,t) from the uniform initial particle distribution p 0 (v) (light blue open circles) to the steady t-Lagrangian PDF P(v) (navy blue solid line) at times t = 0, 10 3 , 10 4 τ c . The results were computed using 10 6 particles. PDF = probability density function.

CHAPTER 3 .

 3 STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIESprocess is stationary. The steady s-and t-Lagrangian statistics are related by flux-weighting according to (3.19).

61

  Figure 1: aaaaa

  ds da, (3.27) where I(•) is 1 if its argument is true and 0 otherwise, and ∆s = ss . The empirical transition probability and the conditional probability density r(v, ss |v ) are related in terms of the stationary joint PDF (3.25) as

Figure 3 . 8 :

 38 Figure3.8: Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for direct particle tracking (open circles) and Markov model (solid lines) simulations from the uniform initial particle distribution p 0 (v) (blue) to the steady s-Lagrangian PDF P(v) (red) at distances s = 0, 8/3, 16/3, 8 p . The results were respectively computed with 10 7 and 10 8 particles for the direct simulation and the model. PDF = probability density function.

Figure 3 . 9 :

 39 Figure 3.9: Evolution of the mean s-Lagrangian velocity v s (s) (circles) for a uniform initial distribution and expression (3.33) (solid line) with c = 2.5 p .

Figure 3

 3 Figure3.10: Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for both direct particle tracking (open circles) and Bernoulli process (solid lines) simulations from the uniform initial particle distribution p 0 (v) (blue) to the steady s-Lagrangian PDF P(v) (red) for distance s = 0, 8/3, 10 p . The results were respectively computed with 10 7 and 5•10 8 particles for the direct simulation and the model. PDF = probability density function.

Figure 3 .

 3 Figure 3.10 compares the evolution of p(v, s) obtained from the direct numerical simulations with the prediction (3.32) of the Bernoulli model for uniform injection conditions.

Figure 3

 3 Figure3.11: Evolution of the ensemble spatial Lagrangian velocity PDF P(v, s) for direct particle tracking (open circles) and Ornstein-Uhlenbeck process (solid lines) simulations from the uniform initial particle distribution p 0 (v) (blue) to the steady s-Lagrangian PDF P(v) (red) at distance s = 0, 4/3, 4, 8, 10 p . The results were respectively computed with 10 7 and 5 • 10 8 particles for the direct simulation and the model. PDF = probability density function.

Figure 3 .

 3 Figure 3.11 compares the results of the velocity Markov model based on the OU process with the data from the direct numerical simulations. The Markov model is capable of predicting the evolution of p(v, s) in every aspect at small, intermediate and large distances from the inlet. The velocity dependence of convergence rates in the OU-based Markov model for v s (s) accurately captures the evolution of p(v, s) for all velocity classes.

CHAPTER 3 .

 3 STOCHASTIC DYNAMICS OF PORE-SCALE VELOCITIES the random time τ = ∆s/v s . s n+1 = s n + ∆s, t n+1 = t n + τ n . (3.42) The Markov property of the s-Lagrangian velocity is transferred to the transition times τ whose distribution evolves in time just like the s-Lagrangian velocity PDF. CTRW formulations that are based on a single transition time distribution ψ(t) are not able to model the impact of non-stationary initial conditions from the s-Lagrangian point of view, or stationary initial conditions from a t-Lagrangian point of view. Note that the process (3.42) describes particle motion along a tortuous streamline. The motion in three-dimensional Cartesian coordinates can be obtained either by an additional characterization of the direction vector ω ω ω(s, a) = v[x(s, a)]/v e [x(s, a)] in equation (3.6) as a stochastic process, or the projection of the streamwise motion on the mean flow direction in terms of the advective

Figure 3

 3 Figure 3.13: Comparison of the streamwise P(v) (orange dashed line), ensemble P(v) (blue circles, computed with 10 6 particles) and mixed Pm (v) (solid black line, computed with 10 4 particles) of the t-Lagrangian velocity PDFs.

  particle position a on the position x(t, a) at time t. The t-Eulerian velocity PDF corresponding to Pt (v,t) is obtained by volumetric sampling of the Eulerian velocity magnitude v e (x) in the subdomain Ω

Figure 3

 3 Figure 3.14: Ensemble t-Lagrangian PDF P(v) (Full circles, computed with 10 6 particles) and s-Lagrangian PDF P(v) (open circles, computed with 10 7 particles) and the fluxweighting relation (3.58) (solid line). PDF = probability density function.

Figure 3 .

 3 Figure3.14 compares the ensemble t-and s-Lagrangian statistics at t = 9 • 10 3 τ c and s = 7 p . Both statistics are in good agreement with respect to relations (3.17) and (3.58), which confirms again that ergodic conditions are attained. Furthermore, note that under ergodic conditions the steady s-Lagrangian PDF P(v) is related to the s-Eulerian PDF by relation (3.50), while the steady t-Lagrangian is equal to the t-Eulerian PDF, see (3.53). Thus the s-and t-Eulerian velocity PDFs are identical under ergodic conditions,

  .63) where we used that ∇ • v(x) = 0. We obtain for the derivative of v s (s, a) = v e [x(s, a)] with respect to s, dv s (s, a) ds = v[x(s, a)] • ∇v s (s, a) , a) = -1 v s (s, a) dv s (s, a) ds J(s, a). (3.65) Integrating the differential Eq. (3.65) for the initial condition J(s = 0, a) = 1 yields J(a, s) = v e (a) v e [x(s, a)]

= Ω 0 δ

 0 (vv e [x(s(t), a)])ρ(a)da, (3.79) since x(t, a) = x(s(t), a). Then we can write it as p(v,t) = Ω 0 δ [ss(t, a)]δ [vv e (x[s(t), a])]ρ(a)dads.

  [START_REF] Porta | Continuum-scale characterization of solute transport based on pore-scale velocity distributions[END_REF] derived a mobile-immobile model to upscale pore-scale transport accounting for information on the pore space and the pore-scale velocity distributions.[START_REF] Meyer | Pore-scale dispersion: Bridging the gap between microscopic pore structure and the emerging macroscopic transport behavior[END_REF] used a Langevin approach to account for the impact of pore-scale velocity heterogeneity on solute dispersion. Due to their central role for transport upscaling from the pore to the Darcy scale, pore-scale particle velocities and their relation to the flow velocity and porous medium structure have been the subject of recent research (de Anna

  .1) where p(x) and ν denote the pressure and the kinematic viscosity of the fluid, respectively. The associated mass conservation equation is ∇ • v(x) = 0. The position vector is denoted by x = (x 1 , x 2 , x 3 ) . The problem is solved by fixing constant pressure at both the inlet and outlet boundaries of the sample and no-slip condition at the void-solid interfaces and at the other physical boundaries of the sample. Details concerning the computations and sample characteristics are given in section 4.2.3. The magnitude of the Eulerian velocity v(x) in the following is denoted by v e (x) = v(x) . The probability density function (PDF) of v e (x) is denoted by p e (v). It can be obtained by spatial sampling over a sampling volume that is representative of flow variability.

  Figure 4.1: Illustration of the subdomain of the Berea sandstone sample under consideration and sample streamlines.

Figure 4 . 2 :

 42 Figure 4.2: Probability distribution function of the s-Lagrangian velocity, p s (v) (black empty circles), the initial velocities for uniform (blue) and flux-weighted (red) injections, p 0 (v), and the Eulerian velocity magnitude p e (v) (black full circles).

Figure 4 . 3 :Figure 4

 434 Figure 4.3: Breakthrough curves for a flux weighted injection obtained from the DNS (circles), the Bernoulli CTRW (dashed line) and OU CTRW (solid line) at planes located at 6 p , 36 p , and 200 p (from dark blue to red).

Figure 4 . 5 :

 45 Figure 4.5: Evolution of the DNS (circles), OU model (solid lines), and Bernoulli model (dashed lines) mean flow direction propagators for (top panel) uniform injection and (bottom panel) flux-weighted injection at times t = 3.5 • 10 -1 τ c (blue), t = 3.5 • 10 2 τ c (orange), and t = 3.5 • 10 4 τ c (red).

Figure 4 .

 4 Figure 4.6 shows the evolution of the mean displacement for uniform and flux-weighted injection conditions. The early time behavior is in both cases linear and given by v 1 t, where

Figure 4 . 6 :

 46 Figure 4.6: Evolution of the (top panel) mean displacement and (bottom panel) displacement variance in time for both uniform (blue) and flux-weighted (red) injections. The DNS and OU models results are respectively displayed with open circles and solid lines. The black lines indicate the t/ ln(t) and t 2 / ln(t) 3 late time scalings.

  Figure5.1: Porosity measured in a cube of increasing side length cube (x-axis) centered in the middle of the Berea sandstone sample studied in Section 5.2. The side length of the full sample is 6 p and its porosity is 0.182 (gray dashed line). p ≈ 1, 5 • 10 -4 m denotes the average pore length.

Figure 5

 5 Figure 5.2: Time series of velocity magnitude of experienced by a particle in the threedimensional digitized Berea sandstone sample shown in Figure 5.1. The intermittency is more visible in the linear plot (top) while the second plot (bottom) gives more information on the low velocity behavior.

Figure 5 Figure 5

 55 Figure 5.3: (Left) The Eulerian velocity PDF computed on cubes centered in the middle of the sample of sizes V = 4 • 10 -8 V T (red solid line), V = 1 • 10 -6 V T (dark orange solid line), V = 10 -5 V T (light orange solid line), V = 0.125V T (light blue solid line) and V T (navy blue circles). (Right) The KL distance between the full sample Eulerian velocity PDF and the Eulerian velocity PDFs computed on growing volumes from V = 0 to V = V T . The x-axis represents the cube side length; the side length of the full sample is 6 p .

Figure 5

 5 Figure5.5: Dispersion in the direction of the mean flow versus normalized time computed using the direct particle tracking simulation (circles) and computed using the onedimensional uncorrelated CTRW model (solid line) for uniform (blue) and flux weighted (red-orange) injection conditions. The theoretical late time scaling t ln 3 (t) is materialized by the black solid line.

  induces a cut-off in the distribution of transition time because the maximum transition time over the length c is τ D = 2 c /D. This is incorporated into the CTRW scheme (5.12) by modeling the transition time τ as the harmonic sum between the purely advective and purely diffusive transition timesτ = τ D τ v τ D + τ v . (5.22) This definition aims to account for both the advective and diffusive impact. It is a function of the Péclet number. This definition follows the time domain random walk methodology

  Figure 5.6: Temporal evolution of the mean flow direction dispersion computed from the direct particle tracking simulation. We display the fully advective case P e = ∞ (blue circles), P e = 700 (red solid line), and P e = 40 (black solid line). The black dashed and red dashed lines represent the stabilization of the dispersion coefficient D predicted by the CTRW theory after diffusive time τ D for P e = 40, and P e = 700 respectively.

Figure 5

 5 Figure 5.7: Evolution of the asymptotic normalized dispersion coefficient (D/D 0 ) in function of the Péclet number. The gray dashed line indicates the low Péclet late time dispersion behavior D ∞ ≈ φ m with m = 2.1. The black solid line represents the dispersion scaling D ≈ P 2 e /(ln(P e ) 3 ) scaling for high Péclet numbers.

Figure 5

 5 Figure 5.8: Average number of pores traveled before reaching asymptotic dispersion in function of the Péclet number. The black solid line indicates a linear scaling with the Péclet number at late times.

Figure 6

 6 Figure 6.1: (Top) Velocity magnitude field (m s -1 ) in the first synthetic heterogeneous porous medium. (Bottom) Velocity magnitude field (m s -1 ) in the second synthetic heterogeneous porous medium.

  locity v = 3.45 × 10 -4 m.s -1 , and a the diffusion coefficient D = 4 × 10 -10 m 2 s -1 , for a resulting Péclet number P e = v L p D = 40. An illustrative zoomed snapshot of the running simulation in the second medium is displayed in Figure 6.2.

Figure 6 . 2 :

 62 Figure 6.2: Illustrative zoomed snapshot of the reactive particle tracking simulation in the second medium. White particles represent the A reactant already in the medium, the red particles depict the B reactant entering the domain, and the green dots symbolize the C species being formed at the interface. The fingering induced by the heterogeneity of the flow field enhances the mixing between the two reactants and therefore accelerates the chemical reaction.

Figure 6

 6 Figure6.3: Evolution of the product mass m C (t) for Pe = 60 from the RPT simulation in the first porous medium (symbols), and from the hydrodynamic dispersion coefficient (black solid line). The analytical solution overpredicts the product total mass in the medium due to the assumption of complete mixing between reactants. Note that similar behaviors have been observed by[START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF].

  .21) which is equivalent to the Langevin equation (6.4). The concentration distribution is represented in terms of the Green function g(x,t|x ) as c

Figure 6 Figure 6

 66 Figure 6.4: Evolution of the concentration distribution g(x,t | y ) evolving from a point injection at y = 3 × 10 -4 m at t = 0.3τ v (top) and t = 68τ v (bottom) for Pe = 60. jection is illustrated in Figure 6.4 for the first medium, which shows the evolution of the Green function g(x,t | y ) at two different times originating from a point source initially

Figure 6

 6 Figure 6.6: Evolution of σ 2 e (t) and σ 2 a (t) from the RPT numerical simulation for the studied case characterized by Pe = 60. The gray solid lines indicate the 2Dt (lower) and 2D h t (upper) behaviors. The insets corresponds to the spatial distribution of the pulse injection at t = 0.26τ v (top-left) and t = 24τ v (top-right). The vertical black dashed lines indicate the times that corresponds to the insets. D h = hydrodynamic dispersion.

Figure 6

 6 Figure 6.6, where the deformation of the plume is visible. Due to the greater heterogeneity of the flow field of the second medium, the plume experiences a larger distortions resulting in an larger overall width.

Figure 6

 6 Figure 6.8: Evolution of the mass m C (t) for Pe = 60 from the RPT simulation (symbols), from the dispersive lamella parameterized by the apparent variance σ 2 a (t) (blue dashedline) and the effective variance σ 2 e (t) (red line) for the first medium.

Figure 6

 6 Figure6.11: Evolution of the mass m C (t) for Pe = 40 from the RPT simulation (green dots), from the dispersive lamella parameterized by the apparent variance σ 2 a (t) (blue solid line) and the effective variance σ 2 e (t) (red solid line) and their respective corrections (blue and red dashed) using the porosity experienced by the plume in time φ (t) for an initial interface between A and B outside of the domain in the second medium.

Figure 8 . 2 :

 82 Figure 8.2: Positive velocity on the x 1 -face, negative on the x 2 -one.

Figure 8

 8 Figure 8.3: Negative velocity on the x 1 -face, positive on the x 2 -one.

8. 1 . 2

 12 Derivation of the Trajectory Equation in Function of The Spatial Situation of the CellIn this section we give the detail of the derivation for the resolution of the trajectory equations in function of the type of situation encountered. The situations where there are no solid voxel or only one nearby the considered mesh cell are discussed in the main text, see Section 2.2.1.2.

Figure 8

 8 Figure 8.4: One solid cell is on the right of the voxel.

  Figure 8.5; and (2) when the solid voxels are located on adjoining faces, which can happen in twelve cases, as shown on Figure8.6.

Figure 8

 8 Figure 8.5: There are two solids next to the empty one, and they are on the same axis.

Figure 8 . 6 :

 86 Figure 8.6: There are two solids next to the empty one, blocking two different directions.

Figure 8 . 8 :

 88 Figure 8.8: Mesh divided by three.

  y 1 )(x 1x)(zz 1 ) + w(x, y, z 1 ). (8.29)

Figure 8

 8 Figure 8.9: There are three solids next to the empty one and two of them block the same direction.

  Figure 8.10: There are three solids near the voxels, each one blocking a different direction.

Figure 8 .

 8 Figure 8.11: There are four solids near the voxel, theses cases can not happen in our simulations.

  we know that f (v 0 ) ≡ δ (tto rewrite Equation (8.63) using (8.64) and (8.66) asψ(τ) = P s (v) x ,tt )ζ (xx )ψ(tt ) dt dx + δ (x)δ (t), (8.75)where ζ is the distribution of lengths.

λ

  and R(x, λ ) = λ p(x, λ ) 1 1ψ(λ ) , (8.76)where the hat denotes Laplace-transformed functions and λ is the Laplace variable. Transforming Eq. (8.75) to the Laplace space we obtainR(x, λ ) = δ (x) + R(xx , λ )ζ (xx ) ψ(λ )dx (8.77)Using relation (8.76) we obtainλ p(x, λ ) 1 1ψ(λ ) = δ (x) + p(xx , λ ) λ ψ(λ ) 1ψ(λ ) ζ (xx ) dx (8.78) Noting that 1 1-ψ(λ ) = 1 + ψ 1-ψ we can write λ p(x, λ ) + λ ψ(λ ) 1ψ(λ ) p(x, λ ) = δ (x) + p(xx , λ ) λ ψ(λ ) 1ψ(λ ) ζ (xx ) dx ,(8.79) factorizing by λ ψ(λ ) 1-ψ(λ ) on the right side, we getλ p(x, λ ) = δ (x) + p(xx , λ )p(x, λ ) λ ψ(λ ) 1ψ(λ ) ζ (xx ) dx .

  La compréhension et la modélisation du transport hydrodynamique de solutés dans les milieux hétérogènes ont récemment connu des progrès significatifs. Cependant, il reste , nous identifierons les informations spatiales nécessaires au paramétrage du modèle de changement d'échelle de type CTRW. Nous paramétrerons le modèle CTRW en une dimension pour la prédiction du transport pré-asymptotique tout d'abord dans un cas purement advectif. Ensuite, nous définirons un volume élémentaire représentatif en termes des statistiques de vitesse eulériennes et lagrangiennes pour la prédiction du transport anormal à plus grandes échelles. Enfin, nous étudierons le comportement du mélange et la quantité de solutés réactifs produite lors d'une réaction simple de type A + B donne C dans les milieux hétérogènes.Le document de thèse est structuré de la façon suivante. Après une introduction (chapitre 1) reprenant en détail les éléments de réflexion et les connaissances acquises ainsi que les

	CHAPTER 1. INTRODUCTION
	Au chapitre 4, nous modélisons la dispersion de l'échelle des pores à l'échelle de
	Darcy dans un cas purement advectif. Nous paramétrons le modèle en termes de distri-
	bution des vitesses eulériennes et de tortuosité advective, qui sont des caractéristiques de
	l'écoulement, et de la longueur moyenne des pores qui est une caractéristique du milieu. Le
	mouvement stochastique des particules est modélisé grâce à une marche aléatoire dans le
	temps de type CTRW: les vitesses sont choisies aléatoirement dans une distribution station-
	naire et conservées pour une distance constante, cela résulte en des transitions temporelles
	également aléatoires. Les vitesses des particules sont modélisées à l'aide de processus
	stationnaires markoviens qui évoluent le long des lignes de courant et qui restent corrélés
	sur une distance de l'ordre d'un pore. Le déplacement le long des trajectoires est pro-
	jeté sur l'axe correspondant à la direction de l'écoulement moyen grâce à la tortuosité du
	milieu. Le modèle stochastique de particules prédit avec précision la dynamique du trans-
	port non-fickéen (pré-asymptotique) obtenue à partir de simulations numériques directes
	de transport de particules dans un échantillon 3D numérisé de grès de Berea. Le mod-
	èle reproduit toutes les caractéristiques du transport et met en lumière la dépendance du
	comportement du transport sur la distribution eulérienne et sur la distribution initiale des
	particules qui sont des attributs essentiels pour la modélisation précise de la dispersion à
	des régimes pré-asymptotiques jusqu'au régime asymptotique (fickéen).
	du pore.
	Ces deux objectifs principaux seront cependant décomposés en objectifs plus spéci-
	fiques adressant les interrogations décrites plus haut. Nous nous attacherons premièrement
	à relier mathématiquement les statistiques de vitesse échantillonnées de manière équidis-
	tante à celles échantillonnées de façon isochrone le long des lignes de courant et à valider
	ces relations par simulations directes effectuées à l'échelle du pore sur des échantillons
	soluté de l'échelle du
	pore à l'échelle de Darcy dans des milieux poreux complexes.

un grand nombre de questions en suspens, par exemple, concernant les relations entre les vitesses (lagrangiennes) du fluide à l'échelle du pore échantillonnées dans le temps et dans l'espace ainsi que leur lien avec les statistiques eulériennes. Des progrès sont également nécessaires pour parvenir à modéliser avec précision le signal spatial des vitesses des particules de solutés et réaliser le changement d'échelle pour obtenir des modèles macroscopiques cohérents avec les processus décrits à l'échelle du pore que le système soit en régime fickéen ou en régime pré-asymptotique (non-fickéen). Dans ce dernier cas, l'utilisation de modèles basés sur des marcheurs aléatoires permettant de simuler un processus de type Continuous Time Random Walk (CTRW) semble appropriée. Les modèles de changement d'échelle reposent cependant sur des paramètres qui sont considérés constant dans un volume élémentaire représentatif (VER). Définir un VER pour les modèles de type CTRW est donc nécessaire. Enfin, il est important de comprendre et de quantifier l'impact de l'hétérogénéité à l'échelle des pores sur le volume de mélange entre les réactifs et donc sur la dynamique du transport réactif qui en résulte. Dans cette optique, nous aborderons cette thèse avec deux objectifs principaux. Le premier objectif sera de comprendre et de caractériser l'impact de l'hétérogénéité des vitesses à l'échelle des pores sur la dispersion et le mélange du soluté. Le deuxième objectif consistera en développer un modèle efficace pour la prédiction du transport à l'échelle de Darcy basée sur des informations à l'échelle numérisés de roches réelles. Dans un deuxième temps nous développerons des processus stochastiques pour la modélisation de l'évolution en espace des vitesses des particules. En-suiteobjectifs résumés ci-dessus, nous décrivons la méthodologie utilisée pour caractériser et modéliser l'écoulement et le transport hydrodynamique dans le chapitre 2. Dans un premier temps, nous détaillons la technique d'acquisition d'images et les étapes successives permettant le calcul de l'écoulement et la détermination du champ de vitesses eulérien. Ensuite, nous présentons la mise en oeuvre de la méthode utilisée pour le suivi de particules: l'algorithme de traçage, le processus de diffusion et le calcul de la réaction.

Le chapitre 3 présente, dans un premier temps, une analyse détaillée de l'évolution des statistiques lagrangiennes et eulériennes et de leur dépendance à l'égard des conditions d'injection dans le but de définir les bases pour, dans un deuxième temps, prédire le transport à l'échelle de Darcy. L'étude présentée ici est basée sur les données de vitesse obtenues à partir de simulations numériques de la dynamique des fluides de l'écoulement de Stokes et du suivi de particules dans la structure poreuse tridimensionnelle d'un échantillon de grès de Berea obtenue à partir d'imagerie microtomographique à haute résolution par rayons X. Nous observons que les séries de vitesses échantillonnées de manière isochrone montrent un comportement intermittent, les séries équidistantes quant à elles, présentent des variations plus régulières. Les deux statistiques évoluent cependant vers des états stationnaires, liés aux statistiques de vitesse eulérienne. La distribution des vitesses lagrangiennes échantillonnées de manière équidistante converge après seulement quelques pores. Ces résultats indiquent que les vitesses équidistantes peuvent être représentées par un processus markovien et ergodique. En conséquence, nous démontrons qu'un modèle markovien de type stochastique est capable de reproduire l'évolution des statistiques de vitesses lagrangiennes. Le modèle est paramétré par la distribution des vitesses eulériennes et par une distance caractéristique de relaxation. Ces informations peuvent être reliées aux propriétés hydrauliques et à la géométrie du milieu. Ces résultats produisent une base solide pour construire une approche stochastique prédictive pour la dispersion de

  Under ergodic conditions, the s-Eulerian PDF is stationary and thus independent from s, P e (v, s) = P e (v) and µ e (s) = v e . The Lagrangian and Eulerian statistics are related by

	3.6. APPENDIX			73
	with the s-Eulerian mean velocity µ e (s). Equation (3.48) means that the s-Lagrangian and
	s-Eulerian velocity PDFs are related through flux weighting. The s-Eulerian mean velocity
	µ e (s) is
	∞
	µ e (s) =	P e (v, s)dv,	(3.49)
	0
	which in general evolves with distance s.
		s)	P e (v, s),	(3.48)

Table 6

 6

	.1.

  (blue dashed line) from the RPT numerical simulation for the case of interest (second medium, Pe = 40). The lower and upper gray solid lines indicate the 2Dt and 2D h t scaling respectively. D = molecular diffusion. D h = hydrodynamic dispersion.

		10 -6
		10 -7	σ 2 a (t) σ 2 e (t)
		10 -8
	(t)
	σ 2	10 -9
	10 -10
	10 -11
				10 0					10 1	10 2
								t/τ v
	Figure 6.7: Evolution of σ 2
	The total mass of C produced is obtained by integrating Equation (6.34) as
			m C (t) = φ	0	L y	dy dx	c 0 2	erfc		| x | 2σ 2 e	,	(6.35)
	which gives
			m C (t) = c 0 L y φ σ e (t)	2 π	.	(6.36)
										)	.	(6.34)

e (t) (red solid line) and σ 2 a (t)

  However, the flow mean on the y = y 1 face is equal to v 1 : Once solved, they allow us to get the path linesx(t) = x 2 -∆x 2 (x 2x p ) ∆x 2 + u 1 (x 2x p )t , (8.20) y(t) = y 1 -v 1 ∆y v 2v 1 + v 1 ∆y + (v 2v 1 )(y py 1 ) (v 2v 1 )∆y 1 + u 1 (x 2x p ) ∆z + (w 2w 1 )(z pz 1 ) (w 2w 1 )∆z 1 + u 1 (x 2x p )

	8.1. PARTICLE TRACKING DETAILS	167
									1 ∆x	x 2 x 1	v(x, y 1 )dx = v 1 .	(8.13)
	We therefore have
							1 ∆x	x 2 x 1	v(x, y 1 )dx = v 1 =	a∆x 2	,	(8.14)
	so
												a =	2v 1 ∆x	,	(8.15)
	and we finally obtain
		v(x, y) =	2(v 2 -v 1 ) ∆x∆y	(x 2 -x)(y -y 1 ) +	2v 1 ∆x	(x 2 -x).	(8.16)
	At this point, we obtain the following equations
		u =	u 1 ∆x 2 (x 2 -x) 2 ,	(8.17)
			v =	2v 1 ∆x	(x 2 -x) +	2(v 2 -v 1 ) ∆x∆y	(x 2 -x)(y -y 1 ),	(8.18)
		w =	2w 1 ∆x	(x 2 -x) +	2(w 2 -w 1 ) ∆x∆z	(x 2 -x)(z -z 1 ).	(8.19)
													∆x 2	t	2∆x(v 2 -v 1 ) u 1 ∆y	,	(8.21)
	z(t) = z 1 -	w 1 ∆z w 2 -w 1	+	w 1 ∆x 2	t	2∆x(w 2 -w 1 ) u 1 ∆z	. (8.22)
													∆y	1	,	(8.10)
	and so we obtain
					v(x, y) =	2(v 2 -v 1 ) ∆x∆y	(x 2 -x)(y -y 1 ) + v(x, y 1 ).	(8.11)
	Integrating v(x, y 1 ) between x 1 and x 2 yields
		1 ∆x	x 2 x 1	v(x, y 1 )dx =	1 ∆x	x 2 x 1	a(x 2 -x)dx =	1 ∆x	a(x 2 -x 1 ) 2 2	=	a∆x 2	.	(8.12)

  = x 1 + (x px 1 ) 1 -2(x px 1 )(y py 1 )(u 2 ∆y + v 2 ∆x) ∆x 2 ∆y 2 = y 1 + (y py 1 ) 1 -2(x px 1 )(y py 1 )(u 2 ∆y + v 2 ∆x) ∆x 2 ∆y 2 = z 1 + w 1 ∆x∆y u 2 ∆y + v 2 ∆x -( w 1 ∆x∆y u 2 ∆y + v 2 ∆x + z 1z p ) 0.5 -2(x px 1 )(y py 1 )[w 1 ∆x∆y(u 2 ∆y + v 2 ∆x) + (z 1z p )(u 2 ∆y + v 2 ∆x) 2 ]

	and
								a =	4w 1 ∆x∆y	.	(8.36)
	Replacing it in (8.29) we obtain
	w(x, y, z) =	4 ∆x∆y	(	u 2 ∆x	+	v 2 ∆y	)(y -y 1 )(x 1 -x)(z -z 1 )
								+	4w 1 ∆x∆y	(x -x 1 )(y -y 1 ). (8.37)
	And if we use the mass equation
								-u 2 ∆x	+	-v 2 ∆y	+	w 1 -w 2 ∆z	= 0,	(8.38)
	which is equivalent to
								w 1 -w 2 ∆z	=	u 2 ∆x	+	v 2 ∆y	,	(8.39)
	we finally obtain
	w(x, y, z) =	4(w +	4w 1 ∆x∆y	(x -x 1 )(y -y 1 ). (8.40)
	The three velocity equations are then
	u =	2u 2 ∆x 2 ∆y	(x -x 1 ) 2 (y -y 1 ),	(8.41)
	v =	2v 2 ∆y 2 ∆x	(x -x 1 )(y -y 1 ) 2 ,	(8.42)
	w =	4w 1					(8.43)
	1 = which allows us to get, after resolution 1 ∆x∆y	x 2 x 1	y 2 y 1	w(x, y, z 1 )dydx,	(8.31)
	= = x(t) -u 2 ∆y 1 ∆x∆y x 2 x 1 y 2 y 1 a(y -y 1 )(x -x 1 )dydx, 1 ∆x∆y x 2 x 1 a(x -x 1 ) u 2 ∆y+v 2 ∆x y 2 (y -y 1 ) 2 2 y 1 dx, = 4∆x∆y . a(y 2 -y 1 ) 2 (x 2 -x 1 ) 2 y(t) -v 2 ∆x u 2 ∆y+v 2 ∆x ,	(8.44) (8.45)	(8.32) (8.33) (8.34)
	and finally z(t) 1 w 1 = a∆y∆x 4 , 2 ∆x 2 ∆y 2	t	(8.35)

2w 1 ) ∆x∆y∆z (yy 1 )(xx 1 )(zz 1 ) ∆x∆y (xx 1 )(yy 1 ) + 4(w 2w 1 ) ∆x∆y∆z (xx 1 )(yy 1 )(zz 1 ),

2

.

(8.46)

  Inversing back to real time gives Since little details have been given in the main text about the Kullback-Leibler (KL) distance

	8.3. NOTES ON THE KULLBACK-LEIBLER DIVERGENCE	179
	8.3 Notes on the Kullback-Leibler Divergence
	8.3.1 Bases
	∂ p ∂t	(x,t) =	0	t	M(t -t )	-c χ	∂ p ∂ x	(x,t -t ) +	2 c 2χ 2	∂ 2 ∂ x 2 p(x,t -t ) dt .	(8.82)

  8.3.2 Using the Kullback-Leibler Divergence as a DistanceSurprisingly, while the KL divergence is often called a distance, it is not. Indeed, it does not fulfill the symmetry assumptions of classical distances since Equation (8.84) is not symmetric (d KL (p 1 , p 2 ) = d KL (p 2 , p 1 )). To tackle this issue the KL divergence can easily be symmetrized (see, for example,[START_REF] Bigi | Using kullback-leibler distance for text categorization[END_REF] asd sKL (p 1 , p 2 ) =where d sKL stands for symmetrical Kullback-Leibler. It results in d sKL (p 1 , p 2 ) being always equal to d sKL (p 2 , p 1 ), allowing for its use as a distance.

	∞ 0	dv[p 2 (v) -p 1 (v)] ln	p 2 (v) p 1 (v)	,	(8.84)

  The normalizedKullback-Leibler then

	reads:	d n sKL (p 1 , p 2 ) =	d sKL (p 1 , p 2 ) d sKL (p 1 , 0)	,	(8.86)

Acknowledgements

This work has been supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ ERC Grant 1. We observe no difference in the results between the two schemes.

PARTICLE TRACKING DETAILS

CHAPTER 6

Prediction of Pore-Scale Reaction in Complex Media

Chapter 5

Is there a Representative Elementary Volume for Anomalous Dispersion?

Introduction

The notion of a representative elementary volume (REV) lies at the heart of macroscopic (continuum) descriptions for systems that exhibit small scale structural and geometric disorder, and phase segregation (solid and void phase, for example), features which are usually referred to as material heterogeneity. In the frame of continuum approaches, the REV is associated to a point of the continuous field where average properties, that are supposed to denote the effective properties of the material, are allocated. For instance the permeability, from which the average fluid velocity is derived, and the hydrodynamic dispersion coefficient are critical properties for modeling steady state flow and solute transport, respectively.

The REV corresponds to the (minimum) volume required to evaluate the effective properties of a heterogeneous material or, in other words, the minimum volume above which the properties are stationary. This is illustrated in Figure 5.1 for the ratio of void to bulk volume. Porosity is defined as the constant limit value of the ratio φ between void and bulk volume

where Ω f denotes the fluid domain, Ω the bulk volume on a scale , V its volume, and I(•)

is the indicator function, which is equal to 1 if its argument is true and 0 otherwise. The length scale at which φ stabilizes defines the REV scale. Accordingly, the REV is clearly definable for two extreme cases: 1) unit volume in a periodic microstructure, and 2) a volume containing a large set of microscale structures displaying homogeneous and ergodic properties. The existence of the REV relies on the Chapter 6

Prediction of Pore-Scale Reaction in Complex Media

Introduction

The study of mixing-controlled chemical reactions in the subsurface has recently attracted significant attention due to its relevance in numerous applications including bioremediation of the subsurface [START_REF] Ding | Elimination of the reaction rate "scale effect": Application of the Lagrangian reactive particle-tracking method to simulate mixing-limited, field-scale biodegradation at the Schoolcraft (mi, usa) site[END_REF][START_REF] Scow | Natural attenuation and enhanced bioremediation of organic contaminants in groundwater[END_REF] and geological CO 2 sequestration [START_REF] Gérard | Miscible viscous fingering induced by a simple a+ b-> c chemical reaction[END_REF]. One of the main complexities in understanding the dynamics of these chemical reactions stems from the fact that natural aquifers are heterogeneous media at all scales. The geometric heterogeneity causes strong modifications of the flow and transport dynamics leading to variations of the mixing interfaces between reactants. Therefore, it alters the chemical reaction behavior. Chemical dynamics in the subsurface are affected by both the degree and the scale dependence of heterogeneity (Dentz et al., 2011b). This starts at the pore scale where the geometrical complexity leads to chemical reaction behavior that does not obey the transport laws established for homogeneous environments (Dentz et al., 2011b;[START_REF] Meile | Scale dependence of reaction rates in porous media[END_REF][START_REF] Li | Upscaling geochemical reaction rates using pore-scale network modeling[END_REF][START_REF] Berkowitz | Anomalous transport in random fracture networks[END_REF].

Mixing is the fundamental process that brings reactants into contact enabling chemical reactions. Mixing limits chemical reaction rates by controlling the time necessary to physically bring reactants into contact in solution, such as fluid-fluid reactions [START_REF] Jiménez-Martínez | Impact of saturation on dispersion and mixing in porous media: Photobleaching pulse injection experiments and shear-enhanced mixing model[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] De Anna | Mixing and reaction kinetics in porous media: An experimental pore scale quantification[END_REF][START_REF] Willingham | Evaluation of the effects of porous media structure on mixing-controlled reactions using pore-scale modeling and micromodel experiments[END_REF][START_REF] Gramling | Reactive transport in porous media: a comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF], or in the solid phase, such as solid-phase reactions [START_REF] Li | Upscaling geochemical reaction rates using pore-scale network modeling[END_REF][START_REF] Wintsch | Fluid-rock reaction weakening of fault zones[END_REF]. In well-mixed systems, chemical reactions occur at a thermodynamically determined rate. However, well-mixed conditions are unusual under natural conditions. In natural porous media, chemical reactions are influenced by the medium heterogeneity and the degree of mixing between reactants, CHAPTER 6. PREDICTION OF PORE-SCALE REACTION to turbulent [START_REF] Lester | Chaotic mixing in three-dimensional porous media[END_REF][START_REF] Villermaux | Mixing as an aggregation process[END_REF]. However, disregarding the eventual merging of lamellae due to transverse diffusion provokes incorrect quantifications of mixing in heterogeneous flows [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF]. Such discrepancies occur when the plume begins to sample the vertical flow contrast. To tackle this issue, [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF] proposed an effective approach they named dispersive lamella to account for the action of transverse diffusion in spatially variable flows. They validated this methodology against analytical solutions and direct reactive particle tracking simulations in a single tube. However, while this methodology was demonstrated to be accurate in a idealized system, the precision in porous media where the heterogeneity of the geometry can induce a complex flow field and therefore a complex deformation of the plume remained to be demonstrated.

To tackle this issue, in this study, we investigate an irreversible instantaneous chemical reaction of the form A + B → C at the pore-scale in two bead pack media characterized by two different degrees of heterogeneity. This simple bimolecular chemical reaction, which can be considered as constituent of more complex reactions, is encountered in many processes, such as the migration of radioactive materials [START_REF] Van Loon | Review of the kinetics of alkaline degradation of cellulose in view of its relevance for safety assessment of radioactive waste repositories[END_REF], metabolic activity of a biofilm [START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF] or ammonia [START_REF] Garg | Mass transfer with instantaneous chemical reaction in finite gas-liquid systems[END_REF]. We use a reactive particle tracking (RPT) model to simulate fluid-fluid reactive transport directly on the pore space of two different synthetic media (Perez et al., 2019a) to capture the quantity of reaction occurring. The methods proves itself to be accurate for complex geometries and its results serve as benchmark for the reaction predictions. We also investigate the lamella methodology [START_REF] Perez | Upscaling of mixing-limited bimolecular chemical reactions in Poiseuille flow[END_REF] for the computation of the interface width between the reactant. We validate this predictive methodology against the direct RPT simulation for heterogeneous flow field situations. These geometries, composed of beads, were selected in order to have configurations that are suitable for generalization. Therefore, even though these structures cannot be said to be representative of every medium, the results apply to more complex structures and to 3D geometries on which the methodology can be generalized easily.

The paper is organized as follows. In Section 2 we introduce the reactive transport problem. We detail how the flow is performed, and explain the different steps that occur during the RPT simulations. Section 3 first discusses the dispersive lamella approach for the prediction of the global mixing evolution in the system before comparing its predictions against the results of the RPT model for two different media. Finally, Section 4 presents the main conclusions of the work.

METHODOLOGY

Methodology

In this section we detail each of the steps of the reactive particle tracking algorithm. We discuss the generation of the geometry, the flow simulation, the streamline tracing algorithm, and the computation of the reaction.

Geometry Generation

We simulate the motion of the reactive species in two different 2-dimensional synthetic heterogeneous porous media. The synthetic porous media consist of a random packing of equally sized circular grains. The representation of granular obstacles in porous media using circular grains provide several advantages compared to real rock samples. It is possible, for example, to test simple geometries that make easier the generalization of results to multiple porous media whereas using experimentally acquired real images might provide representations with a high degree of arbitrariness. We generate the grains by selecting randomly their position from a uniform distribution. We avoid overlapping of grains by rejecting the position of a given grain if it overlaps a grain previously placed.

The algorithm stops when the target porosity is achieved or the maximum number of attempts to place a new grain is exceeded. The binary image of the geometry is composed of regular pixels that represent either void or solid. The mesh is created from a regular hexaedron mesh compatible with OpenFOAM and that exactly respects the geometry made by the pixels of the images. We refine the mesh by dividing each hexaedron by 2 in all directions. The final mesh cells have a size of 1.25 ×10 -6 m in all directions (∆x = ∆y).

This discretization level is selected such that the radius of a grain is divided in 35 cells.

We use this meshing technique to avoid any averaging or smoothing that often occurs in the course of the standard OpenFOAM meshing procedure [START_REF] Bibliography Gjetvaj | Dual control of flow field heterogeneity and immobile porosity on non-Fickian transport in Berea sandstone[END_REF] and to have an accurate representation of the solid beads. The dimensions of the first medium are L x × L y = 7.5 × 10 -3 m × 2 × 10 -3 m. The grain diameter d = 0.93 × 10 -4 m, the average size of pores L p = 3.0671 × 10 -5 m, and porosity φ = 0.5. The resulting discretization for the regular grid consists of 6022 × 1600 cells (corresponding to x and y dimensions respectively) with a porosity φ of 50%. The porosity φ is defined as the ratio between the number of pore cells to the total number of cells. The resulting geometry is displayed together with the flow field in the top panel of Figure 6.1.

The second medium is a 7.3 × 10 -3 by 2 × 10 -3 m rectangle for a resulting porosity of 50%. The average pore length L p is 2.625×10 -5 m. The resulting discretization for the regular grid consists of 5800 × 1600 grid (corresponding to x and y dimensions respectively). We then compare the reaction produced by the direct RPT simulations to the mass predicted using the dispersive lamella parameterized by the average medium porosity. To highlight the dependence of the prediction on the porosity experienced the plume, we also computed the dispersive lamella prediction parameterized by this time varying porosity through Equation (6.36). Figure 6.11 displays the evolution of the mass produced from the RPT and compare to the two predictions of the dispersive lamella: the first is parameterized by the porosity of the medium while the second depends on the porosity of the plume. The former underpredicts the mass created at early times since the porosity that the plume experiences at early times is larger the medium porosity used in Equation ( 6.36), it then slowly converges to the correct mass produced. The latter captures accurately the reaction behavior demonstrating the impact of the porosity on the mixing volume. We aim here at giving a warning for early predictions when the mixing volume porosity is not representative of the medium porosity. This can lead to an overprediction (or underprediction) of the mixing volume and thus result in an overprediction (or underprediction) of the mass produced. There are several situations where this can happen such as the experiment injection we discussed. This could for example happens in small media where the line injection is too small for the plume to experience the full medium heterogeneity at early times. This could also happen at the interface between two media where the porosity experienced by the plume evolves 6.4. SUMMARY AND CONCLUSIONS CHAPTER 7. CONCLUSIONS velocity distribution input for the continuous time random walk approach and second, allow for the validation of the model by comparing transport properties such as breakthrough curves.

• Pore-scale mechanisms such as clogging can bring temporal variations of the flow field. Introducing a temporally varying velocity distribution in CTRW type models could be a solution to account for this changes. A comprehensive analysis of the impact of clogging on the velocity field would be needed.

• The thesis focuses on the upscaling of transport from the pore-to the Darcy scale through the stochastic dynamics of the Lagrangian velocity of fluid particles. A next step might therefore be to investigate the stochastic dynamics of the velocity of fluid particles in three-dimensional Darcy scale media.

• Finally, here is a possibility that would combine every chapter of this thesis into a powerful predictive model. In Chapters 3, 4, and 5, we developed effective CTRW models for the prediction of transport properties, while in Chapter 6, we validated a model that predicts accurately the reaction mass from mixing volume measurements. If the CTRW could be improved to not only predict the displacement variance accurately, but also the effective variance (related to the mixing volume, see Chapter 6), this would translate directly into a predictive and effective CTRW model for the upscaling of chemical reaction. In this case (see Figure 8.1), since the eventual exit would be x 1 , we solve it the same way as before, replacing x(t) by x 1 in Equation (2.10) and solving for t In this case (see Figure 8.2), it is obvious that the particle will not be able to leave the cell on the x axis, that's why we set In that case the normalized KL distance reads as

where, in this scenario, p 0 1 denotes the initial distribution computed on the original volume V 0 , p 2 is our target distribution and p n 1 is the nth velocity distribution computed on the nth volume V n . CHAPTER 8. APPENDIX

Scientific Publications and Presentations

Scientific Publications

• Puyguiraud, A., [START_REF] Puyguiraud | Stochastic dynamics of Lagrangian pore-scale velocities in three-dimensional porous media[END_REF]. Stochastic dynamics of Lagrangian pore-scale velocities in three-dimensional porous media. Water Resources Research, 55.

• Puyguiraud, A., [START_REF] Puyguiraud | Stochastic dynamics of Lagrangian pore-scale velocities in three-dimensional porous media[END_REF]. Upscaling of anomalous pore-scale dispersion. Transport in Porous Media, 128(2), 837-855.

• Puyguiraud, A., Gouze, P. and Dentz. M. (submitted). Is there a representative elementary volume for anomalous dispersion?

• Puyguiraud, A., Perez, L. J., Hidalgo J. J. and Dentz. M. (in preparation). Prediction of pore-scale reaction in complex media.

• Perez, L. J., Puyguiraud, A., Hidalgo J. J., Jimenez-Martinez J., Gouze, P. and Dentz.

M. (in preparation)

. Impact of heterogeneous pore structure on mixing and chemical reactions.

Presentations at Congresses