Simulating reactive transport of selenium coupled with nitrogen in a regional-scale irrigated groundwater system

Journal of Hydrology 515 (2014) 29–46

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Simulating reactive transport of selenium coupled with nitrogen in a regional-scale irrigated groundwater system Ryan T. Bailey ⇑, Timothy K. Gates, Mehdi Ahmadi 1 Department of Civil and Environmental Engineering, Colorado State University, 1372 Campus Delivery, Fort Collins, CO 80523-1372, United States

a r t i c l e

i n f o

Article history: Received 9 May 2013 Received in revised form 24 February 2014 Accepted 14 April 2014 Available online 26 April 2014 This manuscript was handled by Laurent Charlet, Editor-in-Chief, with the assistance of Ewen James Silvester, Associate Editor Keywords: Selenium Groundwater contaminant transport modeling Nitrogen cycling Agricultural groundwater system Marine shale

s u m m a r y The contamination or deficiency of Selenium (Se) has become a critical issue in watershed systems worldwide during the past half-century, with either elevated or deficient Se concentrations in groundwater, surface water, soils, and associated vegetation. Regardless of the reason for concern, there is a basic need for tools for use in assessing baseline concentrations and mass loadings in large-scale environmental systems and for exploring remediation strategies. This paper presents the application of UZF-RT3D, a multispecies variably-saturated reactive transport model for groundwater systems, and its recently developed Se agricultural and chemical reaction module to a regional study site (50,600 ha) along the irrigated Arkansas River Valley in southeastern Colorado. The study region has been monitored for Se contamination during the previous decade, with all segments of the Arkansas River impaired with respect to Se. The Se reaction module accounts for near-surface Se cycling due to agricultural practices, oxidation–reduction reactions, and sorption, and includes a nitrogen cycle and reaction module due to the dependence of Se transformation and speciation on the presence of nitrate (NO3). Of prime importance is the role NO3 plays in oxidation of residual Se from marine shale, which is prevalent throughout the study region, along with its inhibition of Se chemical reduction to less toxic forms. The model is corroborated against groundwater solute concentrations, mass loadings of solutes to the Arkansas River, relationships between solutes within the groundwater system, and overall regional statistics of groundwater solute concentration. Results indicate that the model is successful in replicating the major spatial and temporal patterns in Se and NO3 contamination and transport, hence providing a tool that can be used with confidence to explore remediation schemes. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction The presence of elevated Selenium (Se) concentrations in aquifer and surface water systems has emerged as a serious concern in regions worldwide during the last half-century (Seiler, 1997; Afzal et al., 2000; Mizutani et al., 2001; Zhang et al., 2008; Hudak, 2010). Although an essential micro-nutrient for humans and animals (Combs et al., 1986; Aro et al., 1998), a narrow range exists between dietary deficiency and toxic levels. In recent decades, this ‘‘essential toxin’’ (Stolz et al., 2002; Fernández-Martínez and Charlet, 2009) has caused deformities and death among water fowl and fish populations in surface waters hydraulically connected to contaminated aquifer systems (Flury et al., 1997; Hamilton, ⇑ Corresponding author. Tel.: +1 970 491 5387; fax: +1 970 491 7727. E-mail addresses: [email protected] (R.T. Bailey), [email protected]. edu (T.K. Gates), [email protected] (M. Ahmadi). 1 Present address: Spatial Science Laboratory, Texas Agrilife Research, Texas A&M University, 1500 Research Plaza, College Station, TX 77843-2120, United States. http://dx.doi.org/10.1016/j.jhydrol.2014.04.039 0022-1694/Ó 2014 Elsevier B.V. All rights reserved.

1998; Skorupa, 1998). In many instances, Se contamination occurs in irrigated agricultural aquifer systems underlain by shale from which Se is released via oxidation–reduction processes (Seiler, 1995, 1997; Bailey et al., 2012). Conversely, Se deficiency in soils and crops has been linked to diseases among both animal (Aro et al., 1998; Kishchak, 1998) and human (Alfthan et al., 1995; Peng et al., 1995; Wang and Gao, 2001) populations. Whether the nature of the problem is one of elevated or of deficient concentrations, methodologies and tools need to be developed that can support evaluation of the current level of concentration for a particular groundwater system and exploration of remediation strategies. Due to the spatial variability in land management practices and geologic formations, both of which strongly influence the level of Se contamination (Seiler, 1995), analysis at the regional scale (104–105 ha) is essential. Whereas spot sampling of groundwater and surface water can provide a general assessment of the concentration for a particular region (e.g., Gates et al., 2009; Hudak, 2010), other techniques are required to determine the cause of Se availability and movement

30

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

within the aquifer system, and to identify potential remediation strategies. Due to the complexity of Se fate and transport dynamics in an aquifer system, which include oxidation–reduction (redox) reactions with dissolved oxygen (O2) and nitrate (NO3) (e.g., Weres et al., 1990; White et al., 1991; Bailey et al., 2012), sorption, along with spatio-temporal mass inputs and sinks, the development and use of numerical models is an appealing approach. Numerical models designed for Se fate and transport dynamics thus far have been limited to one-dimensional (1D) soil profile models (Alemi et al., 1988, 1991; Fio et al., 1991; Liu and Narasimhan, 1994; Guo et al., 1999; Mirbagheri et al., 2008) and a small-scale two-dimensional (2D) vertical profile model (Tayfur et al., 2010). Se transport has been simulated in a regional-scale watershed using a three-dimensional (3D) model (Myers, 2013), although Se was treated as a conservative solute. For the 1D studies, Se transport was subject to sorption and redox reactions, while boundary conditions and forcing terms were kept basic (e.g., steady inflow rate and constant concentration of Se within the influent). Mirbagheri et al. (2008) incorporated a more complete set of Se cycling processes for a 1D model, including volatilization, mineralization and immobilization, and plant uptake. Tayfur et al. (2010) applied the same modeling processes to a small-scale (3-m soil profiles) 2D vertical cross-section to simulate Se transport in both saturated and unsaturated soil zones. More recently, the model used by Bailey et al. (2013a) in assessing Se transport in soil profiles at a test site in southeastern Colorado accounted for Se cycling in the root and soil zone and the dependence of redox reactions on the presence of organic carbon (OC). The model also accounted for the reactive transport of O2 and NO3, with the latter governed by Nitrogen (N) cycling. Modeling framework consisted of incorporating Se and N modules into UZF-RT3D (Bailey et al., 2013b), which is linked with MODFLOW-NWT (Niswonger et al., 2011) and simulates the reactive transport of multiple interacting chemical species in variably-saturated groundwater systems. Although chemical reactions governing the release of Se from marine shale have been implemented in the modeling framework, the model has until now not been applied at a large enough scale to require their use. There are two objectives of this paper. The first is to present a numerical modeling tool that can be used to simulate the reactive transport of Se species in a regional-scale irrigated agricultural groundwater system underlain by marine shale. The second objective is to test the tool to establish its credibility in evaluating remediation strategies. The tool was developed by applying UZF-RT3D and the accompanying Se and N reaction modules developed by Bailey et al. (2013a) to a regional study site along the alluvial valley of the Arkansas River in southeastern Colorado, with the groundwater flow pattern and sources/sinks provided by a calibrated and tested groundwater flow model for the same region (Morway et al., 2013). Of particular interest is the influence of marine shale, which underlies the alluvial aquifer and also is present as outcrops and weathered residuum throughout the region, on the source and mobilization of Se in the aquifer. Spatial distributions of irrigation water, canal and river seepage, and crop cultivation, each with associated Se concentration and Se mass residue, also are included. Guided by global sensitivity analysis, the model is calibrated and tested against observed groundwater solute concentration, calculated daily mass loading of solutes from the aquifer and ungauged surface flows to the Arkansas River, observed relationships between solutes in the groundwater, and regional statistics of groundwater solute concentrations. Following a summary of the Se and N reaction modules for UZF-RT3D, the climate, geology, and field monitoring network of the study region is presented, followed by a description of the calibration and testing of the model for the study region.

2. Simulating Se and N reactive transport at the regional scale 2.1. Conceptual model of Se and N reactive transport Se is present in water systems primarily in the four oxidation 2 states of +VI (selenate SeO2 4 ),+IV (selenite SeO3 ), 0 (elemental 0 2 selenium Se ), and -II (selenide Se ). Selenide can occur in forms such as solid Se found in geologic formations in the form of seleno-pyrite (-II) (FeSexS2x), in which Se substitutes for Sulfur (S) in pyrite (FeS2) (Bye and Lund, 1982), organic selenomethionine (-II) (SeMet), and gaseous dimethylselenide (CH3)2Se (-II) (DMSe, a product of the volatilization of SeMet). SeO4, SeO3, and SeMet are soluble, whereas Se0 and other forms of Se2 are insoluble. SeO4 is one of the most toxic Se species, is a weak sorbent (Ahlrichs and Hossner, 1987), and has been reported to account for 90–95% of soluble Se in oxygenated agricultural waters (Masscheleyn et al., 1989; Gates et al., 2009; Gerla et al., 2011), and hence is the primary target in Se contamination remediation. Microbially-mediated chemical reactions (Oremland et al., 1990) reduce SeO4 to SeO3 and SeO3 to either immobile Se0 or mobile SeMet (Ellis and Salt, 2003). SeMet can then be volatilized to non-toxic DMSe (Calderone et al., 1990; Frankenberger and Arshad, 2001). SeO3 also can sorb strongly to soil surface sites (Balistrieri and Chao, 1987). Chemical reductions, however, are inhibited by the presence of O2 and NO3 (Oremland et al., 1989; Weres et al., 1990; White et al., 1991; Zhang and Moore, 1997). A more detailed description of these processes is provided in Bailey et al. (2013a). SeO4 also can be released from shale formations containing seleno-pyrite (Seiler, 1995) via autotrophic reduction of O2 or NO3 (Wright, 1999; Stillings and Amacher, 2010): 2 þ 2FeSex S2x þ 7O2 þ 2H2 O ! 2Fe2þ þ 2xSeO2 4 þ ð4  2xÞSO4 þ 4H

