Quantifying flow and remediation zone uncertainties for partially opened wells in heterogeneous aquifers

This study presents a numerical first-order spectral model to quantify transient flow and remediation zone uncertainties for partially opened wells in heterogeneous aquifers. Taking advantages of spectral theories in solving unmodeled small-scale variability in hydraulic conductivity ( K), the presented nonstationary spectral method (NSM) can efficiently estimate flow uncertainties, including hydraulic heads and Darcy velocities inrand z-directions in a cylindrical coordinate system. The velocity uncertainties associated with the particle backward tracking algorithm are then used to estimate stochastic remediation zones for scenarios with partially opened well screens. In this study the flow and remediation zone uncertainties obtained by NSM were first compared with those obtained by Monte Carlo simulations (MCS). A layered aquifer with different geometric mean of K and screen locations was then illustrated with the developed NSM. To compare NSM flow and remediation zone uncertainties with those of MCS, three different small-scale K variances and correlation lengths were considered for illustration purpose. The MCS remediation zones for different degrees of heterogeneity were presented with the uncertainty clouds obtained by 200 equally likely MCS realizations. Results of simulations reveal that the first-order NSM solutions agree well with those of MCS for partially opened wells. The flow uncertainties obtained by using NSM and MCS show identically for aquifers with small ln K variances and correlation lengths. Based on the test examples, the remediation zone uncertainties (bandwidths) are not sensitive to the changes of small-scale ln K correlation lengths. However, the increases of remediation zone uncertainties (i.e. the Correspondence to: C.-F. Ni (nichuenfa@geo.ncu.edu.tw) uncertainty bandwidths) are significant with the increases of small-scale lnK variances. The largest displacement uncertainties may have several meters of differences when the ln K variances increase from 0.1 to 1.0. Such conclusions are also valid for the estimations of remediation zones in layered aquifers.


