VZJ Download to Citation Manager
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Nolan, B. T.
Right arrow Articles by Kauffman, L. J.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Nolan, B. T.
Right arrow Articles by Kauffman, L. J.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Nolan, B. T.
Right arrow Articles by Kauffman, L. J.
Related Collections
Right arrow Flow
Right arrow Spatial Distribution
Right arrow Ground Water Quality
Published in Vadose Zone Journal 2:677-691 (2003)
© 2003 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH PAPERS

Spatial Variability of Groundwater Recharge and its Effect on Shallow Groundwater Quality in Southern New Jersey

Bernard T. Nolan*,a, Arthur L. Baehrb and Leon J. Kauffmanb

a U.S. Geological Survey, 413 National Center, Reston, VA 20192
b U.S. Geological Survey, West Trenton, NJ
* Corresponding author (btnolan{at}usgs.gov).

1 Use of brand names in this paper does not constitute endorsement by the USGS. Back


Received 27 November 2002.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Point estimates of groundwater recharge at 48 sediment-coring locations vary substantially (-18.5–1840 cm yr-1) in a 930-km2 area of southern New Jersey. Darcian estimates of steady, long-term recharge made at depth in the unsaturated zone were estimated using pedotransfer functions of soil texture and interpolated (mapped) with nonparametric methods to assess aquifer vulnerability in the area. The probability of exceeding the median recharge (29.1 cm yr-1) is low in the southwestern and northeastern portions of the study area and high in the eastern and southeastern portions. Estimated recharge is inversely related to measured percentage clay and positively related to the percentage of well-drained soils near wells. Spatial patterns of recharge estimates, exceedance probabilities, and clay content indicate that sediment texture controls recharge in the study area. Relations with land elevation and a topographic wetness index were statistically insignificant. Nitrate concentration and atrazine (6-chloro-N2-ethyl-N4-isopropyl-1,3,5-triazine-2,4-diamine) percentage detection in samples of shallow groundwater (typically <10 m) are higher for low recharge sites (<=29.1 cm yr-1) than for high recharge sites (>29.1 cm yr-1) in agricultural and urban areas. Differences between high and low recharge sites in these areas are highly significant for NO3 concentration, but not for atrazine concentration.

Abbreviations: CCDF, conditional cumulative distribution function • DEM, digital elevation model • DO, dissolved oxygen • GIS, geographic information system • IGF, Indicative Goodness of Fit • IK, indicator kriging • KED, Kriging with external drift • MLR, multiple linear regression • PTF, pedotransfer function


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
NONPOINT-SOURCE contamination is considered the single greatest threat to water quality (Corwin et al., 1997). Preventing contamination of groundwater is crucial in areas where it is a major source of public and domestic supply. Knowing where an aquifer is vulnerable to surface-derived contaminants would help managers prioritize scarce resources for alternative management practices, monitoring, and cleanup.

Aquifer vulnerability studies at large spatial scales have used index methods, such as DRASTIC and SEEPAGE (Navulur and Engel, 1996), or overlays made with geographic information systems (GISs) (Nolan et al., 1997). Both DRASTIC and SEEPAGE underestimated contamination potential by describing areas with high NO3 concentration as low risk (Navulur and Engel, 1996). Index and overlay methods provide only limited understanding of processes controlling the transport of water and chemicals in the unsaturated zone. Alternatively, deterministic models can simulate water and chemical fluxes, but the spatial variability of sediment properties at field scales and above limits accuracy and imposes large uncertainty on model predictions.

In the current study, we used a combined deterministic–geostatistical approach to assess aquifer vulnerability in the Glassboro, NJ area (Fig. 1). The Darcian method was used to estimate groundwater recharge from water-retention parameters and unsaturated hydraulic conductivity on the basis of sediment texture and moisture content data obtained near the water table at 48 locations in the study area. The recharge estimates were geostatistically analyzed to evaluate the spatial variability of measured sediment properties, to map recharge with respect to land use, and to derive statistical distributions of recharge at specific locations in the study area. The recharge estimates were compared with soils and topographic data to determine whether recharge could be accurately predicted from landscape characteristics. Finally, recharge estimates were compared with concentrations of NO3 and atrazine to evaluate potential effects on the quality of shallow, recently recharged groundwater. In this study, "shallow groundwater" refers to depths <10 m. The depth of the screened interval below water in observation wells in the area is about 3 m. The objectives of the study were to



View larger version (110K):
[in this window]
[in a new window]
  		Fig. 1. Glassboro, NJ study area and outcrop of the Kirkwood-Cohansey aquifer.

 

   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Description of Study Area
The study area (Fig. 1) comprises about 930 km2 within the Coastal Plain Physiographic province of southern New Jersey. Population in the area has increased from about 50 000 people in 1940 to about 250 000 in 2000. Groundwater withdrawals from the surficial, Kirkwood-Cohansey aquifer system, which consists of highly permeable unconsolidated sands and gravels, have recently increased to meet the growing demand for drinking water. As of 1986, the Glassboro area comprised 21% urban land, 26% agricultural land, and 39% undeveloped land (Stackelberg et al., 1997).

The outcrop of the Kirkwood Formation, a confining unit about 30 m thick, underlies the aquifer and forms the northwest boundary of the study area. Aquifer thickness increases to about 75 m at the southeastern boundary (Zapecza, 1989). Unsaturated zone sediment in the study area consists mainly of the Cohansey Sand, which was deposited during the Miocene Age on inner shelf, nearshore, and beach areas during slow retreat of the sea. Sediments in the Cohansey Sand generally are coarser at shallower depths, which is consistent with similarly deposited formations in the New Jersey Coastal Plain (Zapecza, 1989). The Bridgeton Formation overlies the Cohansey Sand in parts of the study area. This formation is generally less than 3 m thick in the study area and consists of coarse, pebbly orange sands deposited during late Tertiary and Quaternary periods.

