Search

Comision Nacional de Investigacion Cientifica y Tecnologica (CONICYT) CONICYT FONDECYT 3170079 3180145 1171032 CONICYT/PIA Project AFB180004

Comision Nacional de Investigacion Cientifica y Tecnologica (CONICYT) CONICYT FONDECYT 3170079 CONICYT/PIA Project AFB180004

Stochastic precipitation generators (SPGs) produce synthetic precipitation data and are frequently used to generate inputs for physical models throughout many scientific disciplines. Especially for large data sets, statistical parameter estimation is difficult due to the high dimensionality of the likelihood function. We propose techniques to estimate SPG parameters for spatiotemporal precipitation occurrence based on an emerging set of methods called Approximate Bayesian computation (ABC), which bypass the evaluation of a likelihood function. Our statistical model employs a thresholded Gaussian process that reduces to a probit regression at single sites. We identify appropriate ABC penalization metrics for our model parameters to produce simulations whose statistical characteristics closely resemble those of the observations. Spell length metrics are appropriate for single sites, while a variogram-based metric is proposed for spatial simulations. We present numerical case studies at sites in Colorado and Iowa where the estimated statistical model adequately reproduces local and domain statistics.

Computer simulation models are continually growing in complexity with increasingly more factors to be identified. Sensitivity Analysis (SA) provides an essential means for understanding the role and importance of these factors in producing model responses. However, conventional approaches to SA suffer from (1) an ambiguous characterization of sensitivity, and (2) poor computational efficiency, particularly as the problem dimension grows. Here, we present a new and general sensitivity analysis framework (called VARS), based on an analogy to “variogram analysis,” that provides an intuitive and comprehensive characterization of sensitivity across the full spectrum of scales in the factor space. We prove, theoretically, that Morris (derivative-based) and Sobol (variance-based) methods and their extensions are special cases of VARS, and that their SA indices can be computed as by-products of the VARS framework. Synthetic functions that resemble actual model response surfaces are used to illustrate the concepts, and show VARS to be as much as two orders of magnitude more computationally efficient than the state-of-the-art Sobol approach. In a companion paper, we propose a practical implementation strategy, and demonstrate the effectiveness, efficiency, and reliability (robustness) of the VARS framework on real-data case studies.

Based on the theoretical framework for sensitivity analysis called “Variogram Analysis of Response Surfaces” (VARS), developed in the companion paper, we develop and implement a practical “star-based” sampling strategy (called STAR-VARS), for the application of VARS to real-world problems. We also develop a bootstrap approach to provide confidence level estimates for the VARS sensitivity metrics and to evaluate the reliability of inferred factor rankings. The effectiveness, efficiency, and robustness of STAR-VARS are demonstrated via two real-data hydrological case studies (a 5-parameter conceptual rainfall-runoff model and a 45-parameter land surface scheme hydrology model), and a comparison with the “derivative-based” Morris and “variance-based” Sobol approaches are provided. Our results show that STAR-VARS provides reliable and stable assessments of “global” sensitivity across the full range of scales in the factor space, while being 1–2 orders of magnitude more efficient than the Morris or Sobol approaches.

We present a Bayesian inversion method for the joint inference of high-dimensional multi-Gaussian hydraulic conductivity fields and associated geostatistical parameters from indirect hydrological data. We combine Gaussian process generation via circulant embedding to decouple the variogram from grid cell specific values, with dimensionality reduction by interpolation to enable Markov chain Monte Carlo (MCMC) simulation. Using the Matern variogram model, this formulation allows inferring the conductivity values simultaneously with the field smoothness (also called Matern shape parameter) and other geostatistical parameters such as the mean, sill, integral scales and anisotropy direction(s) and ratio(s). The proposed dimensionality reduction method systematically honors the underlying variogram and is demonstrated to achieve better performance than the Karhunen-Loeve expansion. We illustrate our inversion approach using synthetic (error corrupted) data from a tracer experiment in a fairly heterogeneous 10,000-dimensional 2-D conductivity field. A 40-times reduction of the size of the parameter space did not prevent the posterior simulations to appropriately fit the measurement data and the posterior parameter distributions to include the true geostatistical parameter values. Overall, the posterior field realizations covered a wide range of geostatistical models, questioning the common practice of assuming a fixed variogram prior to inference of the hydraulic conductivity values. Our method is shown to be more efficient than sequential Gibbs sampling (SGS) for the considered case study, particularly when implemented on a distributed computing cluster. It is also found to outperform the method of anchored distributions (MAD) for the same computational budget.

The fluid flow environment associated with high Glossosoma abundance is predicted by large-eddy simulation of a natural turbulent open-channel flow. The spatial distribution of Glossosoma was depicted by high resolution physical variables described by fluid flow and streambed topography. Variogram analysis of the streambed topography revealed a characteristic length scale of the streambed of the order 0.2 m over which bed roughness height was correlated. Flow simulation output was spatially and temporally averaged over the streambed characteristic length scale and linked to Glossosoma spatial density. A dimensionless scaling relationship between Glossosoma spatial distribution and streamwise velocity averaged in the longitudinal and transverse direction, spatial velocity fluctuation, and spanwise vorticity from the computational fluid dynamics simulation output explained 79% of the variation in observed dimensionless Glossosoma spatial density. The analysis demonstrated that computational fluid mechanics and high resolution bed topography could be instrumental in predicting benthic macroinvertebrate spatial distribution in streams and rivers.

Regional climate models (RCMs) are a standard tool for downscaling climate forecasts to finer spatial scales. The evaluation of RCMs against observational data is an important step in building confidence in the use of RCMs for future projection. In addition to model performance in climatological means and marginal distributions, a model's ability to capture spatiotemporal relationships is important. This study develops two approaches: (1) spatial correlation/variogram for a range of spatial lags, with total monthly precipitation and nonseasonal precipitation components used to assess the spatial variations of precipitation, and (2) spatiotemporal correlation for a wide range of distances, directions, and time lags, with daily precipitation occurrence used to detect the dynamic features of precipitation. These measures of spatial and spatiotemporal dependence are applied to a high-resolution RCM run and to the National Center for Environmental Prediction (NCEP)-U.S. Department of Energy Atmospheric Model Intercomparison Project II reanalysis data (NCEP-R2), which provide initial and lateral boundary conditions for the RCM. The RCM performs significantly better than NCEP-R2 in capturing both the spatial variations of total and nonseasonal precipitation components and the spatiotemporal correlations of daily precipitation occurrences, which are related to dynamic behavior of precipitating systems. The improvements are apparent not only at resolutions finer than that of NCEP-R2 but also when the RCM and observational data are aggregated to the resolution of NCEP-R2.

We develop a stochastic approach to construct channelized 3-D geological models constrained to borehole measurements as well as geological interpretation. The methodology is based on simple 2-D geologist-provided sketches of fluvial depositional elements, which are extruded in the 3rd dimension. Multiple-point geostatistics (MPS) is used to impair horizontal variability to the structures by introducing geometrical transformation parameters. The sketches provided by the geologist are used as elementary training images, whose statistical information is expanded through randomized transformations. We demonstrate the applicability of the approach by applying it to modeling a fluvial valley filling sequence in the Maules Creek catchment, Australia. The facies models are constrained to borehole logs, spatial information borrowed from an analogue and local orientations derived from the present-day stream networks. The connectivity in the 3-D facies models is evaluated using statistical measures and transport simulations. Comparison with a statistically equivalent variogram-based model shows that our approach is more suited for building 3-D facies models that contain structures specific to the channelized environment and which have a significant influence on the transport processes.

