Published online May 26, 2006

Modeling the Impact of Clustered Septic Tank Systems on Groundwater Quality
Liping Pang,* Chris Nokes, Jirka Šimůnek, Heather Kikkert, and Ross Hector
drinking water by the disposal system can be avoided.
Where a number of dwellings are clustered, however,some water supply wells in the community could be
affected by a neighboring disposal system, or they mightexperience the cumulative effect of several upgradientdisposal systems. For the purpose of making policy andoverseeing rural development, there has been a greater
requirement from government authorities to evaluate
the cumulative effects of clustered systems on groundwater quality.
Numerical models are useful tools for the quantitative
assessment of the impact of onsite systems on groundwater quality and the elucidation of the importance offactors that control contaminant concentrations in receiving waters. Studies on the development and applications of numerical models for evaluating the impact ofonsite wastewater systems on environmental quality aresparse (Shutter et al., 1994; MacQuarrie and Sudicky,
2001; Beach and McCray, 2003). Using a variably saturated flow and first-order transport model, Shutter et al.
(1994) simulated the movement of Na and a surfactant
from septic drain field. MacQuarrie and Sudicky (2001)
developed a multicomponent flow and reactive transport model, which couples the most relevant physical,
geochemical, and biochemical processes involved inwastewater plume evolution in sandy aquifers. They experimentally evaluated their model by simulating wastewater migration in a 1-m-long unsaturated column.
Beach and McCray (2003) applied HYDRUS-2D to simulate the influence of soil clogging of wastewater soil
absorption systems on flow regimes. McCray et al. (2005)
presented a critical review of model-input parameters tosimulate the transport of onsite wastewater pollutants.
Most studies reported in the literature, however, havefocused on evaluating the performance of onsite systems
or specific processes that are associated with onsite
systems (e.g., nitrification and denitrification), and little
has been reported on modeling the cumulative impact
of onsite systems on groundwater quality. The need toquantitatively predict potential cumulative effects ofonsite systems on groundwater quality has been indicatedin the review by McCray et al. (2005). They commented
that understanding the cumulative effects of onsite
systems is critical for determining the adverse effects on
nearby drinking water supply wells.
The impact of clustered septic tank systems on groundwater quality in New Zealand has previously been investigated in two field studies performed at a site near
Christchurch, in 1977 (Sinton, 1982) and in 1986 (Closeet al., 1989). Elevated NO32 levels and locally high fecal
coliform levels in groundwater were found in both surveys. The 1977 survey (Sinton, 1982) also found a clear
trend of increasing NO32 concentrations in the downgradient groundwater as the number of upgradient septictank systems increased, but this trend was not found for

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTContamination of groundwater from onsite disposal of septic tankeffluent is an increasing concern. In this study, HYDRUS-2D simulated the impact of clustered disposal systems on NO32 and fecal
coliforms in groundwater in a rural community near Christchurch,
New Zealand. The model included nine disposal boulder pits, embedded 4 m below the surface in alluvial gravel media, in a domain of
3.3 km by 30 m (including both unsaturated and saturated zones).
Water movement between the ground surface and the disposal pits was
simulated using HYDRUS-1D. The performance of the two-dimensional model was evaluated using monitoring data obtained from a
1977 study. Applying the daily climate data for 1974 to 1977, the simulated NO32 and bacteria concentrations in the groundwater weresimilar to those observed. Both observed and simulated results showedthat clustered disposal systems have a significant cumulative impact
on NO32 concentrations in groundwater, but the impact of fecal coliforms from individual systems is localized. This study further demonstrates that groundwater has a limited ability to dilute NO32,
requiring at least 2.9 km for NO32 to be reduced to near backgroundlevels. Therefore, disposal systems must treat effluent efficiently andwater wells downgradient of closely clustered disposal systems must be
deep enough to avoid the adverse health effects of NO32. Sensitivity
analysis suggests that the model results are most sensitive to changes
in hydraulic conductivity, effluent concentrations and discharge rate,
and the removal rate of bacteria in the unsaturated zone. Therefore,
the accurate estimation of these parameters is a fundamental requirement for the model to produce realistic results.

O

are used
worldwide in many rural and urban fringe areas.
These onsite systems are often developed as clusters,with a number of homes adjacent to one another, all on
separate onsite sewage disposal systems. This is particularly true in many parts of New Zealand, where there
is a trend toward subdividing rural land to form “lifestyle
blocks.” Dwellings within these blocks receive neither
reticulated water nor sewerage services, and often the
groundwater used for human consumption and hygieneis drawn from the same areas that are used for the disposal of sewage effluent. The potential threat to publichealth associated with this practice is an increasing
concern as wastewater can contain many potentiallyharmful contaminants, such as disease-causing microorganisms and NO32. Where a dwelling in a lifestyle
block is isolated, and its water supply well is located upgradient of the disposal system, contamination of the
NSITE WASTEWATER TREATMENT SYSTEMS

L. Pang, C. Nokes, H. Kikkert, and R. Hector, Inst. of Environmental
Science & Research Ltd., P.O. Box 29181, Christchurch, New Zealand;
J. Šimůnek, Dep. of Environmental Sciences, Univ. of California,
Riverside, CA 92521. Received 6 Sept. 2005. *Corresponding author(liping.pang@esr.cri.nz).
Published in Vadose Zone Journal 5:599–609 (2006).
Original Researchdoi:10.2136/vzj2005.0108ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

599

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

600

VADOSE ZONE J., VOL. 5, MAY 2006