Groundwater in the study area was sampled with shallow observation wells installed about 3 m below the water table (Stackelberg et al., 1997). Public supply wells and deeper observation wells were sampled in a later study (Stackelberg et al., 2000) that described the migration of compounds to deeper parts of the aquifer, and related compound occurrence to current and historical land-use patterns. Many compounds, including NO3, atrazine, desethylatrazine, simazine, metolachlor, prometon, dieldrin, chloroform, and methyl tert-butyl ether, were detected in these surveys. Groundwater quality in the study area depends on contributing land area, groundwater age, and well type (shallow monitoring, moderate-depth monitoring, or public supply) and related construction characteristics (Kauffman et al., 2001).

The median concentration of NO3 in shallow groundwater beneath agricultural land in the study area is high (13 mg L-1) (Stackelberg et al., 1997). Median NO3 concentration is 2.6 and 3.5 mg L-1 in new and old urban areas, respectively, and is 0.07 mg L-1 in undeveloped areas. The USEPA has established a Maximum Contaminant Level of 10 mg L-1 NO3 as N (USEPA, 1995). Ingestion of NO3 by infants can cause low O2 levels in the blood, a potentially fatal condition known as methemoglobinemia or "blue baby" syndrome (Spalding and Exner, 1993). Additionally, NO3 concentrations of 4 mg L-1 or more have been associated with increased risk of non-Hodgkin's lymphoma in Nebraska (Ward et al., 1996).

Groundwater recharge was estimated at 47 locations where sediment texture and moisture content data were collected during installation of the shallow groundwater observation wells, and at one additional sediment-coring site, for a total of 48 estimates (Baehr et al., 2003). In the current study, we analyzed the spatial variability of the recharge estimates to better understand relations between recharge, spatial scale, landscape characteristics, and the quality of shallow groundwater in the Glassboro area.

Point Recharge Estimates
Groundwater recharge is the rate at which infiltrating water moves across the water table. Greater amounts of infiltrating water might adversely affect water quality by transporting more chemicals to shallow groundwater. Conversely, greater amounts of infiltrating water might attenuate contaminant concentrations through dilution, for a given chemical loading.

Sediment cores were collected at 48 locations during the summer and fall of 1996 to facilitate point estimates of recharge. Well locations were randomly selected within major land-use categories representing agricultural, new urban, old urban, and undeveloped lands. Samples representing prominent sediment horizons were analyzed by optical diffraction to determine percentages of sand, silt, and clay. These percentages were used with Rosetta, a program consisting of hierarchical pedotransfer functions (PTFs) (Schaap et al., 2001), to derive recharge estimates. Pedotransfer functions are commonly used in regional simulations of unsaturated-zone contaminant transport (Görres and Gold, 1996; Petach et al., 1991; Soutter and Pannatier, 1996) because of the difficulty and cost of directly measuring soil hydraulic properties at numerous locations at large spatial scales. Rosetta uses neural network analysis to extract as much information as possible from measured soils data (Schaap et al., 1998). Compared with regression, neural networks require no a priori model concept and optimally relate input and output data using an iterative calibration procedure. Rosetta input consists of sediment texture (percentage sand, silt, and clay) with or without other properties such as bulk density, and the output comprises parameters used in water-retention functions and calculations of partially saturated hydraulic conductivity. Although the PTFs used in this research are calibrated to soil samples outside of the Glassboro study area, textural distributions indicate that the samples on which Rosetta is based are predominantly coarse grained (Schaap et al., 2001), which is consistent with the sandy sediments composing the Kirkwood-Cohansey aquifer system.

Recharge was calculated at each well location as described by Baehr et al. (2003), using Darcy's Law for unsaturated flow:

[1]
where qi is water flux at depth i (cm yr-1), Ki({theta}) is unsaturated hydraulic conductivity at depth i (cm yr-1), {psi}i is matric potential (water pressure or head) at depth i (cm), and zi is depth (cm), which is positive upwards. Although the water fluxes could not be calculated precisely at the water table, they were calculated for layers deep in the unsaturated zone and represent recharge in this study.

The matric potential is defined as a function of the relative saturation (van Genuchten, 1980):

[2]
where {alpha} (1/cm) and n (dimensionless) are curve-fitting parameters. The parameter 1/{alpha} is the air-entry pressure, and the relative saturation (Se) is defined as

[3]
where {theta}r and {theta}s are residual and saturated moisture contents, respectively, and {theta} is measured moisture content at a given sampling depth.

Equation [3] is used with the pore size distribution model of Mualem (1976) to yield the van Genuchten–Mualem model (van Genuchten, 1980) for K({theta}):

[4]
where K0 (the fitted matching point at saturation, cm yr-1) and L (dimensionless) are curve-fitting parameters. The parameters {theta}r, {theta}s, {alpha}, n, Ko, and L were obtained using the Rosetta model for each sediment layer for which texture data were collected (Baehr et al., 2003).

The hydraulic gradient ({psi}i-1 - {psi}i)/( zi-1 - zi) shown in Eq. [1] was calculated based on changes in measured moisture content within a single sediment layer at depth (Baehr et al., 2003). Because sediment texture essentially is uniform within a layer, changes in moisture content reflect the hydraulic gradient and not textural variability. Gravity flow (where K indicates q) was not assumed because uniformity in {psi} is less likely in a layered medium (Nimmo et al., 2002).

