Simulation and validation of concentrated subsurface lateral flow paths in an agricultural landscape

The importance of soil water flow paths to the transport of nutrients and contaminants has long been recognized. However, effective means of detecting concentrated subsurface flow paths in a large landscape are still lacking. The flow direction and accumulation algorithm based on single-direction flow algorithm (D8) in GIS hydrologic modeling is a cost-effective way to simulate potential concentrated flow paths over a large area once relevant data are collected. This study tested the D8 algorithm for simulating concentrated lateral flow paths at three interfaces in soil profiles in a 19.5-ha agricultural landscape in central Pennsylvania, USA. These interfaces were ( 1) the interface between surface plowed layers of Ap1 and Ap2 horizons, ( 2) the interface with subsoil water-restricting clay layer where clay content increased to over 40%, and ( 3) the soil-bedrock interface. The simulated flow paths were validated through soil hydrologic monitoring, geophysical surveys, and observable soil morphological features. The results confirmed that concentrated subsurface lateral flow occurred at the interfaces with the clay layer and the underlying bedrock. At these two interfaces, the soils on the simulated flow paths were closer to saturation and showed more temporally unstable moisture dynamics than those off the simulated flow paths. Apparent electrical conductivity in the soil on the simulated flow paths was elevated and temporally unstable as compared to those outside the simulated paths. The soil cores collected from the simulated flow paths showed significantly higher Mn content at these interfaces than those away from the simulated paths. These results suggest that ( 1) he D8 algorithm is useful in simulating possible concentrated subsurface latCorrespondence to: H. S. Lin (henrylin@psu.edu) eral flow paths if used with appropriate threshold value of contributing area and sufficiently detailed digital elevation model (DEM); (2) repeated electromagnetic surveys can reflect the temporal change of soil water storage and thus is a useful indicator of possible subsurface flow path over a large area; and ( 3) observable Mn distribution in soil profiles can be used as a simple indicator of water flow paths in soils and over the landscape; however, it does require sufficient soil sampling (by excavation or augering) to possibly infer landscape-scale subsurface flow paths. In areas where subsurface interface topography varies similarly with surface topography, surface DEM can be used to simulate potential subsurface lateral flow path reasonably so the cost associated with obtaining depth to subsurface water-restricting layer can be minimized.

Abstract.The importance of soil water flow paths to the transport of nutrients and contaminants has long been recognized.However, effective means of detecting concentrated subsurface flow paths in a large landscape are still lacking.The flow direction and accumulation algorithm based on single-direction flow algorithm (D8) in GIS hydrologic modeling is a cost-effective way to simulate potential concentrated flow paths over a large area once relevant data are collected.This study tested the D8 algorithm for simulating concentrated lateral flow paths at three interfaces in soil profiles in a 19.5-ha agricultural landscape in central Pennsylvania, USA.These interfaces were (1) the interface between surface plowed layers of Ap1 and Ap2 horizons, (2) the interface with subsoil water-restricting clay layer where clay content increased to over 40%, and (3) the soil-bedrock interface.The simulated flow paths were validated through soil hydrologic monitoring, geophysical surveys, and observable soil morphological features.The results confirmed that concentrated subsurface lateral flow occurred at the interfaces with the clay layer and the underlying bedrock.At these two interfaces, the soils on the simulated flow paths were closer to saturation and showed more temporally unstable moisture dynamics than those off the simulated flow paths.Apparent electrical conductivity in the soil on the simulated flow paths was elevated and temporally unstable as compared to those outside the simulated paths.The soil cores collected from the simulated flow paths showed significantly higher Mn content at these interfaces than those away from the simulated paths.These results suggest that (1) the D8 algorithm is useful in simulating possible concentrated subsurface lat-Correspondence to: H. S. Lin (henrylin@psu.edu)eral flow paths if used with appropriate threshold value of contributing area and sufficiently detailed digital elevation model (DEM); (2) repeated electromagnetic surveys can reflect the temporal change of soil water storage and thus is a useful indicator of possible subsurface flow path over a large area; and (3) observable Mn distribution in soil profiles can be used as a simple indicator of water flow paths in soils and over the landscape; however, it does require sufficient soil sampling (by excavation or augering) to possibly infer landscape-scale subsurface flow paths.In areas where subsurface interface topography varies similarly with surface topography, surface DEM can be used to simulate potential subsurface lateral flow path reasonably so the cost associated with obtaining depth to subsurface water-restricting layer can be minimized.

