Search

Abstract: The selective radius shift model was used to relate changes in mineral volume due to precipitation/dissolution reactions to changes in hydraulic properties affecting flow in porous media. The model accounts for (i) precipitation/dissolution taking place only in the water-filled part of the pore space and further that (ii) the amount of mineral precipitation/dissolution within a pore depends on the local pore volume. The pore bundle concept was used to connect pore-scale changes to macroscopic soil hydraulic properties. Precipitation/dissolution induces changes in the pore radii of water-filled pores and, consequently, in the effective porosity. In a time step of the numerical model, mineral reactions lead to a discontinuous pore-size distribution because only the water-filled pores are affected. The pore-size distribution is converted back to a soil moisture characteristic function to which a new water retention curve is fitted under physically plausible constraints. The model equations were derived for the commonly used van Genuchten/Mualem hydraulic properties. Together with a mixed-form solution of Richards’ equation for aqueous phase flow, the model was implemented into the geochemical modelling framework PHREEQC, thereby making available PHREEQC’s comprehensive geochemical reactions. Example applications include kinetic halite dissolution and calcite precipitation as a consequence of cation exchange. These applications showed marked changes in the soil’s hydraulic properties due to mineral precipitation/dissolution and the dependency of these changes on water contents. The simulations also revealed the strong influence of the degree of saturation on the development of the saturated hydraulic conductivity through its quadratic dependency on the van Genuchten parameter α. Furthermore, it was shown that the unsaturated hydraulic conductivity at fixed reduced water content can even increase during precipitation due to changes in the pore-size distribution. [Copyright &y& Elsevier]

Abstract: This paper concerns efficient uncertainty quantification techniques in inverse problems for Richards’ equation which use coarse-scale simulation models. We consider the problem of determining saturated hydraulic conductivity fields conditioned to some integrated response. We use a stochastic parameterization of the saturated hydraulic conductivity and sample using Markov chain Monte Carlo methods (MCMC). The main advantage of the method presented in this paper is the use of multiscale methods within an MCMC method based on Langevin diffusion. Additionally, we discuss techniques to combine multiscale methods with stochastic solution techniques, specifically sparse grid collocation methods. We show that the proposed algorithms dramatically reduce the computational cost associated with traditional Langevin MCMC methods while providing similar sampling performance. [Copyright &y& Elsevier]

• Equivalent hydraulic conductivity is studied using a stochastic approach in 2D and 3D. • A wide range of connectivity scenarios is portrayed using binary media. • Strong changes in connectivity derive in a displacement of the percolation threshold. • Data collapse for all considered cases, including 2D and 3D data. • Useful for applications related to reservoir connectivity. The equivalent hydraulic conductivity (K e q) relates the spatial averages of flux and head gradient in a block of heterogeneous media. In this article, we study the influence of connectivity on K e q of media samples composed of a high conductivity (k +) and a low conductivity (k −) facies. The k + facies is characterized by a proportion p , and also by two connectivity parameters: a connectivity structure type (no, low, intermediate, high), and a correlation integral scale l c. The probability distribution of K e q , and the critical value of p at which percolation occurs (p a v), are studied as a function of these connectivity parameters. The distribution of log (K e q) is Gaussian in all cases, so the results are presented in terms of the geometric mean (〈 K e q 〉) and the variance (σ log (K e q) 2). Both quantities show a data collapse if expressed as a function of p − p a v (for the variance σ log (K e q) 2 , notably, even if 2D and 3D data are plotted together). In 3D, when a connectivity structure exists, K e q is always greater than when no structure exists, and increases (while p a v decreases) as l c increases. The same is observed in 2D, except for the low connectivity structure type (i.e. when the k + facies is disconnected), that shows an unprecedented behaviour: K e q is greater in the absence of structure, and decreases (p a v increases) as l c increases. Our results show that any influence of connectivity on K e q is well accounted for simply by a shift in the percolation threshold p a v , and then, suggest that K e q is controlled mainly by the proximity to percolation. [Display omitted] [ABSTRACT FROM AUTHOR]

Abstract: Unsaturated water flow through soil aggregates is controlled by the contacts between aggregates. The contacts are highly conductive when wet and become bottle-necks for flow when drained. We postulate that the hydraulic conductivity of the contacts is in first place determined by the water-filled contact area. The objective of this study was to measure and model the water-filled contact area and to relate it to the conductivity of a series of aggregates. We performed microscopic tomography of an aggregate pair equilibrated at different water potentials. By means of image analysis and a morphological pore network model, the water-filled contact area was calculated. We found that the aggregate surface is rough and the contact region contains macropores which are rapidly drained. As a consequence the water-filled contact area dramatically decreases as the water potential is diminished. We modeled this process by describing the aggregates as spheres covered by much smaller spheres representing the roughness. The water-filled contact was analytically calculated from this model. Knowing the water-filled contact area we up-scale the hydraulic conductivity of a series of aggregates. This is calculated as the harmonic mean of the contact and aggregate conductivities. The contact conductivity is calculated from the water-filled contact area. Near saturation the conductivity of a series of aggregates is close to the conductivity of a single aggregate, and, when further drained, it rapidly decreases as the water-filled contact area. The model matches the experimental data well. [Copyright &y& Elsevier]

Abstract: Investigations in outcrops of non-homogeneous gravel deposits in north-eastern Switzerland by Jussel et al. [Jussel P, Stauffer F, Dracos T. Transport modeling in heterogeneous aquifers: 1. Statistical description and numerical generation of gravel deposits. Water Resour Res 1994;30:1803–17] revealed distinct sedimentary structures or facies elements appearing as more or less horizontal lenses and layers embedded in extended patches of background gravel material. One of the observed sediment units consists of highly permeable inclined sequences of gravel layers. This inclination manifests itself as a strong anisotropy of the hydraulic conductivity with inclined principal directions. In this paper, the role of highly conductive sediment units with inclined bedding is investigated by a series of three-dimensional numerical flow and advective transport experiments. A typical single lens with inclined principal axes of the hydraulic conductivity tensor is investigated, which is embedded in homogeneous background matrix material. Particle tracking shows pronounced lateral (horizontal and vertical) distortions of the flow path depending on the angles of the bedding of the lens. The results allow better understanding of the flow and transport phenomena that are finally responsible for macrodispersion effects in similar gravel aquifers. [Copyright &y& Elsevier]

