University of Southern California Dissertations and Theses

On the transport dynamics of miscible fluids in porous media under different sources of disorder

(USC Thesis Other)

On the transport dynamics of miscible fluids in porous media under different sources of disorder

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content On the transport dynamics of miscible fluids in porous media under different sources of disorder by Alessandra Bonazzi A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ENVIRONMENTAL ENGINEERING) August 2023 Copyright 2023 Alessandra Bonazzi Table of Contents List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Chapter 2: Mixing in multidimensional porous media - A numerical study of the ef- fects of source configuration and heterogeneity . . . . . . . . . . . . . . . . . . . 5 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.1 Flow and transport model . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.2 Stochastic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.3 Mixing metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.4 Reference solutions for homogeneous media . . . . . . . . . . . . . . . . 12 2.2.4.1 Concentration PDF . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.4.2 Concentration moments . . . . . . . . . . . . . . . . . . . . . . 13 2.2.4.3 Dilution index . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5 Numerical set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.6.1 Concentration mean and variance . . . . . . . . . . . . . . . . . . . . . . 21 2.6.2 Concentration distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6.3 Dilution index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Chapter 3: Relative contributions of permeability heterogeneity and viscosity contrast on solute mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 ii 3.2 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.2 Mixing descriptors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Analysis on the sources of disorder in the flow field . . . . . . . . . . . . . . . . . 43 3.3.1 Viscosity contrast and heterogeneity effects on vorticity . . . . . . . . . . 43 3.3.2 Fingering mechanisms at R̸= 0 . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Numerical Model and Implementation . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Chapter 4: Influence of initial plume shape on miscible porous media flows under density and viscosity contrasts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 Mathematical description of flow and transport . . . . . . . . . . . . . . . . . . . 68 4.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2.2 Boundary and initial conditions . . . . . . . . . . . . . . . . . . . . . . . 70 4.3 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4 Comparison with previous work . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.5.1 Analysis of flow and transport behaviour . . . . . . . . . . . . . . . . . . 76 4.5.2 Mixing and spreading . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.5.2.1 Fluid interface dynamics . . . . . . . . . . . . . . . . . . . . . 81 4.5.2.2 Mixing analysis in the R− Ra parameter space and concentration statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.6 Effects of heterogeneity in the permeability field . . . . . . . . . . . . . . . . . . . 89 4.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Chapter 5: Uncertainty quantification in the presence of viscous fingering and perme- ability heterogeneity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.2.1 Multi-Gaussian Log-Permeability Field . . . . . . . . . . . . . . . . . . . 100 5.2.2 Mathematical Formulation of Flow and Transport . . . . . . . . . . . . . . 101 5.2.3 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2.4 Monte Carlo framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.3.1 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.3.2 Uncertainty in the concentration field . . . . . . . . . . . . . . . . . . . . 106 5.3.3 Uncertainty in early arrival times . . . . . . . . . . . . . . . . . . . . . . . 111 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Chapter 6: Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 iii Chapter 7: List of publications and papers under review . . . . . . . . . . . . . . . . 119 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 A Dimensionless equations for Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . 134 B Dimensional equations for Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . 135 C Details about the verification process in Chapter 4 . . . . . . . . . . . . . . . . . . 136 iv List of Tables 2.1 Input parameters used in the simulations . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 Skewness and kurtosis of the concentration empirical distribution and of the beta distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Relative errors between the fitted concentration CDFs and the empirical CDF . . . 32 3.1 Input parameters used in the simulations . . . . . . . . . . . . . . . . . . . . . . . 50 4.1 Model input parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.1 Input parameters used in the simulations . . . . . . . . . . . . . . . . . . . . . . . 104 v List of Figures 2.1 Comparison with the analytical and semi-analytical solutions . . . . . . . . . . . . 18 2.2 Summary of the considered scenarios . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Concentration field for d= 2 for point and line source . . . . . . . . . . . . . . . . 22 2.4 Concentration field for d= 3 for a vertical planar source . . . . . . . . . . . . . . 23 2.5 Mean concentration in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Concentration variance in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.7 Concentration PDF for distinct source configurations and flow dimensionality . . . 27 2.8 Comparison between concentration CDF and beta distribution model . . . . . . . . 29 2.9 Comparison between the concentration CDF, the Pareto distribution model and the uniform distribution model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.10 Dilution index in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.1 Initial distribution of the concentration field . . . . . . . . . . . . . . . . . . . . . 42 3.2 Concentration and logarithm of the local scalar dissipation rate . . . . . . . . . . . 45 3.3 Difference in the finger tip structures . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4 Single finger regions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5 Longitudinal and transversal velocity fields . . . . . . . . . . . . . . . . . . . . . 48 3.6 Concentration distribution for three different values of R . . . . . . . . . . . . . . 52 3.7 Breakthrough curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.8 Effect of R on first arrival times . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.9 Mean plume concentration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.10 Concentration variance of the solute cloud . . . . . . . . . . . . . . . . . . . . . . 56 vi 3.11 Cumulative density function of the plume concentration . . . . . . . . . . . . . . . 57 3.12 Mean scalar dissipation rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.13 Evolution of the contour length . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.14 Trends in time ofϕ andξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.1 Schematic configuration of the problem studied . . . . . . . . . . . . . . . . . . . 67 4.2 Comparison with previous work . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3 Concentration fields for the slower inlet flux case . . . . . . . . . . . . . . . . . . 75 4.4 Maps of the velocity magnitude for the slower inlet flux case . . . . . . . . . . . . 78 4.5 Concentration fields for the faster inlet flux case . . . . . . . . . . . . . . . . . . . 79 4.6 Maps of the velocity magnitude for the faster inlet flux case . . . . . . . . . . . . . 81 4.7 Scatterplots of the solute plume’sη andσ 2 c for f L = 0.001 . . . . . . . . . . . . . 82 4.8 Scatterplots of the solute plume’sη andσ 2 c for f L = 0.01 . . . . . . . . . . . . . . 82 4.9 Contour plots of the spatial statistics of the concentration for f L = 0.001 . . . . . . 86 4.10 Contour plots of the spatial statistics of the concentration for f L = 0.01 . . . . . . . 86 4.11 Concentration fields for a/b= 1/3 for three different levels of heterogeneity . . . . 91 4.12 Maps of the velocity magnitude for a/b= 1/3 for three levels of heterogeneity . . 92 4.13 Concentration fields for a/b= 3 for three levels of heterogeneity . . . . . . . . . . 93 4.14 Maps of the velocity magnitude for a/b= 3 for three levels of heterogeneity . . . . 94 5.1 Initial distribution of the concentration field . . . . . . . . . . . . . . . . . . . . . 103 5.2 Convergence analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.3 Maps of the concentration variance . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.4 Three-dimensional maps of the concentration variance . . . . . . . . . . . . . . . 108 5.5 Concentration variance (normalized with respect to the initial concentration c 0 ) across horizontal sections of the domain . . . . . . . . . . . . . . . . . . . . . . . 110 5.6 Probability Distribution Function of early arrival times . . . . . . . . . . . . . . . 113 5.7 Cumulative Distribution Function of early arrival times . . . . . . . . . . . . . . . 113 vii Abstract This thesis investigates the effect of multiple sources of disorder in solute transport in porous media through the use of numerical modeling. These sources are 1) the spatial heterogeneity of the permeability field and 2) the viscosity and density contrast between miscible fluids. Firstly, we focus on how degree of heterogeneity of a porous medium, flow dimensionality and source zone configurations impact mixing. Secondly, we examine how permeability heterogeneity and viscosity contrast impact the overall mixing behaviour of a solute. We analyze the degree and rate of mixing, spatial statistics of the concentration field and arrival times at a control plane to characterize spreading and mixing in the domain. Through our analysis, we are able to provide a quantitative separation of the impacts of fingering and heterogeneity. Furthermore, we parameterize the probability distribution function of the concentration field. Thirdly, we focus on viscosity and density contrasts as sources of disorder, and we investigate their impact on the temporal evolution of the spreading and mixing quantities. We show that such impact depends on the initial shape of the source distribution where the solute is injected and on the intensity of the horizontal background flux. We also investigate how a stratified permeability field can interact with these sources of disorder and affect the transport behaviour of the plume. Finally, we study how the uncertainty related to permeability heterogeneity and viscous insta- bilities affect uncertainties of quantities of interest such as early arrival times of a solute at a target location and the concentration distribution in space and time. viii Chapter 1 Introduction Groundwater is the main source of freshwater and accounts for a significant percentage of the water withdrawn for municipal, agricultural and industrial use (Dieter et al., 2018). Because of this, groundwater vulnerability to contamination is an important issue (Enzenhoefer et al., 2012). The very nature of groundwater, hidden from sight, makes it challenging for environmental engineers to monitor groundwater quantity and quality. Understanding the mixing process between the ambient fluid and a solute in a porous medium is of capital importance in both groundwater remediation and risk analysis; in the first case the solute is a remediation fluid that needs to fully mix with contaminated water, while in the second the solute is a contaminant that should not mix with the water eventually withdrawn by extraction wells at sensitive locations. In order to determine the degree of solute mixing in an aquifer we need to characterize flow and transport processes. This task, usually accomplished by using numerical models, is not trivial due to the spatial variability of subsurface’s hydraulic properties such as permeability (Rubin, 2003). Heterogeneity in the permeability field can lead to the formation of preferential paths for the solute plume (Rizzo & de Barros, 2017); permeability heterogeneity also impact the solute plume mixing and dispersive rate (Dentz & de Barros, 2015; Dentz et al., 2011; Le Borgne et al., 2015a; Dentz et al., 2000a; Kapoor & Kitanidis, 1998). In addition to the hydraulic features characterizing the porous formation, fluids properties such as viscosity and density can also affect both flow and transport processes; this aspect however has received less attention by the hydrogeological research community when compared to the medium 1 heterogeneity, and has mostly been studied in non-environmental applications like enhanced oil recovery (Lake, 1989) and CO 2 sequestration (Huppert & Neufeld, 2014). However, viscosity and density contrasts between injected and ambient fluid (i.e. contaminant and water or remediation fluid and contaminated water) are usually present. When the injected fluid is less viscous than the ambient fluid, the fluids’ interface is unstable and leads to a phenomenon known as viscous fingering. This interface instability, denominated as the Saffman-Taylor instability from the authors that firstly described it (Saffman & Taylor, 1958), occurs in both miscible and immiscible fluids. In miscible fluids, the stretching of the interface (i.e. the location where the concentration gradient is higher) caused by the fingers’ development can increase the mixing between two fluids (Jha et al., 2011a); as the fingering caused by dif- ferences in viscosity affects the evolution of the solute plum, the evolving concentration field in turn affects fluid’s properties that affect mobility, such as viscosity and density, since such prop- erties are dependent on concentration. On the other hand viscous fingering can cause channeling, a phenomenon in which the solute plume travels on a preferential path leaving most of the ambi- ent fluid unswept (Jha et al., 2011a), thus reducing mixing, an undesirable effect for instance in enhanced oil recovery applications (Orr & Taber, 1984) and in groundwater remediation. Viscous- fingering induced channeling can also lead to earlier arrival times of the solute plume to a given target location (Jha et al., 2011a). Density contrast between injected and ambient fluid can lead to another type of interface in- stability known as Rayleigh-Taylor instability, that occurs when the more dense fluid is above the less dense one. Some examples are the propagation of dense nonaqueous phase liquids (DNAPLs) plumes (Held & Illangasekare, 1995a) and CO 2 sequestration (Szulczewski et al., 2013; Jenny et al., 2014). This type of instability in porous media leads to the formation of gravity-driven fingers and the onset of a convection-dominated regime when the mixing process reaches its max- imum; as the concentration gradients between the fluids diminish, mixing slows down. This be- haviour has been observed for a homogeneous porous media (Hidalgo & Dentz, 2018a). 2 We can thus identify three main sources of disorder that impact mixing of fluids in a porous medium: permeability heterogeneity, viscosity contrast and density contrast. Comprehension of when and how the considered sources affect solute mixing and spreading is necessary in under- standing and, most importantly, controlling the evolution of a solute plume. For instance in ground- water remediation, where we wish for the injected fluid to sweep most of the domain and maximize mixing with the ambient fluid, our work could help identify the optimal viscosity and density of the remediation fluid and the best injection rate; in risk analysis, in case of aquifer contamination caused by a landfill leakage, a good understanding of how the contaminant plume will move and mix can help managers decide if and how long extraction wells will need to be shut down. In order to investigate the impact of these sources of disorder on mixing, one other feature of solute transport should be taken into account: the geometrical configuration of the solute source zone. In heterogeneous porous media, the source zone configuration has been shown to affect the transport behaviour of a solute plume (Dentz & de Barros, 2013a). In the presence of density and gravity fingering, we might expect one source of disorder being predominant over the other if the initial solute source is oriented perpendicularly to either the direction of gravity or the mean flow direction. The impact of the initial condition on the temporal evolution of mixing metrics has not received much attention despite its importance in applications such as groundwater remediation and CO 2 sequestration. We also perform a stochastic analysis in order to understand the interaction between two sources of disorder (viscosity contrast and permeability heterogeneity) and uncertainty related to lack of characterization of the structure of the permeability field, and the effect of such interac- tion on quantities of interest such as early arrival times at a target location and the concentration distribution. Several methods have been developed to quantify uncertainty of model predictions in groundwater flow and transport analysis Zhang et al. (2010). The most commonly used is the Monte Carlo approach due to its straightforward approach, that consists in running flow and trans- port simulations on a great number of permeability fields realizations due to our lack of knowledge of the actual permeability distribution in space Rubin (2003). 3 The objective of this thesis is to quantify when and how the above mentioned sources of disor- der impact the mixing dynamics of a solute plume. More precisely, we aim to identify conditions where one source of disorder overwhelms the others or when all sources are equally important. Improved understanding of the role of each of these factors, and their interactions, will enable engineers and geoscientists to design more efficient groundwater remediation strategies and to per- form more precise risk analysis assessments. We show how these processes affect several mixing metrics such as dissipation rates, mean and variance of the solute plume, early arrival times and concentration probability distribution functions. Moreover, we show how the contrast in fluid prop- erties (such as viscosity) can modulate the uncertainty in transport observables originating from the lack of characterization of the permeability field. The structure of this document is as follows. Chapter 2 provides a systematic analysis on how flow dimensionality, source zone configuration and degree of heterogeneity impact mixing in porous media in absence of viscosity and density contrasts between two miscible fluids. Chap- ter 3 provides a study on the impact of field heterogeneity and viscosity contrast on the trans- port behavior of an solute in a two-dimensional homogeneous porous medium. Through the use of high-resolution numerical simulations, the evolution in time of several mixing metrics is ana- lyzed. Chapter 4 studies the impact of density and viscosity contrasts on the temporal evolution of spreading and mixing properties, and shows that such impact depends on the initial geometrical configuration of the solute source. Chapter 5 analyzes how permeability heterogeneity and viscous fingering affect uncertainty related to early arrival times at a target location and the concentration distribution. Finally, Chapter 6 provides concluding remarks and an outlook to the possible future areas of research related to the topics here discussed. Chapter 7 provides a list of publications. 4 Chapter 2 Mixing in multidimensional porous media - A numerical study of the effects of source configuration and heterogeneity 2.1 Introduction The spatial variability of the hydraulic conductivity in porous formations leads to complex flow patterns which in turn lead to mixing dynamics that differ from the ones observed under uniform flow conditions. Transport in heterogeneous porous media is characterized by early solute break- through, tailing behavior of the concentration at late times and irregular spreading and mixing rates (Neuman & Tartakovsky, 2009; Dentz et al., 2011; Fiori et al., 2015). The presence of heterogene- ity is also responsible for poorly mixed conditions and spatially dispersed solute plumes (Dentz et al., 2011). The interplay between advective and local-scale dispersive fluxes in heterogeneous porous media results in non-trivial macroscopic transport behavior which is of importance for a broad range of scientific fields and engineering applications (Sahimi, 2011). Examples consist of contaminant migration in hydrology (Dagan & Neuman, 2005), aquifer remediation (Chapman & Parker, 2005), human health risk assessment (Im et al., 2020), geological storage of CO 2 (Hi- dalgo & Carrera, 2009) and safety assessment of waste repositories (Selroos, 1997). Improved understanding of the effects of heterogeneity on mixing is imperative to improve the predictive capabilities of models in the above-mentioned applications. 5 Many works focused on understanding the role of conductivity heterogeneity on temporal scal- ing properties of solute mixing. Mixing in heterogeneous porous media has been analyzed in terms of effective dispersion coefficients (Dentz et al., 2000a; Fiori & Dagan, 2000; Fiori, 2001), entropy-based mixing metrics (Kitanidis, 1994; de Barros et al., 2015) and the temporal decay of the concentration variability (Kapoor & Kitanidis, 1998). Through the use of numerical simu- lations, Le Borgne et al. (2010) showed how moderate-to-strong levels of heterogeneity induced anomalous temporal scaling for the scalar dissipation rate of a solute plume originating from a line source in a two-dimensional porous medium. Other studies investigated the impact of the sequence of fluid deformation events in two-dimensional flows on mixing metrics (de Barros et al., 2012a; Le Borgne et al., 2013, 2015b). The effects of flow focusing on the transverse dilution behavior of steady-state plumes have been topic of numerical (de Barros & Nowak, 2010; Cirpka et al., 2011) and experimental investigations (Rolle et al., 2009; Gueting & Englert, 2013). Kapoor & Kitanidis (1998) analyzed the rate of destruction of the concentration variance by local-scale dis- persion through the use of numerical simulations and approximate analytical solutions. Dentz et al. (2018) showed the importance of the initial condition of the solute plume (i.e. its spatial distribution) on the temporal mixing evolution in two-dimensional porous media. Semi-analytical solutions for the statistical description of the concentration field in spatially heterogeneous porous formations are also reported in the literature (Rubin et al., 1994; Fiori & Dagan, 2000; Fiori, 2001; Tonina & Bellin, 2008; Meyer et al., 2010; Dentz & Tartakovsky, 2010; Dentz, 2012; de Barros & Fiori, 2014). Approximate semi-analytical solutions for the dilution index, introduced by Kitanidis (1994), in two- and three-dimensional porous formations are reported in de Barros et al. (2015) and compared with results from numerical simulators (Boso et al., 2013a; de Barros et al., 2015) and field data (de Barros et al., 2015; Soltanian et al., 2020; de Barros & Fiori, 2021). In general, the above-mentioned semi-analytical solutions are based on perturbation theory, and therefore restricted to low-to-moderate levels of hydraulic conductivity heterogeneity. To ad- dress mixing in porous formations displaying a stronger degree of heterogeneity (i.e. logconduc- tivity variance larger than unity), numerical methods are needed. However, traditional Eulerian 6 grid-based numerical approaches are often plagued by oscillations and numerical dispersion which impact the accuracy of the numerical scheme (Zheng & Bennett, 2002; Ferziger et al., 2002a; Gotovac et al., 2007). Furthermore, in order to capture the effects of small scale heterogeneity on solute mixing, a fine spatial resolution in the numerical model is required which increases the computational burden associated with flow and transport simulations (Ababou et al., 1989b; Bellin et al., 1993). This computational burden is augmented within the context of uncertainty quantifica- tion where a Monte Carlo framework is needed (Moslehi et al., 2015) . Lagrangian-based methods, such as the Random Walk Particle Tracking (RWPT) technique (Salamon et al., 2006), are an ap- pealing alternative since they are globally mass conservative and no subject to artificial oscillation and numerical dispersion. RWPT was used to study dispersion in both two- and three-dimensional heterogeneous porous media (Bellin et al., 1992a; de Dreuzy et al., 2007; Beaudoin & de Dreuzy, 2013). Jankovic et al. (2017) used RWPT to examine the impact of conductivity structure in the mass breakthrough curve in a three-dimensional setting. Sole-Mari et al. (2021), used RWPT sim- ulations to key metrics of transport originating from a source zone occupying the entire inlet area of a bounded three-dimensional porous media displaying non-Gaussian features and a logconduc- tivity variance of unity. Libera et al. (2019) used RWPT to analyze the joint effects of porosity and conductivity variability on both the peak flux-averaged concentration and solute arrival times in three-dimensional porous media for low and high levels of heterogeneity. Using different numer- ical schemes, including RWPT and smoothed particle hydrodynamics (SPH), Boso et al. (2013b) investigated the impact of heterogeneity in transport metrics such as the second central spatial moments of the solute plume and the dilution index. Results obtained from different numerical schemes were compared and reported for a point-like injection and two-dimensional domains for different logconductivity variances (ranging from 0.2 to 10). Despite the benefits associate with RWPT, there are drawbacks such as the presence of local concentration fluctuations and the need to have a significant number of particles to achieve numerical precision of the concentration field, es- pecially when dealing with heterogeneous porous media (Herrera et al., 2009; Boso et al., 2013b). 7 The need to use a large number of particles has implications on the computational costs associated with transport. To overcome this challenge and improve the efficiency associated with Lagrangian-based meth- ods, Rizzo et al. (2019) introduced a GPU-accelerated RWPT, denoted as PAR 2 , that 1) enables the use of a large number of particles and 2) is computationally efficient. PAR 2 has been employed to study hydrogeological connectivity in both Gaussian and non-Gaussian flow fields (Rizzo & de Barros, 2019; Morvillo et al., 2021a), aquifer resilience loss and probabilistic risk analysis (Morvillo et al., 2022), solute transport at the pore-scale (Kamrava et al., 2021) and has been ex- panded to account for chemical reactions (Morvillo et al., 2021b). The code is open source (see details in Rizzo et al. (2019)) and a step-by-step tutorial on its use and how to connect to existing groundwater flow simulation tools can be found in Morvillo et al. (2022). The objective of this Chapter is to provide a systematic numerical investigation of the impact of flow dimensionality, source zone configuration and the degree of heterogeneity on metrics of solute mixing. To achieve our goals, we rely on the RWPT-based simulator PAR 2 (Rizzo et al., 2019) to show the importance of the aforementioned factors in mixing metrics such as the global spatial mean and variance of the concentration field, the dilution index and the concentration probability distribution function. Results are compared to existing analytical solutions for homogeneous and heterogeneous porous media. 2.2 Problem statement 2.2.1 Flow and transport model We consider a d-dimensional porous medium with constant porosityφ and spatially variable (lo- cally isotropic) hydraulic conductivity K(x) with x=[x 1 ,...,x d ] T denoting the Cartesian coordinate system and d= 2 and 3. The flow field is given by 8 ∇· q(x) = 0 (2.1) where the specific discharge q is obtained through Darcy’s law q(x)=− K(x)∇h(x), (2.2) with h representing the hydraulic head. Equation 2.1 is subject to the following boundary condi- tions: prescribed hydraulic head along the longitudinal direction x 1 and no-flow boundary condi- tions along directions x j for j= 2 and 3. These conditions ensures that flow is uniform-in-the-mean along the longitudinal direction. That implies that the mean velocity field is ⟨u⟩=(U,0) for d= 2 and⟨u⟩=(U,0,0) for d= 3. Here the angled bracket represent the ensemble expected value and U denotes the mean longitudinal velocity given by U = K G J/φ where K G is the geometric mean of the hydraulic conductivity andJ is the spatially uniform mean hydraulic gradient. An inert solute is instantaneously released in a domainV o located within the porous formation. Depending on the dimensionality of the flow field, the domain V o can represent a volume or an area. In this Chapter, we will consider different geometrical configurations for the source zone. The geometrical configurations considered are point, line and planar sources in both two- and three dimensional porous formations. Details pertaining the dimensions of the solute injection zone are provided further below in Section 2.5. The spatiotemporal dynamics of the solute plume’s resident concentration c is assumed to be governed by the local advection dispersion equation (ADE): ∂c(x,t) ∂t + u(x)· ∇c(x,t)=∇· [D(x)∇c(x,t)], (2.3) where u= q/φ is the velocity vector and D is the local scale dispersion tensor given by D(x)=(α T |u(x)|+ D m )1+ (α L − α T ) |u(x)| u(x)u(x) T (2.4) 9 with D m denoting the molecular diffusion coefficient and 1 representing the identity matrix. Equa- tion 2.3 is subject to the following initial condition c(x,0)=        c o if x∈V o ; 0 otherwise. (2.5) Due to the spatial randomness of the hydraulic conductivity, the velocity field is spatially vari- able which impacts solute mixing and spreading rates. In the following subsection, we will describe the details regarding the spatial structure of the K field. 2.2.2 Stochastic model To simulate flow and transport in a spatially random porous medium, we assume that the logcon- ductivity field, namely f(x)= lnK(x), is multivariate Gaussian and statistically stationary. There- fore, f is fully characterized by i) its mean value⟨ f⟩= lnK G where K G is the geometric mean of the conductivity field K G and ii) its spatial covarianceC ff (r)=⟨ f(x) f(x ′ )⟩ with r= x− x ′ . The logconductivity variance is given by σ 2 f ≡ C ff (0). In this Chapter, we consider an isotropic cor- relation lengthλ. In the following, we adopt an exponential spatial covariance model forC ff (see Chapter 2 of Rubin, 2003)) such that C ff (r)=σ 2 f exp − |r| λ . (2.6) 2.2.3 Mixing metrics To quantify mixing, we will rely on four descriptors. The first descriptor is the spatially averaged concentration of the plume over volumeϕ d of the d-dimensional domain characterized by c> c ∗ 10 where c ∗ is a low concentration threshold value. The volumeϕ d (t|c ∗ ) represents a line in d = 1, an area for d= 2 and volume for d= 3 and is defined as follows ϕ d (t|c ∗ )= Z Ω H[c(x,t)− c ∗ ]dx (2.7) where H[·] is the Heaviside function andΩ is the flow domain. This means, ϕ d (t|c ∗ ) denotes the volume occupied by the solute, or, in other words, the mixing volume. The spatially averaged concentration is defined by µ c (t)= 1 ϕ d (t|c ∗ ) Z ϕ d (t|c ∗ ) c(x,t)dx. (2.8) Note that for a sufficiently small c ∗ , the integral on the right side may be approximated by 1 due to mass conservation. Thus, the spatial mean concentration is inversely proportional to the mixing volumeϕ d (t|c ∗ ). The second descriptor is the spatial variance of the concentration σ 2 c (t)= 1 ϕ d (t|c ∗ ) Z ϕ d (t|c ∗ ) c(x,t) 2 dx− [µ c (t)] 2 , (2.9) withµ c given in Equation 2.8. We also examine the concentration probability density function (PDF) p(c). The PDF p(c) is obtained from spatially sampling the concentration point values within the domain ϕ d (t|c ∗ ), see Equation 2.7. It can be formally written as p(c)= Z ϕ d (t|c ∗ ) δ[c− c(x,t)]dx, (2.10) whereδ(c) is the Dirac delta distribution. Finally, the fourth mixing metric investigated in this Chapter is the dilution index E introduced by Kitanidis (1994). The dilution index represents a global measure of dilution. It measures the temporal evolution of the volume occupied by the solute plume and allows to quantify the 11 combined effects of advection and local scale dispersive fluxes on mixing. The dilution index is mathematically expressed as E(t) = exp[Λ(t)], with Λ(t) = − Z Ω c(x,t) M o ln c(x,t) M o dx (2.11) where M o = c o V o φ is the total mass injected into the porous formation. We will computeµ c ,σ 2 c , p and E for different values ofσ 2 f , distinct source dimensions (point, line and planar source zones) and for different flow dimensionality ( d= 2 and 3). 2.2.4 Reference solutions for homogeneous media For reference, we give in the following the expected behaviors of these mixing metrics for homoge- neous media in d spatial dimensions. The concentration distribution in response to an instantaneous point source, c(x,t = 0)= c 0 δ(x), in an infinite ν-dimensional domain is given by the Gaussian distribution c(x,t)= c 0 (4πDt) ν/2 exp − (x− ut) 2 4Dt . (2.12) Note that solution for an infinitely extended line source in d = 2 spatial dimensions and for an infinitely extended planar source in d = 3 is given by expression 2.12 for ν = 1, the solution for an infinitely extended line source in d= 3 spatial dimensions is given by expression 2.12 for ν= 2. In general,ν= d− d s , where d is the dimension of space and d s the dimension of the source distribution. The source dimension is d s = 1 for a line source and d s = 2 for a planar source. In the following we rescale c→ c/c 0 . 12 2.2.4.1 Concentration PDF In order to determine the concentration PDF we note that the radius of a concentration isoline measured from the center of mass at ut is related to the concentration c by r(c)= p (4Dt)ln[c m (t)/c], c m (t)= 1 (4πDt) ν/2 . (2.13) Thus, we can write Equation 2.10 for the concentration PDF in spherical coordinates as p(c)= d r(c ∗ ) r(c ∗ ) Z 0 drr ν− 1 δ[c− f(r)], f(r)= exp(− r 2 /4Dt) (4πDt) ν/2 . (2.14) We can evaluate this integral explicitly by using thatδ[1− f(r)]= 1/|d f/dr|δ[r− r(c)]. Thus, we obtain p(c)= ν 2c ln[c m (t)/c] (ν− 2)/2 ln[c m (t)/c ∗ ] ν/2 I(c ∗ < c≤ c m (t)). (2.15) 2.2.4.2 Concentration moments The concentration moments are defined in terms of the concentration PDF by µ (k) (t)= ∞ Z 0 dcc k p(c). (2.16) Using the explicit expression 2.15 for the concentration PDF and rescaling the integration variable by c m (t), we obtain µ (k) (t)= c m (t) k ln[c m (t)/c ∗ ] ν/2 ν 2 1 Z c ∗ /c m (t) dcc k− 1 ln(1/c) (ν− 2)/2 . (2.17) 13 We consider situations for which (c ∗ /c m (t)≪ 1). In this limit we obtain µ (k) (t)= A k c m (t) k +..., (2.18) where the dots denote subleading contributions. The constant A k is defined by A k = ν 2ln(1/c ∗ ) ν/2 1 Z 0 dcc k− 1 ln(1/c) (ν− 2)/2 . (2.19) It depends on the dimension ν of space and the order k of the moment. We denote the mean concentration byµ (1) ≡ µ c . We obtain for the mean concentration in leading order µ c (t)= A 1 c m (t)∼ t − ν/2 , (2.20) whereµ c (t)≡ µ (1) (t). For the concentration variance, we obtain σ 2 c (t)= c m (t) 2 (A 2 − A 2 1 )∼ t − ν . (2.21) 2.2.4.3 Dilution index The expression for the dilution index is obtained by inserting the Gaussian solution 2.12 into the definition 2.11. The explicit analytical solutions are Kitanidis (1994) E(t)= exp(d/2) c m (t) ∼ t ν/2 . (2.22) Note that is scales as 1/µ c (t). This is due to the fact that the mean concentration scales as the inverse mixing volume, while the dilution index describes the evolution of the mixing volume. 14 2.3 Methods For all results presented in Section 2.6, the spatially variable logconductivity field f is randomly generated using the SGeMS tool (Remy et al., 2009). SGeMs is based on a sequential Gaussian simulation model (Rubin, 2003). The flow equation 2.1 is solved using the finite-difference solver MODFLOW (Harbaugh, 2005) with the Python package FloPy (Bakker et al., 2016). The flow domain, L 1 × L 2 (for d = 2) or L 1 × L 2 × L 3 (for d = 3), is discretized into a regular grid with numerical grid block of dimension’s∆. To simulate solute transport (Equation 2.3), we make use of the PAR 2 code developed by Rizzo et al. (2019). PAR 2 is an open source, GPU-accelerated, simulator based on the Random Walk Particle Tracking (RWPT) (Rizzo et al., 2019). A step-by-step tutorial of PAR 2 , and the details of how to link it to an open source Python package, is reported in Morvillo et al. (2022). PAR 2 has been employed to study transport at both field and pore scales (Rizzo et al., 2019; Kamrava et al., 2021). The RWPT is based on the trajectory of the i th solute particle that can be computed using the Itˆ o –Taylor integration scheme (Salamon et al., 2006): X i (t+∆t)= X i (t)+ A(X i (t))∆t+ B(X i (t))· ξ ξ ξ(t) √ ∆t (2.23) whereξ ξ ξ is a normally distributed random variable with zero mean and unit standard deviation,∆t is the time step, and the drift vector A and the displacement matrix B are defined by: A(x) = u(x)+∇· D(x)+ 1 φ D(x)· ∇φ (2.24) 2D(x) = B(x)· B(x) T . (2.25) 15 2.4 Validation In this section we validate the results from the numerical simulators with existing analytical so- lutions for two-dimensional domains. For a homogeneous conductivity field, Kitanidis (1994) derived the following expression for the dilution index for a point source (see also Equation 2.22): E(t)=(4πDt) d/2 exp d 2 (2.26) For a heterogeneous porous medium, de Barros et al. (2015) used perturbation theory to derive a semi-analytical solution for the dilution index. The semi-analytical expression is valid for a point-like instantaneous injection, uniform-in-the-mean flow conditions in the absence of sinks and sources and low-to-mild levels of heterogeneity, i.e. σ 2 f ≲ 1. The semi-analytical solution for the temporal evolution of the dilution index E is given by E(t)=(2π) d/2 exp d 2 d ∏ i=1 p W ii (t) (2.27) whereW ii is the relative particle trajectory covariance along the i th direction given by Fiori (2001) W ii (t)=X ii (t)+ 2Dt− Z ii (t;0), (2.28) whereX ii andZ ii correspond to the one- and two-particle trajectory covariances respectively (Fiori & Dagan, 2000). Note that the “0” present in the left-hand side of the two-particle trajectory covarianceZ ii is to emphasize that the function is evaluated for a point-like source, i.e., the sep- aration distance between two particles originally located within the source zone is approximately zero (details provided in de Barros et al., 2015). The semi-analytical expressions for the parti- cle trajectory covariancesX ii andZ ii were derived in Fiori & Dagan (2000) and are reproduced in Equations 2.29 and 2.30. These expressions are valid for low-to-mild levels of heterogeneity 16 and uniform-in-the-mean flow conditions. For an isotropic and constant local dispersion D and point-like injection, the particle trajectories are X ii (t)= 2 1+d/2 π d/2 Z ∞ 0 Z t 0 (t− τ)cos[k 1 Uτ]e − Dk 2 τ ˆ u ii (k)dτdk, (2.29) Z ii (t;0)= 2 π d/2 Z ∞ 0 Z t 0 Z t 0 cos[k 1 U(t ′ − t ′′ )]e − Dk 2 (t ′ +t ′′ ) ˆ u ii (k)dt ′ dt ′′ dk, (2.30) where mean longitudinal velocity is denoted by U and ˆ u ii (k) corresponds to the Eulerian velocity covariance in Fourier space. The wave number vector is given by k and k 2 =∑ ∞ j=1 k 2 j for j= 1, ..., d. Equations 2.29 and 2.30 are limited to small to mild levels of heterogeneity. The expression for the velocity covariance in Fourier space, i.e. ˆ u ii (k), is (Dagan, 1984): ˆ u ii (k)= U 2 δ 1i − k 1 k i k 2 δ 1i − k 1 k j k 2 ˆ C ff (k), (2.31) with andδ i j is Kronecker’s delta and ˆ C ff (k) is the Fourier transform of the logconductivity covari- ance, Equation 2.6. For the purpose of validation, we will restrict ourselves to a two-dimensional setting. Therefore, the Fourier transform of the isotropic logconductivity exponential covariance function for d= 2 is (Rubin, 2003): ˆ C ff (k)=σ 2 f λ 2 (1+ k 2 1 λ 2 + k 2 2 λ 2 ) − 3/2 . (2.32) Figure 2.1 (left) shows how the flow and transport simulator was able to capture the dilution of a plume in a two-dimensional homogeneous porous medium. The numerical results shows an excellent agreement with the analytical solution, see Equation 2.26. The parameters used for the simulations in the homogeneous setting are reported in Table 2.1. Next, we compare the numerical results with the semi-analytical solutions derived in de Barros et al. (2015) in a heterogeneous setting, see Equation 2.27. The results reported in Figure 2.1 (right) are for a two-dimensional heterogeneous aquifer characterized by σ 2 f = 1. For this comparison, 17 Figure 2.1: Comparison between the results from the numerical flow and transport simulators (red crosses) against the analytical results (continuous line) derived by Kitanidis (1994) for a homo- geneous medium (left) and the semi-analytical solution derived by de Barros et al. (2015) for a heterogeneous medium characterized by σ 2 f = 1 (right). Results depicted for d = 2 and a point source. The advective time scale for the homogeneous medium (left) is given byτ a = L 1 /U where L 1 is the longitudinal dimension of the domain. The advective time scale for the heterogeneous medium (right) isτ a =λ/U withλ as the correlation scale. the spatially random porous medium was generated using HYDRO GEN (Bellin & Rubin, 1996). Given the stochastic nature of the flow field, we averaged the temporal evolution of the dilution index over five random realizations. The results depicted in Figure 2.1 (right) show an excellent agreement with the numerical simulator and the semi-analytical results. This is remarkable given that the perturbation theory based solution (Equation 2.27) is expected to deteriorate forσ 2 f ≳ 1, 2.5 Numerical set-up A description of the scenarios investigated is displayed in Figure 2.2. As shown in Figure 2.2, we will consider both two- and three dimensional flows characterized by two distinct logconductivity variances, namelyσ 2 f = 1 and 4. These values were selected to represent a mildly heterogeneous case and a highly heterogeneous case. In order to capture the effects of heterogeneity of the flow field on mixing, we consider a ratio of λ/∆= 10. The ratioλ/∆ was selected based on the analysis 18 carried out in the literature (Ababou et al., 1989b; Bellin et al., 1993, 1994; de Dreuzy et al., 2007; Moslehi et al., 2015). A detailed numerical error analysis regarding the spatial resolution of the flow field on simulations is reported in figures 5-6 of Leube et al. (2013) and figures 2-6 of Moslehi et al. (2015). Figure 2.2 also provides a graphical representation of the different geometrical configurations for the source zone explored in our study. Details regarding the different scenarios are provide below. Two-dimensional case: For d = 2, we will consider a point source of dimensions 0.1λ× 0.1λ centered at x=(24,100). We will also consider a line source of length 3λ located at x 1 = 24 with transverse dimension x 2 ∈[85,115]. Three-dimensional case: For d= 3, we will consider, analogously to the case for d= 2, a point source of dimensions 0.1λ× 0.1λ centered at x=(24,100,100) and also a line source located at x 1 = 24 and x 3 = 100 with transverse dimension x 2 ∈ [85,115]. We will also consider a plane source located at x 1 = 24 with transverse dimension x 2 ∈ [85,115] and x 3 ∈ [85,115] (i.e., of dimensions 3λ× 3λ). All other parameter values used in the simulations presented in this Chapter are listed in Table 2.1. All the simulations were performed on a desktop computer equipped with a GPU NVIDIA GeForce GTX 745/PCIe/SSE2 necessary to use PAR 2 and with 4GB of internal RAM memory. The time required to complete a high-resolution simulation (i.e.,∆/λ = 0.1) with 10 7 particles ranged from about three hours for the two-dimensional cases to up to six hours for the three-dimensional cases. 2.6 Computational results Figures 2.3 and 2.4 provide examples of the simulated concentration fields in two- and three- dimensional flows for σ 2 f = 1 and 4 at the final simulation time. The spatial concentration plume 19 2 f 2 = 1 f 2 = 4 3 Type of source 2 3 2 3 2 3 1 2 1 2 3 / 3 / 2 / 1 / 2 / 1 / POINT POINT LINE LINE PLANE at 1 / = 2 . 4 at 1 / = 2 . 4 at 1 / = 2 . 4 Figure 2.2: Summary of the considered scenarios. We analyze transport in domains with d = 2 and d= 3 characterized by logconductivity fields with variance σ 2 f = 1 andσ 2 f = 4. In each field we study the impact on solute transport of the source dimension. 20 Table 2.1: Input parameters used in the simulations Parameter Symbol Value Calculated as Correlation length in x 1 and x 2 λ 10 m - Correlation length in x 3 (d= 3 only) λ 10 m - Domain length in x 1 L 1 1000 m 100λ Domain length in x 2 L 2 200 m 20λ Domain length in x 3 (d= 3 only) L 3 200 m 20λ Mesh size in x 1 and x 2 ∆ 1 m λ/10 Mesh size in x 3 (d= 3 only) ∆ 1 m λ/10 Porosity φ 0.25 - Mean hydraulic conductivity K G 1 m/d - Mean longitudinal flow velocity U 0.02 m/d - Dispersivity in x 1 α L 0.1 m 2 0.1∆ Dispersivity in x 2 (except Section 2.4) α T 0.01 m 2 0.01∆ Dispersivity in x 2 (Section 2.4 only) α T 0.1 m 2 0.1∆ Dispersivity in x 3 (d= 3 only) α T 0.01 m 2 0.01∆ Molecular diffusion D m 8.64× 10 − 5 m 2 /d - Number of particles N p 10 7 - shows more complex features with increasing σ 2 f that are accompanied with overall lower con- centration values, which indicates increased dilution. In the following, we study the mixing and dilution dynamics in a more quantitative fashion in terms of the evolution of the mean concentra- tions, concentration variance, concentration PDF and dilution index introduced in Section 2.2.3. 2.6.1 Concentration mean and variance Figures 2.5 and 2.6 report the temporal evolution of the spatial concentration mean and variance defined in Equations 2.8 and 2.9 for d= 2 and d= 3 spatial dimensions for point, line and planar source distributions. The scalings derived in Section 2.2.4 for homogeneous media indicate the dependence of the decay of the concentration mean and variance on the dimensions of space and the source distribution. We found that the mean concentration scales as µ c (t)∼ t − (d− d s )/2 , and the concentration variance asσ 2 c (t)∼ t − (d− d s ) , where d is the space dimension and d s the source dimension. This means, that the concentration mean and variance for a homogeneous medium decay faster for higher spatial dimensions and lower source dimension. We observe qualitatively 21 Point source Line source σ 2 f = 1 σ 2 f = 4 Figure 2.3: Concentration field for d = 2 at the final simulation time t/τ a = 20 for point (left column) and line source (right column), in the conductivity field with σ 2 f = 1 (top row) andσ 2 f = 4 (bottom row). In all four images, the initial position of the solute source is depicted in orange color. similar dependence on space and source dimension for heterogeneous media, while the decay is in general faster for the heterogeneous than for the homogeneous media. This effect is increasing with increasingσ 2 f . Forσ 2 f = 1, we find that the concentration mean and variance for point and line sources in d= 2 and 3 decay only slightly faster than expected for homogeneous media. Forσ 2 f = 4, in contrast we observe a strong acceleration of the decay of concentration mean and variance compared to homogeneous media. For the planar source in d= 3 dimensions, concentration mean and variance decay significantly faster than for a homogeneous medium both for σ 2 f = 1 and 4. Such faster decay can be explained by the fact that a planar source is able sample more fluctuations of the permeability field. This means that higher concentration gradients have the opportunity to be generated, thus leading to faster mixing. As illustrated in Figures 2.3 and 2.4, the concentration distributions are strongly dispersed for both heterogeneity strengths, and the initial distribution evolves into a lamellar structure. This means, spatial heterogeneity creates strong concentration gradients, which are efficiently degraded through local scale dispersion, thus explaining the faster 22 f 2 = 4 f 2 = 1 Figure 2.4: Concentration fields at the final simulation time t/τ a = 20 for d= 3 withσ 2 f = 1 (top box) andσ 2 f = 4 (bottom box) for a vertical planar source. 23 decay of mean concentration and concentration variance in heterogeneous media compared with homogeneous media. Figure 2.5: Mean concentration in time for: point source (top left), line source (top right), 3d field withσ 2 f = 1 (bottom left), and 3d field with σ 2 f = 4 (bottom right). 2.6.2 Concentration distribution Figure 2.7 illustrates the concentration PDF, p(c) for point and line sources in both d = 2 and 3 spatial dimensions at different times. As discussed in the previous section, at increasing times the concentration mean and variance decrease. The maximum observable concentrations are shifted 24 Figure 2.6: Concentration variance in time for: point source (top left), line source (top right), three-dimensional field with σ 2 f = 1 (bottom left), and three-dimensional field with σ 2 f = 4 (bottom right). 25 towards smaller values. At small concentration values, p(c) behaves as predicted for a homoge- neous medium as p(c)∼ c − 1 for all cases under consideration. The solutions for a homogeneous medium are given by Equation 2.15. At increasing concentration values, the behaviors expected for a homogeneous medium are different from the observations for the heterogeneous media. For the point source in d= 2 and the line source in d= 3, we expect p(c)∼ c − 1 for c< c m (t). For the point source in d= 3, the concentration PDF behaves as p(c)∼ c(lnc) 1/2 and for the line source in d= 2 as p(c)∼ c(lnc) − 1/2 . The concentration PDF decreases faster towards smaller maximum values than the corresponding homogeneous solutions. These behaviors can be understood in the lamellar mixing framework of Villermaux (2012) and Le Borgne et al. (2015b) as the result of stretching-enhanced local dispersion. As shown in Figures 2.3 and 2.4, in heterogeneous media, the initially homogeneous solute distribution evolves into a lamellar structure due to the spatial variability of the underlying flow field. Across individual lamellae, the concentration distribution can be approximated by Gaussian distributions character- ized by a maximum distribution that depends on the lamella’s stretching history. For a lamella that has passed through regions of strong stretching, the maximum concentration decreases much faster than for the lamella that remains undeformed, as long as the values of local dispersion large enough to efficiently dissipate concentration gradients. The concentration PDF for an individual lamella is then characterized by 2.15 for a Gaussian concentration distribution. The full concentration PDF is then constructed by the superposition of lamellae with different stretching histories. This explains the behavior of p(c)∼ c − 1 at small concentration values, which is characteristic of the Gaussian concentration profile across lamellae, and the decay behavior at large concentration val- ues, which reflects the stretching and deformation history of the lamellae. In this sense, the mixing behavior at the simulation times under consideration here, which are smaller than the characteristic dispersion time over a correlation length, can be understood by the superposition of independent stretched lamellae. At times larger than the characteristic dispersion time scale, lamellae start co- alescing, which leads to different behaviors of the concentration PDF as discussed in Le Borgne et al. (2015b). However, this regime is beyond the scope of this paper. 26 Figure 2.7: Concentration PDF for distinct source configurations and flow dimensionality. The results in this Figure are for the fields with σ 2 f = 1. Results reported for a point source for d = 2 (top left), point source for d = 3 (top right), line source for d = 2 (bottom left), and line source for d= 3 (bottom right). The dotted line represents the concentration PDF in a homogeneous field according to Equation 2.15. 27 Next, we compare the cumulative distribution function (CDF), namely P(c) = R c 0 p(ψ)dψ, obtained from the numerical simulations with a parametric CDF model, the beta CDF. The beta CDF was originally suggested by Fiori (2001) as a potential distribution to capture the statistics of the concentration field given its flexibility to shift from a bimodal to a unimodal PDF. It has been widely employed to quantify the distribution of the concentration values in heterogeneous domains (Fiorotto & Caroni, 2002; Bellin & Tonina, 2007) and tested against numerical and analytical solutions (Schwede et al., 2008; Boso et al., 2013a; de Barros & Fiori, 2014; Boso & Tartakovsky, 2016). The performance of the beta distribution model was also assessed for mixing of fluids with distinct viscosities (Bonazzi et al., 2021) as well as mixing in non-Gaussian random flow fields (de Barros et al., 2022). The beta CDF model is defined by: P(c)= Γ[a+ b] Γ[a]Γ[b] Z c 0 w a− 1 (1− w) b− 1 dw, (2.33) whereΓ[z] is the Gamma function. The exponents a and b are given in terms of the empirically determined concentration meanµ c and varianceσ 2 c as a= µ c (t) β , b= 1− µ c (t) β β = σ 2 c (t) µ c (t)[1− µ c (t)]− σ 2 c (t) , (2.34) The results displayed in Figure 2.8 show the comparison between the CDF obtained from the numerical simulations and the beta CDF model 2.34 fitted with Equations 2.8 and 2.9, i.e., we calculated the spatial average and spatial variance of the concentration in the mixing volume obtained from the numerical simulations and used those values to estimate the parameters a and b in Equation 2.34. Figure 2.8 illustrates the comparison in both two-dimensional (Figure 2.8, left) and three-dimensional (Figure 2.8, right) domains for distinct source configurations at t/τ a = 12.3. In agreement with previous findings reported in the literature (Fiorotto & Caroni, 2002; Boso et al., 2013a; de Barros & Fiori, 2014), the beta distribution fails to capture the probabilities at low concentrations. This is particularly evident in the case of two-dimensional domains (see Figure 2.8, left). The fitting of the beta distribution improves with an increase of the dimensionality of 28 the porous medium (compare Figure 2.8 left and right). As depicted in Figure 2.8 (right), the agreement between the beta CDF model and the empirical CDF is augmented when the solute source zone is a vertical plane. The discrepancy between the beta model and the data at small concentrations can be traced back to the fact that the mixing behavior at the small concentration fringes of the solute plume is dominated by diffusion. We further explore this discrepancy at low concentrations by analyzing the third (skewness) and fourth (kurtosis) statistical moments of the concentration data at t/τ a = 12.3 for a log-conductivity variance of unity. Table 2.2 provides a comparison between the higher order moments obtained from the fitted beta CDF and the empirical CDF obtained from the high resolu- tion numerical simulations. The higher moments of the beta distributions were calculated from the distribution coefficients in Equation 2.34. As expected, the higher order moments estimated from the fitted beta model are different than the ones characterizing the empirical distribution. Figure 2.8: Comparison between the concentration cumulative distribution function (CDF) and beta distribution model for d = 2 (left) and d = 3 (right). The results in this Figure are for the fields with σ 2 f = 1 at time t/τ a = 12.3. Results computed for different source zone configurations (point, line and vertical planar sources). To overcome the limitations of the beta distribution at capturing the probabilities of low con- centration occurrences, we explore alternative distributions. To do so, we propose to employ two one-to-one variable transforms, namely the power law and logarithm transforms. We transform the 29 Table 2.2: Skewness and kurtosis of the concentration empirical distribution (originated from the numerical simulations) and of the fitted beta distribution at t/τ a = 12.3. The ratios of the higher order moments between the empirical and the beta distributions are also reported. The skewness and kurtosis of the empirical distribution obtained from the numerical simulations are represented byµ (3) andµ (4) whereas the ones for the beta distribution areµ (3) β andµ (4) β . Skewness Kurtosis Simulation µ (3) µ (3) β µ (3) /µ (3) β µ (4) µ (4) β µ (4) /µ (4) β d= 2, point 3.05 4.55 0.67 12.54 30.94 0.41 d= 2, line 1.82 2.92 0.62 5.82 12.58 0.46 d= 3, point 4.15 5.14 0.81 22.94 39.71 0.58 d= 3, line 3.79 4.23 0.90 20.17 26.71 0.76 d= 3, plane 2.28 2.72 0.84 8.85 10.98 0.81 concentration data originated from the numerical simulations as follows: √ c and ln[c]. The results depicted in Figure 2.9 show a good performance, within the low concentration range, of the Pareto type IV distribution for the √ c and a uniform distribution for ln[c]. Table 2.3 reports the relative error between the estimates for the probability that c< c ∗ obtained from the fitted concentration CDFs and the empirical CDF. The relative error is given by ε = 100×| P(c ∗ )− b P(c ∗ )|/P(c ∗ ), where P(c) is the empirical concentration CDF obtained from the numerical simulations and b P(c) is the fitted CDF (i.e. beta, uniform or Pareto type IV). The concentration threshold c ∗ is defined as follows: c ∗ (t)= c min (t)+η[c max (t)− c min (t)], (2.35) where c max and c min are the maximum and minimum concentrations within in the plume at time t. The quantities are defined as: c min (t) = min x∈ϕ d c(x,t); (2.36) c max (t) = max x∈ϕ d c(x,t), (2.37) 30 whereϕ d is defined in Equation 2.7. For the results present in Table 2.3, we compute c ∗ forη = 0.01, 0.05 and 0.1. The results reported in Table 2.3 inform that at low concentrations, both the Pareto type IV and uniform distributions perform well. As expected, the fit of the beta distribution improves as c ∗ increases. Figure 2.9: Comparison between the concentration cumulative distribution function (CDF) and the Pareto distribution model for d= 2 (top left) and d= 3 (top right), and between the concentration CDF and the uniform distribution model for d= 2 (bottom left) and d= 3 (bottom right). Results are for the fields with σ 2 f = 1 at time t/τ a = 12.3. 31 Table 2.3: Relative errorsε between the fitted concentration CDFs and the empirical CDF obtained from the numerical simulations. Results computed for different values ofη, see Equation 2.35. η = 0.01 Simulation c ∗ ε for beta ε for Pareto type IV ε for uniform d= 2, point 9.5× 10 − 5 13.97 2.20 7.60 d= 2, line 6.4× 10 − 4 36.55 6.99 13.27 d= 3, point 1.5× 10 − 5 2.02 2.12 1.59 d= 3, line 1.1× 10 − 4 2.58 6.38 3.32 d= 3, plane 7.8× 10 − 4 9.47 3.87 6.25 η = 0.05 Simulation c ∗ ε for beta ε for Pareto type IV ε for uniform d= 2, point 4.7× 10 − 4 5.80 4.21 5.12 d= 2, line 3.1× 10 − 3 13.67 12.01 15.63 d= 3, point 7.2× 10 − 5 2.29 1.19 0.95 d= 3, line 5.2× 10 − 4 1.75 2.07 1.00 d= 3, plane 3.2× 10 − 3 4.59 11.87 4.21 η = 0.1 Simulation c ∗ ε for beta ε for Pareto type IV ε for uniform d= 2, point 9.5× 10 − 4 2.11 3.69 5.27 d= 2, line 6.3× 10 − 3 3.20 11.55 9.51 d= 3, point 1.4× 10 − 4 0.77 0.90 0.21 d= 3, line 1.0× 10 − 3 1.38 0.42 0.81 d= 3, plane 6.2× 10 − 3 1.83 6.58 1.91 32 2.6.3 Dilution index The temporal evolution of the dilution index defined by Equation 2.11 is shown in Figure 2.10. The evolution of the dilution index mirrors the evolution of the mean concentration. Both quantities are related to the mixing volume. The mean concentration is inversely proportional to the mixing volume, while the dilution index is directly a measure for the mixing volume (Kitanidis, 1994). For a homogeneous medium, the dilution index scales as E(t)∼ t (d− d s )/2 . Thus, it increases faster with increasing dimension d of space and decreasing source dimension d s . While the dependence on space and source dimensions is qualitatively similar for heteroge- neous media, dilution is significantly accelerated compared to homogeneous media. In analogy to the behaviors observed for the mean concentration, dilution increases with increasing variance σ 2 f , this means with increasing level of heterogeneity (see top row, Figure 2.10). This is because the macro-scale spreading of the solute plume increases with an increase inσ 2 f . Larger values of σ 2 f implies more erratically shaped plumes and an increase of the plume’s surface area with the surrounding (i.e., ambient) fluid. As a consequence, diffusive and local-scale dispersive solute mass fluxes are augmented. Comparing, for instance, the plots in the bottom row of Figure 2.10, one can observe that, for all source dimensions, the rate of increase of the dilution index is higher when compared to the homogeneous scaling depicted by the dashed lines. The strongest increase, compared to the behavior for homogeneous media, is observed for the planar source case (Figure 2.10, bottom). This is in agreement with what was observed for the mean concentration results reported in Section 2.6.1. Close inspection of Figure 2.10 (bottom row) also reveals that the source zone mainly impacts the rate of dilution. The rate of dilution is highest for a point-like injection. For this case, the particles that constitute the point-like source, in a three-dimensional setting, will all occupy the same position at t = t 0 . However, at the next considered time, t 1 > t 0 , each particle will have moved in a different direction due the combination of advection and dispersion. This means that the mixing volume at t 1 would have significantly increased (due to spreading in all three directions) when compared to the initial volume at t 0 , when it was just a point. 33 Figure 2.10: Dilution index in time for: point source (top left), line source (top right), 3d field with σ 2 f = 1 (bottom left), and 3d field with σ 2 f = 4 (bottom right). 34 2.7 Summary Modeling solute mixing in spatially heterogeneous porous media flows is a challenging task. The multiscale variability of the hydraulic conductivity impacts both spreading and mixing rates of the solute body. The presence of preferential flow paths (i.e. hydrogeological connectivity), low- conductivity zones and abrupt facies transitions require numerical models to have a high spatial resolution. The need for such fine spatial resolution increases the computational costs associated with flow and transport simulations, especially when dealing with three-dimensional stochastic systems. In this Chapter, we performed high resolution numerical simulations to study the temporal evolution of the mixing behavior of a solute plume in heterogeneous porous media. Mixing is quantified in terms of the solute plume’s concentration mean and variance, probability distribution and the dilution index. Our analysis showed how mixing is affected by different factors such as the dimensionality of the flow domain, the degree of heterogeneity in the conductivity field and the geometrical configuration of the source zone. With few exceptions, most studies up-to-date have focused their attention solely to studying mixing in spatially heterogeneous two-dimensional domains for point and line sources. In this Chapter, we provide a systematic investigation of mixing in both two- and three-dimensional spa- tially heterogeneous porous media flows using a computationally efficient transport simulator de- noted PAR 2 developed in Rizzo et al. (2019). PAR 2 code is an open source, GPU accelerated transport simulator based on the random walk particle tracking technique. We tested the accuracy of the numerical results with analytical and semi-analytical results reported in the literature for both homogeneous and heterogeneous porous media. Our results highlight the importance of the source configuration and flow dimensionality on the temporal evolution of different mixing patterns. In particular, hydraulic conductivity fields with d = 3 exhibit faster decays in time of both µ c and σ 2 c respect to d = 2 fields with the same level of heterogeneity. Moreover, such decays of µ c and σ 2 c happen faster for a point source and slower for a planar source, analogously to the analytical results in a homogeneous field. The source 35 configuration also affects dilution, since the dilution index is higher for a planar source and lower for a point source, although such difference seems to decrease when heterogeneity is increased. These results are in agreement with the theoretical analysis of Dentz & de Barros (2013b) where it was analytically shown that both source zone configuration and flow spatial dimension impact the self-averaging transport behavior of the solute plume. Here we expand this theoretical analysis for different mixing metrics and higher levels of heterogeneity in the conductivity field. Finally, we compared the statistics of the concentration values with the beta PDF model. Using a parametric model to approximate the concentration PDF is appealing given that one can estimate the concentration mean and variance from a limited data set (i.e. due to high costs of data acqui- sition). Previous works showed the suitability of the beta PDF model to capture the main features of the concentration field statistics for both point-like and line sources (Bellin & Tonina, 2007; Boso et al., 2013a). Our results show that the performance of the beta PDF improves when the flow field is three-dimensional and transport originates from a vertical planar source. Given that the performance of the beta distribution deteriorates for low concentration values, we propose two alternative distributions to better capture the probabilities of observing c≪ c o . Our results show that the Pareto type IV distribution performs well after employing a power-law transform to the concentration data. Similarly, the uniform distribution also performs well at low concentrations when a logarithm transform is applied to the concentration data. The analysis carried out in this study shows a promising path in employing a mixture distribution approach to model the statistics of solute concentration field. 36 Chapter 3 Relative contributions of permeability heterogeneity and viscosity contrast on solute mixing 3.1 Introduction The fundamental mechanisms controlling transport phenomena in natural porous formations are significantly affected by the spatially random fluctuations of the permeability field at different scales. The presence of heterogeneity in the permeability field leads to the creation of fast flow channels that impact both the dispersive and mixing rate of the solute cloud (Dagan, 1984; Koch & Brady, 1988; Dentz & Carrera, 2007; de Barros et al., 2012a; Cirpka et al., 2012; Dentz & de Barros, 2015; Ye et al., 2016; Kapoor & Kitanidis, 1998) as well as first passage times (Sahimi et al., 1983; Tyukhova et al., 2016; Rizzo & de Barros, 2017). These heterogeneity features have been investigated in applications such as management and remediation of groundwater contamina- tion (de Barros et al., 2012c, 2013), artificial groundwater recharge (Hartmann et al., 2017), human health risk analysis (Maxwell & Kastenberg, 1999; Henri & de Barros, 2016) and CO 2 storage (Hovorka et al., 2004; Dai et al., 2018). An often neglected but important aspect controlling transport is the viscosity contrast between the solute and the ambient fluid. As reported in the lit- erature (Nicolaides et al., 2015), groundwater contaminants such as nonhalogenated semivolatile compounds and jet fuel, are more viscous than water. On the other hand, contaminants such as 37 halogenated volatiles (amongst others) are less viscous than water. Moreover, the physical proper- ties of a solute, such as viscosity, can vary in space and time due to the dependence of the property on the solute concentration and, therefore, mixing (Jha et al., 2011a). The contrast of viscosity be- tween two fluids leads to the well-known phenomenon of viscous fingering, when the less viscous fluid displaces the more viscous fluid (Homsy, 1987). Viscous fingering is a type of hydrodynamic instability known as the Saffman-Taylor instability that occurs in porous media or Hele-Shaw cells under both miscible and immiscible flow conditions (Saffman, 1986). This phenomenon has re- ceived renewed attention due to its role in applications such as enhanced oil recovery (Jha et al., 2011a; Tran & Jha, 2020), geological CO 2 storage (Van der Meer, 1993), and chromatography separation (Broyles et al., 1998). The effect of viscous fingering on spreading and mixing of fluid slugs has been reported through laboratory experiments and numerical simulations. Experiments have been conducted using Hele- Shaw cells (Kopf-Sill & Homsy, 1988; Petitjeans et al., 1999), glass beads or sand packs (Bacri et al., 1992; Held & Illangasekare, 1995a; Kretz et al., 2003), and naturally consolidated rocks (Davies et al., 1991; Kempers & Haas, 1994a). Numerical simulations have been conducted using higher-order finite difference methods (Christie, 1989; Rubin et al., 1993), particle tracking meth- ods (Tchelepi & Orr Jr, 1994), spectral methods (Tan & Homsy, 1988; Zimmerman & Homsy, 1992; De Wit et al., 2005; Mishra et al., 2008) and compact finite difference-finite volume meth- ods (Chen & Meiburg, 1998; Jha et al., 2011a,b, 2013a). The effect of low levels of heterogeneity on spreading of slugs due to viscous fingering has also been studied (Welty & Gelhar, 1991; Tan & Homsy, 1992; Tchelepi & Orr Jr, 1993; Loggia et al., 1996; De Wit & Homsy, 1997a,b; Chen & Meiburg, 1998; Kempers & Haas, 1994a; Sajjadi & Azaiez, 2013; Nijjer et al., 2019). Nicolaides et al. (2015) showed how the interplay between two sources of disorder (i.e. viscosity contrast and physical heterogeneity of the porous medium) impacted macroscopic features of the solute trans- port behavior, namely solute breakthrough and removal times and mixing, in a flow field induced by injection and extraction wells. 38 This Chapter focuses on further understanding of the joint role of medium heterogeneity and viscous fingering on solute transport under miscible flow conditions. The analysis carried out in the present contribution differs from Nicolaides et al. (2015) by considering additional mixing and statistical metrics such as the concentration statistical distribution. Different than the work of Nicolaides et al. (2015), we consider transport under uniform-in-the-mean flow conditions. Our flow configuration mimics an ambient flow, while in Nicolaides et al. (2015), the flow is controlled by injection and extraction wells. The presence of wells in the flow field impact the solute’s mixing behavior and its concentration statistics (see Libera et al., 2017). We are particularly interested in quantifying the relative importance of the two factors in the overall transport behavior. We achieve this goal in two ways. First, we provide a theoretical analysis on how the sources of disorder, namely viscosity contrast and medium heterogeneity, affect the solute transport process. Second, through the use of high-resolution numerical simulations, we examine the effects of both medium heterogeneity and viscosity contrast on the temporal evolution of statistical descriptors of transport (i.e. concentration mean, variance and probability density function) and solute cloud’s contour length. These metrics provide important information about the solute mixing state (Le Borgne et al., 2015a), and therefore its eventual chemical reactivity (Dentz & de Barros, 2015; de Barros et al., 2015), and assist site managers to estimate the risks associated with contamination (Boso et al., 2013a). For example, the effectiveness of groundwater remediation techniques may be re- lated to the contact surface area between two fluids, which is increased by solute spreading (Zhang et al., 2009; Di Dato et al., 2018). We also show how the early breakthrough, associated with the leading edge of the solute cloud, is impacted by permeability heterogeneity and viscous finger- ing. Different from previous works, we introduce two novel metrics that provide a quantitative separation of the impacts of fingering and heterogeneity. These two metrics allow to characterize and clarify the interplay between heterogeneity and viscosity contrast on the macroscopic trans- port behavior. Finally, a relation between the solute cloud’s concentration cumulative distribution function (CDF) and the beta distribution is reported. 39 3.2 Theoretical background 3.2.1 Governing equations For this Chapter, we consider a two-dimensional (2D) porous medium where x = (x 1 ,x 2 ) is the Cartesian coordinate system. The physical domain has dimensions L 1 × L 2 ={(x 1 ,x 2 )|x 1 ∈ L 1 andx 2 ∈ L 2 }. The permeability of the porous medium, denoted by k(x), is assumed to be lo- cally isotropic and spatially heterogeneous. Furthermore, we assume a constant porosity φ. The spatially heterogeneous structure of the log-permeability f(x)≡ ln[k(x)] is considered to be multi- Gaussian and modelled through a Random Space Function (RSF) (Rubin, 2003). The RSF model for f is characterized by its mean⟨ f⟩, varianceσ 2 f , and a correlation lengthλ. The angled brackets ⟨·⟩ represent the average operator. In this study, we adopt an isotropic Gaussian log-permeability covariance model: C f (r) = σ 2 f e − r 2 /λ 2 (3.1) with r=|x ′′ − x ′ | denoting the lag-distance. Our study neglects the effects of density variations. This assumption is reasonable when deal- ing with 2D planar flow fields and for solutes characterized by density values similar to the sur- rounding fluid. Regarding the latter, the densities of certain contaminants are similar to water. Examples consist of m-Cresol, Chlorobenzene and Benzene (with densities equal to 1.03 g/ml, 1.11 g/ml and 0.88 g/ml respectively). For our study, we designate dimensional variables with the hat symbol (ˆ· ). The governing equation for the flow field is thus given by c t ∂ ˆ p(ˆ x, ˆ t) ∂ ˆ t +∇ ˆ x · ˆ q(ˆ x, ˆ t)= 0, (3.2) 40 where ˆ q is the specific discharge, ˆp denotes the pore fluid pressure and c t represents compressibil- ity. Assuming incompressible flow ( c t ≈ 0) in the absence of both sinks/sources and temporally variable boundary conditions, Equation 3.2 becomes ∇ ˆ x · ˆ q(ˆ x, ˆ t)= 0. (3.3) The heterogeneous permeability field is mapped on the specific discharge through Darcy’s law: ˆ q(ˆ x, ˆ t)=− ˆ k(ˆ x) ˆ µ( ˆ C(ˆ x, ˆ t)) ∇ ˆ x ˆ p(ˆ x, ˆ t), (3.4) where ˆ µ is the ambient fluid mixture’s viscosity, which depends on the concentration ˆ C of the injected solute, and ˆ p denotes the pore fluid pressure. In this study we consider that the fluid carrying the solute has a viscosity different than the viscosity of the ambient fluid and we assume an exponential viscosity model (Jha et al., 2011a): ˆ µ( ˆ C(ˆ x, ˆ t))=µ 0 e − R ˆ C(ˆ x,ˆ t) C o (3.5) and R is the log-viscosity ratio, R= ln[µ 0 /µ 1 ], C o is the inlet concentration of the solute, µ 0 is the viscosity of the pure fluid in the absence of a solute, and µ 1 is the viscosity of the solute. We consider permeameter-like boundary conditions, i.e. prescribed pressures at the inlet (x 1 = 0) and outlet (x 1 = L 1 ) of the flow domain and zero-flux at the layer boundaries x 2 = 0 and x 2 = L 2 . Under these conditions, the spatially heterogeneous flow field is driven by uniform-in-the-mean pressure gradient. An inert solute is instantaneously injected through an areal source zone with dimensionsℓ 1 and ℓ 2 (withℓ 1 ≪ ℓ 2 ) into the divergence-free Darcy flow. The spatiotemporal distribution of the solute concentration is provided by the advection-dispersion equation: φ ∂ ˆ C(ˆ x, ˆ t) ∂ ˆ t + ˆ q(ˆ x, ˆ t)· ∇ ˆ x ˆ C(ˆ x, ˆ t)=φD∇ 2 ˆ x ˆ C(ˆ x, ˆ t), (3.6) 41 Figure 3.1: Initial distribution of the concentration field in the flow domain. The solute is instanta- neously injected along a rectangular source zone in a spatially heterogeneous porous medium. where D corresponds to the local-scale dispersion coefficient (assumed to be constant) and is mod- elled through Scheidegger’s theory (Scheidegger, 1961). The boundary conditions for the transport problem are periodic on the left and right boundaries, as well as the top and bottom of the domain. The initial condition for the concentration distribution, as reported in Figure 3.1, is a zero concen- tration in the whole domain with the exception of the source area, where the concentration is set to be equal to unit value, with a small numerical perturbation on the edges to allow the formation of viscous fingers. Appendix A reports the governing equations, i.e. Eqs. (3.3)-(3.6), in dimensionless form, see Eqs. (A.2)-(A.5). As indicated in the appendix, the ˆ· symbol is removed to denote all dimensionless variables (see dimensionless groups, Equation A.1 in the appendix). The appendix also shows that the dimensionless advection-dispersion equation is expressed in terms of the P´ eclet number. The P´ eclet number is given by Pe≡ Uλ/D. Due to the dependence of q on the concentration field C, Equation (A.5) is non-linear and needs to be solved numerically. 42 3.2.2 Mixing descriptors The degree and rate of mixing can be quantified in the terms of the concentration variance ( σ 2 c ). The temporal evolution ofσ 2 c under periodic boundary conditions (or in absence of any net injec- tion or extraction of the fluids) is governed by an ordinary differential equation that relates σ 2 c to the mean scalar dissipation rateε c (Le Borgne et al., 2010; Jha et al., 2011a), dσ 2 c (t) dt =− 2ε c (t). (3.7) The mean scalar dissipation rate is defined in terms of concentration gradients as ε c (t)≡ ⟨|∇C(x,t)| 2 ⟩ Pe , (3.8) where⟨·⟩ denotes the spatial averaging operator. 3.3 Analysis on the sources of disorder in the flow field 3.3.1 Viscosity contrast and heterogeneity effects on vorticity The viscosity contrast R and spatial heterogeneity in the k field induce fluctuations in the flow field, which give rise to fingering and channeling patterns on the solute interface. These fluctuations can be quantified in terms of the vorticity ω ω ω =∇× q of the flow field. Using Equation 3.4 and µ = e − RC , we obtain (Camhi et al., 2000) ω ω ω(x,t)= R∇C(x,t)× q(x,t)+ 1 k(x) ∇k(x)× q(x,t)=ω ω ω R (x,t)+ω ω ω k (x,t), (3.9) which clearly shows how the two sources of disorder –viscosity contrast and heterogeneity– act on the flow field to generate two types of vorticity, ω ω ω R andω ω ω k , defined as ω ω ω R = R∇C(x,t)× q(x,t) andω ω ω k = 1 k(x) ∇k(x)× q(x,t). In a 2D field, only the out-of-plane component of the vorticity vector 43 is non-zero, i.e. ω ω ω =[0,0,ω]. Then Equation 3.9 can be written in terms of the stream functionψ, which is defined as ∂ψ/∂x 1 =− q x 2 and∂ψ/∂x 2 = q x 1 to obtain ω(x,t)=− R∇C(x,t)· ∇ψ(x,t)− ∇ln[k(x)]· ∇ψ(x,t)=[ϑ R (x,t)+ϑ k (x)]· ∇ψ(x,t), (3.10) where we split the vorticity generation mechanisms into its two parts, ϑ R and ϑ k . The stream function and vorticity are also related to each other as∇ 2 ψ =− ω. Let us analyze the two sources individually. • The magnitude of ϑ R (x,t)=− R∇C(x,t) (3.11) increases with non-zero R and∇C, the latter being a dynamic quantity. The interface between the initial solute cloud and the ambient fluid provides the initial ∇C, which is amplified by R to generate the vorticity near the interface that creates viscous fingers and further growth of the concentration gradients. This shows how the fingering process can grow like an insta- bility. The concentration gradient initially sparks the creation of a finger, and the tip of the finger in turn “pushes” to travel faster than the ambient fluid, thus leading to the existence of a concentration gradient on the front of the interface between scalar and fluid; as the finger moves and progresses, such gradient moves along with it in the domain, causing local growth of the concentration gradient. Such local growth can be seen in more detail in Figure 3.2. • The magnitude of ϑ k (x)=− ∇ f(x) (3.12) (with f = lnk) is not, on the other hand, a dynamic quantity and is constant in time for a given permeability field. These differences between the two sources of disorder, namely R andσ 2 f , give rise to different flow structures depending on whether viscosity contrast or heterogeneity dominates the flow (see 44 C(x,t) ln[ε l c (x,t)] Figure 3.2: Concentration field C(x,t) (left column) and logarithm of the local scalar dissipation rateε l c = |∇C(x,t)| 2 Pe (right column) at four subsequent time steps (from top to bottom) for a viscosity contrast R= 1 in a homogeneous permeability field. Note how the concentration gradient pro- gresses in the orange square, moving from the bottom left quadrant to occupy the whole bottom half in the first two time snapshots. As another finger enters the region defined by the orange square, the concentration gradient locally increases also in the upper half of the considered square. 45 R= 3.5,σ 2 f = 0 R= 0,σ 2 f = 3 Figure 3.3: Difference in the finger tip structures between fingering-dominated ( top row) and heterogeneity-dominated (bottom row) flows. The left column shows the concentration fields in the two cases within a zoom-in view window. The right column shows the concentration con- tours (black lines) superposed on the streamlines (red lines). Fingering-dominated flows produce rounder and better delineated tips than heterogeneity-dominated flows. Fingering also produces tip-splitting that is absent at R= 0. Figure 3.3). When fingering dominates heterogeneity (sufficiently large R compared toσ 2 f ), vor- ticityω R is high immediately behind the sharp interface of a finger tip, which grows and becomes “rounder” as time increases until the tip splits into two nascent fingers (Homsy, 1987; Zimmerman & Homsy, 1991; Jha, 2014a). Roundness and tip-splitting of fingers is absent at R= 0 regard- less of how high σ 2 f is. Instead, fingers at high σ 2 f gradually thins towards the tip and the tip fades away in analogy with the Taylor dispersion (Taylor, 1953) effect (Figure 3.3). In viscous fingering-dominated flow, Taylor dispersion at the tip is suppressed due to the vorticity from ϑ R . This suppression happens because the solute cloud that is creating the finger has the tendency to travel fast, pushing against the front and keeping it sharp. In this study, we consider both R> 0 and R< 0 displacements. Fingering occurs at the leading edge of the solute cloud in case of former and at the trailing edge in case of latter. For R< 0, the fingers are developed in the trailing edge of the plume, opposite to the direction of mean flow because their tips lie at the rear front of the less viscous fluid. The sign of R affects the viscosity 46 Figure 3.4: Single finger regions: straight wall, curved rear and curved tip. contrast-induced vorticity,ω R =− R∇C· ∇ψ, in three important ways. The first effect is obvious because R appears explicitly in the expression ofω R . The second effect is through the∇C vector. Note thatω R is nonzero only near the interface because∇C is zero away from the interface. The direction of∇C vector changes along the entire interface and in particular along the interface of a single finger, which can be split into three regions: the curved tip, curved rear, and straight wall regions of the finger, as shown in Figure 3.4. At the tip, ∇C and q vectors are opposite to each other for R> 0, whereas they point in the same direction for R< 0. The third effect is through the magnitude of q. The transverse velocity q x 2 is highest near the tip and lowest along the wall region (Ghesmat & Azaiez, 2008). Moreover, q x 1 and q x 2 magnitudes are higher at the tip of a R> 0 finger than that of a R< 0 finger, as visible from Figure 3.5. 3.3.2 Fingering mechanisms at R̸= 0 To understand the effect of fingering on the transport metrics of the solute cloud, we need to refer to the physical mechanisms of tip-splitting, channeling, shielding, and merging of fingers that have been observed in viscous fingering simulations and experiments (Homsy, 1987; Tan & Homsy, 1988; Zimmerman & Homsy, 1991; Jha et al., 2011b; Tran & Jha, 2020). As mentioned above, tip-splitting refers to the splitting of the tip of a finger into two smaller, nascent fingers. The nascent fingers compete with each other as well as with other fingers in their vicinity by virtue 47 q x 1 (x,t) q x 2 (x,t) Figure 3.5: Longitudinal and transversal velocity fields ( q x 1 (x,t) and q x 2 (x,t) respectively) for a solute cloud with viscosity contrast R= 1 (top row) and R=− 1 (bottom row) in a homogeneous field. The black lines represent the contour of the solute cloud at a threshold concentration value C= 0.1. of the global nature of the pressure field. This manifests as shielding of a smaller finger by a larger, faster-moving finger and merging of the smaller finger into the larger finger. The shielded finger experiences a stunted growth and either merges with the larger finger or fades away through diffusion. At high R, multiple merging and shielding events can lead to emergence of a single finger as a dominant feature of the flow. This is called channeling, which is often undesirable in mixing applications because it leads to a reduction in transverse mixing of the two fluids compared to the scenario where multiple fingers grow on the interface. However, note that channeling also implies a faster breakthrough at pumping wells which may be desirable during some contaminant removal operations. These fingering mechanisms are nonlinear and arise due to two-way coupling between flow and transport processes. They remain active throughout the displacement period, albeit with different strengths. For example, tip-splitting is more pronounced at early times whereas channeling emerges at later times. Also, the effect of two mechanisms can negate each other, e.g. tip-splitting can negate merging in terms of their net effect on the interface length. 48 To understand the impact of the sign of R on the solute transport behavior, we can analyze how tip-splitting and channeling mechanisms manifest themselves at the leading and trailing edges of the cloud. From Equation 3.10, ω R is larger at larger specific discharges q x 1 and q x 2 . Therefore, we expect tip-splitting to be more pronounced at R> 0 than at R< 0 because the velocities q x 1 and q x 2 at the tip are higher at R> 0 than at R< 0. Channeling is almost exclusively reserved for R> 0 flows because the tip velocity (with respect to the mean flow velocity) is smaller in R< 0 than in R> 0. 3.4 Numerical Model and Implementation We use a second-order accurate finite volume method to solve Equation 3.3 for cell-centered pres- sures. We linearize the equation by using the two-point flux approximation (LeVeque, 2002) of the Darcy flux in Equation 3.4, which expresses the face-centered flux in terms of cell-centered pressures of the two neighboring cells and the face transmissibility computed from the harmonic mean of the cell permeabilities, viscosities, and dimensions. To solve the transport problem ( Equation 3.6), we use a spectral method to discretize the spatial derivatives of cell-centered con- centrations and a third-order explicit Runge-Kutta scheme to integrate time (for details, see Jha et al., 2011a,b). Higher-order accuracy of the spectral method allows us to resolve sharp gradients in the concentration field that arise from viscous fingering. The pressure field is relatively smooth because of the assumption of miscibility of the two fluids and, therefore, a second-order finite volume method provides sufficient numerical accuracy for the flow problem. We use a 2D Cartesian grid with a uniform mesh size. Each grid block has dimensions∆× ∆ (see Table 3.1). The number of cells is N x × N y = 1440× 264 such that the ratio∆/λ = 0.084, in order to capture the effects of the spatially heterogeneous permeability field on the spatiotemporal dynamics of the solute cloud (e.g., Ababou et al., 1989a; Bellin et al., 1992b; Rubin et al., 1999). The simulations were performed on high-performance computing nodes. Details about the input parameters for the numerical model are reported in Table 3.1. 49 Table 3.1: Input parameters used in the simulations Parameter Symbol Value Calculated as Correlation length in the x and y-direction λ 5 m - Domain length in the x 1 -direction L 1 600 m 120λ Domain length in the x 2 -direction L 2 110 m 22λ Mesh size in the x 1 - and x 2 -directions ∆ 0.42 m λ/12 Distance of source from top and bottom boundary s tb 15 m 3λ Distance of source from left boundary s l 50 m 10λ Length of source in the x 1 -direction ℓ 1 15 m 3λ Length of source in the x 2 -direction ℓ 2 80 m 16λ Longitudinal distance of control plane from left boundary L cp 250 m 50λ Mean longitudinal flow velocity U 1 m/d - Local dispersivity α 0.042 m 0.1∆ Local scale dispersion coefficient D 0.042 m 2 /d α U P´ eclet number Pe 120 Uλ/D Concentration at the source zone C o 1 mg/l - Threshold concentration for solute cloud delineation C t 10 − 3 mg/l - 3.5 Results All computational results are reported in dimensionless form. Both longitudinal and transverse directions are normalized byλ, time is normalized by the advective time scaleτ u =λ/U and the concentration is normalized by the inlet concentration C o . In Figure 3.6 we report the results of the simulations performed for a k-field characterized by σ 2 f = 0.5 at different times for the three selected values of viscosity contrast R= -1, 0 and 1. Qualitative analysis of Figure 3.6 illustrates that the mixing behaviour is different for different R values, for example the absence of viscosity contrast (R= 0) reduces the degree of mixing of the cloud. The exact quantification of this behaviour is reported in the following sections. It is also clear that in the case of a positive R, the solute cloud travels faster, while the opposite happens for negative values of R. This holds true also for the fields with σ 2 f = 0 and σ 2 f = 0.25 (not shown in Figure 3.6). We can thus expect earlier first arrival times for cases where R> 0 with respect to transport in the absence of viscosity contrast. Figure 3.6 also shows that clouds with R< 0 will be characterized by late arrival times. This becomes clear in Figure 3.7, which illustrates 50 the impact of the mobility ratio R on the concentration breakthrough curve (BTC) at the control plane for different levels of heterogeneity. Note that Figure 3.7 reports the averaged normalized concentration over the control plane (i.e. transverse direction) located at x 1 = L cp , namely: ⟨C(L cp ,x 2 ,t)⟩ ⊥ = 1 L 2 Z L 2 0 C(L cp ,x 2 ,t)dx 2 . (3.13) The case for a homogeneous porous media is illustrated in Figure 3.7.a. As the level of hetero- geneity increases, the maximum concentration decreases and macro-scale spreading is augmented leading to more tailing effects (compare Figures 3.7.a, 3.7.b and 3.7.c). The effects of the mobility ratio are also clearly depicted in Figure 3.7. For R= 1, the presence of fingering is pronounced in the leading edge of the solute cloud which leads to earlier breakthrough at the control plane as opposed to the case of R=− 1. The asymmetry of the concentration signal at the control plane is amplified when R is different than zero. On one hand, the effects of heterogeneity tend to diminish the impact of the mobility ratio on the maximum concentration (compare Figure 3.7.b with 3.7.c, results for R=− 1 and R= 1). On the other hand, Figure 3.7 clearly illustrates that the mobility ratio impacts the spreading of the solute cloud independent of the log-permeability variance. In all cases, the macro-spreading effects (Dagan, 1984) (the plume spreading induced by fluctuations in the sources of disorder) are larger for R< 0. To evaluate the coupled effect of viscosity contrast and heterogeneity on early arrival times, we define the variable η as η = t e | R̸=0 t e | R=0 , (3.14) where t e is the time at which the normalized concentration reaches a value of C = 0.05. The trend ofη for varyingσ 2 f is reported in Figure 3.8, where the blue curve represents the cases with negative viscosity contrast while the light blue the positive one. From Figure 3.8 is possible to notice that fingers in the front of the cloud travel faster with increasing levels of heterogeneity. The opposite occurs when fingers are developed in the tailing edge of the plume (i.e. the cloud slows down with respect to the case with R= 0). This means that the hydraulic connectivity (i.e. 51 Figure 3.6: Two-dimensional representation of the concentration distribution with k-field charac- terized with σ 2 f = 0.5 for three different values of viscosity ratio R at three different simulation times. 52 Figure 3.7: Breakthrough curves (BTCs) at a control plane located at x 1 = L cp for the three different levels of heterogeneity of the permeability field: (a) homogeneous, (b) heterogeneous with σ 2 f = 0.25 and (c) heterogeneous withσ 2 f = 0.5. early arrival times) of the leading edge of the plume is enhanced for R> 0 whereas it is diminished for R≤ 0. These results show that neglecting viscosity contrast in transport models can have an impact on the estimation of early arrival times which are critical in risk analysis and contaminant site management (e.g., de Barros et al., 2012c; de Barros et al., 2016). Next, we evaluate the spatial mean of the concentration (normalized by C o ) at all time steps (Figure 3.9). The spatial mean is calculated in the region of the domain occupied by the solute cloud, defined as the volume V(t) where the concentration is higher than the threshold value re- ported in Table 3.1, C(x,t)> C t . The dimensionless mean concentration⟨C⟩ is thus calculated as: ⟨C⟩= 1 V(t) Z V(t) C(x,t)dV. (3.15) As shown in Figure 3.9.a, for a homogeneous porous media the temporal evolution of the mean concentration is the same for both R= 1 and R=− 1. Furthermore, the mean concentrations for R̸= 0 cases are lower than the one obtained in the absence of viscous fingering ( R= 0). When R̸= 0, the presence of fingering contributes to an increase of the surface area of the solute cloud with the surrounding ambient fluid, also called interface stretching. It also contributes to a higher concentration gradient along the interface. These two effects combined contribute to an enhance- ment of the diffusive mass flux across the interface, which in turn leads to a more diluted cloud. Furthermore, more vigorous tip-splitting in case of R= 1 increases the length of the fluid-fluid 53 Figure 3.8: Effect of the viscosity ratio R onη for increasing levels of heterogeneity of the perme- ability field. Results computed at x 1 = L cp , where the value of L cp is reported in Table 3.1. interface compared to R=− 1, which explains the larger dilution or smaller mean concentration observed at R= 1. This explains also why, if we consider a scaling relation for the decay of the mean concentration⟨C⟩∼ t b at large times, the values of b are always slightly higher in modulus for R< 0 than for R> 0. In fact, in the latter case the enhanced mixing at earlier times results into smaller magnitudes of the gradients at intermediate times which retards mixing at later times. Both cases however present a decay faster than the case of R= 0. The effects of the mobility ratio are reduced when the level of heterogeneity increases (compare Figures 3.9.a, 3.9.b and 3.9.c). Figure 3.10 depicts the temporal dynamics of the concentration variance (normalized by C 2 o ) within the solute cloud of volume V(t), calculated as: σ 2 c (t)= 1 V(t) Z V(t) (C(x,t)−⟨ C(t)⟩) 2 dV. (3.16) As observed in Figure 3.10, the presence of fingering causes a sudden increase of the concentra- tion variance followed by fast decay when compared to the R= 0 case, regardless of the level of 54 Figure 3.9: Mean plume concentration for the three different levels of heterogeneity of the per- meability field: (a) homogeneous, (b) heterogeneous with σ 2 f = 0.25 and (c) heterogeneous with σ 2 f = 0.5. In the insets the detail of late times plotted in a log-log scale, with the values of the parameter b. 55 Figure 3.10: Concentration variance of the solute cloud for the three different levels of hetero- geneity of the permeability field: (a) homogeneous, (b) heterogeneous with σ 2 f = 0.25 and (c) heterogeneous withσ 2 f = 0.5. 56 heterogeneity. Similar to the mean concentration behavior observed in Figure 3.9, more vigorous tip-splitting in case of R> 0 leads to an increase of the surface area of the solute body with the surrounding fluid which contributes to the dilution of the cloud and therefore, lower variance (com- pare curves for R=− 1 and R= 1 in Figures 3.10.b-c). It is interesting to note the difference in the impact of R between the heterogeneous and homogeneous cases. In fact, for the heterogeneous cases, the curves for R=− 1 and R= 1 become more distinct with respect to the homogeneous one. Figure 3.11: Cumulative density function of the plume concentration in the k-field with σ 2 f = 0.5 at three different simulation times for three different values of viscosity ratio R: (a) R=− 1, (b) R= 0 and (c) R= 1. The cumulative distribution function (CDF) of the concentration field is plotted at different times forσ 2 f = 0.5 and R= -1, 0 and 1 in Figure 3.11. As shown in Figure 3.11, the probability that the concentration is below a certain value (i.e. C= 0.25) is higher for R= 1 (Figure 3.11.c) when compared to the case of R=− 1 (Figure 3.11.a). The tip-splitting mechanism, which is stronger for positive values of R, is the key reason for the increase in the likelihood of observing lower concentrations. As shown in the literature (de Barros & Fiori, 2014), the concavity of the CDF can be used as a measure of the dilution state of the solute body. Based on the evidence published in previous papers in the turbulence, atmospheric and hydrologic communities (Wall et al., 2000; Munro et al., 2003; Bellin & Tonina, 2007; Boso et al., 2013a), the concentration CDF can be approximated by a beta distribution for transport in spatially variable flow fields. Here 57 we examine if the beta distribution can approximate the concentration distribution in the presence of viscous fingering. The beta CDF is given by: P C (c;x,t) = Γ[a o + b o ] Γ[a o (x,t)]Γ[b o (x,t)] Z c 0 χ a o (x,t)− 1 (1− χ) b o (x,t)− 1 dχ, with Γ[z] = Z ∞ 0 ζ z− 1 e − ζ dζ ; a o (x,t)= ⟨C(x,t)⟩ β(x,t) ; b o (x,t)= 1−⟨ C(x,t)⟩ β(x,t) β(x,t) = σ 2 c (x,t) ⟨C(x,t)⟩(1−⟨ C(x,t)⟩)− σ 2 c (x,t) , (3.17) whereΓ[z] is the Gamma function. Figure 3.11 depicts an excellent agreement between the empirical CDF (obtained from the raw numerical data) and the beta CDF model ( Equation 3.17) parameterized by the mean and variance of C, namely⟨C⟩ and σ 2 c . From an application point of view, both quantities, i.e. concentration mean and variance, can be estimated from monitoring wells in a contaminated site (Boso et al., 2013a). The results reported in Figure 3.11 show that the beta-CDF model can be employed to capture the full probabilistic distribution of the solute cloud in the presence of viscous fingering for the range of heterogeneity and mobility ratio explored in this study, i.e. σ 2 f < 1 and − 1≤ R≤ 1. This result can be used to predict the probability that the solute concentration will exceed a threshold value which is of importance in risk analysis and aquifer remediation (Boso et al., 2013a; de Barros & Fiori, 2014). The temporal evolution of the mean scalar dissipation rate ( Equation 3.8) is presented in Figure 3.12 for the homogeneous and heterogeneous flows. As illustrated in Figure 3.12, the presence of viscous fingering leads to a sharp increase in the scalar dissipation rate followed by a decay. The differences between the curves for R= -1, 0 and 1 are attenuated in the heterogeneous scenarios (compare Figure 3.12.a with Figures 3.12.b-c). In all cases, the peak of the scalar dissipation rate is higher for R= 1. Figure 3.13 illustrates the solute cloud’s contour length evolution in time. The contour length is computed based on the normalized concentration C ∗ = 0.05 isoline. As shown in Figure 3.13, 58 Figure 3.12: Mean scalar dissipation rate ε c for the three different levels of heterogeneity of the permeability field: (a) homogeneous, (b) heterogeneous with σ 2 f = 0.25 and (c) heterogeneous withσ 2 f = 0.5. Figure 3.13: Evolution of the length of C= 0.05 contour for the three different levels of hetero- geneity of the permeability field: (a) homogeneous, (b) heterogeneous with σ 2 f = 0.25 and (c) heterogeneous withσ 2 f = 0.5. 59 the effect of the medium’s disorder contributes to an increase of the cloud’s contour length. For the heterogeneous cases, i.e. σ 2 f = 0.25 and 0.5, the contour length is larger for R= 1 than the ones computed for R=− 1. Figure 3.14: Trends in time of the mean scalar dissipation rate ε c normalized with respect to (a) the effect of viscosity contrast, defined by the metric ϕ ( Equation 3.19) and to (b) the effect of heterogeneity, defined by the metric ξ ( Equation 3.18). Next, we quantify the relative contribution of the medium’s disorder level and the viscosity contrast on solute mixing. In order to do so, the effects of heterogeneity and viscosity contrast have been here quantified by normalizing the mean scalar dissipation rate ε c (Figure 3.12) to the cases where those effects are not present. The heterogeneity-induced dissipation is represented by the parameterξ thus calculated as: ξ(t)= ε c (t|R,σ f )− ε c (t|R,σ f = 0) ε c (t|R,σ f = 0) , (3.18) while the finger-induced dissipation is quantified by ϕ as: ϕ(t)= ε c (t|R,σ f )− ε c (t|R= 0,σ f ) ε c (t|R= 0,σ f ) . (3.19) The temporal evolution ofξ andϕ is reported in Figure 3.14 for the cases where both heterogeneity and viscosity contrast are present. It is clear that at early times, the effects induced by the viscosity 60 contrast are more pronounced on the dissipation rate as opposed to the ones induced by the physical heterogeneity, as the values ofϕ in Figure 3.14.a are higher than the values ofξ in Figure 3.14.b. From Figure 3.14.a, we notice that for the same level of heterogeneity, a positive viscosity contrast has a greater impact on ε c compared to a negative R case. At parity of R, the fingering effect on mixing is attenuated by higher levels of heterogeneity, which is as expected. At later times, on the other hand, heterogeneity (Figure 3.14.b) becomes the main contributor to dissipation, though its importance is relatively marginal compared to the viscosity contribution to mixing at early times. The behaviour ofξ at parity of R does not seem to be greatly affected by an increase inσ 2 f , but this may be due to the fact that the considered levels of heterogeneity do not differ much. However, heterogeneity seems to enhance mixing at later times especially for R< 0; this can be explained considering that mixing happens faster for R> 0, so at later times the process is almost over, while for R< 0, heterogeneity has sufficient time to kick in and have a more significant impact. This suggests that, for the simulations considered here, the time scale of fingering-induced mixing mechanism is shorter than that of the heterogeneity-induced mixing mechanism. The two time scales can further be related to the length scales of diffusion (Pe) and permeability correlation length (λ). The prescribed initial condition on the concentration field, i.e. the sharpness of the interface at t= 0, also affects the separation of two time scales in this plot (Tan & Homsy, 1992). 3.6 Summary In this study we investigate the impact of two sources of disorder, namely viscous fingering and permeability heterogeneity, on the temporal evolution of key statistical descriptors of transport. Our numerical simulations indicate that: 1. Heterogeneity reduces the peak concentration observed at the control plane and, while it does not affect the impact of the viscosity ratio, it produces asymmetry in the BTCs and a larger cloud spreading for R< 0. 61 2. A positive viscosity contrast enhances the connectivity of leading edge of the solute cloud which leads to earlier breakthrough, while a negative viscosity contrast delays the arrival of the solute to the control plane; this opposite impact of viscosity on connectivity increases with increasing levels of heterogeneity. 3. At late times, the mean concentration decay is described by a relation of the type ⟨C⟩∼ t b , where the effect of the mobility ratio on the value of b is reduced for higher level of heterogeneity. 4. The concentration CDF, whose concavity can be used as a proxy for dilution (de Barros & Fiori, 2014), can be perfectly described as a beta distribution (parameterized solely by the concentration spatial mean and variance) in the presence of viscous fingering for the range ofσ 2 f and R considered in the present study. This result can be used for predictions in risk analysis, and to the authors’ knowledge has never been shown before for flow in the presence of viscosity contrast. 5. The mean scalar dissipation rate for R̸= 0 presents a sharp increase at early times, with higher peaks for R> 0 because of the more intense tip-splitting. 6. Longer cloud contours (and thus mixing interfaces) are caused by higher heterogeneity levels and viscosity contrast, with R> 0 having a bigger impact on the contour than R< 0. 7. At early times the effect of viscosity contrast on mixing is more significant than the effect of heterogeneity, despite being attenuated for increasing values of σ 2 f . At later times, the impact of heterogeneity becomes the main contributor to mixing, especially for R< 0. The results shown in this Chapter illustrate the importance of viscous fingering and hetero- geneity in controlling the temporal evolution of the statistical descriptors of the solute cloud. The conclusions reported are limited to the range of values explored forσ 2 f and R. 62 Chapter 4 Influence of initial plume shape on miscible porous media flows under density and viscosity contrasts 4.1 Introduction Improved understanding of the mixing mechanisms of fluids in a porous medium is of relevance to a wide range of applications in hydrology, water resources engineering, enhanced oil recovery, geothermal systems and risk analysis of groundwater pollution. Several factors impact the spa- tiotemporal transport dynamics of a solute injected in a saturated porous medium. These consist of properties characterizing both the porous medium and the injected fluid. Another important factor is the initial configuration of the injected solute. In this work, we investigate how the initial configuration of the injected fluid and the mean flow velocity control the relative impact of multi- ple sources of disorder on the mixing evolution of two miscible fluids (i.e., ambient fluid and an injected solute). The main sources of disorder examined in this study are the viscosity and density contrasts between the injected solute and the ambient fluid. We also analyze the effect of a third source of disorder, permeability heterogeneity of the porous medium, in the form of stratification. Much of the previous work focused on the effects of permeability variability on mixing (Le Borgne et al., 2015a; Dentz et al., 2018) and dispersion (Dagan, 1984; de Barros & Rubin, 2011). This is due to the fact that the variability of the permeability field over multiple length scales leads to enhanced solute mixing (Dentz et al., 2011; Bonazzi et al., 2022). However, other mechanisms can 63 impact solute mixing and spreading rates. For example, differences in viscosity between the added solute and the ambient fluid impacts the mixing dynamics in porous media. When an injected fluid is less viscous than the ambient fluid, fluid displacement can lead to the formation of a hydrody- namic interface instability denoted as the Saffman-Taylor instability. This hydrodynamic insta- bility generates viscous fingering (Saffman, 1986; Homsy, 1987) which can potentially augment dilution rates (Jha et al., 2011b, 2013b). Due to the spatiotemporal dependency of the viscosity on the solute concentration, plume mobility and, consequently, mixing, are significantly affected by viscosity contrast (Jha et al., 2011a; Nijjer et al., 2018; Tran & Jha, 2020, 2021; Bonazzi et al., 2021). The interaction between viscosity contrast and permeability heterogeneity and its effect on mixing have been studied for low levels of heterogeneity (Tan & Homsy, 1992; Tchelepi & Orr, 1994; De Wit & Homsy, 1997a,b; Nicolaides et al., 2015; Bonazzi et al., 2021). These studies shed light on how and when solute mixing are affected by the fluctuations of the permeability field. Through the use of high resolution numerical simulations, Bonazzi et al. (2021) showed how key transport metrics, such as the solute arrival times and the temporal evolution of spatial statistical descriptors of the plume, were impacted by a spatially heterogeneous multi-Gaussian random log-permeability field. Other studies investigated the effect of permeability heterogeneity and viscosity contrast on mixing by considering layered porous media (Sajjadi & Azaiez, 2013; Nijjer et al., 2019). Density contrast between two fluids is another important source of disorder. A type of interface instability called Rayleigh-Taylor instability can be triggered when a more dense fluid is located above a less dense one, causing the formation of gravity-driven fingers (Gopalakrishnan, 2020; Slim, 2014; Jenny et al., 2014). Under these conditions, analogous to the case of viscous fingering, density is a function of the solute concentration, thus creating an interdependence between plume mobility, mixing and density distribution (Hidalgo & Dentz, 2018b; Hassanzadeh et al., 2007). Studies that focus on the effect of density-driven fingering are usually presented in the context of CO 2 sequestration in saline aquifers (Van der Meer, 1993; Hidalgo & Carrera, 2009; Riaz et al., 64 2006). Variable-density flow and its impact on solute transport in coastal aquifers is also reported in the literature (Dentz et al., 2006). The combined effects of viscosity and density contrast has received progressively more at- tention given its importance in many environmental and energy-related processes such as CO 2 sequestration and aquifer remediation. Tchelepi & Orr (1994) proposed a hybrid particle tracking- finite difference method to analyze the effects of viscosity and density contrasts in 2D and 3D heterogeneous porous media. Hidalgo et al. (2013) used numerical simulations to examine the dis- solution flux of injected carbon dioxide in a homogeneous medium. Daniel & Riaz (2014) focused on the effect of viscosity contrast on a gravitational unstable interface, and Pramanik & Mishra (2016) studied the instability in a vertical flow of two fluids with different densities and viscosi- ties under distinct mean flow velocities. Both analytical and semi-analytical solutions have been proposed to predict the location of the interface between the formation brine and the CO 2 rich phase while accounting for contrast in viscosity and density (Nordbotten et al., 2005; Vilarrasa et al., 2010). Experimental laboratory studies, aimed at investigating the joint effects of density and viscosity contrasts on fluid mixing, have also been conducted, amongst others, by Held & Illangasekare (1995b) and Liyanage et al. (2020). Some other studies focused on fluid mixing when both viscosity and density contrasts between a solute and the ambient fluid are present in a heterogeneous permeability medium (Kempers & Haas, 1994b). Loggia et al. (1996) focused on transport in a multi-layered medium. The authors showed that for certain levels of density and viscosity contrast, the flow field is sufficiently stable to mask the effect of heterogeneity, causing a multi-layered concentration profile to merge into a single one (see Figure 2 of Loggia et al., 1996). Despite significant efforts, there is still a need to better understand how these sources of dis- order impact the dynamics of a solute plume under different geometrical configurations of the initial solute source. Here, we define this initial condition as the zone where a fluid, characterized with viscosity and density different than the ones of the ambient fluid, will be injected. The ini- tial configuration of the solute plume is particularly relevant in applications such as groundwater management and risk assessment. In these applications, contaminant plumes can originate from 65 leaking tanks and pipes, accidental spills, injection wells or landfills. In the context of aquifer clean-up, remediation fluids (e.g. remedial reagents) are injected through wells. Injection wells are also utilized to inject CO 2 into the subsurface. In all of these cases, solutes enter the geologi- cal formation by means of different spatial configurations. Furthermore, as previously mentioned, most of the injected fluids have viscosity and density that differ from the ambient fluid. Several works showed the importance of the solute source configuration on mixing and spread- ing dynamics in the absence of viscosity and density contrasts (Dentz et al., 2000b; Dentz & de Barros, 2015; Bonazzi et al., 2022). Koch & Nowak (2015) showed how the uncertainty of con- taminant source architecture impacted fate and transport of nonaqueous-phase liquids in aquifers which are characterised by density and viscosity values that differ from groundwater. To the best of our knowledge, a systematic analysis of the impact of the geometrical configuration of the in- jection zone on mixing in the presence of viscosity and density contrasts has not been investigated, despite its importance in applications in subsurface energy, risk analysis, and groundwater remedi- ation. The investigations carried out in this Chapter show that the initial configuration of the solute plume controls the relative importance of these sources of disorder, namely viscosity and density contrasts. Intuitively, one could expect the dominance of viscous fingering when the solute source is positioned perpendicular to the direction of flow. Similarly, more intense density fingering (also known as “fingered tongue”) could be predicted when the solute source is parallel to the direction of the flow. In order to study these effects for different ranges of viscosity and density contrasts between the solute and ambient fluid, we analyze transport for an initial solute source shaped like an ellipse of radii a and b. We examine how the transport behaviour is affected by different values of the ratio a/b, as illustrated in Figure 4.1. We also study how varying the mean flow velocity impacts our results. After analyzing the results for a homogeneous porous medium, we add permeability heterogeneity, in the form of a stratified permeability field, as an ulterior source of disorder. This paper is structured as follows. Section 4.2 provides the equations governing fluid flow and transport in a porous medium. Details regarding the numerical schemes and implementation are 66 Figure 4.1: Schematic configuration of the problem studied. A solute is injected within an elliptical source zone into a confined porous medium. Upper Figure: Transport of an elliptical solute plume of radii a and b across an aquifer layer of thickness H confined between impermeable layers. The Figure shows the vertical cross-section across the layers. The impact of the source zone on transport is examined by varying the ratio a/b. Lower Figure: Different initial shapes of the plume investigated in this study. The solute source zones range from vertically oriented (a/b= 1/3) to horizontally oriented (a/b= 3). reported in Section 4.3. We compare the output of our numerical simulations with existing results from the literature in Section 4.4. The computational analysis is carried out for homogeneous and heterogeneous media in Sections 4.5 and 4.6. The main findings of this Chapter are summarised in Section 4.7. 67 4.2 Mathematical description of flow and transport 4.2.1 Governing equations In this Chapter we consider a 2D porous medium in a dimensionless Cartesian coordinate system x=(x,z). The dimensionless lengths of our domain are L× H. All equations and quantities in this section are normalized and the dimensionless groups employed are reported in Appendix B. We start by considering a homogeneous porous medium, characterised by a constant permeability ˆ k. Assuming incompressible flow, in the absence of temporally variable boundary conditions and sources and sink terms, the governing (dimensionless) equation for the flow field is given by: ∇· q q q(x,t)= 0, (4.1) where q q q denotes the specific discharge, defined by Darcy’s law as: q q q(x,t)=− 1 µ(c(x,t)) [∇p(x,t)− F(c(x,t))e e e z ]. (4.2) The velocity field is given by q/φ, whereφ is the medium porosity. In Equation 4.2, p(x,t) denotes the pressure field, and µ(c) the viscosity of the ambient fluid mixture of concentration c. The unit vector along the vertical z direction is denoted by e e e z . For our work, we adopt an exponential viscosity model (Jha et al., 2011a) to relate µ and c. An exponential viscosity model can capture the concentration-dependence of viscosity of mixtures such as water and glycerol (Petitjeans & Maxworthy, 1996). The functional relationship betweenµ and c used in this Chapter is: µ(x,t)=µ 1 exp[R(1− c(x,t))], (4.3) where µ 1 is the viscosity of the injected solute, and R is the log-viscosity ratio R= ln(µ 0 /µ 1 ), in which µ 0 is the viscosity of the ambient fluid. For reference, in geological carbon sequestration 68 applications, the log-viscosity ratio R ranges between 1.9 and 2.7 for the viscosity values of CO 2 and water brine reported in Frailey & Leetaru (2009). The functionF(c) in Equation 4.2 originates from the density profile model adopted in Daniel & Riaz (2014): ρ(c(x,t))=ρ 0 +∆ρF(c(x,t)), (4.4) whereρ 0 is the ambient fluid density (i.e., the density at c= 0) and∆ρ =ρ m − ρ 0 , in whichρ m is the maximum possible density of the mixture. In our work, we will set F(c)= c, thus adopting a linear density profile (see Daniel & Riaz, 2014; Horton & Rogers Jr, 1945; Cowell et al., 2020). A linear density profile is able to describe the density-dependence on concentration of a potassium permanganate and water mixture; such mixture has been employed in experimental studies such as Slim et al. (2013) as a model for CO 2 in brine. For the linear model employed, the maximum density ρ m will coincide with the density of the injected solute ρ 1 (i.e., the density at c = 1). Therefore Equation 4.2 becomes: q q q(x,t)=− 1 µ(c(x,t)) [∇p(x,t)− c(x,t)e e e z ]. (4.5) The governing dimensionless transport equation for an inert substance is given by: ∂c(x,t) ∂t + q q q(x,t)· ∇c(x,t)− 1 Ra ∇ 2 c(x,t)= 0, (4.6) in which the Rayleigh number Ra is defined as Ra= k c ∆ρgℓ φDµ 1 , (4.7) where k c is the medium characteristic permeability (for the homogeneous medium, k c = ˆ k),ℓ is the domain characteristic length (here,ℓ represents the vertical dimension of the porous formation),φ is the medium porosity and D is the isotropic local-scale dispersion coefficient. Both φ and D are assumed to be constant. 69 4.2.2 Boundary and initial conditions The boundary conditions for the flow field are, as shown in Figure 4.1, a no flux condition on the top and bottom of the domain, a constant inlet flux f L on the left boundary and an assigned constant pressure p R on the right boundary. These are mathematically equivalent to: ∂ p ∂z z=0 = ∂ p ∂z z=H = 0, (4.8) ∂ p ∂x x=0 = f L , (4.9) p x=L = p R . (4.10) The dimensionless inlet flux f L can be understood as a measure of the relative strength of horizontal advection compared to diffusion (equivalent to a P´ eclet number), whereas Ra controls the relative strength of vertical convection compared to diffusion. Therefore, we expect both parameters to control the strength of convective instabilities - viscous fingering and gravitational fingering - with the effect of Ra being limited to the latter. For the transport problem, we adopt a no-flux boundary condition at the top and bottom edges of the computational domain, analogously to the flow problem, which is consistent with the as- sumption of confining impermeable layers, ∂c ∂z z=0 = ∂c ∂z z=H = 0. (4.11) The no flux-boundary condition applies to the left boundary as well: ∂c ∂x x=0 = 0, (4.12) while on the right boundary of the domain we assume a natural outflow boundary condition, i.e. the solute mixture is allowed to exit the domain through convection only, following the direction of the average flow velocity. 70 The initial condition for the concentration field of the injected fluid is a zero concentration in the whole domain except where the initial source is located. Our initial source zone is represented by an ellipse of radii a and b along the longitudinal and transverse directions respectively, as displayed in Figure 4.1. Within the source zone, represented here by the areaS of the ellipse, the value of concentration is unitary at initial time t= t 0 , c(x,t 0 )=        1, for x∈S, 0, for x/ ∈S. (4.13) By changing the ratio a/b, we can explore the effects of source orientation (figure 4.1) on the spatiotemporal dynamics of c and mixing metrics. Note that the area of the ellipseS is maintained constant throughout all the values of the ratio a/b explored. The theory of linear stability analysis of viscosity and density-driven fingering can provide the instability onset time and the wave number vs. growth rate relation for the early-time behavior of the system. The linear analysis, including the non-modal analysis, is valid before the fingers start to interact via non-linear mechanisms such as tip-splitting, shielding, coalescence, coarsening, and channeling. Given our focus on plume spreading and mixing, we use direct numerical simulations, instead of the linear theory, to resolve such nonlinear interactions and their impact on spreading and mixing metrics. 4.3 Numerical implementation The numerical implementation of our model is similar to the one presented in Bonazzi et al. (2021), with the due modifications in the discretised equations to account for the concentration dependence on the solute density. To solve Equation 4.1 numerically for cell-centered pressures, we employ a finite volume method at second order accuracy (Ferziger et al., 2002b). Since Equation 4.1 is expressed in terms of the Darcy flux, i.e. the face-centered flux at the cells interface, we lin- earize Equation 4.1 using the two-point flux approximation (LeVeque, 2002) of the Darcy flux in 71 Parameter Symbol Values Inlet flux at the left boundary f L 1× 10 − 2 , 1× 10 − 3 Aspect ratio of the initial source a/b 1/3,1/2,1,2,3 Log-viscosity ratio R [-3.5 : 0.5 : 3.5] Rayleigh number Ra [1000 : 500 : 10000] Table 4.1: Model input parameters. All parameter values are dimensionless following the di- mensionless groups reported in Appendix B. R and Ra are varied over a range of values that are uniformly spaced between a minimum value (value left of the first colon) and a maximum value (value right of the second colon) with a constant increment or spacing between successive values (middle value between the two colons). Here R= ln(µ 0 /µ 1 ) and Ra is defined in Equation 4.7. Equation 4.2. In this way, we express the flux in terms of the pressure gradient, the concentra- tion distribution, and the cell interface transmissibility, which depends on viscosity. The transport problem, governed by Equation 4.6, is solved explicitly for the concentration field, c(x,t), in terms of the advective and diffusive concentration fluxes computed at the previous time step using a second-order finite difference method. In our numerical model we use a 2D Cartesian grid where each cell has dimensions ∆x× ∆z, with ∆x= 8× 10 − 3 and ∆z= 4× 10 − 3 . The number of cells in both directions, namely N x and N z , is the same, N x = N z = 250. Thus, the domain lengths in the longitudinal and vertical directions are respectively L× H = 2× 1. Our choice for ∆x and ∆z are based on a numerical grid resolution analysis that focused on quantities of interest such as the horizontally averaged concentration profile in the vertical direction and the plume’s concentration variance, to ensure that the relative percentage error with respect to a more refined grid ( N x = N z = 1000) was below 5%. The initial source zone in our model is an ellipse with radii a and b in the horizontal and vertical direction respectively. In our paper, we investigate different ratios a/b (see Figure 4.1 and table 4.1) in order to explore the effects of the initial plume configuration on transport. Although the ratio a/b is varied, the values of a and b are calculated such that the area of the ellipse,S , at t 0 is the same in all the simulations; this way the same mass of solute is initially injected in the domain. For all the ratios explored, we maintain the center of the ellipse at the dimensionless coordinate x c =(0.3,0.5). 72 4.4 Comparison with previous work In this section we compare the results originating from the numerical implementation of the gov- erning equations as described in section 4.3. We test the performance of the developed numerical simulator using the same physical setup as reported in the work of Hewitt et al. (2013). Ap- pendix C provides a brief description of the physical setup and parameter values used by Hewitt et al. (2013). A relevant difference between our numerical simulation and the one in Hewitt et al. (2013) is that we chose a horizontal discretization with ∆x = (1024) − 1 , while in Hewitt et al. (2013) ∆x=(2048) − 1 . This was done in order to reduce the computational burden of the sim- ulation. Another difference is that the discretization in the vertical direction is ∆z=(250) − 1 in Hewitt et al. (2013), and a vertical coordinate transformation was used to ensure high resolution of the dynamics near the fluid-fluid interface. Here, we opt for a constant vertical discretization, ∆z=(750) − 1 . Figure 4.2 (left) shows a comparison between the horizontally averaged concentra- tion profiles obtained from our simulations and the ones reported by Hewitt et al. (2013) at distinct times. We can see that there is a reasonably good agreement between the two results for the con- sidered times, although some discrepancies are present. Such discrepancies are expected given that the hydrodynamic instability, which causes the fingering behavior, is triggered in numerical models by the introduction of a small random numerical perturbation which is amplified in time by the physics of hydrodynamic instability and the properties of the numerical method used. Dif- ferences in the amplitude and structure of the initial perturbation, which are often unreported in the studies, and difference in the accuracy of the numerical discretization method can also cause differences between studies. Figure 4.2 (right) also illustrates a comparison of the temporal evo- lution of the horizontally averaged inlet solute mass flux at the top boundary of the domain. As shown in figure 4.2 (right), our computational results are based on a single realization of the ini- tial perturbation, therefore displaying some noise. The mass flux estimates obtained from Hewitt et al. (2013) is smoother given that it represents an average from multiple realizations of the initial perturbation. 73 Figure 4.2: Left: Comparison between the horizontally averaged concentration profiles, ¯c(z,t), of Figure 4 in Hewitt et al. (2013) (circles) and the ones obtained in our simulation (solid line), at three selected times. Right: Comparison between the mass flux ˙m s evolution of Figure 5b in Hewitt et al. (2013) (dashed line) and the one obtained in our simulation (solid line). All quantities are dimensionless in accord with the normalization provided by Hewitt et al. (2013) and reported in Appendix C. 4.5 Results and discussion In this section we perform a series of numerical experiments with the model input parameter values reported in table 4.1. In Section 4.5.1, we provide a general analysis of the overall behaviour of miscible flow systems for different values of R, Ra, a/b, and two different values for the inlet flow rate ( f L ). We illustrate how these different parameters impact the solute plume configuration and the corresponding flow field. Section 4.5.2 analyzes how the aforementioned parameters impact specific metrics that quantify mixing and spreading of the plume. Although we have performed simulations for several values of the source zone configuration, we limit our illustrations for the cases of a/b= 1/3, 1 and 3 for the sake of brevity (simulations were also carried out for a/b=1/2 and 2 but are not shown). 74 a/b= 1/3 a/b= 3 R= 3 Ra= 1000 R= 0 Ra= 1000 R= 0 Ra= 9000 R=− 3 Ra= 9000 Figure 4.3: Concentration fields for the slower inlet flux case, f L = 0.001, at the simulation time t= 1. The left column shows results for the vertical source zone (a/b= 1/3), and the right column shows results for the horizontal (a/b= 3) source zone. All cases are for a plume heavier than the ambient fluid ( Ra> 0). The top two rows present cases with low density contrast (Ra= 1000), while the bottom two rows present high density contrast cases (Ra= 9000). R= 3 implies a less viscous solute compared to the ambient fluid, R=− 3 implies a more viscous solute, and R= 0 implies no viscosity contrast. 75 4.5.1 Analysis of flow and transport behaviour Figure 4.3 shows the spatial distributions of the concentration field at dimensionless time t = 1 from a set of eight representative simulations characterized by distinct values of R, Ra, and a/b. The results displayed in Figure 4.3 corresponds to f L = 0.001, which we shall denote here as the low inlet flux case. The plumes depicted in Figure 4.3 originate from two different initial source configurations, namely a vertical source ( a/b= 1/3) and a horizontal source (a/b= 3) (see Fig- ure 4.1 for reference). Comparing the plume snapshots in the left and right columns of Figure 4.3, we observe that the initial plume configuration has a significant impact on the spatial distribution of the solute plume. One of the mechanisms driving this dependency on the initial configuration is gravitational fingering, the strength of which depends on the orientation of the source zone, e.g., compare the plumes obtained for Ra= 9000 and R= 0 for both a/b= 1/3 and 3. Since the gravitational instability grows downward from the plume’s interface, a horizontal initial source configuration ( a/b= 3) provides a longer interface for the gravitational fingers to form and grow. The number of these fingers at early time depends on the typical finger width and the horizontal extent of the source. The typical finger width is a function of Ra, as per the linear stability analysis of gravitational fingering (Riaz et al., 2006), which predicts the finite size wavelength of initial fingers. This is best shown by comparing the second and third rows in Figure 4.3 for a/b= 3; an increase in Ra from 1000 to 9000 leads to a decrease in the finger width which means the same interface length can host three fingers instead of one. The case of a vertical initial plume does not provide enough interface length for gravitational fingers to form (i.e., a/b= 1/3, left column, Figure 4.3). Higher Ra means sharper plume interface and slower effective diffusion, as visible from comparing second and third rows. The results depicted in the first row of plots in Figure 4.3 also highlight the effects of R> 0 in rapidly diluting the plume via viscous fingering, which (unlike gravitational fingering) does not have a preferential direction. Viscous fingering takes place in both vertical and horizontal directions, and in the case of vertical direction, the effect of viscous fingering superposes on the effect of gravitational fingering, leading to complex plume shapes such as the one shown in the 76 first row, right column of Figure 4.3. A horizontal initial source under relatively low Ra leads to one gravitational finger growing downward and one viscous finger growing rightward from the plume interface (figure 4.3, first row, right column). Due to viscous instability at the tip of the gravitational finger, the finger travels faster and reaches the bottom boundary earlier than R≤ 0 cases (second, third and fourth rows). On the other hand, due to gravitational instability, the horizontally traveling viscous finger bends downward. Because of the absence of surface tension effects in a miscible flow, the gravitational and viscous fingers remain connected to the main plume body, which stretches and distorts in a manner reminiscent of chaotic mixing (Ottino, 1989). In essence, the two fingering instabilities interact with each other and our results indicate that the initial source configuration modulates the strength of that interaction. The last row of Figure 4.3 shows the combined effects of R< 0 and high Ra in hindering mixing. Higher values of density contrast are known to inhibit diffusive mass transfer compared to lower values. For a/b= 3, gravitational fingering occurs at R≥ 0 but is absent at R=− 3 because at R=− 3 the velocity magnitude inside the plume is below the threshold for generating the instability; see the velocity field in Figure 4.4, second row. This important result has been overlooked in some of the prior studies of buoyant mixing that neglect viscosity contrast between the solute and the ambient fluid. When the initial source is vertical or perpendicular to the mean flow direction, i.e. a/b= 1/3, the low mobility of the plume caused by R< 0 (see Figure 4.4, second row) retards the sinking of plume caused by Ra= 9000, thus preventing the solute from reaching the bottom boundary. A close inspection of the last row in Figure 4.3 for the case of a/b= 1/3 reveals the emergence of backward propagating fingers on the rear interface. This occurs because for R< 0, the rear interface is viscously unstable and the front interface is viscously stable. However, due to the low value of the inlet flux, f L = 0.001, the backward propagating fingers do not grow to become dominant features of the flow. Next, we examine the impact of the inlet flux (or the background flux), namely f L , on both the transport dynamics and the velocity field of the plume. The inlet flux is expected to play a key role in controlling the spreading and mixing behaviour of the plume because it represents a source 77 a/b= 1/3 a/b= 3 R= 3 Ra= 1000 R=− 3 Ra= 9000 Figure 4.4: Maps of the velocity magnitude,|v|, at t = 1 for f L = 0.001. Arrow orientation in- dicates the flow direction and arrow length indicates the magnitude of the velocity vector. The colorbars represent the logarithm of base 10 of the velocity magnitude, log 10 |v|, where the white areas correspond to zones with velocity equal to the mean flow velocity, red-orange areas have velocities faster than the mean, and blue areas have velocities slower than the mean. The black line outlines the position of the plume. of mechanical energy directly inserted into the system. We want to investigate how our findings, discussed in Figure 4.3, change when f L is higher. To investigate this impact, we performed a set of simulations using the same parameter values (R, Ra and a/b) for an inlet flux that is ten times higher than the one in Figure 4.3, i.e., we increased f L from 0.001 to 0.01. As a consequence, the mean flow velocity is ten times higher, and we discuss the resulting eight plumes in Figure 4.5 at dimensionless time t = 0.1 (which corresponds to a tenth of the time considered in Figure 4.3). By doing this, when we compare the results from t = 1 for f L = 0.001 with results from t = 0.1 for f L = 0.01, we are comparing plumes whose center of mass would be in the same position if no sources of disorder (e.g., viscosity and density contrast) were acting on them. As illustrated in the first row in Figure 4.5 ( R= 3 and Ra= 1000), an increase in the hori- zontal flow velocity causes viscous fingering to dominate over density-driven fingering (compare Figures 4.3 and 4.5). We observe that the vertically oriented initial source (a/b= 1/3) provides sufficiently long interface for the viscous fingers to nucleate and grow in the downstream direction. The first row of Figure 4.5 reveals that the effect of density becomes limited to the bottom finger, 78 a/b= 1/3 a/b= 3 R= 3 Ra= 1000 R= 0 Ra= 1000 R= 0 Ra= 9000 R=− 3 Ra= 9000 Figure 4.5: Concentration fields for the faster inlet flux case, f L = 0.01, at t= 0.1. Analogously to Figure 4.3, results are showed for the vertically and horizontally extended source zones (left and right column respectively). The two top rows show results for Ra= 1000, while the two bottom rows present results for Ra= 9000. The two center rows shows results in absence of viscosity contrast, while the top and bottom row presents cases where the solute is respectively less and more viscous than the ambient fluid. 79 which is thicker, horizontally slower, more concentrated with the denser solute and therefore more prone to sinking toward the bottom boundary; see the corresponding velocity magnitude map in Figure 4.6. When a/b= 3, the horizontal source orientation does not provide sufficient vertical interface for viscous fingers to form. Hence, at R= 3, the plume travels as a single finger with high velocity at the finger tip as shown in the kinematic analysis of Figure 4.6. This causes a rapid horizontal stretching of the plume, which triggers gravitational instability on the bottom interface and away from the viscous finger’s tip. This is a novel mechanism arising out of interplay between viscosity and density-driven instabilities; viscous fingering triggers density fingering, however they dominate different parts of the plume, which are determined by the velocity field. The initial shape of the source (horizontal vs. vertical) modulates the strength of the interplay; a horizontal source (a/b= 3) enhances the effect of gravity fingering compared to a vertical source ( a/b= 1/3). This is further supported by the kinematic analysis associated with Figure 4.6. The second and third rows of Figure 4.5 show the solute plumes in the absence of viscosity contrast (R= 0). We observe that an increase in the mean flow velocity, which is horizontally oriented in our study, significantly decreases the effect of density contrast on mixing even in ab- sence of viscous fingering. An increase in Ra from 1000 to 9000 (figure 4.5, second and third rows respectively) does not affect the shape of the plume or, more importantly, the degree of mixing (as seen from the colour scale of the two plumes), for either a/b= 1/3 nor a/b= 3. The flow is predominantly advective due to the higher inlet flux; compare Figures 4.3 and 4.5. When R< 0 (bottom row in Figure 4.5), we observe that a larger value of f L impacts the com- petition between viscous and gravitational fingering; compare with the bottom row of Figure 4.3. The results in Figure 4.5 show that viscous fingering becomes dominant over density-driven finger- ing at a higher inlet flux f L . For a/b= 1/3, the initial source orientation paired with the increased horizontal velocity allows the growth of the backward viscous fingers that could not grow for the same values of R, Ra and a/b in Figure 4.3. As depicted in the last row of Figure 4.6 for R< 0, the low velocity magnitude inside the plume affects the mixing behaviour also in the case for a/b= 3, when viscous fingers are not present due to the source orientation. In this case, the low mobility 80 a/b= 1/3 a/b= 3 R= 3 Ra= 1000 R=− 3 Ra= 9000 Figure 4.6: Maps of the velocity magnitude, |v|, at t = 0.1, for f L = 0.01. Arrow orientation indicates the flow direction and arrow length indicates the magnitude of the velocity vector. The colorbars represent the logarithm of base 10 of the velocity magnitude, i.e. log 10 |v|, where the white areas correspond to zones with velocity equal to the mean flow velocity, red-orange areas indicate velocities faster than the mean, and blue areas indicate velocities slower than the mean. The black line outlines the position of the plume. of the plume prevents it from traveling with the mean flow, thus obstructing spreading. The plume shape resembles the shape of the plumes observed for the cases with R= 0, however it is asymmet- ric between downstream and upstream edges of the plume due to the viscosity contrast between the two fluids. A “wake” emerges downstream of the plume reminiscent of turbulent flow around a cavity or a bluff body (Wu, 1972). This effect, paired with a high Ra that hinders diffusion, causes the plume to be poorly mixed with the ambient fluid. 4.5.2 Mixing and spreading 4.5.2.1 Fluid interface dynamics We are interested in relating the plume’s interface to the concentration spatial variance for different initial plume configurations. The concentration variance σ 2 c can be viewed as a plume mixing indicator (Boso et al., 2013a; Bonazzi et al., 2022). The concentration spatial variance is defined as σ 2 c (t)= 1 Ω(t) Z Ω(t) [c(x,t)−⟨ c(x,t)⟩] 2 dΩ, (4.14) 81 Figure 4.7: Scatterplots of the solute plume’sη andσ 2 c for f L = 0.001 at ten equidistributed time intervals between t= 0.1 (dark colour shade) and t= 1 (light colour shade). Results for Ra= 1000 (top row) and Ra= 9000 (bottom row), and for R=− 3,0 and 3 (left, center and right column respectively). Figure 4.8: Scatterplots of the solute plume’s η and σ 2 c for f L = 0.01 at equidistributed time intervals between t = 0.025 (dark colour shade) and t = 0.5 (light colour shade), or until the time the plume reaches the right boundary of the domain. Results for Ra= 1000 (top row) and Ra= 9000 (bottom row), and for R=− 3,0 and 3 (left, center and right column respectively). 82 where Ω(t) is the area occupied by the solute plume at time t. It is defined as the area where concentration is higher than a threshold value c i , with c i = 10 − 3 , and⟨c(x,t)⟩ is the plume’s spatial mean concentration, defined as: ⟨c(x,t)⟩= 1 Ω(t) Z Ω(t) c(x,t)dΩ. (4.15) Since the solute is neither injected nor withdrawn from the domain during the time period of analysis, we can obtain the evolution equation of σ 2 c from the advection-dispersion equation 4.6 following the steps described in Jha et al. (2011b), dσ 2 c dt =− 2⟨ε⟩, (4.16) which impliesσ 2 c decreases and mixing increases monotonically with time for all parameter values studied here and the rate at which concentration variance decreases is proportional to the mean scalar dissipation rate⟨ε⟩. To better quantify the effect of shape of the initial source on mixing and spreading, we define a plume shape metricη as follows: η(t)= ˆ η(t) η o , (4.17) where ˆ η is calculated as: ˆ η(t)= ℓ c i (t) A c i (t) , (4.18) where ℓ c i is the length of the isoline of concentration c= c i , andA c i is the area of the plume delimited by such isoline. In this study, we set c i = 1× 10 − 3 . We normalize ˆ η by its initial value given byη o ≡ ˆ η(t 0 )=(ℓ c i /A c i )| t=t 0 . The shape metric η is defined to account for the length of the interface between the solute plume and the ambient fluid, and for the area occupied by the plume. Since the interface is where mixing happens, and the plume area is related to spreading, the definition of η includes mecha- nisms of both mixing and spreading. For diffusion-dominated transport of an isotropic (circular) 83 source, e.g., R= 0, Ra≤ 1000 and a/b= 1, the plume radius increases proportionally to √ t, the plume interface isℓ c i ∝ √ t, and the area isA c i ∝ t (Jha et al., 2011a), rendering η ∝ t − 1/2 , i.e., η decreases monotonically in time during diffusive transports. For diffusive systems,σ 2 c also de- cays with time such thatη∝σ 2 c , with the proportionality factor given by the initial source shape. In convection-dominated transports due to viscous fingering or permeability heterogeneity, η can increase at early time due to rapid stretching of the plume interface while the plume area stays relatively constant or grows slower than the interface length. This would imply that mixing is still in its early stages (σ 2 c is relatively high) and the solute has not mixed with the ambient fluid enough to occupy an area significantly larger than the initial one, but mixing will increase rapidly in the future across the stretched interface as implied by highη. As time progresses and the solute becomes more mixed with the ambient fluid, we expect η to decrease as the area of the plume increases and its interface becomes less sinuous. The result is thatη evolves non-monotonically in convective systems. Since the concentration spatial varianceσ 2 c quantifies the degree of mixing (high σ 2 c means lack of mixing, see Bonazzi et al., 2021), we analyze the co-evolution ofη andσ 2 c in Figures 4.7 and 4.8 to learn how the interplay of viscous and gravitational fingering affects spreading and mixing of plumes of different initial shapes. In the sense of nonlinear dynamical systems, theη vs. σ 2 c plot can be seen as a phase space diagram for investigating the effect of fingering and shape parameters on mixing. This is equivalent toσ 2 c vs. ⟨ε⟩ phase diagram in viscous fingering (Jha, 2014 b). We first analyze the evolution behaviour for a diffusive system, i.e., [R,Ra]=[0,1000], to set a baseline for interpreting these plots. For a diffusive system, bothη andσ 2 c decrease monotonically in time, which can be inferred from Figures 4.7 and 4.8 by noting that darker shades denote early time and lighter shades denote late times. Therefore,η vs. σ 2 c trend is also monotonic and agrees with our previously mentioned hypothesis for diffusive transport of a circular source. The effect of vertical or horizontal initial shape onη vs. σ 2 c trend is minimal in diffusive systems except in the presence of boundary condition effects; the vertical source (a/b= 1/3) at slower mean flow ( f L = 0.001) reaches the bottom boundary fast and experiences the effect of no-flux boundary condition, which 84 causes a rapid drop inη. At higher flux f L = 0.01 (figure 4.8), this effect is absent and we recover the baseline trend in η vs. σ 2 c . In convective systems with fingering instabilities, η behaves non-linearly with σ 2 c because of significant changes in the plume shape, e.g., horizontal viscous fingering in Figure 4.8 for [R,Ra,a/b]=[3,1000,3]. The effect of R on plume shape (η) and mixing (σ 2 c ) metrics is amplified when convection dominates the flow, which happens at higher Ra (second row in Figure 4.7) and higher f L (both rows in Figure 4.8). When both density and viscosity-driven fingering occurs, their effects are superposed on the baseline diffusive trend, rendering theη vs. σ 2 c trend complex. In these cases, the initial source shape (a/b) plays a more important role in determining theη vs. σ 2 c trend than in the cases where fingering is absent. The importance of relating an indicator of mixing, such asσ 2 c , to a plume shape metric, such as η, stems from the observation that mixing is more intense where the concentration gradients are higher, i.e., at the interface between an injected solute and the ambient fluid. Engineers and hydrogeologists who may need to maximize or minimize mixing between two fluids can attempt to control the plume shape evolution in time by leaning toward one initial solute source orientation or another. For example, if a fluid is being injected in an aquifer for the purpose of groundwater remediation, it is desirable to ensure maximum mixing between the two fluids. Traditionally, this implies that the remediation fluid selection and its injection are designed in a way to maximize the interface between the two fluids; in other words, a high η has been traditionally associated with good mixing, or low values ofσ 2 c . We have shown in Figures 4.7 and 4.8 that, depending on several factors (the mean flow velocity, the initial shape and orientation of the injected solute source, and the degree of viscosity and density contrast between the two fluids), the same degree of mixing can be achieved for different values of the shape parameterη, that is, for different shapes of the solute plume. Similarly, the same value ofη can be associated with different degrees of mixing; see, for example, cases for R=− 3 and Ra= 9000 in Figure 4.8 for a/b= 3 and a/b= 1/3. 85 a/b= 1/3 a/b= 1 a/b= 3 ⟨c⟩ σ 2 c Figure 4.9: Contour plots of the spatial statistics of the concentration. The mean concentration (top row) and concentration variance (bottom row) of the plume at t= 1 for f L = 0.001 and for the range of R and Ra investigated for different initial shapes: a/b= 1/3, 1 and 3. a/b= 1/3 a/b= 1 a/b= 3 ⟨c⟩ σ 2 c Figure 4.10: Contour plots of the spatial statistics of the concentration. The mean concentration (top row) and concentration variance (bottom row) of the plume at t = 0.1 for f L = 0.01 and for the range of R and Ra investigated for a/b= 1/3, 1 and 3. 86 4.5.2.2 Mixing analysis in the R− Ra parameter space and concentration statistics To understand how mixing evolves globally in the domain as a function of the density and viscosity contrasts, we plot the spatially-averaged mean⟨c⟩ and variance σ 2 c of the concentration field as contour maps in the parameter space of R and Ra for f L = 0.001 and f L = 0.01 and distinct values for a/b (see Figures 4.9 and 4.10 respectively). The spatial statistics of the concentration field can be viewed as global measures of mixing. As discussed in Bonazzi et al. (2022), the spatial mean concentration is related to the dilution index (Kitanidis, 1994) and therefore can be used as a descriptor for the global mixing state of a solute plume. Similarly, the concentration variance is related to the global scalar dissipation rate (Kapoor & Gelhar, 1994). Both Figures present results for three different values of the initial solute source shape (a/b= 1/3, 1 and 3), presenting the contour maps for the concentration spatial mean (top row of Figures 4.9 and 4.10), and variance (bottom row of Figures 4.9 and 4.10). Results for⟨c⟩ and σ 2 c are shown at time snapshots t = 1 (figure 4.9) and t= 0.1 (figure 4.10). Focusing on the low inlet flux case reported in Figure 4.9, we can see that, for R< 0, a vertical initial source (a/b= 1/3, left column) leads to better mixing compared to a horizontal source, in agreement with the analysis carried out for the bottom row images in Figure 4.3. For R> 0, however, a better degree of global mixing (lower⟨c⟩ andσ 2 c ) is achieved for a/b= 3 because the horizontal initial source facilitates the creation of gravity fingers along with viscous fingers. For R≈ 0, values of⟨c⟩ andσ 2 c lie between the ones computed for R> 0 and R< 0. We also observe that if a/b= 1/3, the impact of viscosity contrast between the fluids on the resulting concentration spatial mean and variance is not as strong as in the case of a/b= 3, in which the difference between the upper and the lower part of the parameter space is more noticeable. In regard to the effect of Ra, the results in Figure 4.9 show that for a fixed R, better mixing is achieved for lower values of Ra. This is due to the fact that high density contrasts lead to less diffusive mass transfer. However, the impact of Ra on mixing becomes almost negligible for R> 0: note how on the upper half of the parameter spaces in Figure 4.9, the concentration isolines are almost horizontal. A similar isoline orientation is found for R< 0 when Ra≳ 6000, suggesting that as Ra increases its impact on 87 mixing tends to “stabilize”. The reason is that reduction in mixing due to lower diffusion at higher Ra is balanced by increase in mixing due to a longer interface at higher Ra; a similar observation has been made for fluid mixing from viscous fingering (Jha et al., 2011a). Further investigation involving a broader range of Ra considered would be warranted to support this hypothesis. As expected, the impact of Ra on mixing is higher for a/b= 3 because a horizontally oriented initial source allows gravitational fingering. Next, we analyze Figure 4.10 to understand how the physics of mixing changes when the inlet flux is increased tenfold ( f L = 0.01). The observation that a vertically extended initial source translates to better mixing, compared to a horizontal source, for R< 0 is still valid. However, unlike what was observed in Figure 4.9, when the mean flow velocity increases, for R> 0, a plume originating from a horizontally extended initial source (a/b= 3) does not correspond to the best mixed scenario. Since a higher flow velocity means that viscous fingers are clearly forming in the direction of flow, a/b= 1/3 (perpendicular to the flow direction) is the shape that offers more interface and favors finger formation. As seen in Figure 4.5, for R> 0 and a/b= 3, the plume travels as a single finger, which inhibits convective mixing. Another difference in regard to the discussion of Figure 4.9 is that in Figure 4.10, values of R≈ 0 do not result in a mean concentration that lies in between the ones observed for R> 0 and R< 0. For R≈ 0 we have the highest computed values of⟨c⟩, and this trend is exacerbated for a/b= 3. The reason is that a high f L value enhances viscous fingering of the vertical source, forward fingering for R> 0 and backward fingering for R< 0. Fingering-induced mixing leads to a drop in⟨c⟩. Fingering is absent at R≈ 0, and hence we observe highest⟨c⟩ values at R≈ 0. As mentioned previously, at R= 0, mixing is not significantly affected by the density contrast, whose effect on solute transport is being hampered by the higher mean flow velocity. The solute plume travels along the mean flow direction almost undisturbed. On the other hand, the lower velocity inside the plume due to R< 0 causes downstream stretching of the plume and appearance of the previously mentioned “wake” (see Section 4.5.1, discussion for Figure 4.5), thus resulting in a mean concentration lower than the one computed with R= 0. The difference in⟨c⟩ between R= 0 and R< 0 is not as large as the 88 one between R= 0 and R> 0, where the formation of a single finger that travels faster than the mean flow velocity is sufficient to achieve the best degree of mixing among the ones obtained for a/b= 3. The effect of Ra on both⟨c⟩ andσ 2 c is almost negligible for R> 0, similar to what was observed in Figure 4.9. Recall that these observations are at a specific timestep; Ra does influence the time evolution of mixing, e.g., the evolution ofη in Figure 4.8. If a/b= 1/3, density contrast does not seem to have a significant impact on the mixing metrics considered here for any value of R. For a/b= 3 on the other hand, if R< 0, we observe again that density differences between the fluids seem to have an impact only for Ra≲ 6000, analogous to the case of lower f L in Figure 4.9. 4.6 Effects of heterogeneity in the permeability field Next, we focus on quantifying the effect of permeability heterogeneity on plume mixing and spreading under hydrodynamic instabilities. In order to account for heterogeneity, we consider a stratified porous medium. Stratification or layering is a common feature in many geophysical flows. For example, a layered geological formation leads to depth-dependence in hydraulic prop- erties, e.g., causing the permeability to change significantly from a clean sandstone layer of high permeability to an adjacent silty sandstone layer of low permeability. This is relevant for solute transport because it is expected that the solute plume will move faster and farther in a high per- meability layer compared to a low permeability layer. Therefore, it is intuitive to expect enhanced plume stretching, sharper concentration gradients along the stretched interface, and faster mixing across the fluid interface. We model stratification following the model described in Equation 2.9 in Nijjer et al. (2019): ln(k)=− σ 2 st cos(2πnz)− ln(I 0 (σ 2 st )), (4.19) 89 where σ 2 st is the log-permeability variance, z is the vertical coordinate, I 0 is the modified Bessel function of the first kind, and n is the wavenumber which determines the length scale of the per- meability variations. In this study, we set n= 2, thus obtaining four layers in the domain. We sim- ulate transport for two different values of the log-permeability variance: σ 2 st = 0.05 andσ 2 st = 0.1, with the same average permeability value. We also consider the baseline case of a homogeneous medium (σ 2 st = 0) with a uniform permeability value equal to the average permeability in the het- erogeneous cases. Regarding fluid properties, we set R= 2 (the plume is less viscous than ambient fluid) and Ra= 5000 (the plume is denser than the ambient fluid). We also present results for R= 0 for comparison. For the upcoming results, we set the inlet flux parameter to f L = 0.01 and con- sider two values of the initial source configuration: a/b= 1/3 (vertical orientation) and a/b= 3 (horizontal orientation). Figure 4.11 shows the concentration fields for different levels of heterogeneity and values of R. Results are depicted at dimensionless time t= 0.1 for a/b= 1/3. Increasing heterogeneity in a stratified medium results in the formation of preferential flow paths in higher permeability layers. As σ 2 st increases, the middle layer becomes more permeable, which causes additional stretching of the plume as seen in Figure 4.11. The presence of heterogeneity changes the location of the peak concentration in the plume. For a homogeneous case (σ 2 st = 0), highest concentration is at the bottom of the plume because of gravitational sinking (Ra = 5000). For layered cases, the highest concentration occurs in the high permeability middle layer. This occurs with or without the viscosity contrast, however, viscous fingering at R = 2 creates additional mechanisms that are superposed on the effects due to density contrast and layering. We know from the literature (Sajjadi & Azaiez, 2013; Nicolaides et al., 2015; Nijjer et al., 2019; Bonazzi et al., 2021) that viscous fingering interacts with permeability heterogeneity. For example, during continuous slug injection in layered media (Nijjer et al., 2019), viscous fingering dominates the mixing behavior at early times, and when the fingers have grown to become comparable to the layer thickness scale, heterogeneity begins to dominate mixing. Such investigations of the interplay between fingering and heterogeneity mechanisms are missing for finite-size solute bodies, which we address here 90 R= 2 R= 0 σ 2 st = 0 σ 2 st = 0.05 σ 2 st = 0.1 Figure 4.11: Concentration fields at t= 0.1 for a/b= 1/3, for three different levels of heterogene- ity: no heterogeneity (σ 2 st = 0, top row),σ 2 st = 0.05 andσ 2 st = 0.1 (center and bottom row respec- tively). On the left column are the results for R= 2, while on the right columns are presented the corresponding cases in absence of viscosity contrast (R= 0). All results are for Ra= 5000 and f L = 0.01. The permeability field has been superimposed in transparency to the concentration fields; light shades correspond to higher permeability layers, while darker shades to low perme- ability ones. 91 R= 2 R= 0 σ 2 st = 0 σ 2 st = 0.05 σ 2 st = 0.1 Figure 4.12: Maps of the velocity magnitude, |v| at t = 0.1 for f L = 0.01 and a/b= 1/3. The arrows indicate the velocity vector with both direction and magnitude information. The colorbars represent the logarithm of velocity magnitude, i.e. log 10 |v|. The black line outlines the position of the plume. All results are for Ra= 5000. (see Figure 4.11 and Figure 4.12). The high flow velocity in the middle layer (Figure 4.11, bottom left image) causes tip splitting of the fingers followed by tip blunting, which we explained in our prior work; see Figure 3 in Bonazzi et al. (2021). The result of these interactions is rapid mixing in the high permeability middle layer. At R= 0, these effects are absent, and high concentrations in the middle layer persist for a longer duration. When the initial solute source is horizontally oriented (a/b= 3, Figure 4.13), permeability stratification does not have a significant impact on solute transport in absence of viscous fingering (R = 0, right column). This is because the source zone is not crossing multiple layers at the considered time. However, when viscosity contrast is present (R= 2, left column of Figure 4.13), stratification favors the formation of viscous fingers in the horizontal direction, a feature that was not observed for analogous cases in homogeneous media discussed in Section 4.5. The horizontal 92 R= 2 R= 0 σ 2 st = 0 σ 2 st = 0.05 σ 2 st = 0.1 Figure 4.13: Concentration fields for a/b= 3 at t = 0.1 for three levels of heterogeneity: σ 2 st = 0 (no heterogeneity, top row),σ 2 st = 0.05 (second row), andσ 2 st = 0.1 (third row). All results are for Ra= 5000 and f L = 0.01. 93 R= 2 R= 0 σ 2 st = 0 σ 2 st = 0.05 σ 2 st = 0.1 Figure 4.14: Maps of the velocity magnitude, |v| for a/b = 3 at t = 0.1 for f L = 0.01. The arrows indicate the velocity vector with both direction and magnitude information. The colorbars represent the logarithm of velocity magnitude, i.e., log 10 |v|. The black line outlines the position of the plume. All results are for Ra= 5000. finger formation can be explained by observing Figure 4.14, left column. As the difference in velocity magnitude between layers increases, the upper and lower edges of the plume experience lower velocity magnitudes compared to the center. The horizontal direction of the velocity field causes the center of the plume to travel faster than its edges and split into two fingers, following the tip splitting mechanism. Atσ 2 st = 0.05 (second row, left column of Figure 4.14), the upper finger eventually merges into the dominating lower finger, following the finger coarsening mechanism (Jha et al., 2011b). At σ 2 st = 0.1, heterogeneity is strong enough to sustain the growth of both fingers for a longer time (last row, left column of Figure 4.13). Eventually, the fingers will merge as in the case of a vertically oriented source (a/b= 1/3), and permeability heterogeneity will be the dominant agent affecting solute mixing and transport. 94 4.7 Summary The need for modelling the effects of viscosity and density contrasts in the spreading and mixing dynamics of a solute body has been recognized and addressed in the literature given its importance in contaminant transport, energy recovery and storage in the subsurface environment. In all these applications, solutes are characterized by viscosity and density values that differ from those of the ambient fluid. Furthermore, these solutes are injected into the subsurface through source zones of different spatial configurations. Our work focuses on systematically investigating the impact of the initial configuration of the source zone on the transport behaviour in the presence of hydrodynamic instabilities. The numerical analysis carried out in the present Chapter illustrates the importance of both the initial solute source configuration and the intensity of the inlet (background) flux in controlling the relative significance of R and Ra in the mixing and spreading behaviour. In the following, we highlight the key findings of our analysis. Effect of viscosity contrast on mixing: For a homogeneous medium and low inlet or background flux , gravitational fingering is impeded when the solute is more viscous than the ambient fluid, i.e., R< 0. When R> 0, viscous fingering occurs in both horizontal and vertical directions and in the latter case, viscous fingering and gravitational fingering effects overlap. For R< 0, better mixing is achieved with a vertically oriented initial source. If more mixing is desired under R> 0 conditions, a horizontally oriented solute source is preferable at a low background flux, whereas a vertical source is preferable at high background flux. At higher inlet flux , when R> 0, a horizontally oriented source does not result in the best degree of mixing because the dominance of viscous fingering (as the primary driver of mixing) allows a vertically extended source to form a larger number of fingers. The importance of viscosity contrast in enhancing mixing when the inlet flux is higher means that the degree of mixing is lower when no viscosity contrast is present (R≈ 0). Effect of density contrast on mixing: For a homogeneous medium and low inlet flux , higher values of Ra were also found to allow lower rates of diffusive mass transfer. The impact of Ra 95 on mixing was found to be almost negligible when R> 0. A horizontally oriented initial source zone allows gravitational fingers to form and grow in the vertical direction. The number of fingers increases with Ra in accordance with the linear stability analysis theory. At higher inlet flux , the role of Ra on mixing is minor for a vertically oriented source for any value of R considered. For a horizontally oriented source, Ra appears to affect mixing when R< 0. Effect of background flow velocity on mixing: We established that a higher mean flow veloc- ity hampers gravitational fingering, thus rendering viscous fingering the dominant force driving mixing. As a consequence, the source orientation becomes a significant factor when R< 0. Un- der these conditions, we observed the following: (1) for a/b= 1/3 (vertically oriented source), a higher flow velocity allows the formation of viscous fingers, which will lead to enhanced mixing; while (2) for a/b= 3 (horizontally oriented source) the absence of viscous fingering and the low velocity inside the plume result in poor mixing. Relationship between mixing and plume shape: The numerical analysis presented here allowed us to investigate the relationship between the degree of mixing and the length of the interface between the solute and the ambient fluid (i.e., plume shape). Traditionally, one would expect en- hanced mixing rates for increased elongation of the plume’s interface with the surrounding fluid. In the context of aquifer remediation, hydrologists and engineers aim to maximize the interface between the contaminant and groundwater in order to enhance dilution and reaction rates (Bagt- zoglou & Oates, 2007; Kitanidis & McCarty, 2012; Trefry et al., 2012; Piscopo et al., 2013). As shown in Section 4.5.2.1, the relationship between mixing and plume shape is non-trivial when fingering instabilities are present. The results in Section 4.5.2.1 shed novel insights on the role of source orientation, inlet flux and fingering instabilities (i.e., viscosity and density contrasts) in controlling the plume shape and its mixing dynamics. Our results show that the same degree of mixing can be achieved for different plume shapes. 96 Effect of permeability heterogeneity on mixing: To account for the heterogeneity, we considered a stratified permeability field. We found that, for a vertically oriented solute source, heterogeneity in the porous medium affects the location of the plume’s highest concentration values by over- laying the effect of gravitational instability and favoring higher concentration values in the higher permeability layer. The high flow velocity in the more permeable layer causes tip splitting and results in rapid mixing. For a horizontally oriented solute source, the combined effect of viscosity contrast and permeability heterogeneity results in fingering along the flow direction, a result not observed for the homogeneous media analysis presented in Section 4.5. The simulations carried out in this Chapter highlight the importance of the initial solute source configuration and the intensity of the inlet flux in modulating both the plume’s spreading and mixing behaviour in homogeneous and heterogeneous permeability fields. Additional research is needed to understand these complex interactions in 3D settings, larger computational domains, in absence of interaction with the bottom boundary of the domain, with different injection modes and distinct conceptualizations of the permeability field, including higher levels of permeability heterogeneity. 97 Chapter 5 Uncertainty quantification in the presence of viscous fingering and permeability heterogeneity 5.1 Introduction Obtaining deterministic predictions of the spatiotemporal dynamics of a solute body in the sub- surface environment is an elusive goal given limited site characterization data. Due to the high cost associated with data acquisition, in-field measurements are scarce, and thus a detailed charac- terization of the spatial heterogeneity of hydrogeological properties (such as the permeability) is impractical. Because of this lack of information on all relevant scales, flow and transport models are subjected to uncertainty and the stochastic paradigm is invoked. Under the stochastic ap- proach, the permeability field is conceptualized as a random space function (RSF) (Rubin, 2003; Fiori et al., 2019) and environmental performance metrics (EPM) relevant for risk analysis, such as solute travel times and concentrations, are characterized in terms of their probability density function (PDF). Strategies for propagating the uncertainty in the permeability field to solute travel time and concentration predictions have been discussed in the literature (for an extensive review, see Fiori et al., 2015; Linde et al., 2017, and references therein). Analytical methods, typically based on per- turbation theory, have been applied to compute the statistical moments of solute travel time (Rubin 98 & Dagan, 1992; Cvetkovic et al., 1992; Sanchez-Vila & Guadagnini, 2005) and resident concen- tration (Kapoor & Kitanidis, 1998; Andriˇ cevi´ c, 1998; Fiori & Dagan, 2000; Tonina & Bellin, 2008; de Barros et al., 2022). The (semi)-analytical solutions derived in these studies are typically lim- ited to low levels of heterogeneity in the log-permeability field, e.g., log-permeability variance σ 2 f ≲ 1. To estimate uncertainty in EPMs in aquifers displaying highly heterogeneity, numerical Monte Carlo is often employed (see for example, Moslehi & de Barros, 2017; Libera et al., 2019). Semi-analytical methods that allow to obtain an estimate of solute arrival times have also been developed (Cvetkovic et al., 2016; de Barros et al., 2016; Fiori et al., 2017). Other studies, based on analytical or numerical Monte Carlo approaches, focused on the estimation of the one-point concentration PDF as a function of the RSF properties of the log-permeability field (Dentz, 2012; Cirpka et al., 2011; de Barros & Fiori, 2014; Boso & Tartakovsky, 2016). In absence of local- scale dispersion, Dagan (1982) showed that the concentration PDF is binary, that is, either zero or equal to the initial concentration at the solute source location. Subsequent studies showed that the concentration PDF becomes unimodal if the sampling volume is increased and the effects of local-scale dispersion are taken into consideration (Bellin et al., 1994; Caroni & Fiorotto, 2005; Schwede et al., 2008). For low levels of heterogeneity of the permeability field, the concentration PDF can be approximated by theβ-distribution model as originally suggested in Fiori (2001) (see also Fiorotto & Caroni, 2002; Bellin & Tonina, 2007; Bonazzi et al., 2022). Methods for com- puting the concentration statistics for higher levels of heterogeneity are also available (see Meyer et al., 2010). The above-mentioned studies (amongst many others) have provided insight into the propaga- tion of uncertainty, with respect to scales of fluctuation of the random permeability, on solute travel times and concentration. Given a few exceptions (Welty & Gelhar, 1991; Talon et al., 2004), most of these studies on uncertainty analysis neglect the effects of viscosity contrast in subsurface systems with multiple fluids, e.g., ambient groundwater and an injected solute. In many applica- tions, groundwater contaminants (i.e., nonhalogenated semivolatile compounds and jet fuel) are 99 more viscous than water. Flowers & Hunt (2007) provide a series of examples in environmen- tal engineering and hydrology where viscosity contrast between the ambient groundwater and a contaminant can be found. However, incorporating the effects of viscosity contrast in numerical simulation models can be challenging due to high computational costs associated with the non- linearity of the governing equations. Viscosity contrast between the ambient fluid and the solute can lead to the hydrodynamic instability known as the Saffman-Taylor instability (Saffman, 1986), with formation of viscous fingers at the unstable interface. Due to the viscosity dependence on solute concentration, physical properties of the fluids vary in space and time, thus affecting the flow field (Jha et al., 2011a). In the presence of permeability heterogeneity, viscous fingering can enhance the mobility of a solute plume as shown in Bonazzi et al. (2021). Despite the extensive literature on viscous fingering (Tran & Jha, 2020; Van der Meer, 1993; Christie, 1989; Tchelepi & Orr Jr, 1994; Tan & Homsy, 1988; De Wit et al., 2005; Nicolaides et al., 2015; Bonazzi et al., 2021), additional efforts are needed to investigate how viscosity contrast impact the uncertainty estimates of EPMs when the flow field is spatially random. Given the impact of viscosity contrast on mixing (Bonazzi et al., 2021; Nicolaides et al., 2015; Jha et al., 2011b,a) and early arrival times, and considering that permeability heterogeneity, as previously discussed, heavily affects flow and solute transport, we wish to quantify in this Chapter the uncertainty related to both the concentration distribution and the early arrival times at a target location, and how these estimates are affected by both the level of viscosity contrast and the degree of permeability heterogeneity. 5.2 Methods 5.2.1 Multi-Gaussian Log-Permeability Field The structure of the log-permeability field, denoted here by f(x,ω)≡ ln[k(x,ω)], is modelled through a Random Space Function (RSF) Kitanidis (1997); Rubin (2003) with x∈D andω∈Ω (the probability space). For the purpose of illustration, Y(x,ω) is represented as a multi-Gaussian 100 spatial random field and assumed to be statistically stationary such that it can be fully characterized by its mean⟨ f(x)⟩ and the correlation function of its fluctuations f ′ (x,ω)= f(x,ω)−⟨ f(x,ω)⟩, which is provided byC f (x,x ′ ) = ⟨ f ′ (x,ω) f ′ (x ′ ,ω)⟩. The variance of the f -field is defined as C f (0)=σ 2 f . In this work, we adopt an isotropic exponential covariance model: C f (x,x ′ )=σ 2 f exp   − v u u t 2 ∑ j=1 (x j − x ′ j ) 2 λ 2   , (5.1) whereλ represents the correlation length. 5.2.2 Mathematical Formulation of Flow and Transport Analogously to Chapter 3, we consider a two-dimensional (2D) porous medium described by the Cartesian coordinate system x=(x 1 ,x 2 ), and we consider a domain with dimensions L 1 × L 2 . The hydraulic conductivity of the porous medium, K(x), is set to be locally isotropic and heteroge- neous. The physical equations describing the phenomenon we are investigating are the same we dis- cussed in Chapter 3. Following the dimensionless groups reported in Appendix A, the dimension- less equations we implemented in our numerical model are here presented. Assuming incompressible flow and the absence of sink/source terms, the governing equation for the subsurface flow field is given by: ∇ x · q(x,t)= 0, (5.2) where q(x,t) is the specific discharge from Darcy’s law: q(x,t)=− k(x) µ(c(x,t)) ∇ x p(x,t) (5.3) 101 in which p(x,t) is the pressure field, k(x) is the permeability field, and µ(c(x,t)) is the viscosity as a function of the concentration of solute in the ambient fluid, c(x,t), expressed through the exponential viscosity model: µ(c(x,t))= e − Rc(x,t) (5.4) where R is the log-viscosity ratio between the ambient fluid viscosity µ 0 and the solute viscosity µ 1 , R= ln(µ 0 /µ 1 ). The governing equation for solute transport is the classic advection dispersion equation: ∂c(x,t) ∂t + q(x,t) φ · ∇ x c(x,t)= 1 Pe ∇ 2 x c(x,t), (5.5) whereφ is the medium porosity, and Pe represents the P´ eclet number, defined as Pe ≡ Uλ/D, in whichλ represents the characteristic length scale of the log-permeability field, U is the character- istic velocity, defined as the mean longitudinal velocity, and D is the constant local-scale dispersion coefficient. The boundary conditions for the flow problem are a no-flux condition on the top and bottom boundaries, and prescribed pressure on the left and right boundaries. For the transport problem, we imposed a no-flux condition on the top, bottom and left boundaries of the domain, while on the right boundary, analogously to the set-up presented in Chapter 4, we imposed a natural outflow boundary condition, which meant that the solute can exit the domain through convection only, following the direction of the average flow velocity. The initial condition for the transport problem is a zero concentration in the domain, as shown in Figure 5.1, with the exception of the source area. There, the solute is instantaneously injected and the concentration c 0 is set to be equal to unit value, with a small numerical perturbation on the edges to allow the development of the hy- drodynamic instability caused by the viscosity contrast between the solute and the ambient fluid. In the remaining of this Chapter, we will refer to the solute concentration normalized by its initial value, c/c 0 , as just c for simplicity, given that c 0 = 1. 102 Figure 5.1: Initial distribution of the concentration field in the considered domain. The white dashed line represents the location of the control plane. 5.2.3 Numerical implementation The flow governing equation is solved using the numerical methods described in Section 3.4 and here briefly summarized. Equation 5.2 is solved for cell-centered pressures using a finite-volume method with second-order accuracy. The Darcy flux in Equation 5.3 is linearized using the two- point flux approximation. As for the scalar transport, we use an explicit scheme for the concen- tration field in Equation 5.5. In particular, we employ a second-order finite difference method, and compute the concentration at the next time step in terms of the mass fluxes (advective and diffusive). We discretized the domain with a two-dimensional Cartesian grid with grid blocks of equal size in the x 1 and x 2 directions. Details about the chosen parameters to describe the domain and the transport problem in Equation 5.5 are reported in Table 5.1. 103 Table 5.1: Input parameters used in the simulations Parameter Symbol Value Calculated as Correlation length in the x and y-direction λ 5 m - Domain length in the x 1 -direction L 1 60 m 12λ Domain length in the x 2 -direction L 2 50 m 10λ Mesh size in the x 1 - and x 2 -directions ∆ 0.42 m λ/12 Distance of source from top and bottom boundary s tb 15 m 3λ Distance of source from left boundary s l 10 m 2λ Length of source in the x 1 -direction ℓ 1 5 m λ Length of source in the x 2 -direction ℓ 2 20 m 4λ Longitudinal distance of control plane from left boundary L cp 55 m 11λ Mean longitudinal flow velocity U 1 m/d - Local dispersivity α 0.003 m 0.008∆ Local scale dispersion coefficient D 0.03 m 2 /d α U P´ eclet number Pe 300 Uλ/D 5.2.4 Monte Carlo framework We randomly generated a large number, N r = 7000, of hydraulic conductivity fields using a ran- dom space function model as described in Section 5.2.2. For each realization, we then obtain the permeability field k(x) from the hydraulic conductivity field K(x) using the relation: k(x)= K(x) µ 0 ρ 0 g , (5.6) where µ 0 and ρ 0 are the viscosity and density of the ambient fluid (here assumed to be water), respectively, and g is gravity acceleration. We then solved the flow and transport problem in Equations 5.2-5.5 for the N r realizations. From the output of our simulations, we analyzed two EPMs: solute arrival times and concentration. We compute the early arrival times t e at a control plane, perpendicular to the direction of mean flow, located eight correlation lengths downstream from the forefront edge of the initial solute source. We define the early arrival time t e as the time when the average concentration⟨c⟩ ⊥ at the control plane is equal to 0.01. We also analyze the concentration field c(x,t) and in particular the variance of the concentration in each point of the field across the considered N r realizations, σ 2 c . Note that, unlike in the previous Chapters, here 104 σ 2 c does not denote the spatial variance of the concentration inside the plume. The concentration variance is computed at every grid block over the ensemble of size N r . 5.3 Results 5.3.1 Convergence analysis Figure 5.2 shows how the first two statistical moments of the early arrival times t e at the third control plane vary as a function of sample size N r of the permeability fields. We can see that for all values of permeability heterogeneity and viscosity contrast considered, the mean of the early arrival time at the third control plane (Figure 5.2, left) converges to its reference value (that is, the mean of t e over all 7000 Monte Carlo realizations) for N r ≲ 2000. On the other hand, when we consider the variance of the early arrival times,σ 2 t e (Figure 5.2, right), there is a distinction in the number of realization required to reach convergence, and that distinction is based on the level of heterogeneity of the log-permeability field. For the lower level of permeability heterogeneity considered (σ 2 f = 0.5), for both values of viscosity contrast convergence is reached for N r ≲ 3000, while the higher level of permeability heterogeneity (σ 2 f = 1.5) requires a higher number of re- alizations to reach convergence: N r ≈ 5000 for R= 1.5, and N r ≳ 6000 for R= 0. In particular, the high heterogeneity case with no viscosity contrast (σ 2 f = 1.5, R= 0) requires additional sim- ulations to reach full convergence of σ 2 t e . This is due to the fact that the presence of a viscosity contrast (R> 0) increases the mixing rates of the solute plume thus reducing the sample-to-sample fluctuations of the ensemble. This was observed in the past when examining transport at finite Pe and 0<σ 2 f ≲ 1 (Dentz & de Barros, 2013b). 105 Figure 5.2: Convergence analysis for the meanµ t e (left) and varianceσ 2 t e (right) of the early arrival time t e at the third control plane for the four sets of simulations considered in this study. The dashed grey lines represent the final value of µ t e (left) andσ 2 t e (right) at N r = 7000. 5.3.2 Uncertainty in the concentration field At a generic point in the domain with Cartesian coordinates x P =(x P 1 ,x P 2 ), we define the variance of the concentration over the N r Monte Carlo realizations, at a specific time t, as: σ 2 c (x P ,t)= 1 N+ 1 N ∑ i=1 (c i (x P ,t)−⟨ c(x P ,t)⟩) 2 , (5.7) where c i (x P ,t) is the value of concentration at x P at time t for the i th realization of the permeability field, and ⟨c(x P ,t)⟩ is the ensemble mean concentration at x P and time t over the N r realizations, defined as: ⟨c(x P ,t)⟩= 1 N N ∑ i=1 c i (x P ,t). (5.8) The computation of σ 2 c (x P ,t) in every point of the domain (i.e., x P ∈[L x × L y ]), enables us to compute the maps ofσ 2 c (Figure 5.3). A close-up of the images in the left column of Figure 5.3 is shown in Figure 5.4, to highlight the distribution ofσ 2 c at an early time in the simulations, t= 5. It is clear that, at t = 5, σ 2 c exhibits a bimodal behaviour for the lower level of permeability heterogeneity considered (σ 2 f = 0.5, first and second rows in Figure 5.3, left column), a result in ac- cordance with the observations reported in the literature (Rubin, 1991; de Barros et al., 2012b). The 106 t= 5 t= 20 R= 0 σ 2 f = 0.5 R= 1.5 σ 2 f = 0.5 R= 0 σ 2 f = 1.5 R= 1.5 σ 2 f = 1.5 Figure 5.3: Maps of the concentration varianceσ 2 c in the domain, at t= 5 (left column) and t= 20 (right column), for the low permeability heterogeneity cases (σ 2 f = 0.5, first and second row) and the high permeability heterogeneity cases (σ 2 f = 1.5, third and fourth rows), and for R= 0 (first and third rows), and R= 1.5 (second and fourth rows). The black dashed line outlines the initial location of the solute source. 107 σ 2 f = 0.5 σ 2 f = 1.5 R= 0 R= 1.5 Figure 5.4: Three-dimensional representation of the images in Figure 5.3, for t= 5. higher level of heterogeneity considered, seeσ 2 f = 1.5, at t = 5 (right column in Figure 5.4), also exhibits a bimodal behaviour, but it is less pronounced than for the lower heterogeneity case.This is due to the fact that for higher levels of heterogeneity, the leading edge of the plume is subject to more fingering and channelling, which enhance mixing. This means that there will be less sample- to-sample fluctuations in concentration, thus causing at early times σ 2 c to be lower at the leading edge of the plume compared to its value at the trailing edge. Figure 5.4 indicates that higher per- meability heterogeneity causes a higher uncertainty related to the concentration distribution closer to the initial location of the solute source (black dashed lines in Figure 5.3). This can be explained considering that higher permeability heterogeneity implies that the plume will be more spreaded, with the front traveling faster and mixing with the ambient fluid at a faster rate than the back of the plume. At the same considered time, the trailing edge of the plume has not traveled and mixed as much, thus being subjected to higher fluctuations in local concentration between different Monte Carlo Realizations. 108 At later times we can notice by comparing the higher heterogeneity cases with the lower ones, for the same level of viscosity contrast R (compare the first and third rows, and the second and fourth rows, right column of Figure 5.3), thatσ 2 f = 1.5 cases present higherσ 2 c close to the solute source respect to σ 2 f = 0.5 cases, preserving the trend discussed for earlier times. However, the overall values of σ 2 c are lower for the higher heterogeneity cases (as noticeable from the lighter color shade in the second and fourth rows, right column of Figure 5.3, compared to the first and third rows). This happens because higher permeability heterogeneity enhances mixing, thus in- creasing the dilution of the solute. The effect of viscosity contrast in enhancing dilution at later times is also visible by compar- ing the first and second rows of the right column in Figure 5.3 ( R= 0 and R= 1.5 respectively, σ 2 f = 0.5). The uncertainty related to the case with no viscosity contrast is i) higher in magnitude compare to the uncertainty related to the case with R= 1.5, ii) located closer to the initial solute source, and iii) occupying a smaller area of the domain. These differences highlight the effect of viscous fingering on i) plume dilution, ii) center of mass location, and iii) extension respectively. The same considerations hold for the higher heterogeneity case (third and fourth rows, right col- umn, in Figure 5.3, as well as earlier times (compare first and last rows in Figure 5.4. In particular, in Figure 5.4 we can see that at earlier times, when the behaviour of σ 2 c is bimodal, the effect of viscosity contrast is to “shift” uncertainty toward the front of the solute plume. This can be seen from the reduction in σ 2 c related to the left peak of the bimodal distribution, and the increase in the peak on the right, i.e. at the front of the plume, for σ 2 f = 0.5 (Figure 5.4, left column). This happens because the effects of viscous fingering compound with the effects of fingering caused by permeability heterogeneity. Therefore, more variability in concentration within the ensemble is found in the leading edge of the plume. This impact of R in shifting the location of uncertainty is hampered by the enhanced dilution caused by permeability heterogeneity forσ 2 f = 1.5 (Figure 5.4, right column), where we can see a reduction in the left peak but no corresponding increase in the right peak. 109 σ 2 f = 0.5 σ 2 f = 1.5 R= 0 R= 1.5 Figure 5.5: Concentration varianceσ 2 c along the x 1 -direction, at t= 5 (continuous lines) and t= 20 (dashed lines), across a horizontal section in the middle of the domain (x 2 = 5λ, dark blue lines), and across a horizontal section edging the bottom boundary of the initial solute source (x 2 = 3λ, light blue lines). 110 From the σ 2 c distribution in the domain reported in Figure 5.3, we focus on the values of σ 2 c along two lines parallel to the mean flow direction: one located in the middle, along the x 2 - direction, of the initial solute source (i.e., x 2 = 5λ), and one located in correspondence of the bottom boundary of the solute source (i.e., x 2 = 3λ). The results are reported in Figure 5.5. For σ 2 f = 1.5, when R= 0 (top right image), viscous fingering is not present to enhance solute mo- bility, and that clearly impacts the back of the plume, where dilution is limited and heterogeneity causes fluctuations in the values of concentration. This persists at later times, as indicated by the asymmetry of the dashed (t= 20) blue line in the top right image in Figure 5.5; such asymmetry is not observed for the lower heterogeneity case (Figure 5.5, top left image), where viscosity contrast is present and enhances the mobility of the rear tailing of the plume. The bimodal trend inσ 2 c is visible for all R andσ 2 f considered at early times and at the center of the domain (continuous, dark blue lines). However, it is not present along the edges of the solute source (x 2 = 3λ, light blue lines), nor at later times (dashed lines). We can also observe, by comparing the top and bottom rows in Figure 5.5 (related to R= 0 and R= 1.5, respectively), that viscosity contrast increases uncertainty at the edges of the plume (light blue lines). This is in contrast with the effect of R at the center of the domain, where the enhanced mixing caused by viscous fingering results in a reduction of uncertainty, as discussed for Figure 5.4. Such increase in uncertainty along x 2 = 3λ when a viscosity contrast is present ceases to exist at later times for σ 2 f = 1.5 (Figure 5.5, left columns, dashed light blue lines), since at later times the impact of permeability heterogeneity on mixing becomes dominant over the impact of viscosity contrast. 5.3.3 Uncertainty in early arrival times Next, we consider another EPM, the early arrival times t e at the control plane located at x 1 = L cp (see Table 5.1). In Figure 5.6, in each image are reported the Probability Density Functions (PDFs) of t e for R= 0 (in blue) and R= 0 (in yellow), for σ 2 f = 0.5 (Figure 5.6, left) and σ 2 f = 1.5 (Figure 5.6, right). When viscous fingering is present, the solute reaches the control plane earlier, a result presented, for instance, in Bonazzi et al. (2021), highlighting the effect of viscosity 111 contrast in enhancing the connectivity of the permeability field. This also causes a reduction in the uncertainty associated with t e , especially for the higher level of permeability heterogeneity considered, as visible from the right image in Figure 5.6, where the tails of the PDF for R= 0 span a larger t e interval compared to the PDF for R= 1.5. We can also see from Figure 5.6 that the Gamma distribution is able to capture the PDFs of the early arrival times at a control plane. The Gamma PDF (Γ-PDF) is given by: G(t e ) = β α Γ(α) t α− 1 e e − βt e , with Γ[α] = Z ∞ 0 ζ α− 1 e − ζ dζ (5.9) whereΓ[α] is the Gamma function andα andβ are the shape and scale parameters, defined as: β = µ t e σ 2 t e , α = µ 2 t e σ 2 t e . (5.10) The Gamma PDF has been used to quantify arrival times PDF (e.g., de Barros, 2018), but to the best of our knowledge this had not been previously shown in the presence of viscous contrast between solute and ambient fluid. We can see from Figure 5.6 that there is a reasonably good agreement between the Γ-PDF and the PDF obtained from the Monte Carlo framework, for all levels of viscosity contrast and permeability heterogeneity considered. However, theΓ-PDF is not perfectly able to capture where the peak of the PDF is located, and the probability associated with it tends to be underestimated. In Figure 5.7 is reported the Empirical Cumulative Distribution Function (ECDF) of the early arrival times t e obtained from the N r Monte Carlo realizations. It is clearly visible that earlier breakthroughs at the control plane are to be expected when both high permeability heterogeneity and viscous fingering are present (red continuous line), as both sources of disorder enhance con- nectivity of the permeability field, thus allowing the front of the plume to travel faster. Lowering the permeability heterogeneity (dashed orange line) results in a lower probability P[t<τ], where 112 σ 2 f = 0.5 σ 2 f = 1.5 Figure 5.6: PDFs of early arrival times t e for the low permeability heterogeneity cases (σ 2 f = 0.5, left) and the high permeability heterogeneity cases (σ 2 f = 1.5, right). The histograms represent the probability distribution obtained from the numerical simulations. The red and purple lines represent (for R= 0 and R= 1.5 respectively) the probability distribution of the gamma distribution model. Figure 5.7: Empirical CDFs of early arrival times t e for the four combinations of R andσ 2 f consid- ered. 113 τ is the early arrival time of which we want to know the probability of non-exceedance. However, maintaining σ 2 f = 1.5 but removing the effect of viscous fingering (compare continuous red and green lines in Figure 5.7) results not only in lower probabilityP[t<τ] (in fact, even lower than for the case withσ 2 f = 0.5 and R= 1.5), but also in longer tailing of the ECDF, especially for higher t e , suggesting that the role of viscous fingering in enhancing connectivity might be predominant over the role of permeability heterogeneity, at least for the values of σ 2 f and R here considered. Later early arrival times are to be expected when neither viscosity contrast nor permeability het- erogeneity contribute to enhancing connectivity (σ 2 f = 0.5 and R= 0, dashed light green line). In this case however, less tailing is present when compared with σ 2 f = 1.5 in absence of viscosity contrast (continuous dark green line). 5.4 Summary We investigated how uncertainty related to the spatial structure of the permeability field affects the concentration field and early arrival times at a target location. We randomly generated 7000 realizations for the permeability field, and performed a Monte Carlo simulation to quantify the statistics of the early arrival times and point concentration values. As expected, higher levels of permeability heterogeneity translates into higher uncertainty in concentration distribution close to the initial solute source at early times. We also showed that the presence of a viscosity contrast between fluids reduces the uncertainty associated with c(x,t) in magnitude, but caused the uncer- tainty to be diffused over a larger portion of the domain. At early times, σ 2 c exhibits a bimodal behaviour along the direction of flow but this trend is not observed at the longitudinal edges of the plume. When this bimodal behaviour is observed, the presence of viscous fingering increases un- certainty at leading edge of the plume, but this effect is hampered by higher levels of permeability heterogeneity. The bimodal trend ofσ 2 c eventually disappears at later times. As expected, viscous fingering enhances connectivity of the permeability field, causing ear- lier breakthroughs at a target location (Bonazzi et al., 2021). We found that the presence of a 114 viscosity contrast reduces the uncertainty associated with early arrival times t e for both levels of permeability heterogeneity considered, but the effect is particularly noticeable for the higher level of heterogeneity. We also show that the Γ-PDF is able to capture the PDF of t e , but it tends to underestimate the peak of probability. Analyzing the CDF of t e , we observed that, for the same level of permeability heterogeneity, viscosity contrast reduces the tailing of the CDF, due to the increased solute mobility caused by viscous fingering. On the other hand, for the same level of viscosity contrast, higher permeability heterogeneity increases tailing of the CDF. Further research is needed to investigate the effect of the uncertainty related to the viscosity model that describes the relationship between concentration and fluid viscosity, to determine its importance relatively to the uncertainty stemming from the lack of characterization of the perme- ability field that was discussed in this Chapter. This would help decision makers to decide where to allocate characterizing efforts and resources in order to reduce uncertainty associated with EPMs of interest. 115 Chapter 6 Concluding remarks This work focused on the effect of multiple sources of disorder on mixing of fluids in porous media. We considered, as sources of disorder, the spatial heterogeneity of the permeability field, the viscosity contrast between the solute and the ambient fluid, and the density contrast. We first focused on the impact of permeability heterogeneity on mixing, and how it is depen- dent on the degree of heterogeneity, on the domain dimensionality, and on the initial geometrical configuration of the solute source. In particular, we showed that for the same level of heterogene- ity, three-dimensional fields exhibits faster decays in time of the solute plume’s spatial mean and variance respect to two-dimensional fields. Faster decays are also associated to a point source com- pared to a line or plane source. We also compared the plume concentration PDF with the beta PDF model, and found that the performance of the beta PDF improves for three-dimensional flow fields. Due to the failure of the beta PDF to capture the extreme values of concentration, we proposed a power low transform of the concentration data, and found that the Pareto-type IV distribution was then able to capture the low tail of the transformed concentration PDF. Then, we studied the interplay between two sources of disorder, permeability heterogeneity and viscosity contrast between the solute and the ambient fluid. We found that, at early times, viscosity contrast is the dominant source of disorder acting on mixing, while at later times heterogeneity becomes the dominant mixing force. We also found that if the solute is less viscous than the am- bient fluid the connectivity of the permeability field at the leading edge of the plume is enhanced, leading to earlier arrival times at a target location. On the other hand, if the solute is more viscous 116 than the ambient fluid, the arrival time is delayed. This opposite impact of viscosity contrast on connectivity is enhanced for higher levels of permeability heterogeneity. We also showed that the plume concentration CDF can be described by the beta distribution model even in the presence of viscosity contrast, for the range of permeability heterogeneity and viscosity contrast considered. We then focused on viscosity contrast and density contrast between a solute and the ambient fluid as sources of disorder on mixing. We analyzed how their relative contribution to solute mixing and spreading is controlled by the initial solute source configuration and by the intensity of the background flux. We showed that if the solute is less viscous than the ambient fluid, the impact of density contrast on mixing is almost negligible for the considered range of density and viscosity contrast. For the higher inlet flux studied, the effect of density fingering on mixing was minor if the solute source was oriented perpendicular to the mean flow direction, while if the solute source was parallelly oriented then density contrast was shown to have an impact on mixing if the solute was more viscous than the ambient fluid. Higher inlet flux hampers the effect of gravitational fingering, rendering viscosity contrast the dominant source of disorder affecting mixing. We also briefly analyzed how a stratified permeability field can interact with fingering and density instabilities and impact the solute transport behaviour. Finally, we investigated how uncertainty stemming from lack of characterization in the struc- ture of the permeability field affects the concentration distribution and early arrival times at a target location. We showed that higher heterogeneity causes higher uncertainty in concentration distri- bution closer to the initial solute source, while viscosity contrast reduces the level of uncertainty in concentration distribution, but increases the overall area of the domain affected by uncertainty. Viscosity contrast also reduces the uncertainty related to early arrival times to a control plane, espe- cially for higher levels of permeability heterogeneity. We also showed that the Gamma distribution model is able to describe the PDF of early arrival times. Future efforts should be devoted to study how mixing evolves in time when all three sources of disorder discussed (permeability heterogeneity, viscosity contrast and density contrast) are acting on solute transport. Studies on this topic should also consider the initial solute source orientation 117 as a factor impacting the concentration distribution in time and space, and modulating the rela- tive importance of the sources of disorder on mixing. Noting that we limited the investigation on the effects of viscous and gravitational instabilities on mixing to two-dimensional domains, further investigation is needed to assess how their interplay (and their interaction with permeabil- ity heterogeneity) translates in three-dimensional domains. Furthermore, future studies should also investigate how the interplay between sources of disorder evolves in time when the solute is slowly released in the domain. Noting that we have considered in this work instantaneous solute injections, our results might be affected if the solute is released in the porous medium through a continuous injection. Another relevant development in this field of study should be to assess the impact of multiple sources of uncertainty arising from the three considered sources of disorder, to better understand where to allocate resources in the effort to characterize the physical properties of both the medium and the mixing fluids. Is uncertainty in a certain quantity of interest (e.g., arrival times at a target location) reduced mainly by reducing uncertainty in the structure of the permeability field? Or will precisely determining the viscosity model, which relates viscosity and concentration, yield the least uncertainty? Or is it the density model that needs to receive more attention? And how is this affected by the initial orientation of the solute source? Much work still needs to be done to fully understand the interplay between these three sources of disorder, their impact on solute mixing, and how they affect uncertainty of quantities of interest in groundwater management applications. This thesis is a further step down the road to a more comprehensive knowledge on fluids mixing, and we hope it might be a source of inspiration for future scientific endeavours. 118 Chapter 7 List of publications and papers under review • Bonazzi, A., Jha, B., & de Barros, F.P.J. (2021). Transport analysis in deformable porous media through integral transforms. International Journal for Numerical and Analytical Methods in Geomechanics, 45(3), 307-324. • Morvillo, M., Bonazzi, A., Rizzo, C. B., & de Barros, F.P.J. (2021). Improving the com- putational efficiency of first arrival time uncertainty estimation using a connectivity-based ranking Monte Carlo method. Stochastic Environmental Research and Risk Assessment, 35, 1039-1049. • Bonazzi, A., Morvillo, M., Im, J., Jha, B., & de Barros, F.P.J. (2021). Relative impacts of permeability heterogeneity and viscosity contrast on solute mixing. Physical Review Fluids, 6(6), 064501. • Bonazzi, A., Dentz, M., & de Barros, F.P.J. (2023). Mixing in multidimensional porous media: a numerical study of the effects of source configuration and heterogeneity. Transport in Porous Media, 146(1-2), 369-393. • Bonazzi, A., Jha, B., & de Barros, F.P.J. Influence of initial plume shape on miscible porous media flows under density and viscosity contrasts. Journal of Fluid Mechanics (under re- view). 119 Bibliography ABABOU, R., MCLAUGHLIN, D., GELHAR, L.W. & TOMPSON, A.F.B. 1989a Numerical simu- lation of three-dimensional saturated flow in randomly heterogeneous porous media. Transport in Porous Media 4 (6), 549–565. ABABOU, R., MCLAUGHLIN, D., GELHAR, L. W. & TOMPSON, A. F. B. 1989b Numerical sim- ulation of three-dimensional saturated flow in randomly heterogeneous porous media. Transport in Porous Media 4 (6), 549–565. ANDRI ˇ CEVI ´ C, R. 1998 Effects of local dispersion and sampling volume on the evolution of con- centration fluctuations in aquifers. Water Resources Research 34 (5), 1115–1129. BACRI, J.C., RAKOTOMALALA, N., SALIN, D. & WOUMENI, R. 1992 Miscible viscous fin- gering: Experiments versus continuum approach. Physics of Fluids A: Fluid Dynamics 4 (8), 1611–1619. BAGTZOGLOU, A.C. & OATES, P.M. 2007 Chaotic advection and enhanced groundwater reme- diation. Journal of Materials in Civil Engineering 19 (1), 75–83. BAKKER, M., POST, V., LANGEVIN, C. D., HUGHES, J. D., WHITE, J. T., STARN, J.J. & FIENEN, M. N 2016 Scripting modflow model development using python and flopy. Ground- water 54 (5), 733–739. DE BARROS, F.P.J. 2018 Evaluating the combined effects of source zone mass release rates and aquifer heterogeneity on solute discharge uncertainty. Advances in Water Resources 117, 140– 150. DE BARROS, F.P.J., BELLIN, A., CVETKOVIC, V., DAGAN, G. & FIORI, A. 2016 Aquifer heterogeneity controls on adverse human health effects and the concept of the hazard attenuation factor. Water Resources Research 52 (8), 5911–5922. DE BARROS, F.P.J., DENTZ, M., KOCH, J. & NOWAK, W. 2012a Flow topology and scalar mixing in spatially heterogeneous flow fields. Geophysical Research Letters 39. DE BARROS, F.P.J., EZZEDINE, S. & RUBIN, Y. 2012b Impact of hydrogeological data on mea- sures of uncertainty, site characterization and environmental performance metrics. Advances in Water Resources 36, 51–63. DE BARROS, F.P.J., FERN ` ANDEZ-GARCIA, D., BOLSTER, D. & SANCHEZ-VILA, X. 2013 A risk-based probabilistic framework to estimate the endpoint of remediation: Concentration rebound by rate-limited mass transfer. Water Resources Research 49 (4), 1929–1942. 120 DE BARROS, F.P.J. & FIORI, A. 2014 First-order based cumulative distribution function for solute concentration in heterogeneous aquifers: Theoretical analysis and implications for human health risk assessment. Water Resources Research 50 (5), 4018–4037. DE BARROS, F.P.J. & FIORI, A. 2021 On the maximum concentration of contaminants in natural aquifers. Transport in Porous Media pp. 1–18. DE BARROS, F.P.J., FIORI, A., BOSO, F. & BELLIN, A. 2015 A theoretical framework for modeling dilution enhancement of non-reactive solutes in heterogeneous porous media. Journal of Contaminant Hydrology 175, 72–83. DE BARROS, F.P.J., FIORI, A., BOSO, F. & BELLIN, A. 2015 A theoretical framework for modeling dilution enhancement of non-reactive solutes in heterogeneous porous media. Journal of Contaminant Hydrology 175, 72–83. DE BARROS, F.P.J, GUADAGNINI, A. & RIVA, M. 2022 Features of transport in non-gaussian random porous systems. International Journal of Heat and Mass Transfer 184, 122244. DE BARROS, F.P.J. & NOWAK, W. 2010 On the link between contaminant source release condi- tions and plume prediction uncertainty. Journal of Contaminant Hydrology 116 (1-4), 24–34. DE BARROS, F.P.J. & RUBIN, Y. 2011 Modelling of block-scale macrodispersion as a random function. Journal of Fluid Mechanics 676, 514–545. DE BARROS, F.P.J., SOUHEIL, E. & RUBIN, Y. 2012c Impact of hydrogeological data on mea- sures of uncertainty, site characterization and environmental performance metrics. Advances in Water Resources 36, 51–63. BEAUDOIN, A. & DE DREUZY, J.R. 2013 Numerical assessment of 3-d macrodispersion in het- erogeneous porous media. Water Resources Research 49 (5), 2489–2496. BELLIN, A., RINALDO, A., BOSMA, W. J. P., VAN DER ZEE, S. E. A. T. M. & RUBIN, Y. 1993 Linear equilibrium adsorbing solute transport in physically and chemically heterogeneous porous formations: 1. analytical solutions. Water Resources Research 29 (12), 4019–4030. BELLIN, A. & RUBIN, Y. 1996 Hydro gen: A spatially distributed random field generator for correlated properties. Stochastic Hydrology and Hydraulics 10 (4), 253–278. BELLIN, A., RUBIN, Y. & RINALDO, A. 1994 Eulerian-lagrangian approach for modeling of flow and transport in heterogeneous geological formations. Water Resources Research 30 (11), 2913–2924. BELLIN, A., SALANDIN, P. & RINALDO, A. 1992a Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. Water Re- sources Research 28 (9), 2211–2227. BELLIN, A., SALANDIN, P. & RINALDO, A. 1992b Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. Water Re- sources Research 28 (9), 2211–2227. 121 BELLIN, A. & TONINA, D. 2007 Probability density function of non-reactive solute concentration in heterogeneous porous formations. Journal of Contaminant Hydrology 94 (1-2), 109–125. BONAZZI, A., DENTZ, M. & DE BARROS, F.P.J. 2022 Mixing in multidimensional porous media: A numerical study of the effects of source configuration and heterogeneity. Transport in Porous Media pp. 1–25. BONAZZI, A., MORVILLO, M., IM, J., JHA, B. & DE BARROS, F.P.J. 2021 Relative impacts of permeability heterogeneity and viscosity contrast on solute mixing. Physical Review Fluids 6 (6), 064501. BOSO, F., DE BARROS, F.P.J., FIORI, A. & BELLIN, A. 2013a Performance analysis of statistical spatial measures for contaminant plume characterization toward risk-based decision making. Water Resources Research 49 (6), 3119–3132. BOSO, F., BELLIN, A. & DUMBSER, M. 2013b Numerical simulations of solute transport in highly heterogeneous formations: A comparison of alternative numerical schemes. Advances in Water Resources 52, 178–189. BOSO, F. & TARTAKOVSKY, D.M. 2016 The method of distributions for dispersive transport in porous media with uncertain hydraulic properties. Water Resources Research 52 (6), 4700–4712. BROYLES, B.S., SHALLIKER, R.A., CHERRAK, D.E. & GUIOCHON, G. 1998 Visualization of viscous fingering in chromatographic columns. Journal of Chromatography A 822 (2), 173–187. CAMHI, EMMANUEL, MEIBURG, ECKART & RUITH, MICHAEL 2000 Miscible rectilinear dis- placements with gravity override. Part 2. Heterogeneous porous media. Journal of Fluid Me- chanics 420, 259–276. CARONI, E. & FIOROTTO, V. 2005 Analysis of concentration as sampled in natural aquifers. Transport in porous media 59, 19–45. CHAPMAN, S. W. & PARKER, B. L. 2005 Plume persistence due to aquitard back diffusion fol- lowing dense nonaqueous phase liquid source removal or isolation. Water Resources Research 41 (12). CHEN, C.Y. & MEIBURG, E. 1998 Miscible porous media displacements in the quarter five-spot configuration. part 2. effect of heterogeneities. Journal of Fluid Mechanics 371, 269–299. CHRISTIE, M.A. 1989 High-resolution simulation of unstable flows in porous media. SPE Reser- voir Engineering 4 (03), 297–303. CIRPKA, O.A., DE BARROS, F.P.J., CHIOGNA, G., ROLLE, M. & NOWAK, W. 2011 Stochas- tic flux-related analysis of transverse mixing in two-dimensional heterogeneous porous media. Water Resources Research 47 (6). CIRPKA, O.A., ROLLE, M., CHIOGNA, G., DE BARROS, F.P.J. & NOWAK, W. 2012 Stochastic evaluation of mixing-controlled steady-state plume lengths in two-dimensional heterogeneous domains. Journal of Contaminant Hydrology 138, 22–39. 122 COWELL, S., KENT, J. & TREVELYAN, P.M.J. 2020 Rayleigh–taylor instabilities in miscible flu- ids with initially piecewise linear density profiles. Journal of Engineering Mathematics 121 (1), 57–83. CVETKOVIC, V., FIORI, A. & DAGAN, G. 2016 Tracer travel and residence time distributions in highly heterogeneous aquifers: Coupled effect of flow variability and mass transfer. Journal of Hydrology 543, 101–108. CVETKOVIC, V, SHAPIRO, AM & DAGAN, G 1992 A solute flux approach to transport in hetero- geneous formations: 2. uncertainty analysis. Water Resources Research 28 (5), 1377–1388. DAGAN, G. 1982 Stochastic modeling of groundwater flow by unconditional and conditional prob- abilities: 2. the solute transport. Water Resources Research 18 (4), 835–848. DAGAN, G. 1984 Solute transport in heterogeneous porous formations. Journal of Fluid Mechan- ics 145, 151–177. DAGAN, G. & NEUMAN, S.P. 2005 Subsurface flow and transport: a stochastic approach . Cam- bridge University Press. DAI, Z., ZHANG, Y., BIELICKI, J., AMOOIE, M.A., ZHANG, M., YANG, C., ZOU, Y., AM- POMAH, W., XIAO, T., JIA, W. & OTHERS 2018 Heterogeneity-assisted carbon dioxide storage in marine sediments. Applied Energy 225, 876–883. DANIEL, D. & RIAZ, A. 2014 Effect of viscosity contrast on gravitationally unstable diffusive layers in porous media. Physics of Fluids 26. DAVIES, G.W., MUGGERIDGE, A.H. & JONES, A.D.W. 1991 Miscible displacements in a het- erogeneous rock: detailed measurements and accurate predictive simulation. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. DE WIT, A., BERTHO, Y. & MARTIN, M. 2005 Viscous fingering of miscible slices. Physics of Fluids 17 (5), 054114. DE WIT, A. & HOMSY, G.M. 1997a Viscous fingering in periodically heterogeneous porous media. i. formulation and linear instability. The Journal of Chemical Physics 107 (22), 9609– 9618. DE WIT, A. & HOMSY, G.M. 1997b Viscous fingering in periodically heterogeneous porous media. ii. numerical simulations. The Journal of Chemical Physics 107 (22), 9619–9628. DENTZ, M. 2012 Concentration statistics for transport in heterogeneous media due to stochastic fluctuations of the center of mass velocity. Advances in Water Resources 36, 11–22. DENTZ, M. & DE BARROS, F.P.J. 2013a Dispersion variance for transport in heterogeneous porous media. Water Resources Research 49 (6), 3443–3461. DENTZ, M. & DE BARROS, F.P.J. 2013b Dispersion variance for transport in heterogeneous porous media. Water Resources Research 49 (6), 3443–3461. 123 DENTZ, M. & DE BARROS, F.P.J. 2015 Mixing-scale dependent dispersion for transport in het- erogeneous flows. Journal of Fluid Mechanics 777, 178–195. DENTZ, M., DE BARROS, F.P.J., LE BORGNE, T. & LESTER, D.R. 2018 Evolution of solute blobs in heterogeneous porous media. Journal of Fluid Mechanics 853, 621–646. DENTZ, M. & CARRERA, J. 2007 Mixing and spreading in stratified flow. Physics of Fluids 19 (1), 017107. DENTZ, M., KINZELBACH, H., ATTINGER, S. & KINZELBACH, W. 2000a Temporal behavior of a solute cloud in a heterogeneous porous medium: 1. point-like injection. Water Resources Research 36 (12), 3591–3604. DENTZ, M., KINZELBACH, H., ATTINGER, S. & KINZELBACH, W. 2000b Temporal behavior of a solute cloud in a heterogeneous porous medium: 2. Spatially extended injection. Water Resources Research 36 (12), 3605–3614. DENTZ, M., LE BORGNE, T., ENGLERT, A. & BIJELJIC, B. 2011 Mixing, spreading and reaction in heterogeneous media: A brief review. Journal of Contaminant Hydrology 120, 1–17. DENTZ, M., TARTAKOVSKY, D.M., ABARCA, E., GUADAGNINI, A., S ´ ANCHEZ-VILA, X. & CARRERA, J. 2006 Variable-density flow in porous media. Journal of Fluid Mechanics 561, 209–235. DENTZ, M. & TARTAKOVSKY, D. M. 2010 Probability density functions for passive scalars dis- persed in random velocity fields. Geophysical Research Letters 37 (24). DI DATO, M., DE BARROS, F.P.J., FIORI, A. & BELLIN, A. 2018 Improving the efficiency of 3-d hydrogeological mixers: Dilution enhancement via coupled engineering-induced transient flows and spatial heterogeneity. Water Resources Research 54 (3), 2095–2111. DIETER, C.A., MAUPIN, M.A., CALDWELL, R.R., HARRIS, M.A., IVAHNENKO, T.I., LOVELACE, J.K., BARBER, N.L. & LINSEY, K.S. 2018 Estimated use of water in the United States in 2015. U.S. Geological Survey. DE DREUZY, J-R, BEAUDOIN, A. & ERHEL, J. 2007 Asymptotic dispersion in 2d heterogeneous porous media determined by parallel numerical simulations. Water Resources Research 43 (10). ENZENHOEFER, R., NOWAK, W. & HELMIG, R. 2012 Probabilistic exposure risk assessment with advective–dispersive well vulnerability criteria. Advances in water resources 36, 121–132. FERZIGER, J.H., PERI ´ C, M. & STREET, R.L. 2002a Computational methods for fluid dynamics , , vol. 3. Springer. FERZIGER, J.H., PERI ´ C, M. & STREET, R.L. 2002b Computational methods for fluid dynamics . Springer. FIORI, A. 2001 The lagrangian concentration approach for determining dilution in aquifer trans- port: Theoretical analysis and comparison with field experiments. Water Resources Research 37 (12), 3105–3114. 124 FIORI, A., BELLIN, A., CVETKOVIC, V., DE BARROS, F.P.J. & DAGAN, G. 2015 Stochastic modeling of solute transport in aquifers: From heterogeneity characterization to risk analysis. Water Resources Research 51 (8), 6622–6648. FIORI, A. & DAGAN, G. 2000 Concentration fluctuations in aquifer transport: A rigorous first- order solution and applications. Journal of Contaminant Hydrology 45 (1-2), 139–163. FIORI, A., ZARLENGA, A., BELLIN, A., CVETKOVIC, V. & DAGAN, G. 2019 Groundwater contaminant transport: Prediction under uncertainty, with application to the made transport ex- periment. Frontiers in Environmental Science 7, 79. FIORI, A., ZARLENGA, A., JANKOVIC, I. & DAGAN, G. 2017 Solute transport in aquifers: The comeback of the advection dispersion equation and the first order approximation. Advances in Water Resources 110, 349–359. FIOROTTO, V. & CARONI, E. 2002 Solute concentration statistics in heterogeneous aquifers for finite peclet values. Transport in Porous Media 48 (3), 331–351. FLOWERS, T.C. & HUNT, J.R. 2007 Viscous and gravitational contributions to mixing during vertical brine transport in water-saturated porous media. Water resources research 43 (1). FRAILEY, S.M. & LEETARU, H. 2009 Geological factors affecting CO 2 plume distribution. En- ergy Procedia 1 (1), 3107–3112. GHESMAT, K. & AZAIEZ, J. 2008 Viscous fingering instability in porous media: effect of anisotropic velocity-dependent dispersion tensor. Transport in Porous Media 73 (3), 297–318. GOPALAKRISHNAN, S.S. 2020 On the instability of buoyancy-driven flows in porous media. Jour- nal of Fluid Mechanics 892. GOTOVAC, H., ANDRICEVIC, R. & GOTOVAC, B. 2007 Multi-resolution adaptive modeling of groundwater flow and transport problems. Advances in Water Resources 30 (5), 1105–1126. GUETING, N. & ENGLERT, A. 2013 Hydraulic conditions at the source zone and their impact on plume behavior. Hydrogeology Journal 21 (4), 829–844. HARBAUGH, A. W. 2005 MODFLOW-2005, the US Geological Survey modular ground-water model: the ground-water flow process . US Department of the Interior, US Geological Survey Reston, V A. HARTMANN, A., GLEESON, T., WADA, Y. & WAGENER, T. 2017 Enhanced groundwater recharge rates and altered recharge sensitivity to climate variability through subsurface hetero- geneity. Proceedings of the National Academy of Sciences 114 (11), 2842–2847. HASSANZADEH, H., POOLADI-DARVISH, M. & KEITH, D.W. 2007 Scaling behavior of convec- tive mixing, with application to geological storage of CO 2 . AIChE Journal 53 (5), 1121–1131. HELD, R.J. & ILLANGASEKARE, T.H. 1995a Fingering of dense nonaqueous phase liquids in porous media: 1. experimental investigation. Water Resources Research 31 (5), 1213–1222. 125 HELD, R.J. & ILLANGASEKARE, T.H. 1995b Fingering of dense nonaqueous phase liquids in porous media: 1. experimental investigation. Water Resources Research 31 (5), 1213–1222. HENRI, C.V.AND FERN ` ANDEZ-GARCIA, D. & DE BARROS, F.P.J. 2016 Assessing the joint impact of dnapl source-zone behavior and degradation products on the probabilistic characteri- zation of human health risk. Advances in Water Resources 88, 124–138. HERRERA, P. A., MASSAB ´ O, M. & BECKIE, R. D. 2009 A meshless method to simulate solute transport in heterogeneous porous media. Advances in Water Resources 32 (3), 413–429. HEWITT, D.R., NEUFELD, J.A. & LISTER, J.R. 2013 Convective shutdown in a porous medium at high Rayleigh number. Journal of Fluid Mechanics 719, 551–586. HIDALGO, J.J. & CARRERA, J. 2009 Effect of dispersion on the onset of convection during CO 2 sequestration. Journal of Fluid Mechanics 640, 441–452. HIDALGO, J.J. & DENTZ, M. 2018a Mixing across fluid interfaces compressed by convective flow in porous media. Journal of Fluid Mechanics 838, 105–128. HIDALGO, J.J. & DENTZ, M. 2018b Mixing across fluid interfaces compressed by convective flow in porous media. Journal of Fluid Mechanics 838, 105–128. HIDALGO, J.J., MACMINN, C.W. & JUANES, R. 2013 Dynamics of convective dissolution from a migrating current of carbon dioxide. Advances in Water Resources 62, 511–519. HOMSY, G.M. 1987 Viscous fingering in porous media. Annual Review of Fluid Mechanics 19 (1), 271–311. HORTON, C.W. & ROGERS JR, F.T. 1945 Convection currents in a porous medium. Journal of Applied Physics 16 (6), 367–370. HOVORKA, S.D., DOUGHTY, C., BENSON, S.M., PRUESS, K. & KNOX, P.R. 2004 The impact of geological heterogeneity on CO 2 storage in brine formations: a case study from the Texas Gulf Coast. Geological Society, London, Special Publications 233 (1), 147–163. HUPPERT, H.E. & NEUFELD, J.A. 2014 The fluid mechanics of carbon dioxide sequestration. Annual review of fluid mechanics 46, 255–272. IM, J., RIZZO, C.B. & DE BARROS, F.P.J. 2020 Resilience of groundwater systems in the pres- ence of bisphenol a under uncertainty. Science of The Total Environment 727, 138363. JANKOVIC, I., MAGHREBI, M., FIORI, A. & DAGAN, G. 2017 When good statistical models of aquifer heterogeneity go right: The impact of aquifer permeability structures on 3d flow and transport. Advances in Water Resources 100, 199–211. JENNY, P., LEE, J.S., MEYER, D.W. & TCHELEPI, H.A. 2014 Scale analysis of miscible density- driven convection in porous media. Journal of Fluid Mechanics 749, 519–541. JHA, B. 2014a Flow through porous media: From mixing of fluids to triggering of earthquakes. PhD thesis, Massachusetts Institute of Technology, Cambridge, MA. 126 JHA, B. 2014b Flow through porous media: From mixing of fluids to triggering of earthquakes. PhD thesis, Massachusetts Institute of Technology. JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2011a Fluid mixing from viscous fingering. Physical Review Letters 106. JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2011b Quantifying mixing in viscously unsta- ble porous media flows. Physical Review E 84 (6), 066312. JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2013a Synergetic fluid mixing from viscous fingering and alternating injection. Physical Review Letters 111 (14), 144501. JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2013b Synergetic fluid mixing from viscous fingering and alternating injection. Physical Review Letters 111 (14), 144501. KAMRAVA, S., IM, J., DE BARROS, F.P.J. & SAHIMI, M. 2021 Estimating dispersion coefficient in flow through heterogeneous porous media by a deep convolutional neural network. Geophys- ical Research Letters 48 (18), e2021GL094443. KAPOOR, V. & GELHAR, L.W. 1994 Transport in three-dimensionally heterogeneous aquifers: 1. Dynamics of concentration fluctuations. Water Resources Research 30 (6), 1775–1788. KAPOOR, V. & KITANIDIS, P.K. 1998 Concentration fluctuations and dilution in aquifers. Water Resources Research 34 (5), 1181–1193. KEMPERS, L.J.T.M. & HAAS, H. 1994a The dispersion zone between fluids with different den- sity and viscosity in a heterogeneous porous medium. Journal of Fluid Mechanics 267, 299–324. KEMPERS, L.J.T.M. & HAAS, H. 1994b The dispersion zone between fluids with different den- sity and viscosity in a heterogeneous porous medium. Journal of Fluid Mechanics 267, 299–324. KITANIDIS, P.K. 1997 Introduction to geostatistics: applications in hydrogeology. Cambridge university press. KITANIDIS, P. K. 1994 The concept of the dilution index. Water Resources Research 30 (7), 2011– 2026. KITANIDIS, P. K. & MCCARTY, P. L. 2012 Delivery and Mixing in the Subsurface: Processes and Design Principles for in Situ remediation. Springer Science & Business Media. KOCH, D.L. & BRADY, J.F. 1988 Anomalous diffusion in heterogeneous porous media. Physics of Fluids 31 (5), 965–973. KOCH, J. & NOWAK, W. 2015 Predicting dnapl mass discharge and contaminated site longevity probabilities: Conceptual model and high-resolution stochastic simulation. Water Resources Re- search 51 (2), 806–831. KOPF-SILL, A.R. & HOMSY, G.M. 1988 Nonlinear unstable viscous fingers in Hele-Shaw flows. i. experiments. Physics of Fluids 31 (2), 242–249. 127 KRETZ, V., BEREST, P., HULIN, J.P. & SALIN, D. 2003 An experimental study of the effects of density and viscosity contrasts on macrodispersion in porous media. Water Resources Research 39 (2). LAKE, L.W. 1989 Enhanced oil recovery. Old Tappan, NJ; Prentice Hall Inc. LE BORGNE, T., DENTZ, D., BOLSTER, D., CARRERA, J., DE DREUZY, J.R. & DAVY, P. 2010 Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media. Advances in Water Resources 33, 1468–1475. LE BORGNE, T., DENTZ, M. & VILLERMAUX, E. 2013 Stretching, coalescence, and mixing in porous media. Physical Review Letters 110 (20), 204501. LE BORGNE, T., DENTZ, M. & VILLERMAUX, E. 2015a The lamellar description of mixing in porous media. Journal of Fluid Mechanics 770, 458–498. LE BORGNE, T., DENTZ, M. & VILLERMAUX, E. 2015b The lamellar description of mixing in porous media. Journal of Fluid Mechanics 770, 458–498. LEUBE, P.C., DE BARROS, F.P.J., NOWAK, W. & RAJAGOPAL, R. 2013 Towards optimal allo- cation of computer resources: Trade-offs between uncertainty quantification, discretization and model reduction. Environmental Modelling & Software 50, 97–107. LEVEQUE, R. J. 2002 Finite Volume Methods for Hyperbolic Problems. UK: Cambridge Univer- sity Press. LIBERA, A., DE BARROS, F.P.J. & GUADAGNINI, A. 2017 Influence of pumping operational schedule on solute concentrations at a well in randomly heterogeneous aquifers. Journal of hy- drology 546, 490–502. LIBERA, A., HENRI, C.V. & DE BARROS, F.P.J. 2019 Hydraulic conductivity and porosity het- erogeneity controls on environmental performance metrics: Implications in probabilistic risk analysis. Advances in Water Resources 127, 1–12. LINDE, N., GINSBOURGER, D., IRVING, J., NOBILE, F. & DOUCET, A. 2017 On uncertainty quantification in hydrogeology and hydrogeophysics. Advances in Water Resources 110, 166– 181. LIYANAGE, R., RUSSELL, A., CRAWSHAW, J.P. & KREVOR, S. 2020 Direct experimental obser- vations of the impact of viscosity contrast on convective mixing in a three-dimensional porous medium. Physics of Fluids 32. LOGGIA, D., RAKOTOMALALA, N., SALIN, D. & YORTSOS, Y.C. 1996 Phase diagram of stable miscible displacements in layered porous media. EPL (Europhysics Letters) 36 (2), 105. MAXWELL, R.M. & KASTENBERG, W.E. 1999 Stochastic environmental risk analysis: an inte- grated methodology for predicting cancer risk from contaminated groundwater. Stochastic En- vironmental Research and Risk Assessment 13 (1-2), 27–47. 128 VAN DER MEER, L.G.H. 1993 The conditions limiting CO 2 storage in aquifers. Energy Conver- sion and Management 34 (9-11), 959–966. MEYER, D. W., JENNY, P. & TCHELEPI, H. A. 2010 A joint velocity-concentration pdf method for tracer flow in heterogeneous porous media. Water Resources Research 46 (12). MISHRA, M., MARTIN, M. & DE WIT, A. 2008 Differences in miscible viscous fingering of finite width slices with positive or negative log-mobility ratio. Physical Review E 78 (6), 066306. MORVILLO, M., BONAZZI, A., RIZZO, C.B. & DE BARROS, F.P.J. 2021a Improving the compu- tational efficiency of first arrival time uncertainty estimation using a connectivity-based ranking monte carlo method. Stochastic Environmental Research and Risk Assessment 35 (5), 1039– 1049. MORVILLO, M., IM, J. & DE BARROS, F.P.J. 2022 Visu-hydra: A computational toolbox for groundwater contaminant transport to support risk-based decision making. Frontiers in Earth Science 10. MORVILLO, M., RIZZO, C.B. & DE BARROS, F.P.J. 2021b A scalable parallel algorithm for reactive particle tracking. Journal of Computational Physics 446, 110664. MOSLEHI, M. & DE BARROS, F.P.J. 2017 Uncertainty quantification of environmental perfor- mance metrics in heterogeneous aquifers with long-range correlations. Journal of contaminant hydrology 196, 21–29. MOSLEHI, M., RAJAGOPAL, R. & DE BARROS, F.P.J. 2015 Optimal allocation of computational resources in hydrogeological models under uncertainty. Advances in Water Resources 83, 299– 309. MUNRO, R.J., CHATWIN, P.C. & MOLE, N. 2003 A concentration pdf for the relative dispersion of a contaminant plume in the atmosphere. Boundary-Layer Meteorology 106 (3), 411–436. NEUMAN, S. P. & TARTAKOVSKY, D. M. 2009 Perspective on theories of non-fickian transport in heterogeneous media. Advances in Water Resources 32 (5), 670–680. NICOLAIDES, C., JHA, B., CUETO-FELGUEROSO, L. & JUANES, R. 2015 Impact of viscous fingering and permeability heterogeneity on fluid mixing in porous media. Water Resources Research 51 (4), 2634–2647. NIJJER, J.S., HEWITT, D.R. & NEUFELD, J.A. 2018 The dynamics of miscible viscous fingering from onset to shutdown. Journal of Fluid Mechanics 837, 520–545. NIJJER, J.S., HEWITT, D.R. & NEUFELD, J.A. 2019 Stable and unstable miscible displacements in layered porous media. Journal of Fluid Mechanics 869, 468–499. NORDBOTTEN, J.M., CELIA, M.A. & BACHU, S. 2005 Injection and storage of CO 2 in deep saline aquifers: Analytical solution for CO 2 plume evolution during injection. Transport in Porous Media 58 (3), 339–360. 129 ORR, JR F.M. & TABER, J.J. 1984 Use of carbon dioxide in enhanced oil recovery. Science 224 (4649), 563–569. OTTINO, J.M. 1989 The Kinematics of Mixing: Stretching, Chaos, and Transport. Cambridge University Press. PETITJEANS, P., CHEN, C.Y., MEIBURG, E. & MAXWORTHY, T. 1999 Miscible quarter five- spot displacements in a hele-shaw cell and the role of flow-induced dispersion. Physics of Fluids 11 (7), 1705–1716. PETITJEANS, P. & MAXWORTHY, T. 1996 Miscible displacements in capillary tubes. part 1. ex- periments. Journal of Fluid Mechanics 326, 37–56. PISCOPO, A.N., NEUPAUER, R.M. & MAYS, D.C. 2013 Engineered injection and extraction to enhance reaction for improved in situ remediation. Water Resources Research 49 (6), 3618– 3625. PRAMANIK, S. & MISHRA, M. 2016 Coupled effect of viscosity and density gradients on finger- ing instabilities of a miscible slice in porous media. Physics of Fluids 28. REMY, N., BOUCHER, A. & WU, J. 2009 Applied geostatistics with SGeMS: A user’s guide. Cambridge University Press. RIAZ, A., HESSE, M., TCHELEPI, H.A. & ORR, F.M. 2006 Onset of convection in a gravita- tionally unstable diffusive boundary layer in porous media. Journal of Fluid Mechanics 548, 87–111. RIZZO, C.B. & DE BARROS, F.P.J. 2017 Minimum hydraulic resistance and least resistance path in heterogeneous porous media. Water Resources Research 53 (10), 8596–8613. RIZZO, C.B. & DE BARROS, F.P.J. 2019 Minimum hydraulic resistance uncertainty and the de- velopment of a connectivity-based iterative sampling strategy. Water Resources Research 55 (7), 5593–5611. RIZZO, C.B., NAKANO, A. & DE BARROS, F.P.J. 2019 Par2: parallel random walk particle tracking method for solute transport in porous media. Computer Physics Communications 239, 265–271. ROLLE, M., EBERHARDT, C., CHIOGNA, G., CIRPKA, O.A. & GRATHWOHL, P. 2009 Enhance- ment of dilution and transverse reactive mixing in porous media: Experiments and model-based interpretation. Journal of Contaminant Hydrology 110 (3-4), 130–142. RUBIN, B., BARKER, J.W., BLUNT, M.J., CHRISTIE, M.A., CULVERWELL, I.D. & MANS- FIELD, M. 1993 Compositional reservoir simulation with a predictive model for viscous finger- ing. In SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers. RUBIN, Y. 1991 The spatial and temporal moments of tracer concentration in disordered porous media. Water Resources Research 27 (11), 2845–2854. RUBIN, Y. 2003 Applied stochastic hydrogeology. Oxford University Press. 130 RUBIN, Y., CUSHEY, M.A. & BELLIN, A. 1994 Modeling of transport in groundwater for envi- ronmental risk assessment. Stochastic Hydrology and Hydraulics 8 (1), 57–77. RUBIN, Y. & DAGAN, G. 1992 Conditional estimation of solute travel time in heterogeneous formations: Impact of transmissivity measurements. Water Resources Research 28 (4), 1033– 1040. RUBIN, Y., SUN, A., MAXWELL, R. & BELLIN, A 1999 The concept of block-effective macrodispersivity and a unified approach for grid-scale-and plume-scale-dependent transport. Journal of Fluid Mechanics 395, 161–180. SAFFMAN, P.G. 1986 Viscous fingering in Hele-Shaw cells. Journal of Fluid Mechanics 173, 73–94. SAFFMAN, P.G. & TAYLOR, G.I. 1958 The penetration of a fluid into a porous medium or hele- shaw cell containing a more viscous liquid. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 245 (1242), 312–329. SAHIMI, M. 2011 Flow and transport in porous media and fractured rock: from classical methods to modern approaches. John Wiley & Sons. SAHIMI, M., DAVIS, H.T. & SCRIVEN, L.E. 1983 Dispersion in disordered porous media. Chem- ical Engineering Communications 23 (4-6), 329–341. SAJJADI, M. & AZAIEZ, J. 2013 Scaling and unified characterization of flow instabilities in lay- ered heterogeneous porous media. Physical Review E 88 (3), 033017. SALAMON, P., FERN ` ANDEZ-GARCIA, D. & G ´ OMEZ-HERN ´ ANDEZ, J.J. 2006 A review and nu- merical assessment of the random walk particle tracking method. Journal of Contaminant Hy- drology 87 (3-4), 277–305. SANCHEZ-VILA, X. & GUADAGNINI, A. 2005 Travel time and trajectory moments of conser- vative solutes in three dimensional heterogeneous porous media under mean uniform flow. Ad- vances in water resources 28 (5), 429–439. SCHEIDEGGER, A.E. 1961 General theory of dispersion in porous media. Journal of Geophysical Research 66 (10), 3273–3278. SCHWEDE, R.L., CIRPKA, O.A., NOWAK, W. & NEUWEILER, I. 2008 Impact of sampling vol- ume on the probability density function of steady state concentration. Water Resources Research 44 (12). SELROOS, J.O. 1997 A stochastic-analytical framework for safety assessment of waste reposito- ries: 1. theory. Groundwater 35 (3), 468–477. SLIM, A.C. 2014 Solutal-convection regimes in a two-dimensional porous medium. Journal of Fluid Mechanics 741, 461–491. SLIM, A.C., BANDI, M.M., MILLER, J.C. & MAHADEVAN, L. 2013 Dissolution-driven con- vection in a hele–shaw cell. Physics of Fluids 25 (2), 024101. 131 SOLE-MARI, G., RIVA, M., FERN ` ANDEZ-GARCIA, D., SANCHEZ-VILA, X. & GUADAGNINI, A. 2021 Solute transport in bounded porous media characterized by generalized sub-gaussian log-conductivity distributions. Advances in Water Resources 147, 103812. SOLTANIAN, M.R., BEHZADI, F. & DE BARROS, F.P.J. 2020 Dilution enhancement in hierarchi- cal and multiscale heterogeneous sediments. Journal of Hydrology 587, 125025. SZULCZEWSKI, M.L., HESSE, M.A. & JUANES, R. 2013 Carbon dioxide dissolution in structural and stratigraphic traps. Journal of Fluid Mechanics 736. TALON, L., MARTIN, J., RAKOTOMALALA, N. & SALIN, D. 2004 Stabilizing viscosity contrast effect on miscible displacement in heterogeneous porous media, using lattice bhatnagar–gross– krook simulations. Physics of Fluids 16 (12), 4408–4411. TAN, C.T. & HOMSY, G.M. 1988 Simulation of nonlinear viscous fingering in miscible displace- ment. Physics of Fluids 31 (6), 1330–1338. TAN, C.T. & HOMSY, G.M. 1992 Viscous fingering with permeability heterogeneity. Physics of Fluids A: Fluid Dynamics 4 (6), 1099–1101. TAYLOR, G.I. 1953 Dispersion of soluble matter in solvent flowing slowly through a tube. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 219 (1137), 186–203. TCHELEPI, H.A. & ORR, F.M. 1994 Interaction of viscous fingering, permeability heterogeneity, and gravity segregation in three dimensions. SPE Reservoir Engineering 9 (04), 266–271. TCHELEPI, H.A. & ORR JR, F.M. 1993 The interaction of viscous fingering, permeability het- erogeneity and gravity segregation in 3d. In SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers. TCHELEPI, H.A. & ORR JR, F.M. 1994 Interaction of viscous fingering, permeability heterogene- ity, and gravity segregation in three dimensions. SPE Reservoir Engineering 9 (04), 266–271. TONINA, D. & BELLIN, A. 2008 Effects of pore-scale dispersion, degree of heterogeneity, sam- pling size, and source volume on the concentration moments of conservative solutes in hetero- geneous formations. Advances in Water Resources 31 (2), 339–354. TRAN, M. & JHA, B. 2020 Coupling between transport and geomechanics affects spreading and mixing during viscous fingering in deformable aquifers. Advances in Water Resources 136, 103485. TRAN, M. & JHA, B. 2021 Effect of poroelastic coupling and fracture dynamics on solute trans- port and geomechanical stability. Water Resources Research 57 (10), e2021WR029584. TREFRY, M.G., LESTER, D.R., METCALFE, G., ORD, A. & REGENAUER-LIEB, K. 2012 To- ward enhanced subsurface intervention methods using chaotic advection. Journal of Contami- nant Hydrology 127 (1-4), 15–29. 132 TYUKHOVA, A., DENTZ, M., KINZELBACH, W. & WILLMANN, M. 2016 Mechanisms of anomalous dispersion in flow through heterogeneous porous media. Physical Review Fluids 1 (7), 074002. VILARRASA, V., BOLSTER, D., DENTZ, M., OLIVELLA, S. & CARRERA, J. 2010 Effects of CO 2 compressibility on CO 2 storage in deep saline aquifers. Transport in Porous Media 85 (2), 619–639. VILLERMAUX, E. 2012 Mixing by porous media. C. R. M´ ecanique 340, 933–943. WALL, C., BOERSMA, B.J. & MOIN, P. 2000 An evaluation of the assumed beta probability density function subgrid-scale model for large eddy simulation of nonpremixed, turbulent com- bustion with heat release. Physics of Fluids 12 (10), 2522–2529. WELTY, C. & GELHAR, L.W. 1991 Stochastic analysis of the effects of fluid density and viscos- ity variability on macrodispersion in heterogeneous porous media. Water Resources Research 27 (8), 2061–2075. WU, T. 1972 Cavity and wake flows. Annual Review of Fluid Mechanics 4, 243–284. YE, Y., CHIOGNA, G., CIRPKA, O.A., GRATHWOHL, P. & ROLLE, M. 2016 Experimental in- vestigation of transverse mixing in porous media under helical flow conditions. Physical Review E 94 (1), 013113. ZHANG, D., SHI, L., CHANG, H. & YANG, J. 2010 A comparative study of numerical approaches to risk assessment of contaminant transport. Stochastic Environmental Research and Risk As- sessment 24 (7), 971–984. ZHANG, P., DEVRIES, S.L., DATHE, A. & BAGTZOGLOU, A.C. 2009 Enhanced mixing and plume containment in porous media under time-dependent oscillatory flow. Environmental Sci- ence & Technology 43 (16), 6283–6288. ZHENG, C. & BENNETT, G.D. 2002 Applied contaminant transport modeling, , vol. 2. Wiley- Interscience New York. ZIMMERMAN, W.B. & HOMSY, G.M. 1991 Nonlinear viscous fingering in miscible displacement with anisotropic dispersion. Physics of Fluids A: Fluid Dynamics 3 (8), 1859–1872. ZIMMERMAN, W.B. & HOMSY, G.M. 1992 Three-dimensional viscous fingering: A numerical study. Physics of Fluids A: Fluid Dynamics 4 (9), 1901–1914. 133 Appendices A Dimensionless equations for Chapter 3 We define the following dimensionless groups: x = ˆ x λ ; u= ˆ u U ; t= ˆ t τ u ; C= ˆ C C o ; µ = ˆ µ µ 0 ; p= ˆ p p c ; k= ˆ k k c (A.1) where the (ˆ·) symbol denotes a dimensional variable. Hereλ represents the characteristic length scale (i.e. log-permeability correlation scale), U denotes the characteristic (mean longitudinal) ve- locity andτ u =λ/U corresponds to the characteristic advective time scale. The inlet concentration is C o . Furthermore, we set the ambient fluid viscosity µ 0 as a characteristic viscosity and a unitary value for the characteristic permeability k c = 1m 2 . Finally, the characteristic pressure is defined as p c =µ 0 Uλ/k c . Substituting the dimensionless groups (A.1) into Eqs. (3.3)-(3.6), we obtain: q = − k µ ∇ x p (A.2) µ = e − RC (A.3) ∇ x · q = 0, (A.4) and ∂C ∂t + q φ · ∇ x C = 1 Pe ∇ 2 x C, (A.5) where the P´ eclet number is given by Pe≡ Uλ/D andφ is the porosity of the porous medium. 134 B Dimensional equations for Chapter 4 The governing dimensional equation for the flow field is given by: ∇· ˆ q(ˆ x, ˆ t)= 0, (B.1) where the (ˆ·) symbol denotes a dimensional variable and ˆ q is the specific discharge defined by Darcy’s law as: ˆ q(ˆ x, ˆ t)=− ˆ k(ˆ x) ˆ µ( ˆ c(ˆ x, ˆ t)) (∇ ˆ p(ˆ x, ˆ t)− ∆ρgF( ˆ c(ˆ x, ˆ t)). (B.2) The dimensional governing equation for transport of an inert solute is: φ ∂ ˆ c(ˆ x, ˆ t) ∂ ˆ t + ˆ q(ˆ x, ˆ t)· ∇ ˆ c(ˆ x, ˆ t)=φD∇ 2 ˆ c(ˆ x, ˆ t). (B.3) In a homogeneous field, ˆ k(ˆ x) in Equation B.2 is a constant, ˆ k(ˆ x) = ˆ k, with corresponding dimensionless group k= ˆ k/k c , with k c = ˆ k. We then define the following dimensionless groups: x= ˆ x ℓ ; q= ˆ q Q ; t= ˆ t τ ; c= ˆ c c o ; µ = ˆ µ µ 1 ; p= ˆ p p c (B.4) whereℓ represents the characteristic length of the domain (chosen to be its thickness); Q denotes the characteristic velocity defined as Q= k c ∆ρg/µ 1 ;τ corresponds to the characteristic time scale defined as τ =µ 1 ℓφ/k c ∆ρg; p c is the characteristic pressure defined as p c =ℓ∆ρg; c o is the in- jected solute initial concentration and µ 1 is the viscosity of the injected solute. By substituting Equation B.4 into Equations B.1-B.3, we find the corresponding dimensionless equations in sec- tion 4.2, Equations 4.1, 4.2 and 4.6. 135 In a heterogeneous field, the dimensionless form of the governing equations for flow and trans- port will remain the same as presented in Equations 4.1 and 4.6, respectively; however permeability will not disappear from the dimensionless form of Darcy’s law: q q q(x,t)=− k(x) µ(c(x,t)) [∇p(x,t)− c(x,t)e e e z ]. (B.5) C Details about the verification process in Chapter 4 We test our numerical simulator by comparing its result for a fixed interface case with the result for the same case from Hewitt et al. (2013). The setup for this case consists of a zero mass and buoyancy fluxes condition at the bottom and side boundaries, while the upper boundary has zero vertical velocity and a constant concentration C ∗ m imposed. The initial concentration C ∗ − in the porous medium is constant with C ∗ − < C ∗ m , so that as time evolves the dimensional concentration C ∗ is C ∗ − < C ∗ < C ∗ m . In Hewitt et al. (2013) the dimensionless concentration C is defined in such a way that− 1< C< 0. We rescaled it in our verification process so that 0 < C< 1. Times are normalized with the convective time scale T , with T =φH/U. Lengths are normal- ized with the initial interfacial height H, which in the examined case corresponds to the domain thickness. U is the convective velocity scale, U = k∆ρg/µ, analogous to Q in Appendix A with the notable difference that in Hewitt et al. (2013) the system viscosity µ is a constant. Pressures are normalized with the characteristic pressure p c =µUH/k. Further details about the process to obtain the dimensionless governing equations can be found in section 3.2 of Hewitt et al. (2013). 136

Abstract (if available)

Abstract This thesis investigates the effect of multiple sources of disorder in solute transport in porous media through the use of numerical modeling. These sources are 1) the spatial heterogeneity of the permeability field and 2) the viscosity and density contrast between miscible fluids.
Firstly, we focus on how degree of heterogeneity of a porous medium, flow dimensionality and source zone configurations impact mixing.

Secondly, we examine how permeability heterogeneity and viscosity contrast impact the overall mixing behaviour of a solute. We analyze the degree and rate of mixing, spatial statistics of the concentration field and arrival times at a control plane to characterize spreading and mixing in the domain. Through our analysis, we are able to provide a quantitative separation of the impacts of fingering and heterogeneity. Furthermore, we parameterize the probability distribution function of the concentration field.

Thirdly, we focus on viscosity and density contrasts as sources of disorder, and we investigate their impact on the temporal evolution of the spreading and mixing quantities. We show that such impact depends on the initial shape of the source distribution where the solute is injected and on the intensity of the horizontal background flux. We also investigate how a stratified permeability field can interact with these sources of disorder and affect the transport behaviour of the plume.

Finally, we study how the uncertainty related to permeability heterogeneity and viscous instabilities affect uncertainties of quantities of interest such as early arrival times of a solute at a target location and the concentration distribution in space and time.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Efficient connectivity assessment of heterogeneous porous media using graph theory

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data

PDF

Flow and thermal transport at porous interfaces

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

PDF

Reproducible and rapid computational approaches for assessing contamination in natural aquifers

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Groundwater contaminant transport predictions in a sustainable water management scenario

PDF

Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale

PDF

Investigating the debris motion during extreme coastal events: experimental and numerical study

PDF

Evaluating the role of energy system decarbonization and land cover properties on urban air quality in southern California

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Electrokinetic transport of Cr(VI) and integration with zero-valent iron nanoparticle and microbial fuel cell technologies for aquifer remediation

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Impact of exposure to brine/CO₂ on the mechanical and transport properties of the Mt. Simon sandstone

Asset Metadata

Creator Bonazzi, Alessandra (author)

Core Title On the transport dynamics of miscible fluids in porous media under different sources of disorder

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Environmental Engineering

Degree Conferral Date 2023-08

Publication Date 07/25/2023

Defense Date 05/04/2023

Publisher University of Southern California. Libraries (digital)

Tag groundwater,Hydrology,mixing,OAI-PMH Harvest,permeability heterogeneity,porous media,viscous fingering

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor de Barros, Felipe P.J. (committee chair), Jha, Birendra (committee member), Lynett, Patrick (committee member)

Creator Email alessandrabonazzi92@gmail.com,bonazzi@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113289514

Unique identifier UC113289514

Identifier etd-BonazziAle-12150.pdf (filename)

Legacy Identifier etd-BonazziAle-12150

Document Type Dissertation

Rights Bonazzi, Alessandra

Internet Media Type application/pdf

Type texts

Source 20230726-usctheses-batch-1074 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA

Repository Email cisadmin@lib.usc.edu