The recharge estimates are a snapshot in time and do not incorporate potential seasonal changes or year to year climatic variability. However, temporal variability at depth is dampened by physical averaging over several years such that a single measurement represents constant flux (Nimmo et al., 2002). The median depth of the recharge calculations is 2.4 m and median unsaturated zone thickness is 3.8 m (Baehr et al., 2003). The 2.4-m depth is greater than depths (0.6 m or less) at which Bayless (2000) observed short-term variations in moisture content and soil tension in glacial till in central Indiana. They measured soil tension profiles to assess hydraulic gradients, and soil moisture following significant rainfall. Soil tension was fairly stable at depths of 0.9 to 1.5 m during the period of data collection (May–September 1994), indicating reduced effects from short-term wetting and drying periods. However, low-frequency (e.g., decadal) variations in moisture content can propagate to deeper depths, and these were not considered in the current study.

Water flux estimates from PTFs have considerable uncertainty because K({theta}) and {psi}({theta}) are nonlinear functions of moisture content. We evaluated uncertainty in several ways. First, for selected sites we compared PTF estimates of recharge with those based on direct measurement of K({theta}) and {psi}({theta}) by the steady-state centrifuge method (Nimmo et al., 1987; Conca and Wright, 1998). We also compared the median of the PTF estimates with results from water-budget calculations in the region. Second, Rosetta provides a bootstrap estimate of the standard deviation (uncertainty) of each parameter used in the water-flux calculations. Third, we used indicator kriging to develop probabilistic models of recharge uncertainty at unsampled locations. Finally, kriging provides a map of standard deviations that indicates the reliability of predictions.

The Darcian method determines water flux in the sediment matrix and does not consider preferential pathways caused by, for example, desiccation cracks and root holes. Estimates of K({theta}) and {psi}({theta}), however, are functions of measured moisture content and can reflect focused recharge in topographically low-lying areas. Spatial concentration of recharge can occur in depressions and channels where high water content increases K({theta}) and promotes rapid water flux (Nimmo et al., 2003). Extremes of water flux within the sediment matrix are a type of preferential flow.

Mapping and Spatial Analysis Techniques
Kriging, an interpolation method that considers the spatial dependence of samples, was used to predict recharge at unsampled locations. We emphasized kriging in this study because (i) it is an exact interpolator and, therefore, provides a realistic map (kriging returns the observed value at locations with measured data) and (ii) kriging quantifies prediction uncertainty through the kriging standard deviation. In contrast, other least-squares algorithms (such as multiple linear regression, MLR) yield smoothed trend surfaces that often do not fit environmental data very well because the trend surfaces are pulled toward extreme values. However, MLR is potentially useful in defining relations between recharge and landscape variables and can augment kriging predictions. Multiple linear regression in a kriging context is discussed below after indicator kriging.

Indicator kriging (Isaaks and Srivastava, 1989; Goovaerts, 1997; Deutsch and Journel, 1998) was used to reduce the influence of extreme values of recharge estimates on the variogram and the kriged prediction. The recharge estimates are highly skewed (skewness coefficient = 4.1) and yield erratic variograms that are difficult to interpret. Indicator kriging is a nonparametric geostatistical technique that involves transformation of a variable to a binary response (0, 1) based on cut-off values such as quartiles or quintiles of the statistical distribution of the variable. The binary, transformed recharge data are referred to in this paper as "indicators." The indicator transform is consistent with a categorical description of recharge, used in this study to mitigate the uncertainty of the PTF estimates (see Uncertainty of Recharge Estimates below). Ordinary kriging is performed on the indicators corresponding to each cut-off using a unique variogram for each, and probabilities of exceeding a specified cut-off value are generated for each class of indicators. Median approximation indicator kriging simplifies the calculations by using the same variogram (e.g., corresponding to the median cut-off) for all classes of indicators. The most representative variogram is used, even if based on indicators from a cut-off other than the median (Isaaks and Srivastava, 1989). In the remainder of this paper, ordinary kriging of recharge indicators by the median approximation method is referred to simply as "indicator kriging" (IK).

We calculated experimental variograms of recharge indicators from covariance functions using the following relation (Isaaks and Srivastava, 1989):

[5]
where {gamma}*(h) is the experimental variogram for separation distance (lag vector) h, {sigma}2 is sample variance, and C(h) is the covariance function.

The covariance function is given by (Deutsch and Journel, 1998):

[6]
where N(h) is the number of sample pairs separated by lag vector h, xi is the recharge value at the beginning ("tail") of the lag vector, yi is the recharge value at end ("head") of the lag vector, m-h is the mean of the tail values = {Sigma}N1xi, and m+h is the mean of the head values = {Sigma}N1yi. The covariance function is resistant to erratic values because, unlike the traditional variogram, it accounts for the sample means of lag distance categories (Isaaks and Srivastava, 1989).

We used the spherical variogram model (Deutsch and Journel, 1998) to fit the experimental variograms of recharge indicators:

[7]
where {gamma}(h) is the variogram (probability units)2, c is the sill (probability units), h is separation distance (lag vector) between pairs of samples (m), and a is the variogram range (m).

Model parameters used in IK were obtained by fitting Eq. [7] to the experimental variograms of recharge indicators. The range of the variogram (a) represents the limit of spatial dependence, within which samples are spatially autocorrelated, and beyond which samples are independent. The sill (c) represents the maximum value of the variogram and is approximated by the sample variance. We performed variogram analysis using Variowin (Pannatier, 1996) and kriged the recharge indicators in Surfer (Golden Software, Inc., 1999).1