Introduction
Partially opened wells are common elements in groundwater remediation technologies.Such well systems associated with aquifer heterogeneity can create complex flow dynamics around wells and affect significantly the remediation zones for the wells (Zlotnik, 1997).Determination of well remediation zones provides key information to define an area in an aquifer for developments of remediation systems.Due to the complex nature of aquifer heterogeneity and limited capability for data measurements, the incomplete knowledge of aquifer properties, particularly the hydraulic conductivity (K) or the tramsmissivity (T ), will generally lead to the uncertainties of flows and then propagate to the uncertainties of well remediation zones.To quantify the remediation zone uncertainty caused by data limitation and aquifer heterogeneity, a stochastic approach is usually employed (Bair et al., 1991).
Two common approaches, including Monte Carlo simulation (MCS) and so called first-order methods, are generally employed to define stochastic remediation zones (e.g.Varljen and Schafer, 1991;Franzetti and Guadagnini, 1996;Vassolo et al., 1998;Guadagnini and Franzetti, 1999;Riva et al., 1999;Van Leeuwen et al., 1998, 2000;Kunstmann and Kinzelbach, 2000;Feyen et al., 2003a, b; C.-F. Ni et al.: Quantifying flow and remediation zone uncertainties Lessoff and Indelman, 2004;Indelman et al., 2006;Riva et al., 2006;Kunstmann and Kastens, 2006;Guha, 2008).The MCS is conceptually straightforward for determining stochastic remediation zones in heterogeneous aquifers.Using MCS to delineate remediation zones is based on generating a series of equally likely realizations of the K fields that are characterized by the same statistic structure (i.e. the mean value, covariance function and the associated variance and correlation lengths).These K fields are then used as the input for solving groundwater flow equations, resulting in a series of head distributions.Subsequently, the remediation zones for a specified time are defined based on particle tracking algorithms.Collecting the equally like remediation zones then results in a probability distribution of the remediation zone.However, for problems with realistic complexity and sizes, the convergence criteria and the computation effort remain important issues for MCS to quantify flow and remediation zone uncertainty.Discussions regarding to the limitations of MCS have been made in many previous studies (e.g.Guadagnini andNeuman, 1999, Kunstmann andKinzelbach, 2000;Zhang, 2002, Feyen et al., 2003a, b;Ballio and Guadagnini, 2004;Dagan, 2004;Neuman, 2004;Li et al., 2003, 2004a, b, Ni and Li, 2005, 2006).
The first-order methods provide alternatives to the solutions of MCS.Unlike the MCS to resolve small-scale variability directly, most first-order methods focus on solving the transformed functions that link the relationship between input (i.e. the hydraulic conductivity) and output (i.e. the hydraulic head and seepage velocities) variability (Li et al., 2004a, b;Ni andLi, 2005, 2006;Ni et al., 2010).The transform functions such as the statistical moments, Green function, and sensitivity equation, can be solved either analytically or numerically (e.g.Dagan, 1989;Gelhar, 1993;Zhang, 2002;Rubin, 2003;Li et al., 2004a, b).Recent applications of first-order methods have been extended to the determinations of stochastic well capture zones (e.g.Kunstmann and Kinzelbach, 2000;Stauffer et al., 2002Stauffer et al., , 2004;;Lu and Zhang, 2003;Zhang and Lu, 2004;Lessoff and Indelman, 2004;Bakr and Butler, 2005;Riva et al., 2006;Kunstmann and Kastens, 2006;Indelman et al., 2006).Most studies on the subject dealt with depth-averaged two-dimensional problems (Kunstmann and Kinzelbach, 2000;Stauffer et al., 2002Stauffer et al., , 2004;;Lu and Zhang, 2003;Zhang and Lu, 2004;Bakr and Butler, 2005;Riva et al., 2006;Kunstmann and Kastens, 2006).Only a few studies considered problems in threedimensional porous media (Lessoff and Indelman, 2004;Indelman et al., 2006).These proposed three-dimensional solutions are applicable for problems with fully penetrating wells.The efficient closed form solutions in the studies of Lessoff and Indelman (2004) and Indelman et al. (2006) are available for some specified conditions, including infinite domain for boundary conditions, relatively large aquifer thickness compared with vertical correlation scales, and negligible pore scale dispersion in the transport process.
Applications of capture zone delineations can be in aquifers with partially opened wells, where the well screens are relatively small compared with aquifer thickness.Additionally, the degrees of aquifer heterogeneity may cause significant differences in defining remediation zones.These conditions are important especially for the implementation of a remediation well for a contaminant site with either conservative plumes or NAPLs.Motivated by the needs to delineate well remediation zones for such conditions, a numerical profile model in cylindrical coordinate is required for better interpretation of complex flow dynamics around wells.The objectives of this study are (1) to develop a first-order numerical model for stochastic remediation zones in cylindrical coordinate system, and (2) to quantify how and to what degrees the effect of aquifer heterogeneity, well screen locations, and mean flow behavior of layered aquifer on the remediation zone uncertainties.More specifically, a numerical spectral method is employed to predict flow uncertainties for partially opened wells in heterogeneous aquifers.Based on the transient flow uncertainty evaluated by the developed stochastic model, the concept of direct propagation of uncertainties of particle tracks proposed by Kunstman and Kastens (2006) is employed to estimate the uncertainty bandwidth of a remediation zone.To reduce the number of release particles for simulations, the particle backward tracking method will be used for all the simulation examples.This study will first evaluate the accuracy of the developed model for flow and remediation zone uncertainties by employing numerical MCS.A variety of conditions, including the degrees of aquifer heterogeneity, well screen locations, and layered aquifers , were then be considered to quantify the effect of such conditions on the variation of remediation zones uncertainties.

Flow equations
Assuming transient flow in a heterogeneous confined aquifer, the groundwater flow equations in two-dimensional cylindrical coordinate can be formally expressed as (2) subject to the initial and boundary conditions

