Adaptive-grid simulation of groundwater flow in heterogeneous aquifers

The prediction of contaminant transport in porous media requires the computation of the flow velocity. This work presents a methodology for high-accuracy computation of flow in a heterogeneous isotropic formation, employing a dual-flow formulation and adaptive gridding. The dual equations, describing the hydraulic head and the streamfunction, are numerically solved through finite element approximations. The application of classic finite-element methods requires a rather large number of nodes to represent suitably the flow in high-contrast formations. We present a mesh-adaptive approach that enhances the accuracy of the numerical flow solution for a given computational effort. We rely on an a posteriori error estimator to identify areas where refinements of the finite element mesh are needed or unrefinements are acceptable. We also demonstrate through numerical experiments that the developed methodology efficiently enhances accuracy through successive mesh adaptation.     

Random domain decomposition for flow in heterogeneous stratified aquifers

We study two-dimensional flow in a layered heterogeneous medium composed of two materials whose hydraulic properties and spatial distribution are known statistically but are otherwise uncertain. Our analysis relies on the composite media theory, which employs random domain decomposition in the context of groundwater flow moment equations to explicitly account for the separate effects of material and geometric uncertainty on ensemble moments of head and flux. Flow parallel and perpendicular to the layering in a two-material composite layered medium is considered. The hydraulic conductivity of each material is log-normally distributed with a much higher mean in one material than in the other. The hydraulic conductivities of points within different materials are uncorrelated. The location of the internal boundary between the two contrasting materials is random and normally distributed with given mean and variance. We solve the equations for (ensemble) moments of hydraulic head and flux and analyze the impact of unknown geometry of materials on statistical moments of head and flux. We compare the composite media approach to approximations that replace statistically inhomogeneous conductivity fields with pseudo-homogeneous random fields. This work was performed under the auspices of the US Department of Energy (DOE): DOE/BES (Bureau of Energy Sciences) Program in the Applied Mathematical Sciences contract KC-07–01–01 and Los Alamos National Laboratory under LDRD 98604. This work made use of STC shared experimental facilities supported by the National Science Foundation under Agreement No. EAR-9876800. This work was supported in part by the European Commission under Contract No. EVK1-CT-1999–00041 (W-SAHaRA).     

Dimensionality reduction for efficient Bayesian estimation of groundwater flow in strongly heterogeneous aquifers

We focus on the Bayesian estimation of strongly heterogeneous transmissivity fields conditional on data sampled at a set of locations in an aquifer. Log-transmissivity, Y, is modeled as a stochastic Gaussian process, parameterized through a truncated Karhunen–Loève (KL) expansion. We consider Y fields characterized by a short correlation scale as compared to the size of the observed domain. These systems are associated with a KL decomposition which still requires a high number of parameters, thus hampering the efficiency of the Bayesian estimation of the underlying stochastic field. The distinctive aim of this work is to present an efficient approach for the stochastic inverse modeling of fully saturated groundwater flow in these types of strongly heterogeneous domains. The methodology is grounded on the construction of an optimal sparse KL decomposition which is achieved by retaining only a limited set of modes in the expansion. Mode selection is driven by model selection criteria and is conditional on available data of hydraulic heads and (optionally) Y. Bayesian inversion of the optimal sparse KLE is then inferred using Markov Chain Monte Carlo (MCMC) samplers. As a test bed, we illustrate our approach by way of a suite of computational examples where noisy head and Y values are sampled from a given randomly generated system. Our findings suggest that the proposed methodology yields a globally satisfactory inversion of the stochastic head and Y fields. Comparison of reference values against the corresponding MCMC predictive distributions suggests that observed values are well reproduced in a probabilistic sense. In a few cases, reference values at some unsampled locations (typically far from measurements) are not captured by the posterior probability distributions. In these cases, the quality of the estimation could be improved, e.g., by increasing the number of measurements and/or the threshold for the selection of KL modes.     

Stochastic analysis of bounded unsaturated flow in heterogeneous aquifers: Spectral/perturbation approach

