Search

• A new percolation-based model for unsaturated hydraulic conductivity is developed. • Weibull distribution is used in combination with critical path analysis and percolation theory. • Model validation using measured and simulated data confirms good performance. • Capillary pressure curve slightly outperforms pore-throat size in estimating K (S w). Accurate estimation of unsaturated hydraulic conductivity, K (S w), is crucial for modeling water and solute transport in variably-saturated soils. In this study, a new model for the estimation of K (S w) based on percolation theory and critical path analysis is presented, in which the pore-throat size distribution is described by the three-parameter Weibull probability density function. To evaluate the performance of the proposed K (S w) model, two datasets including more than three hundred samples are analyzed. The first dataset includes 240 pore-network simulations, while the second dataset consists of 101 soil samples from the UNSODA database. For the pore-network simulations, for which both pore-throat size distributions and capillary pressure curves are available, we estimate K (S w) from the capillary pressure curve with RMSLE = 0.54 slightly better than the estimations obtained from the pore-throat size distribution with RMSLE = 0.6. For the soil samples, since only the capillary pressure curve is measured, we estimate K (S w) from the capillary pressure curve and evaluate three previously proposed approaches for determining the exponent, α, in Poiseuille's law. We find that RMSLE = 0.98 for α = 3, RMSLE = 0.87 for α = 2(4 − D) − (3 − D)/(2 D − 3), and RMSLE = 1.06 for α = 4 − 0.74 D , in which D is the fractal dimension characterizing the pore space. Comparing the results with previous studies in which the power-law probability density function was used to characterize capillary pressure curve and pore sizes shows that the proposed three-parameter Weibull distribution, in conjunction with the CPA analysis, may more accurately estimate K (S w) for a wide variety of soils ranging from very fine to very coarse textures. [ABSTRACT FROM AUTHOR]

The incorporation of hard data in geostatistical modeling is crucial for enhancing the accuracy of interpolating or stochastically estimating subsurface spatial features. The hard data at specified points in the model domain serve as a guide in optimizing the unknown parameters to follow the patterns of the hard data. Recently, a novel approach to solving hydrogeologic/reservoir modeling problems has emerged by using deep generative models, specifically generative adversarial networks (GANs), to generate realistic and diverse images of channelized aquifers. This subsequently can be coupled with inverse models to solve parameter estimation problems. This study focused on using an improved GAN, called a progressive growing generative adversarial network (PGGAN), conditioned with hard data to perform parameter estimation of complex facies models by coupling an ensemble smoother with multiple data assimilation (ES-MDA). First, the PGGAN was trained to an image with 128 × 128 resolution. The trained PGGAN was used to generate hydraulic conductivity fields when fed an ensemble of latent variables and hard data. The ES-MDA then was used to update the latent variable with the help of hydraulic head data obtained from the groundwater model. The approach was tested on synthetic hydraulic conductivity data. Results show that this approach was able to perform efficient estimation of an unknown facies model domain. Additionally, the proposed method was applied to a different test case of a facies model exhibiting different statistical characteristics. The results were satisfactory, demonstrating that the method is not constrained to the particular hydraulic conductivity fields introduced in the generator's training. • Coupling deep learning with data assimilation for groundwater modeling. • Both hard data and soft data are conditioned. • Synthetic examples are used to demonstrate the effectiveness. [ABSTRACT FROM AUTHOR]

• A high-speed camera observed gravity-driven unstable finger formation. • At the fingertip, water advances intermittently one pore at a time. • Pore invasions occur at high velocities. • High velocities increase dynamic contact angles. • Greater contact angles reduce the matric potential at the front and cause instability. The imbibition of fluids in porous media has been studied widely. Still, processes of preferential flow under gravity due to instability at the wetting front, crucial in groundwater contamination, have yet to be fully understood. Recent theories using dynamic contact angles could describe unstable flow phenomena but have not been proven experimentally. Therefore, infiltration experiments in small sand-filled chambers were conducted to explore the effect of dynamic contact angles. A high-speed camera recorded pore invasion at the unstable imbibition liquid front. A tensiometer recorded the matric potential. The results show that water moved in milliseconds through a small pore at the wetting front, followed by a stationary period. The maximum observed pore water velocity was 25 cm/s, exceeding the saturated hydraulic conductivity by three orders of magnitude. The relationship between dynamic contact angle and water velocity found by Hoffman in glass tubes could describe which was observed in soils. [Display omitted] [ABSTRACT FROM AUTHOR]

• The variability of solute transport in heterogeneous confined aquifers of variable thickness is quantified by the variance of solute displacement. • Variability of solute displacement is caused by variability of the hydraulic conductivity and the thickness of the aquifer. • Even if the random fields of hydraulic conductivity and aquifer thickness are second-order stationary, the hydraulic head fields could be nonstationary. • The nonstationary hydraulic head fields lead to nonstationary flow fields and thus to nonstationary solute displacement fields. In this work, the variability of regional-scale transport of inert solutes in heterogeneous confined aquifers of variable thickness is quantified by the variance of the displacement of a solute particle. In the traditional stochastic approach to solute transport at the regional scale, the variability of solute displacement is attributed to the variability of the transmissivity fields. In this work, however, variability in solute displacement is attributed to variability in hydraulic conductivity and aquifer thickness fields. A general stochastic methodology for deriving the variance of the displacement of a solute particle based on the convection velocity of solute particles, developed from the relationship between the two-dimensional depth-averaged solute mass conservation equation and the Fokker-Planck equation, is given. Assuming that the fluctuations in log hydraulic conductivity and log thickness of the confined aquifer are second-order stationary processes, a closed-form expression for the solute displacement variance in the mean flow direction is obtained for the case of advection-dominated solute transport. The results show that variation in hydraulic conductivity and aquifer thickness can lead to nonstationarity in the covariance of flow velocity, making longitudinal macrodispersion anomalous and increasing linearly with travel time at large distances. [ABSTRACT FROM AUTHOR]