Mean and perturbation equations
This study considers the variability of K to be solely the source of uncertainty and treats the natural logarithm of hydraulic conductivity (ln K) as stochastic processes.We therefore assume lnK = F + f , where F is the geometric mean of K, and f denotes the perturbations from the mean.The responses of hydraulic head and Darcy velocities to the variation of K are represented by h = H + h , u r = U r + u r , and u z = U z + u z , respectively, where H , U r , and U z are the means and h , u r , and u z represent perturbations.Substituting these stochastic variables (i.e.ln K, h, and u r and u z ) into Eqs.
(1) to (3), neglecting perturbation terms with orders higher than one, and taking expected values of the equations generates the following mean equations (Li and McLaughlin 1991;Gelhar 1993): and the initial and boundary conditions for mean flow are In Eqs. ( 7) to (9), µ r = ∂F (r,z)/∂r, and µ z = ∂F (r,z)/∂z are the gradients of geometric mean K (i.e.K trends) in r-and z-directions, while J r = −∂H (r,z,t)/∂r and J z = −∂H (r,z,t)/∂z are head gradients.Notation Kg = Kg (r,z) is the geometric mean of hydraulic conductivity K.
The mean removed perturbation equations are then given as: subject to the initial and boundary conditions h = 0, ( 16) Note that the assumption that products of fluctuations can be neglected can only be justified when the fluctuation variances of K in aquifers are very small (Dagan, 1989;Gelhar, 1993;Zhang, 2002;Li et al., 2003).Here the perturbations (i.e.Eqs. 13 to 15) describe the linear, nonstationary transformation from f to h to u r and u z .Because the direct solutions of Eqs. ( 13) to ( 15) are unavailable, equations with moment formulas are typically used to analyze the variable correlations forf , h , u r , and u z ( Dagan, 1989;Gelhar, 1993;Zhang, 2002).The Nonstationary Spectral Method (NSM) is a perturbation approach and does not require dependent fluctuations to be stationary.This method differs from other classical perturbation methods primarily in the form of the spectral representation of the output variable fluctuations.The dependent fluctuations are represented as stochastic integrals expanded in terms of sets of unknown complex-valued "transfer functions" such as ψ hf = ψ hf (r,z,k r ,k z ,t) for head fluctuation and ψ u r f = ψ u r f (r,z,k r ,k z ,t) and ψ u z f = ψ u z f (r,z,k r ,k z ,t) for velocity fluctuations, where k r and k z are wave numbers for components r and z, respectively.These fluctuations then have the following Fourier-Stieltjes representation (e.g.Prisetley, 1981;Li andMcLaughlin, 1991, 1995;Li et al., 2004a, b): and C.-F. Ni et al.: Quantifying flow and remediation zone uncertainties where i = (−1) 1/2 and dZ f (k r ,k z ) is the random Fourier increment of f (r,z,t), evaluated at (k r ,k z ) in the spectra domain.The dZ f (k r ,k z )has the following properties Such properties represent that the stochastic process Z f has zero mean and uncorrelation of two different frequencies.The orthogonal increment of Z f will lead to the relationship E[dZ f (k r ,k z ) dZ * f (k r ,k z )] = S ff (k r ,k z )dk r dk z , where S ff (k r ,k z ) is the spectral density function (Priestly, 1981;Gelhar, 1993).The Fourier representation can be viewed as the continuous version of a Fourier series expansion of f .The random Fourier increment at a particular wave number is analogous to the random amplitude of one of the terms in the Fourier integral.The symbols ψ hf , ψ u r f , and ψ u z f are unknown head and velocity transfer functions introduced to account for nonstationary flow transformations.All the transfer functions must be selected such that h , u r , and u z satisfy the governing perturbation equations (i.e.Eqs. 13 to 15).Substituting Eqs. ( 19) to ( 22) into Eqs.( 13) to ( 15) gives the following transfer function equations: ∂ψ hf ∂r and with respect to the initial and boundary conditions Equations ( 23) to (25) are deterministic and complex-valued differential equations.Unlike the classical stationary spectral method, which requires transfer functions to be spatially invariant, the transfer functions introduced here are spatially variants.Three transfer functions ψ hf , ψ u r f , and ψ u z f obtained by solving Eqs. ( 23) to ( 25) can then be used to derive the variances of head and Darcy velocities in the same way as the classical stationary spectral method (e.g.Mizell et al., 1982;Gelhar, 1993): and where the asterisk superscript represents the complex conjugate and S ff (k r ,k z ) is the spectral density function of the log hydraulic conductivity (Priestly, 1981;Gelhar, 1993).Note that the transfer functions obtained from Eqs. ( 23) to ( 25) require a numerical discretization in complex-valued format.
In this study the exponential spectral density function is used for illustrative examples.For specific implementations, a minor revision of the program may be required to involve different spectral density functions.