fecal coliforms. The aim of our study was to establish acontaminant transport model that can predict the effects
of clustered disposal systems on groundwater quality,
using the same field site as Sinton (1982) as an example.
The capability of the contaminant transport model wasthen evaluated using historic monitoring data. The
transport of NO32 and fecal coliforms in the unsaturated
zone and groundwater was simulated, and simulated
concentrations were compared with those monitored ingroundwater. We also examined the relative importanceof the transport parameters to the model results. These
results will help with the design of future experiments sothat more accurate parameter values can be obtained forfurther model refinement.
MATERIALS AND METHODS
Study Site and Characteristics of Subsurface MediaThe site selected for modeling was a 3.3-km length of
Buchanans Road in the Yaldhurst area of Christchurch, New
Zealand. The uppermost layer (1.9 m) of soil is Waimakariri
deep silt loam on sand (interbeds of silt loam, sandy loam,
and sand), and the substrata below are alluvial gravels with asand matrix. No soil information was directly available fromthis site. Water retention data were therefore obtained froma nearby reference site with similar soil profiles. Table 1 lists
the hydraulic properties of the soil layers, determined fromwater retention data using the RETC program (van Genuchten et al., 1991).
In the 1970s and 1980s, this area of Christchurch was asemirural unsewered community. Each property was served byan individual septic tank system and a drinking water well.
Onsite disposal systems have been used in this area since the1950s, and the number of systems has increased over the years.
As a common practice in the past on the Canterbury Plains,
effluent was not disposed of through any soil layers but via aboulder pit (about 4 m deep and 1 m diameter), which was indirect contact with underlying alluvial gravels (Sinton, 1982).
The rate of effluent discharge into the each pit varied between
200 and 600 L d21.
Although direct discharge of septic tank effluent into gravel
was very effective for the rapid disposal of unwanted effluent,
unfortunately it was ineffective in the removal of contaminants. In a baseline groundwater quality survey in the Yaldhurst area performed in 1976, 33% of 120 household water
wells contained fecal coliforms or streptococci (Sinton, 1982).
Septic tank systems were considered to be the major source of
fecal contamination. Subsequently, 25 wells were selected in1977 for monitoring of fecal coliforms and NO32 (Sinton,1982). These 25 wells had indicator bacteria present andhigh levels of electrical conductivity in the 1976 survey. The
results of the 1977 survey showed that the drinking water wellsin this study site had elevated levels of NO32 (up to 25.8 mg
L21) and fecal coliforms [up to 74 cfu (100 mL)21]. The NO32
level in the background groundwater, immediately upgradient of the study area, was approximately 6 mg L21 and it
had zero fecal coliforms. This background NO32 concentration was estimated from the average of five samples taken in
1977 from two wells used by Environment Canterbury (aregional authority) to monitor groundwater quality. The depthof the groundwater table was about 12 to 16 m below groundlevel, with a hydraulic gradient of approximately 0.0012 to
0.0014 estimated from water levels. The groundwater flow
direction in this area was approximately parallel to the Buchanans Road so any cumulative effects resulting from onsite
effluent disposal systems should become apparent in thewells along the road. Figure 1 shows the locations and depthsof disposal pits and drinking water wells along the Buchanans
Road. Although there were many other septic tanks andwells in the same area at a distance from the Buchanans Road,
the cross-section of the two-dimensional model does not embody them.

Model Construction and Inputs
HYDRUS-2D was used to simulate contaminant transport
through the unsaturated zone and groundwater system.
HYDRUS-2D is a Microsoft Windows-based finite-element
model used to simulate water flow and contaminant transportin variably saturated porous media (Šimůnek et al., 1999). The
program solves the Richards equation for saturated–unsaturated water flow, uses a Fickian-based advection–dispersion
equation for contaminant transport, and includes provisions
for linear equilibrium adsorption, zero-order production, andfirst-order reduction. The governing equations are solved
using a Galerkin-type linear finite-element scheme.
To reduce the complexity of the problem and the simulation
time, a two-step approach was used in the modeling. The
model conceptualization is illustrated in Fig. 1. As the effluent
was generally discharged at a depth of 4 m (Sinton, 1982), a
one-dimensional model using HYDRUS-1D (Šimůnek et al.,1998) was first constructed for the unsaturated soil layers
above the 4-m depth. Nine soil layers were included in thismodel (Table 1). The top boundary of this one-dimensionalmodel received daily rainfall and potential evapotranspiration
using local climate data measured from 1974 to 1977. A freedrainage boundary was assigned at its bottom. Only flow
movement was simulated in the one-dimensional model.
Hydraulic conductivities of the soil layers considered in the
one-dimensional model were assumed not to be affected by

Table 1. Soil hydraulic properties used in the HYDRUS models.†
Model dimensions
1
1
1
1
1
1
1
1
1
2
2

Depth
m
0.27
0.38
0.54
1.13
1.24
1.60
1.80
1.90
4.00
12 to 1830.0

Lithologyrb
23

silt loamfine sandy loam
loamy fine sandfine sandcoarse sandsilt loamfine sandy loam
loamy sandgravels with sand matrix
unsaturated gravels with sand matrix
saturated gravels with sand matrix

g cm
1.38
1.38
1.38
1.38
1.38
1.38
1.38
1.38
1.38
2.12
2.12

Ksat
21

md
1.798
0.533
0.621
0.415
1.798
1.798
0.533
0.533
25.92
600.0
600.0

ur

us
3

23

m m
0.000
0.524
0.020
0.493
0.024
0.517
0.000
0.492
0.000
0.524
0.000
0.524
0.020
0.493
0.020
0.493
0.059
0.507
0.059
0.400
0.059
0.300

a

n

21

m
3.583
9.495
5.656
6.527
3.583
3.583
9.495
9.495
6.700
6.700
6.700

1.237
1.278
1.381
1.417
1.237
1.237
1.278
1.278
2.267
2.267
2.267

† rb, bulk density; Ksat, saturated hydraulic conductivity; ur , residual water content; us, saturated water content; a, coefficient in the soil water retentionfunction; n, exponent in the soil water retention function.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

601

Fig. 1. Schematic of the transport domain and conceptual model, including defined boundary conditions.