This paper describes a stochastic analysis of steady state flow in a bounded, partially saturated heterogeneous porous medium subject to distributed infiltration. The presence of boundary conditions leads to non-uniformity in the mean unsaturated flow, which in turn causes non-stationarity in the statistics of velocity fields. Motivated by this, our aim is to investigate the impact of boundary conditions on the behavior of field-scale unsaturated flow. Within the framework of spectral theory based on Fourier–Stieltjes representations for the perturbed quantities, the general expressions for the pressure head variance, variance of log unsaturated hydraulic conductivity and variance of the specific discharge are presented in the wave number domain. Closed-form expressions are developed for the simplified case of statistical isotropy of the log hydraulic conductivity field with a constant soil pore-size distribution parameter. These expressions allow us to investigate the impact of the boundary conditions, namely the vertical infiltration from the soil surface and a prescribed pressure head at a certain depth below the soil surface. It is found that the boundary conditions are critical in predicting uncertainty in bounded unsaturated flow. Our analytical expression for the pressure head variance in a one-dimensional, heterogeneous flow domain, developed using a nonstationary spectral representation approach [Li S-G, McLaughlin D. A nonstationary spectral method for solving stochastic groundwater problems: unconditional analysis. Water Resour Res 1991;27(7):1589–605; Li S-G, McLaughlin D. Using the nonstationary spectral method to analyze flow through heterogeneous trending media. Water Resour Res 1995; 31(3):541–51], is precisely equivalent to the published result of Lu et al. [Lu Z, Zhang D. Analytical solutions to steady state unsaturated flow in layered, randomly heterogeneous soils via Kirchhoff transformation. Adv Water Resour 2004;27:775–84].     

Describing near surface, transient flow processes in unconfined aquifers below irrigated lands: A deforming finite element model for heterogeneous aquifers

Hydrologic models of irrigated lands generally adopt either a basin-scale or a root-zone perspective. While basin-wide macro-scale models rely on the aggregation of important spatial and temporal data across large areas, micro-scale root-zone models depend on the definition of rigid boundaries around the zone of plant–soil–water interaction. In reality, irrigation management decisions are made on a field by field basis and can interact across field boundaries. This paper first describes a shallow water table model, based on deforming finite element (DFE) framework, to characterize the near-surface field-to-field hydrologic response to various irrigation and drainage management regimes along a gently sloping alluvial fan. The model is then enhanced through changing geometry of a fluctuating water table below a series of irrigated fields. Such an enhancement also offers computational flexibility relative to the saturated–unsaturated models commonly used in micro-scale studies. The model is designed with the alluvial fan aquifers of California’s western San Joaquin Valley as reference systems.     

Numerical study on tide-driven submarine groundwater discharge and seawater recirculation in heterogeneous aquifers

There are many factors affecting submarine groundwater discharge (SGD). However, systematic study of the influences of these factors is still limited. In this study, numerical modeling is performed to quantitatively explore the influences of various factors on SGD in a coastal aquifer. In such locations, tidal and terrestrial hydraulic gradients are the primary forces driving fresh and salt water movement. Unlike steady-state flow, dynamic fresh and salt water mixing at the near-shore seafloor may form an intertidal mixing zone (IMZ) near the surface. By constructing a general SGD model, the effects of various model components such as boundary conditions, model geometry and hydraulic parameters are systematically studied. Several important findings are obtained from the study results: (1) Previous studies have indicated there will be a freshwater discharge tube between the classic transition zone and the IMZ. However, this phenomenon may become unclear with the increase of heterogeneity and anisotropy of the medium’s conductivity field. (2) SGD and IMZ are both more sensitive to the vertical anisotropy ratio of hydraulic conductivity (K_x/K_z) than to the horizontal ratio (K_x/K_y). (3) Heterogeneity of effective porosity significantly affects SGD and IMZ. (4) Increase of the storage coefficient decreases fresh water discharge but increases mixing salt water discharge and total SGD. The increase will also change the shape of the IMZ. (5) Variation of dispersivities does not affect SGD, but significantly changes the distributions of the IMZ and the whole mixing zone. These findings will be helpful to the sampling design of field studies of SGD and to the application of dynamic SGD models to field sites for model development and calibration.     

Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: A comparative study