Kriging with external drift (KED) was attempted as an alternative to indicator kriging to incorporate information from densely sampled secondary attributes consisting of landscape variables. Kriging with external drift is appropriate when the local mean of a process (e.g., recharge) varies within the kriging search neighborhood (Goovaerts, 1997). The KED prediction at each grid node is obtained by adding a regression estimate to the kriged estimate of the regression residual. Multiple linear regression was performed to derive the drift component of KED and also was used in a separate exercise to identify variables that significantly influence recharge in the study area. Prior researchers observed that recharge was focused beneath topographic depressions at both upland and lowland locations in an agricultural field (Delin et al., 2000); and in the current study, recharge values >36.5 cm yr-1 in the east and southeast generally occur in low-lying areas with low clay content (Fig. 2).



View larger version (61K):
[in this window]
[in a new window]
  		Fig. 2. Soil coring locations, recharge estimates, and percentage clay in the Glassboro, NJ study area. Percentage clay is from the USDA Soil Survey Geographic database.

 
Landscape variables evaluated in KED include land elevation, depth to groundwater, soil survey attributes, and a topographic wetness index given by ln(a/tanß), where a is the upslope area per unit contour length and tanß is the slope gradient along which drainage occurs. The topographic wetness index was computed from 1:250 000-scale digital elevation models (DEMs) with a resolution of about 90 m (Wolock and McCabe, 1995). The parameter a is based on the total area draining into each DEM grid cell, which was calculated using a single flow direction algorithm described by Wolock and McCabe (1995). We anticipated that topographic wetness index would be negatively correlated with recharge. High values of the index represent near-stream saturated (discharge) areas associated with a rise in water-table elevation. The rise in the water table occurs when infiltrating water migrates laterally in low-lying areas and can result in streamflow generation.

Lastly, we attempted cokriging to exploit potential relations between the primary variable (recharge) and densely sampled secondary variables, such as drainage characteristic from soil survey data. However, cokriging requires reasonable correlation (typically 0.5 or greater) between the primary and secondary variables. Cokriging was unsuitable with these data because Spearman correlations between recharge and soil survey properties are 0.2 or less.

Water Quality Relations
High recharge rate might adversely affect groundwater quality by transporting more chemicals to the water table. Conversely, for a given chemical load, high recharge rate might reduce contaminant concentration through dilution. To test these opposing hypotheses, median NO3 concentration and atrazine percentage detection were compared for high and low recharge categories representing agricultural, new urban, old urban, and undeveloped lands. Recharge values >29.1 cm yr-1 (the median of PTF-derived recharge estimates) were categorized as "high" and values <=29.1 cm yr-1 were categorized as "low." We attempted to account for differences in chemical loading by stratifying data by land-use category. Reliable data on chemical loading for specific locations within the study area are unavailable. State-wide data (New Jersey Department of Environmental Protection, 1996) were compiled in a GIS of the study area to assign wells to land-use categories. Water samples generally were collected from observation wells within two months after well installation (Stackelberg et al., 1997). Nitrite plus nitrate concentration (based on elemental N) is referred to in this paper as "nitrate" because nitrite contribution to nitrite plus nitrate in groundwaters of the area generally is <0.02 mg L-1. Atrazine concentration was determined by solid phase extraction with capillary column gas chromatography–mass spectrometry (Zaugg et al., 1995). Medians, which are resistant to the effects of outliers typical of skewed data sets (Helsel and Hirsch, 1992), were used as the measure of central tendency of NO3 data grouped by land-use (agricultural, new urban, old urban, and undeveloped) and recharge categories. Censored values were set to one-half the detection limit before calculating medians of categories. Atrazine concentration in groundwater commonly is below the detection limit of 0.001 µg L-1 and therefore less amenable to comparisons involving median concentrations; accordingly, percentage detection of atrazine was computed within categories for initial water quality comparisons.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Uncertainty of Recharge Estimates
Recharge estimates in the study area range from -18.5 to 1840 cm yr-1, and median recharge is 29.1 cm yr-1 (Table 1). The sign convention in Eq. [1] was reversed so that positive values indicate downward movement of water to the aquifer, and negative values indicate upward movement of water in response to evapotranspiration. The wide range of values reflects the considerable spatial variability of sediment properties on which the estimates are based, and also the uncertainty of the flow calculations. However, the median value compares favorably with regional estimates of 33.1 to 49.3 cm yr-1 obtained by water-budget calculations (Baehr et al., 2003).


View this table:
[in this window]
[in a new window]
  		Table 1. Summary statistics of recharge estimates for the Glassboro, NJ study area.{dagger}

 
The uncertainty of recharge (q) was evaluated by comparing the PTF estimates with calculations based on measured values of K({theta}) and {psi}({theta}) obtained by the steady-state centrifuge method for selected samples (Baehr et al., 2003). Although the PTF estimates do not precisely match the latter, they are well within an order of magnitude (Table 2). Additionally, both methods predict low q for four sites and high q for two sites, indicating that conversion of the PTF estimates to categorical variables ("low" and "high" in Table 2) is appropriate. Because the value of the categorical variable is unaffected by the method of recharge estimation, it helps compensate for the uncertainty of the PTF estimates.


View this table:
[in this window]
[in a new window]
  		Table 2. Comparison of recharge estimates based on the Rosetta model and on measured values of hydraulic conductivity and matric potential from the steady-state centrifuge method (from Baehr et al., 2003).

 
The standard deviation of bootstrap estimates in Rosetta provides an additional measure of the uncertainty inherent in the recharge calculations. Mean values of Rosetta output parameters ± 1 standard deviation are shown in Table 3 for the same sites that had measured values of K({theta}) and {psi}({theta}) by the centrifuge method. The range of Ko standard deviations [0.239–0.360 log(cm d-1)] is somewhat higher than the range of 0.079 to 0.11 log(cm d-1) reported for saturated hydraulic conductivity (Ks) in predominantly sandy soils (Schaap et al., 2001). We emphasized Ko in this study because it yields better predictions of unsaturated hydraulic conductivity when suctions are at least a few centimeters (Schaap et al., 2001). To enhance comparison with prior work, we also calculated bootstrap standard deviations for Ks [0.072–0.146 log(cm d-1)], which compare favorably with the above values.


