We study the temporal evolution of solute dispersion in three-dimensional porous rocks of different heterogeneity and pore structure. To this end, we perform direct numerical simulations of pore-scale flow and transport in a sand-like medium, which exhibits mild heterogeneity, and a Berea sandstone, which is characterized by strong heterogeneity as measured by the variance of the logarithm of the flow velocity. Solute dispersion is quantified by effective and ensemble dispersion coefficients. The former is a measure for the typical width of the plume, the latter for the deformation, that is, the spread of the mixing front. Both dispersion coefficients evolve from the molecular diffusion coefficients toward a common finite asymptotic value. Their evolution is governed by the interplay between diffusion, pore-scale velocity fluctuations and the medium structure, which determine the characteristic diffusion and advection time scales. Dispersion in the sand-like medium evolves on the transverse diffusion time across a characteristic streamtube diameter, which is the mechanism by which pore-scale flow variability is sampled by the solute. Dispersion in the Berea sandstone in contrast is governed by both the diffusion time across a typical streamtube, and the diffusion time along a pore conduit. These insights shed light on the evolution of mixing fronts in porous rocks, with implications for the understanding and modeling of transport phenomena of microbes and reactive solutes in porous media.

We study the upscaling and prediction of dispersion in two-dimensional heterogeneous porous media with focus on transverse dispersion. To this end, we study the stochastic dynamics of the motion of advective particles that move along the streamlines of the heterogeneous flow field. While longitudinal dispersion may evolve super-linearly with time, transverse dispersion is characterized by ultraslow diffusion, that is, the transverse displacement variance grows asymptotically with the logarithm of time. This remarkable behavior is linked to the solenoidal character of the flow field, which needs to be accounted for in stochastic models for the two-dimensional particle motion. We analyze particle velocities and orientations through equidistant sampling along the particle trajectories obtained from direct numerical simulations. This sampling strategy respects the flow structure, which is organized on a characteristic length scale. Perturbation theory shows that the longitudinal particle motion is determined by the variability of travel times, while the transverse motion is governed by the fluctuations of the space increments. The latter turns out to be strongly anti-correlated with a correlation structure that leads to ultraslow diffusion. Based on this analysis, we derive a stochastic model that combines a correlated Gaussian noise for the transverse motion with a spatial Markov model for the particle speeds. The model results are contrasted with detailed numerical simulations in two-dimensional heterogeneous porous media of different heterogeneity variance.

Hyporheic zone reaction rates are highest just below the sediment-water interface, in a shallow region called the benthic biolayer. Vertical variability of hyporheic reaction rates leads to unexpected reaction kinetics for stream-borne solutes, compared to classical model predictions. We show that deeper, low-reactivity locations within the hyporheic zone retain solutes for extended periods, which delays reactions and causes solutes to persist at higher concentrations in the stream reach than would be predicted by classical approaches. These behaviors are captured by an upscaled model that reveals the fundamental physical and chemical processes in the hyporheic zone. We show how time scales of transport and reaction within the biolayer control solute retention and transformation at the stream scale, and we demonstrate that accurate assessment of stream-scale reactivity requires methods that integrate over all travel times.

We present an upscaled Lagrangian approach to predict the plume evolution in highly heterogeneous aquifers. The model is parameterized by transport independent characteristics such as the statistics of hydraulic conductivity and the Eulerian flow speed. It can be conditioned on the tracer properties, the conductivity data at the injection region, and is able to account for mobile-immobile mass transfer. Thus, the model is transferable to different solutes and hydraulic conditions. It captures the large scale non-Gaussian features for the evolution of the longitudinal mass distribution observed for the bromide and tritium tracer plumes at the macrodispersion experiment (MADE) site (Columbus, Ohio, USA), which are characterized by a slow moving peak and pronounced forward tailing. These large scale features are explained by advective tracer propagation due to a broad distribution of spatially persistent Eulerian flow speeds as a result of spatial variability in hydraulic conductivity.

Mixing along the salt-freshwater interface is critical for geochemical reactions, transport and transformation of nutrients and contaminants in coastal ecosystems. However, the mechanisms and controls of mixing are not well understood. We develop an analytical model, based on the coupling between flow deformation and dispersion, that predicts the mixing dynamics along the interface for steady state flow in coastal aquifers. The analytical predictions are compared with the results of detailed numerical simulations, which show that non-uniform flow fields, inherent to seawater intrusion in coastal aquifer, result in a non-monotonic evolution of mixing width and mixing rates along the interface. The analytical model accurately captures these dynamics over a range of freshwater flow rates and dispersivities. It predicts the evolution of the mixing width and mixing rates along the interface, offering a new framework for understanding and modeling mixing and reaction processes in coastal aquifers.

The freshwater seawater mixing zone is a critical region for chemical activity. Yet, little is known about the influence of ever present spatial heterogeneity on the dynamics of mixing and calcite dissolution, which play a key role in the understanding of Karst development. We analyse the impact of different heterogeneity structures and strengths on the local and global response of mixing and dissolution rates across the saltwater freshwater mixing zone. We find that the initial heterogeneity structure significantly impacts on observed dissolution and mixing patterns, which sheds some new light on Karst propagation in coastal aquifers.

