Published online 17 May 2007
Published in Vadose Zone J 6:387-396 (2007)
DOI: 10.2136/vzj2005.0113
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: SAVANNAH RIVER SITE
Simulating Tritium Fluxes in the Vadose Zone under Transient Saturated Conditions
Karin T. Rebela,d,*,
Susan J. Rihaa,
Derek Karssenbergc and
Jery R. Stedingerb
a Dep. of Earth and Atmospheric Sciences, Cornell Univ., Ithaca, NY 14853
b School of Civil and Environmental Engineering, Cornell Univ., Ithaca, NY 14853
c Dep. of Physical Geography, Utrecht Univ., Utrecht, the Netherlands
d current address: Dep. of Hydrology and Geo-Environmental Sciences, Vrije Universiteit, Amsterdam, the Netherlands
* Corresponding author (karin.rebel{at}falw.vu.nl).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 15 September 2005.
 |
ABSTRACT
|
|---|
A phytoremediation project was established at the U.S. Department of Energy's Savannah River Site (South Carolina) to reduce fluxes of tritium-contaminated groundwater to surface waters. Contaminated groundwater was collected in a pond and applied by spray irrigation to a catchment of mixed forest on Atlantic Coastal Plain soils. The objectives of this research project were to simulate tritium uptake by the vegetation and to determine if subsurface lateral flow at the sand-clay interface impacts tritium uptake by forest vegetation. To simulate water and tritium fluxes within the catchment, we developed a spatially explicit water and solute transport model. Vertical water flow was simulated within a grid using a layered capacitance model with lateral flow between cells dependent on head development and the local slope of the impeding clay layer. Tritium movement was simulated on a daily basis assuming complete mixing within a cell. The model was evaluated by comparing biweekly measurements of soil tritium activity and soil water content in 18 measurement clusters distributed across the irrigated portion of the catchment. Although lateral flow was predicted to occur locally, after 3 yr of tritium irrigation, the model predicted that lateral transfer of tritium was <70 m downslope of the irrigation area. Transient saturated conditions were observed and simulated to occur within parts of the catchment. Including the simulation of subsurface lateral flow of water and tritium under transient saturated conditions, however, did not significantly alter average tritium uptake predictions for the catchment.
Abbreviations: ldd, local drainage direction map • SRS, Savannah River Site • WTB, water and tracer balance.
|
INTRODUCTION
|
|---|
Simulation of solute fluxes in complex terrain typically fails to account for subsurface lateral movement of solutes under transient saturated conditions. This type of solute flux could be important in phytoremediation (the use of plants to remove or neutralize contaminants in polluted soil or water), and more generally in evaluating the fate of solutes applied to soil. Solutes could move laterally and remain within the root zone, and thus could be taken up by the vegetation rather than leaching or moving offsite.
A number of catchment-scale models have been developed that simulate transient saturated conditions in the vadose zone. The Variable Source Area (VSA) approach recognized that some areas within a catchment contribute to surface and subsurface runoff more than others, and that these areas change with time, expanding and contracting as a function of rain events and antecedent moisture (Hewlett and Hibbert, 1967). Simple assumptions have been made to represent VSA processes in models: declining saturated hydraulic conductivity with depth, steady-state catchment water table response, topographically defined water flow paths, and linear wetting and drying from the valley bottom upward to the ridge; examples of such models are the Soil Moisture Routing model (SMR) (Frankenberger et al., 1999) and TOPMODEL (Beven and Kirkby, 1979). The SMR represents the vadose zone using a water balance approach, where subsurface lateral flow is simulated using a simplified Richards equation. TOPMODEL uses geographic similarity (i.e., contributing area and slope of soil surface) as a basis for modeling saturated zone catchment response; Shaman et al. (2002) added simulation of subsurface lateral flow using the same TOPMODEL assumptions. Band et al. (1993) coupled TOPMODEL and the ecosystem model BIOME-BGC to account for the impact of lateral fluxes on ecological processes. In this coupled model, simulated water and solute fluxes occur in the saturated zone rather than in the vadose zone. Wigmosta et al. (1994) introduced DHSVM, a semi-three-dimensional hydrology vegetation model that includes canopy interception, evaporation, transpiration, and snow accumulation. Whereas TOPMODEL was not spatially explicit, DHSVM was one of the first spatially explicit models that simulated transient saturated conditions in the vadose zone.
We developed a semi-three-dimensional hydrology–vegetation model, similar to DHSVM (Wigmosta et al., 1994) but incorporating multiple soil layers, solute movement, and a more detailed water uptake routine. To calculate phytoremediation efficiency (contaminant uptake/contaminant applied), it is important to calculate contaminant activity with depth across the catchment, so activity is known in a spatially explicit fashion within the model domain.
Flow routing routines have been a topic of discussion during the last decade. Wigmosta and Lettenmaier (1999) concluded that explicit cell-to-cell routing of subsurface flow performed better than the statistical dynamic method used in TOPMODEL. In their comparison, the direction of the subsurface flow was simulated (as is frequently done) using the hydrologic gradients derived from the surface topography (Beven et al., 1995; Frankenberger et al., 1999). Recent studies have found that the use of the surface topography, estimated using digital elevation maps, may not be a good indicator of the probable occurrence of lateral subsurface flow pathways (Burt and Butcher, 1985; Buttle and McDonald, 2002; Moore and Thompson, 1996). Freer et al. (1997) suggested that an improvement would be to use subsurface topography of the soil–bedrock interface or other hydrologically impeding layers. Freer et al. (1997) studied two hillslopes in different hydrogeological–climatic settings and related observed patterns of subsurface storm flow to model predictions based on either surface or subsurface topographic indices, i.e., the index of Kirkby, which divides the upslope contributing area by the tangent of the local slope angle (Beven and Kirkby, 1979). They found that, where bedrock topography is distinctly different from surficial topography, the bedrock surface has a considerable influence on local hydrological gradients and therefore the dominant flow path directions.
In this study, we explored alternative lateral flow routing routines and their impact on flow patterns. We propose a new routing routine, where subsurface lateral flow is routed over the impeding layer, but the direction of the flow is adjusted for water buildup above the impeding layer, i.e., head development. This new routing routine would be an improvement over routing along the top of the impeding layer without adjustment for head development, because in the proposed routine, water can flow out of a local depression. In previous routines, during a wet period, water can accumulate in a local depression and can achieve unrealistically high values. To compare this newly proposed routing routine to the two routing routines currently used (routing over the impeding layer topography and routing over the soil surface topography), we used an explicit cell-to-cell routing model (Wigmosta and Lettenmaier, 1999).
An objective of this study was to determine if subsurface lateral flow at a sand–clay interface impacts tritium uptake by mixed forest vegetation. To understand and predict tritium uptake at a tritium phytoremediation site (Hitchcock et al., 2005; Rebel et al., 2005), is it sufficient to apply several one-dimensional models located over space (Rebel et al., 2005), or do we need to apply a spatially explicit, semi-three-dimensional model? To test this hypothesis, we developed and applied a water and solute balance model to simulate both vertical and lateral water and solute fluxes, including water and solute uptake by the vegetation. We ran this model with and without the lateral flow algorithm, and compared the results. The model was applied to the Savannah River Site (SRS) in South Carolina, where a forest has been irrigated with tritium-enriched water (Hitchcock et al., 2005). The spatially explicit water and solute model was used to predict the mass of tritium transpired by the vegetation in comparison to the mass of tritium that leached below the rooting zone (potentially contaminating groundwater), moved laterally offsite (potentially contaminating an uncontaminated area), or was being retained in the soil with time.
|
Materials and Methods
|
|---|
Site Description
The research was conducted at the U.S. Department of Energy's Savannah River Site, a National Environmental Research Park. A phytoremediation project was initiated in 2000 to reduce the amount of tritiated groundwater discharging to surface waters on the SRS. A small dam was installed below the discharging area to collect contaminated groundwater. Contaminated water impounded by the dam was used to irrigate 8.9 ha of the adjacent catchment, which was covered by an approximately 50-yr-old mixed pine–oak forest. Irrigation water was uniformly applied to 28 blocks, using a grid of evenly spaced sprinklers (1 per 30 m2) (Blake, 1999). The phytoremediation project is described in more detail by Hitchcock et al. (2005).
The SRS is located on the Atlantic Coastal Plain of the southeastern USA. Like much of the Atlantic Coastal Plain, the soils are gently sloping to moderately steep and have a sandy surface texture underlain at variable depths (10–250 cm) by a slowly permeable, finer textured layer (Daniels et al., 1984; Rogers and Herren, 1990). Water and tritium percolating rapidly down through the macropores of the sandy surface layer may accumulate on top of the less permeable clay layer and be routed laterally through adjacent sandy layers, while only slowly percolating in the vertical dimension through the smaller pores of the clay layers (Jury et al., 1991; Kung, 1990).
Field Measurements
Six of the 28 irrigation blocks on the phytoremediation catchment were instrumented at three locations (upslope, midslope, and downslope); these 18 locations are referred to as plots (of approximately 5 by 5 m; Fig. 1B and 1C). A number of instruments to monitor soil water content and tritium activity were installed. Each plot contained a piezometer, inserted to a depth of 3 m, to study if transient saturated conditions occurred at the study site. Depth to water in the tubes was measured monthly. As expected, transient saturated conditions were spatially and temporally variable (Fig. 2). During a wet period in 2003, saturated conditions were observed in 13 of the 18 plots. Time-domain reflectometry (TDR; Trime, MESA Systems Co., Medfield, MA) was used to measure volumetric soil water content. Two TDR tubes were placed approximately 1 m apart in each plot. The tubes were inserted to a depth of 140 to 185 cm, depending on ease of installation. Soil water content was measured biweekly at depths of 25, 55, 95, 135, and (where possible) 175 cm. To estimate residual soil tritium, each of the 18 monitoring sites included five soil vapor samplers installed at depths of 25, 55, 135, 205, and 295 cm. Vapor samplers consisted of a polyvinyl chloride well pipe screened over a 15-cm span at the base. Soil pore water was collected by pulling a small vacuum on the well pipe at the surface and condensing the water vapor (Hitchcock et al., 2005). This technique proved to be effective in yielding sufficient pore water for tritium analysis, regardless of the soil moisture conditions. Tritium activity in these water samples was analyzed using liquid scintillation analysis with an estimated detection limit of 740 Bq/L and a counting error below 2% for elevated tritium concentrations (Minaxi Tri-Carb 4000, Packard Instrument Co., Downers Grove, IL).
View larger version (42K):
[in this window]
[in a new window]
|
FIG. 1. Surface topography with (A) irrigated blocks (white lines) and (B) the six instrumented irrigation plots (white lines); and (C) interpolated map of sand layer thickness with the six instrumented irrigation plots (white lines).
|
|
Leaf area index was calculated for each plot from seasonal ceptometer (AccuPar, PAR80, Decagon Devices, Pullman, WA) measurements of direct and diffuse ray penetration, as described by Norman and Campbell (1989). Hourly measurements of radiation, air temperature, vapor density, and wind speed were taken at the SRS meteorological station. Daily precipitation was measured both at the phytoremediation site and at the meteorological tower. Pond tritium activity was measured several times a week and linearly interpolated to predict daily tritium activity in the irrigation water. Time and rate of irrigation applications to the remediation areas were reported daily.
Augering was used in and around the irrigated area to measure the depth to the loamy clay layer, which varied from 25 to 250 cm. In total, 129 locations in 140 ha of the catchment area were sampled and the depth to texture determined by feel (Brady and Weil, 1999). Block kriging, using a block size of 50 by 50 m, was then used to create a continuous depth-to-clay map for the 140-ha catchment.
|
Model Description
|
|---|
To simulate water and tritium fluxes in the vadose zone, we developed a spatially explicit water and tracer balance (WTB) model. Although this study focused on tritium fluxes, the model can be applied to simulate the transport of any nonreactive solute. The simulated catchment area is divided into square grids (5 by 5 m for this study), each with its own vegetation, soil, and topographical properties (Fig. 3). Precipitation, irrigation, depth to clay layer, depth of rooting, and leaf area index can vary among grids. The model is spatially explicit, i.e., water and tritium can be transferred among neighboring grids. Each grid consists of k cells (this study: nine, Table 1), also called voxels (Fig. 3). We used a bucket approach for simulating water fluxes, where the assumption is that water can move through the sandy part of the soil surface in one time step. The adopted daily time step was consistent with this assumption.
Water and Tritium Mass Balance
The water mass balance for each time step and each grid is
 | [1] |