operation of effluent disposal (e.g., clogging) as disposal ofeffluent occurred below these layers.
The second model, using HYDRUS-2D, simulated bothflow movement and contaminant transport. The model
domain was 3.3 km wide and 30 m deep, and included an
8- to 12-m unsaturated zone and an 18- to 22-m saturated zone.
The grid consisted of a total of 6232 nodes and 12 062 finite
elements. Two material layers (the unsaturated zone and saturated zone) were distinguished (Table 1). Nine onsite sewage
disposal systems, located along the Buchanans Road downgradient from each other, were included in the model domain,with effluent discharge rates between 200 and 600 L d21. As
the site was not surveyed, locations of each disposal systemwere estimated from site visits. The grid was spaced verticallyat 0.1 to 0.5 m in the top 5 m (denser at the top) and 1 m
downward, and horizontally at 0.5 m for the effluent pit and its
adjacent nodes and then gradually became sparser by a factorof two for its up- and downgradient areas.
The top boundary of the two-dimensional model received
the daily soil water drainage from the bottom boundary of the
one-dimensional model for the areas without effluent disposal
pits, and constant effluent discharges for the areas with theeffluent pits. The left and right boundaries were assumed tohave constant pressure heads in the saturated zone, calculated
from the hydraulic gradient of the groundwater, and no flow for
the unsaturated zone. No flow was assumed at the bottom
boundary. Conceptualized boundary conditions are shown in
Fig. 1. No surface ponding of effluent and rainwater was allowed at the highly permeable gravel media. The permeabilityof coarse gravel media was much greater than the intensity ofrainfall and effluent discharge. First- and third-type solutetransport boundary conditions were applied for the effluent
pits and other boundaries, respectively, excluding the no-flow
boundary. Drinking water wells (Fig. 1) were treated as observation wells for groundwater quality and the effect ofpumping was not considered in the model.
The simulation period for both one- and two-dimensionalmodels started from winter 1974, about 3 yr before the 1977
monitoring data were obtained, to allow equilibration of thesystem. Winter was chosen as the starting point to allow a
relatively wet soil profile to be initially assigned to the model tofacilitate the stabilization of the simulated systems. To further
stabilize the system, the models were first run for a 20-d
simulation and its final pressure heads were imported as theinitial pressure heads of the models that were run for a 3-yr

simulation. Time discretizations were assigned as follows: initial time step 1 3 1023 d, minimum allowed time step 1 3 1024 d,
and maximum allowed time step 4 3 1023 d.
A longitudinal dispersivity value of 0.4 m was assigned forMaterial Layer 1 (depth 4 m) as dispersivity is about one orderof magnitude smaller than the transport distance (Gelhar
et al., 1992). The longitudinal dispersivity value for Material
Layer 2 was highly uncertain. In a review by Gelhar et al.
(1992), longitudinal dispersivities for a scale of 1 to 3.5 km
varied between 6 and 170 m. We were, however, careful about
using a very large dispersivity value, as the largest scale thathas highly reliable dispersivity data reported in the Gelhar
et al. (1992) survey was only about 250 m. Considering multiple inputs of contaminants, a longitudinal dispersivity valueof 10 m was assigned for Material Layer 2. The rationale forthe choice of this value will be further discussed below. Transverse dispersivity values were assumed to be one-tenth of theirrelevant longitudinal dispersivity values.
Choosing a suitable value of the hydraulic conductivity forheterogeneous coarse gravel media was also difficult. Hydraulic conductivities determined from pumping tests varied between 0.86 and 864 m d21 for clean sand and gravels on the
Canterbury Plains (North Canterbury Catchment Board andRegional Water Board, 1983), while those determined from a
tracer experiment performed at the nearby site at Burnham
varied from 5064 to 14 112 m d21 (Pang et al., 1998). Different
results obtained by the two methods were explained bydifferent segments of aquifers being measured and by thetracer study probably being affected by preferential flow.
While a pumping test deals with the entire pore network within
its zone of influence, a tracer experiment measures only thedomain between sampling wells and has the tendency to overestimate water fluxes since the tracer may move through preferential flow pathways. Considering that the model domain
is fairly large and probably involves both preferential andmatrix flow, we considered that the use of the K (hydraulicconductivity) value of 600 m d21 in our model was appropriate.
This value proved to be the most reasonable value in the sensitivity analysis.
We have assumed that the hydraulic conductivity of thesubsurface media did not change with time due to processes
such as clogging. No experimental data supporting the hypothesis of hydraulic conductivity changes due to operation of
effluent disposal were collected at the site, nor is relevant data
available in the literature for coarse gravel media. The impact

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

602

VADOSE ZONE J., VOL. 5, MAY 2006

of effluent disposal on hydraulic conductivity would be muchless significant in highly permeable coarse gravel media of thesite than in fine porous media (such as sands and silts). Alluvial
gravels on the Canterbury Plains are typically very permeable
and coarse, with an average of 10% finer than 0.90 (0.20–17)mm and 50% finer than 18.3 (12.9–34) mm, and with cobbles
up to 150 mm in size. Although the substrata also contain somesilt and clay, our field experiments suggested that solute transport in heterogeneous coarse gravels is dominantly in preferential flow paths (Pang et al., 1998).
In sewage effluent, N is almost entirely in the form of
dissolved NH41 and NO32 concentrations are insignificant(Wilhelm et al., 1994). Therefore, NH4+ and fecal coliforms
were introduced to the model domain through the flux andconcentration of the sewage discharge. As the actual concentrations of NH41 and fecal coliforms in sewage effluent were
not available for the study site, they were obtained from themean concentrations reported by Auckland Regional Council (2004), and values of 64 mg NH41 L21 and 106 fecal
coliforms per 100 mL were used in the model (Table 2). These
concentration values are typical for sewage effluent. In a review of model-input parameters for transport of onsite wastewater treatment by McCray et al. (2005), the median NH41
concentration was summarized as 75 mg L21 (n 5 37, Table 2).
In a summary of data from a number of studies, Pang et al.
(2003) reported that sewage effluent typically contains
106 fecal coliforms per 100 mL. Considering the relatively
dry climate in the Canterbury Plains, a mean effluent discharge
rate of 200 L person21 d21 was used. This assumes that allwater used by a household is discharged through the wastewater disposal system. This discharge rate is slightly lower than
the rate of 260 L person21 d21 summarized in a review by
McCray et al. (2005).
When dissolved NH41 leaves the anaerobic environment of
the septic tank, it is transformed to NO32 by autotrophic bacteria in the presence of O2 (MacQuarrie and Sudicky, 2001).
Ammonium is nitrified to NO32 via a three-species nitrification chain reaction (NH41 ! NO22 ! NO32) as follows:

NH4 1 1 1:5O2 ! NO2 2 1 H2 O 1 2H1[1]