ð1aÞ 2 5FeSex S2x þ 14NO3 þ 4Hþ ! 5Fe2þ þ 5xSeO2 4 þ ð10  5xÞSO4

þ 7N2 þ 2H2 O

ð1bÞ

Se fate and transport in an irrigated alluvial stream-aquifer system is represented in Fig. 1A. This conceptualization is the result of numerous field studies conducted during the previous decades. Sources of Se in watersheds like this are Se-bearing geologic formations, particularly marine shale (Williams and Byers, 1934; Moxon et al., 1938; Schultz et al., 1980; Kulp and Pratt, 2004). This residual Se in either outcropped or bedrock shale is subject to oxidation by O2 in leaching soil water and groundwater and, more recently, by NO3-laden water resulting from over-fertilization and over-irrigation of cultivated area. Once released from shale, typically in the form of SeO4, Se is transported through the aquifer to either shallow soil zones or to surface water discharge points (Fig. 1A), with transport tempered by possible sorption and redox reactions that lead to immobilized, volatilized, and/or precipitated forms of Se. In regions downstream from Se source locations, Se-laden river water may be diverted into canal networks (Fig. 1A), resulting in Se transport to the root zone and shallow soil zones via infiltrating irrigation water and canal seepage (Fig. 1A). If present in the root zone in cultivated areas, Se is subject to cycling in the plant-soil system similar to other nutrients (Shrift, 1954; Stolz et al., 2002). The conceptual model of Se fate and transport in an agricultural soil system is shown in Fig. 1B. All redox reactions as discussed previously are shown, including the chemical reduction of O2 and NO3 (denitrification), and the nitrification and volatilization of ammonium (NH4). Soil organic matter is composed of litter (fast-decomposing), humus (slow-decomposing), and manure. Solute mass is taken up by crop roots throughout the growing season and incorporated into the plant material. At harvest the dead root mass of the crop, with associated organic

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

31

Fig. 1. (A) Conceptual model of the movement, transformation, and storage of Se and N species in an irrigated agricultural watershed and (B) detailed root zone processes and chemical reactions and transformations.

Se and N, is assimilated into the litter pool, and at plowing the remainder of the root mass and the above-ground crop mass not removed at harvest (i.e., stover) is incorporated into the litter pool, upon which decomposition of the organic matter occurs. Mineralization of organic matter to inorganic species and immobilization of inorganic species to organic matter also occur. For simplicity of notation, NH4–N, NO3–N, SeO4–Se, SeO3–Se will be written as NH4, NO3, SeO4, and SeO3 throughout the remainder of this paper. 2.2. Se and N cycling and reaction modules for UZF-RT3D This section describes the development of the regional-scale Se transport model. Following a brief description of UZF-RT3D and the developed Se and N reaction modules, considerations for applying the model to a large-scale groundwater system are discussed. 2.2.1. Reactive transport of Se and N in agricultural groundwater systems The Se and N reactive transport model (Bailey et al., 2013a) was developed by incorporating Se and N reaction modules into the base reactive transport model UZF-RT3D (Bailey et al., 2013b), which is linked with the Unsaturated-Zone Flow (UZF1) package (Niswonger et al., 2006) of MODFLOW-NWT (Niswonger et al., 2011) and which simulates the reactive transport of multiple interacting species in variably-saturated soil and groundwater systems. UZF1 assumes vertical homogeneity of the unsaturated zone and neglects the diffusive term in the Richards equation, thus demanding less computational effort and data support than models that solve the full Richards equation and hence providing a suitable tool

for simulating variably-saturated flow and transport at the regional scale. UZF-RT3D solves a system of advection-dispersion-reaction (ADR) equations for dissolved-phase and solid-phase species (Bailey et al., 2013b) using the operator-split strategy (Yeh and Tripathi, 1989), with the kinetics of multiple interacting species simulated using an ordinary differential equation solver. Using the form of the ADR equation in UZF-RT3D, the following equations are written for the Se dissolved-phase species:

  @ðC SeO4 hÞ @ @ @C SeO4 RSeO4 ¼  ðhv i C SeO4 Þ þ hDij @t @xi @xi @xj   imm þ qf C fSeO þ F SeO4  U SeO4 þ e r min s;Se  r s;Se 4   het þ h rauto f ;SeO4  r f ;SeO4

ð2aÞ

  @ðC SeO3 hÞ @C SeO3 @ @ ðhv i C SeO3 Þ þ hDij RSeO3 ¼  @xi @xi @xj @t þ qf C fSeO þ F SeO3  U SeO3 3   hetðSes Þ hetðSeMetÞ þ h rhet f ;SeO4  r f ;SeO3  r f ;SeO3   @ðC SeMet hÞ @ @ @C SeMet þ qf C fSeMet ¼ ðhv i C SeMet Þ þ hDij @t @xi @xi @xj   hetðSeMetÞ  U SeMet þ h r f ;SeO3  rhet f ;SeMet

ð2bÞ

ð2cÞ

where C is solute concentration [MfLf3], with f denoting fluid phase; Dij is the hydrodynamic dispersion coefficient [L2T1]; v is

32

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

the pore velocity [LbT1] provided by MODFLOW-UZF1; / is porosity [Lf3Lb3]; h is the volumetric water content [Lf3Lb3]; qf is the volumetric flux of water representing sources and sinks [Lf3T1Lb3] such as irrigation water, canal and river seepage, groundwater discharge to the river, or pumped groundwater, with b denoting the bulk phase; Cf is the concentration of the source or sink [MfLf3]; rf represents the rate of all reactions that occur in the dissolved1 phase [MfL3 ]; e is the volumetric solid content [Ls3Lb3] with s f T denoting the solid phase, and is equal to 1  /; Rj is the retardation factor for species j and is equal to 1 + (qb Kd,j)/h, where qb is the bulk density of the porous media [MbLb3] and is the partitioning coeffi1 cient [Lf3Mb]; F is the inorganic fertilizer application [MfL3 ]; U b T 1 is the potential uptake rate [MfL3 T ]; min and imm signify minerb alization and immobilization, respectively; auto and het represent autotrophic and heterotrophic chemical reduction, respectively. Similar equations are written for NH4, NO3, and O2, with nit indicating nitrification:

  @ðC NH4 hÞ @C NH4 @ @ þ qf C fNH4 ðhv i C NH4 Þ þ hDij RNH4 ¼  @xi @xi @xj @t   imm þ F NH4  U NH4 þ e r min s;N  r s;N   v ol þ h r nit f  rf   @ðC NO3 hÞ @C NO3 @ @ þ qf C fNO ¼ ðhv i C NO3 Þ þ hDij 3 @xi @xi @xj @t   nit het auto þ F NO3  U NO3 þ h r f  r f ;NO3  rf ;NO3   @ðC O2 hÞ @ @ @C O2 þ q f C fO ¼ ðhv i C O2 Þ þ hDij 2 @t @xi @xi @xj   het auto þ h r f ;O2  r f ;O2

    @ðC HSe eÞ dec imm min ¼ e r dec s;SeðL!HÞ  r s;SeðHÞ þ e r s;SeðHÞ  r s;SeðHÞ @t   @ðC MSe eÞ imm min ¼ M Se  erdec s;SeðMÞ þ r s;SeðMÞ  r s;SeðMÞ @t

het rhet f ;O2 ¼ kO2 C O2



C O2 K O2 þ C O2

het rhet f ;NO3 ¼ kNO3 C NO3



 CO2;prod E K CO2 þ CO2;prod



C NO3 K NO3 þ C NO3



ð5aÞ

CO2;prod K CO2 þ CO2;prod



 I O2 E IO2 þ C O2 ð5bÞ

het rhet f ;SeO4 ¼ kSeO4 C SeO4



CO2;prod K CO2 þ CO2;prod



I O2 I O2 þ C O2



 INO3 E INO3 þ C NO3 ð5cÞ

hetðSe Þ

hetðSe Þ

rf ;SeO3s ¼ kSeO3 s C SeO3



CO2;prod K CO2 þ CO2;prod



I O2 IO2 þ C O2



 INO3 E INO3 þ C NO3 ð5dÞ

hetðSeMetÞ

rf ;SeO3

hetðSeMetÞ

¼kSeO3

C SeO3



ð3aÞ

CO2;prod K CO2 þCO2;prod



I O2 IO2 þC O2



 INO3 E INO3 þC NO3 ð5eÞ

het rhet f ;SeMet ¼ kSeMet C SeMet



CO2;prod K CO2 þ CO2;prod



IO2 I O2 þ C O2



 INO3 E INO3 þ C NO3 ð5fÞ

ð3bÞ 1

ð3cÞ

A glossary for all model parameters is included in Bailey et al. (2013a). For both Se and N, the daily mass of uptake during the growing season is calculated using a logistic equation (Johnsson et al., 1987) and is distributed across a profile of grid cells according to the depth and mass distribution of the root system (Neitsch et al., 2011). Fate equations also are written for the Se solid-phase litter LSe, humus, HSe, and manure MSe organic matter:

@ðC LSe eÞ ¼ aRt;Se PRt þ aSt;Se PSt @t   dec dec dec þ e rdec s;SeðH!LÞ þ r s;SeðM!LÞ þ r s;SeðL!LÞ  r s;SeðLÞ   min þ e rimm s;SeðLÞ  r s;SeðLÞ

Rate law expressions for chemical heterotrophic reduction of O2, NO3, SeO4, SeO3, and SeMet using Monod terms are written as (Bailey et al., 2013a):