• An analytical model is developed for inter-aquifer, abandoned-well leakage. • A well injected to the confined aquifer is considered in the leakage model. • A thinner confined aquifer will lead to a higher abandoned well leakage rate. A non-plugged abandoned well may cause unwanted flows and solute transport between aquifers. This research provides a novel analytical model to calculate the leakage rate in an abandoned well between a confined aquifer and an unconfined aquifer under steady-state flow conditions. The abandoned well is treated as a passive preferential flow channel connecting a confined aquifer receiving injection from a separate well and an overlying unconfined aquifer; an impermeable aquiclude separates the two aquifers. A closed-form analytical solution for the steady-state leakage rate through an abandoned well is derived. The concerned abandoned well in this study fully penetrates the receiving unconfined aquifer but only penetrates the ceiling of the injecting confined aquifer with a very small screen length (with respect to the confined aquifer thickness). The injection well is sealed from the receiving unconfined aquifer and fully penetrating in the injecting confined aquifer. The solution developed in this study is compared with a previous solution by Avci (1994) (Avci, 1994. Evaluation of flow leakage through abandoned wells and boreholes. Water Resources Research, 30 (9), 2565–2578). Compared to Avci (1994), this research has two different assumptions: (1) The Avci solution assumed that the leakage effect occurs between two confined aquifers, and leakage concerned in this research occurs between an injecting confined aquifer and a receiving unconfined aquifer. (2) The abandoned well in the Avci (1994) solution fully penetrates both the receiving and injecting confined aquifers, while the abandoned well in this research fully penetrates the receiving unconfined aquifer but only penetrates the ceiling of the injecting confined aquifer with a minor screen length. In both this study and the Avci (1994), the injecting well is isolated from the receiving aquifer and fully penetrating in the injecting aquifer. The results show that the two solutions have similarities in thin confined aquifers but are much different in relatively thick confined aquifers (with thickness greater than 10 m). Using the Sobol method for global sensitivity analysis, we tested controls of the hydraulic conductivity ratio between confined and unconfined aquifers, the transmissivity of the confined aquifer over the injection rate per unit screen length, and the aquifer thickness on model performance. The results show that the model is most sensitive to the transmissivity of the confined aquifer over the injection rate per unit screen length. This research can be applied in multiple applications, such as modeling the leakage effect among two aquifers involving stormwater recharge and using a monitoring well to detect the existence of abandoned wells when the abandoned wells are not visible on the surface or not documented, among others. [ABSTRACT FROM AUTHOR]

• Spatially variable, scale-dependent dispersivities are jointly estimated with hydraulic conductivity. • Joint inversion of tracer and ERT data improves aquifer characterization and transport predictions. • Spatially variable dispersivities outperform uniform dispersivities in transport prediction when adding ERT data. The hydraulic conductivity (K) has been recognized as one of the controlling parameters that significantly impact the transport behavior in groundwater systems. To improve the predictive ability of a transport model, concentration measurements are commonly used in inverse modeling to characterize the K distribution. In addition to K , the spread of the solute plume also depends on dispersion coefficients. The evolution of the plume morphology can only be satisfactorily reproduced with appropriate dispersivities. Dispersivity values or distributions rely on the resolution of the inferred K distribution. Dispersivities were typically treated as spatially uniform in inverse modeling and jointly estimated with K. Although employing spatially variable dispersivities may hold the potential to better capture the solute transport behavior, it may significantly increase the number of unknown parameters for inverse modeling. Solving such an ill-posed inversion problem is challenging, especially when the amount of intrusive borehole-based measurements is limited. As a non-intrusive, cost-effective, and high sampling density method, time-lapse geophysical technique (e.g. , electrical resistivity tomography, ERT) has not yet drawn much attention as a useful source of information for dispersivity estimation. In this study, we developed a joint hydrogeophysical inversion framework to simultaneously estimate spatially variable K and dispersivities by integrating tracer-concentration and ERT-derived potential data. We evaluated the proposed framework by conducting numerical tracer experiments in a highly heterogeneous aquifer. Results show that using limited tracer concentration data alone failed to resolve the fine structure in K field. In this case, there was little difference in the prediction performance of solute transport by considering spatially uniform or variable dispersivities. In contrast, the K distribution was estimated at a significantly higher resolution when integrating tracer-concentration and ERT data, contributing to a remarkable improvement in transport modeling. Powered by the sufficient amount of ERT measurements, jointly estimating spatially variable dispersivities further enhanced the reproduction of tracer plumes and breakthrough curves. [ABSTRACT FROM AUTHOR]

• Preferential pathways at the Darcy scale can be characterized by their bifurcations, and USR formation formed by their tortuous path. • The bifurcation fraction, USR, tortuosity, and fractal dimension scale with the conductivity heterogeneity. • The same power law exponent captures the scaling of bifurcation fraction, USR, tortuosity, and fractal dimension. Transport in porous media on the Darcy scale can be both Fickian and non-Fickian, an outcome dependent on the degree of homogeneity of the hydraulic conductivity pattern, as well as the boundary conditions and flow rate. The non-Fickian manifestation is generally associated with heterogeneous media, which promotes the formation of preferential pathways that funnel the transport. Yet this funneling occurs in both weakly and strongly heterogeneous domains, through a mechanism that is yet to be fully characterized. We model Darcy-scale transport using a particle tracking code (PT) that samples a 2D, lognormally distributed conductivity field with a variance representing a range of homogeneous to heterogeneous domains. We find that the resulting preferential pathways tend to split into more pathways (bifurcations), leaving regions into which particles do not invade, which we refer to as "under sampled regions" (USR), while forming a tortuous path. The fraction of bifurcations decreases downstream, reaching an asymptotic value, with a trend that can be fitted as a power-law of the variance. We show that the same power-law exponent relating the bifurcations to the variance holds true for the USR fraction, tortuosity, and fractal dimension with the same variance. An extension of our work is also presented for varying correlation length of the conductivity spatial distribution. We further expand our analysis to a case of impermeable fraction in a uniform conductivity field and show that the power-law fit still holds. Key Points Preferential pathways at the Darcy scale can be characterized by their bifurcations, and USR formation formed by their tortuous path. The bifurcation fraction, USR, tortuosity, and fractal dimension scale with the conductivity heterogeneity. The same power law exponent captures the scaling of bifurcation fraction, USR, tortuosity, and fractal dimension. [ABSTRACT FROM AUTHOR]

