Claus Haslauer

We analyse two datasets of hydraulic conductivity (K) from the MAcroDispersion Experiment (MADE) site, one measured by direct-push injection logging (DPIL) and the other by flowmeter profiling. The analysis is performed using copula techniques which do not rely on the assumption of multivariate Gaussianity and provide a means to characterise differing degrees of spatial dependence in different quantiles of the K distribution. This characterisation provides better insights into the similarities and differences between the two datasets. In addition to the marginal distributions and the traditional two-point geostatistical measures, copula-based bivariate rank correlation and asymmetry measures are analysed and compared. Furthermore, the parameter estimates obtained by likelihood estimation using n-point theoretical models are analysed. This analysis confirms the similarity of the spatial dependence of K between the two datasets in terms of their marginal distributions and bivariate measures, particularly in the vertical direction. We demonstrate clear indications of the existence of non-Gaussian spatial dependence structures of K at this site. We were able to improve the estimation of the K distribution by taking into account either non-Gaussianity or a censoring threshold, which are expected to lead to a more realistic description of processes that are dependent on K.

The spatial distribution of hydraulic properties in the subsurface controls groundwater flow and solute transport. However, many approaches to modeling these distributions do not produce geologically realistic results and/or do not model the anisotropy of hydraulic conductivity caused by bedding structures in sedimentary deposits. We have developed a flexible object‐based package for simulating hydraulic properties in the subsurface—the Hydrogeological Virtual Realities (HyVR) simulation package. This implements a hierarchical modeling framework that takes into account geological rules about stratigraphic bounding surfaces and the geometry of specific sedimentary structures to generate realistic aquifer models, including full hydraulic‐conductivity tensors. The HyVR simulation package can create outputs suitable for standard groundwater modeling tools (e.g., MODFLOW), is written in Python, an open‐source programming language, and is openly available at an online repository. This paper presents an overview of the underlying modeling principles and computational methods, as well as an example simulation based on the Macrodispersion Experiment site in Columbus, Mississippi. Our simulation package can currently simulate porous media that mimic geological conceptual models in fluvial depositional environments, and that include fine‐scale heterogeneity in distributed hydraulic parameter fields. The simulation results allow qualitative geological conceptual models to be converted into digital subsurface models that can be used in quantitative numerical flow‐and‐transport simulations, with the aim of improving our understanding of the influence of geological realism on groundwater flow and solute transport.

Often, the primary variable of interest can be measured only rarely and/or indirectly. Hence, it is desired to use secondary correlated data to improve the estimation of the primary variable. This paper extends the Gaussian copula-based geostatistical approach for interpolation to include both real-valued primary and secondary data. Even more, multiple types of real-valued secondary data and categorical exhaustive secondary information were incorporated and the associated benefits analysed. The performance of the proposed inclusion of secondary data and information was compared to benchmark geostatistical methods (Ordinary Kriging, co-Kriging, external drift Kriging, and copula-based interpolation of the primary variable alone) by cross-validation and split sampling, with varying compositions of primary and secondary variables. The proposed method was tested for several groundwater quality parameters and the results show that a) the joint copula-based method outperforms co-Kriging in all quality measures; (b) physically impossible estimates (i.e., negative concentrations) never occur; (c) copula-based methods outperform Kriging methods in the quantification of uncertainty and result in more realistic spatial patterns of uncertainty. Joint methods lead to better results in cases when the primary variable is available at fewer locations than the secondary variable (i.e., undersampled case) and when the variables are meaningfully related. In copula-based spatial dependence models that are constructed in probability space, detrimental effects of non-normal (here: bi-modal) marginal distributions on the joint distribution and dissimilar measurement scales of the primary and secondary variables are prevented. The proposed method offers an efficient and non-linear alternative to co-Kriging.

This paper presents a methodology to predict the shape of solute breakthrough curves in heterogeneous aquifers at early times and/or under high degrees of heterogeneity, both cases in which the classical macrodispersion theory may not be applicable. The methodology relies on the observation that breakthrough curves in heterogeneous media are generally well described by lognormal distributions, and mean breakthrough times can be predicted analytically. The log‐variance of solute arrival is thus sufficient to completely specify the breakthrough curves, and this is calibrated as a function of aquifer heterogeneity and dimensionless distance from a source plane by means of Monte Carlo analysis and statistical regression. Using the ensemble of simulated groundwater flow and solute transport realizations employed to calibrate the predictive regression, reliability estimates for the prediction are also developed. Additional theoretical contributions include heuristics for the time until an effective macrodispersion coefficient becomes applicable, and also an expression for its magnitude that applies in highly heterogeneous systems. It is seen that the results here represent a way to derive continuous time random walk transition distributions from physical considerations rather than from empirical field calibration.