Resolution analysis is a crucial appraisal procedure in solving general estimation problems, especially for correctly interpreting the results of spatial analysis schemes. Resolution analyses based on the resolving kernels are typically applied to small inverse problems only when the inverse operators are explicitly accessible. Stochastic simulation schemes have been proposed to extract empirical resolution information for solving large inverse problem. In this study, we generalize the formulation of the empirical resolution length and derive the characteristic length of the point spread function for general estimation methods such as minimum curvature gridding and kriging interpolation schemes that are not equipped with explicitly accessible resolving kernels. The implementation of these resolution diagnostics has not been possible in the past and is demonstrated in this study to facilitate clarifying the advantages and limitations of these widely used methods. In addition, we compare these schemes, based on the resolution appraisal, with a multiscale gridding algorithm in the spatial analysis of the Pacific seafloor heat flow observations. By depicting the pattern of the resolution length variations of both the empirical averaging function and the point spread function for each of the estimated models, we demonstrate that schemes equipped with multiscale capability are more favorable for accommodating sparse, nonuniform data distribution than stationary schemes, such as the kriging method. Furthermore, the empirical resolution pattern constructed in this study facilitates the selection of an appropriate reference function and radii of influence for fitting the variogram, which is difficult but critical when using the kriging method.

[1] Nonstationary random fields such as fractional Brownian motion and fractional Levy motion have been studied extensively in the hydrology literature. On the other hand, random fields that have nonstationary increments have seen little study. A mathematical argument is presented that demonstrates processes with stationary increments are the exception and processes with nonstationary increments are far more abundant. The abundance of nonstationary increment processes has important implications, e.g., in kriging where a translation-invariant variogram implicitly assumes stationarity of the increments. An approach to kriging for processes with nonstationary increments is presented and accompanied by some numerical results.

[1] In mountainous regions, snow accumulation on the ground is crucial for mountain hydrology and water resources. The present study investigates the link between the spatial variability in snowfall and in snow accumulation in the Swiss Alps. A mobile polarimetric X-band radar deployed in the area of Davos (Switzerland) collected valuable and continuous information on small-scale precipitation for the winter seasons of 2009/2010 and 2010/2011. Local measurements of snow accumulation were collected with airborne laser-scanning for the winters of 2007/2008 and 2008/2009. The spatial distribution of snow accumulation exhibits a strong interannual consistency that can be generalized over the winters in the area. This unique configuration makes the comparison of the variability in total snowfall amount estimated from radar and in snow accumulation possible over the diverse winter periods. As expected, the spatial variability, quantified by means of the variogram, is shown to be larger in snow accumulation than in snowfall. However, the variability of snowfall is also significant, especially over the mountain tops, leads to preferential deposition during snowfall and needs further investigation. The higher variability at the ground is mainly caused by snow transport.

[1] 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.

A new approach is presented to condition geostatistical simulation of high-permeability zones in fractured rock to hydraulic-connection data. A simulated-annealing algorithm generates three-dimensional (3-D) realizations conditioned to borehole data, inferred hydraulic connections between packer-isolated borehole intervals, and an indicator (fracture zone or background-K bedrock) variogram model of spatial variability. We apply the method to data from the U.S. Geological Survey Mirror Lake Site in New Hampshire, where connected high-permeability fracture zones exert a strong control on fluid flow at the hundred-meter scale. Single-well hydraulic-packer tests indicate where permeable fracture zones intersect boreholes, and multiple-well pumping tests indicate the degree of hydraulic connection between boreholes. Borehole intervals connected by a fracture zone exhibit similar hydraulic responses, whereas intervals not connected by a fracture zone exhibit different responses. Our approach yields valuable insights into the 3-D geometry of fracture zones at Mirror Lake. Statistical analysis of the realizations yields maps of the probabilities of intersecting specific fracture zones with additional wells. Inverse flow modeling based on the assumption of equivalent porous media is used to estimate hydraulic conductivity and specific storage and to identify those fracture-zone geometries that are consistent with hydraulic test data.

A modification of previously published solutions regarding the spatial variation of hydraulic heads is discussed whereby the semivariogram of increments of head residuals (termed head residual increments HRIs) are related to the variance and integral scale of the transmissivity field. A first-order solution is developed for the case of a transmissivity field which is isotropic and whose second-order behavior can be characterized by an exponential covariance structure. The estimates of the variance σY2 and the integral scale λ of the log transmissivity field are then obtained via fitting a theoretical semivariogram for the HRI to its sample semivariogram. This approach is applied to head data sampled from a series of two-dimensional, simulated aquifers with isotropic, exponential covariance structures and varying degrees of heterogeneity (σY2 = 0.25, 0.5, 1.0, 2.0, and 5.0). The results show that this method provided reliable estimates for both λ and σY2 in aquifers with the value of σY2 up to 2.0, but the errors in those estimates were higher for σY2 equal to 5.0. It is also demonstrated through numerical experiments and theoretical arguments that the head residual increments will provide a sample semivariogram with a lower variance than will the use of the head residuals without calculation of increments.

We investigate from a theoretical basis the impacts of the number, location, and correlation among measurement points on the quality of an estimate of the semivariogram. The unbiased nature of the semivariogram estimator γˆ(r) is first established for a general random process Z(x). The variance of γˆZ(r) is then derived as a function of the sampling parameters (the number of measurements and their locations). In applying this function to the case of estimating the semivariograms of the transmissivity and the hydraulic head field, it is shown that the estimation error depends on the number of the data pairs, the correlation among the data pairs (which, in turn, are determined by the form of the underlying semivariogram γ(r)), the relative locations of the data pairs, and the separation distance at which the semivariogram is to be estimated. Thus design of an optimal sampling program for semivariogram estimation should include consideration of each of these factors. Further, the function derived for the variance of γˆZ(r)) is useful in determining the reliability of a semivariogram developed from a previously established sampling design.

Permeability upscaling is physically investigated by making over 31,000 permeability measurements on a meter-scale block of volcanic tuff. The experiments are made possible by a specially adapted minipermeameter test system. Here we present and analyze 5185 permeability values, corresponding to five different sample supports (i.e., sample volumes) collected from one of the six block faces. The results show that the measured spatial permeability patterns, bimodal permeability distribution, and semivariogram structure/length scales are closely related to the strong textural contrast characterizing the tuff sample (i.e., highly porous pumice fragments embedded in a tight rock matrix). Each of the summary statistics shows distinct and consistent trends with increasing sample support (i.e., upscaling). As the sample support increases, the mean and variance decrease according to a power law relation, and the semivariogram range increases linearly, while the general structure of the semivariogram (isotropic, spherical model) remains unchanged. Interpretation of these results is pursued from two very different points of view; one addresses upscaling of the ensemble statistics, while the second examines upscaling from a local or pointwise perspective. We find the general upscaling trends exhibited by the ensemble statistics (given above) to be consistent with the basic concepts of volume averaging, albeit nonlinear volume averaging. The bimodal characteristics of the tuff sample and the nonuniform flow conditions imparted by the minipermeameter contribute to the nonlinearity. The local analysis reveals strong variability in permeability upscaling from point to point throughout the sampling domain. Specifically, the permeability upscaling exhibited by zones rich in pumice is very different from zones dominated by matrix, unless the averaging volume is significantly larger than the spatial correlation scale.