Introduction
Contribution of concentrated subsurface lateral flow in soils to rapid transport of nutrients and chemicals has been well recognized (e.g., Tsukamoto and Ohta, 1988;Elliot et al., 1998).Therefore, generating three-dimensional (3-D) scheme of subsurface flow paths in a landscape can help nutrient management and pollution control.However, limited means are available for detecting (especially nondestructively) subsurface flow paths in a large landscape.In addition, most studies on concentrated subsurface lateral flow reported in the literature have been conducted in forested catchments (e.g., Kitahara et al., 1994;Sidle et al., 2001;Lin et al., 2006), with much fewer studies conducted in agricultural landscapes.
The soil-bedrock interface has been recognized in a number of recent studies as an important concentrated subsurface Q. Zhu and H. S. Lin: Simulation and validation of concentrated subsurface lateral flow lateral flow path.For example, Freer et al. (1997) reported a positive correlation between total flow volume and the contributing area calculated from a digital elevation model (DEM) of the soil-bedrock interface (instead of the soil surface).Noguchi et al. (1999) demonstrated through dye tracing that bedrock topography was important in contributing to preferential flow in a forested hillslope.Buttle and McDonald (2002) found that water flow at bedrock surface occurred in a thin saturated layer.Haga et al. (2005) demonstrated that saturated subsurface flow above the soil-bedrock interface was dominant subsurface runoff.Fiori et al. (2007) also reported that the principal mechanism for the stream flow generation was subsurface flow along the soil-bedrock interface.
Significant changes in texture, structure, and bulk density can often be observed across the boundary of two adjacent soil horizons in a soil profile.Therefore, soil horizon interface can alter water flow direction and pattern (e.g., Kung, 1990Kung, , 1993;;Ju and Kung, 1993;Gish et al., 2005;Lin et al., 2006).Several studies have reported water accumulation and subsequent lateral preferential flow above a high clay content and low hydraulic conductivity B horizon (called argillic horizon) (e.g., Haria et al., 1994;Perillo et al., 1999;Heppell et al., 2000).Slowly-permeable fragipans in many soils have been recognized to develop seasonal perched water table and thus trigger lateral preferential flow (e.g., Palkovics and Peterson, 1977;McDaniel et al., 2008).Because of compaction caused by farming equipments, plowpan (Ap2 horizon) underneath plowed layer (Ap1 horizon) can potentially generate lateral seepage, especially in rice paddy soils (e.g., Chen et al., 2002;Sander and Gerke, 2007).Sidle et al. (2001) also observed lateral flow at organic horizon-mineral soil interface in forested hillslopes.
Although concentrated subsurface lateral flow at the interfaces between soil horizons and between soil and underlying bedrock are important to water flow and chemical transport across a landscape, methods for effectively determining where and when concentrated subsurface lateral flow occurs remain very limited.In recent years, flow direction and accumulation simulations based on DEM have been implemented in Geographic Information System (GIS) hydrologic modeling tools (e.g., Maidment, 2002).These simulations are based on single-direction flow algorithm (D8) (also called steepest gradient algorithm) (O'Callaghan and Mark, 1984).The D8 algorithm has been widely used in the simulation of surface flow paths (e.g., Marks et al., 1984;Jones, 2002;Schäuble et al., 2008).However, it has not been widely used to simulate subsurface flow paths at different interfaces, although some studies have used other algorithms (e.g., MD8 and Dinf) to predict groundwater level thus groundwater flow (e.g., Zinko et al., 2005;Sørenson et al., 2006).Gish et al. (2005) have used the D8 modeling tool to identify concentrated subsurface lateral flow paths above the clay layer in an agricultural watershed, which were confirmed by ground penetration radar (GPR) investigations.Bakhsh and Kanwar (2008) reported that flow accumulation generated from the D8 method contributed significantly to discriminate subsurface drainage clusters.
However, the D8 method only allows one of eight flow directions, which constrains the representation of flow path variability (Fairfield and Leymarie, 1991).Flow path simulated using the D8 tends to be concentrated to distinct, often artificially straight lines, as reported by Seibert and McGlynn (2007).In addition, Kenny et al. (2008) pointed out that the D8 algorithm can not yield good simulation results in low relief areas or areas with poor DEMs.Efforts to alleviate these drawbacks have focused on introducing models with multiple-flow directions, also called dispersive algorithms.For example, the algorithm proposed by Quinn et al. (1991) (MD8) distributes flow to all neighboring downslope cells weighted according to slope.However, dispersive algorithms produce numerical dispersion from a DEM cell to all neighboring cells with a lower elevation, which may be inconsistent with the physical definition of upstream drainage area (Orlandini et al., 2003).Tarboton (1997) proposed a nondispersive algorithm (Dinf) that assigns flow direction angle between 0 and 2 π radian and allows an infinite number of possible flow directions.However, a certain degree of dispersion still remains in this method (Orlandini et al., 2003).Seibert and McGlynn (2007) incorporated and extened the MD8 and Dinf models and generated a new triangular multiple flow direction algorithm (MD∞) that takes the advantages of both models.According to Paik (2008), dispersive algorithms cannot define specific flow paths (nondispersive); therefore they are not suitable for investigating the transport of nutrient, pollutant, and water through channel corridors.In this respect, nondispersive algorithms (e.g., D8) are preferable.However, other studies also showed that the dispersive algorithms can be used to define dispersive flow paths (e.g., Bogaart and Troch, 2006).A few studies have suggested that the D8 method can yield good results in areas of substantial relief using a high resolution DEM (e.g., 3-5 m resolution DEM) (Guo et al., 2004;Kenny et al., 2008;Paik, 2008;Wu et al., 2008).
Soil apparent electrical conductivity (ECa) readings from electromagnetic induction (EMI) surveys are affected by soil properties such as clay content, soil water status, organic matter content, salinity, and depth to bedrock (Rhoades et al., 1976;Auerswald et al., 2001;Corwin and Lesch, 2005).Previous studies have used soil ECa from EMI surveys to represent soil water content.For example, Sherlock and Mc-Donnell (2003) reported that soil ECa from EM38 (Geonics, Mississauga, Canada) vertical dipole mode could explain over 70% of gravimetrically determined soil-water variance.Reedy and Scanlon (2003) used the same sensor to explain 80% of the averaged volumetric water content in the soil profile.Although soil ECa values are affected by many soil properties, most of them (e.g., clay content, depth to bedrock, and organic matter content) are temporally stable over a relatively short period of time.Therefore, the temporal variation Hydrol.Earth Syst.Sci., 13,[1503][1504][1505][1506][1507][1508][1509][1510][1511][1512][1513][1514][1515][1516][1517][1518]2009 www.hydrol-earth-syst-sci.net/13/1503/2009/ of soil ECa measured repeatedly within a short period of time (e.g., within a few months) should reflect the temporal variation in soil water content, after temperature correction (Zhu et al., 2009).Soil morphological features (e.g., redox features, soil structure, macropores, and many others) are indicative of soil water movement (Lin et al., 2005).For example, previous studies have suggested that soil manganese (Mn) content is a good indicator of water movement in soil profiles.This is because Mn can be easily reduced and mobilized with moving water, and then oxidized and re-deposited when soil dries and O 2 reenters the soil (Patrick and Henderson, 1981).Yaalon et al. (1977) found that soil Mn content was topography and drainage related in three catenas.McDaniel and Buol (1991) and Walker and Lin (2008) reported greater soil Mn content at footslope and concave landscape positions because of water accumulation.Cassel et al. (2002) reported a positive relationship between subsurface flow paths and dissolved Mn along hillslopes.However, such simply indicator has not been utilized in larger scale studies (e.g., catchment and farm scales) to interpret landscape-scale subsurface water movement.
The objective of this study was to investigate the reliability of high-resolution DEM-derived flow direction and accumulation algorithm (D8) implemented in ArcGIS 9.2 (ESRI, Redlands, CA, USA) for simulating concentrated subsurface lateral flow paths at three interfaces in a large agricultural landscape.The interfaces investigated included (1) the interface between the surface layers of Ap1 and Ap2 horizons, (2) the interface between the upper soil profile and subsoil clay layer where clay content increased to over 40%, and (3) the interface between soil and the underlying bedrock.Three field indicators were then used to validate the simulated flow paths, including field soil moisture monitoring, EMI surveys, and soil Mn content observation at these interfaces.

Study site
This study was conducted in an agricultural landscape typical of the valley in the Northern Appalachian Ridges and Valleys Physiographic Region in the USA (Fig. 1).The 19.5ha study area is located on The Pennsylvania State University's Kepler Farm in Rock Springs, PA.Typical crops grown on this farm are corn, soybean, and winter wheat.Elevation ranges from 373 m at the footslope in the northeastern corner to 396 m at the ridge top located in the middle portion of the field (Fig. 1).Depth to bedrock ranges from less than 0.25 m on the summit to more than 3 m on the footslope based on our field investigations.According to the second-order soil survey (Soil Survey Division Staff, 1993), five soil series have been identified in this landscape: the Hagerstown, Opequon, Murrill, Nolin, and Melvin soil series (Fig. 1).There are some transition zones among these soil series that we have identified on a refined soil map, including the Opequon-Hagerstown variant, Hagerstown-Murrill variant, Hagerstown-Nolin variant, and Nolin-Melvin variant (Fig. 1).The dominant soil series are the Hagerstown silt loam (fine, mixed, semiactive, mesic Typic Hapludalfs) and the Opequon silty clay loam (clayey, mixed, active, mesic Lithic Hapludalfs).These are well-drained soils derived from limestone residuum, with the Hagerstown solum over 1.0 m thick and the Opequon solum <0.5 m thick.The Murrill series (fine-loamy, mixed, semiactive, mesic Typic Hapludults) consists of deep, well-drained soils formed in sandstone colluvium with underlying residuum weathered from limestone.The Melvin silt loam (fine-silty, mixed, active, nonacid, mesic Fluvaquentic Endoaquepts) and the Nolin silt loam (fine-silty, mixed, active, mesic Dystric Fluventic Eutrudepts) are deep soils formed in alluvium washed from surrounding uplands with limestone lying underneath the alluvium.The Nolin series is well-drained while the adjacent Melvin series is poorly drained (closer to a nearby stream).

Subsurface flow paths simulation
The DEMs of the three interfaces (the Ap1 to Ap2 interface, the interface with the clay layer, and the soil-bedrock interface) were generated by subtracting the land surface DEM (3-m resolution) by the Ap1 horizon thickness, depth to clay layer, and depth to bedrock, respectively.The 3-m resolution surface DEM was developed from 2-m contours interpreted from aerial photos.The DEMs of the three subsurface interfaces were dominated by the variation in land surface elevation since the maximum difference in surface elevation was 23 m in the study area, while the largest difference in depth to the Ap2 horizon or the clay layer or the bedrock was <2 m.All spatial operations, including interpolations described below, were implemented using ArcGIS.
A 1.1-m long intact soil core (0.038-m in diameter) was collected from 145 soil moisture monitoring sites when we installed PVC access tubes for soil moisture monitoring in this landscape (Fig. 1).These 145 sites covered all of the landforms, soil series, and depth to bedrock ranges in the study area (Table 1).The Ap1 horizon thickness was recorded from these 145 soil cores.Seventy out of these 145 soil cores were selected for determining clay content of each horizon and depth to clay layer where clay content increased to over 40% (Table 1 and Fig. 1).The clay content was analyzed using a simplified method proposed by Kettler et al. (2001).In this method, the sand particles are filtered out from the soil slurry using the 0.053-mm mesh, the slurry passed through the 0.053-mm mesh (silt + clay) were settled for 2-6 h to separate clay (in the suspension) and silt (settled down) particles.For the Nolin and Melvin series, the horizon with >40% clay was not observed (27-29% clay in their B horizons); however, a restrictive horizon with greater density (>1.6 g/cm 3 ) was presented at the depth range of 0.6-1.0 m.For simplicity, we used the depth to this restrictive horizon in the Nolin and Melvin series (which only occupied a small portion of the overall landscape; see Fig. 1) to approximate their equivalent depth to clay layer.Ordinary kriging was used to generate the maps of the Ap1 horizon thickness and the depth to the clay layer.Both of these two soil properties had a small ratio of sample spacing over spatial correlation range (<0.5) and a medium to strong spatial structure (nugget over sill ratio <0.6).According to Kravchenko (2003) and Zhu and Lin (2009), these two soil properties can be reliably interpolated with ordinary kriging.Depth to bedrock was obtained from 77 point observations (Table 1 and Fig. 1), which showed a strong correlation with auxiliary variables such as terrain indices and ECa values obtained from EMI surveys (R 2 >0.82).According to Kravchenko and Robertson (2007) and Zhu and Lin (2009), depth to bedrock can be better interpolated with regression kriging that incorporate the information of the auxiliary variables.2).This provides a grid of flow directions from one cell to its steepest downslope using the D8 algorithm.The flow accumulation determines the accumulated water from all cells in contributing area that flow into each downslope cell.Wu et al. (2008), via comparing different DEM resolutions (10, 30, 60, 90, 150, and 200 m), found that finer DEM resolution led to more accurate D8 simulation.The 3-m resolution DEM used in this study is finer than many previous studies reported in the literature (e.g., 5-and 10-m resolution DEMs in the studies of Erskine et al., 2006 andThompson et al., 2006, respectively).
A threshold of contributing area was used to determine whether a cell was involved in a concentrated lateral flow path.A smaller threshold indicates more cells participating in the flow.Thus, a smaller threshold is better for simulating flow paths under wet conditions, while a larger threshold is better suited for dry conditions.To date, no definitive model has emerged that provides clear criteria for selecting such threshold value.Besides, concentrated flow initiation mechanisms are likely to vary depending on local characteristics of climate, geology, soils, relief, and vegetation (e.g., Kirkby, 1994;Vogt et al., 2003).In the study of Gish et al. (2002), a threshold of 100 m 2 was suggested, while in the study of Bakhsh and Kanwar (2008), a threshold of 420 m 2 was used.In our study, instead of using a single threshold, we compared three thresholds of contributing area: 1000, 500, and 100 m 2 .The output from each flow simulation in ArcGIS was a raster file, which was converted to a vector file for generating three buffer zones of 0-5, 5-10, and 10-15 m away from the simulated flow paths for comparing with the EMI survey data.
After flow path simulations, 145 monitoring sites were superimposed to determine whether a site was on or off the simulated flow paths (Fig. 2).If a monitoring site was in the cell of the simulated flow path or in the cell adjacent to the simulated path, it was considered as on the flow path; otherwise,  it was considered as off the flow path.Details about determining whether a site is on or off the simulated flow path are illustrated in Fig. 2d.Gish et al. (2005) used a similar criteria (<5 m away from the simulated flow path) to determine whether a cell was on or off predicted lateral flow paths.