Determinations of stochastic remediation zones
To determine the uncertainty bandwidth of a remediation zone, this study employs the concept of direct propagation of uncertainties of particle tracks proposed by Kunstman and Kastens (2006).The propagation of particles in the mean flow field U r (r,z,t) and U z (r,z,t) can be formally expressed by: and where t is the time for particle tracking and r(t) and z(t) indicate the location of a particle at a specified time.Because the transient mean velocities U r (r,z,t) and U z (r,z,t) are known values at grid points over the entire modeling area, the position of the particle along its flow path can be calculated by using fourth-order Range-Kutta method (e.g.Zheng and Bennet, 2000;Bakr and Bultler, 2005) where n is the number of total tracking steps.For each particle, a bilinear interpolation algorithm was used in this study to calculate the values of velocity standard deviations in Eqs. ( 35) to (37).When the displacement uncertainties (i.e.σ r , σ rz , and σ ) z of a particle is obtained from Eqs. ( 35) to (37), the uncertainty bandwidth of the particle location can then be approximately calculated by minus and plus one standard deviation (i.e.σ R ) from the mean particle displacement (i.e.R = r 2 + z 2 ), where r and z indicate the mean particle location in r-and z-directions.Estimations of inner and outer bounds of displacement standard deviation for each particle are based on the solutions obtained from NSM velocity uncertainties(i.e.Eqs.35 to 37).In this study, the σ R is estimated based on the following formula (Liu and Zhang, 2003): where symbols W r and W z stand for the referenced point to calculate R for each particle.This study uses the center point of the opened well screen to be the referenced point.
Figure 1 illustrates the concept to calculate the uncertainty bandwidths for a capture zone of a partially opened well.Note that particles are only released along the screen portion of the well.To obtain the uncertainty bandwidths for transition zones (marked by downward diagonal lines in Fig. 1), the locations and displacement uncertainties of the first and the last particles over tracking times need to be recorded.For a specified tracking time, connect the mean trace line of the first particle, particles other than the first and last particles, and the trace line of the last particle, one can obtain the mean remediation zone.Based on the particle locations along the mean remediation zone, the R and σ R values are calculated based on Eq. ( 38) and the associated displacement uncertainties.