• We solve low-velocity non-Darcian flow to a confined and partially penetrating well in a leaky aquifer system. • An analysis is conducted to examine the temporal and spatial changes in confined groundwater head differences between Darcian and non-Darcian flow conditions. • Non-Darcian flow in the aquitard needs to be considered under small confined aquifer hydraulic conductivity and large threshold pressure gradient of the aquitard. Previous studies on the flow dynamics of leaky aquifer systems have primarily relied on Darcy's law, assuming fully penetrating wells in confined aquifers. However, when dealing with aquitards dominated by clay, the flow behavior often deviates from Darcy's law. Moreover, in the majority of cases, this entails partially penetrating well pumping. This paper introduces a mathematical model for non-Darcian flow in a confined and partially penetrating well system within a leaky aquifer system. The model is based on the low-velocity non-Darcian flow equation, incorporating the threshold pressure gradient. In the confined aquifer, the flow exhibits two-dimensional Darcian flow behavior, while in the aquitard, the flow is characterized by one-dimensional non-Darcian flow in the vertical direction. In the shallow aquifer, the flow is one-dimensional Darcian flow in the radial direction. The mathematical model is solved using the finite difference method, and the obtained results are compared with calculations based on traditional Darcian flow theory. The findings indicate that the confined groundwater head difference, as predicted by the Darcian and non-Darcian flow theories, diminishes radially as the distance from the pumping well increases and vertically as one moves away from the top of the confined aquifer. As pumping progresses, the confined groundwater head difference initially rises and subsequently stabilizes gradually. Moreover, the magnitude of the confined groundwater head difference is more pronounced when the threshold pressure gradient, vertical hydraulic conductivity of the aquitard, pumping rate, hydraulic conductivity and specific yield of shallow aquifer are larger, or when the hydraulic conductivity and specific storage of the confined aquifer are smaller. The effect of well screen length on the confined groundwater head difference is minimal. [ABSTRACT FROM AUTHOR]

• The immiscible fluid flow with zero-velocity and continuity interface are simulated at pore scale. • A modified scaling model of non-wetting fluid conductance with zero-velocity interface is given. • The viscous coupling effect with fluid viscosity, wettability and wetting-film length is addressed. • An accurate scaling model of hydraulic conductance with viscous coupling effect is proposed. Viscous coupling effect on pore-scale flow capacity and its impact on macroscopic flow properties in porous media remain subjects of ongoing debate. This study employs the Stokes equation in conjunction with the zero-velocity interface and continuity interface to investigate viscous coupling effect during two-phase flow through a square capillary tube. Our findings substantiate the significance of the viscous coupling effect in angular pore multiphase flows. Enhanced fluid conductance is linked to fluid viscosity ratios, contact angles, and wetting-film lengths at the pore scale. Proposing a novel upscaling approach, we formulate hydraulic conductance in a square capillary tube for multiphase flows by incorporating a viscous coupling term, and its validation is accomplished through comparison with results from lattice Boltzmann method simulations. Our scaling model predicts hydraulic conductance with a mean relative error of less than 1%, outperforming prior viscous coupling models with errors reaching up to 19%. The derived scaling model, incorporating the viscous coupling effect, holds potential for integration within pore-network models, offering an efficient and precise simulation method for characterizing two-phase flow through porous media at a representative scale. [ABSTRACT FROM AUTHOR]

• Physics Informed Neural Networks (PINNs) is a promising technique. • Current implementation of PINNs face convergence issues in heterogeneous domains. • A new implementation is proposed based on the mixed formulation of the equations. • The new implementation simulate flow, regardless the level and type of heterogeneity. • It provide high accurate results without additional learning time. Current implementations of Physics Informed Neural Networks (PINNs) can experience convergence problems in simulating fluid flow in porous media with highly heterogeneous domains. The objective of this work is to develop an accurate implementation of PINNs to model fluid flow in heterogeneous porous media that avoids alternative approaches in the literature that tend to impose unphysical continuity of the hydraulic conductivity field. The proposed implementation is inspired by the mixed formulation of the governing equations, which separates the mass continuity equation and Darcy's law. This methodology has been used with the mixed finite element method and is known to provide accurate pressure and velocity fields in highly heterogeneous domains. In this work, a similar methodology is applied to PINNs, where the training loss function is based on the decoupled continuity equation and Darcy's law. The separation of the continuity equation and Darcy's law allows for calculating the residual term in the learning loss function without evaluating the spatial derivatives of the discontinuous hydraulic conductivity and associated non-differentiable functions. This approach provides more accurate automatic differentiation of the neural networks for pressure and velocity fields. The new implementation of PINNs can be applied to simulate flow in porous media, regardless of the type and level of heterogeneity with a discontinuous hydraulic conductivity distribution. A variety of structures of the new implementation are investigated. Numerical experiments show that the structure based on different neural networks for the pressure and each component of the velocity field provides the highest accuracy with equivalent training time as other structures. The proposed methodology presents a novel approach to broaden the scope of PINNs in modeling fluid flow in porous media, while preserving the precise representation of the domain's discontinuous features. [ABSTRACT FROM AUTHOR]

Contaminant source and aquifer characterization (CSAC) is critical in groundwater pollution evaluation and remediation. The ensemble smoother with multiple data assimilation (ES-MDA) is utilized to jointly identify contaminant source information and hydraulic conductivities by assimilating time-lapse electrical resistivity tomography (ERT) data. In a synthetic profile with a time-varying release history in a heterogeneous aquifer, we verify the performance of the proposed data assimilation framework. The results show that the CSAC problem could be handled by the proposed approach. The time-varying release history and the high permeability area can be identified with adequate time-lapse ERT measurements. Further reproduction of the evolution of the plume after CSAC also shows consistency with the reference plume. The poor conditioning inversion caused by the filter inbreeding is analyzed by comparing four scenarios with different apparent resistivity measurements. Furthermore, we also evaluate the impact of uncertainties in the petrophysical properties and geophysical observations on our data assimilation framework. The results show that the proposed ES-MDA data assimilation framework could provide a convincing inversion of the time-varying release history and hydraulic conductivities. • We propose a data assimilation framework including the coupling of the simulators of MODFLOW, MT3DS and ResIPy into the Ensemble Smoother with Multiple Data Assimilation (ES-MDA) implementation. • The time-lapse electrical resistivity tomography (ERT) measurements are first employed for the joint identification of contaminant source information and hydraulic conductivities. • We analyze the influence of different bipole–bipole electrode array configurations on the ES-MDA data assimilation framework. [ABSTRACT FROM AUTHOR]