Abstract: The hydraulic conductivity of heterogeneous porous media depends on the distribution function and the geometry of local conductivities at the smaller scale. There are various approaches to estimate the effective conductivity K eff at the larger scale based on information about the small scale heterogeneity. A critical geometric property in this ‘upscaling’ procedure is the spatial connectivity of the small-scale conductivities. We present an approach based on the Euler-number to quantify the topological properties of heterogeneous conductivity fields, and we derive two key parameters which are used to estimate K eff. The required coefficients for the upscaling formula are obtained by regression based on numerical simulations of various heterogeneous fields. They are found to be generally valid for various different isotropic structures. The effective unsaturated conductivity function K eff (ψ m) could be predicted satisfactorily. We compare our approach with an alternative based on percolation theory and critical path analysis which yield the same type of topological parameters. An advantage of using the Euler-number in comparison to percolation theory is the fact that it can be obtained from local measurements without the need to analyze the entire structure. We found that for the heterogeneous field used in this study both methods are equivalent. [Copyright &y& Elsevier]

Abstract: Formation properties are one of the key factors in numerical modeling of flow and transport in geologic formations in spite of the fact that they may not be completely characterized. The incomplete knowledge or uncertainty in the description of the formation properties leads to uncertainty in simulation results. In this study, the ensemble Kalman filter (EnKF) approach is used for continuously updating model parameters such as hydraulic conductivity and model variables such as pressure head while simultaneously providing an estimate of the uncertainty through assimilating dynamic and static measurements, without resorting to the explicit computation of the covariance or the Jacobian of the state variables. A two-dimensional example is built to demonstrate the capability of EnKF and to analyze its sensitivity with respect to different factors such as the number of realizations, measurement timings, and initial guesses. An additional example is given to illustrate the applicability of EnKF to three-dimensional problems and to examine the model predictability after dynamic data assimilation. It is found from these examples that EnKF provides an efficient approach for obtaining satisfactory estimation of the hydraulic conductivity field with dynamic measurements. After data assimilation the conductivity field matches the reference field very well, and different kinds of incorrect prior knowledge of the formation properties may also be rectified to a certain extent. [Copyright &y& Elsevier]

Abstract: An existing capillarity correction for free surface groundwater flow as modelled by the Boussinesq equation is re-investigated. Existing solutions, based on the shallow flow expansion, have considered only the zeroth-order approximation. Here, a second-order capillarity correction to tide-induced watertable fluctuations in a coastal aquifer adjacent to a sloping beach is derived. A new definition of the capillarity correction is proposed for small capillary fringes, and a simplified solution is derived. Comparisons of the two models show that the simplified model can be used in most cases. The significant effects of higher-order capillarity corrections on tidal fluctuations in a sloping beach are also demonstrated. [Copyright &y& Elsevier]

The total travel distance and the hydraulic conductivity in steady-state groundwater flow systems are modeled with asymmetric gamma distribution functions. The ratio of the total travel distance over the hydraulic conductivity, scaled by a factor involving the hydraulic gradient and the aquifer porosity, equals the residence time of groundwater. The residence time measures the period elapsed between recharge and discharge of groundwater from an aquifer flow system. The residence time is a random variable whose probability density function (p.d.f.) is determined by those of the total travel distance and the hydraulic conductivity. The p.d.f. of the residence time was derived and from that, applying basic mass-balance considerations, followed that of the groundwater age. The age is the time elapsed since the recharge of groundwater in transit through an aquifer flow system. It is shown that the p.d.f.s of the residence time and the groundwater age are related by a first-order, linear, differential equation. The cumulative distribution functions of the residence time and the groundwater age and their first and second moments are also derived, and so is the output concentration of a solute that enters a steady-state groundwater flow system dissolved in recharge. A computational example illustrates the advantage of treating the residence time and the groundwater age statistically and constraining them with Darcy''s law and mass-balance principles. Solute output shows a complex structure with stable and radioactive solutes in deceptively simple steady-state groundwater flow regimes. [Copyright &y& Elsevier]

The Boussinesq equation appears as the zeroth-order term in the shallow water flow expansion of the non-linear equation describing the flow of fluid in an unconfined aquifer. One-dimensional models based on the Boussinesq equation have been used to analyse tide-induced water table fluctuations in coastal aquifers. Previous analytical solutions for a sloping beach are based on the perturbation parameter, ϵN=αϵcotβ (in which β is the beach slope, α is the amplitude parameter and ϵ is the shallow water parameter) and are limited to tan−1(αϵ)≪β⩽π/2. In this paper, a new higher-order solution to the non-linear boundary value problem is derived. The results demonstrate the significant influence of the higher-order components and beach slope on the water table fluctuations. The relative difference between the linear solution and the present solution increases as ϵ and α increase, and reaches 7% of the linear solution. [Copyright &y& Elsevier]

• Finite particle method (FPM) simulates groundwater solute transport. • FPM outperforms smoothed particle hydrodynamics (SPH) with irregular particles. • FPM improves SPH by using a modified kernel gradient for concentration gradient. When a smoothed particle hydrodynamics (SPH) method, a Lagrangian and meshfree numerical scheme, is used to solve the advection-dispersion equation (ADE), SPH solutions may not be accurate when particles are irregularly distributed, i.e., the distance between two neighbor particles varies irregularly in a simulation domain. Particle irregularity may be caused by nonuniform groundwater flow in a heterogeneous field of hydraulic conductivity. This study explores for the first time whether the Finite Particle Method (FPM) can provide more accurate ADE solutions than SPH does for irregularly distributed particles. FPM is similar to SPH in theory, but uses a modified kernel gradient to construct a SPH approximation of solute concentration gradients. Performance of SPH and FPM with irregularly distributed particles is evaluated by using two numerical cases. The first case considers only diffusive transport, and has an analytical solution for the evaluation. The second case considers both advection and dispersion, and uses a numerical solution as a reference for the evaluation. For each of the two cases, several numerical experiments are conducted using multiple sets of irregularly distributed particles with different levels of particle irregularity due to different levels of heterogeneity of hydraulic conductivity. Numerical results indicate that, for the numerical experiments of this study, FPM outperforms SPH to yield more accurate ADE solutions. However, FPM solutions are still subject to numerical errors, and the errors increase when the level of heterogeneity of hydraulic conductivity increases. Further improvement of FPM is warranted in a future study. [ABSTRACT FROM AUTHOR]