Data collections for validating simulated flow paths
Three sets of data were collected in the field to validate the simulated flow paths.These were: (1) soil moisture monitoring (including volumetric soil water content and matric potential), (2) EMI surveys, and (3) observable soil Mn content at the three interfaces (Table 2).Soil water content at multiple depths was monitored at 145 locations distributed throughout the study area (Fig. 1).Our procedure followed that of Lin et al. (2006).Briefly, a portable TRIME-FM Time Domain Reflectomery (TDR) Tube Probe (IMKO, Ettlingen, Germany) was used to determine volumetric soil water content while being placed at specific depth interval in a PVC access tube installed at each site.Readings were taken at six depth intervals of 0-0.2, 0.1-0.3,0.3-0.5, 0.5-0.7,0.7-0.9, and 0.9-1.1 m (representing soil water content at 0.1, 0.2, 0.4, 0.6, 0.8, and 1.0 m depth, respectively).If the depth to bedrock at a monitoring site was not sufficiently deep to allow all six depth measurements, fewer readings were taken.For example, the actual numbers of subsoil water content observations were 110 and 96 for the 0.3-0.5 m and 0.7-0.9m depth intervals, respectively.Entire study area's soil water content at these 145 sites was collected for 12 times from 2005 to 2007 (Table 3).
Seventy-four out of these 145 monitoring sites also had tensiometers installed (Fig. 1).These 74 locations were selected based on landforms and soil series in the study area (Table 1).Nested tensiometers were installed at five depths of 0.1, 0.2, 0.4, 0.8, and 1.0 m at each of these 74 sites.They were 0.15 m away from the TDR access tubes.Soil matric potentials along with soil water contents at these 74 locations were measured for 14 times from 2006 to 2008 (Table 3).
Seventy out of these 145 soil cores were selected for profile description (Fig. 1), including the quantity and size of visible pores, roots, and Mn mottles in each horizon (including that at the three interfaces).These morphological features were estimated visually following the standard soil survey procedures (Soil Survey Division Staff, 1993).
Soil ECa values were collected with EM38 sensor on four dates of 16 January, 10 March, 30 April, and 4 June of 2008.The EM38 sensor operated at a frequency of 13.2 KHz and provided effective theoretical measurement depths of 1.5 m when operated in vertical dipole mode.These EMI surveys were conducted with the same sample spacing of about 3×8 m (3-m spacing between two consecutive readings along each traverse line, and 8-m apart between traverse lines across the study area).We assumed that only soil moisture was changed during the period of our EMI surveys from January to June 2008 while other soil properties (e.g., clay content, depth to bedrock, and organic matter content) remained unchanged.Although temperature was also a changing factor, all EMI readings were corrected to a standard temperature of 25 • C. Ordinary kriging was used to generate the ECa maps for the entire study area based on its spatial structure (Kravchenko, 2003;Zhu and Lin, 2009).