NO2 2 1 0:5O2 ! NO3 2[2]

Under anoxic conditions and in the presence of an organic C
source, the reverse process occurs and NO32 is reduced to N2gas. The denitrification chain reaction isNO3 2 1 2CH3 OH ! 6NO2 2 1 2CO2 1 4H2 O[3]

6NO2 2 1 3CH3 OH ! 3N2 1 3CO2 1 3H2 O 1 6 OH2[4]These processes can be considered in the HYDRUS model.
Although an intermediate product, NO22, is involved in athree-species chain reaction, the conversion between NH41
and NO22 is so rapid (McCray et al., 2005) that the reaction canbe simplified to two single steps (NH41 ! NO32 or NO32 ! N2gas). The coupled first-order reaction rates for nitrification anddenitrification are

][NH1
4 ]
5 2l[NH1
4 ]
]t

[5]

][NO2
2
3 ]
5 l[NH1
4 ] 2 m[NO3 ]]t[6]where l and m are the nitrification and denitrification rate
coefficients [T21], respectively, and t is the time [T].
Denitrification was not considered in our study. Althoughgroundwater NH41 was not analyzed in the 1977 survey, it was
examined later in a survey in 1986 (Close et al., 1989). The
NH41 to NO32 ratio was estimated from the 1986 data of Closeet al. (1989) to be usually ,2% (i.e., NH41–N/NO32–N ratio#1%). This indicates that almost all inorganic N wasconverted to the oxidized form. Since groundwater was at
fairly oxidized conditions, denitrification was thus negligible
for the system investigated here. We believe that oxidation
also predominated over reduction in the overlying permeable
unsaturated gravels. Groundwater samples taken from asimilar field site 15 km nearby contained dissolved O2 at aconcentration of 11.6 mg L21 (Pang and Close, 1999), and anorganic C content of 0.04% (Pang et al., 2005). These findings
further support our inference of an unfavorable environment
for denitrification in coarse gravel aquifers.
Reference values for the first-order nitrification rate (l) andlinear adsorption coefficients of NH41 (Kd) are reported inTable 2. Parameters identified to be sensitive to the model results and their reasonable values.
Parameter

Reasonable range

Aquifer properties
Hydraulic conductivity for coarse gravels, m d21
Longitudinal dispersivity, m, for a scale of 1–3.5 km
Nitrogen1
21
NH4 concentration in septic tank effluent, mg L
Nitrification rate, d

1

21

NH4 adsorption coefficient, L kg

21

Fecal coliform
21
Removal rate in saturated gravels (d )21
Removal rate in unsaturated gravels, d )21Concentration in septic tank effluent, cfu (100 mL)

0.86–864

Reference

5,064–14,112 (n 5 15)
6–170

North Canterbury Catchment Board &Regional Water Board (1983)†
Pang et al. (1998)‡
Gelhar et al. (1992)

22–230 (n 5 37, mean 75)
51–77 (poor operation, mean 64)
26–39 (good operation, mean 33)
0.08–211 (n 5 19, mean 2.9)
0.02–0.5 (arable soils)
0.23–0.43 (soils)
0.15–0.25 (unsaturated soils),0.15 (coarse gravels)
3–4 (arable soils)
1.5 (soils)
3.5 (unsaturated soils),1.5 (coarse gravels)

McCray et al. (2005)
Auckland Regional Council (2004)
Auckland Regional Council (2004)
McCray et al. (2005)
Lotse et al. (1992)
Selim and Iskandar (1981)
Ling and Al-Kadi (1998)
Considered in this study
Lotse et al. (1992)
Selim and Iskandar (1981)
Ling and Al-Kadi (1998)
Considered in this study

1.14
.1.14
5
6
4 3 10 – 8.4 3 10

Sinton et al. (1997)§
Considered in this study
Pang et al. (2003)† Pumping test for coarse clean sand and gravel on Canterbury Plains.
‡ Tracer experiment at a nearby site 15 km away.
§ Determined from a tracer experiment in sewage-contaminated gravel aquifer at a nearby site.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

some studies (Table 2); however, most of these reported values
are derived from natural soils, not onsite wastewater systems
specifically (McCray et al., 2005). In addition, reported values
are predominately derived from cropped soils within the rootzones, where microbial activity, O2 levels, and the presence oforganic matter favor nitrification and adsorption of NH41 onto
soil media. In contrast, in our study septic tank effluent wasdischarged at a depth of 4 m, where O2 levels, microbial activity, and levels of organic matter in the coarse gravels are
relatively low. In addition, coarse gravel media are highly
permeable and heterogeneous. Taking all these into account,values of l 5 0.12 d21 and Kd 5 0.01 L kg21 were chosen in ourstudy. While these values are lower than most values reported
in the literature, we consider them to be reasonable for thesystem investigated. The impact of these parameter values onmodel results will be further discussed below. The same transformation rate was applied to both the liquid and solid phases.
Although an application of effluent would introduce microbes
and organic matter into the system, thus affecting nitrification rates, since no experiments and monitoring were performed for the disposal systems in the study area, we assumed
that these parameter values were constant.
The attenuation and transport of indicator bacteria in the
groundwater of alluvial gravels have been studied at two
nearby field sites, Templeton (Sinton et al., 1997) and Burnham (Pang et al., 1998; Sinton et al., 2000; Pang et al., 2005).
These studies showed that bacteria are not retarded in the
groundwater of coarse gravel aquifers. Therefore, the adsorption coefficient for the bacteria was set to zero in the model. A
removal rate of 1.14 d21 for fecal coliforms was adopted from
Sinton et al. (1997) for the sewage-contaminated aquifer. Thisrate encompasses the effects of all irreversible processes (e.g.,
die-off, filtration, and irreversible sorption). Pang et al. (2005)
demonstrated that the removal rate of microbes is generally
one order of magnitude lower in sewage-contaminated aquifers than it is in clean aquifers. The value of the removal
rate for the unsaturated zone was the most uncertain, compared with all the other model parameters, and it was expected
to be greater than that for the saturated zone. Since there
were no reference values available in the literature for bacterial removal rates in unsaturated coarse gravel media, theremoval rate was manually calibrated by trial and error to
obtain a reasonable result from the simulations. A removal
rate of 3 d21 was found to be the optimal value for fecal
coliforms, and this was demonstrated in the sensitivity analysis. An input concentration of 106 fecal coliforms per 100 mL
was used.