Three-dimensional (3D) natural convection (NC) processes in heterogeneous porous media and associated energy losses and mixing processes are still poorly understood. Studies are limited to two-dimensional domains because of computational burden, worsened by heterogeneity, which may demand grid refinement at high permeability zones for accurate evaluation of buoyancy forces. We develop a meshless Fourier series (FS) solution of the natural convection problem in a porous enclosure driven by thermal or compositional variations. We derive the vector potential formulation of the governing equations for vertical and horizontal heterogeneity of hydraulic conductivity and implement an efficient method to solve the spectral system with an optimized number of Fourier modes. 3D effects are induced either by heterogeneity or variable boundary conditions. The developed FS solution is verified against a finite element solution obtained using COMSOL Multiphysics. We evaluate entropy generation (viscous dissipation and mixing) indicators using FS expansions and assess how they are affected by heterogeneity. We define a large-scale Rayleigh number to account for heterogeneity by adopting an arithmetic average effective permeability. The FS solution is used to investigate the effect of the large-scale Rayleigh number and level of heterogeneity on NC processes and energy losses. Results show that increasing the Rayleigh number intensifies fluid flow, thus enhancing convective transfer, which causes a dramatic increase in total entropy generation. Both viscous dissipation and mixing (and thus chemical reactions in the solute transport case) increase. The third dimension effect, which also enhances flow and entropy indicators, is more pronounced at high Rayleigh numbers. Surprisingly, entropy variation indicators remain virtually unchanged in response to changes in heterogeneity, for fixed Rayleigh number, which we attribute to the arithmetic average permeability being indeed appropriate for NC in 3D. This study not only explores the effect of Rayleigh number and heterogeneity on natural convection processes and the associated entropy generation and mixing processes, but also provides a highly accurate solution that can be used for codes benchmarking. [ABSTRACT FROM AUTHOR]