• Transient tidal overheight represented by phase-averaged tidal boundary condition. • Planar beach model successfully predicts overheight for complex coastal morphology. • Boundary condition performance is sensitive to coastal morphology representation. • Phase-resolved simulations and observed data can be replicated reasonably well. • Considerable reduction of coastal groundwater flow simulation runtimes achieved. Ocean tides cause water table overheight near the seaward boundary of coastal aquifers which can have a large impact on groundwater flows and saltwater intrusion. Despite this, regional-scale coastal groundwater flow models often neglect the effects of ocean tides or alternatively assume that the influences of ocean tides and sea level are constant over time. These simplifications are often required because simulation of the phase-resolved tidal signal at a coastal boundary is computationally demanding. Our objective was to derive a phase-averaged tidal boundary condition for application along gently sloping coastlines that enables simulation of real, complex tidal signals combined with sea-level changes due to meteorological effects. The boundary condition extends an existing analytical solution of tidal overheight by including an empirical correction function for conditions where the assumptions of the existing solution are violated. The correction function was developed by conducting a parametric study on an idealized two-dimensional cross-sectional coastal aquifer model and considering the effects of parameters including horizontal hydraulic conductivity, vertical anisotropy, specific yield, aquifer thickness, and beach slope. The performance of the new phase-averaged tidal boundary condition was assessed using a three-dimensional groundwater flow model of the island of Spiekeroog, Northwest Germany, which comprises complex coastal morphology with non-planar beach slopes. Simulations applying the new boundary condition were compared to simulations that adopted a phase-resolved tidal boundary condition and to observed groundwater levels. Accounting for time-varying changes in the tidal signal and sea level was found to be critical to simulate observed data and to adequately reproduce the transient tidal overheight simulated by the phase-resolved model. The performance of the new phase-averaged boundary condition in the regional-scale model varied depending on how the coastal morphology was represented in the boundary condition with local morphological features increasingly important for lower intertidal vertical infiltration capacity (low isotropic hydraulic conductivity or high vertical anisotropy) or specific yield. In conclusion, the new boundary condition presented overcomes current limitations in simulating transient tidal overheight in regional-scale groundwater models. [ABSTRACT FROM AUTHOR]

The connectivity of the more conductive hydrofacies strongly determines flow and transport in heterogeneous media. Here we study solute transport in 3D binary isotropic samples with a proportion p of a high hydraulic conductivity facies (k +), and (1 − p) of a low (k −) one. The k + facies is characterized by two connectivity parameters: a connectivity structure type (no, low, intermediate and high), that controls how well the k + facies is connected, and an integral scale I b , that controls the heterogeneity characteristic lengthscale. Under ergodic conditions, and in the asymptotic Fickian regime that arises only very far from the injection plane, we analyze two transport quantities: the normalized mean solute arrival time 〈 t a ∗ 〉 , and the longitudinal dispersivity α L. As p reaches the percolation threshold p c (p c depends on the connectivity parameters), k + channels spanning the sample along the mean flow direction appear, giving rise to fast flow pathways. A sharp decrease of 〈 t a ∗ 〉 , and a sharp increase of α L , occur when p → p c . As p exceeds p c , a subsequent minimum of 〈 t a ∗ 〉 and a maximum of α L are observed. This result is in contrast with previous ones by other authors that found a maximum of α L at p = p c . On the other hand, p kept fixed, α L decreases as the connectivity of the k + facies increases. We conclude that the connectivity features sampled by the solute particles during their trajectories are retained in the transport quantities even after the asymptotic regime is attained. Also, that connectivity mainly affects α L through a shift or displacement of p c. Finally, the existence of a spatial connectivity structure may imply early, but also late, arrival times, compared with the absence of structure. • Solute transport in percolative binary media is studied, with focus on connectivity. • Media have a proportion p of a high conductivity facies and (1-p) of a low one. • The connectivity sampled by the solute particles is kept in the transport quantities. • Influence of connectivity is accounted for by a shift of the percolation threshold p c. • Solute dispersivity is not maximal at p = p c , but for p > p c . [ABSTRACT FROM AUTHOR]

• First to use superparamagnetic DNA tagged colloids (SiDNAmag) for hydraulic parameter estimation of a 3D sand aquifer, using monte carlo modelling approach at different ionic strengths. • Hydraulic parameter distributions estimated from sidnamag breakthrough curves through Monte Carlo simulation approach were statistically similar to salt tracer. • SiDNAmags are good candidate for determining hydraulic parameters of saturated 3D sand aquifer at short travel distance (50 cm). We investigated the applicability of Silica encapsulated, superparamagnetic DNA particles (SiDNAmag) in determining aquifer hydraulic parameters at different ionic strengths (1 mM, 5 mM, and 20 mM phosphate buffer) of injection suspension. Thereto, in a homogeneous, unconsolidated sand tank we pulse - injected two uniquely sequenced SiDNAmag at two injection points. At 0.5 m and 0.8 m downstream from the injection points, we measured the concentration of SiDNAmags at three vertically distributed and two horizontally distributed sampling locations. We estimated the hydraulic parameter distributions from the SiDNAmag breakthrough curves through a Monte – Carlo approach and compared the parameter distributions with salt tracer breakthrough curves. Our results indicated that at all the ionic strengths, the times of peak concentrations, and the shapes of the breakthrough curves were similar to the salt tracer. As compared to the salt, a 1 – 3 log units reduction in the maximum effluent concentration of SiDNAmag was due to kinetic attachment. The attachment rate reduced from 1 mM to 5 mM phosphate buffer possibly due to competitive adsorption of phosphate onto the favourable attachment sites. SiDNAmag attachment rate further increased in 20 mM buffer suspension, possibly due to the compression of electric double layer and reduction in energy barrier for attachment. The parameter distributions of hydraulic conductivity (k), effective porosity (n e), longitudinal dispersivity (α L), vertical transverse dispersivity (α TV /α L) and horizontal transverse dispersivity (α TH /α L) estimated from the SiDNAmag and the salt tracer breakthrough curves were statistically similar. Our work contributes to the applicability of colloidal SiDNAmags for determining hydraulic parameters at different ionic strength conditions. [ABSTRACT FROM AUTHOR]

• A multi-fidelity Gaussian process model was implemented to estimate the hydro-geological properties by different sources of data. • The accuracy of the model prediction is dependent on the locations and the distribution of both high- and low-fidelity data, especially when data points are sparse. • A Bayesian experimental design algorithm was coupled with the multi-fidelity Gaussian process model to determine future sampling locations. Enhanced water management systems depend on accurate estimation of subsurface hydraulic properties. However, geologic formations can vary significantly, so information from a single source (e.g., widely spaced boreholes) is insufficient in characterizing subsurface aquifer properties. Therefore, multiple sources of information are needed to complement the hydrogeology understanding of a region. This study presents a numerical framework in which information from different measurement sources is combined to characterize the 3D random field in a multi-fidelity prediction model. Coupled with the model, a Bayesian experimental design was used to determine the best future sampling locations. The Upper Sangamon watershed in east-central Illinois was selected as the case study site, where the multi-fidelity Gaussian process model was used to estimate the hydraulic conductivity in the region of interest. Multi-source observation data were obtained from electrical resistivity and borehole pumping tests. The accuracy of the model prediction is dependent on the locations and the distribution of both high- and low-fidelity data. Furthermore, the multi-fidelity model was compared with the single-fidelity model. The uncertainties and confidence in the measurements and parameter estimates were quantified and used to design future cycles of data collection to further improve the confidence intervals. [ABSTRACT FROM AUTHOR]