It has been shown by Di Federico and Neuman [1997, 1998a, b] that observed multiscale behaviors of subsurface fluid flow and transport variables can be explained within the context of a unified stochastic framework, which views hydraulic conductivity as a random fractal characterized by a power variogram. Such a random field is statistically nonhomogeneous but possesses homogeneous spatial increments. Di Federico and Neuman [1997] have demonstrated that the power variogram and associated spectra of a statistically isotropic fractal field can be constructed as a weighted integral from zero to infinity (an infinite hierarchy) of exponential or Gaussian variograms and spectra of mutually uncorrelated fields (modes) that are homogeneous and isotropic. We show in this paper that the same holds true when the field and its constituent modes are statistically anisotropic, provided the ratios between principal integral (spatial correlation) scales are the same for all modes. We then analyze the effect of filtering out (truncating) modes of low, high, and intermediate spatial frequency from this infinite hierarchy in the real and spectral domains. A low-frequency cutoff renders the truncated hierarchy homogeneous. The integral scales of the lowest- and highest-frequency cutoff modes are related to length scales of the sampling window (domain) and data support (sample volume), respectively. Taking the former to be proportional to the latter renders expressions for the integral scale and variance of the truncated field dependent on window and support scale (in a manner previously shown to be consistent with observations in the isotropic case). It also allows (in principle) bridging across scales at a specific locale, as well as among locales, by adopting either site-specific or generalized variogram parameters. The introduction of intermediate cutoffs allows us to account, in a straightforward manner, for lacunarity due to gaps in the multiscale hierarchy created by the absence of modes associated with discrete ranges of scales (for example, where textural and structural features are associated with distinct ranges of scale, such as fractures having discrete ranges of trace length and density, which dissect the rock into matrix blocks having corresponding ranges of sizes). We explore mathematically and graphically the effects that anisotropy and lacunarity have on the integral scale, variance, covariance, and spectra of a truncated fractal field. We then develop an expression for the equivalent hydraulic conductivity of a box-shaped porous block, embedded within a multiscale log hydraulic conductivity field, under mean-uniform flow. The block is larger than the support scale of the field but is smaller than a surrounding sampling window. Consequently, its equivalent hydraulic conductivity is a random variable whose variance and spatial autocorrelation function, conditioned on a known mean value of support-scale conductivity across the window, are given explicitly by our multiscale theory.

The Washita '92 experiment, which took place in the Little Washita River watershed in Oklahoma, generated maps describing the volumetric surface soil moisture field during an 8-day drying period. The Little Washita is part of the Red River basin, a major tributary of the Lower Mississippi. A semivariogram analysis revealed that the surface soil moisture field of this 18 km × 25 km study region as a whole was spatially nonstationary and exhibited certain fractal characteristics. However, for smaller subareas within the study region, it was determined that the moisture field could be considered stationary for distances smaller than 5 km, with ranges of approximately 1 km over which surface soil moisture was spatially correlated. A semivariogram analysis for each of the individual soil types over the entire Washita '92 region revealed that surface soil moisture was stationary within each soil type, again with a range of correlation of approximately 1 km. This indicates that soil moisture variability is strongly affected by soil type.

A methodology for the choice of optimal spatial sampling designs is proposed. The problem is to find an optimal finite set of space locations where a random field has to be sampled, in order to minimize the variability of the parametric variogram estimator. Under the hypothesis that the random field is Gaussian second-order stationary and using a first-order approximation for the linearization of the parametric variogram with respect to its parameters, we are able to express the covariance matrix of the parameter estimators as a function of the sampling design. An optimization algorithm based on a generalized least squares approach is proposed in order to reduce the variability of these estimators through the minimization of the determinant of their covariance matrix. In order to validate this approach, a practical case study is conducted for two different variogram models. It shows that compared to random sampling designs, the benefit of the optimization procedure is somewhat limited for a variogram model without nugget effect. For a variogram model with nugget effect, random sampling designs are associated with much higher variability on the parameter estimators than optimized designs because of the typical lack of information offered by random sampling designs for small distances between locations. The performance of random sampling designs compared to optimized or alternative regular designs is shown to be poor with respect to parameter estimation, especially when a nugget effect is included in the variogram model. The parameters that benefit the most from the optimization procedure are the variance and the nugget effect, whereas the improvement for the range parameter estimation is limited. The optimization algorithm provides yardstick results, yielding reference values for the selection of an alternative regular design. The applicability of the algorithm is very wide, and it can greatly help the user to understand the way the parameters are influenced by the choice of a set of sampling locations instead of another.

Verification of a stochastic model which models the phenomena of flow and transport in a nonstationary conductivity field [Rubin and Seong, 1994] is attempted with field data from the Waste Isolation Pilot Plant (WIPP) site near Carlsbad, New Mexico. The nonstationarity is manifested as a spatial linear trend in the mean log conductivity field. Analysis of the log conductivity and the hydraulic head data shows that the site corresponds to Rubin and Seong's α2 model, where the mean flow direction and the gradient of mean log conductivity are normal to each other. The experimental semivariogram of the log conductivity obtained from field data strongly suggests that it is nonstationary and can be fit with a Gaussian covariance function. The α2 model is successful in predicting the experimental semivariogram of the head obtained from the field data, especially in capturing the finite sill characteristics. The stationary model predicts an unbounded semivariogram at large distances. Solute travel times are also investigated by comparing the travel time cumulative distribution function (CDF) resulting from the stationary and the nonstationary hydraulic conductivity models. The stationary model predicts a higher CDF at early times and hence overestimates the very early travel times to the regulatory plane of compliance than the α2 model.

This study presents a centroid-based method of regional analysis employing the Kriging technique. A case study in streamflow estimation using Kriging was carried out in the Choushui Chi Basin in central Taiwan. A variable (runoff depth) distributed in space was viewed as a regionalized variable, and the representative value of the runoff depth was allocated on the centroid of the area considered. The results show that the spatial structure of the measurements in the experimental variogram is not purely random and the stationarity hypothesis is appropriate in estimating streamflow.