• The PCR-GLOBWB model was calibrated for 10-year rolling periods for 1901 to 2010. • The optimal parameter values changed notably through time in the Rhine-Meuse basin. • Changes in optimal parameter values are larger at the upstream parts of the basin. • Parameter changes are highly correlated with land use change and climate change. • Systemic change has occurred in the upstream parts of the Rhine-Meuse basin. In hydrological modelling, traditionally one calibration was performed over a certain calibration period before the model is used to study the hydrological system. This implies that a constant model structure and parameterization are assumed. However, if the catchment system is subject to changes that are not incorporated in the model, the parameter values found in a calibration period may not be optimal for other periods, which is called systemic change. The aim of this study was to identify systemic change and its possible causes with the PCR-GLOBWB hydrological model in the Rhine-Meuse basin, by performing a brute-force calibration for multiple periods for five calibration locations between 1901-2010. Systemic change was studied for the main model components, by selecting a key parameter from each component (minimum soil depth fraction, saturated hydraulic conductivity, groundwater recession coefficient, degree day factor, Manning's n). These parameters were calibrated for 10-year rolling periods between 1901-2010. The results showed that at the downstream locations, the changes in optimal parameter values were small, while at the upstream locations, the optimal values of most parameters changed considerably over the different rolling calibration periods, signifying systemic change. Especially the degree day factor showed large variations, varying over time between 0.5 and 2.5 times its default value at Basel and Maxau (upstream and middle part of the Rhine basin). Based on correlation analysis, it was found that climate change as well as changes in land use and river structure are possible causes of changes in optimal parameter values through time. [ABSTRACT FROM AUTHOR]

• Geostatistical inversion of head and flux data recorded under steady-state conditions. • For few observations, flux data is more informative than head data. • For many observations, the inversion results are comparable regardless of data type. • For equal number of observations, individual and joint inversion perform similarly. • The pumping borehole boundary condition affects individual, but not joint inversion. Hydraulic tomography is a state-of-the-art method for inferring hydraulic conductivity fields using head data. We employed geostatistical inversion using synthetically generated head and flux data individually and jointly in a steady-state experiment. We designed 96 inversion scenarios to better understand the relative merits of each data type. For the typical case of a small number of observation points, we find that flux data provide a better resolved hydraulic conductivity field compared to head data when considering data with similar signal-to-noise ratios. This finding is further confirmed by a resolution analysis. When considering a high number of observation points, the estimated fields are of similar quality regardless of the data type. In terms of borehole boundary conditions, the best setting for flux and head data are constant head and constant rate, respectively, while joint inversion results are insensitive to the borehole boundary type. When considering the same number of observations, the joint inversion of head and flux data does not offer advantages over individual inversions. When considering the same number of observation points and, hence, twice as many observations, the joint inversion performs better than individual inversions. The findings of this paper are useful for future planning and design of hydraulic tomography tests comprising flux and head data. [ABSTRACT FROM AUTHOR]

• Novel methodology for physically and dynamically constraining Kosugi hydraulic parameters to reduce non-uniqueness. • Development of a parsimonious bimodal Kosugi hydraulic model. • Derive bimodal Kosugi hydraulic parameters exclusively from θ (ψ) and K s data. Hydraulic parameters define the water retention, θ (ψ), and the unsaturated hydraulic conductivity, K (θ), functions. These functions are usually obtained by fitting experimental data using inverse modelling. The drawback of inverting the hydraulic parameters is that they suffer from non-uniqueness and the optimal hydraulic parameters may not be physical. To reduce the non-uniqueness, it is necessary to invert the hydraulic parameters simultaneously from observations of θ (ψ) and K (θ) , and ensure the measurements cover the full range of θ from saturated to oven dry. The challenge of using bimodal θ (ψ) and K (θ) compared to unimodal functions is that it requires double the number of parameters, one set for the matrix and another set for the macropore domain. The objective of this paper is to address this shortcoming by deriving a procedure to reduce the number of parameters to be optimized to obtain a unique physical set of bimodal soil Kosugi hydraulic parameters from inverse modelling. To achieve this, we (1) derive residual volumetric soil water content from the Kosugi standard deviation parameter of the soil matrix, (2) derive macropore hydraulic parameters from the water pressure head threshold between macropore and matrix flow, and (3) dynamically constrain the Kosugi hydraulic parameters of the soil matrix. The procedure successfully reduces the number of optimized hydraulic parameters and dynamically constrains the hydraulic parameters without compromising the fit of the θ (ψ) and K (θ) functions, and the derived hydraulic parameters are more physical. The robustness of the methodology is demonstrated by deriving the hydraulic parameters exclusively from θ (ψ) and K s data, enabling satisfactory prediction of K (θ) even when no additional K (θ) data are available. [ABSTRACT FROM AUTHOR]

• We propose a deep neural network to predict 2D subsurface property fields from 1D measurement data (e.g., time series) thereby offering a new way to perform inversion. • The method is illustrated using inversion of synthetic (1) first-arrival GPR travel times and (2) time series of transient hydraulic heads from a multi-pumping experiment. • The method generalizes by producing more appropriate models when fed with data than those found in the training set. • Uncertainty is assessed using a deep ensemble. We present vec2pix , a deep neural network designed to predict categorical or continuous 2D subsurface property fields from one-dimensional measurement data (e.g., time series), thereby offering an alternative approach to solve inverse problems. The performance of the method is investigated through two types of synthetic inverse problems: (a) a crosshole ground penetrating radar (GPR) tomography experiment with GPR travel times being used to infer a 2D velocity field, and (2) a multi-well pumping experiment within an unconfined aquifer with time series of transient hydraulic heads being used to retrieve a 2D hydraulic conductivity field. For each type of problem, both a multi-Gaussian and a binary channelized subsurface domain with long-range connectivity are considered. Using a training set of 20,000 examples (implying as many forward model evaluations), the method is found to recover a 2D model that is in much closer agreement with the true model than the closest training model in the forward-simulated data space. Further testing with smaller training sample sizes shows only a moderate reduction in performance when using 5000 training examples only. Although these 5000 to 20,000 forward runs can be performed offline in parallel, the associated computational expense may still be prohibitive for very demanding forward solvers. In addition, re-training is required for each new measurement configuration. Even if the recovered models are visually close to the true ones, the data misfits associated with their forward responses are generally larger than the noise level used to contaminate the true data. Uncertainty of the inverse solution is partially assessed using deep ensembles, in which the network is trained repeatedly with random initialization. Overall, this study advances understanding of how to use deep learning to infer subsurface models from indirect measurement data. More work is needed to evaluate the suitability of our proposed approach for real data, for which model errors and uncertainty in the geologic scenario will bring additional complexities. [ABSTRACT FROM AUTHOR]

• Flow and transport equations for porous media are solved by new random walk schemes. • Global random walk (GRW) algorithms are designed as L-schemes for Richards equation. • The GRW L-schemes solve nonlinearities and model unsaturated/saturated transitions. • The coupled flow-transport solvers are robust and second-order convergent in space. • GRW schemes for both flow and transport are practically free of numerical diffusion. In this article, we present new random walk methods to solve flow and transport problems in saturated/unsaturated porous media, including coupled flow and transport processes in soils, heterogeneous systems modeled through random hydraulic conductivity and recharge fields, processes at the field and regional scales. The numerical schemes are based on global random walk algorithms (GRW) which approximate the solution by moving large numbers of computational particles on regular lattices according to specific random walk rules. To cope with the nonlinearity and the degeneracy of the Richards equation and of the coupled system, we implemented the GRW algorithms by employing linearization techniques similar to the L -scheme developed in finite element/volume approaches. The resulting GRW L -schemes converge with the number of iterations and provide numerical solutions that are first-order accurate in time and second-order in space. A remarkable property of the flow and transport GRW solutions is that they are practically free of numerical diffusion. The GRW solvers are validated by comparisons with mixed finite element and finite volume solvers in one- and two-dimensional benchmark problems. They include Richards' equation fully coupled with the advection-diffusion-reaction equation and capture the transition from unsaturated to saturated flow regimes. [ABSTRACT FROM AUTHOR]

• Log-conductivity (Y) is modeled as a generalized sub-Gaussian (GSG) field. • Average longitudinal spreading is smaller for GSG than for Gaussian Y fields. • Plume stretching in GSG fields is enhanced with respect to Gaussian fields. • Local dispersion masks identification of impacts of non-Gaussianity on transport. • The GSG nature of Y can be difficult to identify relying only on transport metrics. There are increasing evidences that probability distributions and associated statistical moments of a variety of hydrogeological and soil science variables and their spatial increments display distinctive scale-dependent features that are not captured by a typical Gaussian model. A Generalized Sub-Gaussian (GSG) model is able to capture key aspects of this pattern. We present the results of a suite of computational analyses set in a Monte Carlo framework and aimed at assessing the impact of a GSG structure of log hydraulic conductivity (Y) on transport of a conservative solute through a three-dimensional bounded porous medium under steady-state saturated Darcy flow. Our results indicate that the longitudinal spreading of a plume is on average significantly smaller for Sub-Gaussian than for Gaussian Y fields. Otherwise, the velocity field arising from a Sub-Gaussian Y field induces enhanced plume stretching with respect to what can be observed in a Gaussian Y setting, this aspect potentially influencing the strength of solute mixing within these two types of conductivity domains. We also find that, in some cases, it may be difficult to identify the nature of the underlying conductivity field relying solely on observations of solute concentrations migrating within the system. In this regard, we show that the action of local dispersion tends to mask the influence of Sub-Gaussianity on major transport metrics. [ABSTRACT FROM AUTHOR]

To consider the non-Darcian nature of groundwater flow and the permeability reduction owing to well and aquifer clogging, an Izbash's law-based analytical model with a time-dependent hydraulic conductivity function is proposed for constant-rate recharge and constant-head recharge. By means of linearization approximations, variable substitutions and the Laplace transform, approximate analytical solutions for the model are derived. Comparison of developed approximate analytical solutions with numerical solutions indicates that the approximate solutions are highly accurate except during the short period of initial recharge stage. The proposed model provides a theoretical basis for managed aquifer recharge involving clogging-related permeability reduction and non-Darcian flow. [ABSTRACT FROM AUTHOR]

• Geostatistical K-field properties inferred from electrical conductivity time-series. • Approximate Bayesian computation assessed with Kullback–Leibler (KL) measure. • Largest KL measure for time-series of horizontal equivalent electrical conductivity. • The variance of the log-transformed K-field is the best constrained parameter. We use Approximate Bayesian Computation and the Kullback–Leibler divergence measure to quantify to what extent horizontal and vertical equivalent electrical conductivity time-series observed during tracer tests constrain the 2-D geostatistical parameters of multivariate Gaussian log-hydraulic conductivity fields. Considering a perfect and known relationship between salinity and electrical conductivity at the point scale, we find that the horizontal equivalent electrical conductivity time-series best constrain the geostatistical properties. The variance, controlling the spreading rate of the solute, is the best constrained geostatistical parameter, followed by the integral scales in the vertical direction. We find that horizontally layered models with moderate to high variance have the best resolved parameters. Since the salinity field at the averaging scale (e.g., the model resolution in tomograms) is typically non-ergodic, our results serve as a starting point for quantifying uncertainty due to small-scale heterogeneity in laboratory-experiments, tomographic results and hydrogeophysical inversions involving DC data. [ABSTRACT FROM AUTHOR]

• A physically-based model is presented to calculate hydraulic conductivity of stony soils. • Rock fragments can be permeable and effective saturated and unsaturated hydraulic conductivity of a stony soil can be calculated for all matric potentials. • The model considers the effect of shape and orientation of rock fragments. • The model is also valid for high volumetric rock fragment contents where existing approaches are not applicable. Despite a longstanding interest in understanding the effects of rock fragments on soil hydraulic properties, physically-based models to describe the effective hydraulic conductivity of stony soils are still rare. The conductivity of stony soils is mostly described by simple scaling approaches, assuming impermeable rock fragments. We present a model based on the general effective medium (GEM) theory to calculate the effective hydraulic conductivity of stony soils. The model accounts for the volumetric content, shape, and orientation of rock fragments in the main direction of water flow. Rock fragments can either be permeable or impermeable in the model which allows to calculate the effective unsaturated hydraulic conductivity curve of stony soils over the full range of soil matric potentials. Existing models assume rock fragments to be isolated and their application is therefore limited to low to moderately stony soils. In GEM model interaction between rock fragments is allowed, which makes it applicable in stony soils with high contents of rock fragments. The mathematical simplicity of the model makes it suitable for practical applications and provides a strong theoretical basis for improving predictions of flow and transport in hydrological models applied in forested and agricultural watersheds and mountainous areas. [ABSTRACT FROM AUTHOR]

• Variations in the support/observation scale of hydraulic conductivity induce changes in the information content of the associated Darcy fluxes. • Propagation of information about conductivity onto Darcy flux through the flow and mass balance equation can be quantified via information theory. • Flux data at large scales provide more information about fluxes at small scales than what can be observed with reference to conductivities. We rest on an Information Theory perspective and assess (i) the average information content and (i) the information shared between Darcy flux fields associated with fields of hydraulic conductivity (K) characterized by differing support (or measurement) scale. We treat hydraulic conductivity as a spatial random field, characterized by a given distribution and correlation structure. The latter is modeled through a truncated power law variogram (TPV), which explicitly takes into account a characteristic length scale of the support volume of K through a lower cutoff scale. We then frame our study in a numerical Monte Carlo context where groundwater flow is evaluated across a collection of realizations of hydraulic conductivity characterized by different values of the TPV parameters and subject to uniform in the mean flow. We quantify information through the Shannon entropy of the probability mass functions of the Darcy flux components as well as the mutual information and the multivariate mutual information respectively shared by pairs and triplets of Darcy flux components related to hydraulic conductivity fields evaluated at diverse scales and associated with various levels of heterogeneity. Partitioning of multivariate mutual information according to unique, redundant and synergetic contributions is also quantified. We found consistent trends (i) in the variation of the average information content with respect to the size of lower cutoff scale and (ii) in the way information is shared between pairs and triplets of Darcy flux components associated with diverse support scales of the underlying conductivities. [ABSTRACT FROM AUTHOR]

• A general mathematical model describing non-Darcian flow in a confined aquifer for a single well circulation system was developed. • This study provided an extension of the analytical solution based on Darcy's Law to non-Darcian flow conditions for a single well circulation system. • An approximate analytical solution for non-Darcian flow in the Laplace domain and the steady-state analytical solution in the time domain were obtained. • The characteristics of non-Darcian flow and the influences of different parameters on drawdown were investigated. An analytical model is proposed in this study to describe transient drawdown induced by non-Darcian flow in a confined aquifer with a single well circulation groundwater heat pump system. The Izbash equation and a linearization method are employed to describe non-Darcian flow in the horizontal direction of a confined aquifer and to approximate the nonlinear term in the governing equation, respectively. By applying a combination of the Laplace and Fourier cosine transforms, an approximate analytical solution in the Laplace domain is obtained, which is numerically inverted to obtain transient drawdown in the time domain using the Stehfest algorithm method. The results of the derived analytical solution for the special case of Darcian flow (m = 1) correspond well with the existing solution derived using Darcy's law. The steady-state analytical solution in the time domain is obtained by applying the Fourier cosine transform. Moreover, the sensitivity analysis is performed to investigate the influence of selected parameters, such as the power index m , the radial hydraulic conductivity K r , the aquifer specific storage S , and the length of the sealed section d 2 , on drawdown. The results show that each parameter has its influence period on drawdown, and that drawdown is more sensitive to the power index m compared to other parameters. [ABSTRACT FROM AUTHOR]

Hydraulic barriers using either two or three injection wells are studied analytically to evaluate their performance on improving groundwater extraction in a coastal aquifer under the threat of seawater intrusion. The multiple well system is operated by extracting fresh groundwater through a landward well and partially reinjecting into the aquifer through two or three wells located parallel to the coastline and between the interface toe and the extraction well. We found that with appropriate well locations and pumping rates, the multiple-well system could outperform a single extraction well in terms of net extraction rate, which is up to 90% higher for an aquifer with a thickness of 20 m, hydraulic conductivity of 10 m/d, and freshwater influx of 0.24 m2/d. Sensitivity analyses reveal that the monotonicity of the maximum net extraction rate is determined by the total injection rate. In specific, when the total injection rate is fixed, the maximum net extraction rate is a non-monotonic function of the distance between the injection well and the coastline and the distance between two injection wells, resulting in a global maximum. When the total injection rate is not fixed, by contrast, the maximum net extraction rate depends on the total injection rate, and a larger total injection rate leads to a larger maximum net extraction rate. By using different injection rates for side and middle injection wells of the three-injection-well hydraulic barrier, the maximum net extraction rate could be further improved, depending on the distance between the two injection wells. Analytical analyses can provide initial guidance for the design of a multiple-well system in coastal aquifers with the purpose of obtaining enhanced extraction under the threats of seawater intrusion. [ABSTRACT FROM AUTHOR]

• Represents local heterogeneity at micro-scale level with an average error variance less than 0.02. • The FSM shows consistent results at well node locations with changing grid dimensions. • ROM is computationally efficient and replicates FSM with desired accuracy. • Both FSM and ROM shows good mesh convergence. Groundwater, being a vital component of the natural water resource system, needs continuous monitoring and dynamic management strategies. That said, we require computationally inexpensive groundwater flow models for repetitive solutions with desirable accuracy under budgetary limitation(s). Natural aquifer systems inherit strong heterogeneity at local scales. In this work, we have proposed ordinary kriging-based sequential algorithm for generating replicates of randomly distributed heterogeneous hydraulic conductivity field (Monte Carlo method-based algorithm) conditioned by field values from sampled locations in an irregular-unstructured grid system. Finite Volume method-based groundwater models often encounter difficulties with the representation of point source/sink terms operating within the domain. In this paper, we have proposed an irregular-unstructured grid Finite Volume discretization technique for overcoming the singularity of point source/sink term to yield a consistent output with different grid dimensions. Furthermore, full-system groundwater models often come with a substantial computational burden. Hence, reduction in model order cuts down the computational expenses (in terms of CPU time and usage) to a significant level. We have also put forth a model order reduction methodology for three different illustrative pumping tests. The proposed framework for the model order reduction projects the governing groundwater flow equation onto a set of identified patterns or orthonormal basis functions, applying the Galerkin Projection method to compute a vector of time-dependent coefficients. We have performed pattern identification by Singular Value Decomposition (SVD) of snapshots of full-system model simulation data at selected time instants within the pumping test time domain. The numerical results of the proposed reduced-order models show a good approximation of the full-system models at a comparatively lesser computational time. The accuracy and efficiency of the models attempt to ensure their potential applicability for identifying groundwater dynamics. [ABSTRACT FROM AUTHOR]

• An integro-differential equation for upscaling reactive transport in a correlated velocity field is derived. • The proposed approach provides a way to predict reactive-transport based on the geostatistical properties of the media. • Breakthrough curves of a non-linearly decaying solute from the upscaled model are compared with the ones from a 2D model. While commonly non-local transport models have been shown to reproduce breakthrough curves resulting from transport in heterogeneous aquifers successfully, open questions include the formal link between the upscaled governing equations and the sub-scale heterogeneity, and the ability to account for the effect of heterogeneity on effective chemical reaction rates in the presence of non-linear multi-component reactions. Time-domain random walk approaches based on velocity Markov models provide a framework to resolve these issues by incorporating the spatial correlation of velocity of a solute particle in consecutive locations along its trajectory. These approaches often rely on particle tracking approaches, which, however, can be computationally burdensome especially when non-linear reactions are sought to be modeled. In this paper, an integro-differential equation is proposed for upscaling multi-component reactive transport in heterogeneous media that relies on copulas for representing the velocity correlation structure. For this purpose, we express concentration or flux of solutes as a distribution over their velocity. We then derive the integro-differential equation that governs the evolution of concentration distribution over a quantity defined as velocity-rank. In this way, the spatial evolution of breakthrough curves away from the source is predicted based on ergodic cross-sectional velocity distributions and a parameterized copula function which expresses the correlation between velocity-ranks of a solute particle along its trajectory. We demonstrate the validity of the proposed model by comparing breakthrough curves for conservative and non-linearly reacting solutes based on realizations of hydraulic conductivity fields to the results of the upscaled model. [ABSTRACT FROM AUTHOR]

• Steady source flow in aquifer of lognormal permeability. • The equivalent conductivity K eq quantifies heterogeneity. • The dependence of Keq on distance and properties is determined. • K eq tends rapidly to the effective conductivity in uniform flow. • K eq expression can be conveniently used in the inversion of hydraulic tomography data. The effective conductivity of heterogeneous aquifers K efu relates the mean flux to the mean head gradient in mean uniform (natural gradient) flow. While K efu depends on the hydraulic conductivity K statistical moments (K G , the geometric mean and σ Y 2 , the logconductivity variance for a given structure model) it is not a local property for strongly nonuniform well flow. In hydraulic tomography, flow to a short well segment which is considered here, can be approximated by that toward a source in a three-dimensional heterogeneous medium. We investigate the equivalent conductivity K eq , which unlike K efu , is a measurable property, directly related to the head H (or its mean) and the discharge. The dependence of K eq on K statistics and distance r from the source is determined for steady source flow in a formation of random stationary and statistically isotropic K. For the first time K eq is derived by numerical simulations for σ Y 2 ≤ 4. The derivation of the transition of K eq from the local exact value for r / I << 1 (I is the integral scale) to K efu for r / I ⪆1 is one of the main results of the study. Application to a synthetic example of hydraulic tomography to derive K G , σ Y 2 and I , from H measurements, demonstrates the promise of the method and calls for its further development and validation. [ABSTRACT FROM AUTHOR]

• Flume experiments and numerical models of solute transport in a layered bedform were conducted. • A combination of mechanisms produces irregular spatial patterns of solute transport in the hyporheic zone. • Disconnected transportation was formed under a discontinuous low-permeability layer. Hyporheic exchange is affected by bedform geometry, which induces complex flow paths within the bedform. Additional factors that influence flow and solute transport in the hyporheic zone are layered profile sediments and density-driven flow. This study explored the combined effects of these factors on hyporheic exchange through laboratory experiments and numerical simulations involving infiltrating solute displacing less-dense resident water in a layered bedform with a low permeability layer (LPL). The bedform consisted of three horizontal layers, in which the hydraulic conductivity of the middle layer (LPL) was less than that of the top (TL) and bottom layers (BL). The results demonstrated that a previously unexplored combination of mechanisms (density effects and layered bedform) produces irregular spatial patterns of solute transport in the hyporheic zone. For instance, the width of solute plume within the bottom layers becomes narrowed compared with tracer transport. With increasing density contrast between infiltrating solute and resident water, the solute plume becomes much narrower, forming fingers. Numerical modeling further shows that the hydraulic conductivity contrast (HCC) and relative thickness (RT) of the hyporheic zone layers also affect the spatial solute transport patterns. As the hydraulic conductivity contrast or relative thickness increases, the plume becomes much narrower. Horizontal ambient flow (HAF) dominated in the bottom layers, and lateral solute spreading and mixing intensified with a higher hydraulic conductivity contrast and relative thickness. Furthermore, the vertical solute plume was detached by the horizontal ambient flow in the bottom layers with a discontinuous low permeability layer, forming a discontinuous zone of vertical solute transport. Image, graphical abstract [ABSTRACT FROM AUTHOR]

• Standard DNN methods might lead to highly uncertain and inaccurate predictions. • Multiphysics-informed neural network (MPINN) method is proposed. • Multiphisics constraints improve predictive ability of DNNs. • MPINN allows assimilating multiphysics data. • MPINN reduces data requirements for parameter estimation. Data assimilation for parameter and state estimation in subsurface transport problems remains a significant challenge because of the sparsity of measurements, the heterogeneity of porous media, and the high computational cost of forward numerical models. We present a multiphysics-informed deep neural network machine learning method for estimating space-dependent hydraulic conductivity, hydraulic head, and concentration fields from sparse measurements. In this approach, we employ individual deep neural networks (DNNs) to approximate the unknown parameters (e.g., hydraulic conductivity) and states (e.g., hydraulic head and concentration) of a physical system. Next, we jointly train these DNNs by minimizing the loss function that consists of the governing equations residuals in addition to the error with respect to measurement data. We apply this approach to assimilate conductivity, hydraulic head, and concentration measurements for the joint inversion of these parameter and states in a steady-state advection–dispersion problem. We study the accuracy of the proposed data assimilation approach with respect to the data size (i.e., the number of measured variables and the number of measurements of each variable), DNN size, and the complexity of the parameter field. We demonstrate that the physics-informed DNNs are significantly more accurate than the standard data-driven DNNs, especially when the training set consists of sparse data. We also show that the accuracy of parameter estimation increases as more different multiphysics variables are inverted jointly. [ABSTRACT FROM AUTHOR]

• Multiscale characterization of the effective conductivity distributions. • Improved analytical results using the energy dissipation formulation. • Efficient numerical estimator valid at all scales using that formulation. • Compare scale dependence of full pdf, mean and variance with two other estimators. In large scale heterogeneous aquifer simulations, determining the appropriate coarsening scale λ to define an effective hydraulic conductivity K eff is a challenging task, that involves a trade-off between accuracy and cost. Efficiently adjusting the scale λ is then key, in particular for uncertainty quantification. In this paper, we obtain improved analytical results for the variance of K eff , valid at any scale, in the context of energy dissipation formulation. Using this formulation, we then derive an efficient K eff numerical estimator, and compare it with those of the potential-flow average and permeameter formulations in 2D, for lognormal and binary media, over a wide range of λ and of heterogeneity. We analyze the probability density function (pdf), mean, and variance, of these estimators, comparing them with the analytical results. In the lognormal case, the pdf's are rather similar for the three estimators, and remain lognormal at all scales. In the binary case, slow convergence to an asymptotic regime is observed close to the percolation threshold. Image, graphical abstract [ABSTRACT FROM AUTHOR]

• An effective porosity approach was proposed for describing permeability variation during microbially induced calcite precipitation. • The performance of using different porosity-permeability relationships to simulate the permeability change during microbially induced calcite precipitation has been discussed. • The predictive ability of the model has been estimated through simulations of laboratory tests under different treatment strategies. • The effects of the reaction rate on the spatial and temporal evolution of the precipitated calcite were analyzed. Microbially induced calcite precipitation (MICP) provides the potential for developing innovative and environmentally friendly techniques to improve the engineering properties of soil through reduction of permeability in soil and increase of soil stiffness and strength. In the present research work, coupled bio-chemo-hydraulic modelling was developed to enhance the understanding of the coupled processes involved in MICP and to predict the MICP performance in permeability reduction. In the model, an overall kinetically controlled model is adopted to describe the biochemical reactions, where the reaction rate is dependent on the concentration of both bacteria and chemical reactants. Specifically, an effective porosity concept was proposed and implemented in the model to capture the effects of pore throat blockage on permeability reduction caused by precipitated calcite. Correspondingly, the Kozeny-Carman equation was modified to describe the permeability variation during MICP. This model has been applied to simulate two laboratory experiments. The observed changes in chemical components and hydraulic conductivity can be well reproduced in the model. Through comparisons of the effective porosity concept to other commonly used porosity-permeability relationships, it is estimated that effective porosity should be considered to simulate the permeability reduction observed in MICP. Furthermore, a sensitivity analysis of the reaction-related parameters was conducted. The results of the sensitivity analysis indicate that the maximum urease rate has a strong influence on the biochemical hydraulic responses in MICP. [ABSTRACT FROM AUTHOR]

• A convolutional neural network (CNN) is built to estimate macro-dispersivity from hydraulic conductivity field directly. • The validation of the CNN shows that most estimations fall within ± 20% error range for conductivity fields with moderate heterogeneity. • For a given size of training datasets, the performance of estimating drops with an increase in the heterogeneity of conductivity field. • The CNN trained by conductivity fields with high heterogeneity can achieve superior performance in estimations for fields with less heterogeneity. Building a quantitative relation between the spatial heterogeneity of the hydraulic conductivity fields and the macroscale behavior of solute transport is fundamental for groundwater environment problem. In this work, the deep learning technique is explored to build the functional mapping between the hydraulic conductivity field and the longitudinal macro-dispersivity. We examine the capability of the deep neural network in estimating macro-dispersivities of conductivity fields with different variances. The universality of the trained deep neural network is investigated. Comparisons of the neural network results and the reference values (macro-dispersivities from transport simulation) suggest the promising potential of deep learning technique in porous media with moderate heterogeneity. For a given size of training datasets, the deep neural network produces better macro-dispersivity estimation for the conductivity field with smaller variance. The trained neural network by conductivity fields with larger variance has stronger universality for macro-dispersivity estimation. This study demonstrates that deep neural network can be an effective alternative for estimating macroscale behavior of solute transport by directly interpreting hydraulic conductivity fields. [ABSTRACT FROM AUTHOR]

• This paper introduces the concept of Traveling Pilot Points (TRIPS) for solving the categorical inverse problem • The TRIPS methodology combined with subspace methods is used to identify channel structures in a categorical aquifer system while honoring calibration and prior constraints • Tests on a synthetic example suggest that the TRIPS methodology has promise and could be tested on more complicated problems Categorical parameter distributions are common-place in hydrogeological systems consisting of geologic facies/categories with distinct properties, e.g., high-permeability channels embedded in a low-permeability matrix. Parameter estimation is difficult in such systems because the discontinuities in the parameter space hinder the inverse problem. Previous research in this area has been focused on the use of stochastic methods. In this paper, we present a novel approach based on Traveling Pilot points (TRIPS) combined with subspace parameter estimation methods to generate realistic categorical parameter distributions that honor calibration constraints (e.g., - measured water levels). In traditional implementations, aquifer properties (e.g., hydraulic conductivity) are estimated at fixed pilot point locations. In the TRIPS implementation, both the properties associated with the pilot points and their locations are estimated. Tikhonov regularization constraints are incorporated in the parameter estimation process to produce realistic parameter depictions. For a synthetic aquifer system, we solved the categorical inverse problem by combining the TRIPS methodology with two subspace methods: Null Space Monte Carlo (NSMC) and Posterior Covariance (PC). A posterior ensemble developed with the rejection sampling (RS) method is compared against the TRIPS ensembles. The comparisons indicated similarities between the various ensembles and to the reference parameter distribution. Between the two subspace methods, the NSMC method produced an ensemble with more variability than the PC method. These preliminary results suggest that the TRIPS methodology has promise and could be tested on more complicated problems. [ABSTRACT FROM AUTHOR]

• Heterogeneous aquifers are modeled with log-normal random conductivity fields. • Numerical schemes are tested on two-dimensional benchmark flow problems. • Code verification performed by comparison with analytical manufactured solutions. • Validation done through comparisons between MC statistics and first-order results. • Limited computational feasibility for exponentially correlated log-conductivity. This article presents numerical investigations on accuracy and convergence properties of several numerical approaches for simulating steady state flows in heterogeneous aquifers. Finite difference, finite element, discontinuous Galerkin, spectral, and random walk methods are tested on two-dimensional benchmark flow problems. Realizations of log-normal hydraulic conductivity fields are generated by Kraichnan algorithms in closed form as finite sums of random periodic modes, which allow direct code verification by comparisons with manufactured reference solutions. The quality of the methods is assessed for increasing number of random modes and for increasing variance of the log-hydraulic conductivity fields with Gaussian and exponential correlation. Experimental orders of convergence are calculated from successive refinements of the grid. The numerical methods are further validated by comparisons between statistical inferences obtained from Monte Carlo ensembles of numerical solutions and theoretical first-order perturbation results. It is found that while for Gaussian correlation of the log-conductivity field all the methods perform well, in the exponential case their accuracy deteriorates and, for large variance and number of modes, the benchmark problems are practically not solvable with reasonably large computing resources, for all the methods considered in this study. [ABSTRACT FROM AUTHOR]

• Model for effective conductivity of highly heterogeneous random media. • Computational studies for conductivity of random fields. • Effects of spatial variability of flow fields. • Upscaling the effects of percolation on the macro scale conductivity of porous media. We develop a phenomenological model for the effective conductivity of porous media that consist of two distinct materials characterized by different values of hydraulic conductivity. Our focus is two-fold: first, to relate the terms and parameters of macroscale conductivity to the spatial variability of local flow and conductivity observed in computational experiments, and second to develop a reduced order model for effective conductivity based on those observations. Darcys law is assumed to hold at local (mesoscopic) and large (macroscopic) scales. At the mesoscale a composite medium is a configuration of irregular subdomains consisting of two different materials. The resulting heterogeneities force fluid to flow along irregular paths, and this produces spatial variability in the magnitude and the orientation of the Darcy flux. This variability, along with the macroscopic effective conductivity, depends on the proportion of the total volume allocated to each material, the ratio of the two conductivity values, and the spatial connectivity of material subvolumes. Computational experiments indicate that the effects of heterogeneity are most pronounced when the two conductivity values are very different, and when the volume fraction of the more conductive material is near a percolation threshold. The percolation threshold is the critical volume fraction of the more conductive material at which the medium is no longer traversed by connected paths through that material. The percolation threshold determines the existence of three regimes in effective conductivity, two where the effective conductivity obeys power laws and is dominated by one of the materials, and a third intermediate regime that interpolates between the power laws. [ABSTRACT FROM AUTHOR]

• New model for solute transport in deformable medium with dynamic hydraulic conductivity and degree of saturation • The Dynamic hydraulic conductivity model (Model Kp) results in less excess pore pressure build up, faster soil deformation and slightly slower solute transport. • The Dynamic degree of saturation model (Model Srp) does not affect to peak excess pore pressure but would accelerate contaminant migration. • Dynamic K and dynamic Sr are coupled effects and can not be simply added up. To correctly consider dynamic effects for both simultaneously, Model Kp+Srp is applicable. The effects of the pore pressure related dynamic hydraulic conductivity and dynamic degree of saturation on the consolidation-induced solute transport in deformable unsaturated soils are investigated. The storage equation and solute transport equations are revised to account for the dependence of these two soil parameters on pore pressures. Three dynamic models were conducted for dynamic hydraulic conductivity, dynamic degree of saturation, and both. Compared with the conventional model, the simulation results showed that both hydraulic conductivity and degree of saturation increased near the soil surface where pore pressure exceeded air-entry value. Dynamic hydraulic conductivity results in a slightly slower solute transport while dynamic degree of saturation accelerates the migration of contaminants. Including both dynamic effects produced limited differences in solute concentration, while consolidation results were affected significantly. Although dynamic degree of saturation has some influences on consolidation-induced solute transport, the correlation of the dynamic path is less significant. The air-entry value is an important parameter to determine when soil parameters become dynamic, and it affects the soil consolidation process in certain extent. This study provides new methods to incorporate in-homogeneous soil parameters, which are not only spatiotemporal variables, but also dynamic variations with pore pressure. [ABSTRACT FROM AUTHOR]

• This setup can describe design of convergent tracer tests or containment systems. • Second-order characterization of K suffices to capture the main features of the BTC. • FOA solution provides a reasonable BTC prediction even for high heterogeneity. This paper studies the impact of the hydraulic conductivity structure on solute transport in a radially convergent heterogeneous flow field, which is generated by a well pumping a confined aquifer. The solute is injected instantaneously in a fully penetrating well at a known distance from the pumping well and the analysis is performed numerically by considering four conductivity structures sharing the same two-point statistics, but displaying widely different higher-order moments and connectivity features. The overall objective of the present work is to assess the impact of three-dimensional hydraulic conductivity structures on the Breakthrough Curve (BTC) at a pumping well as a function of a few structural parameters: the logconductivity variance the integral scales along the horizontal and vertical directions, the local-scale dispersivity and the distance of the injection well. The numerical simulations are also compared with two analytical solutions from literature based on the First Order Approximation and the Self Consistent Approximation, respectively. The results indicate that the principal features of the BTC are weakly dependent on the hydraulic conductivity structure, similar to the findings of Janković (2017) for transport in uniform flow. Also, the simple First-Order solution provides a good approximation of the BTC in most of the scenarios considered, and thus it can be a useful and effective tool for screening purposes in applications dealing with aquifers' contamination. [ABSTRACT FROM AUTHOR]

• Comparison of the micro-CT images obtained under the spreading and nonspreading conditions allowed us to gain a significantly improved understanding of the governing displacement physics and the matrix-fracture interactions in fractured rocks. • During secondary gas injection under the spreading condition, the oil spread and formed thick and stable layers between brine and gas. These spreading layers were present not only in the matrix pores but also in the fracture conduit. • The oil phase in the matrix established a hydraulic connectivity through the connected and stable spreading layers with the oil in the fracture, thereby allowing the oil in the matrix to drain into the fracture, which resulted in higher oil recovery. • However, infrequent oil layers observed under nonspreading condition were not stable and collapsed easily as capillary pressure was increased during gas injection. This meant that oil had poor connectivity that hindered oil movement from the matrix to the fracture, leading to lower oil recovery compared to that of the spreading system. This paper presents the results of a detailed experimental study performed to examine fluid flow in a water-wet fractured sandstone rock. Using high-resolution X-ray micro computed tomography technique, we systematically investigate the pore-scale displacement mechanisms and governing interactions between the matrix and fracture during gas injection. We perform two sets of flow experiments using brine/Soltrol 170 (spreading oil)/nitrogen and brine/decalin (nonspreading oil)/nitrogen fluid systems to probe the possible beneficial role of spreading phenomena in transferring oil from the matrix to the fracture during gas injection. Gas injection after primary oil drainage was used to generate a wide range of oil saturations and pore fluid arrangements with both fluid systems. Direct visualization of the pore fluid occupancies at different oil saturations in the medium reveals the significant role that the spreading oil layers play in maintaining the hydraulic conductivity of the oil phase between the matrix and the fracture. This mainly takes place at low remaining oil saturations under the spreading condition where layer drainage, in the presence of stable and connected spreading oil layers, is the dominant displacement mechanism. In the case of high remaining oil saturations, it is observed that the oil movement is primarily governed by piston-like displacements as well as the ability of the gas phase to access the pore elements adjacent to the fracture. Under the nonspreading condition, oil cannot maintain its connectivity due to the absence of the spreading layers, leading to higher residual oil saturations in the matrix and lower ultimate oil recovery. [ABSTRACT FROM AUTHOR]