• Hydraulic conductivity is a key control of the way precipitation signals are filtered in the subsurface. • Hydrologic controls of the signal filtering depend on the timescale being investigated. • The filtering that occurs in streamflow, unsaturated storage, and saturated storage signals have varying physical controls. Watersheds act as low-pass filters, damping and attenuating climatic signals as water moves through the surface and subsurface. This is a well observed phenomenon; however, the ways in which watershed properties control the nature of this filtering are less well documented, especially with respect to groundwater surface water interactions. Here, we use a physically based groundwater surface water model to simulate idealized hillslope ensembles with varying watershed properties (hillslope slope, hydraulic conductivity, and precipitation magnitude) to quantitively explore the impact of watershed configuration on temporal filtering in both the surface and subsurface. To limit the complexities of this system an idealized titled-v domain is used. Multi-decadal simulations (95 years) are run, and then power spectral densities and transfer functions are used to quantify the temporal dynamics and damping of each simulation. Overall, we show that the degree of filtering and the degree of signal transformation is controlled by the total time spent in the subsurface and the degree of groundwater surface water exchanges. The ratio of precipitation to hydraulic conductivity controls the partitioning between infiltration and runoff. Greater infiltration results in less filtering in the subsurface and more filtering in streamflow. For a given precipitation conductivity ratio, deeper water table depths lead to greater streamflow filtering for periods less than 5 years. For time periods greater than 5 years the streamflow filtering is most strongly related to hydraulic conductivity which controls the baseflow dynamics. The majority of the input signal is filtered in the subsurface for short periods less than one year. For longer time scales, hydraulic conductivity is found to be the primary control of filtering and power shift taking place in the subsurface with larger conductivities correlated to less filtering and less of a signal transformation. Deeper water table depths lead to more signal transformation in saturated storage but are not correlated to filtering in unsaturated storage. This is likely due to counteracting effects of higher conductivity (which decreases filtering) and deeper water table depths (which increase filtering). [ABSTRACT FROM AUTHOR]

• An analytical model is developed to handle horizontal aquifer/aquitard heterogeneity. • A high- K aquitard window tends to dampen tide-induced aquifer head fluctuations. • Aquifer head fluctuations may be dampened by a large-size high- S s aquitard window. • Sandbox experiments validated the effects of a high- K aquitard window. Low-permeability aquitards are important sedimentary units that protect the underlying confined aquifers from contamination by seawater and anthropogenic sources in coastal leaky aquifer systems. Therefore, the integrity of the aquitard is of crucial importance. However, the aquitard integrity may be vulnerable to horizontal heterogeneities in the form of aquitard windows. In this study, an innovative analytical model was developed to describe groundwater flow in a coastal leaky aquifer system with horizontal heterogeneities in hydraulic properties of either the aquitard or the confined aquifer. Consequently, the effects of aquitard windows, which manifest as horizontally local heterogeneities in hydraulic conductivity and/or specific storage, on the tidal wave propagation can be considered in the analysis. Analytical solutions were derived for tidal-head or constant-flux inland boundary conditions. The derivation process avoided previous methodologies that involved the assembly of large-size matrices. The newly developed analytical model can cover multiple existing analytical models that are associated with the assumptions of homogeneity in hydraulic properties of the aquifer and the aquitard, infinite horizontal extent of the coastal leaky aquifer system, impermeability of the aquitard, or ignorance of aquitard storativity. Analytical investigations revealed that high-permeability aquitard windows acting as preferential-flow conduits tended to dampen the tide-induced head propagation within the underlying confined aquifer. Moreover, an aquitard with a higher storativity (but the same permeability) and of considerable size can result in an apparent propagation bias. The newly developed analytical model in this study was further corroborated through sandbox experiments. Both piezometer measurements and parameter estimates indicated the non-negligible effects of high-permeability aquitard windows on the tidal wave propagation. The sandbox investigation also highlighted the importance of geostatistical inverse modeling in characterizing local heterogeneities within the confined aquifer and in detecting the location and size of aquitard windows, when measurement data of tide-induced aquifer head fluctuations were provided. [ABSTRACT FROM AUTHOR]

• Clogging risk assessment at managed aquifer recharge (MAR) sites. • Site characterization techniques to estimate model parameters for physical clogging. • The fines remobilization affects the formation of clogging patterns in the field. • Field spatial factors are included to estimate the evolution in the site's infiltration capacity. Managed aquifer recharge (MAR) techniques are in demand to cope with water scarcity challenges posed by climate change and groundwater overexploitation. One of the long-lasting technical issues associated with MAR systems is physical clogging. The intrusion and deposition of external fines during water recharge reduce the infiltration capacity of the site over time. Operation and maintenance (O&M) costs are experienced directly at the site to restore the original efficiency of infiltration rates. Thus, investors need reliable estimations of the risk of clogging during the planning of the site. As a rule, in MAR design, the main parameter of concern for physical clogging is the total suspended solids (TSS), and most clogging models rely on experiment calibrations in 1D sand columns. However, secondary processes can control the development and spatial distribution of physical clogging in field conditions. The proposed work aims to detect key clogging factors directly in the field and to model these processes for reproducibility at other sites. The fieldwork is conducted at the two-stage infiltration basin in Suvereto (Tuscany, Italy). Spatial factors are included in the analysis (i.e. basin topography) to explain clogging patterns in the field altered by erosion processes. The observed clogging profiles at two sampled locations exhibiting clogging are replicated by a mathematical model. Based on the computation of annual erosion rates in the pond and fines' redistribution, the exceeding fines' contents over depth are validated with an RMSE of 2.53% and 12.53%. The infiltration capacity of the site is estimated to reach a stable value of 90% of the initial infiltration capacity over 20 years, given the Suvereto basin features. The model's parameterisation from field measurements represents a great advantage over existing clogging models due to its transferability to other MAR sites. The assessment of the risk of clogging supported by field characterization and numerical modelling is cost-effective and assists the deduction of O&M schemes for MAR sites. [ABSTRACT FROM AUTHOR]