Data analysis
For the 74 sites with both soil water content and matric potential measurements, soil water retention curves (SWRC) at different depths in each site were fitted with the model of van Genuchten (1980): where α and n are empirical parameters; θ s and θ r are saturation and residual water contents, respectively; h is matric potential and θ(h) is volumetric water content under h.Examples of fitted SWRC for typical soil series, texture classes, and horizons in the study area are shown in Fig. 3. Texture class was one of the main factors affecting the shape of the SWRC.For example, in the Ap horizon, as the texture class changed from silty clay loam to silt loam, θ s decreased from 0.43 to 0.37 m 3 m −3 (Fig. 3a, d, g).For the same texture class (e.g., silt loam) and soil series (e.g., the Murrill), θ s decreased from 0.37 to 0.32 m 3 m −3 as the soil horizon changed from Ap to Bt2 (Fig. 3g, h, i).Because of plowing and root growth, surface soils had lower bulk density, more pore space, and thus greater θ s than subsurface soils.Volumetric soil water contents at the field capacity (0.33 kPa) and saturation (0 kPa) for each depth and each site were estimated through the fitted curve.The estimated water contents at field capacity ranged from 0.2 to 0.3 m 3 m −3 , while the estimated water contents at saturation ranged from 0.35 to 0.45 m 3 m −3 for the entire study area.The ratio of field capacity over saturation ranged from 0.60 to 0.65 (with a mean of 0.63).These estimated soil water contents at field capacity and saturation of all depths and all 74 sites were grouped according to their soil horizons (Ap, Bt1, and Bt2)  and textural classes (silt loam, silty clay loam, and silty clay) (Fig. 4).The difference between Ap1 and Ap2 horizons was not considered here for two reasons: (1) the Ap1 horizon varied in thickness from 0.08 to 0.13-m in the study area; therefore, tensiometers installed at 0.1-m depth were located in the transition zones of Ap1 to Ap2 horizons; and (2) the TRIME-FM TDR probe was 0.18-m long, thus soil water content of the Ap1 (0-0.1 m below the ground surface) and Ap2 (0.1-0.3 m) could not be clearly separated.Although the Ap2 horizon was denser than the Ap1 horizon, such density contrast was less strong as compared to that between Ap and Bt horizons in our study area.
For volumetric soil water content collected at each specific depth of each monitoring site, relative degree of saturation (RS) was calculated by dividing volumetric soil water content by the estimated saturated water content of this horizon and textural class (Fig. 4).A RS close to 1 suggests near saturation.At a specific soil horizon interface, 95% confidence intervals for the RS values between sites on and off the simulated flow paths were calculated using SAS (SAS Institute Inc., Cary, NC, USA).These confidence intervals were then compared using one-way ANOVA to determine whether significant differences in the RS existed between sites on Q. Zhu and H. S. Lin: Simulation and validation of concentrated subsurface lateral flow and off the simulated flow paths (Table 2).Similarly, 95% confidence intervals of the RS values at, below, and above a specific interface (Ap1-Ap2, clay layer, or soil-bedrock) were also compared to determine whether significant statistical differences existed (Table 2).
At each interface, temporal stability of the RS values of all 145 monitoring sites was analyzed using the approach proposed by Vachaud et al. (1985): where R j is the arithmetic mean of RS at all sites in day j ; R ij is the RS of a particular interface at site i in day j ; N is the number of monitoring sites (in this study, N = 145); δ i is the arithmetic mean of the relative difference of RS at site i; M is the number of times that the whole study area's soil water content was measured (in this study, M = 12); and S δ is the standard deviation of δ i .Positive or negative δ i suggests that at a particular interface, site i is wetter or dryer, respectively, than the average condition of the entire study area.
The S δ depicts the magnitude of temporal stability of RS at a particular interface at site i.Higher S δ indicates a more dynamic change (i.e., temporally unstable) in soil moisture.Temporal stability of the RS values between sites on and off the simulated flow paths was also compared at each interface (Table 2) Following the same procedure, we conducted the temporal stability analysis of ECa values collected from four EMI surveys.Temporal stability of ECa in the three buffer zones (0-5, 5-10, and 10-15 m away from the simulated flow paths) and the rest of the study area were statistically compared with each other through t-test in SAS (p<0.05)(Table 2).The temporal changes in ECa reflected the change in soil water content.Therefore, the magnitude of ECa temporal stability can represent the degree of change in soil water storage in the three buffer zones of the predicted flow paths vs. the rest of the study area.
At the three interfaces, Mn contents estimated from soil cores were also statistically compared between sites on and off the simulated paths through t-test in SAS (p<0.05)(Table 2).The Mn mass observed in the soil is an indicator of soil water movement as demonstrated in other studies (e.g., McDaniel et al., 2008;Walker and Lin, 2008).