where kj is the base rate constant for species j [T ]; Kj is the Monod half-saturation constant for species j [MfLf3]; IO2 and INO3 are the O2 and NO3 inhibition constants [MfLf3] signifying the species concentration at which lower-redox species can undergo appreciable rates of reduction; CO2;prod is the total mass of CO2 produced during organic matter decomposition and is used as an indicator of available OC for microbial consumption (Birkinshaw and Ewen, 2000); and E [–] is an environmental reduction factor that accounts for soil moisture h and soil temperature T and acts to temper the reaction rates based on microbial activity (Birkinshaw and Ewen, 2000; Bailey et al., 2013a). Reaction rates thus are tempered by the species concentration, the concentration of other reactant species, the presence of electron (e) donors, the presence of higher-redox species, and microbial activity. With the assumption that FeS2 is in unlimited supply in any adjacent shale material, rate law expressions for autotrophic reduction of O2 and NO3 can be written as:

ð4aÞ auto rauto f ;O2 ¼ kO2 C O2



C O2 K O2 þ C O2

ð4bÞ auto rauto f ;NO3 ¼ kNO3 C NO3

ð4cÞ

where rs represents the rate of all reactions that occur in the solid1 phase [MsL3 ]; PRt and PSt are the application rates of root and s T after-harvest stover mass, respectively; aRt,Se and aSt,Se the portions of the root and stover mass attributed to Se; dec signifies organic matter decomposition; and L, H, and M represent the litter, humus, and manure pool, with the arrow representing the direction of mass flow. Similar equations for N cycling, following the pattern established by Johnsson et al. (1987) and Birkinshaw and Ewen (2000), also are implemented but not shown. Mathematical expressions for organic matter decomposition, mineralization, and immobilization are implemented using first-order kinetics (van Veen and Paul, 1981), and are contained in Birkinshaw and Ewen (2000) for N and in Bailey et al. (2013a) for Se.





C NO3 K NO3 þ C NO3

ð6aÞ 

I O2 I O2 þ C O2

 ð6bÞ

The mass of Se released during these two reactions is dependent on the portion of autotrophically-reduced O2 or NO3 that contributes to the production of SeO4, and is dependent on the stoichiometry of Eq. (1) and the ratio of S to Se in the shale material: auto rauto f ;SeO4 ¼ r f ;O2

    1 1 Y Se:O2 þ r auto Y Se:NO3 f ;NO3 n n

ð7Þ

where Y Se:O2 is the mass of Se produced for mass of O2 consumed in Eq. (1a), Y Se:NO3 is the mass of Se produced for mass of NO3 consumed in Eq. (1b), and n is the ratio of S to Se in the shale material. This term n is included since autotrophic reduction of O2 and NO3 release both S and Se from shale. Referring to Eq. (1), Y Se:O2 is equal to 1.41 (315.84 g/224.0 g) and Y Se:NO3 is equal to 4.03 (789.6 g/196.0 g).

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

2.2.2. Applying the model to regional-scale agricultural systems Eqs. (2)–(4) were solved for the regional model domain using a three dimensional (3D) finite-difference grid. Input data requirements to solve the ADR equations include parameters to calculate values of N and Se cycling (fertilizer loading F, uptake U, organic matter decomposition and mineralization/immobilization); rate constants and other chemical reaction parameters for the rate law expressions of Eqs. (5)–(7); and the solute concentration of sources and sinks such as applied irrigation water, canal and river seepage, and groundwater pumping. For the first set of requirements, each vertical column of cells in the 3D grid is assigned a set of crop parameters according to the portions of cultivated fields residing within the spatial area of the grid cell. Crop parameters include the planting day, harvest day, plowing day, fertilizer loading and timing, maximum seasonal uptake values Seup and Nup [MLb2], maximum rooting depth drt,max [L], mass of roots and stover, and constants defining the root growth and daily uptake rate. Each crop type cultivated in the region is assigned a unique set of crop parameter values. Seasonal deficits of uptake are tracked and subtracted from the amount of Se and N placed back into the soil through root and stover material, to maintain mass balance (Wriedt and Rode, 2006). For perennial crops, incorporation of root and stover mass occurs only at the end of the perennial cycle. For cells encompassing portions of more than one field, model subroutines apply a weighting scheme that calculates composite crop parameter values, using the portion of each grid cell occupied by each crop type for each year of the simulation. Rate constants and other chemical reaction parameters are assigned either uniformly across the model domain or on a cell-by-cell basis. Each cell in the 3D grid also is assigned a material type to indicate the presence of FeS2 and FeSexS2x. During the model simulation, autotrophic reduction proceeds according to Eq. (6) within a given cell if any adjacent cell contains FeS2, with the mass of SeO4 released calculated according to Eq. (7). A reaeration term also is included that supplies O2 mass to the saturated groundwater zone via gaseous diffusion from the ambient atmosphere through the unsaturated zone (Bailey et al., 2013b). Since the model does not simulate solute transport in surface water, the concentration of O2, NH4, NO3, SeO4, and SeO3 in canal and river seepage water is provided by field data, with concentration values specified for specific dates during the model simulation and linear interpolation used to calculated daily concentration values between two sampling dates. The mass of each species brought into or removed from the subsurface via groundwater-surface water interaction is calculated for the cells designated in the River package of the MODFLOW model. Determining solute concentration for the infiltrating water at the ground surface of cultivated fields requires more complex accounting, due to the possibility of portions of multiple fields residing in a single grid cell, and also to infiltrating water based from canal-derived irrigation water, groundwater-pumped irrigation water, and rainfall. Composite concentration values for each grid cell are calculated by using the fraction of infiltrating water derived from each of the three sources and the solute concentrations for each water source. Solute concentration in canal-derived water is determined using the specified canal concentration values derived from field data. Solute concentration in pumped groundwater is determined by first, identifying the pump that supplies water to the cultivated field; and second, retrieving the current model-simulated groundwater solute concentrations from the grid cell associated with the screened depth of the pump. If the screened depth spans more than one grid cell, than an average concentration value can be calculated.

33

3. Model application to lower Arkansas River Valley, Colorado 3.1. Description of study region The study region (Fig. 2) is located within the Lower Arkansas River Valley (LARV) in southeastern Colorado. The boundary of the model domain is shown with a black line, and encompasses an area of 50,600-ha (125,000 acres), of which 26,400 ha (65,300 acres) are irrigated. The site is chosen because (i) it is designated as a seleniferous river basin (Seiler, 1997); (ii) all river segments of the LARV are designated as ‘‘water quality impaired’’ with respect to Se and placed on the current Clean Water Act 303(d) list for TMDL development (Colorado Water Quality Control Commission, 2012); (iii) it has been monitored for both hydrologic and chemical species components during the last decade, hence yielding an extensive dataset of measured groundwater and surface water solute concentrations; and (iv) a calibrated and tested groundwater flow model has been implemented for the region. The climate is semi-arid, with average monthly temperatures and month precipitation ranging from 1 °C and 0.7 cm during the winter months, respectively, to 25 °C and 5.0 cm during the summer months. The LARV, containing a total of about 109,000 ha (270,000 ac) of irrigated land, has been intensively irrigated for more than 100 years, and has served as one of the most productive agricultural areas for the state of Colorado. The Arkansas River (Fig. 2), originating from snowmelt in the upper Arkansas River Valley northwest of Pueblo, CO, has a flow pattern typical of snowmelt-derived rivers with high flows in the spring and low flows in the late summer and fall. The irrigation season commences mid to late-March and ends in early November, with un-lined irrigation canals receiving water from the Arkansas River during the period of March 15th to November 15th. Irrigation water in the study region is derived from one of six principal irrigation canals (Rocky Ford Highline, Catlin, Otero, Rocky Ford, Fort Lyon, and Holbrook) or groundwater pumps, as shown in Fig. 2. The designation of the study region as seleniferous is likely due to the preponderance of shale. The alluvial aquifer in the LARV, which ranges from approximately 4–34 m in thickness (Fig. 3A), is underlain by Cretaceous Shale (Pierre Shale, Niobrara Shale, Carlisle Shale, Graneros Shale) (Scott, 1968; Sharps, 1976) in both solid and weathered form. In the study region, the aquifer is underlain principally by Niobrara Shale, followed by Greenhorn Limestone, Carlile Shale, and Graneros Shale. Shale formations, predominantly Niobrara and Carlisle shale, residing within approximately 2 m of the ground surface also have been delineated (Fig. 2). Amount and distribution of cultivated crops in the study region is representative of the LARV. In terms of spatial extent, alfalfa is the dominant crop, followed by sorghum, corn, grass/pasture, wheat, melons, onion, oats, sunflower, and soybeans. Melons and onions, as the principal cash crops, receive the most irrigation water per unit cultivated area. The crop type cultivated for each field in 2006, provided by the Farm Service Agency in Rocky Ford, CO, is shown in Fig. 3B. Due to over-irrigation, canal seepage, and poor drainage resulting in salinization and waterlogging, the region has experienced a decrease in crop productivity (Burkhalter and Gates, 2005; Morway and Gates, 2012). The groundwater flow model of the region was constructed, calibrated, and tested as described in Morway et al. (2013) using the UZF1 package (Niswonger et al., 2006) of MODFLOW. Horizontal cell dimensions are 250 m by 250 m, resulting in grid cells with spatial areas comparable to dimensions of cultivated fields. The alluvial aquifer is divided into two layers, with a third layer constituting the bedrock. Applied irrigation water is based on crop type and weekly volumes of canal diversions from the Arkansas River, with timing of water application based on a hierarchy of crop

34

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Fig. 2. Location and hydrologic features of the study region in the Lower Arkansas River Valley in southeastern Colorado, showing the cultivated fields, irrigation canals, groundwater pumping wells, and extent of near-surface shale (within 2 m of ground surface).

Fig. 3. (A) Aquifer thickness (between ground surface and marine shale bedrock), (B) crop type for each cultivated field for the year 2006, (C) Calculated N fertilizer loading for 2006, and (D) calculated after-harvest stover mass plowed into the top soil layers for the year 2006. Seasonal fertilizer and stover loading are calculated using spatial distribution of crop type and stochastic perturbations based on realistic ranges of values.