A deep learning technique called Physics Informed Neural Networks (PINNs) is adapted to study steady groundwater flow in unconfined aquifers. This technique utilizes information from underlying physics represented in the form of partial differential equations (PDEs) alongside data obtained from physical observations. In this work, we consider the Dupuit–Boussinesq equation, which is based on the Dupuit–Forchheimer approximation, as well as a recent, more complete model derived by Di Nucci (2018) as underlying models. We then train PINNs on data obtained from steady-state analytical solutions and laboratory based experiments. Using PINNs, we predict phreatic surface profiles given different input flow conditions and recover estimates for the hydraulic conductivity from the experimental observations. We show that PINNs can eliminate the inherent inability of the Dupuit–Boussinesq equation to predict flows with seepage faces. Moreover, the inclusion of physics information from the Di Nucci and Dupuit–Boussinesq models constrains the solution space and produces better predictions than training on data alone. PINNs based predictions are robust and show a little effect from added noise in the training data. Furthermore, we compare the PINNs solutions obtained via the Di Nucci and Dupuit–Boussinesq flow models to examine the effects of higher order flow terms that are included in the Di Nucci formulation but are neglected by the Dupuit–Boussinesq approximation. Lastly, we discuss the effectiveness of using PINNs for examining groundwater flow. • PINNs predict hydraulic conductivity and water table heights using experimental data. • PINNs eliminate inability of Dupuit–Boussinesq equation when predicting seepage face. • Inclusion of physics improves PINNs predictions compared to plain neural networks. [ABSTRACT FROM AUTHOR]

• Cross-gradients constraint inversion is fused into hydraulic tomography. • The method permits deterministic optimization of structural similarities. • Limited improvements are obtained from the cross-gradients constraint. The inclusion of a priori information is helpful to obtain high resolution parameter identifications. A common practice is to exploit structural information in the search for detailed hydrogeological characterizations of subsurface. In this study, we proposed to incorporate cross-gradients as a constraint in nonlinear inverse problems. The resultant inversion scheme is successfully applied to three synthetic experiments and one laboratory rock block experiment with structural similarities in hydraulic parameters and also compared to estimated results with a traditional method. Our results show similar estimated hydraulic parameters whether the objective function has the cross-gradients operator or not. The improvements are only enhanced when assuming one type of parameter is completely known. The limited improvement is proved by the small contribution from the cross-gradients term in the objective function. Therefore, the cross-gradients constraint is a weak structural constraint. The main reason is that drawdown data has already contained correlated hydrogeological parameter information. Estimated results substantiate the rationalization for a simplified assumption that hydraulic conductivity and specific storage are treated as mutual independent variables. [ABSTRACT FROM AUTHOR]

We consider non-linear Bayesian inversion problems targeting the geostatistical hyperparameters of a random field describing hydrogeological or geophysical properties given hydrogeological or geophysical data. This problem is of particular importance in the non-ergodic setting as there are no analytical upscaling relationships linking the data to the hyperparameters, such as, mean, standard deviation, and integral scales. Full inversion of the hyperparameters and the local properties of the field (typically involving many thousands of unknowns) brings substantial computational challenges, such that simplifying model assumptions (e.g., homogeneity or ergodicity) are typically made. To prevent the errors resulting from such simplified assumptions while also circumventing the burden of high-dimensional full inversions, we use a pseudo-marginal Metropolis–Hastings algorithm that treats the random field as latent variables. In this random effects model, the intractable likelihood of observing the data given the hyperparameters is estimated by Monte Carlo averaging over realizations of the random field. To increase the efficiency of the method, low-variance approximations of the likelihood ratio are obtained by using importance sampling and by correlating the samples used in the proposed and current steps of the Markov chain. We assess the performance of this correlated pseudo-marginal method by considering two representative inversion problems involving diffusion-based and wave-based physics, respectively, in which we infer the hyperparameters of (1) hydraulic conductivity fields using apparent hydraulic conductivity data in a data-poor setting and (2) fracture aperture fields using borehole ground-penetrating radar (GPR) reflection data in a more data-rich setting. For the first test case, we find that the correlated pseudo-marginal method generates similar estimates of the geostatistical mean as classical rejection sampling, while an inversion assuming ergodicity provides biased estimates. For the second test case, we find that the correlated pseudo-marginal method estimates the hyperparameters well, while rejection sampling is computationally unfeasible and a simplified model assuming homogeneity leads to biased estimates. [Display omitted] • Hyperparameters enable general descriptions of target property fields. • No analytical relation linking hyperparameters and data in the non-ergodic setting. • Correlated pseudo-marginal method accounts for the influence of sample realizations. • It avoids a high-dimensional target space and errors from simplified assumptions. [ABSTRACT FROM AUTHOR]

Precise information about the spatial distribution of hydraulic conductivity (K) in an aquifer is essential for the reliable modeling of groundwater flow and transport processes. In this study, we present results of a new inversion procedure for induced polarization (IP) data that incorporates petrophysical relations between electrical and hydraulic parameters, and therefore allows for the direct computation of K. This novel approach was successfully implemented for the Bolstern aquifer analog by performing synthetic IP experiments with a combined surface and cross-borehole setup. From these data, the distribution of K was retrieved with high accuracy and resolution, showing a similar quality compared to images achieved by hydraulic tomography. To further improve the quantitative estimates of K , we use synthetic pumping test data to inform two novel calibration strategies for the IP inversion results. Both calibrations are especially helpful for correcting a possible bias of the IP inversion, e.g., due to resolution limitations and/or to bias in the underlying petrophysical relations. The simulation of tracer experiments on the retrieved tomograms highlights the accuracy of the inversion results, as well as the significant role of the proposed calibrations. • New IP inversion procedure for direct computation of hydraulic conductivity. • IP inversion results have similar quality as hydraulic tomography. • Hydraulic calibration of IP results to improve transport simulations. [ABSTRACT FROM AUTHOR]