Test examples and numerical considerations
Our objective of this study is to develop a spectral first-order method to quantify flow uncertainties and delineate stochastic remediation zone in the cylindrical coordinate system.The illustrative examples here may not cover all the scenarios for partially opened wells, but we aim to present the accuracy and capability of the developed NSM for possible applications to problems with realistic complexity and sizes.Here a synthetic example with modeling area of 80 m by 20 m is employed to illustrate the developed NSM for estimating flow and remediation zone uncertainties in heterogeneous aquifer systems.We assume that a well with partially opened screen is installed in a confined aquifer.Figure 2 shows the conceptual model of the test example.Depending on the problems to be discussed, the locations of the opened well screens are either in the central (8 m to 12 m), upper (14 m to 18 m), or lower (2 m to 6 m) portions of the well.The initial head H 0 for all the simulation scenarios are 10m in modeling areas.The aquifer top and bottom boundaries are specified with no flow boundary conditions.The left boundary is assumed to connect with the well and the portions without well screens are specified with no flow boundary conditions.At the right side of the modeling area, the constant head boundary condition with 10 m is specified.Such boundary condition implies that the distance is sufficient large and the head changes induced by the pumping well are not significant at the boundary.Within the well screen interval, we use a constant head value of 0 m to be the boundary condition for groundwater flowing toward a well screen.Such constant head condition can produce flow rates within the screen interval proportional to the hydraulic conductivity K values along the well screen.Due to the random nature of the K property, one can say that flow in this interval is driven by a source whose strength is a random space function (Severino et al., 2008).Similar conditions were considered by previous investigations (e.g.Dagan, 1989;Indelman, 1996Indelman, , 2002Indelman, , 2004;;Severino et al., 2008 ).The comparison of different boundary types for well locations can be found in the study of Indelman and Dagan (2004).To analyze the effect of aquifer heterogeneity on the predictions of flow and remediation zone uncertainties, the small-scale fluctuation is modeled stochastically by an exponential spectral density function with the ln K variances of 0.1, 0.5, and 1.0, while the correlation lengths in r-direction (λr) are selected to be 1, 5, and 10 m, respectively.For all the simulation scenarios, the correlation lengths in z-direction (λz) are fixed to 1 m.In this study, all the MCS solutions for flow uncertainties are based on 10 000 equally likely realizations of K fields.Such random fields are generated by using the spectral random field generation algorithm (Ni andLi, 2005, 2006).The grid spaces used for NSM simulations are assigned to be 1m in both r-and z-directions, while the grid spaces for MCS simulations are fixed to 0.25 m for better resolution of small-scale K variability.
To conduct stochastic remediation zones for different small-scale K variances and anisotropic scenarios, 20 particles are released along the opened well screens for both NSM and MCS.To model the transient flows for the examples, the time steps are assigned to be variable.In the early simulation time (from 0 to 0.1 day) the time step is 0.001 day.The time step is 0.01 in the period of 0.1 to 1.0 day.After 1.0 day, the time step is fixed to 0.2 day through the rest of the simulation.Note that the NSM requires only one flow and particle tracking simulation for delineating remediation zones, while the MCS requires a number of flows and particle tracking simulations based on different K realizations.For comparison purpose, the results of NSM remediation zones for different ln K variances and anisotropic scenarios will be overlapped on top of the particle clouds created by 200 MCS realizations.

Results and discussion
The first-order method (i.e. the NSM) to delineate transient remediation zones relies on the solutions of velocity variances and cross variance in r-and z-directions in modeling areas.In this study we first assess the accuracy of flow uncertainties estimated by NSM.Then the stochastic remediation zones are delineated based on the flow uncertainties obtained from NSM.The flow and remediation zone uncertainties obtained by using NSM are compared with the corresponding MCS solutions for different small-scale ln K variances and anisotropic scenarios.Based on the verified NSM, the remediation zone uncertainties for different locations of opened  well screens are investigated in a layered and heterogeneous aquifer system.

