# Characterizing mercury concentrations and fluxes in a Coastal Plain watershed: Insights from dynamic modeling and data

H. E. Golden, 1 

### C. D. Knightes, 2 

### P. A. Conrads, 3 

### G. M. Davis, 2 

### T. D. Feaster, 4 

### C. A. Journey, 3 

### S. T. Benedict, 4 

### M. E. Brigham, 5 

### and P. M. Bradley 3

Received 8 July 2011; revised 9 November 2011; accepted 11 November 2011; published 26 January 2012. 

### [ 1 

### ] 

### Mercury (Hg) is one of the leading water quality concerns in surface waters of the United States. Although watershed-scale Hg cycling research has increased in the past two decades, advances in modeling watershed Hg processes in diverse physiographic regions, spatial scales, and land cover types are needed. The goal of this study was to assess Hg cycling in a Coastal Plain system using concentrations and fluxes estimated by multiple watershed-scale models with distinct mathematical frameworks reflecting different system dynamics. We simulated total mercury (Hg T 

### , the sum of filtered and particulate forms) concentrations and fluxes from a Coastal Plain watershed (McTier Creek) using three watershed Hg models and an empirical load model. Model output was compared with observed in-stream Hg T 

### . We found that shallow subsurface flow is a potentially important transport mechanism of particulate Hg T 

### during periods when connectivity between the uplands and surface waters is maximized. Other processes (e.g., stream bank erosion, sediment re-suspension) may increase particulate Hg T 

### in the water column. Simulations and data suggest that variable source area (VSA) flow and lack of rainfall interactions with surface soil horizons result in increased dissolved Hg T 

### concentrations unrelated to DOC mobilization following precipitation events. Although flushing of DOC-Hg T 

### complexes from surface soils can also occur during this period, DOC-complexed Hg T 

### becomes more important during base flow conditions. TOPLOAD simulations highlight saturated subsurface flow as a primary driver of daily Hg T 

### loadings, but shallow subsurface flow is important for Hg T 

### loads during high-flow events. Results suggest limited seasonal trends in Hg T 

### dynamics.

Citation: Golden, H. E., C. D. Knightes, P. A. Conrads, G. M. Davis, T. D. Feaster, C. A. Journey, S. T. Benedict, M. E. Brigham, and P. M. Bradley (2012), Characterizing mercury concentrations and fluxes in a Coastal Plain watershed: Insights from dynamic modeling and data, J. Geophys. Res.

, 117, G01006, doi:10.1029/2011JG001806.



### 1. Introduction

[ 2 

### ] Mercury (Hg) is one of the leading water quality concerns in surface waters of the United States [ U.S. Environmental Protection Agency , 2005]. Atmospheric deposition is a primary source of Hg to watersheds. As such,transformations of atmospherically derived Hg on the landscape and subsequent fluxes to surface waters are wellrecognized drivers of Hg concentrations in aquatic habitats [ Bradley et al. , 2011; Brigham et al. , 2009; Glover et al.

, 2010; Grigal, 2002; Hurley et al. , 1995]. The consequences of increased Hg in surface waters results from microbial communities in anoxic environments (e.g., benthic sediments and wetlands) that methylate inorganic Hg to producemethylmercury (MeHg), a potent neurotoxin that bioaccumulates and biomagnifies in aquatic food webs [Bloom

, 1992; Krabbenhoft et al. , 1999], and subsequent ingestion of Hg-contaminated fish by humans [ Clarkson and Magos

, 2006; Mergler et al. , 2007] and piscivorous wildlife [e.g., Burgess and Meyer, 2008; Evers et al. , 2008].

[ 3 

### ] The number of watershed-scale Hg cycling studies has substantially increased in the past two decades, providing important insights on Hg inputs, outputs, and processes through empirically based approaches. Previous studiesfocused primarily on small forested catchments of northern climates or boreal regions [ Grigal, 2002] or specific catchment-scale or surface water processes including, for example, Hg mobilization and watershed transport [Bushey et al. , 2008; Demers et al. , 2010; Mason and Sullivan

, 1998] and in-stream cycling [ Brigham et al. , 2009]. Recent trends in Hg exposure science and risk management, however, emphasize a comprehensive understanding of Hg cycling 1 Ecological Exposure Research Division, Office of Research and Development, U.S. Environmental Protection Agency, Cincinnati, Ohio, USA.
2 Ecosystems Research Division, Office of Research and Development, U.S. Environmental Protection Agency, Athens, Georgia, USA.
3 South Carolina Water Science Center, U.S. Geological Survey, Columbia, South Carolina, USA.
4 South Carolina Water Science Center, U.S. Geological Survey, Clemson, South Carolina, USA.
5 Minnesota Water Science Center, U.S. Geological Survey, Mounds View, Minnesota, USA.

This paper is not subject to U.S. copyright. Published in 2012 by the American Geophysical Union. JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 117, G01006, doi:10.1029/2011JG001806 , 2012



### G01006 1 of 17

### processes for a variety of physiographic regions, spatial scales, and land cover types. Research in this area requires an approach that can both spatially and temporally link the information gained from such studies focusing on different components of the watershed Hg cycle and predict future changes in these processes.

[ 4 

### ] Watershed models are primary tools for ecological riskmanagement; however, watershed Hg modeling is a relatively recent research focus area. There are currently few spatially explicit watershed-scale models that focus on a wide range of Hg cycling processes from the landscape to surface water habitats, particularly for stream systems in thegeographically extensive Coastal Plain region. While valuable and abundant, most available literature on Hg cycling in watersheds explore empirical relationships between Hg and various chemical, hydrologic, and geomorphic correlates. Spatially explicit numerical watershed models, however, can(1) improve spatial and temporal linkages between environmental processes and subsequent water quality whenobservational studies are limited and (2) provide a mathematical approach for predicting the impact of current andfuture changes in climate, land cover, or resource management on watershed Hg cycling [National Research Council

, 2007]. Likewise, empirical watershed models directly link important landscape variables to watershed fluxes using data to provide a strong modeling foundation.

[ 5 

### ] Applying a variety of models, which emphasize different watershed Hg processes, allows simultaneous assessment of the reasonableness of each model ’ s underlying conceptual framework and the quality of the respective mathematical simulations. Thus, employing multipleprocess-based and empirical approaches concurrently promotes understanding of overall watershed Hg dynamics above that achievable by individual modeling efforts or monitoring alone. This is especially critical for Coastal Plain systems, where elevated fish Hg burdens are common [ Bradley et al. , 2010; Brumbaugh et al. , 2001; Glover et al.

, 2010; Scudder et al. , 2009] but where limited understanding of watershed Hg cycling exists.