Model Sensitivity Analysis
As the large two-dimensional model required a long CPU(central processing unit) time (.2 d on a 3-GHz computer)and thus inverse modeling was practically impossible; only
direct modeling could be completed. A sensitivity analysis was
performed since the values for a number of input parameters
were highly uncertain due to the wide range of possible values(Table 2). The following parameters with the most uncertain
values in Table 2 were examined in the sensitivity analysis:
hydraulic conductivity, the adsorption coefficient for NH41, therate coefficient for the transformation of NH41 to NO32, longitudinal dispersivity, and NH41 input concentration. In addition, the discharge rate of the effluent was also included in thesensitivity analysis. Model simulations were performed within
the possible ranges of parameter values (Table 2) to examine
their influence on the model’s output. One parameter wastested at a time while all other parameters were fixed. The best
manually calibrated values were used as the baseline for com-603

parison of the model results, together with the change in maximum NO32 concentration.
The sensitivity analysis not only assists in understanding
where weaknesses in the modeling results may lie, but it also
assists with data interpretations, and provides guidance in
identifying which parameter values must be refined by future
work, using field measurements or laboratory experiments.
Model runs were also performed to demonstrate the effect ofthe density of disposal systems on groundwater quality and to
determine the separation distance between two disposal systems in which NO32 concentrations in downgradient groundwater can be reduced to a near natural level.
Since neither the observed nor model-simulated data showeda significant cumulative impact of fecal coliform contaminationin groundwater under the influence of clustered disposal systems (see below), the focus of the sensitivity analysis was mainly
on simulations of NO32. Sensitivity analysis for the bacteria
was only performed for the removal rate in the unsaturated
gravels, as that value was the most uncertain compared with allother model parameters, as mentioned above.

RESULTS AND DISCUSSIONS
Impact of Clustered Septic Tank Systems on
Groundwater Quality
NitrateFigure 2 shows simulated NO32 plumes in the verticalplane at 1000 d, which suggests a clear cumulative impact of clustered septic tank systems on the NO32 levelin groundwater. This cumulative effect can also be seenclearly from the NO32 concentrations at selected depths
across the model domain (Fig. 3). Figure 3 shows that, as
the number of disposal systems increases, the NO32 concentrations in downgradient groundwater increases
cumulatively. At 1 m below the water table, NO32 concentrations rise sharply in association with each disposal
system, followed by a sharp decrease with distance due tothe dilution effect, before the next discharge point. As
the depth increases, however, this sharp change in NO32
concentrations becomes much less apparent. At 10 m
below the water table, NO32 concentrations show a steep
increase in association with a disposal system, but remain
relatively constant between disposal systems.
The density impact of disposal systems on NO32
accumulation in groundwater is further demonstrated in
two hypothetical cases given in Fig. 4. Figure 4 showsthat when there are three upgradient disposal systems,
NO32 concentrations in downgradient groundwater at
the same distance (e.g., 3.3 km) are higher than when
there are only two upgradient disposal systems. Therefore, we can conclude that, with an increase in the densityof clustered onsite systems, a greater drinking water welldepth is required to obtain less contaminated groundwater. We need to note that this conclusion assumes thesame discharge rate for all systems. With fewer systems
in place, a greater well depth is also required when loading of the system is increased. We should also point out
that obtaining less contaminated groundwater by increasing the well depth is a passive approach. The bestsolution is to increase the efficiency of effluent treatment
within septic tanks before the effluent is discharged into
disposal pits.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

604

VADOSE ZONE J., VOL. 5, MAY 2006

Fig. 2. Hydrus-2D simulated NO32 plumes developed under the impact of clustered septic tank systems (results at 1000 d).

An explanation for NO32 accumulation in groundwater is that NO32 is a conservative solute and its
concentrations can only be reduced by dilution when
denitrification is negligible in the investigated system;
however, the capacity of groundwater to dilute NO32
concentrations is very limited, and a very large distanceis needed for NO32 levels to reduce to the backgroundlevel. This is illustrated in Fig. 5, which shows that only
when the distance is greater than 2.9 km can NO32 indowngradient groundwater be reduced to near background levels. Such a large separation distance is not
feasible, and thus it is critical for the disposal systems totreat effluent efficiently. We note here that since we used
a two-dimensional model that cannot account for
dilution in the transverse direction (the direction thatis perpendicular to the transport domain), our results
are more conservative than they would be in reality. A
review by McCray et al. (2005) also discussed the limited
dilution of NO32 in groundwater for onsite wastewater
systems. They indicated that although dilution may playa part in reducing N concentrations in groundwaterin the short term, it is not a practical long-term solution because as the area within which unsewered development increases along with the density of onsite
disposal systems, so the dilution capacity of the aquifer
media is diminished.
Model-simulated NO32 concentrations in groundwater, using the best manually calibrated parameter values,

Fig. 3. Simulated NO32 concentrations in groundwater 1, 4, and 10 m
below the water table across the model domain (results at 1000 d).

compared reasonably well with those observed in the
1977 study (Sinton, 1982), as shown in Fig. 6. Both observed and simulated datasets show a trend of increasing
NO32 concentrations with distance. It should be notedthat to make an unbiased comparison, the data presented in Fig. 6 were only for the nine wells that had thesame sampling dates (1 July and 1 Aug. 1977), and thesame depths were matched when comparing observedwith predicted data. Although the model was only manually calibrated, the good agreement between predicted
and observed NO32 concentrations suggests that themodel we have constructed can simulate the problem
under investigation and that the input parameter values
are reasonable. Figure 6 shows that simulated results
using K 5 600 and 1000 m d21 match equally well withobserved results (with the same r 2 5 0.72 for the 1 July
1977 data and r2 = 0.57 for the 1 Aug. 1977 data for linear
regression relationships between simulated and predicted values). It should be noted that, as mentioned
above, our two-dimensional model cannot account for
dilution (although relatively small) in the direction perpendicular to the transport domain. Therefore, the good
agreement between concentrations predicted with thetwo-dimensional model and those observed in the threedimensional field means that the hydraulic conductivity chosen for the two-dimensional model is probably
higher than that in the field.