This paper describes the first major attempt to compare seven different inverse approaches for identifying aquifer transmissivity. The ultimate objective was to determine which of several geostatistical inverse techniques is better suited for making probabilistic forecasts of the potential transport of solutes in an aquifer where spatial variability and uncertainty in hydrogeologic properties are significant. Seven geostatistical methods (fast Fourier transform (FF), fractal simulation (FS), linearized cokriging (LC), linearized semianalytical )LS), maximum likelihood (ML), pilot point (PP), and sequential self-calibration (SS)) were compared on four synthetic data sets. Each data set had specific features meeting (or not) classical assumptions about stationarity, amenability to a geostatistical description, etc. The comparison of the outcome of the methods is based on the prediction of travel times and travel paths taken by conservative solutes migrating in the aquifer for a distance of 5 km. Four of the methods, LS, ML, PP, and SS, were identified as being approximately equivalent for the specific problems considered. The magnitude of the variance of the transmissivity fields, which went as high as 10 times the generally accepted range for linearized approaches, was not a problem for the linearized methods when applied to stationary fields; that is, their inverse solutions and travel time predictions were as accurate as those of the nonlinear methods. Nonstationarity of the “true” transmissivity field, or the presence of “anomalies” such as high-permeability fracture zones was, however, more of a problem for the linearized methods. The importance of the proper selection of the semivariogram of the log10 (T) field (or the ability of the method to optimize this variogram iteratively) was found to have a significant impact on the accuracy and precision of the travel time predictions. Use of additional transient information from pumping tests did not result in major changes in the outcome. While the methods differ in their underlying theory, and the codes developed to implement the theories were limited to varying degrees, the most important factor for achieving a successful solution was the time and experience devoted by the user of the method.

In a previous paper we offered an interpretation for the observation that the log hydraulic conductivity of geologic media often appears to be statistically homogeneous but with variance and integral scale which grow with the size of the region (window) being sampled. We did so by demonstrating that such behavior is typical of any random field with a truncated power (semi)variogram and that this field can be viewed as a truncated hierarchy of mutually uncorrelated homogeneous fields with either exponential or Gaussian spatial autocovariance structures. The low- and high-frequency cutoff scales λl and λu are related to the length scales of the sampling window and data support, respectively. We showed how this allows the use of truncated power variograms to bridge information about a multiscale random field across windows of different sizes, either at a given locale or between different locales. In this paper we investigate mean uniform steady state groundwater flow in unbounded domains where the log hydraulic conductivity forms a truncated multiscale hierarchy of Gaussian fields, each associated with an exponential autocovariance. We start by deriving an expression for effective hydraulic conductivity, as a function of the Hurst coefficient H and the cutoff scales in one-, two-, and three-dimensional domains which is qualitatively consistent with experimental data. We then develop leading-order analytical expressions for two- and three-dimensional autocovariance and cross-covariance functions of hydraulic head, velocity, and log hydraulic conductivity versus H, λl and λu; examine their behavior; and compare them with those corresponding to an exponential log hydraulic conductivity autocovariance. Our results suggest that it should be possible to bridge information about hydraulic heads and groundwater velocities across windows of disparate scales. In particular, when λl ≫ λu, the variance of head is infinite in two dimensions and grows in proportion to λl2+2H in three dimensions, while the variance and longitudinal integral scale of velocity grow in proportion to λl2H and λl, respectively, in both cases.

The process of designing a piezometer network for estimating the sample variogram of the hydraulic gradient is used to illustrate one strategy for designing a sampling network which is both appropriate for the estimation of spatial statistics and efficient in terms of the measurement instrument. Constraints which impacted this design can be broken into three categories: (1) constraints discussed by previous authors and related to increasing confidence in the estimate of the sample variogram, (2) new constraints related to the spatial distribution of sample location for the purpose of estimating the sample variogram, and (3) constraints related to the measurement instrument. After assessing a number of sampling strategies, including uniform grids, randomized grids, various geometric patterns, and scaled sampling patterns, an acceptable design for the present application was developed from a scaled sampling pattern which was modified to utilize the instrument more efficiently. In this situation the final design could not have been established except through consideration of both the statistics to be derived from the sampling effort and the instrument used to obtain the measurements.

Results of three meter-scale outcrop studies suggest a quantifiable relationship between fluvial bounding surfaces and the permeability correlation structure. Mapping of bounding surfaces and in situ measurements of permeability were obtained in three types of deposits: a low-energy channel deposit, a paleosol, and a high-energy channel deposit. For the low-energy channel deposit only second-order surfaces separating lithofacies were mapped. The variogram exhibits anisotropic exponential behavior, and an analysis of variance (ANOVA) suggests that the within-facies variation dominates over between-facies variation. The paleosol contained no internal bounding surfaces and exhibited a “bellshaped” variogram consistent with a smoothly varying style of heterogeneity. Both firstand second-order bounding surfaces were mapped in the high-energy channel deposit. Results of the ANOVA suggest that the between-facies variation dominates over the within-facies variation. A model of spatial correlation using superposition of independent processes acting over different scales provides a good fit to the data. These three studies suggest that fluvial bounding surfaces provide a geological basis for conceptualizing and modeling heterogeneity in alluvial deposits.

The purpose of this work is to describe, evaluate, and demonstrate a laboratory-based method for physically investigating permeability upscaling. The method makes use of a simple instrument, the gas permeameter, to acquire rapid, precise, and nondestructive permeability measurements from heterogeneous blocks of dry rock. Critical to investigating permeability upscaling is the ability to acquire data at multiple sample supports subject to consistent boundary conditions and flow geometry. Such measurements, spanning almost 4 orders of magnitude on a per volume basis, are made with the gas permeameter by simply varying the size of the permeameter tip seal. The precision and consistency of measurements made in this way were evaluated using a suite of data collected from blocks of three relatively homogeneous materials: Berea Sandstone and two synthetic rocks. Results suggest that measurement error is small (approximately ±1% of the measured permeability) and consistent, and measurements made at different sample supports are free from systematic bias. To demonstrate the ability of this method to measure and quantify upscaling processes, limited data sets were collected with four different-sized tip seals from the Berea Sandstone block. Analysis reveals distinct and consistent trends diagnostic of permeability upscaling relating the sample mean (increased), variance (decreased), and semivariogram to increasing sample support.

An interpretation is offered for the observation that the log hydraulic conductivity of geologic media often appears to be statistically homogeneous but with variance and integral scale which grow with domain size. We first demonstrate that the power (semi)variogram and associated spectra of random fields, having homogeneous isotropic increments, can be constructed as weighted integrals from zero to infinity (an infinite hierarchy) of exponential or Gaussian variograms and spectra of mutually uncorrelated homogeneous isotropic fields (modes). We then analyze the effect of filtering out (truncating) high- and low-frequency modes from this infinite hierarchy in the real and spectral domains. A low-frequency cutoff renders the truncated hierarchy homogeneous with an autocovariance function that varies monotonically with separation distance in a manner not too dissimilar than that of its constituent modes. The integral scales of the lowest- and highest-frequency modes (cutoffs) are related, respectively, to the length scales of the sampling window (domain) and data support (sample volume). Taking each relationship to be one of proportionality renders our expressions for the integral scale and variance of a truncated field dependent on window and support scales in a manner consistent with observations. The traditional approach of truncating power spectral densities yields autocovariance functions that oscillate about zero with finite (in one and two dimensions) or vanishing (in one dimension) integral scales. Our hierarchical theory allows bridging across scales at a specific locale, by calibrating a truncated variogram model to data observed on a given support in one domain and predicting the autocovariance structure of the corresponding multiscale field in domains that are either smaller or larger. One may also venture (we suspect with less predictive power) to bridge across both domain scales and locales by adopting generalized variogram parameters derived on the basis of juxtaposed hydraulic and tracer data from many sites.