demand. Pumping from groundwater wells (see Fig. 2) and canal seepage also are simulated. The model was tested against water table depth, groundwater discharge to the Arkansas River, and

several other variables. An example of model-calculated spatial distribution of water table elevation and soil water content is shown in Fig. 4.

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

35

Fig. 4. (A) Temporally-averaged water table elevation (m) and (B) Volumetric water content for June 2006 within the top soil layers, as simulated by the groundwater flow model of Morway et al. (2013).

Fig. 5. Location of groundwater observation wells (maroon targets) and surface water sampling points (red triangles) in the study region, from which water samples were taken and analyzed for species concentrations. Also shown is the division of the study region into canal command areas, with fields receiving irrigation water from the Highline, Otero, Catlin, Rocky Ford, Holbrook, and Fort Lyon canals. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3.2. Field data for model corroboration An extensive data set is available to test the model for the study region. Water samples have been collected and analyzed for concentrations of O2, NO3, SeO4, SeO3, and other salt ions since the growing season of 2006. Samples were taken routinely from 45 groundwater observation wells, 10 locations along the Arkansas River, and four locations in tributaries during nine sampling events between June 2006 and July 2009 (Gates et al., 2009). Samples from an additional 61 observation wells were taken for several of the sampling events. Samples were collected and analyzed according to accepted field and laboratory protocol, as described thoroughly in Gates et al. (2009). Fig. 5 shows the location of the observation wells and surface water sampling points within the study region. Model testing was conducted in a manner advocated by Konikow (2011), with the objective of reproducing major trends and spatio-temporal statistics rather than time series of concentrations for individual observation wells. With this objective, model results were tested against observed spatio-temporal averages of C SeO4 , C NO3 , and C O2 , the observed C NO3 —C SeO4 relationship (Bailey et al., 2012), the portion of dissolved Se that is accounted for by SeO4, and the global frequency distributions of C SeO4 and C NO3 . The spatio-temporal averages were calculated by command area (i.e., the collection of cultivate fields that receive irrigation water

from the same canal, see Fig. 5) and across multiple sampling events. Tables showing these values are presented in Supplementary Data. Spatial averaging by command area was performed due to the unique water rights priority, irrigation, and cultivation histories of each command area. As is common in regional modeling studies, the observation scale and model simulation scale were dissimilar. In comparing observed and simulated solute concentration values, it is recognized that point values from observation wells have uncertainty when compared with model results, i.e. the observed solute concentration averaged over a spatial area covered by a computational grid cell is not known. In an effort to evaluate the uncertainty in using point measurements as comparators for areal values, use was made of measurements of electrical conductivity (EC), as specific conductance at 25 °C, in groundwater monitoring wells at multiple points within each of eight cultivated fields in the study region. Between 7 and 13 wells were installed within each field, and approximately 12 sampling events occurred during each year. Three fields were sampled for two years, whereas five fields were sampled only for one year. For each sampling event, the average EC value and corresponding standard deviation were used to calculate the coefficient of variation (CV) for that sampling event, corresponding to the degree of variability of the EC of groundwater underyling the field. A more complete description of the analysis,

36

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

including data values for each sampling event for each field, is contained in Supplementary Data. Taking into account the CV for each field and for each year, the average CV of EC within a field is 0.42. Assuming that C SeO4 and C NO3 have a similar variability (statistically significant Pearson correlation coefficients of 0.38 and 0.29, respectively, with EC), this value was assigned to describe uncertainty in the spatio-temporal averages of observed solute concentration. In addition to using groundwater solute concentrations for model testing, Se and NO3 mass-balance calculations were performed for the surface water system (Arkansas River and its tributaries) to determine the approximate daily mass loadings of Se and N in the Arkansas River which are not accounted for by loading from the tributaries. It was assumed that a substantial portion of the unaccounted-for mass loading, especially during the non-irrigation season, can be attributed to mass loading from the aquifer to the Arkansas River, and thereby provides an additional test for model results. A complete description of the surface water mass balance calculations are provided in Supplementary Data. 3.3. Model set-up The period of model simulation was January 1 2006 through October 31 2009, with the period of January 1 2006 through March 31 2008 designated as the calibration period, and April 1 2008 through October 31 2009 designated as the testing period. Daily time steps were used. The finite-difference surface grid is the same as for the groundwater flow simulation, with grid cell dimensions of 250 m by 250 m in the horizontal directions. The vertical cell discretization was modified to use 7 grid layers to simulate localized chemical reactions and physical processes, particularly in the vicinity of the root zone (e.g., root growth and corresponding solute uptake, incorporation of dead root mass and stover mass). The top four layers are comprised in the first grid layer of the flow model, the next two layers are contained in the grid second layer, and the seventh layer, representing the bedrock shale, is the same as the third grid layer (see Fig. S-2 in Supplementary Data). The top two layers, comprising the root zone, each have a thickness of 0.5 m, followed by a third layer with a thickness of 1.0 m. The top three layers therefore correspond generally to the unsaturated zone, whereas layers 4–6 correspond to the saturated zone. Layer 4 corresponds to the depth at which groundwater samples typically were taken from the observation wells, therefore model results from layer 4 were compared with the observed groundwater concentrations during model calibration and model testing. To

preserve the groundwater flow field as established through the calibration procedure of Morway et al. (2013), volumetric flow rate and volumes of water sources and groundwater sinks were mapped from the three-layer grid to the seven-layer grid. A complete description of the mapping scheme is presented in Supplementary Data. Key crop parameters governing management, growth, solute uptake, and death are provided for each crop type in Table 1. Average values were supplied either by Michael Bartolo at the Arkansas Valley Research Center (AVRC) in Rock Ford, CO (Bartolo, personal communication) or by literature (e.g., Johnsson et al., 1987; Birkinshaw and Ewen, 2000). Typically, 40% of the annual fertilizer load is applied one week before the planting day, within the remaining 60% applied six weeks after planting. Values of selected parameters (planting day, harvest day, N fertilizer loading, root mass, stover mas, Se root content, Se stover content) were perturbed stochastically to establish uncertainty in parameter values and spatial heterogeneity in cultivation practices throughout the study region. CV values of 0.02, 0.01, 0.05, 0.10, 0.10, 0.05, and 0.05 are assigned to the average parameter values, respectively, with the CV values chosen to provide a reasonable range of parameter values. Resulting frequency distributions for planting day, N fertilizer loading, and stover mass are shown in Fig. S-9 in Supplementary Data. The spatial distribution of annual N fertilizer loading (kg/ha) and applied stover mass (kg/ha) in the model domain for the 2006 growing season are shown in Fig. 3C and D. Notice the high values of N fertilizer loading and stover mass in areas of corn cultivation (compare to Fig. 3B). A list of parameters for the chemical reactions and associated values is presented in Table 2, with parameters grouped according to the species they affect: organic matter decomposition, O2, N, and Se. Parameter values were derived from the published literature (e.g., Johnsson et al., 1987; Tokunaga et al., 1997; Guo et al., 2000; Lee et al., 2006) and are used in pre-calibration simulations as discussed in Section 3.4. In pre-calibration simulations, values are treated as spatially-constant. However, they were varied spatially during model calibration according to canal command area. The measured daily average air temperature at the AVRC climatic station, using in the model to calculate soil temperature for moderating chemical reaction rates, is shown in Supplementary Data. Also shown in Supplementary Data is the percentage of maximum microbial activity as a function of percent soil saturation for mineralization, denitrification, and general chemical reactions. For simulation of autotrophic reduction and release of SeO4 in the

Table 1 Agricultural management and crop parameter values for the model application to the study region in the Arkansas River Valley in southeastern Colorado. Crop type Units

Planting day –

Harvest day –

Plow day –

PSt kg ha1

drt,max m

CNRT –

F NH4 kg ha1

Nup kg ha1

aSt,Se

aRt,Se

–

–

Seup g ha1

Alfalfa Bean Corn Melon Onion Pasture Pumpkin Sorghum Spring Grain Squash Sunflower Vegetable Winter Wheat

30-April 20-May 1-May 15-May 20-March 30-August 1-June 20-May 1-April 20-May 1-June 25-April 30-September

30-September 30-September 25-October 10-August 15-September 30-September 30-September 15-October 15-July 25-July 10-October 30-August 5-July

20-October 20-October 14-November 30-August 5-October 20-October 20-October 4-November 4-August 14-August 30-October 19-September 25-July

561.6 561.6 5616 561.6 561.6 0 561.6 1684.8 1684.8 561.6 561.6 561.6 1684.8

1.83 0.91 1.22 1.22 0.46 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.91

25 25 70 25 25 70 25 70 70 25 25 25 70

22.4 140 252 112 140 140 140 112 112 140 140 140 112

22.4 84.2 224.6 112.3 78.6 112.3 84.2 112.3 112.3 84.2 84.2 84.2 112.3

0.0000 0.0000 0.0010 0.0000 0.0000 0.0030 0.0000 0.0010 0.0010 0.0000 0.0000 0.0000 0.0010

0.0002 0.0001 0.0022 0.0001 0.0001 0.0053 0.0001 0.0022 0.0022 0.0001 0.0001 0.0001 0.0022

1.3 0.3 10.8 1.3 2.7 10.7 2.7 10.4 8.7 2 0.1 2.7 8.7

dpw (depth of plowing) is 1.0 m for all crops except beans (0.8 m). PRt (seasonal mass of root mass) is 500 kg ha1 for all crop types. CNST (carbon:nitrogen ratio in stover mass) is 50 for all crop types. kup,Se (rate constant of Se uptake) is 0.01 d1 for all crop types. cup (SeO4–SeO3 crop uptake ratio) is 10 for all crop types.

37

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46 Table 2 Parameters for chemical reactions involving organic matter decomposition, Dissolved Oxygen, Nitrogen species, and Selenium species. Org. Matter Decomp.

Dissolved oxygen

Nitrogen

Selenium

Param.

Value

Unit

Param.