Fig. 4. Comparison of simulated NO32 concentrations in groundwater
4 m below the water table in the presence of two and three upgradient disposal systems (results at 1000 d).

605

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 5. The effect of separation distance between two adjacent disposal systems on NO32 accumulation in groundwater (results at 1000 d,
4 m below the water table).

Although the cumulative impact of clustered disposal
systems is shown by both observed and model-simulateddata, the NO32 concentrations in groundwater are still
below the World Health Organization drinking water
guideline of 50 mg L21 (WHO, 2004), even where a highinput effluent NH41 concentration was used. The WHO
guideline of 50 mg L21 for NO32 in drinking water is also
the maximum acceptable value allowed in New Zealand(Ministry of Health, 2000). When drinking wells are 10 m
below the water table, the NO32 concentrations in groundwater can be reduced to a level of ,23 mg L21 (Fig. 3), andtherefore the adverse health effect of NO32 in drinking
water can be significantly reduced.
Fecal Coliforms
In contrast to the results of NO32, model predictions
do not show cumulative effects in relation to the concentrations of fecal coliforms in groundwater as a result ofclustered disposal systems—there is a localized impact
on bacterial concentrations (Fig. 7). This is consistent
with Sinton’s (1982) finding that there was no correlation
between microbial contamination and distance in the
direction of groundwater flow. This can be explained by
the efficiency of bacterial removal processes (e.g., straining, attachment, and die-off) in the subsurface media,

Fig. 6. Comparison of NO32 concentrations in groundwater simulated
in this study (at hydraulic conductivity K values of 600 and 1000 md21) and observed (Obs) in 1977 by Sinton (1982) for the nine wells
that were sampled on the same dates (1 July and 1 Aug. 1977).

Fig. 7. Simulated concentrations of fecal coliforms in groundwater 1 m
below the water table (results at 1000 d).

preventing the bacteria from being carried far downstream and reducing the likelihood of cumulative effects
becoming apparent.
Simulated bacteria concentrations decrease withdepth; however, to avoid the diagram being cluttered,
only results at 1 m below the water table are shown in
Fig. 7. Simulated fecal coliform concentrations in groundwater are mostly ,16 cfu (100 mL)21. These results agreewith those observed in 1977 [mean 2.9 cfu (100mL)21,
ranging 0–74 cfu (100mL)21] by Sinton (1982). Unfortunately, fecal coliforms at individual sampling wells measured in the 1977 survey were not given in Sinton (1982)
and, thus, we were unable to make any further comparison between predictions and observations. Nevertheless, the good match in the magnitude of observed
and predicted bacteria concentrations suggests that themodel we have constructed and its inputs describe reasonably well the transport of fecal coliforms in the gravel
media investigated.