View this table:
[in this window]
[in a new window]
  		Table 3. Uncertainty of Rosetta output parameters used to calculate matric potential and partially saturated hydraulic conductivity, for selected sites.

 
Recharge values in excess of the 60th percentile (36.5 cm yr-1) in the east and southeast generally coincide with areas having low percentage clay in the upper 1.8 m of soil (Fig. 2). Percentage clay in Fig. 2 is based on Soil Survey Geographic (SSURGO) data (USDA-NRCS, 1995). Areas with <12.5% clay predominantly occur near major drainages in the eastern and south-central portions of the study area. High-recharge sites in these areas generally are located in coarse-textured soils near streams. Streamflow in the study area is derived in part from young, recently recharged groundwater (<25 yr) entering near the edges of the stream (Kauffman et al., 2001). Several values of recharge exceed the annual precipitation rate of 109 cm yr-1. Very high values of recharge (>177 cm yr-1) that cluster in the southeast might indicate focused recharge in local, topographically low-lying areas with coarse-textured sediments. Because of the extreme variability in estimated recharge, categorical descriptions of recharge ("low" and "high" as described above) are emphasized in this study and the transformed estimates are scaled by indicator kriging, a nonparametric method of spatial interpolation.

Application of Geostatistical Techniques
Indicator kriging was performed at the median cut-off by transforming the recharge estimates to 0 if less than or equal to the median (29.1 cm yr-1) and 1 if greater than the median. Additional indicator data sets were created and analyzed based on the quintiles of the distribution of recharge values to evaluate variograms for several cut-offs. A spherical variogram model corresponding to 60th percentile cut-off indicators was used to krig all of the indicator data sets (median cut-off and 20th, 40th, 60th, and 80th percentile cut-offs) because it fit the experimental data well and best represented anisotropic patterns in the data. The models were manually fitted to the experimental variograms, and fit was assessed with the Indicative Goodness of Fit (IGF) statistic (Pannatier, 1996); values approaching zero indicate good fit. The IGF statistic measures how well the variogram model fits the experimental data based on the number of lag classes, number of sample pairs, mean and maximum separation distances, and the differences between observed and fitted variogram values. Figure 3a indicates that omnidirectional variograms corresponding to the 20th and 80th percentile cut-offs are somewhat lower than for the 40th to 60th percentile cut-offs. Experimental variograms of all cut-offs, however, indicate similar range (approximately 5000–7000 m).



View larger version (29K):
[in this window]
[in a new window]
  		Fig. 3. (a) Omnidirectional, experimental variograms of binary recharge indicators for all cut-offs and (b) anisotropic, spherical variogram models corresponding to the 60th percentile cut-off (36.5 cm yr-1).

 
The variogram represents variation between samples as a function of separation distance and indicates the degree of spatial dependence (autocorrelation) between samples. The omnidirectional experimental variograms of 60th percentile recharge indicators increase to about 6000 m and then level off, indicating spatial dependence up to this distance (Fig. 3a). Beyond about 6000 m, the recharge indicators are spatially uncorrelated. There is only one data point near the origin, which creates uncertainty about the variogram at small separation distances. This point is based on 66 sample pairs, however, which is reasonable for variogram analysis. Use of smaller classes of h to yield more points on the experimental variogram was not feasible with the overall number of samples (48). More samples closer than 6000 m would improve definition of variogram behavior near the origin.

Variograms of 60th percentile recharge indicators were analyzed in all directions (0°, 20°, and so on up to 160°) to identify potential anisotropy. Zonal anisotropy causes the variogram range to differ depending on direction, precluding use of isotropic models. Variogram models obtained with Eq. [7] indicate maximum and minimum ranges in the 40° and 130° directions, respectively. The fit of the models is reasonable because they pass through the experimental values near the origin, which helps define the nugget, and the IGF statistic (0.052) is reasonably close to zero. The 40° direction is about parallel with the current shoreline (Fig. 1). Both variogram models have nugget = 0.13 and sill = 0.24, but the range of the 40° variogram is 7600 m, and the range of the 130° variogram is 4560 m (Fig. 3b). The 40° direction is referred to as the axis of major anisotropy, and the 130° direction is designated the axis of minor anisotropy. The nugget (y-intercept) represents sampling error and small-scale variability within the smallest sampling interval, and the sill represents the maximum variation between samples.

The 40° variogram range (7600 m) appears related to cyclical patterns of sediment texture in the study area, rather than patterns of recharge and discharge, because distances separating water bodies are substantially less. Soils with <12.5% clay indicate broad, low-lying areas in the east and south (Fig. 2). The width of these areas is about 8000 m in the northeast–southwest direction. Similarly, the northeast–southwest distance separating areas in the west with clay content in excess of 16.4% is about 8000 m. Recharge estimates in this study are based in part on measured sediment texture, which is correlated with the SSURGO data mapped in Fig. 2. The Spearman correlation between surficial, measured clay content and SSURGO percentage clay at sediment coring locations is 0.453 (p < 0.001).

The axis of major anisotropy (40°) is consistent with the depositional history of the Kirkwood-Cohansey aquifer system, which generally is aligned in the northeast–southwest direction in the study area (Fig. 1). The Cohansey Sand composes most of the unsaturated zone in the study area and was deposited in inner shelf, nearshore, and beach areas during slow retreat of the sea during the Miocene Age (Zapecza, 1989). The depositional history causes sediment properties to be more alike in the northeast–southwest direction than in the perpendicular direction because it parallels the ancestral shoreline. This similarity is manifested as increased spatial continuity of estimated recharge in the 40° or northeasterly direction.