Value

Unit

Param.

Value

Unit

Param.

Value

Unit

kL

0.25

d1

khet O2

2.0

g mf3

HC/N

12.0

–

HC/Se

1.75  105

–

kH

0.003

d1

kauto O2

0.1

g mf3

BC/N

8.0

–

BC/Se

1.23  105

–

fe fh

0.5 0.2

– –

K O2

1.0

g mf3

IO2 knit

1.0 0.8

g mf3 d1

INO3

0.50 0.02

g mf3 d1

K CO2

0.75

g mf3

kvol

0.1

d1

0.02

d1

0.1

d

1

0.02

d1

d

1

0.02

d1

3000 0.1

– –

khet NO3 kauto NO3

0.01

K NO3 K d;NH4

presence of shale, cells intersecting with the surface shale (see Fig. 2) and with the bedrock layer of shale in the model were designated as containing FeS2 and FeSexS2x. During field sampling it was recognized that a yellow clay material located in certain portions of the aquifer also contained residual Se. Grid cells intersecting the spatial distribution of yellow clay hence also were designated as S and Se sources. The ratio of S mass to Se mass contained in the shale material (n in Eq. (7)) was given a value of 3000, representative of numerous published S:Se ratios in marine shale material (see Fig. S-7 in Supplementary Data). For canal solute concentrations, which were used for canal seepage and canal-derived irrigation water solute mass loadings to the system, the species concentration at the canal diversion point on the Arkansas River were used, with the species concentration assumed to be constant along the length of the canal. The species concentrations for each canal for each of the 9 sampling events are presented in Supplementary Data. For Arkansas River solute concentrations, data are available for each of the 9 sampling events at 9 sampling points along the reach of the river within the model domain. For most Arkansas River tributaries (Patterson Hollow, Timpas Creek, Crooked Arroyo, Anderson, and Horse Creek, see Fig. 2), one measurement of species concentration is available for each of the 9 sampling events. 3.4. Estimation of model parameter values Estimation of model parameters involved analyzing solute concentrations and mass loadings during a 20-year spin-up simulation and during the 2006–2008 calibration period. The spin-up simulation was included to achieve steady fluctuations of solute concentrations and to provide initial conditions for the calibration period simulation. The parameter estimation methodology includes both manual calibration and automated calibration using the PEST

10.0 3.5

g mf –

3

khet SeO4 hetðSe Þ kSeO3 s hetðSeMetÞ kSeO3 khet SeMet n K d;SeO3

(Parameter ESTimation) software (Doherty, 2007). An iterative approach using both the spin-up simulation and the 2006–2008 simulation was used in the model calibration procedure due to the dependence of initial concentrations on the parameters used in the spin-up simulation. The overall procedure is represented in Fig. 6, with the following elements: (i) Establish an initial set of initial conditions for the 2006– 2008 simulation using a 20-year spin-up simulation (ii) Use global sensitivity analysis for the 2006–2008 simulation to identify influential parameters (iii) Use PEST to provide an estimate of influential parameter values, with parameter values assigned and groundwater concentrations averaged according to command area (iv) Re-run the 20-year spin-up with the new parameter values to establish new set of initial conditions, with comparisons made between Se and NO3 mass loading to the Arkansas River. The spin-up simulation was re-run using parameter values until the mass loadings at the end of the spin-up simulation matched well with the range of observed daily mass loadings during the 2006–2008 time period. For (i), the flow and cropping patterns for the year 2006 were repeated for 20 years using the crop and chemical reaction parameters in Tables 1 and 2. Initial concentrations for each grid layer for the spin-up simulation are summarized in Supplementary Data. Using the set of concentrations from the end of the spin-up simulation as initial conditions, a sensitivity analysis was conducted on the 2006–2008 simulation using the FAST (Fourier Amplitude Sensitivity Test) (Cukier et al., 1973) method to determine the relative influence of selected parameters and system sources/sinks on model results. Selected source/sink terms include NH4 fertilizer, seasonal N uptake, seasonal Se uptake, and C NO3 and C SeO4 in

Fig. 6. Framework for model simulation initialization and performing estimation of model parameter values.

38

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

surface water. Selected chemical reaction parameters include the rates of litter pool decomposition, kL; humus pool decomposition, kH; autotrophic reduction of O2 in the presence of shale, kauto O2 ; autotrophic reduction of NO3 in the presence of shale, kauto NO3 ; nitrification, knit; NH4 volatilization, kvol; and heterotrophic reduction of het NO3, khet NO3 , and of SeO4, kSeO4 . A total of 1053 simulations were run to analyze the influence of all 13 parameters and source/sink terms. For each simulation the parameter values were sampled from the parameter space based on the probability distributions, with the average values contained in Tables 1 and 2. All parameters were given a CV of 0.20 except for the surface water concentration of NO3 and SeO4, which was assigned a CV of 0.30. CV values were selected by comparing the resulting spread of parameter values to values found in the literature and from field data in the LARV. Chemical reaction parameters were perturbed according to a lognormal distribution, whereas all others were perturbed according to a normal distribution. Resulting C NO3 and C SeO4 in groundwater from each of the 1053 simulations were spatially-averaged globally (across the entire model domain), according to command area, and according to crop type. The latter two were included to determine if parameter influences varied according to irrigation and cultivation patterns. Limited results are shown here. Results are presented using global sensitivity plots (Saltelli et al., 2008), with the first-order sensitivity index (the main effect of the parameter on the change in model output), Si, along the abscissa and the interaction effect (effect of the parameter on the change in model output due to interaction with other parameters), Sti, along the ordinate. Fig. 7 shows the results for C NO3 and C SeO4 in groundwater averaged across the entire model domain for 2006 (Fig. 7A and C), with knit, C NO3 of surface water, and kauto NO3 governing simulated C NO3 in groundwater and kauto dominating the resulting simulated C SeO4 . The influence of O2 parameters, however, varies locally within the model domain. auto For example, kauto NO3 has a stronger influence than kO2 on both C NO3 and C SeO4 in groundwater within the Catlin canal command area (Fig. 7B and D). auto Values of knit, kauto O2 , and kNO3 were assigned to each command area and estimated using PEST, with the command area average values of C NO3 and C SeO4 in groundwater used as targets. The objective function was defined as a minimization of the sum of squared deviations between simulated and observed values. Using parameter values derived from the PEST simulation, the 20-year spin-up

simulation was re-run to allow these new parameter values to impact the initial conditions of the 2006–2008 simulation. It was recognized that kauto and khet O2 SeO4 strongly influences the Se loading to the Arkansas River, hence values were modified during step (iv). Several iterations were required to achieve acceptable simulated values of Se mass loading due to the interplay between kauto O2 (releasing SeO4 into the groundwater from shale, inhibiting chemical reduction of SeO4 in groundwater, and influencing the resulting mass loading to surface water) and khet SeO4 (chemically reducing SeO4 in the riparian zone). Initial parameter values for each command area, as well as parameter values estimated after the completion of steps (iii) and (iv), are shown in Table 3. Based on the range of published values of first-order denitrification rates (see Fig. S-8 in Supplementary Data), the values reported in Table 3 are reasonable and physically realistic.

Table 3 auto het Initial (Pre-PEST), PEST-derived, and final values of kauto O2 , kNO3 , knit , and kSeO4 for each command area within the study region. The final values were determined by modifying the PEST-derived values during re-runs of the 20-year spin-up simulation. The unit for all values is d1. Command area

kauto O2

RF Highline

0.10

0.703

0.450

Catlin Rocky Ford Fort Lyon Holbrook Outside

0.10 0.10 0.10 0.10 0.10

3.00 0.00012 0.0001 0.108 0.492

3.00 0.00012 0.0001 0.040 0.010

kauto NO3

a

Pre-PEST

PEST

Finala

Parameter

RF Highline

0.01

0.0003

0.0003

Catlin Rocky Ford Fort Lyon Holbrook Outside

0.01 0.01 0.01 0.01 0.01

1.00 0.0002 0.0003 0.303 0.01

1.00 0.0002 0.0003 0.303 0.01

knit

RF Highline Catlin Rocky Ford Fort Lyon Holbrook Outside

0.80 0.80 0.80 0.80 0.80 0.80

0.80 0.80 0.80 0.80 0.80 0.80

1.50 2.00 1.00 0.80 0.20 0.40

khet SeO4

–

0.06

0.06000

0.02000

Based on additional runs of the 20-year spin-up simulation.

Fig. 7. Global sensitivity plots for (A) model domain-averaged groundwater C NO3 , (B) groundwater C NO3 within the Catlin canal command area, (C) model domain-averaged groundwater C SeO4 , and (D) groundwater C SeO4 within the Catlin canal command area.

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

4. Results and discussion 4.1. Spin-up simulation Results from the final spin-up simulation are shown in Figs. 8 and 9. Fig. 8 shows the times series of daily calculated solute concentration during the spin-up simulation for C SeO4 (Fig. 8A) and C NO3 (Fig. 8C) for selected cells in layer 4 that underlie various crop types. A times series of the simulated daily mass loading (kg) of NO3 and Se to the Arkansas River is presented in Supplementary Data (Fig. S-9). Both Figs. 8 and S-9 demonstrate the dynamic equilibrium reached after 20 years of simulation. Fig. 9 shows contour plots of spin-up simulation results. The contour plots of species concentration correspond to initial concentrations used for the 2006–2008 simulation. Fig. 9A shows C SeO4 , with high concentration values corresponding typically to areas adjacent to near-surface shale. Fig. 9B shows the cumulative mass of chemically reduced SeO4, demonstrating the correspondence of SeO4 reduction to riparian zones. The distribution of C SeO3 , the product of SeO4 reduction, is shown in Fig. 9C. Fig. 9D shows the spatial distribution of C NO3 . For C O2 (Fig. 9F), a plot of water table depth for the same simulation time also is shown (Fig. 9E) to demonstrate the influence of the reaeration term. Areas of shallow water tables typically have high values of C O2 , and vice versa, since areas with deep and shallow water tables are supplied via the reaeration term with small and large amounts of O2 from the ground surface, respectively. Areas of high water table elevation also typically correspond to areas of cultivation and associated application of irrigation water, with O2 mass also supplied in the irrigation water. 4.2. Calibration and testing periods 4.2.1. Groundwater solute concentrations An example of cell-by-cell results for solute concentration is presented in Fig. 8. Fig. 10B shows the daily-calculated C SeO4 during 2006–2009 for the same cells as presented for the spin-up simulation in Fig. 8A, with initial conditions in 2006 established by the final concentrations of the spin-up simulation. Similar results are presented for C NO3 in Fig. 8D. Fig. 8B and D demonstrate the