This paper demonstrates a maximum likelihood (ML)-based approach to derive representative (“best guess”) contaminant concentrations from data with censored values (e.g., less than the detection limit). The method represents an advancement over existing techniques because it is capable of estimating the proportion of measurements that are true zeros and incorporating varying levels of censorship (e.g., sample specific detection limits, changes through time in method detection). The ability of the method to estimate the proportion of true zeros is validated using precipitation data. The stability and flexibility of the method are demonstrated with stochastic simulation, a sensitivity analysis, and unbiasedness analysis including varying numbers of significant digits. A key aspect of this paper is the application of the statistical analysis to real site rock core contaminant concentration data collected within a plume at two locations using high resolution depth-discrete sampling. Comparison of the representative values for concentrations at each location along the plume center-line shows a larger number of true zeros and generally lower concentrations at the downgradient location according to the conceptual site model, leading to improved estimates of attenuation with distance and/or time and associated confidence; this is not achievable using deterministic methods. The practical relevance of the proposed method is that it provides an improved basis for evaluating change (spatial, temporal, or both) in environmental systems.

This paper demonstrates quantitative reasoning to separate the dataset of spatially distributed variables into different entities and subsequently characterize their geostatistical properties, properly. The main contribution of the paper is a statistical based algorithm that matches the manual distinction results. This algorithm is based on measured data and is generally applicable. In this paper, it is successfully applied at two datasets of saturated hydraulic conductivity (K) measured at the Borden (Canada) and the Lauswiesen (Germany) aquifers. The boundary layer was successfully delineated at Borden despite its only mild heterogeneity and only small statistical differences between the divided units. The methods are verified with the more heterogeneous Lauswiesen aquifer K data-set, where a boundary layer has previously been delineated. The effects of the macro- and the microstructure on solute transport behaviour are evaluated using numerical solute tracer experiments. Within the microscale structure, both Gaussian and non-Gaussian models of spatial dependence of K are evaluated. The effects of heterogeneity both on the macro- and the microscale are analysed using numerical tracer experiments based on four scenarios: including or not including the macroscale structures and optimally fitting a Gaussian or a non-Gaussian model for the spatial dependence in the micro-structure. The paper shows that both micro- and macro-scale structures are important, as in each of the four possible geostatistical scenarios solute transport behaviour differs meaningfully.

The spatial variability of hydraulic conductivity is known to have a strong impact on solute spreading and mixing. In most investigations, its local anisotropy has been neglected. Recent studies have shown that spatially varying orientation in sedimentary anisotropy can lead to twisting flow enhancing transverse mixing, but most of these studies used geologically implausible geometries. We use an object‐based approach to generate stacked scour‐pool structures with either isotropic or anisotropic filling which are typically reported in glacial outwash deposits. We analyze how spatially variable isotropic conductivity and variation of internal anisotropy in these features impacts transverse plume deformation and both longitudinal and transverse spreading and mixing. In five test cases, either the scalar values of conductivity or the spatial orientation of its anisotropy is varied between the scour‐pool structures. Based on 100 random configurations, we compare the variability of velocity components, stretching and folding metrics, advective travel‐time distributions, one and two‐particle statistics in advective‐dispersive transport, and the flux‐related dilution indices for steady state advective‐dispersive transport among the five test cases. Variation in the orientation of internal anisotropy causes strong variability in the lateral velocity components, which leads to deformation in transverse directions and enhances transverse mixing, whereas it hardly affects the variability of the longitudinal velocity component and thus longitudinal spreading and mixing. The latter is controlled by the spatial variability in the scalar values of hydraulic conductivity. Our results demonstrate that sedimentary anisotropy is important for transverse mixing, whereas it may be neglected when considering longitudinal spreading and mixing.

Droughts are serious natural hazards, especially in semi-arid regions. They are also difficult to characterize. Various summary metrics representing the dryness level, denoted drought indices, have been developed to quantify droughts. They typically lump meteorological variables and can thus directly be computed from the outputs of regional climate models in climate-change assessments. While it is generally accepted that drought risks in semi-arid climates will increase in the future, quantifying this increase using climate model outputs is a complex process that depends on the choice and the accuracy of the drought indices, among other factors. In this study, we compare seven meteorological drought indices that are commonly used to predict future droughts. Our goal is to assess the reliability of these indices to predict hydrological impacts of droughts under changing climatic conditions at the annual timescale. We simulate the hydrological responses of a small catchment in northern Spain to droughts in present and future climate, using an integrated hydrological model calibrated for different irrigation scenarios. We compute the correlation of meteorological drought indices with the simulated hydrological time series (discharge, groundwater levels, and water deficit) and compare changes in the relationships between hydrological variables and drought indices. While correlation coefficients linked with a specific drought index are similar for all tested land uses and climates, the relationship between drought indices and hydrological variables often differs between present and future climate. Drought indices based solely on precipitation often underestimate the hydrological impacts of future droughts, while drought indices that additionally include potential evapotranspiration sometimes overestimate the drought effects. In this study, the drought indices with the smallest bias were the rainfall anomaly index, the reconnaissance drought index, and the standardized precipitation evapotranspiration index. However, the efficiency of these drought indices depends on the hydrological variable of interest and the irrigation scenario. We conclude that meteorological drought indices are able to identify years with restricted water availability in present and future climate. However, these indices are not capable of estimating the severity of hydrological impacts of droughts in future climate. A well-calibrated hydrological model is necessary in this respect.