This work deals with a comparison of different numerical schemes for the simulation of contaminant transport in heterogeneous porous media. The numerical methods under consideration are Galerkin finite element (GFE), finite volume (FV), and mixed hybrid finite element (MHFE). Concerning the GFE we use linear and quadratic finite elements with and without upwind stabilization. Besides the classical MHFE a new and an upwind scheme are tested. We consider higher order finite volume schemes as well as two time discretization methods: backward Euler (BE) and the second order backward differentiation formula BDF (2). It is well known that numerical (or artificial) diffusion may cause large errors. Moreover, when the Péclet number is large, a numerical code without some stabilising techniques produces oscillating solutions. Upwind schemes increase the stability but show more numerical diffusion. In this paper we quantify the numerical diffusion for the different discretization schemes and its dependency on the Péclet number. We consider an academic example and a realistic simulation of solute transport in heterogeneous aquifer. In the latter case, the stochastic estimates used as reference were obtained with global random walk (GRW) simulations, free of numerical diffusion. The results presented can be used by researchers to test their numerical schemes and stabilization techniques for simulation of contaminant transport in groundwater.     

Positive solution of two-dimensional solute transport in heterogeneous aquifers

The transport of contaminants in aquifers is usually represented by a convection-dispersion equation. There are several well-known problems of oscillation and artificial dispersion that affect the numerical solution of this equation. For example, several studies have shown that standard treatment of the cross-dispersion terms always leads to a negative concentration. It is also well known that the numerical solution of the convective term is affected by spurious oscillations or substantial numerical dispersion. These difficulties are especially significant for solute transport in nonuniform flow in heterogeneous aquifers. For the case of coupled reactive-transport models, even small negative concentration values can become amplified through nonlinear reaction source/sink terms and thus result in physically erroneous and unstable results. This paper includes a brief discussion about how nonpositive concentrations arise from numerical solution of the convection and cross-dispersion terms. We demonstrate the effectiveness of directional splitting with one-dimensional flux limiters for the convection term. Also, a new numerical scheme for the dispersion term that preserves positivity is presented. The results of the proposed convection scheme and the solution given by the new method to compute dispersion are compared with standard numerical methods as used in MT3DMS.     

Spatial moments for colloid-enhanced radionuclide transport in heterogeneous aquifers

We consider colloid facilitated radionuclide transport by steady groundwater flow in a heterogeneous porous formation. Radionuclide binding on colloids and soil-matrix is assumed to be kinetically/equilibrium controlled. All reactive parameters are regarded as uniform, whereas the hydraulic log-conductivity is modelled as a stationary random space function (RSF). Colloid-enhanced radionuclide transport is studied by means of spatial moments pertaining to both the dissolved and colloid-bounded concentration. The general expressions of spatial moments for a colloid-bounded plume are presented for the first time, and are discussed in order to show the combined impact of sorption processes as well as aquifer heterogeneity upon the plume migration. For the general case, spatial moments are defined by the aid of two characteristic reaction functions which cannot be expressed analytically. By adopting the approximation for the longitudinal fluid trajectory covariance valid for a flow parallel to the formation bedding suggested by Dagan and Cvetkovic [Dagan G, Cvetkovic V. Spatial Moments of Kinetically Sorbing Plume in a Heterogeneous Aquifers. Water Resour Res 1993;29:4053], we obtain closed form solutions.     

Cross-borehole flowmeter tests for transient heads in heterogeneous aquifers

Cross-borehole flowmeter tests have been proposed as an efficient method to investigate preferential flowpaths in heterogeneous aquifers, which is a major task in the characterization of fractured aquifers. Cross-borehole flowmeter tests are based on the idea that changing the pumping conditions in a given aquifer will modify the hydraulic head distribution in large-scale flowpaths, producing measurable changes in the vertical flow profiles in observation boreholes. However, inversion of flow measurements to derive flowpath geometry and connectivity and to characterize their hydraulic properties is still a subject of research. In this study, we propose a framework for cross-borehole flowmeter test interpretation that is based on a two-scale conceptual model: discrete fractures at the borehole scale and zones of interconnected fractures at the aquifer scale. We propose that the two problems may be solved independently. The first inverse problem consists of estimating the hydraulic head variations that drive the transient borehole flow observed in the cross-borehole flowmeter experiments. The second inverse problem is related to estimating the geometry and hydraulic properties of large-scale flowpaths in the region between pumping and observation wells that are compatible with the head variations deduced from the first problem. To solve the borehole-scale problem, we treat the transient flow data as a series of quasi-steady flow conditions and solve for the hydraulic head changes in individual fractures required to produce these data. The consistency of the method is verified using field experiments performed in a fractured-rock aquifer.     

Non-Darcian groundwater flow in leaky aquifers

Abstract