39

day-to-day and year-to-year fluctuation of concentrations due to variations in fertilizer loading, applied irrigation water, rainfall, crop uptake, etc. Observed and simulated command area averages of C SeO4 , C NO3 , and C O2 in groundwater are shown in Fig. 10 for both the calibration and testing periods. Results for the Rocky Ford canal command area as well as for the Outside area should be treated with caution due to the lower number of samples (8 and 10 samples, respectively, available for the calibration period). The large discrepancy between the simulated (93.8 lg L1) and observed (10.9 lg L1) C SeO4 values for the Rocky Ford canal command area (Fig. 10A) could be an artifact of not having samples in areas of high C SeO4 . In general, for command areas with a large number of samples (see Supplementary Data), the simulated values of C SeO4 for the calibration period compare favorably with the observed values (Fig. 10A and B), particularly if the uncertainty in the observed values are taken into account. Excluding the Rocky Ford canal values, and treating the reported observed value as the true value and hence neglecting uncertainty, the sum of the absolute differences between observed and simulated is 92.6 lg L1 and 125.4 lg L1 for the calibration and testing periods, respectively. As shown in Fig. 10C and D, favorable matches occur for C NO3 for all command areas except for the Catlin command area during the Calibration period (11.9 mg L1 observed vs. 2.2 mg simulated mg L1). The high observed average, however, is a result of two groundwater samples with measured C NO3 values of 66.0 mg L1 and 45.6 mg L1, measured at Well 12 (see Fig. 4) on October 6 2007 and Well 41 on May 23 2007, respectively. These concentrations constitute two of the four highest measured values of C NO3 during the Calibration period throughout the region, and also are much higher than other measured concentrations at these same wells during the calibration period (10.3, 11.8, and 26.0 mg L1 for Well 12 and 14.4 mg L1 for Well 41). Hence, the high average of C NO3 for the Catlin command area is a result of unusually high single-occurrence measurement concentrations that were not captured by the model. Overall, the sum of the absolute differences between observed and simulated for C NO3 is 17.9 mg L1 and 10.0 mg L1 for the Calibration and Testing periods, respectively. Favorable matches also occur for C O2 (Fig. 10E and F) except for the Holbrook command area where the simulated value

Fig. 8. Fluctuation of C SeO4 in groundwater for selected individual grid cells underyling various cultivated crops for the (A) spin-up and (B) 2006–2009 simulations. Fluctuation of C NO3 in groundwater for selected individual grid cells underyling various cultivated crops for the (C) spin-up and (D) 2006–2009 simulations.

40

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Fig. 9. Contour plots of system response variables at the end of the 20-year spin-up simulation for (A) C SeO4 , with regions of near-surface shale depicted, (B) total mass of reduced SeO4 during the 20-year period, (C) C SeO3 , (D) C NO3 , (E) water table depth, and (F) C O2 .

(6.7 mg L1) is much higher than the observed value (2.5 mg L1) during the calibration period. A similar discrepancy occurs for the testing period (5.7 mg L1 compared to 2.0 mg L1). As seen in Fig. 9E and F, the shallow water tables in the Holbrook command area result in high C O2 in the saturated zone. Hence, the reaeration term likely is overestimating the amount of O2 supplied to the saturated zone in this area. Overall, the sum of the absolute differences between observed and simulated for C O2 is 9.9 mg L1 and 8.8 mg L1 for the calibration and testing periods, respectively. In regards to the high simulated C SeO4 values in the Rocky Ford Highline canal command area during the testing period (Fig. 10B), the average C SeO4 in the groundwater did not decrease during the testing period due to the presence of high-concentration areas in the vicinity of surface shale. Although several observation wells in the Rocky Ford Highline command area are located near surface shale, the average observed C SeO4 likely would be higher if more samples were collected near shale. For example, Well 20 is located near shale and had an average of 266.8 lg L1 over the 9 sampling events with a maximum concentration of 1350 lg L1 measured for a sample collected on June 19, 2006. In addition, groundwater sampled from a piezometer well in the Rocky Ford Highline canal command area (Bailey et al., 2012) had an average C SeO4 value of 2508 lg L1 for seven samples taken from August 2009 to December 2010 (values of 1620, 1520, 2870, 2650, 2860, 3300, and 2740 lg L1). The piezometer well had a screen located within

yellow clay. High model-simulated averages of C SeO4 as compared to observed values is not a definite proof that the model yields inaccurate results, i.e., a higher average value of C SeO4 in the Rocky Ford Highline canal command area may exist in the actual aquifer system and would be observed if additional observation wells were installed. Time-series comparisons of observed and simulated C SeO4 and C NO3 for several observation wells are shown in Fig. 11, for a ‘‘good’’ match (Well 204, Fig. 11A), a ‘‘poor’’ match (Well 65, Fig. 11B), and a ‘‘satisfactory’’ match (Well 93, Fig. 11C). Although not a target for calibration, these point-by-point comparisons are presented to demonstrate that although poor matches can result due to inadequate representation of actual system stresses (Fig. 11B), in general the model is able to produce cell-by-cell results that are similar in magnitude and range to observed point values. 4.2.2. Mass loading to the Arkansas River Simulated mass loadings of Se and NO3 and observed unaccounted-for mass loadings during the 2006–2009 period are shown in Fig. 12. Average spatial distribution of loadings (kg) are shown along the Arkansas River for both Se (Fig. 12B) and NO3 (Fig. 12D) to indicate the areas of high loading to the river. For Se loadings (Fig. 12A), the simulated daily mass loadings are similar to observed values in regards to fluctuation pattern and general range of values. As indicated in Section 3.2, the observed loadings

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

41

Fig. 10. Comparison between simulated and observed average C SeO4 , C NO3 , and C O2 for each command area, for the Calibration period (A, C, E) and for the testing period (B, D, F). The whisker bars associated with the observed values represent the average values ± one standard deviation of an assumed normal distribution with CV estimated as described in Section 3.2.

provide a measure of all mass not accounted for by loading from the gauged and sampled tributaries. This includes mass loaded to the river by ungauged tributaries and drains and by direct surface runoff, which is not simulated in the model. As these surface flows occur mainly during the irrigation season, this may account for the smaller discrepancy in Fig. 12A and C during the non-irrigation season. This is especially true for NO3, due to the high concentration of NO3 in the surface runoff water due to the fertilizer loading, and likely accounts for the larger discrepancy in the NO3 groundwater mass loading (Fig. 12C) than for the Se mass loading.

4.2.3. Additional model corroboration Further comparison between simulated and observed values is performed to assess model corroboration. Fig. 13A and B shows log–log plots of C NO3 vs. C SeO4 for observed and simulated values for both the Calibration and Testing periods. For observed values, 141 data values are available for the four sampling events of the Calibration period and 199 data values are available for the five sampling events of the Testing period. For simulated results, values of C NO3 and C SeO4 were taken from each grid cell in layer 4 for the simulation times corresponding to the dates of field sampling. This resulted in 30,374 and 38,140 simulated values for the Calibration and Testing periods, respectively. Also shown in Fig. 13A and B are the C NO3 and C SeO4 detection limits for the groundwater samples analyzed in the laboratory (Gates et al., 2009). Data values shown at these limits likely have actual concentrations lower than the limit values, and hence would correspond more closely with simulated values.

As seen in Fig. 13A and B, the simulated values capture the main positive correlation between C NO3 and C SeO4 . For the Calibration period, the correlation coefficient for the observed and simulated values is 0.56 and 0.36, respectively, and for the Testing period the correlation coefficient for the observed and simulated values is 0.71 and 0.42. For both the Calibration and Testing periods, however, it can be seen in Fig. 13A and B that simulated C SeO4 values can be high (>102 lg L1) in grid cells where C NO3 values are close to 0.0 mg L1 (<101 mg L1), contrary to the expected positive relationship. This also occurs for the observed data and is seen particularly for the Testing period results in Fig. 13B, wherein many data points have values of C NO3 at the detection limit, indicating that high values of C SeO4 occur in groundwater samples that have low C NO3 values. These high C SeO4 values likely are due to the presence of O2 in the groundwater, which similar to NO3 is able to oxidize residual Se from shale as well as inhibit SeO4 reduction. This possibility is explored using a log–log plot of C NO3 vs. C O2 (Fig. 13C and D) for the calibration and testing periods. As seen in the figure, high values of C O2 (>100.5 mg L1) occur in some grid cells where C NO3 values are close to 0.0 mg L1 (<101 mg L1), and hence may be the cause of high C SeO4 values in grid cells where C NO3 values are close to 0.0 mg L1. The second global measure is the portion of dissolved Se mass in groundwater accounted for by SeO4, with remaining mass attributed to SeO3. This measure determines if the model is able to adequately capture the observed relationship between SeO4 and SeO3, with the latter species a product of SeO4 chemical reduction. The measure thus also analyzes the ability of the model to simulate accurately the occurrence and rate of SeO4 reduction. For analysis,

42

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Fig. 11. Comparison of simulated and observed C NO3 and C SeO4 at (A) Well 204, (B) Well 65, and (C) Well 93, demonstrating a variety in model results when comparing with point measurements. Simulated values are taken from grid cells in layer 4 of the model. Error bars are included to indicate the variability of the observed value for a representative cultivated field.