Study region The Lerma catchment, a small (7.3 km2) sub-catchment of the Ebro Basin in northern Spain. Study focus The Lerma catchment underwent a monitored transition to irrigated agriculture, using water from outside the catchment, between 2006 and 2008. This transition has successfully been simulated using the partial-differential-equation-based model HydroGeoSphere, simulating coupled evapotranspiration, surface water, and groundwater flow in the catchment. We use the calibrated model to study how irrigation practices influence the response of the Lerma catchment to the climate change projected for northern Spain. We consider four different irrigation scenarios: no irrigation, present irrigation, climate-adapted irrigation with current crops, and adapted irrigation for crops requiring less water. The climate scenarios are based on four regional climate models and two downscaling methods. New hydrological insight The simulated catchment responses to climate change show clear differences between the irrigation scenarios. In future climate, groundwater levels and base flows decrease more when irrigation is present than without irrigation, because groundwater levels and base flow in present climate are already at low levels without irrigation. In contrast, annual peak discharges increase more in non-irrigated cases than in irrigated cases. Irrigation increases water availability and an associated rise in potential evapotranspiration results in higher actual evapotranspiration during summer. In non-irrigated scenarios, by contrast, actual evapotranspiration in summer is controlled by precipitation and thus decreases in future climate.

Hydraulic tomography (HT) is a method for resolving the spatial distribution of hydraulic parameters to some extent, but many details important for solute transport usually remain unresolved. We present a methodology to improve solute transport predictions by combining data from HT with the breakthrough curve (BTC) of a single forced‐gradient tracer test. We estimated the three dimensional (3D) hydraulic‐conductivity field in an alluvial aquifer by inverting tomographic pumping tests performed at the Hydrogeological Research Site Lauswiesen close to Tübingen, Germany, using a regularized pilot‐point method. We compared the estimated parameter field to available profiles of hydraulic‐conductivity variations from direct‐push injection logging (DPIL), and validated the hydraulic‐conductivity field with hydraulic‐head measurements of tests not used in the inversion. After validation, spatially uniform parameters for dual‐domain transport were estimated by fitting tracer data collected during a forced‐gradient tracer test. The dual‐domain assumption was used to parameterize effects of the unresolved heterogeneity of the aquifer and deemed necessary to fit the shape of the BTC using reasonable parameter values. The estimated hydraulic‐conductivity field and transport parameters were subsequently used to successfully predict a second independent tracer test. Our work provides an efficient and practical approach to predict solute transport in heterogeneous aquifers without performing elaborate field tracer tests with a tomographic layout.

Partial-differential-equation based integrated hydrological models are now regularly used at catchment scale. They rely on the shallow water equations for surface flow and on the Richards’ equations for subsurface flow, allowing a spatially explicit representation of properties and states. However, these models usually come at high computational costs, which limit their accessibility to state-of-the-art methods of parameter estimation and uncertainty quantification, because these methods require a large number of model evaluations. In this study, we present an efficient model calibration strategy, based on a hierarchy of grid resolutions, each of them resolving the same zonation of subsurface and land-surface units. We first analyze which model outputs show the highest similarities between the original model and two differently coarsened grids. Then we calibrate the coarser models by comparing these similar outputs to the measurements. We finish the calibration using the fully resolved model, taking the result of the preliminary calibration as starting point. We apply the proposed approach to the well monitored Lerma catchment in North-East Spain, using the model HydroGeoSphere. The original model grid with 80,000 finite elements was complemented with two other model variants with approximately 16,000 and 10,000 elements, respectively. Comparing the model results for these different grids, we observe differences in peak discharge, evapotranspiration, and near-surface saturation. Hydraulic heads and low flow, however, are very similar for all tested parameter sets, which allows the use of these variables to calibrate our model. The calibration results are satisfactory and the duration of the calibration has been greatly decreased by using different model grid resolutions.

A transformation of a Gaussian random field results in a change of the spatial dependence. In simple terms, the range and the shape of a variogram are different between an underlying Gaussian and a transformed Gaussian field. There is a relationship between the underlying and the transformed spatial dependence. This paper shows how this relationship can be expressed analytically. Rather than generating multidimensional random fields and conducting spatial statistical analysis, we develop an accurate and efficient approach based on a unique mapping of the correlation coefficients of the original multi‐Gaussian fields to the transformed correlation coefficients to evaluate the spatial correlation of transformed non‐Gaussian random fields for any type of geostatistical parameterization. Such a mapping can also yield accurate estimation of the spatial correlation of the underlying Gaussian field given the spatial correlation of the transformed field. Results indicate that (1) a nonlinear transformation of spatially distributed fields usually changes the spatial dependence structure; and (2) the relationship between the dependence structures of the underlying and the transformed field can be expressed analytically, and it is sufficient to do this in one dimension. We use the developed approach to investigate the change of correlations of connected random fields generated by the absolute‐value transformation. Results show (1) the correlation lengths of the underlying Gaussian fields are 1.67 and 2.64 times of those of transformed non‐Gaussian fields for Gaussian and exponential covariance models, respectively; (2) the anisotropic ratio does not change; and (3) the anticorrelation in hole‐effect correlation models disappear.