Simulated subsurface flow paths
The spatial patterns of potential lateral flow paths at the three interfaces simulated with different thresholds of contributing area (i.e., 100, 500, and 1000 m 2 ) are illustrated in Fig. 2a-c for the soil-bedrock interface, where the patterns using the thresholds of 1000 and 500 m 2 were close to each other but quite different from that using the threshold of 100 m 2 .Because of the topography of the study area, very few locations (<8% of the entire area) had a contribution area > 500 m 2 .Therefore, the simulated flow paths using the threshold contribution areas of 1000 and 500 m 2 were sparse and similar.In comparison, 25% cells of the entire study area had contribution area >100 m 2 .In the subsequent analysis, we focus on comparing the simulated flow paths obtained with 500 and 100 m 2 thresholds.In Fig. 2a-c, water moved laterally across the landscape through the soil-bedrock interface in three main areas: the north-east corner, the mid-west depressional area, and the mid-south portion.When using a smaller threshold (100 m 2 ), more areas participated in the flow, leading to 61% of the soil water monitoring sites (total 88 sites) being identified as on the simulated flow paths (Fig. 2c).In contrast, during drier conditions using a larger threshold of 500 m 2 , only 35% of the 145 monitoring sites were identified as on the simulated flow paths (Fig. 2b).