where 
Sw is the change in water storage, P is precipitation, I is irrigation, Qin is lateral inflow, A is surface accumulation from the previous time step, D is drainage, AT is the actual transpiration, and Qout is lateral outflow, all in millimeters per day.
The water mass balance for each cell k in a grid is
 | [2] |
in which Sw,k is the change in water storage in a cell, Dk–1 is water flowing vertically into cell k, Dk is water flowing vertically out of cell k, Qin,k is the lateral inflow of water from the neighboring upslope cells, Qout,k is the lateral outflow of water to neighboring downslope cells, ATk is the actual root water uptake or transpiration from the cell, all in millimeters per day, and D0 is the sum of water entering the surface cell (k = 1) from rain, irrigation, and any surface accumulation from the previous time step.
When tritium enters a cell, it is assumed to mix completely with the water and tritium already present in that cell. Each of the water fluxes (in mm) described below are multiplied by tritium activity S (activity m–3) of the cell where the flux is moving from and divided by water density 
(kg m–3) to calculate the flux of tritium (activity m–2).
Water Uptake
Daily potential evapotranspiration (PET, mm d–1) is calculated using the Priestley–Taylor equation (Priestly and Taylor, 1972). The value of PET was calculated using both Priestley–Taylor and Penman–Monteith for a period of 3 yr, but on a daily basis there was very little difference between the two approaches and Priestley–Taylor was the preferred approach since it is easier to parameterize. To obtain potential transpiration (PT, mm d–1), we multiplied PET by a vegetation factor (Vf, 0–1), which is dependent on solar radiation intercepted by the canopy (Landsberg, 1986; Norman and Campbell, 1989):
 | [3] |