Recent studies have shown that a useful scheme to characterize subsurface heterogeneity is to use fractional Brownian motion (fBm) and fractional Gaussian noise (fGn). Commonly used techniques to detect these fractal distributions include rescaled range, spectral, and semivariogram analyses. Although rescaled range analysis is considered to be superior to spectral or semivariogram techniques, it cannot be used to clearly discriminate fGn and fBm from a given data set in some cases. On the basis of the fact that both adjusted range and rescaled range analyses, for a true fGn data set, should yield the same slope value in the relevant log-log plots, a simple scheme is developed to discriminate fGn and fBm using a combination of the two analyses. The usefulness of the proposed scheme is evaluated with one-dimensional computational experiments and by an application to a data set of observed air permeability logs. Lack of clear discrimination may have contributed to inconsistent past results, and it is recommended that old data sets be reanalyzed.

This study evaluates commonly used geostatistical methods to assess reproduction of hydraulic conductivity (K) structure and sensitivity under limiting amounts of data. Extensive conductivity measurements from the Cape Cod sand and gravel aquifer are used to evaluate two geostatistical estimation methods, conditional mean as an estimate and ordinary kriging, and two stochastic simulation methods, simulated annealing and sequential Gaussian simulation. Our results indicate that for relatively homogeneous sand and gravel aquifers such as the Cape Cod aquifer, neither estimation methods nor stochastic simulation methods give highly accurate point predictions of hydraulic conductivity despite the high density of collected data. Although the stochastic simulation methods yielded higher errors than the estimation methods, the stochastic simulation methods yielded better reproduction of the measured In (K) distribution and better reproduction of local contrasts in In (K). The inability of kriging to reproduce high In (K) values, as reaffirmed by this study, provides a strong instigation for choosing stochastic simulation methods to generate conductivity fields when performing fine-scale contaminant transport modeling. Results also indicate that estimation error is relatively insensitive to the number of hydraulic conductivity measurements so long as more than a threshold number of data are used to condition the realizations. This threshold occurs for the Cape Cod site when there are approximately three conductivity measurements per integral volume. The lack of improvement with additional data suggests that although fine-scale hydraulic conductivity structure is evident in the variogram, it is not accurately reproduced by geostatistical estimation methods. If the Cape Cod aquifer spatial conductivity characteristics are indicative of other sand and gravel deposits, then the results on predictive error versus data collection obtained here have significant practical consequences for site characterization. Heavily sampled sand and gravel aquifers, such as Cape Cod and Borden, may have large amounts of redundant data, while in more common real world settings, our results suggest that denser data collection will likely improve understanding of permeability structure.

The ability of passive microwave remote sensing (PMRS) to assess soil moisture has been amply demonstrated. Building on the expertise, this research is a first attempt at obtaining near-surface soil hydraulic characteristics from PMRS. We show that relative change in surface soil moisture can be determined accurately from PMRS measurements of brightness temperature, given the commonly observed linear relationship between the two variables. Relative, rather than absolute, change within a series of PMRS images, obtained during a period of drying, is also a very useful statistic for filtering out the drydown heterogeneity caused by spatial variability in initial wetness rather than soil heterogeneity. We develop a test of soil heterogeneity based on the detection of deviations of relative change from the average rate in a series of images by more than the instrument error. With the test, the existence of soil heterogeneity for a pair of images showing a drydown during the 1992 Hydrologic Atmospheric Pilot Experiment (HAPEX-Sahel) experiment can be proven. Areas of relatively fast or slow drying can be mapped. These areas coincide with geomorphological features and soil-type patterns that are expected to show the observed moisture dynamics. The variogram of relative change in the HAPEX images differs from those of brightness temperature, corroborating the rejection of the soil homogeneity hypothesis. Furthermore, the variogram indicates the existence of an effective correlation length much larger than that commonly observed in ground-based soil surveys. These findings encourage the use of PMRS for assessing soil hydraulic characteristics that are valid at a scale appropriate for hydrometeorological models.

Past rescaled range analyses of porosity and hydraulic conductivity (K) distributions have indicated the presence of long-range correlations in the data typical of the related stochastic functions known as fractional Gaussian noise and fractional Brownian motion [Hewett, 1986; Molz and Boman, 1993]. New K data analyzed herein lend further support to this notion. Horizontal processes that mimic fBm will display a power-law variogram. The Mandelbrot-Weierstrass random fractal function is introduced as an analytical model for fBm and used to illustrate several concepts. With the exception of the Hurst coefficient value (H), our analysis supports the existence of fractal-like K fields similar to those visualized by Neuman [1990, 1994]. Past studies which produced H values greater than or less than 0.5 appear to differ mainly because of the underlying model, fractional motion or fractional noise, that was assumed in the respective analyses.

A new methodology for solution of the inverse problem in groundwater hydrology is proposed and applied to a site in southeastern New Mexico with extensive hydrogeologic data. The methodology addresses the issue of nonuniqueness of the inverse solutions by generating an ensemble of transmissivity fields considered to be equally likely, each of which is in agreement with the measured transmissivity and pressure data. It consists of generating a selected number of conditionally simulated transmissivity fields and then calibrating each of the fields to match the measured steady state or transient pressures, in a least squares sense. The calibration phase involves an iterative implementation of an automated pilot point approach coupled with conditional simulations. Pilot points are the parameters of calibration. They are synthetic transmissivity data which are added to the transmissivity database to produce a revised conditional simulation during calibration. Coupled kriging and adjoint sensitivity analysis is employed for the optimal location of pilot points, and gradient search methods are used to derive their optimal transmissivities. The pilot point methodology is well suited for characterizing the spatial variability of the transmissivity field in contrast to methods using zonation. Pilot points are located where their potential for minimizing the objective function is the highest. This minimizes the perturbations in the transmissivities which are optimally assigned to the pilot point and results in minimal changes to the covariance structure of the transmissivity field. The calibrated fields honor the transmissivity measurements at their locations, preserve the variogram, and match the measured pressures in a least squares sense. Since there are numerous options in the execution of this methodology, computational experiments have been conducted to identify the most efficient among them. The method has been applied to the Waste Isolation Pilot Plant (WIPP) site, in southeastern New Mexico, where the U.S. Department of Energy is conducting probabilistic system assessment for the permanent disposal of transuranic nuclear waste. The resulting calibrated transmissivity fields are input to a Monte Carlo analysis of the total system performance. The present paper, paper 1 of a two-paper presentation, describes the methodology. Paper 2, a companion paper, presents the methodology's application to the WIPP site.