A simplified method has been developed for solving leaky aquifer non-Darcian flow hydraulics. The principle of volumetric approach is combined with the confined-aquifer, time-dependent drawdown equation in an observation well. The groundwater flow in the leaky aquifer is assumed to obey a non-Darcian flow law of exponential type. The results are obtained in the form of type-curve expressions from which the necessary bundles of curves are drawn for a set of selective non-Darcian flow aquifer parameters. Although application of the methodology appears as rather limited but it provides a scientific contribution and extension of leaky aquifer theory towards nonlinear flow conditions. The methodology developed herein is applied to some actual field data from the eastern sedimentary basin in the Kingdom of Saudi Arabia.     

Monte Carlo evaluation of microbial-mediated contaminant reactions in heterogeneous aquifers

Monte Carlo simulations are conducted to evaluate microbial-mediated contaminant reactions in an aquifer comprised of spatially variable microbial biomass concentrations, aquifer hydraulic conductivities, and initial electron donor/acceptor concentrations. A finite element simulation model is used that incorporates advection, dispersion, and Monod kinetic expressions to describe biological processes. Comparisons between Monte Carlo simulations of heterogeneous systems and simulations using homogeneous formulation of the same two-dimensional transport problem are presented. For the assumed set of parameters, physical aquifer heterogeneity is found to have a minor effect on the mass of contaminant biodegraded/transformed when compared to a homogeneous system; however, it noticeably changes the dispersion, skewness, and peakness of contaminant concentration distributions. Similarly, for low microbial growth rate, given favorable microbial growth characteristics, biological heterogeneity has minor effect on the mass of contaminant biodegraded/transformed when compared to a homogeneous system. On the other hand, when higher effective growth rates are assumed, biological heterogeneity and spatial heterogeneities in essential electron donor/acceptors reduce the efficiency of biotic contaminant reactions; consequently, model simulations derived from heterogeneous biomass distributions predict remediation time scales that are longer than those simulated for homogeneous systems. When correlations between physical aquifer and biological heterogeneities are considered, the assumed correlation affects predicted mean and variance of contaminant concentration and biomass distributions. For example, an assumed negative correlation between hydraulic conductivity and the initial biomass distribution produces a plume where less efficient biotic contaminant reactions occur at the leading edge of the plume; this is consistent with less degradation/transformation occurring over regions of higher groundwater velocities. However, the presence and absence of these correlations do not appear to affect the efficiency of microbial-mediated contaminant attenuation.     

Comparison of formation and fluid-column logs in a heterogeneous basalt aquifer

Deep observation boreholes in the vicinity of active production wells in Honolulu, Hawaii, exhibit the anomalous condition that fluid-column electrical conductivity logs and apparent profiles of pore-water electrical conductivity derived from induction conductivity logs are nearly identical if a formation factor of 12.5 is assumed. This condition is documented in three boreholes where fluid-column logs clearly indicate the presence of strong borehole flow induced by withdrawal from partially penetrating water-supply wells. This result appears to contradict the basic principles of conductivity-log interpretation. Flow conditions in one of these boreholes was investigated in detail by obtaining flow profiles under two water production conditions using the electromagnetic flowmeter. The flow-log interpretation demonstrates that the fluid-column log resembles the induction log because the amount of inflow to the borehole increases systematically upward through the transition zone between deeper salt water and shallower fresh water. This condition allows the properties of the fluid column to approximate the properties of water entering the borehole as soon as the upflow stream encounters that producing zone. Because this condition occurs in all three boreholes investigated, the similarity of induction and fluid-column logs is probably not a coincidence, and may relate to aquifer response under the influence of pumping from production wells.     

Tide-induced head fluctuations in coastal aquifers of variable thickness