• A novel framework IGSM for heterogeneous steady-state parameter inversion problems. • Coupling of group evolutionary algorithm and reduced-order model. • IGSM is fast, simple, reliable and requires no posterior error estimation. • Results of IGSM are less affected by the density and parameter variation range. • IGSM can quickly solve the parameter inversion of large-scale steady state problems. To solve the time-consuming parameter inversion problem based on swarm evolution algorithm, this study proposes an improved greedy sampling method (GSM) (IGSM) based on model reduction technology to rapidly estimate the hydraulic conductivity of steady-state seepage problem. Different from the traditional GSM, IGSM neither needs to estimate the posterior error bound of the reduced-order model (ROM) nor does it require a training sample set to train the ROM. Instead, the observation data are directly applied to the training process of ROM. Each iterative process of IGSM includes a parameter inversion process to obtain the parameters that best match the observation information and use it to update the current ROM. The improved differential evolution algorithm (e.g., MMRDE) is used as the search algorithm, and the resulting algorithm is denoted as IGSM-MMRDE. After comparative testing, IGSM-MMRDE was found to have higher accuracy and lower calculation cost than GSM-MMRDE. In addition, IGSM-MMRDE is also capable of high-dimensional parameter inversion problems and parameter inversion with a wider search range. After applying IGSM-MMRDE to the inversion of the rock mass hydraulic conductivity of the natural seepage field in the dam site area of a hydropower station, IGSM-MMRDE converges quickly and takes much less time than ADINA-MMRDE based on the full-order model (approximately 0.31% of ADINA-MMRDE). The performance results are even better than the ADINA-MMRDE using the fixed evolutionary generation. Therefore, the advantages of using IGSM in the rapid inversion of the initial seepage field of large-scale projects are extremely evident. [ABSTRACT FROM AUTHOR]

For strongly non-linear and high-dimensional inverse problems, Markov chain Monte Carlo (MCMC) methods may fail to properly explore the posterior probability density function (PDF) given a realistic computational budget and are generally poorly amenable to parallelization. Particle methods approximate the posterior PDF using the states and weights of a population of evolving particles and they are very well suited to parallelization. We focus on adaptive sequential Monte Carlo (ASMC), an extension of annealed importance sampling (AIS). In AIS and ASMC, importance sampling is performed over a sequence of intermediate distributions, known as power posteriors, linking the prior to the posterior PDF. The AIS and ASMC algorithms also provide estimates of the evidence (marginal likelihood) as needed for Bayesian model selection, at basically no additional cost. ASMC performs better than AIS as it adaptively tunes the tempering schedule and performs resampling of particles when the variance of the particle weights becomes too large. We consider a challenging synthetic groundwater transport inverse problem with a categorical channelized 2D hydraulic conductivity field defined such that the posterior facies distribution includes two distinct modes. The model proposals are obtained by iteratively re-simulating a fraction of the current model using conditional multiple-point statistics (MPS) simulations. We examine how ASMC explores the posterior PDF and compare with results obtained with parallel tempering (PT), a state-of-the-art MCMC inversion approach that runs multiple interacting chains targeting different power posteriors. For a similar computational budget, ASMC outperforms PT as the ASMC-derived models fit the data better and recover the reference likelihood. Moreover, we show that ASMC partly retrieves both posterior modes, while none of them is recovered by PT. Lastly, we demonstrate how the power posteriors obtained by ASMC can be used to assess the influence of the assumed data errors on the posterior means and variances, as well as on the evidence. We suggest that ASMC can advantageously replace MCMC for solving many challenging inverse problems arising in the field of water resources. • Adaptive sequential Monte Carlo (ASMC) approximates the posterior and the evidence. • ASMC is combined with sequential geostatistical resampling to handle complex priors. • ASMC outperforms parallel tempering considering the same computational budget. • Intermediate ASMC results are interpreted as results for larger assumed data errors. [ABSTRACT FROM AUTHOR]

A semi-analytical solution of the Richards Equation posed on a two-layered one-dimensional soil supplied with various boundary conditions is derived under a constraint that the constitutive relations are exponentially dependent on the pressure head. It allows for a transformation of the Richards Equation into a linear parabolic partial differential equation that governs a spatial–temporal function that represents the hydraulic conductivity. The procedure is proceeded with expressing this function as a linear combination of a set of eigenfunctions associated with a novel Sturm—Liouville problem that reflects the layer system and an auxiliary function that depends only on the spatial variable and the pressure head at the interface at the time of interest. All the relevant coefficients in the representation satisfy a nonlinear differential–algebraic system gathered from imposing the continuity of the pressure head and its flux at the interface. Two different approximations of the derivatives yield algebraic systems to be solved by the Newton method of iteration. Several pertinent numerical experiments demonstrating the approach are discussed and compared with the standard finite volume method. • Derivations of a semi-analytical solution is constructed. • Two approximations are considered to construct the solution of the problem. • Comparison with a standard finite volume indicates a close similarity of the solutions. [ABSTRACT FROM AUTHOR]

Nitrate contamination of aquifers is threatening the sustainability of water resource production around the globe. Denitrification is an ubiquitous process in groundwater systems and needs to be represented in managerial and predictive efforts. Yet, fully understanding the temporal and spatial dynamics of all biogeochemical and physical properties controlling nitrate reduction at the basin scale remains a challenging task. Conceptual simplifications such as homogenization are used, but their implications on model performance are often not well assessed. This paper (1) analyzes the importance of denitrification spatial variability on the fate of non-point source nitrate contamination, and (2) assesses the performance of upscaled denitrification rate coefficient values. Heterogeneity in the hydraulic conductivity of a typical unconsolidated alluvial aquifer system and its associated uncertainty is systematically represented in our 3D numerical simulations of flow and transport. The spatial variability in the nitrate reduction capacities of the aquifer is accounted for by 3 different conceptual models: a positive and a negative correlation with the aquifer hydraulic properties and a uniform random distribution. Our results indicate that heterogeneity in the denitrification rate coefficient – as well as the model describing spatial variability – has a significant impact on concentration statistics observed at a series of extraction wells. However, we also show that, under the modeling setting used in this study (e.g., stationarity of reactivity), upscaled values of denitrification rate coefficients allow a virtually perfect reproduction of nitrate concentration statistics, a promising outcome for employing parameter estimation (calibration) of effective rate coefficients of denitrification: A detailed description of the spatial variability of the nitrate reduction capacity of the alluvial aquifer is not necessary. Yet, these homogeneous, upscaled denitrification rate coefficients result from a complex interaction between the biochemistry and the physics of the aquifer, which makes their direct estimation a challenging task. [Display omitted] • Spatial variability structure of denitrification rate coefficients influence nitrate dynamics. • Upscaled, effective denitrification rate coefficients reproduce well concentration statistics. • Upscaling denitrification rates is complex due to interplays between key processes. [ABSTRACT FROM AUTHOR]