[1] This is an application of a spatial copula model that is fitted to a real world data set. The copula model allows modeling of pure spatial dependence independently of the marginal distribution. Using non-Gaussian copula models it is demonstrated that the spatial dependence structure of the Borden aquifer is significantly non-Gaussian—despite the fact that the Borden aquifer is commonly thought of as a relatively homogeneous porous medium with a small variance of hydraulic conductivity. In addition to evaluating the spatial dependence structure of the Borden hydraulic conductivity data set using copulas, goal of this study is to explore if the structure of the hydraulic conductivity field influences a physical property, such as plume evolution as evaluated by second spatial central moments of concentration fields. For this comparison, two types of hydraulic conductivity fields were fitted to the Borden hydraulic conductivity data set: one with a Gaussian and the other with a non-Gaussian type of dependence. These two types of hydraulic conductivity fields were constructed such that their second-order spatial moments are identical, and hence they cannot be distinguished by semivariogram-based geostatistics. This paper illustrates that the spatial dependence structure of the Borden hydraulic conductivity data set is significantly non-Gaussian. Despite the fact that Borden is a relatively homogeneous porous medium, and despite the fact that both types of spatial fields are not distinguishable by their variograms, the solute transport characteristics based on these two types of isotropic fields differ significantly in two-dimensional settings. The difference is less pronounced in three-dimensions with anisotropy. It is postulated that non-Gaussian spatial dependence of hydraulic conductivity and a more skewed marginal distribution of hydraulic conductivity will have significant implications in the other more heterogeneous aquifers.

[1] In developing a reliable approach for inferring hydrological properties through inverse modeling of tracer data, decisions made on how to parameterize heterogeneity (i.e., how to represent a heterogeneous distribution using a limited number of parameters that are amenable to estimation) are of paramount importance, as errors in the model structure are partly compensated for by estimating biased property values during the inversion. These biased estimates, while potentially providing an improved fit to the calibration data, may lead to wrong interpretations and conclusions and reduce the ability of the model to make reliable predictions. We consider the estimation of spatial variations in permeability and several other parameters through inverse modeling of tracer data, specifically synthetic and actual field data associated with the 2007 Winchester experiment from the Department of Energy Rifle site. Characterization is challenging due to the real-world complexities associated with field experiments in such a dynamic groundwater system. Our aim is to highlight and quantify the impact on inversion results of various decisions related to parameterization, such as the positioning of pilot points in a geostatistical parameterization; the handling of up-gradient regions; the inclusion of zonal information derived from geophysical data or core logs; extension from 2-D to 3-D; assumptions regarding the gradient direction, porosity, and the semivariogram function; and deteriorating experimental conditions. This work adds to the relatively limited number of studies that offer guidance on the use of pilot points in complex real-world experiments involving tracer data (as opposed to hydraulic head data).

[1] From the Permian through to the modern day, stalagmites are an important archive of environmental change. Annually laminated stalagmites provide both a precise chronology and a paleoclimate proxy. The rate of annual vertical growth of stalagmites is recorded in changes of calcite fabric, annual fluxes of fluorescent organic matter or annual variations in trace element composition. The processes governing stalagmite growth are the flux of water, the CO2 saturation of drip water relative to the cave atmosphere, and the temperature. Although these processes are well understood, they depend on the specific hydrogeological flow routing of individual stalagmites. Therefore, although past climates are recorded in the vertical growth lamina thickness, the climatic signal is perturbed by noise related to local hydrologic factors. To separate local from global factors, we used geostatistical tools to analyze annual growth rate data from eleven stalagmites located on four continents. Variogram analyses permit the quantification of the signal content contained within the growth rate records. The information content ranges from 23 to 87%. Analysis of the growth derivative shows a negative correlation at a 1 year lag, meaning that acceleration in growth rate tends to be systematically followed by deceleration in growth rate and vice versa. We call this behavior “flickering” growth, and argue that it is related to the size of the store feeding the stalagmite. Variogram analysis and flickering are used to screen which types of signals can potentially be recorded in a given speleothem.