Simulations of flow uncertainties
On the basis of the screen opened in the central portion of the well, Fig. 3 shows the simulated mean flow pattern (Fig. 3a) and the magnitude profiles (Fig. 3b).To simplify the comparisons of flow and remediation zone uncertainties, here a constant geometric mean K of 1.0 (approximately 2.718 m day −1 ) is assigned for the entire modeling area.The mean head profiles (Fig. 3b) show that the head distributions along center line profile (z = 10 m) reach a steady state in the early simulation time (approximately by 0.2 day).Figures 4 to 6 show, respectively, the selected head and Darcy velocity uncertainties for ln K variance = 0.5 and different anisotropic scenarios by using NSM.In Fig. 4 the head STDs show that the high head uncertainty occurs near the well location, however, the patterns are very different depending on the values of the r-direction correlation lengths.The high values of the head STD increase with the increased r-correlation lengths.With the increase of r-correlation lengths, the hydraulic effect (i.e. the hydraulic gradient) of well screen on the head variations will propagate longer distance from the well locations.Additionally, the increase of r-correlation lengths will lead to significant change of STD patterns near the well location (see Fig. 4b and c).
Figures 5 and 6 show the selected velocity uncertainties (also plotted with STDs) in r-and z-directions for lnK variance = 0.5 and different correlation lengths in r-direction.In  intervals the velocity STDs in r-direction are small for all the anisotropic scenarios.Similar to the solutions of velocity STDs in r-direction, in Fig. 6 the high values of velocity STD in z-direction also located near the well screen and such high velocity STD areas are limited in 2 to 3 m from the well location.On the basis of the algorithm to delineate stochastic remediation zones, the insignificant difference of velocity STDs (both in r-and z-directions) for different anisotropic scenarios may not lead to significant differences of stochastic remediation zones.The results in Fig. 6 also show that the increase of the r-direction correlation length can restrict the propagation of velocity uncertainty in z-direction.Here the small-scale ln K variances are varied from 0.1 to 1.0 and the correlation lengths in r-direction are varied from 1 to 10 m.The results show that NSM solutions for flow uncertainties agree well with those obtained by MCS.In general, the changes of r-direction correlation lengths do not influence much the accuracy of NSM head and velocity STDs (Fig. 7).For isotropic medium, the solutions of velocity STDs for NSM and MCS show identically.Small ln K variance will lead to more accurate estimations of flow uncertainties by using NSM (Fig. 8).Such result is consistent with the assumption of first-order approximation used in the NSM.Note that the velocity uncertainties at boundaries do not reach to zero (Figs.7b, c and 8b, c).This is because of that the values of hydraulic conductivity at boundaries are uncertain.