Figure 4 indicates the probability of exceeding median recharge (29.1 cm yr-1) in the study area. Indicator kriging was performed on median recharge indicators using the anisotropic variogram model corresponding to the 60th percentile recharge indicators. An elliptical search neighborhood with a maximum distance of about 19 000 m was used to incorporate the zonal anisotropy. Probabilities are low in the southwestern and northeastern portions of the study area and are high in the eastern and southeastern portions.



View larger version (82K):
[in this window]
[in a new window]
  		Fig. 4. Predicted probability of exceeding median recharge (29.1 cm yr-1) in the Glassboro, NJ study area. Conditional cumulative distribution functions (CCDFs) in Fig. 5 were developed for the two locations shown.

 


View larger version (27K):
[in this window]
[in a new window]
  		Fig. 5. Conditional cumulative distribution functions (CCDFs) corresponding to low and high recharge locations in the Glassboro, NJ study area. The recharge locations corresponding to the CCDFs are shown in Fig. 4.

 
The IK map identifies an area where recharge is expected to be less than the regional median (Baehr et al., 2003). The low probabilities in the southwest coincide with the Cohansey River Basin, where average annual recharge is estimated to be 37.1 cm yr-1 (Charles et al., 2001). This value is the 25th percentile of the range of estimates (33.1–49.3 cm yr-1) from water-budget calculations for the region. High probabilities in the east and southeast generally coincide with areas with <12.5% clay (Fig. 2). The spatial patterns of recharge estimates, IK probabilities, and SSURGO clay content, and direct correlation between measured and mapped clay percentages indicate that sediment texture controls recharge in the study area.

In this study, variograms were developed for noninundated land areas, and the variogram range represents large-scale changes in sediment texture. Therefore, IK maps shown here do not apply to discharge areas (i.e., where the water table is at or above the land surface).

Assembling the IK-predicted probabilities for all cutoffs at a single location yields a least-squares estimate of the conditional cumulative distribution function (CCDF) of recharge. The CCDF is a probabilistic model of uncertainty (Deutsch and Journel, 1998) that indicates exceedance probabilities of recharge at unsampled locations. For example, the predicted probability of exceeding the 20th percentile cutoff (0.5 cm yr-1) is 0.51 at a specific IK grid node on the recharge map, the probability of exceeding the 40th percentile cutoff (16.8 cm yr-1) is 0.39, and so on. Indicator kriging is nonparametric in the sense that the CCDF is estimated directly from the data and not from the mean and variance of a statistical distribution of recharge values at a particular location.

Kriging estimates of indicators corresponding to recharge quintiles (20th, 40th, 60th, and 80th percentiles) were assembled into CCDFs at selected IK grid nodes representing high and low recharge areas (Fig. 5). The recharge quintiles represent cutoffs ranging from 0.5 to 177 cm yr-1. The resulting CCDFs represent the statistical distributions of recharge probabilities at unsampled locations and are shaped differently depending on whether they are in high or low areas of recharge as predicted by IK. The CCDF for the high recharge areas generally indicates higher exceedance probabilities at a given recharge rate. For example, the probability of exceeding the 60th percentile value (36.5 cm yr-1) is 0.63 in the high recharge area and is 0.22 in the low recharge area.

Potential Causes of Recharge Variability
We evaluated KED as an alternative to IK to incorporate information from densely sampled secondary attributes. SSURGO attributes screened in regression (the drift component of KED) include individual soil drainage categories, such as "moderately well drained" and "poorly drained," and lumped SSURGO categories representing soil drainage characteristic and soil hydrologic group (Table 4). The latter variable describes the infiltration characteristics of soils. Measured sediment texture was evaluated in a separate MLR exercise but not in KED, which requires secondary information throughout the study area (i.e., at kriging grid nodes).


View this table:
[in this window]
[in a new window]
  		Table 4. Explanatory variables evaluated in kriging with external drift and multiple linear regression.

 
Percentage of moderately drained soils is highly significant (F-ratio p = 0.003) in the KED model, and the poorly drained soils variable approaches statistical significance at the 0.05 level (Table 4). However, KED predicts local extreme values of recharge that are mostly unrelated to nearby observed values. The extreme values are from regression model prediction or the drift component of KED, which is added to the kriged residual. In contrast, IK reduces the influence of observed, extreme values of recharge by transforming data to indicators (0, 1). The resulting IK map of probabilities represents a smoothly varying surface devoid of local extremes. Indicator kriging implicitly incorporates soils information because the recharge estimates are based on measured sediment properties. As a result, areas with low IK probabilities in Fig. 4 generally coincide with areas with high percentage clay in Fig. 2. For these reasons, IK is emphasized in this paper.

Multiple linear regression was used in a separate exercise to identify variables that significantly influence recharge in the study area. Surficial, percentage measured clay was highly significant (F-ratio p = 0.045) in the MLR model, and percentage well-drained soils was significant at the 0.10 level (p = 0.088 in Table 4). Coefficient signs indicate that recharge decreases with increasing clay content and increases with increasing percentages of well-drained soils. This is consistent with the spatial patterns of PTF-estimated recharge and clay content shown in Fig. 2, where fourth and fifth quintile values of recharge cluster in the east and southeast.