In this work, a new analytical solution to describe tide-induced head fluctuations in aquifers of variable thickness is presented. The proposed model assumes a finite and confined aquifer with a thickness that increases or decreases quadratically with the distance to the coast. A closed-form analytical solution is obtained by solving a boundary-value problem with both a separation of variables method and a change of variables method. This solution is a generalization of the solution obtained by Cuello et al., Hydrogeological Journal, 2017, 25, 1509–1515. The analytical solution is expressed in terms of the wedging parameter, a parameter that depends on the length and thicknesses at the coast and at the inland edge of the aquifer. Positive values of the wedging parameter describe aquifers with increasing thickness towards land and negative values describe aquifers with a decreasing thickness in the inland direction. The comparison of the new solution and the solution for a finite aquifer with constant thickness indicates that the sign of the wedging parameter enhances or decreases the amplitude of the tide-induced signal. However, the differences in time-lag between both solutions are negligible near the coast. The slope factor, which quantifies the inconsistencies between aquifer diffusivities estimated from attenuation and time-lag data, is computed and analysed. Near the coast, slope factor values greater than one are obtained for negative wedging parameters while slope factor values less than one are obtained for positive wedging parameters. The analysis of the new solution also indicates that more reliable estimates of the hydraulic diffusivity can be obtained from time-lag data.     

Assessing radial transmissivity variation in heterogeneous aquifers using analytical techniques

Pumping tests are one of the most commonly used in situ testing techniques for assessing aquifer hydraulic properties. Numerous researches have been conducted to predict the effects of aquifer heterogeneity on the groundwater levels during pumping tests. The objectives of the present work were as follows: (1) to predict drawdown conditions and to estimate aquifer properties during pumping tests undertaken in radially symmetric heterogeneous aquifers, and (2) to identify a method for assessing the transmissivity field along the radial coordinate in radially symmetric and fully heterogeneous transmissivity fields. The first objective was achieved by expanding an existing analytical drawdown formulation that was valid for a radially symmetric confined aquifer with two concentric zones around the pumping well to an N concentric zone confined aquifer having a constant transmissivity value within each zone. The formulation was evaluated for aquifers with three and four concentric zones to assess the effects of the transmissivity field on the drawdown conditions. The specific conditions under which aquifer properties could be identified using traditional methods of analysis were also evaluated. The second objective was achieved by implementing the inverse solution algorithm (ISA), which was developed for petroleum reservoirs to groundwater aquifer settings. The results showed that the drawdown values are influenced by a volumetric integral of a weighting function and the transmissivity field within the cone of depression. The weighting function migrates in tandem with the expanding cone of depression. The ability of the ISA to predict radially symmetric and log‐normally distributed transmissivity fields was assessed against analytical and numerical benchmarks. The results of this investigation indicated that the ISA method is a viable technique for evaluating the radial transmissivity variations of heterogeneous aquifer settings. Copyright © 2013 John Wiley & Sons, Ltd.     

Numerical simulations of transport of non-ergodic solute plumes in heterogeneous aquifers

Transport of non-ergodic solute plumes by steady-state groundwater flow with a uniform mean velocity, μ, were simulated with Monte Carlo approach in a two-dimensional heterogeneous and statistically isotropic aquifer whose transmissivity, T, is log-normally distributed with an exponential covariance. The ensemble averages of the second spatial moments of the plume about its center of mass, <S _{i i} (t)>, and the plume centroid covariance, R _{i i} (t) (i=1,2), were simulated for the variance of Y=log T, σ _Y ²=0.1, 0.5 and 1.0 and line sources normal or parallel to μ of three dimensionless lengths, 1, 5, and 10. For σ _Y ²=0.1, all simulated <S _{i i} (t)>−S _{i i} (0) and R _{i i} (t) agree well with the first-order theoretical values, where S _{i i} (0) are the initial values of S _{i i} (t). For σ _Y ²=0.5 and 1.0 and the line sources normal to μ, the simulated longitudinal moments, <S ₁₁(t)>−S ₁₁(0) and R ₁₁(t), agree well with the first-order theoretical results but the simulated transverse moments <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) are significantly larger than the first-order values. For the same two larger values of σ _Y ² but the line sources parallel to μ, the simulated <S ₁₁(t)>−S ₁₁(0) are larger than but the simulated R ₁₁ are smaller than the first-order values, and both simulated <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) stay larger than the first-order values. For a fixed value of σ _Y ², the summations of <S _{i i} (t)>−S _{i i} (0) and R _{i i} , i.e., X _{i i} (i=1,2), remain almost the same no matter what kind of source simulated. The simulated X ₁₁ are in good agreement with the first-order theory but the simulated X ₂₂ are significantly larger than the first-order values. The simulated X ₂₂, however, are in excellent agreement with a previous modeling result and both of them are very close to the values derived using Corrsin's conjecture. It is found that the transverse moments may be significantly underestimated if less accurate hydraulic head solutions are used and that the decreasing of <S ₂₂(t)>−S ₂₂(0) with time or a negative effective dispersivity, defined as , may happen in the case of a line source parallel to μ where σ _Y ² is small.     