Simulations of remediation zone uncertainties
Figures 9 and 10 show the delineated stochastic remediation zones by using NSM (shown with lines) and MCS (shown with symbols) for different r-direction correlation lengths and lnK variances.Here the first 200 realizations of MCS solutions are plotted with particle clouds for better presentation.For each realization, a total of 20 particles are released along the opened well screen and are recorded at the end of t = 80 day.The remediation zone clouds in Figs. 9 and 10 are then obtained by collecting all particle locations from the 200 MCS realizations.Figures 9 and 10 show that the NSM remediation zones and the associated uncertainty bandwidths match well with the remediation zone clouds obtained by using MCS.The longer correlation length in r-direction   will lead to a wider uncertainty bandwidth, i.e. the large displacement uncertainty.However, with the small difference of velocity variances, the NSM remediation zones show slight differences for different anisotropic scenarios (Fig. 9). Figure 10 shows the cases with fixed r-correlation length of 5m and different ln K variances.Results show that the increases of ln K variances will lead to increase of uncertainty bandwidths.The largest value of displacement uncertainty (σ R ) for different ln K variance cases will vary from 3 m (ln K variance = 0.1) to 5m (ln K variance = 1.0).
It is worth to mention here the computational efficiency of the developed NSM to delineate the stochastic remediation zones.Based on our workstation with Intel i7 CPU, the computational time to obtain the NSM solution is in the order of minute.However, the computational time for MCS solution based on 10 000 realizations and statistical calculations is in the order of hour.Note that the presented example here is relatively small, which involves a total of 1600 cells for NSM and a total of 25 600 cells for MCS.For most practical problems, the computational domain can be in the order of hundreds of meters to several kilometers.The computational cost for MCS will be very expensive.For such large-scale problems, the developed NSM can provide efficient approximations to quantify flow uncertainties and estimate stochastic remediation zones.2, the aquifer here is divided into two layers with different values of geometric mean K.The geometric mean of K is kept 1.0 for the lower layer (from z = 0 to 10 m).However, we assign a geometric mean K of 3.0 for the upper layer (from z = 10 to 20 m), in which the K value is approximately one order of magnitude greater than the one in the lower layer.Depending on the problems to be discussed, the locations of the well screens are opened either in the central (8 m to 12 m), upper (14 m to 18 m), or lower (2 m to 6 m) portions of the well.
Figure 11 shows the mean head distribution at t = 80 day and the delineated stochastic remediation zones by using NSM.The results in Fig. 11 indicate that the mean flow patterns are influenced by the screen locations and the mean ln K of aquifer layers.Such local flow patterns lead to differences of the patterns of mean remediation zones and the associated uncertainty bandwidths.The well screen in the central and upper portions of the well show similar largest traveling distances of particles in high mean ln K layers but the fronts of the mean and uncertainty bandwidths are slightly different  in patterns (Fig. 11a and c).In Fig. 11b the traveling distances of particles and the patterns of remediation zones in the high mean ln K layer are away from two other scenarios (i.e. the well screen in central and upper portions).The difference is about 30 m based on the solutions at t = 80 day.In the low mean ln K layer we found similar traveling distances for the scenario shown in Fig. 11a.However the patterns of remediation zone uncertainties in low mean ln K layers are different in both the fronts and the trace lines of the first particles.Due to the strong stress created by well screen in the lower portion of the well, the remediation zone in Fig. 11b covers larger area than that in Fig. 11a near the well in the low mean ln K layer.However, the additional area is very small compared with the situations shown in Figs.11a.Note that the abrupt changes of zone uncertainties near the interfaces of the high and low mean ln K layers may be caused by limited particles near the interfaces(Fig.11a and b).In summary, the fronts of remediation zone uncertainties depend highly on the statistical structure of the small-scale K variability, mainly by the variances of ln K variations.The overall patterns of stochastic remediation zones are still controlled by the mean flow behavior.Here such mean flow behavior is generated by different locations of well screens and the mean ln K values in different layers.
We have developed a first-order spectral method to quantify transient flow and remediation zone uncertainties for partially opened wells in heterogeneous aquifers.The developed NSM employs the concept of traditional spectral method and introduce a transfer function in spectral domain to account for aquifer nonstationarity.Based on the velocity uncertainties evaluated by NSM, the concept of direct propagation of uncertainties of particle tracks is then used to calculate stochastic remediation zones for two-dimensional cylindrical coordinate system.In this study, the solutions of developed NSM were first assessed by comparing the solutions of flow uncertainties with the corresponding numerical solutions of MCS.Three ln K variances and anisotropic conditions were considered for the illustrative examples.Based on the velocity uncertainties obtained by using NSM, the first-order stochastic remediation zones were then delineated approximately.The developed model was then employed to estimate remediation zone uncertainties in a layered aquifer under conditions with three screen locations of a well.
The simulation results show that the flow uncertainties obtained by using NSM agree well with the MCS solutions.For aquifers with small ln K variances and correlation lengths, the velocity uncertainties obtained from NSM and MCS show identically.On the basis of velocity uncertainties from first-order solutions, the delineated stochastic remediation zones show reasonably well when compared those first-order remediation zones with the corresponding MCS results.Our illustrative examples involve partially opened well screens and the screen locations are specified with constant head conditions.Under the condition that the screen in the central portion of the well, the velocity uncertainties show slightly differences for different anisotropic scenarios.The NSM remediation zones for different anisotropic scenarios show that the uncertainty bandwidth increases slightly with the increase of correlation lengths in r-direction.However, the increases of remediation zone uncertainties are significant with the increases of small-scale ln K variances.The remediation zone bandwidths may have several meters of differences when the ln K variances increase from 0.1 to 1.0.
The stochastic remediation zones obtained by using NSM in layered aquifer show that the mean flows control the patterns of mean remediation zones and the associated uncertainty bandwidths.The fronts of remediation zone uncertainties depend highly on the statistical structure of small-scale K heterogeneity, mainly by the variances of ln K variations.The location of the well screen plays an important role for the largest length of a remediation zone in the high mean ln K layer.With the well screen entirely or partly in the high mean ln K layer (as shown in Fig. 11a and c), the stochastic remediation zones are similar in high ln K layers.When the well screen is solely opened in the low mean ln K layer (as shown in Fig. 11b), the remediation zone in high mean ln K layer are significantly smaller than those for well screens opened in high mean ln K layers.
In this study we have put our effort on the development of first-order spectral model for transient two-dimensional cylindrical coordinate system.The proposed NSM method has taken the advantages of spectral theories and provided an opportunity to include stochastic theories in practical groundwater modeling problems.The illustrated examples used here for illustrations are synthetic and the hydrogeologic conditions are well defined in advance.For applications of realistic problems, the modeling domain and hydrogeologic conditions can be adjusted to meet conditions on sites.