Sensitivity Analysis
Considering NO32 concentrations at a specific location (horizontal distance 1907 m, 4 m below water table)
and time (1000 d), the sensitivity of the model results tothe input parameters within the range of their reasonable values is shown in Table 3 and Fig. 8. It is very clearthat K, the NH41 input concentration (Co), and thedischarge rate (Q) have the most significant influence on
the model results. In contrast, the nitrification rate coefficient (l), longitudinal dispersivity (ax), and adsorption coefficient for NH41 (Kd) have much less influenceon the model results. Figure 8 illustrates that the simulated NO32 concentrations are best described with apower function for the impact of K, Kd, ax, and l (witha relative importance of K . ax . Kd . l and with alinear function for the impact of Co (and probably alsoQ). The simulated NO32 concentrations are positively
related to Co, Q, and l but inversely related to K, Kd,
and ax. The relative importance of parameters Co . Q .
K . ax . Kd . l is also shown in Table 3, and is consistent with the results shown in Fig. 8.
The impact of input parameters on model results isfurther demonstrated in Fig. 9, which compares simulated NO32 concentrations 4 m below the water table
across the model domain. Again, K, Co, and Q are
shown to be the most sensitive parameters, and l, ax,

606

VADOSE ZONE J., VOL. 5, MAY 2006

centrations. Therefore our choice of K 5 600 m d21 as
the most likely K input value is reasonable. Similar results were obtained with K 5 5000 and 10 000 m d21, withpredicted NO32 concentrations close to the backgroundlevel. As simulated NO32 concentrations obtained with
K 5 100 m d21 were unrealistically high (up to 70 mg
L21) and their inclusion would have cluttered the plots,
we have excluded them from Fig. 9a. The impact of theK value on model results is reflected in both the magnitude and shape of the NO32 concentration profile as
K is directly proportional to the groundwater flux. As
K (and thus flux) increases, NO32 concentrations ingroundwater decrease, and the peak concentrations associated with the disposal systems become less pronounced. In contrast, at low K (low flux), NO32 from theleached effluent is not rapidly flushed away, which results in a substantial increase in the NO32 concentration.
The Kd coefficient for NH41 does not affect the model
predictions no matter which value is selected from a
relatively large interval of tested values. Therefore ourchoice of Kd 5 0.01 L kg21 as the best input Kd value
does not significantly affect model results, or theconclusions based on them. The insensitivity of themodel results to the Kd value is a consequence of a rapidtransformation of NH41 to NO32. The transformation of
NH41 to NO32 is a relatively fast process compared withother processes occurring in the system as all tested
l values (0.06–211 d21) correspond with a half-life of,12 d. It can thus be expected that the transformation
is mostly completed within the unsaturated zone. For
the same reason, the model results are also insensitive
to l, even though the range of tested l values is reasonably wide. When l , 2.9 d21, the model results are very
similar no matter which value is selected. Therefore ourchoice of l 5 0.12 d21 for the model input would also
have no major effect on the model results.

Table 3. Parameter values used in the sensitivity analysis and achange in the maximum concentration (Cmax).
FactorValue changed

Parameter

Cmax

Change in Cmax

21

1

NH4

21

21

concentration, mg L

321Discharge rate Q, m dNitrification rate, d

21

1

NH4 adsorption coefficient,
L kg21Longitudinal dispersivity, m

100
0.17
300
0.50
600
1.00
1 000
1.67
5 000
8.33
10 000
16.67
22
0.34
32‡
0.50
64§
1.00
75
1.17
230
3.59
1/2 Q
0.50
Q¶
1.00
0.06
0.50
0.08
0.67
0.12
1.00
2.9
24.2
211
1758
0.005
0.50
0.01
1.00
0.1
10.0
1.5
150
4
400
5
0.50
10
1.00
50
5.00
100
10.0

mg L
56.5
45.1
33.6
24.0
1.89
1.89
15.2
19.5
33.6
37.6
97.9
20.7
33.6
32.8
33.1
33.6
33.6
34.2
33.6
33.6
33.0
27.5
27.5
33.0
33.6
22.6
20.4

%
68.2
34.2
0.00
228.6
294.4
294.4
254.6
242.0
0.00
11.90
191
238.4
0.00
22.38
21.5
0.00
0.00
1.79
0.00
0.00
21.79
218.2
218.2
21.79
0.00
232.7
239.3

† Italic numbers are the best values used in the model.
1‡ Good operation (Auckland Regional Council, 2004): NH4 5 25.72–
21
38.57 mg L .
1§ Poor operation (Auckland Regional Council, 2004): NH4 5 51.43–
21
77.14 mg L .
21 21¶ Q 5 no. of residents 3 200 L person d .

and Kd have little impact on the model results. Of thethree most sensitive parameters, possible K values forheterogeneous gravel media have the widest range.
Figure 9a shows that K values between 600 and 1000 m
d21 give the most reasonable predication of NO32 con120

120
hydraulic conductivity

90

40

+

input NH4 concentration-0.87

y = 0.40x +7.33
2
R = 1.0

30

0

10

0
0

4000

8000

0
0

12000

y = 25.8x + 7.8
2
R = 1.0

20

60

30

effluent discharge rate30

90

y = 5432.5x
R2 = 0.90

60

NO3- (mg/L)

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Hydraulic conductivity, m d

100

K (m/day)

200

300

0

Co (mg/L)

40

40

30

30

0.5

1

1.5

Factor change [-]40

NH4+ adsorption coefficientlongitudinal dispersivity
30

nitrification rate
20

0.004

y = 33.485x2
R = 0.77

10

80

160

λ (day-1)

24020

y = 28.79x-0.03
2
R =0.89

10

0
0

20

y = 46.75x-0.18
2
R = 0.94

10
0

0
0

1

2

3

Kd (ml/g)

4

5

0

40

α x (m)80

120

Fig. 8. Sensitivity of NO32 concentrations simulated using the range of parameter values for a specific location (distance 1907 m, 4 m below watertable) and time (1000 d).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

607

Fig. 9. Sensitivity of NO32 concentrations simulated using a range of parameter values for the cross-section 4 m below the water table at 1000 d.

Similarly, model results change only slightly for
different values of the longitudinal dispersivity (ax 5
5–100 m were considered). Our choice of ax 5 10 m
would thus have not impacted the model results in any
significant way. We believe that the insensitivity of themodel results to the ax parameter arises because thephysical dispersion is overshadowed by the numericaldispersion involved in large-scale models where manyhorizontal grid spaces are much greater than the dispersivity values used in the sensitivity analysis. Therefore,
the difference in the simulated NO32 concentrations
when using different dispersivity values might well bewithin the error caused by numerical dispersion. Thiserror caused by numerical dispersion was inevitable forthe large model domain that we constructed. If thehorizontal discretization was finer than the dispersivity
value of 10 m, the model simulation time would havedramatically increased.
As mentioned above, sensitivity analysis for the bacteria predictions was only performed for the removal

rate in the unsaturated zone. As expected, the level of
fecal coliforms in groundwater is very sensitive to theremoval rate in the unsaturated zone (Fig. 10). Thisemphasizes the important role that the unsaturated zone
plays in reducing bacterial contamination in groundwater. Of the removal rate coefficients examined(1.5–11.4 d21), a removal rate of 3 d21 for the unsaturated zone seemed to yield the most satisfactory resultsin comparison with observed results.

CONCLUSIONS AND IMPLICATIONS
The large-scale transport model that was constructedbased on HYDRUS-2D was satisfactorily evaluatedusing the historic groundwater monitoring data. Bothmodel simulations and field observations have demonstrated that, while the clustered septic tank systems havea cumulative impact on NO32 concentrations in groundwater, they have a localized impact on fecal coliform
concentrations in groundwater. These contrasting re-608

VADOSE ZONE J., VOL. 5, MAY 2006

sitive to the NH41 adsorption coefficient, nitrification
rate, and longitudinal dispersivity. The relative importance of model parameters to model outputs provides
useful information for the design of our future experiments and subsequent model refinement.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ACKNOWLEDGMENTS
We wish to thank Environmental Canterbury and thehouseholders within the study area for their help in providing
information for this project. This study was funded by the New
Zealand Foundation for Research, Science, and Technology(Contract CO3X0303).

Fig. 10. Simulated concentrations of fecal coliforms (on a logarithmic
scale) at 1, 2, 3, and 4 m below the water table (mb.w.1) at a distanceof 2005 m. The sensitivity analysis was performed with different
removal rates in the unsaturated zone (results at 1000 d).

sults for NO32 and fecal coliforms reflect the fact that
NO32 is a conservative solute that can only be attenuated by dilution (as denitrification is negligible), whilebacteria are effectively removed in subsurface media(especially through the unsaturated zone) by processes
such as filtration and die-off. The predicted levels of
NO32 and fecal coliforms in groundwater also agreed
well with those observed, suggesting that the constructed
model and its inputs are reasonable for describing nitrification and bacteria transport in alluvial gravel media
under the influence of clustered septic tank systems.
The results of the model simulations suggest that
NO32 accumulates in groundwater as the density of
upgradient clustered onsite systems increases, and that
reduction of NO32 in groundwater by dilution is limited.
For NO32 levels to reduce to a near background level,
the separation distance required between two adjacentdisposal systems would have to be at least 2.9 km in thecoarse gravel aquifers investigated. It is impractical touse such a large distance for regulatory purposes. Therefore, to minimize contamination of groundwater, it isessential for the waste disposal systems to treat effluent efficiently. In addition, drinking water wells downgradient of the clustered waste disposal systems should
be at a sufficient depth, ideally .10 m below the watertable, so that NO32 concentrations can be reduced to thelevel that has a minimal health impact. In spite of theelevated NO32 levels in groundwater within the study
area, the NO32 concentrations observed and simulated
were still below the WHO drinking water guidelines of
50 mg L21, even when the NH4+ concentration in theeffluent was high.
Simulated NO32 results were most sensitive to the
input values for the hydraulic conductivity, effluent concentrations, and the discharge rate. This indicates thatboth the aquifer properties and the operation of thewaste treatment and disposal systems are critical to thepotential impact of disposal systems on groundwater
quality. Therefore, an accurate estimation of these parameters is the fundamental requirement to enable themodel to make reliable predictions. For the ranges of
selected parameter values, the model results are insen-

REFERENCES
Auckland Regional Council. 2004. Onsite wastewater systems: Design
and management manual. 3rd ed. ARC Tech. Publ. no. 58. ARC,
Auckland, New Zealand.
Beach, D.N.H., and J.E. McCray. 2003. Numerical modeling of unsaturated flow in wastewater soil absorption systems. Ground
Water Monit. Rem. 23(2):64–72.
Close, M.E., L.R. Hodgson, and G. Tod. 1989. Field evaluation offluorescent whitening agents and sodium tripolyphosphate as indicators of septic tank contamination in domestic wells. N.Z. J. Mar.
Freshwater Res. 23:563–568.
Gelhar, L.W., C. Welty, and K.R. Rehfeldt. 1992. A critical review ofdata on field-scale dispersion in aquifers. Water Resour. Res. 28:
1955–1974.
Ling, G., and A.I. Al-Kadi. 1998. A lumped parameter model for
nitrogen transformation in the unsaturated zone. Water Resour.
Res. 34:203–212.
Lotse, E.G., J.D. Jabro, K.E. Simmons, and D.E. Baker. 1992. Simulation of nitrogen dynamics and leaching from arable soils. J. Contam. Hydrol. 10:183–196.
MacQuarrie, K.T.B., and E.A. Sudicky. 2001. Multicomponent simulation of wastewater-derived nitrogen and carbon in shallow unconfined aquifers I. Model formulation and performance. J.
Contam. Hydrol. 47:53–84.
McCray, J.E., S.L. Kirkland, R.L. Siegrist, and G.D. Thyne. 2005.
Model parameters for simulating fate and transport of on-sitewastewater nutrients. Ground Water 43:628–639.
Ministry of Health. 2000. Drinking-water standards for New Zealand
2000. Available online at http://www.moh.govt.nz/moh.nsf/by1unid/
70727DB605B9F56A4C25696400802887?Open (verified 12 Apr.
2006). Ministry of Health, Wellington, New Zealand.
North Canterbury Catchment Board and Regional Water Board. 1983.
Interim report on the groundwater resource of the central plains.
North Canterbury Catchment Board and Regional Water Board,
Christchurch, New Zealand.
Pang, L., M.E. Close, and M. Noonan. 1998. Rhodamine WT and
Bacillus subtilis transport through an alluvial gravel aquifer. Ground
Water 36:112–122.
Pang, L., and M.E. Close. 1999. Attenuation and transport of atrazine
and picloram in an alluvial gravel aquifer. N.Z. J. Mar. Freshwater
Res. 33(2):279–291.
Pang, L., M. Close, M. Goltz, L. Sinton, H. Davies, C. Hall, and G.
Stanton. 2003. Estimation of septic tank setback distances based ontransport of E. coli and F-RNA phages. Environ. Int. 29:907–921.
Pang, L., M. Close, M. Goltz, M. Noonan, and L. Sinton. 2005. Filtration and transport of Bacillus subtilis spores and the F-RNA
phage MS2 in a coarse alluvial gravel aquifer: Implications in theestimation of setback distances. J. Contam. Hydrol. 77:165–194.
Selim, H.M., and I.K. Iskandar. 1981. Modeling nitrogen transport andtransformations in soils: 2. Validation. Soil Sci. 131:303–312.
Shutter, S.B., E.A. Sudicky, and W.D. Robertson. 1994. Chemical fate
and transport in a domestic septic system: Application of a variably
saturated model for chemical movement. Environ. Toxicol. Chem.
13(2):223–231.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1998. The
HYDRUS-1D software package for simulating the one-dimensional
movement of water, heat, and multiple solutes in variably-saturated

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

media. Version 2.0. IGWMC TPS-70. Int. Ground Water Modeling
Ctr., Colorado School of Mines, Golden.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiple solutes in variably saturated
media. Version 2.0. IGWMC TPS-53. Int. Ground Water Modeling
Ctr., Colorado School of Mines, Golden.
Sinton, L. 1982. A groundwater quality survey of an unsewered, semirural area. N.Z. J. Mar. Freshwater Res. 16:317–326.
Sinton, L.W., R.K. Finlay, L. Pang, and D.M. Scott. 1997. Transport ofbacteria and bacteriophages in irrigated effluent into and throughan alluvial gravel aquifer. Water Air Soil Pollut. 98(1–2):17–48.

609

Sinton, L.W., M.J. Noonan, R.K. Finlay, L. Pang, and M.E. Close. 2000.
Transport and attenuation of bacteria and bacteriophage in an alluvial gravel aquifer. N.Z. J. Mar. Freshwater Res. 34:175–186.
van Genuchten, M.Th., F.L. Leij, and S.R. Yates. 1991. The RETC
Code for quantifying the hydraulic functions of unsaturated soils.
U.S. Environ. Protection Agency, R.S. Kerr Environ. Res. Lab.,
Ada, OK.
World Health Organization. 2004. Guidelines for drinking-water quality: Vol. 1. Recommendations. 3rd ed. WHO, Geneva, Switzerland.
Wilhelm, S.R., S.L. Schiff, and W.D. Robertson. 1994. Chemical fate
and transport in a domestic septic system: Unsaturated and saturated zone geochemistry. Environ. Toxicol. Chem. 13(2):193–203.