After building the MLR model, we entered the remaining variables one at a time and checked F-ratio p values for significance at the 0.05 level. Variables representing high and low infiltration soils yielded highly significant p values (Table 4), but the coefficient signs are difficult to interpret. For example, the sign of the high infiltration coefficient is negative, which suggests that recharge decreases with higher percentages of well-drained soils. We expected the reverse relation because these soils by definition have moderate to high rates of water transmission even when thoroughly wetted (USDA-NRCS, 1995). The high-infiltration areas might contain some fine-grained sediments at depth. Whereas SSURGO variables describe the upper 1.8 m of soil, recharge was estimated at depths closer to the water table (median depth of recharge calculation = 2.4 m). Alternatively, the high-infiltration areas might contain localized deposits of fine-grained sediments in closed depressions or as overbank flood deposits along streams that are too small to be mapped at the SSURGO scale. Such deposits would restrict recharge to the aquifer.

Topographic wetness index, land elevation, and depth to water are insignificant in regression; they have p values >0.4 when introduced into the MLR model (Table 4). The study area is in the Coastal Plain and is relatively flat. Topographic relief in the area likely is insufficient to significantly affect recharge at the regional scale. Lack of topographic relief is reflected in the low coefficients of variation (expressed as percent) of land elevation (21%) and topographic wetness index (25%). Sediment properties and recharge vary considerably more in the study area. The CV of measured moisture content is 56%. The CVs of SSURGO soil drainage variables are as high as 140%, and the CV of the recharge estimates is 220%. The wide ranges of these variables facilitate correlation in that high values of recharge generally coincide with high percentages of well-drained soils. The high CV of recharge is typical of soil hydraulic properties and reflects both the spatial variability of sediment properties and the uncertainty of PTF parameters. In prior research, CVs of hydraulic properties in field soils were as high as 124% for saturated hydraulic conductivity and 194% for ponded solute velocity (Corwin et al., 1997).

Relation of Recharge to Shallow Groundwater Quality
Both median NO3 concentration and atrazine percentage detection were lower for the high recharge category corresponding to developed lands, suggesting that higher recharge dilutes contaminant source concentrations for a given land use (Table 5). For example, median NO3 concentration is 6.6 mg L-1 for the low recharge category corresponding to old urban land and 2.7 mg L-1 for the high recharge category. Corresponding atrazine percentage detections are 75.0 and 37.5%, respectively. In undeveloped areas both high and low recharge categories have the same median NO3 concentration (0.07 mg L-1) and atrazine percentage detection (20%). A scatterplot of the middle 50% of recharge estimates indicates generally decreasing NO3 concentration with increasing recharge rate in developed areas (Fig. 6). The fitted lines in Fig. 6 are from linear regression.


View this table:
[in this window]
[in a new window]
  		Table 5. Median NO3 concentration and atrazine percentage detection for land uses and recharge categories in the Glassboro, NJ study area.

 


View larger version (33K):
[in this window]
[in a new window]
  		Fig. 6. Relation between NO3 concentration and estimated recharge in the Glassboro, NJ study area (middle 50% of recharge estimates is shown).

 
The dilution effect is consistent with results of prior researchers, who reported declines of about 30% or more in NO3, electrical conductivity, and silica concentrations in middle and deep (0.8 and 1.2 m) unsaturated zone lysimeters from flushing by low-residence time water in the Palouse River Basin of southeastern Washington State (Allen-King et al., 2001). The lysimeters were installed in a no-till agricultural field, and the Basin comprises loess deposits over basalt. Corroborating information includes a large precipitation event during above-freezing temperatures coincident with increases in water table height and soil temperature and a decrease in the {delta}18O content of deep lysimeter (1.2 m) pore water. These data indicate vertical migration to the deeper subsurface of warmer water coinciding with a winter thaw.

Comparatively low NO3 concentration in high recharge areas might also indicate denitrification caused by frequent and longer periods of saturation. Unsaturated zone NO3 concentrations, such as might be obtained with lysimeters, were unavailable in this study. However, we sampled NO3 and dissolved oxygen (DO) using shallow groundwater probes (0.6–0.9 m below water level) paired with 15 observation wells in the study area to assess denitrification potential in recently recharged groundwater. Compared with the observation wells, the shallow probes better characterize the quality of water entering the aquifer at the well location. Median NO3 concentration in groundwater samples from the probes is 3.3 mg L-1, and median DO is 6.7 mg L-1. The elevated NO3 concentration and well-oxygenated groundwater are inconsistent with denitrification, which is promoted by anaerobic conditions. Nitrate and DO concentrations are very similar for probes and paired observation wells, indicating a lack of vertical stratification of water quality in the upper part of the aquifer. Median concentrations of these constituents in samples from the 15 paired wells are 3.6 and 6.7 mg L-1, respectively. Compared with concentrations of NO3 and DO from the probes, these differences are statistically insignificant (p = 0.983 and 0.818, respectively).

Two-factor ANOVA on ranked NO3 and atrazine concentrations from observation wells was performed to test the effects of recharge category and type of developed land (agricultural, new urban, old urban) on water quality. Atrazine concentration (instead of percentage atrazine detection by group) was used in the ANOVA to facilitate computation of means of treatment groups. Considering only developed areas, differences in ranked NO3 concentration were significant for variables representing land use (p < 0.001) and recharge category (p = 0.024) (Fig. 7a and 7b), but differences in ranked atrazine concentration were not (p = 0.370 and 0.206, respectively) (Fig. 7c and 7d). Median atrazine concentration, however, was greater in the low recharge category (Fig. 7d), which is consistent with the dilution hypothesis.



View larger version (24K):
[in this window]
[in a new window]
  		Fig. 7. Distributions of (a, b) NO3 and (c, d) atrazine concentrations in relation to land-use and recharge categories (excludes undeveloped lands). ANOVA p values are based on ranked concentrations.

 