We present an upscaled model to predict the plume evolution in highly heterogeneous alluvial aquifers. The model is parameterized exclusively by the mean, variance and correlation length of the logarithm of hydraulic conductivity, porosity and the mean hydraulic gradient. It can be conditioned on the tracer and conductivity data at the injection region. The model predicts the evolution of the longitudinal mass distribution observed at the MADE site, which is characterized by strongly non-Gaussian plume shapes with a localized peak and pronounced forward tail. The proposed model explains these features by the conductivity heterogeneity at the injection region, and tracer propagation due to a broad distribution of spatially persistent Eulerian flow speeds.

Unsaturated fractured aquifer systems offer a domain for complex gravity-driven flow dynamics leading to the development of preferential flow along fracture networks that often strongly contributes to rapid mass fluxes. This behaviour is difficult to recover by volume-effective modeling approaches (e.g. Richards equation) due to the non-linear nature of free-surface flows and mass partitioning processes at unsaturated fracture intersections. The application of well-controlled laboratory experiments enables to isolate single aspects of the mass redistribution process that ultimately affects travel time distributions across scales. We use custom-made acrylic cubes (20 cm x 20 cm x 20 cm) in analogue percolation experiments to create simple fracture networks with single or multiple horizontal fractures. A high precision multichannel dispenser produces gravity-driven free surface flow (droplets; rivulets) at flow rates ranging from 1 ml/min to 5 ml/min. Hereby, total inflow rates are kept constant while the fluid is injected via 15 (droplet flow) or 3 inlets (rivulet flow) to reduce the impact of erratic flow dynamics. Normalized fracture inflow rates (Q_f/Q_0) are calculated and compared for aperture widths d_f of 1 mm and 2.5 mm. A higher efficiency in filling an unsaturated fracture by rivulet flow observed in former studies can be confirmed. The onset of a capillary driven Washburn-type flow is determined and recovered by an analytical solution. In order to upscale the dynamics and enable the prediction of mass partitioning for arbitrary-sized fracture cascades a Gaussian transfer function is derived that reproduces the repetitive filling of fractures, where rivulet flow is the prevailing regime. Results show good agreement with experimental data for all tested aperture widths.

Infiltration processes in fractured-porous media remain a crucial, yet not very well understood component of recharge and vulnerability assessment. Under partially-saturated conditions flows in fractures, percolating fracture networks and fault zones contribute to the fastest spectrum of infiltration velocities via preferential pathways. Specifically, the partitioning dynamics at fracture intersections determine the magnitude of flow fragmentation into vertical and horizontal components and hence the bulk flow velocity and dispersion of fracture networks. In this work we derive an analytical solution for the partitioning processes based on smoothed particle hydrodynamics simulations and laboratory studies. The developed transfer function allows to efficiently simulate flow through arbitrary long wide aperture fracture networks with simple cubic structure via linear response theory and convolution of a given input signal. We derive a non-dimensional bulk flow velocity ($\widetilde{v}$) and dispersion coefficient ($\widetilde{D}$) to characterize the system in terms of dimensionless horizontal and vertical time scales $\tau_m$ and $\tau_0$. The dispersion coefficient is shown to strongly depend on the horizontal time scale and converges towards a constant value of $0.08$ within reasonable ranges for the fluid and geometrical parameters, while the non-dimensional velocity exhibits a characteristic $\widetilde{v} \sim \tau_m^{-1/2}$ scaling. Given that hydraulic information is often only available at limited places within (fractured-porous) aquifer system, such as boreholes or springs, our study intends to provide a rudimentary analytical concept to potentially reconstruct internal fracture network geometries from external boundary information, e.g., the dispersive properties of discharge (groundwater level fluctuations).

Six conceptually different models of steady groundwater flow and conservative transport are applied to the heterogeneous MADE aquifer. Their predictive capability is assessed by comparing the modelled and observed longitudinal mass distributions at different times of the plume in the MADE-1 experiment, as well as at a later time. The models differ in their conceptualization of the heterogeneous aquifer structure, computational complexity, and use of permeability data obtained from various observation methods (DPIL, Grain Size Analysis, Pumping Tests and Flowmeter). Models depend solely on aquifer structural and flow data, without calibration by transport observations. Comparison of model results by various measures, i.e. peak location, bulk mass and leading tail, reveals that the predictions of the solute plume agree reasonably well with observations if the models are underlined by a few parameters of close values: mean velocity, a parameter reflecting log-conductivity variability and a horizontal length scale related to conductivity spatial correlation. From practitioners perspective the robustness of the models is an important and useful property. The model comparison provides insight into relevant features of transport in heterogeneous aquifers. After further validation by additional field experiments or by numerical simulations, the results can be used to provide guidelines for users in selecting conceptual aquifer models, characterization strategies, quantitative models and implementation for particular goals.