[ 6 

### ] The purpose of this study was to explore Hg cycling within a small Coastal Plain watershed (McTier Creek) by simulating daily and seasonal total Hg (Hg T

### ) concentrationsand fluxes using multiple watershed-scale models with distinct mathematical frameworks that emphasize different system dynamics. Modeling Hg T 

### processes is an important first step toward advancing the science of watershed MeHg dynamics. As such, we aim to explore the value in using amultiple model and data approach for understanding watershed Hg cycling. Our assessments are based on the range of available data; thus, we do not attempt to extend model simulations beyond these boundaries. Further, identifying the limitations of each model for Hg T 

### simulations (i.e., knowing which Hg dynamics are not represented by the simulations) was considered a fundamental step towardimproved understanding Hg cycling in Coastal Plain watersheds. A secondary goal was to identify current needs in the science of watershed-scale mercury modeling.

[ 7 

### ] We applied two spatially explicit mechanistic models of watershed-scale Hg cycling, an empirical daily Hg load model, and a regression-based water quality flux estimator model for seasonal variations in Hg load estimates, to better understand Hg processes in this Coastal Plain system. The Grid-Based Mercury Model (GBMM) is a spatially explicitmechanistic that emphasizes surface runoff, particulate production, and sediment delivery as primary drivers of Hg T 

### fluxes. Therefore, GBMM simulations provided insight onthe particulate fraction of Hg, suggesting that shallow subsurface flow and subsequent discharge are potentially important controls on the mobilization of particulate Hg T 

### during periods of high upland-to-stream connectivity. The Visualizing Ecosystems for Land Management Assessment(VELMA) model for Hg (VELMA-Hg) is the second spatially explicit mechanistic we employed. VELMA-Hg associates Hg T 

### fluxes with multisoil layer hydrologic transport along with carbon, nitrogen, and mercury cycling. As a result of the carbon cycling module housed in its structure, the VELMA-Hg model afforded insights on DOC-bound Hg cycling in the study watershed. For example, VELMA-Hg simulations revealed that DOC-complexed Hg T 

### becomes more important during base flow conditions compared to high-flow events. We also applied a recently developed empirical Hg load model (TOPLOAD), which is linked tothe TOPMODEL hydrologic framework, in order to simulate daily and seasonal Hg T 

### loadings from the watershed and identify flow components contributing to Hg in the stream. TOPLOAD simulations suggested that saturated subsurface flow is the primary component driving daily HgT loadings, but shallow subsurface flows are important contributors to HgT loads during high-flow events. We used observed data to further understand watershed-scale daily Hg T 

### dynamics and to gauge the accuracy of the models ’ daily Hg T 

### flux estimates. This approach allowed a comparison of statistical relationships and conceptualizations of Hg T 

### dynamics among existing empirical and mechanistic models and data, which are used as tools to characterize Hg cycling in a mesoscale, Coastal Plain watershed. Finally, we compared seasonal load estimates at the watershed outlet for each ofthe three models to outputs from S-LOADEST, a regressionbased water-quality constituent flux estimator, which suggested that dissolved fraction of Hg T 

### is dominant compared to the particulate fraction of Hg T 

### , both seasonally and daily.



### 2. Methods
2.1. Site Description

[ 8 

### ] McTier Creek watershed falls within the Sand Hills Region of the Upper Coastal Plain physiographic province in South Carolina (Figure 1) and is described in detail in Feaster et al. [2010]. The study site covers the approximately 79 km2 

### drainage area above U.S. Geological Survey stream gage 02172305 (McTier Creek near New Holland, South Carolina), and modeling outputs were assessed at the gage station. The climate of the McTier Creek watershed issubtropical humid, and the average annual mean air temperature is approximately 18°C (http://www.dnr.sc.gov/climate/ sco/ClimateData/countyData/county_aiken.php/), based on climate data for period from 1945 to 2008. Land cover within the study area is approximately 49% forest, 21% grassland/ herbaceous, 16% agriculture, 8% wetland, 5% developed, and about 1% open water, using the 2001 National Land Cover Database [Homer et al. , 2004].

[ 9 

### ] The Sand Hills is a topographically distinctive region of the Upper Coastal Plain and is characterized by features ranging from broad, nearly flat ridges to rolling hills (up to GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

2 of 17

### 600 ft above sea level) [ Marshall, 1993]. The two majorgeologic units in the Sand Hills, the Cretaceous age Middendorf Formation and the Eocene age Pinehurst Formation are derived from depositional activity from fluvial systems dating 90 mya, and Pinehurst sediments are marine in origin dating approximately 45 mya, respectively [Bennett and Patton , 2008].

[ 10 

### ] The hydrology of McTier Creek watershed is described in detail elsewhere [ Bradley et al. , 2010; Feaster et al. , 2010]. It is important to note, however, that the coarse-grained sandy sediments of the Coastal Plain exhibit efficient vertical recharge and low surface runoff [ Atkins et al. , 1996; Campbell and Coess , 2010]. Discharge from the shallow flow system represents 72

– 100% of the total groundwater discharge to South Carolina Coastal Plain streams [ Atkins et al. , 1996; Campbell and Coess , 2010]. The gradient and the direction of shallow groundwater flow are toward the stream channel during normal to low-flow conditions, with active groundwater/surface water exchange limited to wetlands and channel margins [Bradley et al.

,2010]. The gradient and the direction of shallow groundwater flow remain toward the stream channel during flood conditions, but with increased area for groundwater/surface water exchange [Bradley et al. , 2010]. Figure 1. Study site: McTier Creek Watershed, South Carolina, USA. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

3 of 17

### 2.2. Data Compilation [ 11 

### ] Observed streamflow data for the study were from the U.S. Geological Survey gage 02172305, McTier Creek near New Holland, South Carolina, starting from the inception of the gage (13 June 2007) through 30 September 2009. These data are used as a baseline for runoff estimates from each model. Observed Hg T 

### (filtered and particulate) (analysis based on DeWild et al. [2002]; Olson et al. [1997]; Olund et al. [2004]), dissolved organic carbon (DOC) (analysis based on Aiken et al. [1992]), and suspended sediment concentrations (SSC) (analysis based on Shreve and Downs [2005]) were sampled 41 times at the stream gage during the June 2007 to September 2009 study period using collection and analyses procedure as described in Lewis and Brigham [2004].

[ 12 

### ] Of the samples collected at the study site, 20% were quality assurance quality control (QAQC) samples. The QAQC samples included seven field blanks for Hg and DOC species related to field collection procedures. Filtered HgT 

### concentrations in field blanks ranged from 0.07 to 0.29 ng L

 1 

### and averaged 0.13 ng L

 1 

### (  0.07 standard deviation). These levels were considered negligible because they represented less than 5% of the environmental sample concentrations, on average. Particulate HgT 

### concentrations in field blanks consistently were below the laboratory reporting level (0.07 ng L  1 

### ). DOC concentrations in field blanks ranged from 0.2 to 0.7 mg L  1 

### and averaged 0.4 mg L

 1 

### (  0.2 standard deviation). These levels were considered negligible because theyrepresented less than 7% of the environmental sample concentrations, on average. Two replicates had relative percentage differences (RPDs) of 10% and 11% for HgT 

### . 2.3. Watershed Hg Modeling

[ 13 

### ] The analysis period for each model ’ s simulation was from 13 June 2007 through 30 September 2009, based on the beginning of the period of record for McTier Creek at New Holland, South Carolina. Table 1 provides a descriptive list of each model applied in this study. We used a 30 m grid cell resolution for the two spatially explicit models, GBMMand VELMA-Hg, which totals over 86,000 grid cells contributing to the watershed-scale analysis. Each model was calibrated using the available observed flow and collected Hg T 

### data during the period of analysis. Final calibrated output is discussed in section 3. 2.3.1. Grid-Based Mercury Model (GBMM)

[ 14 

### ] The Grid-Based Mercury Model was designed tosimulate the daily fluxes and mass balances of water, sediment, and Hg from each GIS raster grid cell to the watershedoutlet using a spatially explicit, process-based model structure [ Dai et al. , 2005]. Spatial input layers to the modelinclude a digital elevation map (DEM), soil types and distributions, and land cover data. Daily fluxes of water, sediment, and Hg from each grid cell are routed through the watershed tributary networks to assessment points along stream channels.

[ 15 

### ] GBMM performs a simple water balance and estimates a modified National Resource Conservation Service (NRCS) curve number (NRCS-CN) for each grid cell on a daily basis. The mass balance of Hg at the watershed outlet is simulated using the equation:

dC s dt 1⁄4 L
V

s  K r þ K l þ K ro þ K e ð Þ * C s ð 1 Þ

L 1⁄4 L p for pervious surfaces ð 2 Þ

L 1⁄4 L f þ L d for forested areas ð 3 Þ

V

s 1⁄4 A c * zd ; ð 4 Þ



### where C s 

### is the concentration of Hg in watershed soils ( mg m  3 

### ), L is the Hg load (m g d

 1 

### ), L p 

### is the Hg atmospheric deposition load on pervious land ( mg d

 1 

### ), L f 

### is the Hg atmospheric deposition load on forest land ( mg d

 1 

### ), L d 

### is the litter decomposition Hg load on forestland ( mg/d), K r 

### is the reduction rate constant (d

 1 

### ) where reduced Hg is assumed to immediately volatilize and is considered a loss from the watershed, K l 

### is the leaching loss constant (d

 1 

### ), K

ro 

### is the runoff loss constant (d

 1 

### ), K e 

### is the erosion loss constant (d

 1 

### ), V s 

### is the watershed soil volume (m 3 

### ), A c 

### is the grid area (m2 

### ), and zd 

### is the watershed soil mixing depth (m). Table 1. Summary of Models Applied in McTier Creek Watershed During the Simulation Period, 13 June 2007

–30 September 2009

Model Model Type Key Drivers of Hg Transport Primary Advantages Primary Scientific Limitations

TOPLOAD-H Semi-distributed load model Flow components of TOPMODEL and mass balance daily Hg accounting Simple, empirical framework; links to TOPMODEL hydrology to identify flow components contributing to in-stream Hg No simulations of watershed Hg cycling and transformations

VELMA-Hg Spatially explicit, process-based Dissolved organic carbon (DOC) within a four-layer soil framework with lateral and vertical movement of water and chemicals Spatially explicit; captures some Hg transformations; multilayer soil framework for hydrology; coupled biogeochemical cycles (Hg, DOC, N) Limited in particulate Hg simulations

GBMM Spatially explicit, process-based Surface soil erosion and in-stream deposition with spatially explicit modified curve number hydrology Spatially explicit; captures some Hg transformations; simulates particulate Hg fate and transformations Limited in dissolved Hg simulations

S-LOADEST Empirical load model Fitted regression curve Simple, empirical framework; models seasonal load variations Limited accurate estimates of daily loads in this watershed; no simulations of Hg watershed cycling and transformations GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

4 of 17

### [ 16 

### ] The current version of GBMM is structured to partition Hg T 

### between its solid and aqueous forms in the soil system. Thus, by calibrating the partition coefficient (Kd

### ) in the model (i.e., the ratio of the equilibrium concentration inthe soil particle to the concentration in soil water), GBMMsimulated Hg T 

### dynamics are largely associated with particulate production and movement from soils to surface waters.

[ 17 

### ] Details of GBMM ’ s equations and sensitivity analysis are previously published [ Golden and Knightes , 2011]. GBMM has been applied and validated in multiple settings in the southeastern United States, including the Piedmont Physiographic Province [ Golden and Knightes , 2011; Golden et al. , 2010] and Coastal Plain Province [ Dai et al.

, 2005; Feaster et al. , 2010]. Final input parameters for McTier Creek watershed are found elsewhere for the hydrology module [Feaster et al. , 2010]. Parameters for the sediment and Hg modules are detailed in the auxiliary material of this paper. 1 

### 2.3.2. The Visualizing Ecosystems for Land Management Assessment Model for Hg (VELMA-Hg)

[ 18 

### ] VELMA is a spatially distributed ecohydrologicalmodel that simulates soil water infiltration and redistribution, evapotranspiration, surface and subsurface runoff, carbon (C) and nitrogen (N) cycling in plants and soils, and the transport of dissolved organic carbon (DOC), dissolved inorganic nitrogen (DIN), and dissolved organic nitrogen (DON) from the terrestrial landscape to streams [ Abdelnour et al. , 2011]. Recently, a Hg cycling submodule was addedto the existing VELMA framework, resulting in a new version of VELMA: VELMA-Hg. Following previous studies that associate Hg T 

### with DOC movement [e.g., Brigham et al. , 2009], particularly in forested settings, we here closely associate VELMA-Hg mercury dynamics in the model with the fate and transport of DOC in the watershed.

[ 19 

### ] VELMA-Hg uses a distributed soil column framework to simulate the lateral and vertical movement of water, heat, and nutrients within the soil (see Figure S1 in the auxiliary material). The modeling domain of VELMA-Hg isa three-dimensional matrix covering the topographical surface (x

– y) and four soil layers (z). The soil column model consists of three coupled submodels: (1) a hydrological model that simulates vertical and lateral movement of water within the soil, including a variable source area (VSA) flow component similar to TOPMODEL [ Wolock, 1993] andlosses of water from the soil and vegetation to the atmosphere, (2) a soil temperature model that simulates daily ground soil layer temperatures from surface air temperature (SAT) and snow depth and, (3) a biogeochemistry model that simulates C, N, and Hg dynamics. The soil column model isthen placed within a catchment framework to create a spatially distributed model applicable to watersheds and landscapes. Adjacent soil columns interact with each other through the downslope lateral transport of water. Surface and subsurface lateral flow are routed using a multiple flowdirection method. A DEM is used to determine flow direction and compute flow contribution area. Required input data includes: air temperature, precipitation, soil texture, soil depth, and DEM. We used air temperature and precipitation values from Feaster et al. [2010]. [ 20 

### ] The net equations governing changes in Hg concentration are solved for each soil column (see below) within the catchment using a forward Euler finite difference approximation. The change in soil mercury areal density, C

Hg

II 

### or C MeHg 

### (g m

 2 

### ) for a given time step is:

dC Hg

II dt 1⁄4 L T
; Hg II þ ks dm þ J in  ks r  ks m  J out ð 5 Þ

dC MeHg dt 1⁄4 L T
;MeHg þ ks m þ J in  ks dm  J out ; ð 6 Þ



### where total L T

, Hg II 

### is the Hg(II) load (g m

 2 

### d

 1 

### ], ks dm 

### is the demethylation source (equation (5)) or sink (equation (6)) (g m

 2 

### d

 1 

### ), J in 

### is the daily soil Hg(II) (equation (5)) or MeHg (equation (6)) influx (g m

 2 

### d

 1 

### ), ks
r 

### is the reduction sink (g m

 2 

### d

 1 

### ), ks
m 

### is the methylation sink (equation (5)) or source (equation (6)) (g m

 2 

### d

 1 

### ), J out 

### is the daily soil Hg(II) (equation (5)) or MeHg (equation (6)) outflux (g m

 2


### d

 1 

### ), and L T

,MeHg 

### is MeHg deposition (g m

 2 

### d

 1 

### ). Further details of the VELMA-Hg model, including model equations and a list of other final calibrated input parameters and initial concentrations, are given in the auxiliary material. 2.3.3. TOPLOAD

[ 21 

### ] TOPMODEL is a physically based, semi-distributed watershed model that simulates the hydrologic fluxes in a watershed [Beven and Kirkby , 1979]. Application of TOPMODEL to the McTier Creek watershed is discussed elsewhere [ Feaster et al. , 2010]. TOPMODEL does notinclude a mass balance algorithm for evaluating waterquality loads. Therefore, we developed a mass balance algorithm and applied it to the surface and subsurface flow components simulated by TOPMODEL, resulting in a loadmodel that assesses the daily mass fluxes for a given constituent (Hg T

### ). Details of the TOPLOAD model are also available elsewhere [Benedict et al. , 2011]. In this investigation, we apply the TOPLOAD-H variation of our TOPMODEL load model, which utilizes simulated surface and subsurface flow components from TOPMODEL and assigns each flow component a Hg T 

### concentration. TOPLOAD-Halso includes a groundwater partitioning algorithm, specifically, upper subsurface (unsaturated) flow and lower subsurface (saturated) flow [Hornberger et al. , 1994]. Thus, daily Hg T 

### flux estimates from TOPLOAD-H provide a simple method for users to quickly evaluate the relative influence of different flow contributions of surface water Hg T 

### on the basis of TOPMODEL-simulated daily flow. The general equation for both TOPLOAD-H is as follows:

LOAD 1⁄4 Q sub 1 C sub 1 þ Q sub 2 C sub 2
... þ Q sub n C sub n þ Qqinf C qin

f
þ Qqimp C qimp þ Qqsrip C qsrip þ Qqof C qof ; ð 7 Þ



### where, LOAD represents the estimated watershed load of a given constituent; Q sub1 

### , Q sub2 

### , ... , and Q subn 

### , represent the TOPMODEL subsurface flow as distributed across the number of soil zones, n ; Qqinf 

### , Qqimp 

### , Qqsrip 

### , and Qqof 

### represent the flow associated with the respective TOPMODEL surface flow components; and C sub1 

### , C sub2 

### , ... , C subn 

### , C qinf 

### , C

qimp 

### , C qsrip 

### , and C qof 

### represent the constituent concentration associated with the respective flow component. TOPLOAD-H is essentially a mixing model and the fluxes 1 Auxiliary materials are available in the HTML. doi:10.1029/ 2011JG001806. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

5 of 17

### can be evaluated by assigning different concentrations to the various flow components. Here, we assign a daily 5 ng L

 1 

### Hg T 

### concentration to each TOPMODEL flow component for TOPLOAD-H load calculations. This value represents the average daily concentration for all flows (i.e., low, high, and base flow) on the basis of field measurements from the 41 sampling days during the simulation period and expert judgement. While each basin will have unique considerations in the assignment of constituent concentrations associated with the various surface and subsurface flow components, the studies by Hornberger et al. [1994] and Robson et al. [1992] provide practical insights on how this assignment can be accomplished. 2.3.4. LOAD ESTimator (S-LOADEST) Model

[ 22 

### ] We used a regression-based modeling program, S-LOADEST, to model seasonal estimates of Hg T 

### loads. The strength of the S-LOADEST estimated Hg T 

### flux (particulate + filtered) regression model for the period of record is relatively strong ( R 2 

### = 0.89, p < 0.0001), based on the 41 day sampling period for our study analysis. S-LOADEST is a U.S. Geological Survey “plug-in” version of LOADEST (LOAD ESTimator) (David Lorenz, U.S. Geological Survey, personal communication, 2010) for S+® software (TIBCO Spotfire Co., Palo Alto, Calif.). The original LOADEST is aFORTRAN program that has been used extensively for estimating constituent loads in streams and rivers [Runkel et al.

, 2004]. The rating curve approach takes into account flow and seasonal dependence of constituent concentrations, which is not always adequately accounted for by other approaches (e.g., linear interpolation between sampling dates).

[ 23 

### ] Instantaneous constituent loads are computed by the following equation:

LHg 1⁄4 C Hg * Q i * C l :; ð 8 Þ



### where LHg 

### is the Hg species (or other constituent of interest) load at the time of sampling in milligrams per day, C Hg 

### is the concentration of the Hg species (or other constituent of interest) in ng/L; Q i 

### is the instantaneous stream discharge at the time of sampling in cubic feet per second; and C l 

### is a unit conversion factor (2.447). The model selected for this study accounted for streamflow and seasonality [ Cohn et al. , 1992; Helsel and Hirsch , 1992] with the daily load as:

L 1⁄4 b
o þ b 1 LnQ þ b2 LnQ2 þ b 3 sin 2 p T ð Þþ b4 cos 2 p T ð Þð 9 Þ



### where L is the natural logarithm ( Ln ) of the estimated load in milligrams per day; Q , log of the mean daily stream discharge in cubic feet per second; T, centered time in decimal years; sin, sine; cos , cosine; p , pi; b
n 

### , estimated coefficients for each variable. For the constituent of interest, the formulated regression model was used to estimate loads over a selected time interval (estimation period) of October 2004 to September 2009. Mean load estimates, standard errors, and 95% confidence intervals were developed on a seasonal and annual basis. Validation of model output included examination of model residuals to ensure linear fit, uniform scatter around the fit, normality of distribution, and linearity with all explanatory variables [Helsel and Hirsch , 1992]. This rating curve approach has recently been applied in multiple surface water systems [Brigham et al. , 2009; Journey et al. , 2011].



### 3. Results

3.1. Runoff

[ 24 

### ] Direct comparisons between modeled GBMM and TOPMODEL runoff, including calibration and validationresults at a nested stream gage within McTier Creek watershed, have been presented elsewhere [ Feaster et al. , 2010].The addition of VELMA-Hg generated a new set of comparison statistics (Table 2). The mean and standard deviation of simulated daily flows for the study period were relatively similar among observed data and the three models. Some slight variations were evident. VELMA-Hg ’ s multilayered hydrologic model represented a minor improvement over GBMM ’ s flow simulations; however, the Nash-Sutcliffe for VELMA hydrology is substantially lower than that of TOPMODEL and GBMM (Table 2), which potentially reflects VELMA-Hg ’ s over predictions during high-rainfall
events.

[ 25 

### ] Simulations suggest that daily flow from each of the three models is fairly well matched to the distribution inobserved flow conditions (Figure 2). Differences in the performances of the hydrologic models can be seen in theextremes of the flow distribution curve. VELMA underestimated flows except for high flows greater than 1.0 m3 

### s  1



### , which are overestimated. GBMM accurately simulated the occurrence of average flows 10%

–80% of the time butunderestimated both high and low flows. TOPMODEL accurately simulated flows 5%

–80% of the time but overestimated low flows. VELMA simulations regularly peaked above observed flow during high-flow periods (Figure 3). However, at or below base flow (approximately 0.4 m3 

### s  1 

### ), VELMA daily flows were well below observed flow and exhibited rapid Table 2. Summary Statistics for Daily Hydrologic Simulations in McTier Creek Watershed During the Simulation Period, 13 June 2007

–30 September 2009 a

Simulated Runoff (m 3 s  1 )

Observed GBMM TOPMODEL VELMA
Mean (m 3 s  1 ) 0.61 0.62 0.61 0.59 Maximum (m 3 s  1 ) 4.62 6.57 3.68 6.55 Minimum (m 3 s  1 ) 0.07 0.04 0.14 0.08 Range (m 3 s  1 ) 4.54 6.54 3.54 6.47 Standard Deviation (m 3 s  1 ) 0.55 0.51 0.49 0.76 RMSE

— 0.35 0.29 0.16 R2

— 0.54 0.64 0.57 NS

— 0.50 0.64 0.20

aRMSE is the root mean squared error, and NS is the Nash-Sutcliffe efficiency index [Nash and Sutcliffe , 1970]. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

6 of 17

### groundwater recession (Figure 3). Groundwater recession for both TOPMODEL and GBMM was often slower than observed stream recession (Figure 3).
3.2. Hg T 

### 3.2.1. Daily Hg T 

### Concentrations

[ 26 

### ] Average observed filtered Hg (Hg F,OBS 

### ) concentrations comprised over two-thirds of the average Hg T,OBS

### (particulate + filtered) concentrations (Table 3). An exception to this trend involved some peak flow events, such as those in August 2008 and October 2008 (Figure 3), when observed particulate Hg (Hg P,OBS 

### ) concentrations exceeded those of Hg F,OBS 

### . All observed Hg T 

### concentrations (Hg F,

OBS 

### , Hg P,OBS 

### , and Hg T,OBS 

### ) increased rapidly during January and February 2009, a period at or near base flow conditions preceded by three strong peaks in the hydrograph in December 2008 (Figure 3). GBMM Hg T 

### (Hg T,GBMM

### ) and VELMA Hg T 

### (Hg T,VELMA

### ) daily concentrations for the sampling days (n = 41) were approximately 93% of the Hg T,

OBS 

### concentrations. However, this estimate varies widely ranging from 13% to 367% of observed Hg T 

### concentrations on any given day. Among observed data (n = 41), Hg T,OBS

### concentrations exhibited statistically significant relationships with streamflow, observed dissolved organic carbon (DOC OBS 

### ) concentrations, and observed total suspended sediment concentrations (TSS OBS 

### ; Table 4). Hg T,OBS 

### concentrations were also correlated with HgF ,OBS 

### (r = 0.88) and Hg P,OBS 

### (r = 0.62) concentrations (both p < 0.0001 ; see Figures S2

– S5 in the auxiliary material for scatterplots of all statistical relationships). Hg P,OBS 

### concentrations were significantly related to DOC OBS 

### and TSS OBS 

### concentrations (Table 4), and Hg F,OBS 

### concentrations were significantly related to flow (Table 4).

[ 27 

### ] Average Hg T,GBMM 

### concentrations and Hg P,OBS 

### concentrations were similar on the days sample data were collected (1.18 ng L

 1 

### and 1.27 ng L

 1 

### , respectively, n = 41; Table 3); however, the standard deviation of Hg T,GBMM 

### concentrations was larger than that of Hg P,OBS 

### . Hg T,GBMM 

### concentrations were closest to total (filtered (Hg F,OBS 

### ) + Hg P,OBS 

### ) observed daily Hg T 

### (Hg T,OBS 

### ) concentrationswhen sampling was conducted during peaks in the hydrograph (e.g, August, October, and December 2008 and May and July 2009; Figure 3). Daily Hg T,GBMM 

### concentrations were significantly correlated with simulated sediment (TSS GBMM

### ), TSS OBS 

### , and Hg P,OBS 

### concentrations in addition to streamflow (Table 4; see also Figures S2

– S5 in the auxiliary material).

[ 28 

### ] Average Hg T 

### , VELMA 

### and Hg T,OBS 

### concentrations were similar (3.54 ng L

 1 

### and 3.92 ng L

 1 

### , respectively) and much higher than Hg T 

### , GBMM 

### and Hg P,OBS 

### concentrations (Table 3 and Figure 3). Surprisingly, Hg T 

### , VELMA 

### concentrations often matched Hg P,OBS 

### concentrations rather than Hg F,OBS 

### (Figure 3), which is supported by the statistically significant relationship between DOC OBS 

### and Hg P 

### , OBS 

### (Table 4) and leads to multiple interpretations elucidated in section 4.2 of the Discussion. Hg T 

### , VELMA 

### concentrations were significantly correlated with DOC OBS 

### concentrations, Hg P,OBS 

### concentrations, and simulated streamflow (Table 4). Hg T 

### , VELMA 

### concentrations were inversely related to VELMA-Hg DOC (DOC VELMA

### ) concentrations (Table 4); Figure 2. Flow duration curves of simulated (GBMM, VELMA-Hg, and TOPMODEL) and observed flow conditions during the simulation period (13 June 2007

–30 September 2009). GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

7 of 17

### however, this relationship was strongest when streamflow was above base flow conditions (r =

 0.87, p < 0. 0001). In fact, at or below base flow conditions the relationship between DOC VELMA 

### and Hg T,VELMA 

### concentrations is relatively weak and positive (r = 0.25, p < 0. 0001). Clusters of simulated DOC and Hg T 

### pairs showing this change in relationship further suggest a change in the DOC-Hg T 

### dynamics in the model under different flow regimes (see scatterplots in Figures S2

– S5 in the auxiliary material). The relationship of simulated Hg T 

### to streamflow resulted in overall increases in Figure 3. Simulated daily runoff and Hg T 

### concentrations. OBS, observed streamflow; Hg T,OBS 

### , total (filtered + particulate) observed mercury; Hg F,OBS 

### , observed filtered mercury; and Hg P,OBS 

### , observed particulate mercury.

Table 3. Mean and Standard Deviation of Observed and Simulated Runoff, Total Mercury (Hg T) Concentrations, and Hg T Fluxes for the Simulation Period, 13 June 2007

–30 September 2009

Runoff GBMM (m 3 s  1 ) Runoff VELMA-Hg (m 3 s  1 ) Runoff TOPMODEL (m 3 s  1 ) Runoff OBS (m 3 s  1 )

Average a 0.56 0.68 0.70 0.72 Std Dev a 0.40 0.55 0.51 0.79

Hg T,GBMM (ng L

 1 ) Hg T,VELMA (ng L

 1 ) Hg T,TOPLOAD (ng L

 1 )b Hg P,OBS (ng L

 1 ) Hg F,OBS (ng L

 1 ) Hg T,OBS (ng L

 1 ) Average a 1.18 3.54 5.00 1.27 2.64 3.92 Std Dev a 2.01 2.75 0.00 0.44 0.93 1.17

Hg T,GBMM (mg km

 2 d

 1 ) Hg T,VELMA (m g km

 2 d

 1 ) Hg T,TOPLOAD (m g km

 2 d

 1 ) b Hg P,OBS (m g km

 2 d

 1 ) HgF,OBS ( mg km

 2 d

 1 ) Hg T,OBS (mg km

 2 d

 1 )

Average a 434 4440 3740 1200 2520 3660 Std Dev a 862 7310 2700 1840 3940 5310

a Statistics calculated from days when observed data are available, n = 40 and 41.
b TOPLOAD model is assigned a constant daily concentration for each of TOPMODEL ’ S flow components. Hg P,OBS , observed particulate Hg; Hg F,OBS , observed filtered Hg; Hg T,OBS , Hg P,OBS + Hg F,OBS . GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

8 of 17

### HgT,VELMA 

### concentrations with increased discharge (Figure 3). However, at streamflow greater than 1.5 m3 

### s  1 

### , minimal increases occurred.

[ 29 

### ] Concentration values for the TOPLOAD-H model did not vary across the simulation period; a single daily Hg T 

### concentration was assigned for each TOPMODEL flow component (5 ng l

 1 

### ). 3.2.2. Daily Hg T 

### Fluxes

[ 30 

### ] Hg T,OBS 

### fluxes followed trends in the hydrograph relatively closely (Figure 4); however, increases in fluxes in January and February 2009 correspond with variations in HgT 

### concentrations rather than flow conditions. HgF,OBS 

### fluxes are higher than HgP,OBS 

### fluxes during every measurement period except January and February 2009. Among observed data (n = 41), Hg T,OBS 

### , Hg F,OBS 

### , and Hg P,OBS 

### fluxes are each significantly correlated with streamflow, DOC OBS 

### , and TSS OBS 

### (Table 5). Hg T,OBS 

### fluxes are also significantly correlated with Hg F,OBS 

### (r = 0.96) and Hg P,OBS 

### (r = 0.81), both
p < 0.0001. Scatterplots of the statistical relationships among the variables are shown in Figures S2

– S5 of the auxiliary material.

[ 31 

### ] Average Hg T,GBMM 

### flux simulations (434 mg km

 2


### d

 1 

### ) were low compared to average Hg T 

### flux estimates from the other models (Table 3), which resulted from GBMM underestimates of low flow conditions (Figure 2). Although average daily Hg T,GBMM 

### fluxes are typically lower than observed Hg T 

### fluxes, Hg T,GBMM 

### fluxes peak above observed fluxes in December 2009 and March 2009 (Figure 4) and follow the trend of simulated flow conditions during that period. Hg T,GBMM 

### fluxes exhibited a strong positive statistically significant relationship with Hg P,OBS 

### fluxes

, 

### TSS GBMM 

### fluxes, and streamflow and were correlated to a lesser extent with Hg T,OBS 

### and Hg F,OBS 

### fluxes (Figure 4).

[ 32 

### ] Average daily Hg T,VELMA 

### fluxes (4438 mg km

 2 

### d

 1 

### ) were higher than all other modeled estimates (Table 3). VELMA-Hg fluxes peaked much higher and more regularly than Hg T,GBMM 

### and TOPLOAD Hg T 

### (Hg T,TOPLOAD 

### ) estimates and typically surpassed Hg T,OBS 

### as well (Figure 4). Spikes in the Hg T,VELMA 

### fluxes typically corresponded with peaks in the hydrograph (Figure 4). Hg T,VELMA 

### fluxes were significantly related to streamflow, DOC OBS 

### , and Hg T,OBS 

### , Hg P,OBS 

### , and to a lesser extend Hg F,OBS 

### (Table 5).

[ 33 

### ] The TOPLOAD-H model produced Hg T 

### fluxes higher than HgT,GBMM 

### for the simulation period but lower than HgT,

VELMA. 

### A salient feature and the goal in this TOPLOAD-H study application, however, was its use as a tool to assist in determining which hydrologic components contribute to daily HgT 

### fluxes. The Hg T,TOPLOAD 

### simulations suggested that total daily Hg T 

### fluxes are largely driven by saturated subsurfaceflow conditions (TLOAD-lower, Figure 5). However, infiltration excess overland flow (TLOAD-infiltration, Figure 5), precipitation on water surfaces (TLOAD-riparian, Figure 5), TOPMODEL overland flow (flow just below the surface soil), and, in a more limited capacity, flow from impervious surfaces, regulated peaks in daily Hg T 

### fluxes. TOPLOAD-H simulated these peaks in Hg T 

### fluxes with greater accuracy when HgT 

### concentrations were increased (from 5 ng L

 1 

### to 15 ng L  1 

### ) for the flow component that represents direct precipitation on water (TLOAD-riparian) (see Figure S6 in the auxiliary material). 3.2.3. Seasonal Hg T 

### Concentrations

[ 34 

### ] Average Hg T,OBS 

### concentrations did not exhibit a strong seasonal trend, although the two winter seasons had the highest average concentrations (Figure 6). LOADEST total and filtered Hg T 

### concentrations were typically higher than LOADEST particulate (Hg P,LOADEST

### ) and Hg T,GBMM 

### concentrations for each season. Total LOADEST Hg T 

### (Hg T,LOADEST

### ) concentrations averaged across seasons were strongly driven by trends in Hg F,OBS 

### (Figure 6) and often matched the average Hg T,OBS 

### as well, the latter which reflectsthe structure of the load model (i.e., using observed concentration and streamflow to estimate fluxes). As a result, LOADEST produced the highest Hg T 

### concentrations in the winter, similar to observed concentrations. VELMA-Hg simulations were variable and showed limited similarity to the trends of the other modeled or observed data. GBMM Table 4. Correlation Coefficients and Significance Levels Between Observed Total Hg (Hg T,OBS ) Concentrations, Observed Particulate Hg (HgP,OBS ) Concentrations, Observed Filtered Hg Concentrations (Hg F,OBS ), GBMM-Simulated Hg T Concentrations (Hg T,GBMM), VELMA-Hg Hg T Concentrations (Hg T,VELMA), and Other Observed and Simulated Variables a,b

Hg T,OBS HgF,OBS Hg P,OBS Flow OBS DOC OBS TSS OBS

Flow OBS 0.35* 0.31* 1.00*** 0.67*** DOC OBS 0.33* 0.61** 0.67*** 1.00*** 0.54** TSS OBS 0.35* 0.84*** 1.00***

Hg T,OBS HgF,OBS Hg P,OBS Hg T,GBMM TSS OBS TSS GBMM Flow OBS Flow GBMM Hg T,GBMM 0.65*** 1.00*** 0.60*** 0.87*** 0.30* 0.36*** Flow GBMM 0.35* 0.31* 0.23*** 0.74*** 1.00*** TSS GBMM 0.35* 0.84*** 0.87*** 0.48** 1.00*** 0.56*** 0.23***

Hg T,OBS HgF,OBS Hg P,OBS Hg T,VELMA DOC OBS DOC VELMA Flow OBS FlowVELMA Hg T,vELMA 0.58*** 1.00*** 0.62***

 0.55** 0.31* 0.64*** FlowVELMA 0.44** 0.64*** 0.64*** 0.88*** 0.60*** 0.64*** 1.00*** DOC VELMA  0.46**

 0.32*

 0.49**

 0.55**

 0.54** 1.00***

 0.39* 0.60***

aDOC OBS , observed dissolved organic carbon concentrations; TSS OBS , observed total suspended solids concentration; FLOWGBMM , GBMM-simulated streamflow; TSS GBMM , GBMM-simulated total suspended solids; FLOWVELMA , VELMA-simulated streamflow; DOC VELMA , VELMA-simulated dissolved organic carbon concentrations.

b Significance levels ( p values) are indicated as follows: a single asterisk (*) is p < 0.05, a double asterisk (**) is p < 0.01, and a triple asterisk (***) is
p < 0.0001. All relationships that include observed data, n = 41; all relationships among modeled data only, n = 841. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

9 of 17

### Figure 4. Simulated daily runoff and Hg T 

### fluxes. OBS, observed streamflow; Hg T,OBS 

### , total (filtered + particulate) observed mercury; Hg F,OBS 

### , observed filtered mercury; and Hg P,OBS 

### , observed particulate
mercury.

Table 5. Correlation Coefficients and Significance Levels Between Observed Total Hg (Hg T,OBS ) Fluxes, Observed Particulate Hg (Hg P,OBS ) Fluxes, Observed Filtered Hg Fluxes (HgF,OBS ), GBMM-Simulated Hg T Fluxes (Hg T,GBMM), VELMA-Hg Hg T Fluxes (Hg T,VELMA), and Other Observed and Simulated Variables a,b

Hg T,OBS Hg F,OBS Hg P,OBS Flow OBS DOC OBS TSS OBS

Flow OBS 0.95*** 0.94*** 0.71*** 1.00*** 0.64*** 0.51** DOC OBS 0.61*** 0.58*** 0.51** 0.64*** 1.00*** 0.40* TSS OBS 0.63*** 0.41** 0.92*** 0.51** 0.40* 1.00***

Hg T,OBS Hg F,OBS Hg P,OBS Hg T,GBMM TSS OBS TSS GBMM Flow OBS Flow GBMM Hg T,GBMM 0.62*** 0.41** 0.89*** 1.00*** 0.78*** 0.84*** 0.55*** 0.75*** Flow GBMM 0.63*** 0.62*** 0.46** 0.75*** 0.30* 0.72*** 0.74*** 1.00*** TSS GBMM 0.76*** 0.63*** 0.81*** 0.84*** 0.57*** 1.00*** 0.78* 0.72***

Hg T,OBS Hg F,OBS Hg P,OBS Hg T,VELMA DOC OBS DOC VELMA Flow OBS FlowVELMA Hg T,vELMA 0.59*** 0.38* 0.86*** 1.00*** 0.49** 0.51** 0.97*** FlowVELMA 0.69*** 0.50** 0.88*** 0.97*** 0.71*** 0.17*** 0.47** 1.00*** DOC VELMA 0.58*** 0.70*** 0.38* 1.00*** 0.76*** 0.17***

aDOC OBS , observed dissolved organic carbon fluxes; TSS OBS , observed total suspended solids fluxes; FLOWGBMM , GBMM-simulated streamflow; TSS GBMM , GBMM-simulated total suspended solids fluxes; FLOWVELMA , VELMA-simulated streamflow; DOC VELMA , VELMA-simulated dissolved organic carbon fluxes.

b Significance levels ( p values) are indicated as follows: a single asterisk (*) is p < 0.05, a double asterisk (**) is p < 0.01, and a triple asterisk (***) is
p < 0.0001. All relationships that include observed data, n = 41; all relationships among modeled data only, n = 841. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

10 of 17

### seasonally averaged Hg T 

### concentrations were consistently low, but generally followed the trend of Hg P,LOADEST 

### concentrations from Autumn 2007 through Winter 2009. 3.2.4. Seasonal HgT Fluxes

[ 35 

### ] Trends in observed seasonal Hg T,OBS 

### fluxes were similar to that of seasonally averaged Hg T,OBS 

### concentrations; however, Hg T 

### fluxes were elevated during autumn 2008 when streamflow conditions were highest during the simulation period (Figure 6). Tukey Studentized Range tests determined that winter Hg T,OBS 

### and Hg F,OBS 

### fluxes were significantly higher (p < 0.05) than other seasons. The Hg F,

LOADEST 

### fluxes drive the curve of the Hg T,LOADEST 

### fluxes; however, Hg T,TOPLOAD 

### followed a similar trend, with the exception of Spring 2009 when TOPMODEL overestimates runoff conditions (Figure 4). The Hg T,LOADEST 

### and Hg T,

TOPLOAD 

### estimates exhibited the highest seasonal fluxes among modeled Hg T 

### estimates. Hg T,GBMM 

### and Hg P,LOADEST 

### fluxes followed a similar pattern during each season of our analysis period. The Hg T,VELMA 

### seasonal flux simulations exhibited a similar trend as all other observed and simulated Hg T 

### fluxes from Spring 2007 through Winter 2008. Following this, Hg T,VELMA 

### fluxes matched observed fluxes relatively well in Summer 2008 and Spring 2009 but exhibited minimal seasonal trends during the simulation period.



### 4. Discussion

[ 36 

### ] Applying three models (GBMM, VELMA-Hg, TOPMODEL/TOPLOAD-H) with differing underlying structures and representative processes, in addition toobserved data, affords a greater number of tools for understanding of Hg cycling and loadings beyond the dynamics and assumptions of single model simulations or monitoring alone. The particular emphasis of this discussion focuses on what can be learned about Hg fate and transport and Hgmodeling at the watershed scale in McTier Creek by variations in simulation output from these diverse models anddata. Insight is gained not only when models perform satisfactorily but also when simulations fail to align closely with observed data or produce reasonable results despite a potentially mischaracterized conceptualization of particular system dynamics in the model equations. In this way, although the structure of mechanistic models is generally not modified when new data become available (i.e., unlike Figure 5. Daily TOPLOAD output for the simulation period. Hg T,OBS 

### , Hg F,OBS 

### , and Hg P,OBS 

### represent observed combined (filtered + particulate) Hg T 

### , observed filtered Hg, and observed particulate Hg, respectively. TLOAD-total is the cumulative predicted daily Hg T 

### output. Source contributions include unsaturated shallow subsurface flow (TLOAD-lower); saturated shallow subsurface flow (TLOAD-upper);flow from riparian areas (TLOAD-riparian); total subsurface flow (TLOAD-subsfc); flow from impervious surfaces (TLOAD-impervious); infiltration excess overland flow, above and/or over saturated soils (TLOAD-infiltration); and subsurface flow within fine upper layer of wetted soil (TLOAD-overland). GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

11 of 17

### empirical models), they provide valuable information and approaches to evaluate different conceptualizations of watershed processes. Therefore, this discussion adheres to the following guidelines: (1) observed data and mechanisticmodel results are mutually informative, (2) model parameters are not forced beyond realistic values in an effort to match simulation output to observed values (i.e., somepotential model output error is allowed to maintain reasonable parameterization of the model), (3) the potential contribution of processes that are not explicitly represented in an individual model is always acknowledged.
4.1. Daily Streamflow as the Transport Mechanism for Hg T



### [ 37 

### ] Daily runoff estimates from each model were compared to elucidate Hg T 

### transport characteristics at McTier Creek watershed. This complementary modeling approach highlighted multiple explanations for water transport and interacting dynamics within the watershed ’ s hydrologic system, particularly during peak flow periods. For example, GBMM runoff simulations matched the daily observed hydrograph, particularly peaks in the hydrograph, relativelywell (Figures 2 and 3). Therefore, if GBMM were implemented alone to understand the hydrologic and Hg dynamics of the watershed, simulations would suggest that during high-rainfall events, the system responds rapidly viainfiltration-excess overland flow, the hydrological mechanism driving GBMM. This interpretation would, in turn, affect the analysis of Hg transformation and transport dynamics in the watershed. However, when GBMM resultswere analyzed alongside VELMA and TOPMODEL simulations, the three models exhibited similar trends in the daily hydrograph and matched daily observed runoff relatively well despite their diverse underlying structures (Figure 2; Table 2). Moreover, observed and simulated hydrographs of all three models reflected hydroperiods typical of Coastal Plain streams (Figure 2), exhibiting lower flows from June through October and higher flows from November through March [Hupp , 2000].

[ 38 

### ] Similar results for three models with different conceptualizations of primary rainfall-runoff mechanisms (e.g., infiltration-excess overland flow with GBMM and VSA contributions using VELMA and TOPMODEL) suggest a

“modified” interpretation of peak flow hydrology

— and potential Hg transport

— in McTier Creek watershed. The infiltration-excess overland flow of GBMM and the VSA contributions of VELMA and TOPMODEL (along with analogous flow components in numerous other models), Figure 6. Simulated seasonal Hg T 

### concentrations and flux estimates. Hg T,OBS 

### is the seasonally averaged observed Hg T 

### , and the error bars represent one standard deviation of the mean. Hg T,LOADEST 

### , Hg F,LOADEST 

### , and Hg P,LOADEST 

### represent simulations from S-LOADEST combined (filtered + particulate) Hg T 

### , S-LOADEST filtered Hg, and S-LOADEST particulate Hg, respectively. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

12 of 17

### originated as conceptualizations of water transport in upland systems [Hupp , 2000]. Coastal Plain hydrologic dynamics are typically not well represented in current hydrological models [e.g., Lu et al. , 2006]. However, a rapid runoff response in McTier Creek following rainfall events also may be attributed to shallow groundwater flow path(s) with a kinetic (temporal) response similar to that of higher gradient water transport in upland systems. The flashy character (rapid rising limb of the hydrograph) typically associated with lower permeability, higher gradient upland systems alsocould occur in permeable Coastal Plain sediments if subsurface water storage is restricted by a shallow subsurface impermeable zone (e.g., bedrock). This “modified” conceptualization of the hydrologic framework is consistent with empirical observations of shallow bedrock in the upper reaches of McTier [ Feaster et al. , 2010] along with a highly permeable surficial aquifer, rapid vertical infiltration, limited surface runoff and continuous shallow subsurface discharge of groundwater toward the stream throughout the McTier watershed.
4.2. Daily Hg T 

### : Insights for Characterizing Hg T 

### Processes in Coastal Plain Watersheds

[ 39 

### ] Given that existing watershed Hg T 

### models represent diverse conceptualizations of drivers that influence in-stream Hg concentrations and fluxes, applying multiple models concurrently can improve the understanding of overall watershed Hg dynamics above that achievable by individual modeling efforts or monitoring alone. As such, each modelprovides insights into a particular set of Hg cycling processes, specifically: (1) GBMM elucidated, in part, particulate Hg T 

### fate and transport in the watershed, (2) VELMA-Hgsimulations represented what can be learned about the dissolved fraction of Hg T 

### based on current conceptualizations in watershed Hg T 

### modeling, (3) TOPLOAD model outputs reflected hydrologic sources in the system that contribute to Hg T 

### at the watershed outlet, and (4) LOADEST revealed patterns (or lack thereof) of seasonal variations in Hg T 

### at the watershed outlet. Details of interpretations gained from each model and their relationship of simulations to observed data follow.

[ 40 

### ] Two particular findings from GBMM simulations partially explained patterns in particulate Hg T 

### in the watershed. First, Hg T,GBMM 

### flux simulations exhibited a statistically significant relationship with simulated and observed streamflow, especially for flows above 1.5 m 3 

### s  1 

### (Figure 4). This pattern is consistent with the above “modified” conceptualization of the hydrologic framework in which a highly responsive (flashy), shallow groundwater flow path(s) becomes increasingly important to Hg T 

### transport under highflow conditions. This finding is also in line with a previously published conceptual model of Coastal Plain hydrology depicting increased floodplain to stream connectivity during high-flow events [Bradley et al. , 2010].

[ 41 

### ] Second, GBMM simulations revealed the importance of erosional processes for Hg transport during high-flow events and the need for additional interpretations of in-stream particulate Hg dynamics. For example, in McTier Creek watershed, Hg F 

### , OBS 

### typically comprised 80

– 85% of the Hg T,

OBS 

### concentrations during base flow conditions, as expected in watersheds with a high proportion of forested land [ Munthe et al. , 2007]. During peak flow periods, however, the particulate contribution to total Hg T 

### concentrations often exceeded that of Hg F,OBS 

### . Within these periods, GBMM simulations of daily Hg T 

### fluxes followed the trend of Hg P,

OBS 

### , and both Hg T,GBMM 

### and Hg P 

### , OBS 

### (Figure 4) were highly correlated with TSS OBS 

### and TSS GBMM 

### . However, GBMM fluxes were also three times lower than those of observed data (Table 3). This finding reflects simulated low flows fromGBMM during precipitation events and suggests that mobilization of Hg via erosion of surface sediment and subsequent deposition in surface waters, primary transport mechanisms in GBMM, are important; however, these are not the primary drivers of particulate Hg T 

### dynamics in McTier Creek during high-flow conditions. Other processes, such as stream bank erosion and benthic sediment re-suspension, may potentially lead to increased particulate Hg T 

### in the water column but are not explicitly modeled in the current version of GBMM and are not distinguishable by the current data set. GBMM modeling results, therefore, point to the need for (1) further monitoring to determine the relative importance of in-stream and landscape sources of particulate Hg T 

### during high-flow events and (2) inclusion of stream bank erosion and benthic sediment re-suspension processes in watershed Hg models.

[ 42 

### ] VELMA-Hg relies strongly on the interactions between Hg T 

### and DOC to drive dissolved Hg T 

### dynamics. While GBMM places limited emphasis on DOC as a key factor influencing Hg T 

### concentrations and fluxes, VELMAHg simulations highlight the relevance of DOC in McTier Creek watershed. The link between Hg T 

### and DOC is particularly important here, as previous studies suggest that DOC dynamics are a primary regulator of Hg T 

### concentrations and fluxes in Coastal Plain watersheds [Barringer et al. , 2010; Guentzel, 2009]. As such, the statistical relationship between Hg T,VELMA 

### concentrations and DOC OBS 

### concentrations versus VELMA-simulated DOC (DOC VELMA

### ) concentrationsreveal two different, yet complementary, insights about dissolved Hg T 

### fate and transport in McTier Creek and highlight that the DOC-Hg T 

### relationship is much more complex than a simple linear model could convey.

[ 43 

### ] As Hg T,VELMA 

### increases, DOC VELMA 

### is predicted todecrease. This is especially apparent for simulations of highflow events, as evidenced by (1) the strengthened inverse relationship between DOC VELMA 

### concentrations and Hg T 

### above base flow conditions of approximately 0.4 m 3 

### s  1 

### (i.e., R2 

### changes from 0.30 with combined flow conditions to 0.75 at flows above base flow) and (2) a distinct trend in small clusters of these two modeled constituents at DOC concentrations above approximately 15 mg L

 1 

### (seeFigure S6 in the auxiliary material). Because of the variations in the Hg T 

### -DOC relationship, the correlation between Hg T,VELMA 

### and DOC VELMA 

### fluxes was not significant (Figure 4).

[ 44 

### ] This statistically inverse relationship between DOC VELMA 

### and Hg T 

### , VELMA 

### is imposed by (1) low DOC concentrations and significant Hg T 

### concentrations in rainfall and (2) the model structure, namely the dominance of VSA flow in the model during high-precipitation events. Hg is more concentrated in precipitation than in surface waters. According to VELMA-Hg, during precipitation events the Hg-concentrated rain falls onto the watershed, and VSA runoff in the VELMA-Hg model becomes the dominant flow process. As originally conceived, this flow mechanism would rapidly flush precipitation falling on saturated areas GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

13 of 17

### from the landscape to surface waters, affording minimal infiltration and limited contact with the surface soil system. This rapid runoff process, combined with direct precipitation of rainwater on open water surfaces, would effectively bypass potential interactions with DOC and flushing of DOC-bound Hg in the soil system. As such, Hg T 

### , VELMA 

### concentrations could increase without a corresponding increase in DOC VELMA 

### concentrations (Figure 3). However, Hg T 

### , VELMA 

### , Hg T 

### , OBS 

### , and DOC OBS 

### concentrations increased concomitantly, suggesting a direct relationship between (1) streamflow and Hg T 

### , VELMA 

### concentrations, in addition to streamflow and Hg T 

### , VELMA 

### fluxes (Figure 4), and (2) DOC OBS 

### and Hg T 

### , VELMA 

### concentrations (Figure 3), in addition to DOC OBS 

### and Hg T 

### , VELMA 

### fluxes (Figure 4). The increase in Hg T 

### , VELMA 

### concentrations was particularly strong when runoff increased slightly above base flow (Figure 3).

[ 45 

### ] The differences in the relationship between Hg T 

### , VELMA 

### concentrations and DOC VELMA 

### versus DOC OBS 

### suggest that there are additional Hg dynamics at work in the watershed that are not captured by the VELMA-Hg model. For example, DOC removal from the upper soil horizons often increases when connectivity between the uplands and surface waters is maximized [Boyer et al. , 1997]. The collective increase of streamflow, DOC OBS 

### , Hg T,OBS 

### , and Hg T,

VELMA 

### indicates DOC flushing continues to occur in the system even during high-precipitation events. This disparity between model prediction and observation is consistent with the “modified” conceptualization of the VELMA-Hg VSA flow component in Coastal Plain environments and a highly responsive shallow subsurface flow component, which would continue to supply DOC from the sediment to the stream even during high flow. Increased DOC flushing from highly organic floodplain surface soils near the stream due to vertical groundwater discharge under flood conditions also is consistent with continued DOC supply at high flow [ Bradley et al. , 2010]. The good agreement between VELMA-Hg and observed Hg and DOC data at low flows suggests that the supply of DOC-complexed Hg T 

### from the shallow groundwater system predominates under base flow conditions, when rapid runoff from VSAs and DOC flushing potential from upland soils is minimized.

[ 46 

### ] The TOPLOAD model provides insights on source flow contributions of Hg T 

### . Based on TOPLOAD simulations, saturated subsurface flow drives the daily curve of the HgT loading (TLOAD-lower; Figure 5). This is supported by earlier interpretations of McTier Creek watershed ’ s Coastal Plain system [Bradley et al. , 2010]. However, above base flow TOPLOAD attributes Hg T 

### loads primarily to shallow subsurface flows (TLOAD-upper under non-flood conditions; TLOAD-upper and TLOAD-infiltration under flooded conditions) (Figure 5).

[ 47 

### ] TOPLOAD simulations were also applied to evaluate the hypothesis, based on VELMA simulations, suggesting that when Hg T 

### concentrations are increased in the model flow component that represents direct precipitation on water ( Qrip 

### ; results in higher concentrations in precipitation compared to surface waters), peaks in Hg T,TOPLOAD 

### are higher than those from simulations where direct rainfall Hg T 

### concentrations are equal that of all other subsurface components(see Figure S6 in the auxiliary material). The close relationship of Hg T 

### , VELMA 

### concentrations and fluxes with Hg P 

### , OBS 

### (Figure 3 and 4), for example, supports this conceptualization of ready transport of precipitation derived dissolved Hg T 

### to the stream during rainfall events.
4.3. Insights From Seasonal Variations in Hg T 

### Concentrations and Fluxes

[ 48 

### ] Strong seasonal trends in the simulations and data were not evident for the two-year study period with the exception of higher than average Hg T,OBS 

### and Hg F,OBS 

### based on observed data only. LOADEST and TOPLOAD simulations exhibited the most explicit seasonal pattern of Hg T 

### concentrations or fluxes among the four simulations; however, these LOADEST estimates were based on a linear statistical relationship between observed and predicted fluxes with sine curves assigned to the model for estimating fluxes. Thus, the seasonal pattern may partially represent the structure of the regression model rather than the “true ” Hg T 

### concentration and flux dynamics. Because Hg T,LOADEST 

### concentrations mirrored trends in Hg F,LOADEST 

### concentrations, seasonal results reflected daily results, suggesting that the dissolved fraction of Hg T 

### is largely regulating average seasonal Hg T 

### concentrations in McTier Creek watershed. Further, Hg T 

### , VELMA 

### seasonal simulations are typically slightly higher than Hg P,LOADEST 

### concentrations and seasonal Hg T 

### , GBMM 

### concentrations, which also confirmed that the dissolved fraction of Hg T 

### remains most important in this watershed at a seasonal time scale.
4.4. Advancing Current Watershed Modeling Approaches for Understanding Hg T 

### Fate and Transport

[ 49 

### ] Multiple model and data results suggest that based onthe current state of the science of watershed scale Hg modeling, a multimodel, complementary approach to understanding the dissolved HgT dynamics in McTier Creek watershed is desirable. Ultimately, we suggest development of a spatially explicit watershed Hg model that incorporates, at minimum, both erosion and DOC complexation to explicitly simulate Hg T 

### loadings to surface waters. This is consistent with empirical models developed in other stream ecosystems [e.g., Brigham et al. , 2009; Yin and Balogh

, 2002]. However, if only one of these processes dominates a system, then GBMM or VELMA-Hg may perform well as a stand-alone model. We suggest that (1) GBMM would be most reliably applied as the primary watershed fate and transport model in catchments with highly erodible land, because Hg T 

### transport would most likely be highly associated with sediment transport in these settings [Hurley et al.

, 1995; Munthe et al. , 2007], and (2) VELMA-Hg would be more suitable for forested catchments where DOC is the primary driver of Hg T 

### in-stream concentrations and watershed fluxes.

[ 50 

### ] Additional processes, however, are required in the models to fully capture Hg T 

### dynamics using watershed-scale simulation and include: (1) Hg methylation processes for estimating MeHg, (2) sulfate (SO 4 2

 

### ) cycling, which is known to influence methylation of Hg and fluxes of MeHg from watersheds [ Morel et al. , 1998], (3) variables that increase the availability of other Hg species (e.g., pH,quality and size of organic matter pool, and iron availability), (4) physical factors such as watershed size, topography, and land cover, which regulate the amount of Hg transport GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

14 of 17

### when watershed connectivity to the stream is high [ Munthe et al. , 2007], (5) wetland Hg processing, including DOC sinks and methylation potentials, (6) temperature-varied rate constants, and (7) in-stream processing, including additional particulate dynamics that influence Hg T 

### concentrations andfluxes (such as wetland particulate transport, benthic sediment re-suspension, and stream bank erosion). Models of Hg T 

### dynamics are also complicated by imperfections inhydrological models. Therefore, improvements in hydrological mechanisms in the models are also needed, including: (1) external links to groundwater models, particularly for Coastal Plain systems [ Lu et al. , 2006]; and (2) recent advances in hydrologic modeling, which could potentially decrease variability in flux estimates due to errors in runoff simulation. Finally, temporal gaps in data may not fully capture Hg T 

### concentration-flow dynamics, which informs modeling efforts. For example, our current data set is based on synoptic sampling (n = 41) during a variety of flow conditions. While this level of sampling may be desirablefor statistical methods (i.e., due to serial correlation), highfrequency sampling is beneficial for interpreting complex watershed processes. Our results therefore highlight the benefit of additional monitoring of in-stream, out-of-channel, and upland Hg T 

### processing in Coastal Plain systems and inclusion of these processes in watershed Hg modeling.



### 5. Summary and Conclusions

[ 51 

### ] We present one of the first studies using multiplemodeling techniques and empirical data analysis concurrently to explore Hg cycling and Hg modeling in a Coastal Plain watershed. This study expands beyond the dynamicsrepresented in a single model or data set alone by simultaneously analyzing two spatially explicit watershed Hg models (GBMM, VELMA-Hg), an empirical model driven by TOPMODEL hydrology (TOPLOAD), a seasonal load estimator model (LOADEST), and observed data to exploreHg dynamics in a diverse, mesoscale, Coastal Plain watershed. Using this complementary, value-added modeling approach, we (1) explore various conceptualizations of watershed Hg cycling in McTier Creek Watershed and begin to elucidate the important factors

— represented in the combined multiple watershed Hg models and data

— influencing watershed Hg fluxes and concentrations at the basin outlet and (2) highlight research needed to improve modeling and characterization of watershed Hg dynamics.

[ 52 

### ] We found similar rainfall-runoff results for three models with different conceptualizations of upland watershed hydrologic dynamics, highlighting the need for a “modified” interpretation of peak flow hydrology in this Coastal Plain system. Specifically, this modified interpretation considered the additional importance of shallow groundwater flow paths that respond similarly to high gradient transport in upland systems during and following precipitation events. Consistentwith this conceptualization, results suggest that shallow subsurface flow and subsequent discharge, and to a limited extent overland flow, are potentially important controls on the mobilization of particulate HgT 

### during periods of high uplandto-stream connectivity. Other processes, such as stream bank erosion and benthic sediment re-suspension, may potentially lead to increased particulate HgT 

### in the water column. VELMA-Hg simulations explained a portion of dissolved aqueous HgT 

### , particularly that HgT 

### can be directly transportedfrom VSAs in the watersheds to surface waters during highrainfall events. Observed data suggested that flushing of DOCHgT 

### complexes from surface soils may also be important during this period, which aligns with the “modified” hydrological interpretation of watershed dynamics that focuses on responsive VSA and shallow subsurface contributions of HgT 

### . DOC-bound HgT 

### becomes much more dominant during baseflow conditions, when DOC is the primary mechanism transporting Hg T 

### to the stream. TOPLOAD results suggest that HgT 

### in the stream during base flow conditions is predominately from saturated subsurface flows, or groundwater sources.However, during flows peaking above base flow, contributions from infiltration excess overland flow and unsaturated subsurface flows create peak Hg T 

### load conditions. Finally, strong seasonal trends in the simulated and observed data werenot evident for the two-year simulation period; however, seasonal analysis confirmed that the dissolved fraction of HgT 

### remains most important in this watershed.

[ 53 

### ] Based on collected data and model simulations, we also suggest that a spatially explicit watershed Hg model thatincorporates, at minimum, both erosion and DOC complexation is needed to reliably simulate Hg T 

### loadings to Coastal Plain surface waters. However, improvements tocurrent watershed Hg models are needed including additional biogeochemical processes (e.g, sulfate), wetland cycling, links to groundwater models, and advancements in hydrological modeling.

[ 54 

### ] Our findings demonstrate advancement in watershed Hg research and modeling in Coastal Plain watersheds. The gaps in knowledge we identify provide a potential stepping stone for important questions concerning Hg science and modeling at the watershed scale.

[ 55 

### ] Acknowledgments. We thank Alex Abdelnour, Marc Steiglitz, Bob McKane, and Fei Fei Pan for assistance in developing the VELMA-Hgmodel. We appreciate helpful suggestions from Mark Gabriel and two anonymous reviewers. This paper has been reviewed in accordance with the U.S. Environmental Protection Agency ’ s and the U.S. Geological Survey ’ s peer and administrative review policies and approved for publication. Mention oftrade names or commercial products does not constitute endorsement or recommendation by the U.S. Government. Statements in this publication reflect the authors ’ professional views and opinions and should not be construed to represent any determination or policy of the U.S. Environmental Protection Agency but do represent the views of the U.S. Geological Survey.



### References
Abdelnour, A., M. Stieglitz, F. Pan, and R. McKane (2011), Catchment

hydrological responses to forest harvest amount and spatial pattern, Water Resour. Res.

, 47, W09521, doi:10.1029/2010WR010165. Aiken, G. R., D. M. McKnight, K. A. Thorn, and E. M. Thurman (1992),

Isolation of hydrophilic organic acids from water using nonionic macroporous resins, Org. Geochem.

, 18 , 567

–573, doi:10.1016/0146-6380(92) 90119-I.Atkins, J. B., C. A. Journey, and J. S. Clarke (1996), Estimation of groundwater discharge to streams in the central Savannah River basin of Georgia and South Carolina, U.S. Geol. Surv. Water Resour. Invest. Rep.

, 96

–4179

, 41 pp. [Available at http://pubs.usgs.gov/wri/1996/4179/report.pdf.] Barringer, J. L., M. L. Riskin, Z. Szabo, P. A. Reilly, R. Rosman, J. L.

Bonin, J. M. Fischer, and H. A. Heckathorn (2010), Mercury and methylmercury dynamics in a coastal plain watershed, New Jersey, USA, Water Air Soil Pollut.

, 212 , 251

–273, doi:10.1007/s11270-010-0340-1. Benedict, S. T., P. A. Conrads, T. D. Feaster, C. A. Journey, H. E. Golden,

C. D. Knightes, G. M. Davis, and P. M. Bradley (2011), Data visualization, time-series analysis, and mass-balance of hydrologic and waterquality data for McTier Creek, South Carolina, 2007

–2009, U.S. Geol. Surv. Sci. Invest. Rep.

, 2011 -, in press. Bennett, D. G., and J. C. Patton (Eds.) (2008), A Geography of the Carolinas , 266 pp., Parkway, Boone, N. C. GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

15 of 17Beven, K. J., and M. J. Kirkby (1979), A physically based, variable contributing area model of basin hydrology, Hydrol. Sci. Bull.

, 24 , 43

–69, doi:10.1080/02626667909491834. Bloom, N. S. (1992), On the chemical form of mercury in edible fish and

marine invertebrate tissue, Can. J. Fish. Aquat. Sci.

, 49 , 1010

– 1017, doi:10.1139/f92-113. Boyer, E. W., G. M. Hornberger, K. E. Bencala, and D. M. McKnight

(1997), Response characteristics of DOC flushing in an alpine catchment, Hydrol. Processes

, 11 , 1635

– 1647, doi:10.1002/(SICI)1099-1085 (19971015)11:12<1635::AID-HYP494>3.0.CO;2-H. Bradley, P. M., C. A. Journey, F. H. Chapelle, M. A. Lowery, and P. A.

Conrads (2010), Flood hydrology and methylmercury availability in Coastal Plain rivers, Environ. Sci. Technol.

, 44 (24), 9285

–9290, doi:10.1021/es102917j. Bradley, P. M., D. A. Burns, K. Riva-Murray, M. E. Brigham, D. Button,

L. Chasar, M. Marvin-DiPasquale, M. A. Lowery, and C. A. Journey (2011), Spatial and seasonal variability of dissolved methylmercury in two stream basins in the eastern United States, Environ. Sci. Technol.

, 45 (6), 2048

–2055, doi:10.1021/es103923j. Brigham, M. E., D. A. Wentz, G. R. Aiken, and D. P. Krabbenhoft (2009),

Mercury cycling in stream ecosystems. 1. Water column chemistry and transport, Environ. Sci. Technol.

, 43 (8), 2720

– 2725, doi:10.1021/es802694n. Brumbaugh, W. G., D. P. Krabbenhoft, D. R. Helsel, J. G. Wiener, and

K. R. Echols (2001), A national pilot study of mercury contamination of aquatic ecosystems along multiple gradients: Bioaccumulation in fish, U.S. Geol. Surv. Biol. Sci. Rep., USGS/BRD/BSR-2001-0009, 32 pp. [Available at http://www.cerc.usgs.gov/pubs/center/pdfdocs/BSR2001-0009.pdf.] Burgess, N., and M. Meyer (2008), Methylmercury exposure associated

with reduced productivity in common loons, Ecotoxicology

, 17 (2), 83

–91, doi:10.1007/s10646-007-0167-8. Bushey, J. T., C. T. Driscoll, M. J. Mitchell, P. Selvendiran, and M. R.

Montesdeoca (2008), Mercury transport in response to storm events from a northern forest landscape, Hydrol. Processes

, 22 (25), 4813

–4826, doi:10.1002/hyp.7091. Campbell, B. G., and A. L. Coess (Eds.) (2010), Groundwater availability

in the Atlantic Coastal Plain of North and South Carolina, U.S. Geol. Surv. Prof. Pap.

, 1773 , 241 pp. [Available at http://pubs.usgs.gov/pp/ 1773/.] Clarkson, T. W., and L. Magos (2006), The toxicology of mercury and its

chemical compounds, Crit. Rev. Toxicol.

, 36(8), 609

–662, doi:10.1080/ 10408440600845619. Cohn, T. A., D. L. Caulder, E. J. Gilroy, L. D. Zynjuk, and R. M. Summers

(1992), The validity of a simple statistical model for estimating fluvial constituent loads

—An empirical study involving nutrient loads entering Chesapeake Bay, Water Resour. Res.

, 28 (9), 2353

–2363, doi:10.1029/ 92WR01008. Dai, T., R. B. Ambrose, K. Alvi, T. Wool, H. Manguerra, M. Chokshi,

H. Yang, and S. Kraemer (2005), Characterizing spatial and temporal dynamics: Development of a grid-based watershed mercury loading model, in Managing Watersheds for Human and Natural Impacts: Engineering, Ecological, and Economic Challenges , edited by G. E. Moglen, Am. Soc. of Civ. Eng., Reston, Va., doi:10.1061/40763(178)56.Demers, J. D., C. T. Driscoll, and J. B. Shanley (2010), Mercury mobilization and episodic stream acidification during snowmelt: Role of hydrologic flow paths, source areas, and supply of dissolved organic carbon, Water Resour. Res.

, 46, W01511, doi:10.1029/2008WR007021. DeWild, J. F., M. L. Olson, and S. D. Olund (2002), Determination of

methyl mercury by aqueous phase ethylation, followed by gas chromatographic separation with cold vapor atomic fluorescence detection, U.S. Geol. Surv. Open File Rep.

, 01-445 , 19 pp. [Available at http://pubs. usgs.gov/of/2001/ofr-01-445/pdf/ofr01445v1.pdf.] Evers, D., et al. (2008), Adverse effects from environmental mercury loads

on breeding common loons, Ecotoxicology

, 17(2), 69

– 81, doi:10.1007/ s10646-007-0168-7. Feaster, T. D., H. E. Golden, K. R. Odom, M. A. Lowery, P. A. Conrads,

and P. M. Bradley (2010), Simulation of streamflow in the McTier Creek watershed, South Carolina, U.S. Geol. Surv. Sci. Invest. Rep.

, 2010-5202

, 61 pp. [Available at http://pubs.usgs.gov/sir/2010/5202/pdf/sir2010-5202. pdf.] Glover, J., M. Domino, K. Altman, J. Dillman, W. Castleberry, J. Eidson,

and M. Mattocks (2010), Mercury in South Carolina fishes, USA, Ecotoxicology

, 19 (4), 781

–795, doi:10.1007/s10646-009-0455-6. Golden, H. E., and C. D. Knightes (2011), Simulated watershed mercury

and nitrate flux responses to multiple land cover conversion scenarios, Environ. Toxicol. Chem.

, 30 (4), 773

–786, doi:10.1002/etc.449. Golden, H. E., C. D. Knightes, E. J. Cooter, R. L. Dennis, R. C. Gilliam,

and K. M. Foley (2010), Linking air quality and watershed models for environmental assessments: Analysis of the effects of model-specific precipitation estimates on calculated water flux, Environ. Model. Softw.

, 25 (12), 1722

– 1737, doi:10.1016/j.envsoft.2010.04.015.Grigal, D. (2002), Inputs and outputs of mercury from terrestrial watersheds: A review, Environ. Rev.

, 10(1), 1

–39, doi:10.1139/a01-013.Guentzel, J. L. (2009), Wetland influences on mercury transport and bioaccumulation in South Carolina, Sci. Total Environ.

, 407(4), 1344

– 1353, doi:10.1016/j.scitotenv.2008.09.030. Helsel, D. R., and R. M. Hirsch (1992), Statistical methods in water

resources, in Hydrologic Analysis and Interpretation

, U.S. Geol. Water Resour. Invest. , Book 4, Chap. A3, 512 pp. [Available at http://pubs. usgs.gov/twri/twri4a3/html/pdf_new.html.]Homer, C., C. Huang, L. Yang, B. K. Wylie, and M. Coan (2004), Development of a 2001 National Landcover Database for the United States, Photogramm. Eng. Remote Sens.

, 70 (7), 829

– 840. Hornberger, G., K. E. Bencala, and D. M. Mcknight (1994), Hydrological

controls on dissolved organic carbon during snowmelt in the Snake River near Montezuma, Colorado, Biogeochemistry

, 25 , 147

– 165, doi:10.1007/ BF00024390. Hupp, C. R. (2000), Hydrology, geomorphology and vegetation of Coastal

Plain rivers in the south-eastern USA, Hydrol. Processes

, 14 (16

– 17), 2991

–3010, doi:10.1002/1099-1085(200011/12)14:16/17<2991::AIDHYP131>3.0.CO;2-H. Hurley, J. P., J. M. Benoit, C. L. Babiarz, M. M. Shafer, A. W. Andren, J. R.

Sullivan, R. Hammond, and D. A. Webb (1995), Influences of watershed characteristics on mercury levels in Wisconsin rivers, Environ. Sci. Technol.

, 29 (7), 1867

– 1875, doi:10.1021/es00007a026. Journey, C. A., A. W. Caldwell, T. D. Feaster, M. D. Petkewich, and P. M.

Bradley (2011), Concentrations, loads, and yields of nutrients and suspended sediment in the South Pacolet, North Pacolet, and Pacolet rivers, northern South Carolina and soutwestern North Carolina, October 2005 to September 2009, U.S. Geol. Surv. Sci. Invest. Rep.

, 2010-5252

, 79 pp. [Available at http://pubs.usgs.gov/sir/2010/5252.] Krabbenhoft, D. P., J. G. Wiener, W. G. Brumbaugh, M. L. Olson, J. F.

DeWild, and T. J. Sabin (1999), A national pilot study of mercury contamination of aquatic ecosystems along multiple gradients, in U.S. Geological Survey Toxic Substances Hydrology Program: Proceedings of the Technical Meeting, Charleston, South Carolina, March 8

– 12, 1999

, edited by D. W. Morganwalp and H. T. Buxton, WRI Rep. 99-4018B

, pp. 147

– 160, U.S. Geol. Surv., West Trenton, N. J. Lewis, M. E., and M. E. Brigham (2004), Low-level mercury, in National

Field Manual for the Collection of Water-Quality Data

, U.S. Geol. Surv. Tech. Water Resour. Invest. , Book 9, Chap. A5, 26 pp. [Available at http://water.usgs.gov/owq/FieldManual/chapter5/pdf/5.6.4.B_v1.0.pdf.] Lu, J., G. Sun, D. M. Amatya, S. V. Harder, and S. G. McNulty (2006),

Understanding the hydrologic response of a coastal plain watershed to forest management and climate change in South Carolina, U.S.A., in Hydrology and Management of Forested Wetlands: Proceedings of the International Conference, April 8

– 12, 2006, New Bern, North Carolina

, pp. 231

–239, Am. Soc. of Agric. and Biol. Eng., St. Joseph, Mich. Marshall, W. D. (Ed.) (1993), Assessing change in the Edisto River Basin:

An ecological characterization, South Carolina Water Resour. Comm. Rep.

, 177, 149 pp. [Available at http://www.gpo.gov/fdsys/pkg/CZIC-qh105-s6- a87-1993/pdf/CZIC-qh105-s6-a87-1993.pdf.] Mason, R. P., and K. A. Sullivan (1998), Mercury and methylmercury

transport through an urban watershed, Water Res.

, 32 (2), 321

–330, doi:10.1016/S0043-1354(97)00285-6. Mergler, D., H. A. Anderson, L. Hing Man Chan, K. R. Mahaffey, M. Murray,

M. Sakamoto, and A. H. Stern (2007), Methylmercury exposure and health effects in humans: A worldwide concern, Ambio

, 36 (1), 3

– 11, doi:10.1579/ 0044-7447(2007)36[3:MEAHEI]2.0.CO;2. Morel, F. M. M., A. M. L. Kraepiel, and M. Amyot (1998), The chemical

cycle and bioaccumulation of mercury, Annu. Rev. Ecol. Syst.

, 29, 543

–566, doi:10.1146/annurev.ecolsys.29.1.543. Munthe, J., R. A. Bodaly, B. A. Branfireun, and C. T. Driscoll (2007),

Recovery of mercury-contaminated fisheries, Ambio

, 36(1), 33

–44, doi:10.1579/0044-7447(2007)36[33:ROMF]2.0.CO;2.Nash, J. E., and J. V. Sutcliffe (1970), River flow forecasting through conceptual models, Part I

—

A discussion of principles, J. Hydrol.

, 10 , 282

–290, doi:10.1016/0022-1694(70)90255-6. National Research Council (2007), Models in Environmental Decision

Making, 287 pp., Natl. Acad. Press, Washington, D. C. Olson, M. L., L. B. Cleckner, J. P. Hurley, D. P. Krabbenhoft, and T. W.

Heelan (1997), Resolution of matrix effects on analysis of total and methyl mercury in aqueous samples from the Florida Everglades, Fresenius, Fresenius ’ J. Anal. Chem.

, 358(3), 392

–398, doi:10.1007/s002160050435. Olund, S. D., J. F. DeWild, M. L. Olson, and M. T. Tate (2004), Methods

for the preparation and analysis of solids and suspended solids for total mercury, in Laboratory Analysis

, U.S. Geol. Surv. Tech. Methods , Book 5, GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

16 of 17Chap. A8, 31 pp. [Available at http://pubs.usgs.gov/tm/2005/tm5A8/pdf/ TM5A-8.pdf.] Robson, A., K. Beven, and C. Neal (1992), Towards identifying sources of

subsurface flow: A comparison of components identified by a physically based runoff model and those determined by chemical mixing techniques, Hydrol. Processes

, 6, 199

–214, doi:10.1002/hyp.3360060208. Runkel, R. L., C. G. Crawford, and T. A. Cohn (2004), Load Estimator

(LOADEST): A Fortran program for estimating constituent loads in streams and rivers, in Hydrologic Analysis and Interpretation

, U.S. Geol. Surv. Tech. Methods , Book 4, Chap. A5, 75 pp. [Available at http://pubs. usgs.gov/tm/2005/tm4A5/pdf/508final.pdf.] Scudder, B. C., L. C. Chasar, D. A. Wentz, N. J. Bauch, M. E. Brigham,

P. W. Moran, and D. P. Krabbenhoft (2009), Mercury in fish, bed sediment, and water from streams across the United States, 1998

–2005, U.S. Geol. Surv. Sci. Invest. Rep.

, 2009-5109 , viii + 74 pp. [Available at http://pubs.usgs.gov/sir/2009/5109/pdf/sir20095109.pdf.]Shreve, E. A., and A. C. Downs (2005), Quality-assurance plan for the analysis of fluvial sediment by the U.S. Geological Survey Kentucky Water Science Center Sediment Laboratory, U.S. Geol. Surv. Open File Rep.

, 2005-1230 , 34 pp. [Available at http://pubs.usgs.gov/of/2005/1230/ ofr20051230.pdf.] U.S. Environmental Protection Agency (2005), Regulatory impact analysis

of the Clean Air Mercury Rule, final report, Rep. EPA-452/R-05

– 003 , U. S. Environ. Prot. Agency, Research Triangle Park, N. C. Wolock, D. M. (1993), Simulating the variable-source-area concept of streamflow generation with the watershed model TOPMODEL, U. S. Geol. Surv. Water Resour. Invest. Rep.

, 93-4124 , 39 pp. [Available at http://pubs.usgs.gov/wri/1993/4124/report.pdf.] Yin, X., and S. J. Balogh (2002), Mercury concentrations in stream and

river water: An analytical framework applied to Minnesota and Wisconsin (USA), Water, Air, Soil Pollut.

, 138 (1

–4), 79

– 100, doi:10.1023/A: 1015525729671.

S. T. Benedict and T. D. Feaster, South Carolina Water Science Center, U.S. Geological Survey, Clemson Field Office, Ste. 200, 405 College Ave., Clemson, SC 29631, USA.

P. M. Bradley, P. A. Conrads, and C. A. Journey, South Carolina Water Science Center, U.S. Geological Survey, Columbia Ste. 129, 720 Gracern Rd., Columbia, SC 29210, USA.

M. E. Brigham, Minnesota Water Science Center, U.S. Geological Survey, 2280 Woodale Dr., Mounds View, MN 55112, USA.

G. M. Davis and C. D. Knightes, Ecosystems Research Division, Office of Research and Development, U.S. Environmental Protection Agency, 960 College Station Rd., Athens, GA 30605, USA.

H. E. Golden, Ecological Exposure Research Division, Office of Research and Development, U.S. Environmental Protection Agency, 26 W. Martin Luther King Dr., Cincinnati, OH 45268, USA. (Golden.Heather@epa.gov) GOLDEN ET AL.: COASTAL PLAIN WATERSHED MERCURY DYNAMICS G01006 G01006

17 of 17