A GIS was used to compare agricultural, urban, and undeveloped lands with areas where recharge likely exceeds the median rate of 29.1 cm yr-1. Areas where the IK-predicted probability of exceeding median recharge is >0.5 were designated "high probability," and areas with IK-predicted probability <=0.5 were designated "low probability." High probability areas generally occur in the northwestern, northern, and southeastern portions of the study area (Fig. 8). Land use in the north and southeast is mainly urban and undeveloped. The low probability areas primarily coincide with agricultural and undeveloped lands in the southwest. The cross-hatched areas in Fig. 8 show where, within a given land use, contamination potential is lower because of possible dilution.



View larger version (52K):
[in this window]
[in a new window]
  		Fig. 8. Geographic information system overlay of land use and areas with high recharge probability (likely >29.1 cm yr-1) in the Glassboro, NJ study area.

 
Verification of Recharge Map
We attempted to validate the recharge map for a subset of 31 well locations not used in IK (those lacking measured sediment texture data and a corresponding recharge estimate) by comparing NO3 concentration in areas with low (<=0.5) and high (>0.5) predicted probabilities of exceeding median recharge. We refer to these 31 wells as validation sites because they are unsampled with respect to recharge. A total of 78 wells were sampled for NO3 and other compounds during previous water-quality surveys in the study area (Stackelberg et al., 1997, 2000). Sediment samples (used in recharge estimates) were collected at 47 of these wells and at one additional site not associated with an observation well. Twenty-four of the remaining 31 wells are in new and old urban areas. We performed a Wilcoxon Rank Sum test to compare NO3 data from areas with low and high recharge probabilities, using the 24 validation wells in urban areas. Lumping the new and old urban wells into a single category is reasonable because NO3 concentration is similar for the two groups (Fig. 7a). Stratification of the 31 validation wells by both land use and recharge category was impractical because of the small number of observations (typically 2 to 3) in each subcategory. Considering the 24 validation wells in urban areas, median NO3 concentration is 1.4 mg L-1 in the low probability category and 2.0 mg L-1 in the high probability category. These differences are statistically insignificant (p = 0.401), indicating that the strong relation between recharge estimates and NO3 concentration (ANOVA result in Fig. 7b) could not be validated in unsampled areas without accounting for differences in land use. Although we attempted to compensate for differences in loading by considering only urban wells, N load likely differs from town to town and even from house to house. Variation in N load likely obscures the effect of recharge on water quality at validation sites. The high p value might also be related to the small number of NO3 observations (11–13) in the recharge subcategories.

We repeated the Wilcoxon Rank Sum test using all 50 wells in urban areas (new and old combined) with water-quality data, including those with a recharge estimate. This attempted a general verification of probabilities predicted by IK, but is not validation per se. Median NO3 concentration is 3.5 mg L-1 in the low probability category and 2.5 mg L-1 in the high probability category in urban areas, suggesting dilution by recharge. These differences, however, are statistically insignificant (p = 0.915). Although median NO3 concentration is different for the two recharge probability categories, other percentiles of the distributions are similar.

In addition to variations in N load, IK prediction uncertainty might obscure differences between recharge probability categories. Equation [7] indicates that the greater the separation distance h between samples, the greater the variation in recharge indicators. Figure 9 shows the locations of urban wells with and without a recharge estimate and the standard deviations associated with the IK map of the probability of exceeding median recharge (Fig. 4). The IK standard deviations indicate that prediction uncertainty is greater at validation sites, where recharge data are sparse.



View larger version (91K):
[in this window]
[in a new window]
  		Fig. 9. Distribution of urban wells with and without (validation) recharge estimates in relation to indicator kriging standard deviations associated with kriged map of median recharge indicators.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
Estimates of groundwater recharge varied substantially (-18.5–1840 cm yr-1) over the 930-km2 study area near Glassboro, NJ. The wide range of values reflects both the spatial variability of sediment properties and the uncertainty of PTF estimates of recharge. However, the estimates are well within an order of magnitude of estimates based on the steady-state centrifuge method, and the median value compares favorably with estimates from a regional water budget. We mitigated the uncertainty of the PTF method by using categorical descriptions of recharge. Directional variogram analysis showed that variation in recharge indicators is anisotropic and coincides with the general northeast–southwest orientation of the Kirkwood-Cohansey aquifer. Variogram range in this direction (7600 m) apparently is related to cyclical patterns of sediment texture corresponding to the widths of low-lying areas. Indicator kriging reduced the influence of extreme values in the data set and yielded a map of the probability of exceeding median recharge (29.1 cm yr-1) in the study area. Low probabilities in the southwest coincide mainly with agricultural land over sediments with high clay content, and high probabilities in the east and southeast generally coincide with coarse grained deposits in low-lying areas. Percentages of measured clay and well-drained soils are significant in a multiple linear regression model; however, land elevation and topographic wetness index, which might have indicated topographic effects on recharge, are insignificant. Regression modeling results and the spatial patterns of recharge estimates, IK probabilities, and clay content indicate that sediment texture controls recharge in the study area. Median NO3 concentration and atrazine percentage detection in shallow groundwater are higher for the low recharge category (<=29.1 cm yr-1) than for the high recharge category (>29.1 cm yr-1) in agricultural, old urban, and new urban areas. In the case of NO3, these differences are significant at the 0.05 level. Results suggest that, for a given land use, NO3 and atrazine in shallow groundwaters of the study area are diluted by higher volumes of recharging water. Water quality–recharge relations, however, could not be validated using IK-predicted probabilities because of the small number of validation sites and potential variation in N load between sites.


   ACKNOWLEDGMENTS
 
We thank Kerie J. Hitt for preparing figures in geographic information systems. We also thank Timothy R. Green, Randall T. Hanson, Michael G. Rupert, Chester Zenone, and three anonymous reviewers for providing insightful reviews that substantially improved this paper.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow