On the role of mixing in controlling transport of aqueous species in heterogeneous formations

The fate of reactive solutes in groundwater is largely determined by mixing, since dilution and reactions are controlled by mixing rates. By mixing we refer to the overlap of solute bodies with a different composition, which makes possible the encounter between reacting molecules. Therefore the quantification of mixing has an important role in contamination and risk assessment and remediation technology, when they rely on processes of natural attenuation, biodegradation or chemical delivery. As porous formations are ubiquitously heterogeneous, and heterogeneity features, besides being deterministically unknown, belong to a hierarchy of scales, the description of transport processes has to deal with two main issues: epistemic uncertainty and reference scale. While the heterogeneous nature of porous media interferes with physical and chemical processes (which are inherently related to the quantification of mixing and mixing-controlled processes), the choice of the reference scale is related to the means of modeling the phenomena. In order to have an accurate representation of mixing at the continuum scale, we develop a few numerical tools, all belonging to the Lagrangian framework, and compare them with classic Eulerian and Eulerian-Lagrangian schemes. Typical transport scenarios are characterized by highly fingered plumes and sharp fringes, and pose several numerical problems (e.g. artificial diffusion and spurious oscillations). In particular, artificial diffusion can in some cases overcome the actual local dispersion, thereby possibly determining gross overestimations of reaction rates. Our numerical tests provide a set of guidelines for a conscious choice of the numerical scheme according to the objectives of the investigation and to the heterogeneity level, highlighting the drawbacks of the numerical schemes on both the evaluation of dilution and of the overall effect of reactions. Under the assumption of complete mixing at the Darcy scale, we model both instantaneous and kinetically-controlled reactive transport on synthetic bi-dimensional hydraulic conductivity fields in order to investigate the complex interplay among velocity non-uniformities, local dispersion and reaction rates at increasing levels of physical heterogeneity. We also compare the effects of different local dispersion models and injection modes (uniform vs non-uniform), still analyzing the results on a single-realization basis. Realizations share the same log-conductivity structure but are characterized by variances ranging from low (0.2) to high (10). Resorting to single-realization analysis is uncommon in the literature, unless when ergodicity conditions are fulfilled. On the other hand, ensemble analysis is insensitive to local features and does not often offer a reliable representation of actual field phenomena, especially in non-ergodic conditions. Hence single-realization scenarios can be used for understanding the key processes and their interaction, or for grasping aggregated information on the whole solute body behavior. Under simplified conditions, that is, limiting the investigation to low heterogeneity fields, these numerical results are compared to simplified Lagrangian semianalytical relations aiming at reproducing plume-averaged quantities. This Lagrangian theory provides relevant information relying on a limited amount of information, i.e. low-order geostatistical properties of the formation, aquifer's geometry, reactive parameters and problem forcings (e.g. initial and boundary conditions for the flow field and the concentration of the involved species). The match between empirical and theoretical global moments is very good in all tested conditions (two different Peclet numbers, a few heterogeneity levels up to log-transmissivity variance equal to 2 and three different source sizes), and also Beta Cumulative Frequency Distributions (CFDs) with shape parameters obtained by substituting the theoretical global moments compare well with the numerical CFDs. As expected, coherent estimates of peak concentration are not equally good, because of an inherently different nature of this quantity as opposed to plume-scale concentration moments. The a-priori information expressed by statistical analysis both at the global scale and at the local scale for a conservative tracer z can be transferred to reactive species in case of very fast kinetics. Given this useful property of equilibrium reactions, we develop explicit semianalytical relations for the moments and the probability distribution functions of the concentration of chemical species reacting according to a bimolecular equilibrium homogeneous reaction. We assume that the conservative tracer probability distribution function, both at the local scale and at the global scale, can be modeled with a Beta distribution, fully characterized by the mean and the variance of z. Rigorous numerical testing on highly heterogeneous velocity fields confirms that this assumption holds. A few illustrative cases shed some light on the role of the reaction on the time evolution of (local and global) concentration for the different reactive species, and on the different quality of information contained in local statistics as opposed to global statistics. The Beta distribution is a powerful predicting tool for the space and time evolution of passive concentration and, by extension, also for reactive species in particular chemical conditions. Analytical procedures are needed for predicting the z moments, as for example the Lagrangian ones used in the present work, which are limited to weakly heterogeneous formations. Finally we explore, analytically and numerically, the upscaling from the pore scale to the Darcy scale. Via multiple scale analysis we identify a homogenizability region, in terms of the dimensionless numbers regulating a multicomponent precipitation/dissolution reactive problem, where Darcy-scale (upscaled) transport equations can be used, regardless of sub-Darcy scale inhomogeneities.