Fig. 1 .Fig. 1 .
Fig. 1.The concept to calculate the bandwidths of a stochastic remediation zone.

Fig. 2 .
Fig. 2. The conceptual model for illustrated examples in this study.

Fig. 2 .
Fig. 2. The conceptual model for illustrated examples in this study.

Fig. 3 .
Fig. 3.The mean head distribution for the screen opened in the central portion of the well: (a) the r and z two-dimensional head distribution for t = 80 day, and (b) the snapshots of head centerline profiles along z=10m for different simulation times.

Fig. 3 .
Fig. 3.The mean head distribution for the screen opened in the central portion of the well: (a) the r and z two-dimensional head distribution for t = 80 day, and (b) the snapshots of head centerline profiles along z = 10 m for different simulation time.

Fig. 7 .Fig. 7 .Fig. 8 .
Fig. 7.The center line profiles of flow uncertainties that are obtained by using first-order NSM (lines) and MCS (symbols) at t= 80 day for fixed lnK variance of 0.5 and different r-direction correlation lengths.

Fig. 8 .Fig. 9 .
Fig. 8.The center line profiles of flow uncertainties that are obtained by using NSM (lines) and MCS (symbols) at t = 80 day for fixed r-direction correlation length of 5 m and different ln K variances.

Fig. 11 .
Fig. 11.Stochastic remediation zones obtained by using first-order NSM at t =80 days(flooded contours: mean head distribution; solid lines: mean remediation zone; dashed lines: mean-R, and dash-dotted line: mean+R ) for fixed r-direction correlation length of 5m and lnK variance of 0.5 in a layered aquifer: (a) the screen opened in the central portion of the well (8 to 12m), (b) the screen opened in the lower portion of the well (2 to 6m), and (c) the screen opened in the upper portion of the well (14 to 18m).

Fig. 11 .
Fig. 11.Stochastic remediation zones obtained by using first-order NSM at t = 80 day (flooded contours: mean head distribution; solid lines: mean remediation zone; dashed lines: mean-σ R , and dashdotted line: mean+σ R ) for fixed r-direction correlation length of 5 m and lnK variance of 0.5 in a layered aquifer: (a) the screen opened in the central portion of the well (8 to 12 m), (b) the screen opened in the lower portion of the well (2 to 6 m), and (c) the screen opened in the upper portion of the well (14 to 18 m).
is the hydraulic conductivity, S s = S s (r,z) is the specific storage, and u r = u r (r,z,t) and u z = u z (r,z,t) [L/T] are Darcy velocities in r-and z-directions for the aquifer system.H 0 is the initial head in the simulation domain , H D is the prescribed head on Dirichlet boundary D , and Hydrol.Earth Syst.Sci., 15, 2291-2301, 2011 www.hydrol-earth-syst-sci.net/15/2291/2011/ where h = h(r,z,t)[L] is the hydraulic head, K = K(r,z)[L/T]