Analysis of the impact of injection mode in transport through strongly heterogeneous aquifers

Large-scale advective transport through highly heterogeneous 3D formations is investigated using highly resolved numerical simulations and simple analytic models. Investigations are focused on impacts of two types of contaminant injection on transport through isotropic formations where flow conditions are uniform in the average. Transport is quantified by analyzing breakthrough curves for control planes at various distances from the injection zone. In flux-proportional injection mode local mass in injection zone is proportional to local groundwater flux; this setup models many practical cases such as contaminant injection through wells. In resident concentration mode local concentration in injection zone is constant. Results show that impacts of injection mode on breakthrough curves and their moments are strong and they persist for hundreds of correlation scales. The resident concentration mode leads to a fatter tails of the breakthrough curves, while the peaks are generally underpredicted. For a synthetic porous medium with logconductivity variance of 8, dispersivity computed using resident concentration mode at control plane 100 integral scales away from the injection zone was about 10 times larger than corresponding one for flux-proportional mode. Hence, injection mode impacts on transport through highly heterogeneous formations are strong and they persist for large distances from the injection zone.     

Numerical simulations of transport of non-ergodic solute plumes in heterogeneous aquifers

Transport of non-ergodic solute plumes by steady-state groundwater flow with a uniform mean velocity, μ, were simulated with Monte Carlo approach in a two-dimensional heterogeneous and statistically isotropic aquifer whose transmissivity, T, is log-normally distributed with an exponential covariance. The ensemble averages of the second spatial moments of the plume about its center of mass, <S _{i i} (t)>, and the plume centroid covariance, R _{i i} (t) (i=1,2), were simulated for the variance of Y=log T, σ _Y ²=0.1, 0.5 and 1.0 and line sources normal or parallel to μ of three dimensionless lengths, 1, 5, and 10. For σ _Y ²=0.1, all simulated <S _{i i} (t)>−S _{i i} (0) and R _{i i} (t) agree well with the first-order theoretical values, where S _{i i} (0) are the initial values of S _{i i} (t). For σ _Y ²=0.5 and 1.0 and the line sources normal to μ, the simulated longitudinal moments, <S ₁₁(t)>−S ₁₁(0) and R ₁₁(t), agree well with the first-order theoretical results but the simulated transverse moments <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) are significantly larger than the first-order values. For the same two larger values of σ _Y ² but the line sources parallel to μ, the simulated <S ₁₁(t)>−S ₁₁(0) are larger than but the simulated R ₁₁ are smaller than the first-order values, and both simulated <S ₂₂(t)>−S ₂₂(0) and R ₂₂(t) stay larger than the first-order values. For a fixed value of σ _Y ², the summations of <S _{i i} (t)>−S _{i i} (0) and R _{i i} , i.e., X _{i i} (i=1,2), remain almost the same no matter what kind of source simulated. The simulated X ₁₁ are in good agreement with the first-order theory but the simulated X ₂₂ are significantly larger than the first-order values. The simulated X ₂₂, however, are in excellent agreement with a previous modeling result and both of them are very close to the values derived using Corrsin's conjecture. It is found that the transverse moments may be significantly underestimated if less accurate hydraulic head solutions are used and that the decreasing of <S ₂₂(t)>−S ₂₂(0) with time or a negative effective dispersivity, defined as , may happen in the case of a line source parallel to μ where σ _Y ² is small.     

Stochastic analysis of velocity spatial variability in bounded rectangular heterogeneous aquifers

Transport of inert solutes in two-dimensional bounded heterogeneous porous media is investigated in a stochastic framework. After adopting a first-order approximation of the flow equations, analytical expressions are derived for the velocity covariances. Effects of the boundary conditions and aquifer size upon the statistical moments are analyzed. While the size of the domain is shown to have small influence on the covariances in most cases, the solutions are considerably modified by the boundaries. The results are compared with analytical solutions on infinite domains, and several discrepancies are demonstrated. For example, while the velocity variances on infinite domains are homogeneous, the present results are strongly non-stationary. Finally, the problem is solved numerically by the Monte Carlo simulation method. The results, including the behavior near the boundaries, are shown to be in close agreement with analytical solutions.