where k
is the extinction coefficient for diffuse ray penetration and L is leaf area index (m2 m–2).
Water uptake within a grid is proportioned among the cells containing roots within the grid, according to the fraction of roots in each cell (Rf) relative to the total roots in the grid. Water uptake (Wuk, mm d–1) within cell k can occur until a lower limit of plant-extractable water, the permanent wilting point (PWPk, m3 m–3) in cell k is reached, at which point unmet demand of cell k is added to cell k+1 (Adk+1, mm):
 | [4] |
where Wk is the current volumetric water content of cell k (m3 m–3), dk is the depth of cell k (m), is the water density (kg m–3), and the function min assures that the minimum value of either the uptake demand or the difference between W and PWP is taken as the actual water uptake. This algorithm is implemented sequentially, from top to bottom, through the root zone. The value of Wuk is then summed for all cells containing roots to obtain the actual transpiration (AT, mm d–1). If the cells containing roots cannot meet the PT demand, AT will be less than PT. This algorithm could underestimate AT if cells at the surface are above PWP and could meet the entire PT demand, while deeper cells containing roots are at PWP. This condition did not arise in the irrigated forest systems that we were simulating, but could occur where soils supporting perennial vegetation are exposed to long periods of drought followed by small or infrequent precipitation events.
Downward Movement of Water
Downward movement of water, Dk, is simulated using a "tipping bucket" or a capacitance model approach based on the CERES-maize model (Buttler and Riha, 1992; Jones and Kiniry, 1986; Riha and Rossiter, 1993, p. 102). The drainage coefficient of cell k, dck (0–1), is used to limit the amount of water above field capacity (fc, m3 m–3) that can move between cells in one time step:
 | [5] |