• Darcian seepage into circular and quasi-circular tunnels is modeled analytically and by HYDRUS. • Extremely high pore pressures, hydraulic gradients, and Riesenkampf forces are detected in the vicinity of even grouted tunnels' contours. • The seepage flow rates into the tunnels can be determined either by the sink-source approximation or numerically. • Seepage refraction for grouted tunnels evidences essentially 2D flow nets both in the soil and liner. Darcian, 2-D flows to subterranean holes are studied analytically (by the methods of complex analysis) and numerically (by HYDRUS). For flow towards two circular or quasi-circular tunnels, reconstructed as isobars generated by two sinks under a ponded homogeneous soil surface, the flow nets, the velocity vector fields, and Riesenkampf's seepage force vectors are found. The position of the two tunnels is optimized using a criterion of the total area of the empty space comprised by the isobars with the admissible seepage flow rate as a constraint and the locus of the tunnels as a control variable. The case of a partially-filled tunnel, the contour of which is composed of two conterminous isobaric and equipotential arcs, is also analyzed. For a grouted tunnel, the refraction problem for two potential fields in two subdomains of the half-plane of the seepage flow domain is solved for an arbitrary contrast between the hydraulic conductivities of the liner and ambient soil. The tunnel grouting is, generally, eccentric with respect to the tunnel contour. High hydraulic gradients in the grouting are detected, which is a long-term seepage-harbinger for any liner materials used by geotechnical engineers. [ABSTRACT FROM AUTHOR]

• A low-permeability surface mulch (LPSM) above drains was proposed to improve the leaching efficiency. • The LPSM effectiveness was examined using sand tank experiments and modeling, with an analytical solution used to predict leaching periods. • The LPSM reduced base case leaching periods and freshwater application rates by up to 60% and 80%, respectively. • The LPSM approach outperforms the previously proposed subsurface barrier scheme. Surface ponding is often used in combination with subsurface drainage systems to leach excessive salts from the saline soils of reclaimed coastal areas. However, soils that are not immediately proximal to subsurface drains are typically subjected to reduced leaching efficiency (LE). This study evaluates an approach for increasing the LE through the use of a low-permeability surface mulch (LPSM) above the drain. Sand tank experiments and 2D numerical simulations were conducted to examine the effectiveness of LPSM configurations under continuous surface flooding. The results show that the addition of a LPSM leads to more spatially uniform, accelerated salt leaching towards subsurface drains. As a result, base case leaching periods and freshwater application rates were reduced by up to ∼60% and ∼80%, respectively. The sensitivity analysis based on field-scale models revealed optimal values (i.e., leading to maximum LE) of the LPSM hydraulic conductivity and length. Moreover, an approximate method is applied to estimate the leaching period, and to thereby assist in the design of drainage-LPSM schemes. We conclude that the addition of a LPSM is an effective strategy to enhance salt leaching under surface ponding, adding to existing methods for ameliorating saline soils. [ABSTRACT FROM AUTHOR]

• New analytical solutions for analysing pumping test near an anisotropic fault. • Fault zone draining aquifers with differing properties. • Unconventional application of the well-image theory. • Satisfactory comparison with numerical modelling in all aquifer compartments. • Useful for analysing pumping tests in fractured media. We present a new set of analytical solutions for analysing pumping tests in a well near a vertical fault zone, dyke or vein. We consider a fault zone with finite width, storativity and hydraulic conductivity, which can be anisotropic, and without limitation of the diffusivity contrast between the three aquifer domains (pumped aquifer, fault zone and opposite-side aquifer). This configuration, not described before, also considers flow transience within the fault zone. Drawdown solutions were developed for the three-aquifer domains based on an unconventional application of the well-image theory. Their comparison with numerical modelling results was very satisfactory. Drawdown behaviour at the pumping well is discussed for contrasted properties between the three-aquifer domains. Original analytical solutions for transient flow along both sides of the fault zone, and net transient flow from the fault itself, are proposed and discussed, with particular attention to flow reversal, from the pumped compartment to the fault zone, which occurs where the fault zone is the most transmissive. Finally, we extend the solutions to the case of a well intersecting and pumping a vertical fracture of finite length located near the same vertical geological discontinuity. [ABSTRACT FROM AUTHOR]

• The ability to verify predicted groundwater discharge patterns at scale is limited. • Discharge zones mapped with thermal infrared are used to evaluate predicted patterns. • Predicted and observed discharge patterns agree well in large but not small streams. • Improved modeling approaches are needed to capture headwater discharge processes. Predicting baseflow dynamics, protecting aquatic habitat, and managing legacy contaminants requires explicit characterization and prediction of groundwater discharge patterns throughout river networks. Using handheld thermal infrared (TIR) cameras, we surveyed 47 km of stream length across the Farmington River watershed (1,570 km2; CT and MA, USA), mapping locations of bank and waterline groundwater discharges based on their thermal signature. Using the observed groundwater discharge locations and predicted groundwater discharge rates from 6 variations of a numerical groundwater-flow model (MODFLOW-NWT), we compared 1) predicted groundwater-discharge rates in areas with and without observed groundwater discharge, 2) spatial patterns of observed and predicted groundwater discharge locations, and 3) density of observed groundwater discharge locations with predicted discharge rates. Five of six models reasonably predicted the spatial patterns of discharge locations along the 5th order mainstem, but fewer models predicted groundwater discharge patterns in smaller streams. Our results highlight 1) the feasibility of using TIR observations to evaluate groundwater models, 2) model parameters that influence discharge prediction accuracy (riverbed sediment and bedrock hydraulic conductivity and river-aquifer connections), and 3) current strengths and future opportunities for improved modeling of groundwater-discharge patterns. [ABSTRACT FROM AUTHOR]

The Horizontal Reactive Media Treatment (HRX) well is a novel technology for in situ treatment of contaminated groundwater. The most significant technical performance risk associated with field-scale implementation of the HRX well is the potential for water to bypass the well untreated. Therefore, the groundwater capture simulation plays a very important role in the design of the HRX well. We build a model based on a cell-based smoothed radial point interpolation method (CS-RPIM) for simulating groundwater flow captured by the HRX well in both 3D test pit pilot scale and field scale. This method is compared with the finite element method and measured values. Capture width and groundwater velocity obtained by CS-RPIM are consistent with the measured values. We further develop MATLAB software modules to analyze the effects of the hydraulic conductivities, dimension of well and other parameters on the groundwater capture zone, velocity and residence time in the reactive media. The results show that increase of well diameter has a significant effect on enlargement of the capture zone and the extension of media treatment time. When well diameter increases by 100%, horizontal and vertical capture width increase by 57% and 50% respectively, and media treatment time can be extended by 68%. • A realistic 3D model was built for simulating groundwater captured by the HRX well. • The simulation results outperform FEM or the average capture width formula results. • MATLAB software modules are developed for hydraulic capture zone simulation. • Impact of various factors on the capture zone, groundwater flow velocity and treatment time are analyzed. [ABSTRACT FROM AUTHOR]

• 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]

• 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]

• 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]

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]

• 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]

• 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]

• 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]