Fig. 12. Comparison of observed unaccounted-for daily mass loading and simulated mass loading from the aquifer to the Arkansas River for (A) Se and (C) NO3. Shaded areas represent the non-irrigation season (November 1 to March 15). The spatial distribution of temporally-averaged simulated mass loading along the length of the Arkansas River within the study region also is shown for (B) Se and (D) NO3.

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

43

Fig. 13. Log–log plots of C NO3 vs. C SeO4 for observed (red dots) and simulated (blue dots) values (from layer 4) for both the (A) Calibration and (B) Testing periods, and log–log plots of C NO3 vs. C O2 for observed and simulated values for the (C) calibration and (D) testing periods. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Frequency distributions of observed and simulated values of C SeO4 for the (A) calibration and (B) testing periods, and fitted Pearson 5 distributions for the (C) calibration and (D) testing periods.

the portion of dissolved Se mass accounted for by SeO4 was calculated for both the observed and simulated values, with the simulated values of C SeO4 and C SeO3 taken from each grid cell in layer 4 of the model. The frequency distributions of these values for both the Calibration and Testing periods are presented in Supplementary Data (Fig. S-10). A close match occurs between the frequency distributions for the observed and simulated values, with the vast

majority of the observed and simulated values having very high (>0.99) SeO4 portions. The third global comparison is the frequency distribution of C SeO4 , C NO3 , and C O2 , with distributions constructed using every observed value in the region and every simulated value from layer 4 of the model. Whereas previous comparisons between observed and simulated values were performed according to the mean of the

44

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Fig. 15. Frequency distributions of observed and simulated values of C NO3 for the (A) calibration and (B) testing periods, and fitted Pearson 5 distributions for the (C) calibration and (D) testing periods.

concentration values (see Section 4.2.1), this comparison analyzes the ability of the model to represent accurately the spread of the concentration values as compared to observed values and the relative frequency of concentration values within a given concentration interval, i.e., to verify that the model captures the regional spatio-temporal statistics of solute concentration. Fig. 14 shows the frequency distributions and the fitted statistical distributions for C SeO4 for both the calibration and testing periods, and Fig. 15 shows the same for C NO3 . For statistical fitting, the Pearson 5 distribution was found to best fit the frequency distributions, and parameter values for the Pearson 5 distribution (a, b, shift) are shown in Figs. 14C,D and 15C,D. As can be seen from Figs. 14 and 15, excellent matches occur between the distributions of the observed and simulated values, particularly for the Testing period, indicating that the model produces the spread of C SeO4 and C NO3 values as observed in the aquifer system. 5. Summary and conclusions The UZF-RT3D model, with its accompanying Selenium (Se) and Nitrogen (N) cycling and chemical reaction modules, was applied to a regional-scale study site (50,600 ha) along the Arkansas River Valley in southeastern Colorado to determine its worth in assessing large-scale systems. The reaction modules include all relevant sources and sinks (fertilizer, crop uptake, species mass in irrigation water, canal seepage, crops roots and stover) for agricultural areas and all pertinent chemical reactions (organic matter decomposition, mineralization/immobilization, nitrification, volatilization, heterotrophic chemical reduction, autotrophic chemical reduction) for simulating the reactive transport of Se and N species in agricultural soil and groundwater systems. The model was used to simulate the fate and transport of Se and N species in the soil and groundwater systems during the 2006–2009 time period, with a 20-year spin-up simulation used to establish initial conditions. Model results are tested against spatio-temporal averages of selenate (SeO4), nitrate (NO3), and dissolved oxygen (O2) groundwater concentrations, daily mass loadings of Se and NO3 from the aquifer to the Arkansas River, inter-species relationships in the groundwater, and overall global

statistics of solute groundwater concentration. Results indicate that the model is successful in replicating the major spatial and temporal patterns in Se and NO3 contamination and transport. Of particular mention is the ability of the model to reproduce to a satisfactory degree the positive relationship between NO3 and SeO4 in groundwater. Of prime importance is the inclusion of autotrophic chemical reduction of O2 and NO3 in the presence of sulfur and Se-bearing marine shale, a dominant geologic formation throughout the Arkansas River Valley. Whereas fertilizer loading governs the presence and contamination of NO3 in the groundwater system, model results, corroborated with field data, indicate that the location of shale dictates to a large extent the level of Se contamination in the aquifer and the resultant mass loadings of Se to the Arkansas River. Whereas this paper presents the calibration and testing of the model for the study region, ongoing and future work include an assessment of land- and water-management strategies on the remediation of Se groundwater concentrations and mass loadings to the Arkansas River. Furthermore, as marine shale is present in numerous aquifer systems in the western United States and worldwide, and since typically these aquifer systems underlie cultivated, irrigated landscapes, many opportunities exist to employ the UZF-RT3D model for Se assessment and remediation investigation. Acknowledgements The interested cooperation of more than 120 landowners and growers in the Colorado’s Lower Arkansas River Valley made this study possible and is much appreciated. We appreciate Dr. Michael Bartolo for providing ranges of values for parameter and model inputs. Financial support was provided by grants from the Colorado Agricultural Experiment Station and the Nonpoint Source Program of the Colorado Department of Public Health and Environment, the Colorado Agricultural Experiment Station, the United States Bureau of Reclamation, the United States Department of Agriculture (USDA), the United States Geological Survey (USGS), the Southeastern Colorado Water Conservancy District, the Lower Arkansas Valley Water Conservancy District, the Colorado Water

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Institute, the Colorado Water Conservation Board, and the Colorado Division of Water Resources (CDWR). Valuable cooperative assistance was provided by the USDA Natural Resources Conservation Service, the District 2 Office of the CDWR, the Pueblo Subdistrict Office of the USGS, and the USDA Farm Services Agency. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the opinions or policies of the U.S. government. Mention of trade names or commercial products does not constitute their endorsement by the U.S. government. We also appreciate the thorough and helpful comments by the reviewers. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.jhydrol.2014. 04.039. References Afzal, S., Younas, M., Ali, K., 2000. Selenium speciation sutides from Soan-Sakesar Valley, Salt Range, Pakistan. Water Int. 25, 425–436. Ahlrichs, J.S., Hossner, L.R., 1987. Selenate and selenite mobility in overburden by saturated flow. J. Environ. Qual. 16, 95–98. Alemi, M.H., Goldhamer, D.A., Nielsen, D.R., 1988. Selenate transport in steady-state, water-saturated soil columns. J. Environ. Qual. 17 (4), 608–613. Alemi, M.H., Goldhamer, D.A., Nielsen, D.R., 1991. Selenate transport in steady-state, unsaturated soil columns. J. Environ. Qual. 20, 89–95. Alfthan, G., Wang, D., Aro, A., Soveri, J., 1995. The geochemistry of selenium in groundwaters in Finland. Sci. Total Environ. 162, 93–103. Aro, A., Alfthan, G., Ekholm, P., Varo, P., 1998. Effects of selenium supplementation of fertilizers on human nutrition and selenium status. In: Frankenberger, W.T., Engberg, R.A. (Eds.), Environmental Chemistry of Selenium. Marcel Dekker, New York, New York, USA, pp. 81–98. Bailey, R.T., Hunter, W.J., Gates, T.K., 2012. The influence of nitrate on selenium in irrigated agricultural groundwater systems. J. Env. Qual. 41 (3), 783–792. Bailey, R.T., Gates, T.K., Halvorson, A.D., 2013a. Simulating variably-saturated reactive transport of selenium and nitrogen in agricultural groundwater systems. J. Contam. Hydrol. 149, 27–45. Bailey, R.T., Morway, E.D., Niswonger, R., Gates, T.K., 2013b. Modeling variably saturated multispecies reactive groundwater solute transport with MODFLOWUZF and RT3D. Groundwater 51 (5), 752–761. Balistrieri, L.S., Chao, T.T., 1987. Selenium adsorption by Goethite. Soil Sci. Soc. Am. J. 51, 1145–1151. Birkinshaw, S.J., Ewen, J., 2000. Nitrogen transformation component for SHETRAN catchment nitrate transport modelling. J. Hydrol. 230, 1–17. Burkhalter, J.P., Gates, T.K., 2005. Agroecological impacts from salinization and waterlogging in an irrigated river valley. J. Irrig. Drain. Eng. 131 (2), 197–209. Bye, R., Lund, W., 1982. Determination of selenium in pyrite by the atomic absorption-hydride generate technique. Fresenius Z. Anal. Chem. 313, 211–212. Calderone, S.J., Frankenberger, W.T., Parker, D.R., Karlson, U., 1990. Influence of temperature and organic amendments on the mobilization of selenium in sediments. Soil Biol. Biochem. 22, 615–620. Colorado Water Quality Control Commission, 2012. Colorado’s Section 303(D) List of Impaired Waters and Monitoring and Evaluation List, Regulation No. 93, Colorado Dept. Public Health and Environ., Denver CO. Combs Jr., G.F., Su, Q., Liu, C.H., Combs, S.B., 1986. Effect of dietary selenite, copper and zinc on tissue trace mineral levels in chicks. Biol. Trace Elem. Res. 11, 51– 64. Cukier, R., Fortuin, C.M., Schuler, K.E., Petschek, A.G., Schaibly, J.H., 1973. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. i theory. J. Chem. Phys. 59, 3873–3878. Doherty, J., 2007. PEST User Manual, fifth ed. Watermark Numerical Computing, Brisbane, Queensland, Australia. Ellis, D.R., Salt, D.E., 2003. Plants, selenium and human health. Curr. Opin. Plant Biol. 6, 273–279. Fernández-Martínez, A., Charlet, L., 2009. Selenium environmental cycling and bioavailability: a structural chemist point of view. Rev. Environ. Sci. Bio/ Technol. 8, 81–110. Fio, J.L., Fujii, R., Deverel, S.J., 1991. Selenium mobility and distribution in irrigated and nonirrigated alluvial soils. Soil Sci. Soc. Am. J. 55, 1313–1320. Flury, M., Frankenberger, W.T., Jury, W.A., 1997. Long-term depletion of selenium from Kesterson dewatered sediments. Sci. Total Environ. 198, 259–270. Frankenberger, W.T., Arshad, M., 2001. Bioremediation of selenium-contaminated sediments and water. BioFactors 14, 241–254. Gates, T.K., Cody, B.M., Donnelly, J.P., Herting, A.W., Bailey, R.T., Mueller-Price, J., 2009. Assessing selenium contamination in the irrigated stream-aquifer system of the Arkansas River, Colorado. J. Environ. Qual. 38, 2344–2356.