When W + Sw,k (Eq. [2]) is greater than saturation, water can bypass cell k and move into cell k + 1. The amount of water bypassing is the newly updated water content in layer k (WNk, m3 m–3) minus the saturated water content of layer k (WSk, m3 m–3). When cell k + 1 is a restricting cell (in this study dck+1 
0.25), however, water cannot bypass cell k, which can result in oversaturation of cell k.
Lateral Movement of Water and Tritium
The lateral drainage coefficient ldck (0–1) functions similarly to dck. If the soil is isotropic (e.g., a sandy soil), dck and ldck are equal in value. If the soil is anisotropic (e.g., a loamy clay soil), ldck is smaller than dck (Jury et al., 1991).
Lateral flow is calculated for each cell k, with the input variable Qin,k (mm d–1) being the amount of water moving into cell k from the surrounding upstream cells:
 | [6] |
where the subscript i indicates the different cells in grid space in the same layer k, and where cell k,i + 1 is the downslope neighbor of cell k,i (Fig. 3). The variable Qout,k,i is water moving laterally out of cell k,i. The direction of Qout,k,i is dependent on the local drain direction map (ldd) and the magnitude of Qout,k,i is dependent on the downstream slope (Eq. [9]). The ldd is created using the eight-point-pour algorithm with flow directions from each cell to its steepest downslope neighbor (Burrough, 1986; O'Callaghan and Mark, 1984). Both the slope (Sl, Eq. [8]) and the ldd are calculated from the topography of the restricting layer (Freer et al., 1997; McDonnell, 2003), adjusted for head development (HD, mm, Eq. [7]) by adding the accumulated water above field capacity (fc) above the restricting layer within a grid to the topography of the restricting layer (Tr, m). Both the slope and the ldd are updated daily.
 | [7] |
where rl is cell k where the restricting layer first occurs and Wk is the updated water content (m3 m–3).
 | [8] |
Lateral drainage Qout,k,i (mm d–1) is determined as the amount of water above field capacity (Shaman et al., 2002) that can drain from cell k,i to cell k,i + 1:
 | [9] |
Vertical Redistribution of Water and Tritium
If cell k is oversaturated (newly updated water content is larger than saturated water content, WNk > WSk), the difference in volumetric water content between WNk and WSk (m3 m–3) is transferred to cell k – 1, along with the tritium present in the transported water. If water reappears at the surface, it is added to Dk,0 at the next time step (A in Eq. [1]). Similarly, if tritium reappears at the surface, it will be added to the irrigation influx in the next time step.
|
Model Implementation
|
|---|
Order of Calculation within a Time Step
The model uses a daily time step. Within a time step, there is an order to the calculations. Water and tritium uptake in cells where roots are active are calculated first, then downward movement of water and tritium, followed by lateral movement of water and tritium. Lastly, we redistribute water and tritium of oversaturated cells upward.
Model Parameterization
The spatial resolution of the model is 5 by 5 m and the total area modeled is 140 ha. The simulation area is larger than the irrigation area (8.9 ha), both because the irrigation area is not contiguous and to prevent errors from boundary conditions. At the boundaries of the simulation area, water and tritium can flow out, but cannot flow in (Fig. 1A), as is reasonable given the character of the boundaries.
In a previous study, the 18 plots were simulated using a one-dimensional model (without lateral flow; Rebel et al., 2005). We adopted, for the new three-dimensional model, the soil and vegetation parameters developed in the earlier study. In the one-dimensional parameterization effort, fc for each soil layer at each site (total of 36) was estimated from soil water content data several days after large rainfall events and from late winter water content values in deeper layers. The PWP for the surface soil layer at each site was estimated from soil water content data when the soil was dry and from late summer or early fall water content values in deeper layers. The drainage coefficients dck were estimated for the 36 sites when soils were at or above fc. For each layer at each site, we manually found the dc that minimized the difference between the predicted and measured soil water content for the first year of data.
The soil for the entire simulation area was divided into nine soil layers (Table 1). The surface layer was 0.1 m and the underlying eight layers were between 0.3 and 0.4 m. Soil layers were parameterized as either sand or clay (Table 2), depending on the extrapolated depth to clay measurements (Rebel, 2005). The dc used to calculate drainage out of cell k is the numerically smallest of dck or dck+1. We assumed sand was isotropic, with dc equal to ldc, and lateral flow determined by the slope of the subsurface layer (Eq. [9]: horizontal slope, no lateral flow; vertical slope, lateral flow equals vertical flow). We assumed clay to be anisotropic, with ldc << dc (Table 2).
Root distribution and depth were estimated for the 18 sites by inverse modeling using first-year soil water content data (Rebel et al., 2005). We first assumed all sites to have a root distribution that exponentially declined with depth. Then we adjusted the root distribution, based on visual judgment of the best fit of measured and simulated soil water content for 1 yr and five depths, for each site individually, to best simulate soil water declines (Bouten et al., 1992; Musters and Bouten, 1999; Zuo and Zhang, 2002). In this implementation of the model, the root distribution was assumed to be spatially homogeneous. We used the root distribution that was most frequently determined for the 18 plots (Table 1). Jackson et al. (1996) studied global root distribution and estimated that the root distribution of all temperate and tropical trees was 26% in the top 10 cm and 60% in the top 30 cm of the soil profile, which is similar to our estimated root distribution (Table 1).
Data from a nearby weather station were used to calculate PET. The Vf was determined from seasonal measurements of leaf area index for the 18 monitoring plots. The 18 monitoring plots represented seven vegetation types (Table 3). The Vf differed between winter and summer. In summer, the average Vf for the 18 plots was 0.95 with a standard deviation of 0.02. Therefore, we used 0.95 for the entire simulation area. In winter, the average Vf for the 18 plots was 0.82 with a standard deviation of 0.08, so values were averaged for each vegetation type (Table 3) and assigned to the corresponding grids of the simulated area. Outside the irrigation plots, no detailed vegetation map was available and we assumed the average winter value for Vf (0.82).
Software
The WTB model was created using the PCRaster software, a raster-based environmental modeling system for construction of spatial dynamic simulation models (Karssenberg, 2002; de Jong, 2005). Gstat (Pebesma and Wesseling, 1998), a statistical package incorporated in PCRaster, was used to analyze the spatial autocorrelation among sample data and to generate the digital elevation model of the clay layer.
Assessing Lateral Flow
In the WTB model, we used a more complicated lateral flow routine than is currently incorporated in most water balance models: routing water over the subsurface impeding layer adjusted for head development. It is difficult to obtain knowledge of impeding layers and it costs computing time to recalculate the slope and the ldd every time step. To justify the necessity of such a flow routine, we evaluated this lateral flow routine by comparing it with two more frequently used lateral flow routines, using a flow ratio (fr). The frk is the cumulative lateral flow across time and space out of every cell within the layer k divided by the cumulative vertical flow across time and space into every cell within the layer k:
 | [10] |
In the first scenario, we calculated the slope and ldd adjusted for head development above the clay layer (as described above) in each cell. In the second scenario, we only used the depth of the clay layer to calculate the slope and ldd, and in the third scenario we used the surface topography to calculate the slope and ldd. In the third scenario the clay layer was assumed to be at a depth of 1.15 m below the surface, half way between the minimum and the maximum clay layer depth measured.
|
Results and Discussion
|
|---|
Depth to the Clay Layer
The depth of the sandy surface horizon of the soils (e.g., depth to the clay layer) at the remediation site and adjacent areas varies between 0.25 and 2.5 m. Using the 129 augering data, we created an exponential semivariogram (Fig. 4). The semivariogram has a large nugget (1170 cm2), which implies either a substantial small-scale variability or a large measurement error (Cressie, 1991), but probably a combination of the two. The range of the semivariogram is 250 m and the sill is 3500 cm2 (Shaw et al., 2001).
View larger version (12K):
[in this window]
[in a new window]
|
FIG. 4. Semivariogram for the depth of the clay layer. The number of pairwise comparisons varied from 39 at a distance of 25 m to >300 at a distance of 150 to 400 m.
|
|
Using the semivariogram and block kriging, we interpolated the measured depth to the clay layer. A section of the resulting map of depth to the clay layer for the tritium remediation site is presented in Fig. 1C.
Assessing Lateral Flow Routines
Lateral flow routines were evaluated for Soil Layers 2, 3, and 4, where lateral flow importance is described using fr (Eq. [10]). The fr reaches a quasi-steady-state value after the first few months of simulation. The values for the fr for Scenario 1 (lateral flow on top of subsurface clay topography adjusted for head development), Scenario two (lateral flow over subsurface clay topography), and Scenario three (lateral flow over surface topography) at the end of a 3-yr simulation are shown in Table 4. In the third scenario, where all clay is below Layer 4, and Layers 2, 3, and 4 are sandy, the flow ratio is 0.08 or less. In Layers 2 and 3, the flow ratio is low because vertical flow dominates over lateral flow. In Layer 4, which is underlain by the clay layer, we might expect lateral flow to be more important. At this depth, however, vertical water fluxes have been greatly reduced due to plant water uptake, so saturated conditions rarely develop. In contrast, in Scenarios 1 and 2, 30% of the cells in Layer 2, 60% in Layer 3, and 90% in Layer 4 are underlain with clay. Whether or not we include head development above the clay layer to determine slope and ldd had no impact on the fr. In both cases, the fr for Layer 3 was greater than for Layer 2 or Layer 4 (Table 4). While some cells within Layer 3 (30%) were clay, 60% of the cells underlying Layer 3 were clay. Even though 90% of the cells underlying Layer 4 were clay, the flow ratio decreased for this layer. This was because 60% of the cells in Layer 4 were clay, and therefore had a very low conductivity for lateral flow (Table 1). In Layer 3, which represents a depth of 40 to 75 cm, lateral flow relative to vertical flow was at a maximum (14% of the vertical flow).
The impeding layer topography and the soil surface topography differ, with the clay layer thickness being highly variable in the simulated area (Fig. 1B and 1C). Freer et al. (1997) concluded that, where impeding layer topography is distinctly different from surficial topography, the restrictive surface has a considerable influence on local hydrological gradients and therefore the dominant flow path directions. On the other hand, Hutchinson and Moore (2000) found that lateral flow processes were controlled by the impeding layer in dry times, but in very wet periods, lateral routing characteristics followed the surface topography. This type of behavior was not evident in our simulation, possibly because our site was characterized by a sandy soil layer in which water moved quickly to the subsurface layer.
There were no differences in the fr between the algorithms using the head development and solely the clay layer for calculating the slope and the ldd. The values for fr reported in Table 4, however, are averages for the entire simulation area. There were local differences in the flow direction due to head development resulting in the ldd changing with time. Therefore, when lateral flow was simulated in other components of our study, we used the algorithm in which the clay layer was adjusted for head buildup.
Impact of Lateral Flow on Tritium Fluxes
To evaluate the impact of lateral flow on tritium fluxes, we simulated water and tritium fluxes at the phytoremediation site for 3 yr. In those 3 yr, in addition to 4075 mm of precipitation, 1653 mm of irrigation water with an average tritium activity of 410 552 Bq/L were applied. The total tritium activity in the irrigated plots was dynamic over time and space. At the end of the 3-yr irrigation simulation, the highest simulated tritium activity in a grid was 277 Bq/m2. Tritium enriched water was predicted by WTB to move laterally away from the irrigation sites. However, the simulated tritium activity had decreased to <3.7 x 10–5 Bq/m2 in all grids further than 70 m downslope from the irrigation sites.
Shaw et al. (2001) and Blume et al. (1987) studied lateral flow in similar Coastal Plain soils in Georgia. Shaw et al. (2001) used bromide as a tracer to evaluate if argillic horizons act as impeding layers that may induce lateral flow. Analysis of bromide concentrations suggested that a very slight lateral redistribution occurred. Blume et al. (1987) also used bromide as a tracer to study lateral flow over a plinthite layer. Their results are very similar to ours; they found that lateral flow was occurring and after 1 yr they found detectable bromide 60 m downslope of the bromide injections, whereas we predicted that tritium activity would be elevated 40 m downslope from the irrigation plots in the first year of simulation.
To assess the impact of incorporating lateral flow into the WTB model on tritium activity predictions with depth, we ran the WTB model with and without implementing the lateral flow algorithm. The output of these two scenarios, the average (of 18 plots) simulated tritium activity with depth and time, along with average (of 18 plots) measured values for these plots are shown in Fig. 5 and 6. When lateral flow was included (Fig. 5), the average standard deviation of the simulated tritium activity for all layers was larger (17 131 Bq/L) than when lateral flow was not simulated (9102 Bq/L). Differences in standard deviation were greatest at or below 135 cm (Fig. 5 and 6). This was expected; using the fr, we had determined that most lateral flow occurs in Layers 3 and 4, which is between 55 and 135 cm (Table 4). On average, incorporating lateral flow did not improve the tritium predictions with depth, but locally, tritium activities differed when lateral flow was incorporated.
View larger version (33K):
[in this window]
[in a new window]
|
FIG. 5. Mean (of 18 plots) of both simulated tritium activity using the spatially distributed model with lateral flow (black line) ± SD (gray lines) and actual measured tritium activity (squares) ± SD (bars).
|
|
View larger version (31K):
[in this window]
[in a new window]
|
FIG. 6. Mean (of 18 plots) of both simulated tritium activity using the spatially distributed model without lateral flow (black line) ± SD (gray lines) and actual measured tritium activity (squares) ± SD (bars).
|
|
The error terms for the WTB model tritium predictions including and excluding lateral flow are shown in Tables 5 and 6. The mean absolute error in tritium prediction with depth is similar whether or not lateral flow was included in the WTB model (Table 5). With and without lateral flow, the mean bias in tritium prediction was always positive (measured values greater than predicted values) and was similar at 25 and 55 cm (Table 6). Below 55 cm, however, mean bias was greater when lateral flow was included (Table 6). Whether or not lateral flow was included, simulated tritium activity at 135 cm was less than measured from April to December 2003. In the simulation without lateral flow, between January 2002 and January 2003, predicted tritium activity at 135 cm was higher than measured, which was not the case when lateral flow was included (Fig. 5 and 6). Therefore, the mean bias for the entire simulation period was lower at this depth when lateral flow was not included (Tables 5 and 6). The same pattern, but less extreme, occurred at 205 and 295 cm.
View this table:
[in this window]
[in a new window]
|
TABLE 5. Mean absolute error calculated as  it|mi,t – pi,t|/nit for soil depth k, where mi,t is measured and pi,t is predicted tritium activity, t is time, i is the cells in a layer at depth k, and n is the number of cells in a layer at depth k.
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 6. Mean bias calculated as it|mi,t – pi,t|/nit for soil depth k, where mi,t is measured and pi,t is predicted tritium activity, t is time, i is the cells in a layer at depth k, and n is the number of cells in a layer at depth k.
|
|
The goal of this study was to predict the amount of tritium lost from the site by evapotranspiration relative to the amount applied in irrigation water, and to evaluate if the transpired amount of tritium differs by including lateral flow in the simulations. We defined efficiency as cumulative tritium transpired divided by cumulative tritium applied. This value reached a quasi-steady state after approximately 18 mo (Rebel, 2005). We compared the efficiency of all cells in the irrigated area simulated with and without lateral flow after 3 yr of simulation (Table 7). The spatially averaged efficiency is similar whether or not lateral flow was simulated. There are differences, however, in the spatial variability of the efficiency (Table 7). Efficiencies simulated with lateral flow were more spatially variable than efficiencies simulated without lateral flow. This indicates that while locally it may be important to consider lateral flow, to obtain an average efficiency it is more important to represent the spatial variability of the area using several one-dimensional models than to apply a complete three-dimensional model including lateral flow routines.
View this table:
[in this window]
[in a new window]
|
TABLE 7. Efficiency (tritium transpired/tritium applied) after a 3-yr simulation using the Water and Tracer Balance model with and without lateral flow across the complete irrigation area and standard deviation of efficiency across cells.
|
|
|
Conclusions
|
|---|
We applied a spatially explicit water and tritium balance model to the texture-contrasted soils at the Savannah River Site to quantify tritium evapotranspiration that resulted from a tritium phytoremediation effort. Our specific modeling objectives were twofold: to determine if knowledge of the impeding layer topography is important to simulate subsurface lateral flow and to determine if subsurface lateral flow at a sand–clay interface impacts tritium uptake by mixed forest vegetation.
We conclude that knowledge of the impeding layer topography is important for simulating spatial patterns in daily lateral flow. The fraction of water and tritium flowing laterally, relative to the water and tritium flowing vertically into a cell, decreased from a maximum of 0.14 using impeding clay layer topography, to a maximum of 0.08 when substituting the digital elevation model for the impeding layer. Spatially averaged, we did not find a difference between routing subsurface lateral flow over the impeding layer or over the impeding layer adjusted for head development, but locally, dynamic flow directions were observed when the subsurface topography was adjusted for head development.
In the simulations, lateral flow did occur. After a 3-yr simulation of water and tritium fluxes, tritium activity stayed within 70 m downslope of the area where tritium was applied. Within the phytoremediation area, simulation with and without lateral flow did not result in significantly different tritium activity with depth when averaged across the 18 monitoring plots. The standard deviation of the tritium activity with depth increased when simulation of lateral flow was included, indicating that the local variability increased. The amount of tritium taken up by the mixed forest vegetation, simulated for the complete irrigation area, does not differ whether lateral flow is incorporated or not. Again, the standard deviation of the simulated transpired tritium increased, indicating more local variability. We conclude that, at the Savannah River Site, on the southeastern U.S. Coastal Plain soils with sand overlying clay, on an area average, it is not important to include lateral flow in the simulation to predict tritium uptake by the vegetation. It is more important to capture the spatial variability in depth to clay layer and vegetation differences, which can be done using multiple one-dimensional models. Locally, it can be important to incorporate subsurface lateral flow, routed over the subsurface impeding layer and adjusted for head buildup.
|
ACKNOWLEDGMENTS
|
|---|
Acknowledgments
We would like to acknowledge John Seaman, Chris Barton, Julian Singer, Dan Hitchcock, Susan Bell, John Blake, Stephanie Smith, Robbie Williams, and Jason McRee. This research is funded by the Department of Energy-Savannah River Operations Office through the U.S. Forest Service Savannah River under Interagency Agreement DE-AI09-00SR22188 (Cooperative Agreement no. 01-CA-11083600-001), and by the Cornell NSF-RTG biogeochemistry program.
|
REFERENCES
|
|---|
- Band, L.E., P. Patterson, R. Nemani, and S.W. Running. 1993. Forest ecosystem processes at the watershed scale: Incorporating hillslope hydrology. Agric. For. Meteorol. 63:93–126.[CrossRef]
- Beven, K.J., and M.J. Kirkby. 1979. A physically based, variable contributing area model of basin hydrology. Hydrol. Sci. B 24:43–69.
- Beven, K.J., R. Lamb, P. Quinn, R. Romanowicz, and J. Freer. 1995. TOPMODEL. p. 627–668. In V.P. Singh (ed.) Computer models of watershed hydrology. Water Resour. Publ., Highlands Ranch, CO.
- Blake, J.R. 1999. Potential impacts of operating an irrigation system to enhance evapotranspiration of tritiated water as an interim action to reduce tritium flux to Fourmile Branch. SRI-99–09-R. U.S. For. Serv. Savannah River Inst., Aiken, SC.
- Blume, L.J., H.F. Perkins, and R.K. Hubbard. 1987. Subsurface water movement in an upland Coastal Plain soil as influenced by plinthite. Soil Sci. Soc. Am. J. 51:774–779.[Web of Science]
- Bouten, W., T.J. Heimovaara, and A. Tiktak. 1992. Spatial patterns of throughfall and soil water dynamics in a Douglas fir stand. Water Resour. Res. 28:3227–3233.[CrossRef]
- Brady, N.C., and R.R. Weil. 1999. The nature and properties of soils. Prentice Hall, Saddle River, NJ.
- Burrough, P.A. 1986. Principles of geographical information systems for land resources assessment. Clarendon Press, Oxford, UK.
- Burt, T.P., and D.P. Butcher. 1985. Topographic controls of soil moisture distributions. J. Soil Sci. 36:469–486.[CrossRef]
- Buttle, J.M., and D.J. McDonald. 2002. Coupled vertical and lateral preferential flow on a forested slope. Water Resour. Res. 38(5):1060, doi:10.1029/2001WR000773.[CrossRef]
- Buttler, I.W., and S.J. Riha. 1992. Water fluxes in Oxisols: A comparison of approaches. Water Resour. Res. 28:221–229.[CrossRef]
- Cressie, N.A.C. 1991. Statistics for spatial data. John Wiley & Sons, New York.
- Daniels, R.B., H.J. Kleiss, S.W. Buol, and J.A. Phillips. 1984. Soil systems in North Carolina. Bull. 467. N.C. Agric. Res. Serv., Raleigh.
- de Jong, K. 2005. PCRaster display manual. Available at pcraster.geo.uu.nl/download/doc/PCRasterManual18Feb2005.pdf (access date September 2005; verified 8 Dec. 2006). Dep. of Physical Geogr., Utrecht Univ., Utrecht, the Netherlands.
- Frankenberger, J.R., E.S. Brooks, M.T. Walter, and M.F. Walter. 1999. A GIS-based variable source area hydrology model. Hydrol. Processes 13:805–822.[CrossRef]
- Freer, J., J. McDonnell, K.J. Beven, D. Brammer, D. Burns, R.P. Hooper, and C. Kendal. 1997. Topographic controls on subsurface storm flow at the hillslope scale for two hydrologically distinct small catchments. Hydrol. Processes 11:1347–1352.[CrossRef]
- Hewlett, J.D., and A.R. Hibbert. 1967. Factors affecting the response of small watersheds to precipitation in humid areas. p. 275–291. In W.E. Sopper and H.W. Lull (ed.) Forest hydrology. Pergamon Press, New York.
- Hitchcock, D.R., C.D. Barton, K.T. Rebel, J. Singer, J.C. Seaman, J.D. Strawbridge, S.J. Riha, and J.I. Blake. 2005. A containment and disposition strategy for tritium-contaminated groundwater at the Savannah River Site, South Carolina, United States. Environ. Geosci. 12:17–28.[Abstract/Free Full Text]
- Hutchinson, D.G., and R.D. Moore. 2000. Throughflow variability on a forested hillslope underlain by compacted glacial till. Hydrol. Processes 14:1751–1766.[CrossRef]
- Jackson, R.B., J. Canadell, J.R. Ehleringer, H.A. Mooney, O.E. Sala, and E.D. Schulze. 1996. A global analysis of root distributions for terrestrial biomes. Oecologia 108:389–411.[CrossRef][Web of Science]
- Jones, C.A., and J.R. Kiniry. 1986. CERES-Maize: A simulation model of maize growth and development. Texas A&M Univ. Press, College Station.
- Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil physics. John Wiley & Sons, New York.
- Karssenberg, D. 2002. The value of environmental modelling languages for building distributed hydrological models. Hydrol. Processes 16:2751–2766.[CrossRef]
- Kung, K.-J.S. 1990. Preferential flow in a sandy vadose zone: 1. Field observation. Geoderma 46:51–71.[CrossRef][Web of Science]
- Landsberg, J.J. 1986. Physiological ecology of forest production. Academic Press, London.
- McDonnell, J.J. 2003. Where does water go when it rains? Moving beyond the variable source area concept of rainfall–runoff response. Hydrol. Processes 17:1869–1875.[CrossRef]
- Moore, R.D., and J.C. Thompson. 1996. Are water table variations in a shallow forest soil consistent with the TOPMODEL concept? Water Resour. Res. 32:663–669.[CrossRef]
- Musters, P.A.D., and W. Bouten. 1999. Assessing rooting depths of an Austrian pine stand by inverse modeling soil water content maps. Water Resour. Res. 35:3041–3048.[CrossRef]
- Norman, J.M., and G.S. Campbell. 1989. Canopy structure. p. 301–325. In R.W. Pearcy et al. (ed.) Plant physiological ecology: Field methods and instrumentation. Chapman and Hall, New York.
- O'Callaghan, J.F., and D.M. Mark. 1984. The extraction of drainage networks from digital elevation data. Comput. Vision Graph. Image Process. 28:328–344.
- Pebesma, E.J., and C.G. Wesseling. 1998. Gstat: A program for geostatistical modelling, prediction and simulation. Comput. Geosci. 24:17–31.[CrossRef]
- Priestly, C.H.B., and B.J. Taylor. 1972. On the assessment of surface heat flux and evaporation using large-scale parameters. Mon. Weather Rev. 100:81–92.[CrossRef]
- Rebel, K.T. 2005. Using trees to remediate tritium contaminated groundwater: A modeling and tracer study. Ph.D. diss. (Diss. Abstr. 3163276.) Cornell Univ., Ithaca, NY.
- Rebel, K.T., S.J. Riha, J. Seaman, and C. Barton. 2005. The use of dynamic modeling in assessing tritium phytoremediation. Environ. Geosci. 12:243–250.[Abstract/Free Full Text]
- Riha, S.J., and D.G. Rossiter. 1993. GAPS: General-purpose atmosphere–plant–soil simulator. Dep. of Soil, Crop and Atmos. Sci., Cornell Univ., Ithaca, NY.
- Rogers, V.A., and E.C. Herren. 1990. Soil survey of Savannah River Plant area, parts of Aiken, Barnwell, and Allendale Counties, South Carolina. U.S. Gov. Print. Office, Washington, DC.
- Shaman, J., M. Stieglitz, V. Engel, R. Koster, and C. Stark. 2002. Representation of subsurface storm flow and a more responsive water table in a TOPMODEL-based hydrology model. Water Resour. Res. 38(8):1156, doi:10.1029/2001WR000636.[CrossRef]
- Shaw, J.N., D.D. Bosch, L.T. West, C.C. Truman, and D.E. Radcliffe. 2001. Lateral flow in loamy to sandy Kandiudults of the Upper Coastal Plain of Georgia (USA). Geoderma 99:1–25.[CrossRef][Web of Science]
- Wigmosta, M.S., and D.P. Lettenmaier. 1999. A comparison of simplified methods for routing topographically driven subsurface flow. Water Resour. Res. 35:255–264.[CrossRef]
- Wigmosta, M.S., L.W. Vail, and D.P. Lettenmaier. 1994. A distributed hydrology–vegetation model for complex terrain. Water Resour. Res. 30:1665–1679.[CrossRef]
- Zuo, Q., and R.D. Zhang. 2002. Estimating root-water-uptake using an inverse method. Soil Sci. 167:561–571.[CrossRef]
TOP
ABSTRACT
INTRODUCTION
by soil type.