Spatial variability of surface hydraulic properties and the extrinsic (e.g., traffic, cropping, etc.) and intrinsic (e.g., soil type, pore size distribution, etc.) factors associated with these properties are important for infiltration and runoff processes in agricultural fields. Disc infiltrometers measured infiltration at 296 sites arranged on two parallel transects. To examine and differentiate the factors contributing to spatial structure under different field conditions these measurements were made in the corn rows, no-track interrows, and wheel track interrows of the field using four different soil water tensions * (0, 30, 60, and 150 mm). Unsaturated hydraulic conductivity (K) and saturated hydraulic conductivity (K,) were maximum in the corn rows and minimum in wheel track inter-rows, with no-track interrows intermediate. Exponents (CZ parameters) of K, and K relationships (K = K, exp -aV) for corn rows and no-track interrows were not significantly different from each other but were significantly different from Q for the wheel track interrows a t P = 0.01 level. Spatial variability o f K and K, values showed some pseudoproportional effect in nugget variance for all three field conditions. No-track interrows clearly showed an inverse trend for semivariogram o f K with changing tension (q) values, whereas differences were found for corn rows and wheel traffic interrows. The spatial structure of (Y for all three field conditions were mostly white noise. Under corn rows, in addition to random variation, a small five-row periodic variation at th e P = 0.20 level, matching the five-row traffic configuration, was discovered. The spatial structure of a was influenced by soil type for the no-track interrows. Spatial structure was absent in wheel track interrows, indicating the destruction of pore structure due to compaction.

The permeabilities and dispersivities of geologic media are known to vary with the scale of observation. Particularly well documented is the consistent increase in apparent longitudinal dispersivity with the mean travel distance of a tracer. This has been previously interpreted by the author to imply that the permeabilities of many geologic media scale, on the average, according to the power-law semivariogram γ(s)=c ν s where c is a constant and s is distance. Tracer test data support this conclusion indirectly at least over scales from 10 cm to 3,500 m. The present paper cites evidence for such behavior over scales from 10 cm to 45 km based directly on permeability and transmissivity data

[1] We analyze theoretically the ability of model quality (sometimes termed information or discrimination) criteria such as the negative log likelihood NLL, Bayesian criteria BIC and KIC and information theoretic criteria AIC, AICc, and HIC to estimate (1) the parameter vector of the variogram of hydraulic log conductivity (Y = ln K), and (2) statistical parameters and proportional to head and log conductivity measurement error variances, respectively, in the context of geostatistical groundwater flow inversion. Our analysis extends the work of Hernandez et al. (2003, 2006) and Riva et al. (2009), who developed nonlinear stochastic inverse algorithms that allow conditioning estimates of steady state and transient hydraulic heads, fluxes and their associated uncertainty on information about conductivity and head data collected in a randomly heterogeneous confined aquifer. Their algorithms are based on recursive numerical approximations of exact nonlocal conditional equations describing the mean and (co)variance of groundwater flow. Log conductivity is parameterized geostatistically based on measured values at discrete locations and unknown values at discrete “pilot points.” Optionally, the maximum likelihood function on which the inverse estimation of Y at pilot points is based may include a regularization term reflecting prior information about Y. The relative weight assigned to this term and its components and , as well as are evaluated separately from other model parameters to avoid bias and instability. This evaluation is done on the basis of criteria such as NLL, KIC, BIC, HIC, AIC, and AICc. We demonstrate theoretically that, whereas all these six criteria make it possible to estimate , KIC alone allows one to estimate validly and (and thus ). We illustrate this discriminatory power of KIC numerically by using a differential evolution genetic search algorithm to minimize it in the context of a two-dimensional steady state groundwater flow problem. We find that whereas , and the integral scale of Y can be estimated on the basis of a zero-order mean flow equation, the sill of the Y-variogram is estimated more accurately by a second-order approximation of flow. This notwithstanding, KIC prefers the simpler zero-order moment over the more complex second-order version.

[1] This paper evaluates the use of field data on the spatial variability of snow water equivalent (SWE) to guide the design of distributed snow models. An extensive reanalysis of results from previous field studies in different snow environments around the world is presented, followed by an analysis of field data on spatial variability of snow collected in the headwaters of the Jollie River basin, a rugged mountain catchment in the Southern Alps of New Zealand. In addition, area-averaged simulations of SWE based on different types of spatial discretization are evaluated. Spatial variability of SWE is shaped by a range of different processes that occur across a hierarchy of spatial scales. Spatial variability at the watershed-scale is shaped by variability in near-surface meteorological fields (e.g., elevation gradients in temperature) and, provided suitable meteorological data is available, can be explicitly resolved by spatial interpolation/extrapolation. On the other hand, spatial variability of SWE at the hillslope-scale is governed by processes such as drifting, sloughing of snow off steep slopes, trapping of snow by shrubs, and the nonuniform unloading of snow by the forest canopy, which are more difficult to resolve explicitly. Subgrid probability distributions are often capable of representing the aggregate-impact of unresolved processes at the hillslope-scale, though they may not adequately capture the effects of elevation gradients. While the best modeling strategy is case-specific, the analysis in this paper provides guidance on both the suitability of several common snow modeling approaches and on the choice of parameter values in subgrid probability distributions.

[1] This study investigates connectivity in a small portion of the extremely heterogeneous aquifer at the Macrodispersion Experiment (MADE) site in Columbus, Mississippi. A total of 19 fully penetrating soil cores were collected from a rectangular grid of 4 m by 4 m. Detailed grain size analysis was performed on 5 cm segments of each core, yielding 1740 hydraulic conductivity (K) estimates. Three different geostatistical simulation methods were used to generate 3-D conditional realizations of the K field for the sampled block. Particle tracking calculations showed that the fastest particles, as represented by the first 5% to arrive, converge along preferential flow paths and exit the model domain within preferred areas. These 5% fastest flow paths accounted for about 40% of the flow. The distribution of preferential flow paths and particle exit locations is clearly influenced by the occurrence of clusters formed by interconnected cells with K equal to or greater than the 0.9 decile of the data distribution (10% of the volume). The fraction of particle paths within the high-K clusters ranges from 43% to 69%. In variogram-based K fields, some of the fastest paths are through media with lower K values, suggesting that transport connectivity may not require fully connected zones of relatively homogenous K. The high degree of flow and transport connectivity was confirmed by the values of two groups of connectivity indicators. In particular, the ratio between effective and geometric mean K (on average, about 2) and the ratio between the average arrival time and the arrival time of the fastest particles (on average, about 9) are consistent with flow and advective transport behavior characterized by channeling along preferential flow paths.

[1] The ensemble Kalman filter (EnKF) has recently been proposed as a promising parameter estimation approach for constraining the description of rock flow properties, such as permeability and porosity, to reproduce flow measurements that are modeled as nonlinear functions of these properties. One of the key factors that strongly affect the performance of the EnKF is the quality or representativeness of the prior ensemble of property fields used to initialize the EnKF assimilation procedure. The initial ensemble is commonly constructed by assuming a known geological continuity model such as a variogram. However, geologic continuity models are derived from incomplete information and imperfect modeling assumptions, which can introduce a significant level of uncertainty into the produced models. Neglecting this important source of uncertainty can lead to systematic errors and questionable estimation results. In this paper, we investigate the performance of the EnKF under varying levels of uncertainty in the variogram model parameters. We first attempt to directly estimate variogram model parameters from flow data and show that the complex and nonunique relation they have with the flow data provides little sensitivity for an effective inversion with the EnKF. We then assess the performance of the EnKF for estimation of permeability values under uncertain and incorrect initial variogram parameters and show that any bias in specifying variogram parameters tends to persist throughout the EnKF analysis even though locally reasonable permeability updates may be obtained near observation points. More importantly, we show that when variogram parameters are specified probabilistically to account for the full range of structural variability in the initial permeability ensemble, the EnKF update results are quite promising. The results suggest that under uncertain geologic continuity, the EnKF tends to perform better if a very diverse set of property fields is used to form the initial ensemble than when a deterministic and potentially erroneous variogram model is used. Therefore, in applying the EnKF to model calibration problems, it is preferable to overestimate the uncertainty in geologic continuity and to initialize the EnKF procedure with a wide range of variability in property description than to overlook the variogram uncertainty at the risk of introducing systematic bias that cannot be corrected by the EnKF updates. The practical implications of the results in this paper are significant for designing the EnKF for realistic ensemble model calibration problems where the level of uncertainty in the initial ensemble is usually not known a priori.

[1] Because of the extensive computational burden and perhaps a lack of awareness of existing methods, rigorous uncertainty analyses are rarely conducted for variable-density flow and transport models. For this reason, a recently developed null-space Monte Carlo (NSMC) method for quantifying prediction uncertainty was tested for a synthetic saltwater intrusion model patterned after the Henry problem. Saltwater intrusion caused by a reduction in fresh groundwater discharge was simulated for 1000 randomly generated hydraulic conductivity distributions, representing a mildly heterogeneous aquifer. From these 1000 simulations, the hydraulic conductivity distribution giving rise to the most extreme case of saltwater intrusion was selected and was assumed to represent the “true” system. Head and salinity values from this true model were then extracted and used as observations for subsequent model calibration. Random noise was added to the observations to approximate realistic field conditions. The NSMC method was used to calculate 1000 calibration-constrained parameter fields. If the dimensionality of the solution space was set appropriately, the estimated uncertainty range from the NSMC analysis encompassed the truth. Several variants of the method were implemented to investigate their effect on the efficiency of the NSMC method. Reducing the dimensionality of the null-space for the processing of the random parameter sets did not result in any significant gains in efficiency and compromised the ability of the NSMC method to encompass the true prediction value. The addition of intrapilot point heterogeneity to the NSMC process was also tested. According to a variogram comparison, this provided the same scale of heterogeneity that was used to generate the truth. However, incorporation of intrapilot point variability did not make a noticeable difference to the uncertainty of the prediction. With this higher level of heterogeneity, however, the computational burden of generating calibration-constrained parameter fields approximately doubled. Predictive uncertainty variance computed through the NSMC method was compared with that computed through linear analysis. The results were in good agreement, with the NSMC method estimate showing a slightly smaller range of prediction uncertainty than was calculated by the linear method.

This paper describes a geostatistical method, known as factorial kriging analysis, which is well suited for analyzing multivariate spatial information. The method involves multivariate variogram modeling, principal component analysis, and cokriging. It uses several separate correlation structures, each corresponding to a specific spatial scale, and yields a set of regionalized factors summarizing the main features of the data for each spatial scale. The method is applied to springwater solute contents (Ca, Sr, Mg, K, Na, SO4, NO3 Cl, pH, electrical conductivity, alkalinity) measured in 86 springs situated in the Dyle River catchment area, Belgium. Two scales of spatial variation (1 and 9 km) are identified and interpreted. These correspond to local sources of contaminants due to human activities and regional changes in the geologic characteristics of the aquifer. At each scale, the correlation structure is analyzed and the regionalized factors are estimated by cokriging and then mapped. The results are compared with those of a classical multivariate approach, i.e., principal component analysis of the correlation matrix, which takes no account of the spatial location of the observations.

Received 3 August 2009; revised 12 July 2010; accepted 21 July 2010; published 4 November 2010. [1] We present novel equations for the exact sensitivity matrix of the (ensemble) mean hydraulic head under steady state groundwater flow conditions. These equations are embedded in a geostatistical inverse procedure to condition approximations of stochastic moment equations of flow on measured hydraulic conductivities and heads. Our formulation allows considerable improvement of the methodology proposed by Hernandez et al. (2003, 2006) and renders the inversion of moment equations feasible for a large number of unknown hydraulic parameters. The spatial distribution of the natural logarithm, Y, of conductivity is parameterized within the pilot points framework. Whereas prior values of Y at pilot points are obtained by a variant of kriging, posterior estimates at pilot points are obtained through a maximum likelihood fit of computed to measured heads. The maximum likelihood function also includes a regularization term. By means of a synthetic example and upon adopting formal model information criteria we explore the influence of (1) the number of pilot points and (2) the order of approximation of the governing mean flow equation on our ability to properly estimate the log conductivity and head fields and identify the relative weight of the regularization term and the parameters of the underlying Y variogram. We find that none of the adopted information criteria can identify the optimum number of pilot points and the plausibility weight and variogram parameters values can be determined by the Kashyap’s Bayesian measure.

[1] The goal of simulation of aquifer heterogeneity is to produce a spatial model of the subsurface that represents a system such that it can be used to understand or predict flow and transport processes. Spatial simulation requires incorporation of data and geologic knowledge, as well as representation of uncertainty. Classical geostatistical techniques allow for the conditioning of data and uncertainty assessment, but models often lack geologic realism. Simulation of physical geologic processes of sedimentary deposition and erosion (process-based modeling) produces detailed, geologically realistic models, but conditioning to local data is limited at best. We present an aquifer modeling methodology that combines geologic-process models with object-based, multiple-point, and variogram-based geostatistics to produce geologically realistic realizations that incorporate geostatistical uncertainty and can be conditioned to data. First, the geologic features of grain size, or facies, distributions simulated by a process-based model are analyzed, and the statistics of feature geometry are extracted. Second, the statistics are used to generate multiple realizations of reduced-dimensional features using an object-based technique. Third, these realizations are used as multiple alternative training images in multiple-point geostatistical simulation, a step that can incorporate local data. Last, a variogram-based geostatistical technique is used to produce conditioned maps of depositional thickness and erosion. Successive realizations of individual strata are generated in depositional order, each dependent on previously simulated geometry, and stacked to produce a fully conditioned three-dimensional facies model that mimics the architecture of the process-based model. We demonstrate the approach for a typical subsea depositional complex.

In Part 1 of this series we presented a methodology and a code for modeling the hierarchical sedimentary architecture in braided channel belt deposits. Here, in Part 2, the code was used to create a digital model of this architecture, and the corresponding spatial distribution of permeability. The simulated architecture was compared to the real stratal architecture observed in an abandoned channel belt of the Sagavanirktok River, Alaska by Lunt et al. (2004). The comparisons included assessments of similarity which were both qualitative and quantitative. From the qualitative comparisons we conclude that a synthetic deposit created by the code has unit types, at each level, with a geometry which is generally consistent with the geometry of unit types observed in the field. The digital unit types would generally be recognized as representing their counterparts in nature, including cross stratasets, lobate and scroll bar deposits, channel fills, etc. Furthermore, the synthetic deposit has a hierarchical spatial relationship among these units which represents how the unit types are observed in field exposures and in geophysical images. In quantitative comparisons the proportions and the length, width, and height of unit types at different scales, across all levels of the stratal hierarchy compare wellmore » between the digital and the natural deposits. A number of important attributes of the channel belt model were shown to be influenced by more than one level within the hierarchy of stratal architecture. First, the high-permeability open-framework gravels percolated at all levels and thus formed preferential flow pathways. Open framework gravels are indeed known to form preferential flow pathways in natural channel belt deposits. The nature of a percolating cluster changed across different levels of the hierarchy of stratal architecture. As a result of this geologic structure, the percolation occurs at proportions of open-framework gravels below the theoretical percolation threshold for random infinite media. Second, when the channel belt model was populated with permeability distributions by lowest-level unit type, the composite permeability semivariogram contained structures that were identifiable at more than one scale, and each of these structures could be directly linked to unit types of different scales existing at different levels within the hierarchy of strata. These collective results are encouraging with respect to our goal that this model be relevant as a base case in future studies for testing ideas in research addressing the upscaling problem in aquifers and reservoirs with multi-scale heterogeneity.« less

[1] The spatial variability of hydraulic conductivity in a shallow unconfined aquifer located at North Bay, Ontario, composed of glacial-lacustrine and glacial-fluvial sands, is examined in exceptional detail and characterized geostatistically. A total of 1878 permeameter measurements were performed at 0.05 m vertical intervals along cores taken from 20 boreholes along two intersecting transect lines. Simultaneous three-dimensional (3-D) fitting of Ln(K) variogram data to an exponential model yielded geostatistical parameters for the estimation of bulk hydraulic conductivity and solute dispersion parameters. The analysis revealed a Ln(K) variance equal to about 2.0 and 3-D anisotropy of the correlation structure of the heterogeneity (λ1, λ2, and λ3 equal to 17.19, 7.39, and 1.0 m, respectively). Effective values of the hydraulic conductivity tensor and the value of the longitudinal macrodispersivity were calculated using the theoretical expressions of Gelhar and Axness (1983). The magnitude of the longitudinal macrodispersivity is reasonably consistent with the observed degree of longitudinal dispersion of the landfill plume along the principal path of migration. Variably saturated 3-D flow modeling using the statistically derived effective hydraulic conductivity tensor allowed a reasonably close prediction of the measured water table and the observed heads at various depths in an array of piezometers. Concomitant transport modeling using the calculated longitudinal macrodispersivity reasonably predicted the extent and migration rates of the observed contaminant plume that was monitored using a network of multilevel samplers over a period of about 5 years. It was further demonstrated that the length of the plume is relatively insensitive to the value of the longitudinal macrodispersivity under the conditions of a steady flow in 3-D and constant source strength. This study demonstrates that the use of statistically derived parameters based on stochastic theories results in reliable large-scale 3-D flow and transport models for complex hydrogeological systems. This is in agreement with the conclusions reached by Sudicky (1986) at the site of an elaborate tracer test conducted in the aquifer at the Canadian Forces Base Borden.