45

Gerla, P.J., Sharif, M.U., Korom, S.F., 2011. Geochemical process controlling the spatial distribution of selenium in soil and water, west central South Dakota. Environ. Earth Sci. 62, 1551–1560. Guo, L., Frankenberger, W.T., Jury, W.A., 1999. Evaluation of simultaneous reduction and transport of selenium in saturated soil columns. Water Resour. Res. 35, 663–669. Guo, L., Jury, W.A., Frankenberger, W.T., Zhang, Y., 2000. Characterizing kinetics of sequential selenium transformation in soil. J. Environ. Qual. 29, 1041–1048. Hamilton, S.J., 1998. Selenium effects on endangered fish in the Colorado River basin. In: Frankenberger, W.T., Engberg, R.A. (Eds.), Environmental Chemistry of Selenium. Marcel Dekker, New York, New York, USA, pp. 297–314. Hudak, P.F., 2010. Nitrate, arsenic and selenium concentrations in The Pecos Valley Aquifer, West Texas, USA. Int. J. Environ. Res. 4, 229–236. Johnsson, H., Bergström, L., Jansson, P., Paustian, K., 1987. Simulated nitrogen dynamics and losses in a layered agricultural soil. Agric. Ecosyst. Environ. 18, 333–356. Kishchak, I.T., 1998. Supplementation of selenium in the diets of domestic animals. In: Frankenberger, W.T., Engberg, R.A. (Eds.), Environmental Chemistry of Selenium. Marcel Dekker, New York, New York. Konikow, L.F., 2011. The secret to successful solute-transport modeling. Ground Water 49 (2), 144–159. Kulp, T.R., Pratt, L.M., 2004. Speciation and weathering of selenium in Upper Cretaceous chalk and shale from South Dakota and Wyoming, USA. Geochim. Cosmochim. Acta 68 (18), 3687–3701. Lee, M.-S., Lee, K.-K., Hyun, Y., Clement, T.P., Hamilton, D., 2006. Nitrogen transformation and transport modeling in groundwater aquifers. Ecol. Mod. 192, 143–159. Liu, C.W., Narasimhan, T.N., 1994. Modeling of selenium transport at the Kesterson reservoir, California, USA. J. Cont. Hydrol. 15, 345–366. Masscheleyn, P.H., Delaune, R.D., Patrick, J.W.H., 1989. Transformations of selenium as affected by sediment oxidation–reduction potential and pH. Environ. Sci. Technol. 24, 91–96. Mirbagheri, S.A., Tanji, K.K., Rajaee, T., 2008. Selenium transport and transformation modelling in soil columns and ground water contamination prediction. Hydrol. Process. 22 (14), 2475–2483. Mizutani, T., Kanaya, K., Osaka, T., 2001. Map of selenium content in soil in Japan. J. Health Sci. 47, 407–413. Morway, E.D., Gates, T.K., 2012. Regional assessment of soil water salinity across an intensively irrigated river valley. J. Irrig. Drain. Eng. 138 (5), 393–405. Morway, E.D., Gates, T.K., Niswonger, R.G., 2013. Appraising options to enhance shallow groundwater table and flow conditions over regional scales in an irrigated alluvial aquifer system. J. Hydrol. 495, 216–237. Moxon, A.L., Olson, O.E., Searight, W.V., Sandals, K.M., 1938. The stratigraphic distribution of selenium in the cretaceous formations of South Dakota and the selenium content of some associated vegetation. Am. J. Bot. 25 (10), 794–809. Myers, T., 2013. Remediation scenarios for selenium contamination, Blackfoot watershed, southeast Idaho, USA. Hydrogeol. J. 21, 655–671. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., 2011. Soil and Water Assessment Tool. Theoretical Documentation. Grassland, Soil and Water Research Laboratory, Agricultural Research Service, Temple, Texas. Niswonger, R.G., Prudic, D.E., Regan, R.S., 2006. Documentation of the Unsaturatedzone flow (UZF1) Package for Modeling Unsaturated flow Between the Land Surface and the Water Table with MODFLOW-2005, U.S. Geological Survey Techniques and Methods 6-A19. Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 2011. MODFLOW-NWT, A Newton Formulation for MODFLOW-2005: U.S. Geological Survey Techniques and Methods 6-A37, 44p. Oremland, R.S., Hollibaugh, J.T., Maest, A.S., Presser, T.S., Miller, L.G., Culbertson, C.W., 1989. Selenate reduction to elemental selenium by anaerobic bacteria in sediments and culture: biogeochemical significance of a novel, sulfateindependent respiration. Appl. Environ. Microbiol. 55, 2333–2343. Oremland, R.S., Steinberg, N.A., Maest, A.S., Miller, L.G., Hollibaugh, J.T., 1990. Measurement of in situ rates of selenate removal by dissimlatory bacterial reduction in sediments. Environ. Sci. Technol. 24, 1157–1164. Peng, A., Wang, Z., Whanger, P.D., Combs, G.F., Yeh, J.Y., 1995. In: Environmental Bioinorganic Chemistry of Selenium. C.E.S. Technology, Beijing, China. Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008. Global Sensitivity Analysis. The Primer. John Wiley and Sons. Schultz, L.G., Tourtelot, H.A., Gill, J.R., Boerngen, J.G., 1980. Geochemistry of the Pierre Shale and Equivalent Rocks of the late Cretaceous Age. USGS, Geological Survey Professional Paper 1064-A. Scott, G.R., 1968. Geologic and Structure Contour Map of the La Junta Quadrangle, Colorado and Kansas. U.S.G. Survey, Reston, Virginia. Seiler, R.L., 1995. Prediction of areas where irrigation drainage may induce selenium contamination of water. J. Environ. Qual. 24, 973–979. Seiler, R.L., 1997. Methods to Identify Areas Susceptible to Irrigation-Induced Selenium Contamination in the Western United States. US Dept. of the Interior, USGS, Washington, DC, p. 4, Fact Sheet FS-038-97. Sharps, J.A., 1976. Geologic Map of the Lamar Quadrangle, Colorado and Kansas. U.S.G. Survey, Reston, Virginia. Shrift, A., 1954. Sulfur-Selenium Antagonium. I. Antimetabolite Action of Selenate on the Growth of Chlorella vulgaris. Am. J. Bot. 41 (3), 8. Skorupa, J.P., 1998. Selenium poisoning of fish and wildlife in nature: Lessons from twelve real-world experiences. In: Frankenberger, W.T., Engberg, R.A. (Eds.), Enviromental Chemistry of Selenium. Marcel Dekker, New York, New York, USA, pp. 315–396.

46

R.T. Bailey et al. / Journal of Hydrology 515 (2014) 29–46

Stillings, L.L., Amacher, M.C., 2010. Kinetics of selenium release in mine waste from the Meade Peak phosphatic shale, Phosphoria formation, Wooley Valley, Idaho, USA. Chem. Geol. 269, 113–123. Stolz, J.F., Basu, P., Oremland, R.S., 2002. Microbial transformation of elements: the case of arsenic and selenium. Int. Microbiol. 5, 201–207. Tayfur, G., Tanji, K.K., Baba, A., 2010. Two-dimensional finite elements model for selenium transport in saturated and unsaturated zones. Environ. Monit. Assess. 169 (1–4), 509–518. Tokunaga, T.K., Brown, J.G.E., Pickering, I.J., Sutton, S.R., Bajt, S., 1997. Selenium redox reactions and transport between ponded waters and sediments. Environ. Sci. Technol. 31, 7. van Veen, J.A., Paul, E.A., 1981. Organic carbon dynamics in grassland soils. (1) Background information and computer simulation. Can. J. Soil Sci. 61, 185–201. Wang, Z., Gao, Y., 2001. Biogeochemical cycling of selenium in Chinese envrionments. Appl. Geochem. 16, 7. Weres, O., Bowman, H.R., Goldstein, A., Smith, E.C., Tsao, L., 1990. The effect of nitrate and organic matter upon mobility of selenium in groundwater and in a water treatment process. Water Air Soil Pollut. 49, 251–272.

White, A.F., Benson, S.M., Yee, A.W., Wollenberg, J.H.A., Flexser, S., 1991. Groundwater contamination at the Kesterson Reservoir, California 2. Geochemical parameters influencing selenium mobility. Water Resour. Res. 27, 1085–1098. Williams, K.T., Byers, H.G., 1934. Occurrence of Selenium in Pyrites. Ind. Eng. Chem., 296–297. Wriedt, G., Rode, M., 2006. Modelling nitrate transport and turnover in a lowland catchment system. J. Hydrol. 328, 157–176. Wright, W.G., 1999. Oxidation and mobilization of Selenium by Nitrate in Irrigation Drainage. J. Environ. Qual. 28, 1182–1187. Yeh, G.T., Tripathi, V.S., 1989. A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components. Water Resour. Res. 25 (1), 93–108. Zhang, Y., Moore, J.N., 1997. Reduction potential of selenate in wetland sediment. J. Environ. Qual. 26, 910–916. Zhang, H.H., Wu, Z.F., Yang, C.L., Xia, B., Xu, D.R., Yuan, H.X., 2008. Spatial distributions and potential risk analysis of total soil selenium in Guangdong Province, China. J. Environ. Qual. 37, 780–787.