Validating the simulated flow paths through soil hydrologic monitoring
At each of the three interfaces, the RS values between the monitoring sites on and off the simulated flow paths were statistically compared (Table 2).In relatively dry conditions (average volumetric soil water content − θ at the interfaces with the clay layer or bedrock was smaller than 0.28 and 0.31 m 3 m −3 , respectively), the RS values of the sites on the simulated flow paths (500 m 2 threshold) were significantly greater (p<0.05)than those sites off the paths (Fig. 5b,  c).However, this was not the case at the Ap1-Ap2 interface (Fig. 5a).In relatively wet conditions ( − θ > 0.28 and 0.31 m 3 m −3 at the interfaces with the clay layer or bedrock, respectively), significant difference (p<0.05) in the RS values also existed between the sites on and off the simulated flow paths (100 m 2 threshold) at the interfaces with the clay layer or bedrock (Fig. 5b, c), but not at the Ap1-Ap2 interface (Fig. 5a).
To further verify the water accumulation at the clay layer interface, the RS values at this interface were statistically compared with the RS values right below this interface and 0.2-m above this interface (Fig. 6a).In relatively dry conditions ( − θ <0.28 m 3 m −3 ), the RS values at this interface for the sites on the simulated flow paths (500 m 2 threshold) were significantly greater (p <0.05) than that above and below it.In 32  relatively wet conditions ( − θ >0.28 m 3 m −3 ), such significant difference (p<0.05) in the RS values also existed between the sites on and off the flow paths simulated with the threshold of 100 m 2 .For the sites off the simulated flow paths, differences in the RS values between the clay layer interface and that above or below the interface were not significant under either dry or wet conditions (Fig. 6a).
For sites on the flow paths simulated with the threshold of 100 m 2 , the RS values at the soil-bedrock interface were significantly greater (p<0.05)than those 0.2-m above it, regardless of the wetness condition (Fig. 6b).In comparison, for sites on the flow paths simulated with the threshold 500 m 2 , significant difference in the RS values between the soil-bedrock interface and 0.2-m above it could only be observed in relatively dry conditions (    Sites with high standard deviations of the relative difference were temporally unstable (i.e., more dynamic).
The temporal stability of the RS values between the sites on and off the simulated flow paths (100 m 2 ) at all three interfaces were compared in Fig. 7.At the clay layer and soil-bedrock interfaces, sites on the simulated flow paths had greater relative differences of RS and standard deviations than the sites off the flow paths (Fig. 7b, c), suggesting a more dynamic (and thus unstable) moisture status over time for the sites on the flow paths.At the clay layer interface, 70% of the sites on the flow paths had δ i >0 and 75% of them had standard deviation of δ i >0.1, while only 13% and 36% of the sites off the flow paths had δ i >0 and standard deviation of δ i >0.1, respectively.At the soil-bedrock interface, percentages of the sites on the flow paths with positive δ i value and high standard deviation (>0.1) were even greater, 85% and 90%, respectively, while the corresponding percentages were only 12% and 36% for the sites off the flow paths.At the Ap1-Ap2 interface, the sites on and off the simulated flow paths had no distinct differences in δ i (Fig. 7a).When the flow path threshold was changed to 500 m 2 , results similar to that shown in Fig. 7 were observed (data not shown).

Validating the simulated flow paths through repeated EMI surveys
The temporal stability of ECa values in different buffer zones of the simulated flow paths is shown in Fig. 8, with a threshold of 500 m 2 .As the distance from the simulated flow paths increased, both the mean and standard deviation of the relative difference in ECa decreased, suggesting that the soils  closer to the simulated flow paths tended to have elevated and more dynamic ECa values.Additionally, the positive relative differences in ECa in the three buffer zones (0-5, 5-10, and 10-15 m) imply that their ECa values were greater than the overall average of the entire landscape.The greater standard deviations further suggest that the soil ECa values in areas closer to the simulated flow paths had higher temporal variations than the soils further away from these flow paths (Fig. 8).

Validating the simulated flow paths through soil Mn distribution
The simulated flow paths were further validated through the spatial variation in soil Mn contents at the clay layer and soilbedrock interfaces (Table 2 and Fig. 9).For the sites on the simulated flow paths (threshold of 500 m 2 ), soils Mn content at these two interfaces were generally greater than 1% and reached as high as 5-10% at some locations.For the sites off the simulated flow paths, almost no Mn was observed at these two interfaces.Therefore, locations with greater soil Mn content are expected to be on or closer to subsurface flow paths in the landscape studied.

Discussion
Our hydrological monitoring suggested that concentrated lateral flow paths at the interfaces with the clay layer or the bedrock were reasonably simulated with the D8 algorithm.First, soils on the simulated flow paths at these two interfaces were closer to saturation than those off the paths (Fig. 5b, c).Second, for sites on the simulated flow paths, soils at these two interfaces were also closer to saturation than soils below or above these interfaces (Fig. 6).Third, the temporal stability analysis indicated that the interfaces with the clay layer or the bedrock on the simulated flow paths were generally wetter and temporally less stable (Fig. 7b, c), implying generally more water movement through these interfaces.De Lannoy et al. ( 2006) also documented that concentrated subsurface lateral flow resulted in high temporal variability of soil water content at the clay layer interface.Our hydrological monitoring also indicated that lateral subsurface flow did not exist at the Ap1-Ap2 interface.No significant difference in the RS values was observed between the sites on and off the simulated flow paths at the Ap1-Ap2 interface (Fig. 5a).In addition, the temporal stability between sites on and off the simulated flow paths had no distinct differences at this interface (Fig. 7a).The interfaces of different soil horizons can trigger lateral subsurface flow as reported in some previous studies (e.g., Kung, 1990Kung, , 1993;;Ju and Kung, 1993;Gish et al., 2005).However, as shown in this study, not all interfaces are effective in generating lateral subsurface water flow in agricultural landscapes.When describing soil cores in this study, we noticed that the surface horizons (e.g., Ap1 and Ap2) were biologically active and had many macrospores (earthworm holes and root channels).These macrospores could have transported water from the Ap1 to Ap2 horizon without much restriction during wet periods and thus prevented water from accumulating at the Ap1-Ap2 interface.Another possible reason for the lack of lateral flow between the Ap1 and Ap2 horizons was the possible limitation of the monitoring devices used to determine soil moisture in these two horizons (i.e., the 0.18-m long TDR probe overlapped the Ap1 and Ap2 horizons' readings and some tensiometers were located in the transition zones between the Ap1 to Ap2 horizons).In comparison, clayenriched argillic horizon or dense fragipan have been widely recognized as important lateral flow paths in subsoils (e.g., Heppell et al., 2000;Gish et al., 2005;Lin et al., 2008).The clay layer interface was generally deeper than 0.4 m in our study area.At such depth, fewer macrospores were observed and thus the vertical water percolation could be more restricted.
The repeated EMI surveys suggested that the spatial pattern of subsurface lateral flow paths simulated with GIS was reasonable and compared favorably with the ECa data.We assumed that the change in soil water content was the main control of the temporal variation in ECa values during our repeated EMI surveys from January to June 2008 (note that all temperatures were corrected to a standard value).This is consistent with some previous studies that reported strong correlation between soil ECa and soil water storage (e.g., Sherlock and McDonnell, 2003;Reedy and Scalon, 2003).Thus, the higher and less stable ECa values in areas closer to the simulated flow paths (Fig. 8) suggest that the wetness there was increased and moisture dynamics was more significant.However, because we used approximately 8-m line spacing in our EMI surveys, specific subsurface lateral flow paths could not be clearly identified on the EMI maps.To do so would require a denser line spacing (e.g., ≤1 m) and perhaps also shorter time intervals (e.g., right before and after a large rainstorm) for repeating EMI surveys.
Soil morphological features (Mn distribution in soils) further justified the existence of subsurface lateral flow at the interfaces with the clay layer and the bedrock as such flow paths simulated with GIS matched with the observed Mn content distribution in the soils in the study area.Previous studies have suggested that soil Mn content is a good indicator of water movement in soil profiles (e.g., Yaalon et al., 1977;McDaniel and Buol, 1991;Cassel et al., 2002;Walker and Lin, 2008).In our study, high Mn content was observed Hydrol.Earth Syst.Sci., 13,[1503][1504][1505][1506][1507][1508][1509][1510][1511][1512][1513][1514][1515][1516][1517][1518]2009 www.hydrol-earth-syst-sci.net/13/1503/2009/ at the interfaces with the clay layer and the bedrock on the simulated flow paths, whereas soil Mn content off the simulated flow paths was almost zero (Fig. 9).During the relatively dry period, more drainage areas were required to initiate the lateral subsurface flow upon rainfall inputs.Thus, a larger threshold (e.g., 500 m 2 ) simulated better potential lateral flow paths at the clay layer interface (Fig. 5b).In contrast, during the wet period, a smaller threshold (e.g., 100 m 2 ) was better to simulate potential lateral flow paths at this interface (Fig. 5b).In our study area, the soil water content generally increased with depth and often reached the highest value at the soil-bedrock interface.Even during relatively dry period, the soil at the soil-bedrock interface might still be wet and free water lateral movement could occur after large rainstorms.In such case, a smaller threshold of 100 m 2 could still reasonably simulate concentrated lateral subsurface flow paths at the soil-bedrock interface.This is supported by significant differences in the RS between the sites on and off the simulated flow paths (threshold 100 m 2 ) during the dry period, as shown in Fig. 5c.In addition, significant differences between soil water content at the soilbedrock interface and 0.2 m above it can also be observed for sites on the simulated flow paths (threshold 100 m 2 ) during the dry period (Fig. 6b).This further suggests that a smaller threshold (100 m 2 ) works better to simulate flow paths at the soil-bedrock interface under both dry and wet conditions.
It is interesting to note that the flow paths simulated with the DEMs of the subsurface interfaces and the DEM of the land surface were nearly identical in this study, with over 90% of the simulated flow paths being the same (Fig. 10).This suggested that the flow path simulations were not improved much with the consideration of the subsurface interface DEMs in this particular landscape.This was attributed to two reasons: First, the topography of the interfaces was dominated by the variation in land surface elevation.In our study area, the maximum difference in surface elevation was 23 m between the lowest point in the footslope and the highest point at the ridge top.However, the largest differences in the Ap1 horizon thickness and the depth to the clay layer or the bedrock for the entire landscape were less than 2 m (i.e., <8.7% of the surface elevation change).Second, the fine scale variation of Ap1 horizon thickness and the depth to the clay layer or the bedrock might not have been well captured through our 145 or less point-based observations and their spatial interpolations over the entire landscape of 19.5-ha in size.Similar finding was reported by Birkhead et al. (1996), in which the bedrock topography derived from GPR image was shown to be closely related to the surface topography (elevations of bedrock and ground surface decreased simultaneously for about 1.5 m in a 90-m transect).
However, other studies observed that the topography of subsurface interfaces could be quite different from that of land surface and thus simulation based on subsurface interface topography would yield better results.For example, Freer et al. (2002) found that terrain attributes of land sur- face did not capture their observed spatial patterns of hillslope (30×60 m and 34 • average slope) hydrologic response, while bedrock surface topography significantly improved the interpretation of flow spatial variation.Burns et al. (1998) documented that accumulated area simulated with bedrock surface DEM explained better the base cation concentrations in the subsurface flow along a 20×50-m hillslope.These studies have a common feature, that is, intensive depth to bedrock surveys were made in a comparatively small area (e.g., >180 observations within an area of <0.2 ha).Thus, fine-scale variation in depth to bedrock was reasonably obtained in these studies.
Consequently, for areas where subsurface interface topography is dominated by surface DEM (e.g., our study area), surface DEM can be sufficiently used to simulate the subsurface concentrated lateral flow paths.Otherwise, simulation based on subsurface interface DEM is more desirable.Hence, it is important to find ways to predict depth to subsurface water-restricting layers (such as the clay layer or the bedrock) so that the cost-effectiveness of the GIS modeling could be better realized.
Through validation by soil hydrologic monitoring, EMI surveys, and soil morphological observations, it was apparent that concentrated subsurface lateral flow occurred at the interfaces with the clay layer (or water-restrictive layer) and the underlying bedrock in the agricultural landscape studied, but not at the interface between the surface plowed layers of Ap1 and Ap2 horizons (because of considerable biological activities at that interface and/or the likely limitation of the monitoring devices for clearly separating these two horizons' soil moisture dynamics).The ArcGIS hydrologic modeling (the D8 algorithm) did a reasonable job in simulating potential concentrated lateral flow paths at the interfaces in soil profiles.Such simulated subsurface lateral flow paths were temporally dynamic as they varied with the wetness condition of the landscape.Hence, using different thresholds of contributing area for the GIS hydrologic simulation would be needed to obtain expected results under different moisture conditions (e.g., 500 m 2 for relatively dry conditions and 100 m 2 for relatively wet conditions in this study).A sufficiently detailed DEM is also needed to ensure that the GIS flow algorithm performs with lower uncertainty.We suggest additional testing of this cost-effective means of predicting likely subsurface concentrated flow paths in other landscapes in order to establish a solid protocol for simulating subsurface hydrologic flow paths in different watersheds.However, some costs would incur in obtaining necessary data for such simulation to be effective, mainly in areas where subsurface interface topography differs significantly from land surface topography.Thus, finding means of predicting depth to subsurface water-restricting layers (such as the clay layer or the bedrock) will enhance the cost-effectiveness of the GIS modeling.In areas where subsurface interface topography does vary similarly with surface topography, surface DEM can be used to approximate subsurface interface topography to obtain similar results.Repeated EMI surveys provide another low-cost and nondestructive means of detecting potential concentrated subsurface flow paths; however, dense sample spacing and frequently repeated surveys would be needed if specific subsurface flow paths are to be identified via repeated EMI surveys.While soil morphological features such as Mn distribution in soil profiles also serve as useful simple indicators of subsurface water flow paths, it does require soil sampling (by excavation or augering) with sufficient numbers of observations to possibly infer landscape-scale subsurface flow paths.

Figure 1 .
Figure 1.The study area and the spatial distribution of monitoring and observation sites for soil moisture (TDR), matric potential (tensiometers), soil cores, and depth to bedrock (bedrock observation) at the Kepler Farm located in central Pennsylvania, USA.The inset at the lower right corner shows a 3D rendering of the study area.The background map is soil series distribution.

Fig. 1 .
Fig. 1.The study area and the spatial distribution of monitoring and observation sites for soil moisture (TDR), matric potential (tensiometers), soil cores, and depth to bedrock (bedrock observation) at the Kepler Farm located in central Pennsylvania, USA.The inset at the lower right corner shows a 3-D rendering of the study area.The background map is soil series distribution.

•
Sites on vs. off the simulated flow paths • Soils at vs. above or below the interface • Temporal stability of soil ECa values • Three buffer zones (0-5, 5-10 and 10-15 m) away from the simulated flow paths vs. the rest of study area • Sites on vs. off the simulated flow paths Potential concentrated lateral flow paths at the interfaces of Ap1-Ap2, clay layer, and soil-bedrock were simulated using the flow direction and accumulation algorithm implemented in ArcGIS 9.2 hydrologic modeling tool (Table

29Figure 2 .
Figure 2. Simulated concentrated lateral flow paths at the soil-bedrock interface using three thresholds of contribution area: (a) 1,000 m 2 , (b) 500 m 2 , and (c) 100 m 2 .An illustration of determining whether a soil moisture monitoring site is on or off the simulated flow path is shown in (d) (see the text for details).

Fig. 2 .
Fig. 2. Simulated concentrated lateral flow paths at the soil-bedrock interface using three thresholds of contribution area: (a) 1000 m 2 , (b) 500 m 2 , and (c) 100 m 2 .An illustration of determining whether a soil moisture monitoring site is on or off the simulated flow path is shown in (d) (see the text for details).

Figure 4 .
Figure 4. Average volumetric soil water contents and their standard deviations at field capacity and saturation based on field observed soil water content and matrix potential for different textural classes and soil horizons in the study area.

Fig. 4 .
Fig. 4. Average volumetric soil water contents and their standard deviations at field capacity and saturation based on field observed soil water content and matrix potential for different textural classes and soil horizons in the study area.

Figure 5 .
Figure 5.Comparison of the means and 95% confidence intervals of relative saturation (RS) at the monitoring sites on and off the simulated flow paths with thresholds of 500 m 2 and 100 m 2 for (a) the Ap1-Ap2 interface, (b) the clay layer interface, and (c) the soil-bedrock interface.Dash lines separate the relatively dry and wet conditions.Bars labeled with asterisks (*) indicate statistically significant difference at p<0.05 level between sites on and off the simulated paths.

Fig. 5 .
Fig. 5. Comparison of the means and 95% confidence intervals of relative saturation (RS) at the monitoring sites on and off the simulated flow paths with thresholds of 500 m 2 and 100 m 2 for (a) the Ap1-Ap2 interface, (b) the clay layer interface, and (c) the soil-bedrock interface.Dash lines separate the relatively dry and wet conditions.Bars labeled with asterisks (*) indicate statistically significant difference at p<0.05 level between sites on and off the simulated paths.

−Figure 6 .
Figure 6.Comparison of the means and 95% confidence intervals of the relative saturation (RS) at the monitoring sites on and off the simulated flow paths with thresholds of 500 m 2 and 100 m 2 for (a) the interface with the clay layer and (b) the interface with the bedrock.Within each graph, the comparison is for RS just above or below the specified interface and 0.2 m above the interface.Dash lines separate the relatively dry and wet conditions.Bars labeled with asterisks (*) indicate statistically significant difference at p<0.05 level.

Fig. 6 .
Fig. 6.Comparison of the means and 95% confidence intervals of the relative saturation (RS) at the monitoring sites on and off the simulated flow paths with thresholds of 500 m 2 and 100 m 2 for (a) the interface with the clay layer and (b) the interface with the bedrock.Within each graph, the comparison is for RS just above or below the specified interface and 0.2 m above the interface.Dash lines separate the relatively dry and wet conditions.Bars labeled with asterisks (*) indicate statistically significant difference at p<0.05 level.

Fig. 7 .
Fig. 7. Temporal stability of relative saturation (RS) at the monitoring sites on and off the simulated flow paths (using the threshold of 100 m 2 ) at (a) the Ap1-Ap2 interface, (b) the clay layer interface, and (c) the soil-bedrock interface.Sites with positive relative difference were wetter than the overall average of the entire study area.Sites with high standard deviations of the relative difference were temporally unstable (i.e., more dynamic).

36Figure 8 .
Figure 8. Temporal stability of apparent electrical conductivity (ECa) in areas with different distances away from the simulated flow paths (using the threshold of 500 m 2 ).Areas with positive relative difference in ECa indicate a higher ECa value than the overall mean of the entire study area.Areas with high standard deviation of the relative difference in ECa indicate temporally more dynamic (or unstable).Bars with the same letter are not statistically significantly different from each other at p<0.05 level.

Fig. 8 .
Fig. 8. Temporal stability of apparent electrical conductivity (ECa) in areas with different distances away from the simulated flow paths (using the threshold of 500 m 2 ).Areas with positive relative difference in ECa indicate a higher ECa value than the overall mean of the entire study area.Areas with high standard deviation of the relative difference in ECa indicate temporally more dynamic (or unstable).Bars with the same letter are not statistically significantly different from each other at p<0.05 level.

ervedFig. 9 .
Fig. 9. Observed Mn contents in the soil profiles at (a) the clay layer interface and (b) the soil-bedrock interface in relation to the simulated concentrated flow paths (using a threshold of 500 m 2 ).In the insets, Mn contents on and off the simulated flow paths are compared.Bars with different letters are statistically significantly different at p<0.05 level.

38Figure 10 .Fig. 10 .
Figure 10.Comparison of concentrated flow paths simulated using surface DEM with those using a) the clay layer interface DEM and b) the soil-bedrock interface DEM.The threshold used in this figure is 500 m 2 .

Table 1 .
Distribution of the number of soil moisture monitoring sites and soil core descriptions among different slope classes, depth to bedrock ranges, and soil series in the study area.

Table 2 .
Methods used to simulate and validate the concentrated subsurface lateral flow paths in the study area.D8: deterministic 8 single-flow algorithm; RS: relative degree to saturation; EMI: electromagnetic induction; ECa: apparent electrical conductivity.

Table 3 .
The time table of soil water content and matric potential data collections in this study.