Spatially-distributed influence of agro-environmental factors 1 governing nitrate fate and transport in an irrigated stream-2 aquifer system

Elevated levels of nitrate (NO 3 ) in groundwater systems pose a serious risk to human populations and natural ecosystems. As part of an effort to remediate NO 3 contamination in irrigated stream–aquifer systems, this study elucidates agricultural and environmental parameters and processes that govern NO 3 fate and transport at the regional (500 km 2 ), local (50 km 2 ), and field scales ( 2 ). Specifically, the revised Morris sensitivity analysis method was applied to a finite-difference nitrogen cycling and reactive transport model of a regional-scale study site in the lower Arkansas River valley in southeastern Colorado. The method was used to rank the influence of anthropogenic activities and natural chemical processes on NO 3 groundwater concentration, NO 3 mass leaching, and NO 3 mass loading to the Arkansas River from the aquifer. Sensitivity indices were computed for the entire study area in aggregate as well as each canal command area, crop type, and individual grid cells. Results suggest that fertilizer loading, crop uptake, and heterotrophic denitrification govern NO 3 fate and transport for the majority of the study area, although their order of influence on NO 3 groundwater concentration and mass leaching varies according to crop type and command area. Canal NO 3 concentration and rates of autotrophic denitrification, nitrification, and humus decomposition also dominate or partially dominate in other locations. Each factor, with the exception of O 2 reduction rate, is the dominating influence on NO 3 groundwater concentration at one or more locations within the study area. Results can be used to determine critical processes and key management actions for future data collection and remediation strategies, with efforts able to be focused on localized areas.


Introduction
During recent decades, an elevated concentration of nitrate (NO 3 )C NO 3 (expressed as equivalent N concentration in this paper) in groundwater systems and at points of groundwater discharge to surface water bodies has become a serious environmental issue due to its adverse effects on human populations and natural ecosystems (Spalding and Exner, 1993).Specific problems associated with high C NO 3 include methemoglobinemia for infants (Fan and Steinberg, 1996) and eutrophication in aquatic systems, which induces depletion of dissolved oxygen (O 2 ) (hypoxia) due to increased biological activity.In addition, high C NO 3 can lead to elevated concentrations of sulfate and selenium (Se) via oxidation of pyrite (FeS 2 ) and seleno-pyrite (FeSe 2 ) from marine shale (Frind et al., 1990;Jørgensen et al., 2009;Bailey et al., 2012).NO 3 also has been shown to mobilize uranium via oxidation (Wu et al., 2010).Recent studies have revealed that certain rock formations can yield nitrogen (N) in response to a variety of biogeochemical processes (Holloway and Dahlgren, 2002;Montross et al., 2013).In most cases, however, elevated concentrations result from excessive loadings of organic or inorganic N fertilizer, inducing NO 3 leaching to the saturated zone of the aquifer (Korom, 1992;Spalding and Exner, 1993).
Published by Copernicus Publications on behalf of the European Geosciences Union.

R. T. Bailey et al.: Environmental factors governing nitrate transport
To combat NO 3 contamination, numerous field and modeling studies have been performed to quantify NO 3 fate and transport processes in soil-groundwater systems, identify baseline conditions of N sources and transport patterns, and investigate potential remediation strategies.For the latter, simulation models typically are used to predict the effect of land use and best-management practices such as reduction in fertilizer loading (Chaplot et al., 2004;Almasri and Kaluarachchi, 2007;Lee et al., 2010), reduction in applied irrigation water (Ma et al., 1998;Rong and Xuefeng, 2011), and implementing or enhancing riparian buffer zones (Hefting and de Klein, 1998;Spruill, 2000;Vaché et al., 2002;Sahu and Gu, 2009) on overall C NO 3 and on NO 3 mass loading to and within streams.These studies have been conducted at various scales (Ocampo et al., 2006), ranging from the soil profile and field scale (Johnsson et al., 1987;Ma et al., 1998;Rong and Xuefeng, 2011), to the catchment scale (Birkinshaw and Ewen, 2000;Conan et al., 2003;Wriedt and Rode, 2006;Lee et al., 2010), to the regional-scale watershed or river basin scale (Chaplot et al., 2004;Almasri and Kaluarachchi, 2007;Bailey et al., 2015), and include a variety of fate and transport processes such as soil N cycling, leaching, groundwater transport, and overland transport.
Besides assessing baseline conditions and predicting domain-scale effects on spatial concentrations and loadings, numerical models also can be used in NO 3 remediation to determine the system inputs, parameters, and processes (i.e., model factors) that govern these concentrations and loadings.In general, identifying the most influential processes on resulting C NO 3 and mass loading can assist in establishing optimal remediation strategies.Additional benefits of the analysis include guiding effective field sampling strategies by focusing on influential system variables or inputs; facilitating model calibration and testing by focusing on the identified key factors (Sincock et al., 2003;Almasri and Kaluarachchi, 2007); identifying factors that require additional research to improve model performance (Hall et al., 2009); and detecting non-influential parameters or processes that possibly could be eliminated to simplify the model (Saltelli et al., 2008).
An appealing approach to determine the influence of model factors is sensitivity analysis (SA), which relates changes in model output variables (e.g., concentration, mass loading) to prescribed changes in model factor input values (e.g., initial conditions, system stresses, system parameters).For studies assessing NO 3 fate and transport in groundwater systems using physically based spatially distributed groundwater models, sensitivity analysis typically is performed in a simple fashion due to model complexity and computational cost.For example, Almasri and Kaluarachchi (2007) increased values of selected parameters (e.g., denitrification rate, longitudinal dispersivity, initial concentration, soil mineralization rate, soil nitrification rate, fertilizer loading) by 50 % to determine their influence on simulated C NO 3 in a watershed in Washington, USA, and Ehteshami et al. (2013), using the LEACHN (N module of the Leaching Estimation and Chemistry Model) model, investigated the influence of low and high values of rainfall and initial C NO 3 for two soil types on soil C NO 3 .In a field study using the nitrate transport model RISK-N (nitrate transport) model, Oyarzun et al. (2007) modified values of soil initial N, C NO 3 in irrigation water, fertilizer, N crop uptake, crop evapotranspiration (ET), and soil properties by 50, 70, 100, 125, and 150 % to investigate their influence on NO 3 vadose zone mass flux and C NO 3 in the groundwater.Also, Hartmann et al. (2013) used SA to estimate the influence of model parameters on the time lag between spring discharge and NO 3 at several karst aquifer sites across Europe.Whereas global effects of the model factor on system-response variables can be assessed, local and interaction effects cannot be quantified using this approach.
A more rigorous SA method is global sensitivity analysis (GSA), which searches the entire parameter space to identify the importance of model parameters and interactions thereof.Such methods include the elementary effects (EE) method (Morris, 1991;Cacuci, 2003), a screening method that identifies the most important model factors and is wellsuited for large models (Campolongo and Braddock, 1999), and variance-based methods that quantitatively decompose the variance of model output into fractions that are attributed to model factors (Saltelli et al., 2008).A number of hydrologic modeling studies have used GSA methods for assessing model factor influence on overall watershed nutrient and sediment processes (White and Chaubey, 2005;Arabi et al., 2007;Sun et al., 2012;Ahmadi et al., 2014), flooding and hydraulic characteristics (Hall et al., 2005(Hall et al., , 2009)), in-stream water quality (Cox and Whitehead, 2005;Deflandre et al., 2006;Liu and Zou, 2012;Bailey and Ahmadi, 2014), and in-stream solute transport (Kelleher et al., 2013).
Sensitivity analysis is commonly used in hydrologic and water quality modeling to identify the influence of model parameters on an aggregated measure of model responses such as average annual stream discharge or contaminant loads.A few studies have assessed how the results of SA vary in time.For example Reusser et al. (2011) used hydrologic catchment models to investigate the temporal-varying influence of model factors on a variety of watershed response variables for catchments in Ecuador and Germany.However, the spatial variability of sensitivity indices has been largely neglected.Specifically regarding this study, no studies have quantified the spatial-varying influence of factors on solute concentrations in large-scale groundwater systems.Such information could be valuable in terms of implementing sitespecific remediation strategies, facilitating model calibration for specific model domain regions, and identifying system variables that require additional field data collection, particularly for NO 3 due to its ubiquitous presence in groundwater systems worldwide.
This study aims to identify the spatially varying influence of system factors on NO 3 fate and transport in a regionalscale (506 km 2 ) irrigated hydro-agricultural system.Specifically, the factors' influence on NO 3 groundwater concentra- tions, NO 3 leaching below root zone, and NO 3 groundwater mass loading to the stream network will be quantified for a variety of scales (cultivated field, canal command area, region).A calibrated and tested N fate and transport groundwater model is used for the assessment, with the modified Morris method used for the sensitivity analysis.

Methods
A comprehensive SA method was applied to a regionalscale, intensively irrigated 506 km 2 groundwater system in the lower Arkansas River valley (LARV) in southeastern Colorado to identify the spatially varying influence of system factors on NO 3 concentrations in groundwater, NO 3 mass leaching in the shallow soil zone, and NO 3 mass loading to the Arkansas River.The model used is UZF-RT3D (Reactive Transport in 3 Dimensions in saturated and Unsaturated Zone Flow systems) (Bailey et al., 2013a, b), a MOD-FLOW (Harbaugh, 2005)-based (Niswonger et al., 2011), finite-difference model designed for N fate and transport at the regional scale and recently applied to the study area (Bailey et al., 2014).The model accounts for major agricultural inputs (fertilizer, canal seepage, and irrigation water), processes (N cycling in the root and soil zone, leaching, threedimensional transport, heterotrophic and autotrophic denitrification), and outputs (mass loading to the stream network).
As identifying the relative importance of parameters and processes in space is the objective of this study, and since computational costs of UZF-RT3D are extremely high (runtime of approximately 3.5 h for a single simulation using an Intel ® Core ™ i7-3770 CPU @ 3.40 GHz desktop computer), the SA method used is an improved variant (Cam-polongo et al., 2007) of the Morris method (Morris, 1991) rather than variance-based SA methods such as Sobol (Sobol, 1993) or FAST (Fourier amplitude sensitivity test) (Cukier et al., 1973).Nine model factors are included in the assessment, with their overall influence on NO 3 fate and transport evidenced in a previous study in the region (Bailey et al., 2014).In conjunction with the SA methodology, model results are processed to determine the dominant model factors globally (i.e., averaged for the entire model domain), for each irrigation canal command area, for each crop type (i.e., the set of model grid cells associated with each crop type), and for each grid cell, thereby elucidating parameter influence at varying spatial scales.For the latter, spatial contour maps depicting model sensitivity to individual model factors are shown.Due to the dependence of N fate and transport on the presence of O 2 , the influence of the nine model input factors on C O 2 also is calculated and presented.

Study area
The semi-arid LARV in Colorado extends from the outlet of the Arkansas River from Pueblo Reservoir eastward across southeastern Colorado to the border with Kansas (Fig. 1), with the Arkansas River fed primarily by snowmelt from the mountainous regions of the upper Arkansas basin.In total, the valley supports approximately 109,000 irrigated ha (270 000 ac), and is one of Colorado's most productive agricultural areas.Approximately 14 000 fields are cultivated, with the majority using flood irrigation methods and a small minority using sprinklers or drip irrigation methods.Major crops include alfalfa, corn, grass hay, wheat, sorghum, dry beans, cantaloupe, watermelon, melons, and onions.The region of the LARV focused on in this study is shown in Fig. 1.The boundary of the study area is shown with a black line, and encompasses an area of 50 600 ha (125 000 ac), of which 26 400 ha (65 300 ac) are irrigated.The fields receiving water from each of six main irrigation canals (i.e., canal command areas) are shown in Fig. 2a, with crop type cultivated in 2006 for each field shown in Fig. 2b.Due to over-irrigation and poor subsurface drainage, high water table elevations have been established in recent decades, with water table depth below ground surface often between 1 and 3 m (Morway and Gates, 2012).These high water tables have resulted in salinization and waterlogging, in addition to substantial rates of groundwater return flows (i.e., discharge) to the Arkansas River and its tributaries (Morway et al., 2013).The thickness of the alluvial aquifer ranges from 4 to 34 m (Fig. 4a), and is underlain by Cretaceous shale (Scott, 1968;Sharps, 1976) in both solid and weathered form.
In addition to salinization and associated decrease in crop productivity (Morway and Gates, 2012), elevated groundwater C No 3 has been observed, presumably due to overfertilization on cultivated fields.In a similar irrigated region of the LARV, located about 67 km upstream, Zielinski et al. (1997) examined δ 15 N isotopic signatures to conclude that NO 3 was derived primarily from fertilizer and crop waste, not from proximate geologic sources.To assess the C NO 3 in the study region, groundwater and surface water samples were collected (see locations in Fig. 2a) during 10 sampling events over the period 2006-2009(Gates et al., 2009)).For groundwater, samples were taken routinely from 52 observation wells, with groundwater from 37 additional observation wells sampled non-routinely (aperiodic).Surface water samples were taken from 10 locations along the Arkansas River and 5 locations in tributaries.Detailed results of the monitoring scheme are shown in the Supplement.In summary, for groundwater the 85th percentile values of C NO 3 −N were at or in excess of the 10 mg L −1 (85th percentile) EPA drinking water standard for the first three sample trips.The maximum measured value was 66 mg L −1 .The means for the samples gathered from the Arkansas River and its tributaries were 1.53 and 1.95 mg L −1 , respectively.The annual median values of the Arkansas River samples were 0.95, 1.20, 1.10, and 2.20 mg L −1 for each of the successive years within the period 2006-2009, compared to the Colorado interim standard of 2 mg L −1 (CDPHE, 2012) for total inorganic N concentration (C NO 3 -N + C NO 2 -N + C NH 4 -N ) C NO 3 -N .The concentration of C NO 3 -N C NO 3 -N exceeded 2 mg L −1 in about 25 % of the samples gathered in the river over this period and exceeded 2.5 mg L −1 in about 12 % of the samples, signifying the growing concern about N pollution in the river.Analysis of 22 river samples and 15 tributary samples in 2013 revealed that C NO 3 -N C NO 3 -N made up greater than 80 % of total dissolved N in the river and about 76 % of total dissolved N in the tributaries.

UZF-RT3D N reaction module and baseline application
UZF-RT3D simulates the reactive transport of multiple interacting chemical species in variably saturated porous media using groundwater flow rates, water content, and a variety of groundwater sources and sinks (e.g., applied irrigation water, pumping, canal seepage, or groundwater-surface water interactions) simulated by a MODFLOW-NWT (Newton) model using the UZF1 package.The N cycling and reaction module add-on package (Bailey et al., 2013b) was designed for model application in an irrigated agricultural groundwater system, and accounts for the major hydrologic, chemical, and land management processes that govern N fate and transport in an irrigated stream-aquifer system.Also, due to the A schematic of the fate and transport of N species and O 2 as simulated by the N reaction module of UZF-RT3D is depicted in Fig. 3a.N mass (NO 3 or NH 4 ) enters the subsurface via fertilizer loading (single application or split application), canal seepage, infiltrating irrigation water (either from canal water or pumped groundwater), or seepage from the stream network (Arkansas River and its tributaries).N mass exits the subsurface via groundwater discharge to the stream network.N cycling occurs in the root and soil zone, with organic N and carbon (C) added to soil organic matter (manure M N , fast-decomposing litter L N , flow-decomposing humus H N ) via after-harvest plowing or decaying root mass and subsequently mineralized to NH 4 , which can be volatilized, nitrified to NO 3 , or taken up with NO 3 into crop roots during the growing season.The timing of land management actions, e.g., fertilizer loading (40 %, 60 % split application), irrigation events, harvesting, and plowing, adopted in the module is shown in Fig. 3b.NH 4 is sorbed readily to soil surface sites, whereas NO 3 is transported by one-dimensional transport in the unsaturated zone and three-dimensional transport in the saturated zone, subject to heterotrophic denitrification in near-surface areas and autotrophic denitrification in the presence of FeS 2 -bearing marine shale (see Fig. 1).O 2 is also subject to heterotrophic and autotrophic chemical reduction.
UZF-RT3D solves a system of advection-dispersionreaction (ADR) equations for interacting dissolved-phase and solid-phase species using the finite-difference approach.Including ADR processes and source/sink terms as depicted, the following mass conservation equations are written for the dissolved-phase species (NO 3 , NH 4 , O 2 ) in the N reaction module: where with "f" denoting fluid phase; v is the pore velocity [L b T −1 ], provided by MODFLOW-UZF1; θ is the volumetric water con-

R. T. Bailey et al.: Environmental factors governing nitrate transport
, also provided by MODFLOW-UZF1; D ij is the hydrodynamic dispersion coefficient [L 2 T −1 ]; q f is the volumetric flux of water representing sources and sinks [L 3 f T −1 L −3 b ] such as irrigation water, canal and river seepage, groundwater discharge to the river, or pumped groundwater, with b denoting the bulk phase; C f is the concentration of the source or sink with s denoting the solid phase, and is equal to 1 − φ, where φ is porosity [L 3 f L −3 b ]; r f represents the rate of all reactions that occur in the dissolved-phase [M f L −3 f T −1 ]; min, imm, nit, and vol signify mineralization, immobilization, nitrification, and volatilization, respectively; and auto and het represent autotrophic and heterotrophic chemical reduction, respectively.ε is included for the min and imm reactions to denote a mass transfer between the solid and dissolved phases.For NH 4 , which is subject to sorption, R is the retardation factor and is equal to 1 . The daily mass of potential N crop uptake during the growing season is determined using a logistic equation (Johnsson et al., 1987) and is distributed across the vertical column of grid cells encompassing the crop rooting depth according to the mass density of the root system.Mass conservation equations (not shown) for solidphase organic N (and C) species L N , H N , and M N also are implemented.
The rate of chemical reactions r f included in Eqs.
(1)-( 3) is governed by the dependence of the chemical reaction on soil temperature T , θ , and the presence of O 2 and C.These rates are simulated using first-order Monod kinetics.For example, the following rate law expression represents the process of heterotrophic denitrification, with others contained in Bailey et al. (2015): where λ is the base rate constant for the reaction [T f ] signifying the species concentration at which lower-redox species can undergo appreciable rates of reduction; CO 2,prod is the total mass of CO 2 produced during organic matter decomposition and is used as an indicator of available organic carbon (OC) for microbial consumption (Birkinshaw and Ewen, 2000); and E [-] is an environmental reduction factor that accounts for θ and T and acts to temper microbial activity rates (Birkinshaw and Ewen, 2000;Bailey et al., 2013b).Nitrification, mineralization, and denitrification each have uniquely specified relationships between θ and microbial activity.
The UZF-RT3D model used in this study is the same as that described in Bailey et al. (2014).The model uses output from a calibrated and tested MODFLOW-NWT (Niswonger et al., 2011) model of the study region (Morway et al., 2013), which uses the UZF1 unsaturated-zone flow package (Niswonger et al., 2006).The flow model uses weekly estimates of irrigation water, precipitation, canal seepage, and crop ET to estimate groundwater level and groundwater-surface water interactions for the 1999-2009 time period.Figure 4b-d show the finite-difference grid, the simulated water content of the soil in June 2006, and the average simulated water table elevation (m) during the 1999-2009 time period, respectively.
The UZF-RT3D model uses the same model domain and finite-difference grid as the flow model (see Fig. 4b).The model has seven vertical layers, with layers 1-2 (0.5 m each) corresponding to the root zone, layer 3 (1.0 m) corresponding to the leaching zone, layers 4-6 to the saturated zone, and layer 7 to the shale bedrock formation.Thickness of layers 4, 5, and 6 varies according to saturated thickness, with layer thickness ranging from 2.8 to 12.6 m.Each vertical column of cells in the three-dimensional grid is assigned a set of crop parameter values according to the portions of fields within the grid cell area.Crop parameters, with values shown in Table 1 for each crop type in the study area, include planting day; harvest day; plowing day; mass of stover plowed into the soil P St (kg ha −1 ) after harvest; maximum rooting depth d rt,max (m), which controls N uptake; C-N ratio of root mass CN RT ; fertilizer loading F NH 4 (kg ha −1 ); maximum seasonal uptake values of N N up (kg ha −1 ); depth of plowing d pw (m); mass of decaying roots P Rt (kg ha −1 ); C-N ratio of stover mass CN ST ; and constants defining root growth and daily uptake rate U .Chemical reaction parameter values are shown in Table 2, with an asterisk ( * ) indicating the mean value of all the grid cells.C NO 3 and C O 2 of canal water and irrigation water were based on observed data.The model was run for the 2006-2009 and tested against spatiotemporal averages of groundwater C NO 3 and NO 3 mass loadings from the aquifer to the Arkansas River.

Morris SA methodology
The Morris screening method for global SA is based on an individually randomized one-at-a-time (OAT) design that provides information regarding (i) the main effect of each input parameter on model output responses and (ii) the overall effects including interactions between parameters.For example, consider a model M with a vector of k parameters (ω i , i = 1, . . ., k) within the feasible parameter space, , that simulates m response vectors of the system (S j,j = 1, . . ., m): [S 1 , . .., S m ] = M (ω 1 , . .., ω k ) .
(  Similar to any standard SA practice, parameters are drawn from their predefined probability distributions, with each model input parameter ω i varied across p discrete values (Saltelli et al., 2008).Generally, results of SA are not sensitive to the choice of distribution from which values are sampled.After running model M for the given parameter sets, the local sensitivity measure (also referred to as the elementary effect, EE) is then computed for each parameter i for model response j as follows: where is a value in the predefined increments (i.e., [1/(p − 1), . . ., 1 − 1/(p − 1)][1/(p − 1), . . ., 1 − 1/(p − 1)]) and ω = ω 1 , . . ., ω k is a random sample in the parameter space so that the transformed point θ 1 , . . ., θ i−1 θ i + , . . .θ k (ω 1 , . . ., ω i−1 , ω i + , . . .ω k ) is still within the parameter space (Saltelli et al., 2008).The - * Indicates mean value, with specific values assigned to each command area according to the values reported in Bailey et al. (2014).
resulting distribution EE i associated with each parameter ω i is then analyzed to determine µ, the mean of the distribution that assesses the overall importance of the parameter on the model output, and σ , the standard deviation of the distribution, which indicates nonlinear effects and/or interactions (Campolongo et al., 2007).
To determine sensitive and insensitive values, it is recommended to evaluate a graphical representation of σ vs. µ.However, for non-monotonic models, some EE values with opposite signs may cancel out when µ is calculated, and hence Campolongo and Saltelli (1997) proposed the use of µ * , the sample mean of the distribution of absolute values of EE. µ * includes all types of effects that parameters can have on output responses and, therefore, is a global measure of output sensitivity to the parameters (Campolongo et al., 2007).µ * i,j is defined as the mean of absolute values of the computed elementary effects EE i,j .The total computational cost of the Morris experiment is n = r(k + 1) runs, where r is the selected size of each sample.
As noted above, an important objective of SA is to determine the most influential model input parameters.Hence, it is important to measure the level of agreement between results of SA experiments with an emphasis on the high-ranked parameters.Campolongo and Saltelli (1997) suggested the use of the Savage score to facilitate comparison of results from different SA experiments (see next section).The Savage score is defined as follows (Iman and Conover, 1987): where i is the rank assigned to the ith model parameter based on the Morris µ * ; for example, the highest ranked variable would have a score of 1/1 + 1/2 + 1/3 + . . .+ 1/k, the second ranked variable would have a score of 1/2 + 1/3 + . . .+ 1/k, etc. Savage scores typically are pre-ferred because they place higher emphasis on the agreement of the key drivers (i.e., higher ranked parameters), rather than the overall agreement.The Savage score can be used in aggregating the results from different SA methods.

Model input factors analyzed
In applying the SA method to the UZF-RT3D model of the study area, nine model input factors were analyzed for impact on model results: F NH 4 , N up , C NO 3 in canal water Canal NO 3 , rate of litter pool decomposition λ L , rate of humus pool decomposition λ H , rate of autotrophic reduction of O 2 in the presence of shale λ auto O 2 , rate of autotrophic reduction of NO 3 in the presence of shale λ auto NO 3 , rate of nitrification λ nit , and rate of heterotrophic denitrification λ het NO 3 .Canal NO 3 conveys NO 3 mass into the subsurface system via applied irrigation water as well as seeped canal water.For each simulation, separate values of F NH 4 and N up were generated for each crop type, separate values of Canal NO 3 were generated for each of the six canal command areas, and separate values of λ auto O 2 , λ auto NO 3 , and λ nit were generated for each command area.The mean of each parameter value is derived from the baseline simulation (see Tables 1 and 2), with the mean values of λ auto O 2 , λ auto NO 3 , and λ nit for each command area estimated during the calibration phase (Bailey et al., 2014).
Setting the number of replications r and levels p of the Morris scheme to 20 and 10, respectively, a total of 280 simulations were run.Parameter values were perturbed using a coefficient of variation (CV) of 0.2 for all parameters except for Canal NO 3 , which was perturbed with a CV of 0.1 based on variance in observed canal water concentrations.Perturbation for the reaction rates (λ L , λ H , λ auto O 2 , λ auto NO 3 , λ het NO 3 , λ nit ) was performed using log values since statistically these rates typically conform to a lognormal distribution (Parkin and Robinson, 1989;McNab Jr. and Dooher, 1998) 5, with averages of 250 kg ha −1 , 1.055 × 10 −4 day −1 , and 2.6 g m −3 , respectively.The values shown in Fig. 5a are for grid cells that contain corn, and the values shown in Fig. 5b and c are for the grid cells within the Rocky Ford Highline canal command area (canal feeding the gray-shaded fields in Fig. 2a).
For each of the 280 simulations, the model was run for a 2-year spin-up period, followed by the 2006-2009 period.Model results were processed to determine the influence of the nine targeted model input factors on groundwater C NO 3 , NO 3 mass leached from the root zone, and total NO 3 mass loading to the Arkansas River from the aquifer.Post-processing was implemented to determine this influence (i) globally for the entire study area, i.e., averaging values from all grid cells; (ii) for individual crop types, i.e., averaging values from all grid cells corresponding to a given crop type; (iii) for individual canal command areas, i.e., averaging values from all grid cells within a given command areas; and (iv) for individual grid cells.As total NO 3 mass loading to the Arkansas River occurs along the entire reach of the river within the study area, parameter influence is assessed only for (i).Values of average concentration, average leaching, and total mass loading were processed from the final year of the model simulation (i.e., 2009).For groundwater C NO 3 , concentration values were taken from layer 4 of the model, which corresponds to the depth of observation well screens in the study area.For NO 3 leaching, values were taken from layer 3 (i.e., the mass leached from layer 3 to layer 4).For parameter influence on C NO 3 for individual grid cells (item iv), the Savage score as calculated by Eq. ( 7) was used for presentation of results.Also for (iv), the parameter influence on C O 2 was presented.

General model results
Model results from 1 of the 280 simulations is shown in Fig. 6, with spatial distribution of C O 2 and C NO 3 shown in Fig. 6a and b, respectively for 22 July 2009, and the spatial distribution of NO 3 mass loading shown for 1 week during the winter (2 December 2006, Fig. 6c) and 1 week during the summer (10 August 2008, Fig. 6d).Mass loadings from the aquifer to the stream network (discharge) are displayed in red, whereas loadings from the stream network to the aquifer (seepage) are displayed in green.For concentrations in groundwater, values of C O 2 range from 0.0 to 10.3 mg L −1 , with an average value of 2.7 g m −3 for the 7776 active grid cells.Values of C NO 3 range from 0.0 to 78.3 mg L −1 , with an average value of 1.84 mg L −1 .
Hotspots occur for both C O 2 and C NO 3 , with those of C NO 3 typically occurring in locations of corn cultivation due to (1/day) and (c) nitrate concentration of canal water Canal NO 3 (mg L −1 ) for the Rocky Ford Highline canal command area, for each of the 280 simulations in the revised Morris SA scheme.the higher loading of F NH 4 as compared to other crop types.NO 3 mass loadings occur along the Arkansas River and the tributaries, with discharge and seepage both occurring along the length of the canals during the summer (Fig. 6d).The spatiotemporal average value of C NO 3 in groundwater for each command area during the entire 2006-2009 time period is shown in Fig. 7 for each of the 280 simulations.The average value for all grid cells in non-cultivated areas is also shown.Average C NO 3 across all simulations for each command area are (average of observed field values are in parentheses) Rocky Ford Highline: 2.0 mg L −1 (3.1 mg L −1 ); Catlin: 1.4 mg L −1 (6.1 mg L −1 ); Rocky Ford: 1.5 mg L −1 (3.8 mg L −1 ); Fort Lyon: 3.7 mg L −1 (1.6 mg L −1 ); Holbrook: 1.9 mg L −1 (3.5 mg L −1 ); and non-cultivated areas: 3.5 mg L −1 (4.2 mg L −1 ).Average values correspond closely to results from the tested baseline model (Bailey et al., 2014).

Parameter influence on global concentration, leaching and loading of NO 3
The global influence of the nine model input factors on NO 3 fate and transport in the study area is shown in Fig. 8. Global sensitivity plots are used, with nonlinear effects and/or interactions σ plotted against mean µ * .The influence of the  factors on C NO 3 in layer 1 (top 0.5 m of the root zone), C NO 3 in layer 4 (shallow saturated zone), NO 3 leaching from layers 3 to 4 L NO 3 Lay3→4 (generally from the unsaturated zone to the saturated zone), and total NO 3 mass loading to the Arkansas River Load NO 3 are shown in Fig. 8a-d, respectively.As seen in Fig. 8a, C NO 3 in the root zone is governed principally by fertilizer loading (F NH 4 ) and seasonal N up and to a smaller degree by λ het NO 3 and λ nit .In the shallow saturated zone (Fig. 8b), where NO 3 mass is received from the upper soil zone via leaching, F NH 4 and N up still are dominant, but the concentration of NO 3 in the canals (Canal NO 3 ) has a stronger direct impact than λ het NO 3 .The rate of humus decomposition (λ H ) and autotrophic denitrification (λ auto NO 3 ) also have a slight impact.NO 3 leaching is also governed by F NH 4 , N up , λ het NO 3 , Canal NO 3 , and λ H (Fig. 8c), as higher F NH 4 , lower N up , lower λ het NO 3 , and higher Canal NO 3 increase the mass of NO 3 leached, and vice versa.Load NO 3 is governed by F NH 4 , N up , and λ het NO 3 (Fig. 8d), with λ het NO 3 influencing not only how much NO 3 is leached to the water table and carried to the stream network via groundwater flow, but also how much NO 3 undergoes denitrification in the riparian areas of the stream network.
The high σ values for N up , F NH 4 , λ het NO 3 and Canal NO 3 shown in Fig. 8 signify the large spread in EE values for these parameters, indicating that their influence on C NO 3 , NO 3 leaching, and NO 3 mass loading is strongly dependent on the values of other parameters.For example, in reference to C NO 3 in the shallow saturated zone (Fig. 8b), the value of µ * for N up signifies the average effect of N up on C NO 3 , but some values of EE for N up are much smaller and larger than µ * .Smaller values of EE indicate that the combined influence of other parameter values produced a small effect of crop uptake on C NO 3 , such as a lower N fertilizer loading and higher rates of denitrification, whereas larger values indicate that other parameters produced a larger effect of crop uptake on C NO 3 , such as a higher N fertilizer loading and lower rates of denitrification.Also, higher values of Canal NO 3 increase the influence of crop uptake on C NO 3 , as more NO 3 mass is brought into the soil zone via canal seepage and infiltrating irrigation water.

Parameter influence on C NO 3 and leaching for each crop type
The influence of each of the nine parameters on C NO 3 in the shallow groundwater zone and on NO 3 leaching for each crop type in the study area is summarized in Tables 3 and 4, respectively, using values of µ * .The µ * values of the three most influential parameters for each crop type are in bold.
For the majority of crop types, C NO 3 the shallow groundwater zone is governed by N fertilizer loading (F NH 4 ), seasonal crop N uptake (N up ), and heterotrophic denitrification λ het NO 3 (Table 3), similar to the global analysis of C NO 3 in the shallow soil layers as presented in Section 3.2.For example, µ * for F NH 4 , N up , and λ het NO 3 is 0.94, 0.72, and 0.30, respectively, for corn-cultivated areas, and 0.84, 0.81, and 0.28 for sorghum-cultivated areas.The exception is areas that cultivate onion, in which Canal NO 3 (µ * = 0.45) ranks in the top three behind F NH 4 (1.21) and N up (0.99).This is due to the fact that onions receive irrigation water (from the canals) more frequently than do the other crop types.Hence, the C NO 3 of the canal water has a stronger influence on groundwater C NO 3 underlying onion-cultivated fields than for the other crop types.For many of the crops, λ H and λ nit have a small to moderate influence, whereas litter pool decomposition rate (λ L ), autotrophic reduction of O 2 (λ auto O 2 ), and autotrophic denitrification (λ auto NO 3 ) have a negligible to small influence on C NO 3 .The influence of the nine parameters on NO 3 mass leaching to the shallow saturated zone (Table 4) follows the same pattern as for their influence on C NO 3 , with fertilizer N loading, uptake, and denitrification dictating the amount of NO 3 leached to the water table (values in boxes) and canal concentration, nitrification, and humus and litter pool decomposition having small to moderate values of µ * .For corncultivated areas, the average effect µ * of F NH 4 , N up , and

Parameter influence on C NO 3 and leaching in individual canal command areas
Summaries of the influence of each of the nine parameters on C NO 3 in the shallow groundwater zone and on NO 3 leaching for each canal command area also are provided in Tables 3 and 4. The results show important differences between the command areas, with a mixture of F NH 4 , N up , λ nit , λ het NO 3 , λ auto NO 3 , and Canal NO 3 providing noteworthy impacts on C NO 3 and NO 3 mass leaching.For influence on C NO 3 (Table 3), the top three influential parameters within the Catlin command area are N up (µ * = 0.26), λ nit (0.16), and F NH 4 (0.12), whereas the top three for the Rocky Ford command area are Canal NO 3 (0.51), λ auto NO 3 (0.20), and N up (0.15), with the strong influence of λ auto NO 3 due to the presence of outcropped shale in the command area and hence locations of autotrophic denitrification.λ auto NO 3 also has a strong influence in the Holbrook command area, with the third highest value of µ * (0.11).Canal NO 3 is ranked third or higher in terms of µ * in three of the six command areas (Rocky Ford, Otero, Rocky Ford Highline).F NH 4 , N up , and λ het NO 3 govern NO 3 mass leaching for each of the command areas (Table 4) except for the Catlin command area, in which λ nit is ranked second (µ * = 38.0)and the Rocky Ford, in which Canal NO 3 is ranked first ( µ * = 30.3).

Spatial distribution of parameter influence on C NO 3 and C O 2
Cell-by-cell plots of Savage scores for the parameters according to their ranking in influencing C NO 3 in shallow groundwater are shown in Fig. 9. Plots are presented for each of the targeted nine parameters except for λ auto O 2 due to the negligible influence of O 2 autotrophic reduction on C NO 3 .The value for each cell represents the ranking (1-9) and associated Savage score for the given parameter.High ranking is displayed in maroon-red coloring, whereas low ranking is displayed in blue.As seen in the plots, the ranking of each parameter in its influence on groundwater C NO 3 is highly spatially variable.For example, the locations where canal NO 3 concentration (Canal NO 3 ) has the strongest influence (maroon coloring) (Fig. 9b) are scattered throughout the region, with entire local areas (encompassed by circles in Fig. 9b) governed by this parameter.For the cultivated areas, the dominant inputs/processes are fertilizer loading (Fig. 9a), crop N uptake (Fig. 9d), and heterotrophic denitrification (λ het NO 3 ) (Fig. 9e), with humus decomposition (Fig. 9g) having a moderate influence and litter decompo-sition (Fig. 9h) having a small influence.Whereas fertilizer loading and N uptake have the most influence on C NO 3 in most of the cultivated areas, some areas are governed principally by heterotrophic denitrification and humus decomposition (cells colored in maroon in Fig. 9e and g).Denitrification is particularly important in riparian areas along tributaries and the Arkansas River (Fig. 9e), where dense vegetation provides a natural filter of NO 3 before being loaded to surface water.Values of humus decomposition (λ H ) and litter decomposition (λ L ) control the rate of organic C and organic N decomposition and hence the availability of C for heterotrophic denitrification to proceed.
No area has λ L being the dominant influence on C NO 3 .Nitrification rate has a strong impact on C NO 3 in the Holbrook command area (red-pink cell coloring in Fig. 9c), with small impact elsewhere in the study area.Autotrophic denitrification is the dominant parameter in areas along the Arkansas River and several of the tributaries (Fig. 9f) that are adjacent to shale formations (see Fig. 1).However, it is interesting to note that there are many locations in the study area adjacent to outcropped shale in which λ auto NO 3 is not the dominant parameter.These locations are indicated by circles in Fig. 9f.In these areas, other system inputs and processes such as F NH 4 , N up , λ het NO 3 , and λ H are the governing influences on C NO 3 , demonstrating that knowledge of shale locations alone cannot be used to determine where C NO 3 will be affected the most by autotrophic denitrification.
Similar cell-by-cell plots of parameter Savage scores are shown in Fig. 10 for influence on C O 2 in shallow groundwater.λ H and λ L govern C O 2 in the cultivated areas (Fig. 10c  and d), with F NH 4 (Fig. 10b), N up (Fig. 10e) and Canal NO 3 (Fig. 10a) exhibiting small to moderate influence on C O 2 in the cultivated areas.The strong influence of λ H and λ L occurs due to their control of the rate of organic C decomposition, and hence the availability of C for heterotrophic reduction of O 2 .The rate of autotrophic reduction of O 2 (λ auto O 2 ) is dominant in localized areas where shale is present (see maroonshaded cells in Fig. 10f) with small influences in other areas of the study region, mainly in areas down gradient of the shale areas.

Discussion of results
Results provide information regarding the system inputs and processes that control NO 3 fate and transport generally (across the entire study region), by crop type, by canal command area, and by local regions.For the entire study region, detailed field sampling and observation of N fertilizer loading, N crop uptake, heterotrophic denitrification in the shallow soil layers, and concentration of NO 3 in canal water must be performed as often as possible to provide accurate model input data.NO 3 in canal water not only seeps through the perimeter of the earthen irrigation canals into the aquifer, but also is loaded to cultivated fields via applied irrigation wa- ter.In addition, results indicate these inputs and processes must be controlled via implemented management practices if NO 3 groundwater concentration, NO 3 leaching, and NO 3 mass loading to the river network are expected to decline in future decades, whereas other processes (rates of organic N decomposition, nitrification of NH 4 ) are not critical target factors.
These results agree with other previous studies from regions worldwide, which indicated that key controls on NO 3 fate and transport in groundwater and watershed systems, and hence targets for management action, include N fertilizer application (Chaplot et al., 2004;Botter et al., 2006;Al-masri and Kaluarachchi, 2007;Arabi et al., 2007;Bailey et al., 2015) and rate of denitrification (Wriedt and Rode, 2006;Almasri and Kaluarachchi, 2007;Schilling et al., 2007), with the order of their influence varied depending on the study region.However, these studies did not analyze the influence of NO 3 in canal irrigation water or the influence of crop N uptake.Molénat and Gascuel-Odoux (2002) demonstrated the strong influence of NO 3 leaching on in-stream NO 3 concentration, similar to our assessment of N uptake and denitrification (which influence NO 3 leaching) on NO 3 loading from the aquifer to the stream network.The same system parameters that govern NO 3 fate and transport at the regional scale also govern NO 3 for each individual crop type.N fertilizer loading (less), N crop uptake (more), and heterotrophic denitrification (more) typically must be controlled to decrease groundwater NO 3 concentration and NO 3 leaching, with NO 3 concentration in canal water controlled to lower these values for onion-cultivated areas.For canal command areas, N fertilizer loading and N uptake must be managed to decrease groundwater NO 3 concentration and NO 3 mass leaching in the majority of command areas.However, nitrification of NH 4 is an important control for the Catlin command area, NO 3 concentration in canal water is important for the Rocky Ford Highline, Otero, and Rocky Ford command areas, heterotrophic denitrification is important for each command area except Catlin and Rocky Ford Ditch, and autotrophic denitrification is important for only the Holbrook and Rocky Ford command areas.These reaction rate parameters must be focused on in field data monitoring schemes and in model parameter estimation.Results demonstrate that targeted inputs/outputs and processes vary depending on command area.
Similarly, different targets are required for controlling NO 3 fate and transport in localized areas throughout the study region.In reference to Fig. 9, each system parameter, with the exception of litter pool decomposition, is the most influential in controlling NO 3 fate and transport in at least several areas within the study region.N fertilizer loading is the dominant parameter in the majority of cultivated areas, although N uptake, heterotrophic denitrification, and NO 3 concentration in canal water also are the most influential in much of the study area.The rate of autotrophic denitrification (λ auto NO 3 ) is influential in many of the areas adjacent to outcropped marine shale.However, it is interesting to note that there are many locations in the study area adjacent to outcropped shale in which λ auto NO 3 is not the dominant parameter.These locations are indicated by circles in Fig. 9f.In these areas, other system inputs and processes are dominant, demonstrating that knowledge of shale locations alone cannot be used to determine where groundwater NO 3 concentration will be affected the most by autotrophic denitrification.
Whereas other studies (Chaplot et al., 2004;Botter et al., 2006;Wriedt and Rode, 2006;Almasri and Kaluarachchi, 2007;Arabi et al., 2007;Schilling et al., 2007;Bailey et al., 2015) have focused on the response of the entire groundwater and/or watershed system, the novelty of this study is the assessment of NO 3 transport control in localized areas within a region.Almasri and Kaluarachchi (2007) stated that the importance of denitrification in controlling NO 3 in groundwater may differ from location to location.In this study we quantify this difference spatially for denitrification and for each of the other eight targeted parameters (see Fig. 9).

Summary and concluding remarks
This study used a 506 km 2 regional-scale N fate and transport numerical model to examine the influence of forcing terms (fertilizer loading, crop N uptake, N concentration of applied irrigation water and canal seepage Canal NO 3 ) and chemical processes (litter and humus organic N decomposition; nitrification of NH 4 to NO 3 ; heterotrophic and autotrophic reduction of NO 3 , with the latter occurring in the presence of pyrite-bearing marine shale; and autotrophic reduction of O 2 , also occurring in the presence of shale) on NO 3 concentration in groundwater C NO 3 , NO 3 leaching from the unsaturated zone to the saturated zone of the aquifer, and NO 3 mass loading from the aquifer to the Arkansas River via groundwater discharge.The influence of each of the nine model factors was computed using the revised Morris method for sensitivity analysis, with results processed to determine parameter influence globally for the entire study region and specific to crop type, canal command area (i.e., the group of fields receiving irrigation water from a given canal), and individual grid cells.For the latter, spatial plots of sensitivity indices are presented to display the spatial distribution of influence for each model factor.
Results indicate that, generally, fertilizer loading, crop N uptake, and heterotrophic denitrification governed NO 3 mass transport, particularly in cultivated areas.However, their order of influence on C NO 3 and NO 3 mass leaching varies according to crop type and command area, and several command areas are influenced more, or at least to a significant degree, by nitrification, autotrophic denitrification, and Canal NO 3 .Spatial plots of cell-by-cell sensitivity indices further enhance the understanding of localized model factor influence, with each factor except for rate of heterotrophic O 2 reduction having the dominant influence over C NO 3 at one or more locations within the study area.Results also indicate that the concentration of O 2 in groundwater C O 2 is governed by rates of organic matter decomposition, which releases CO 2 and hence enhances heterotrophic reduction of O 2 .
In general, the procedure followed in this study provides key information regarding overall NO 3 fate and transport in an agricultural groundwater system, guidance for future data collection and monitoring programs, an indication of which parameters should be targeted during model parameter estimation, and guidance for implementing bestmanagement practices (BMPs) for NO 3 remediation, i.e., decreasing groundwater concentrations and mass loading to the stream network.For example, fertilizer loading, crop N up-take, and Canal NO 3 should be targeted in field data collection and observation, with Canal NO 3 monitored for each irrigation canal as often as possible, whereas first-order kinetic rate constants for nitrification, denitrification, and organic matter decomposition should be targeted during parameter estimation efforts.Furthermore, the procedure followed in this study also allows for data collection, management practice implementation, and parameter estimation to be performed on location-specific basis.For example, results suggest that a specific BMP (e.g., reduction in N fertilizer loading) may be optimal for several of the command areas but not for others, or that decreasing Canal NO 3 or the amount of NO 3 denitrified in shale outcrop locations will help remediate NO 3 only in a few specific locations within the study area.Also, data collecting points for specific model factors can be restricted to sub-region areas, either to a given command area or, with the use of the spatial plots of sensitivity indices, to even more localized sites.
The Supplement related to this article is available online at doi:10.5194/hess-19-4859-2015-supplement.

Figure 1 .
Figure 1.Location and hydrologic features of the study region in the lower Arkansas River valley in southeastern Colorado, showing the Arkansas River and tributaries (red), cultivated fields (yellow), irrigation canals (light blue), groundwater pumping wells (black dots), and the extent of near-surface shale (within 2 m of the ground surface) (green).

Figure 2 .
Figure 2. Features of the cultivation and data collection of the study region, including (a) canal command areas and location of groundwater observation wells, with a command area defined as the collection of fields receiving irrigation water from the same canal, and (b) the spatial distribution of crop cultivation during the 2006 growing season.

Figure 3 .
Figure 3. Depiction of the main processes simulated by the N reaction module of the UZF-RT3D model, with (a) conceptual model of the fate and transport of O 2 and N species in an irrigated soil-aquifer system wherein fertilizer, irrigation, and canal seepage bring solute mass into the subsurface environment, and (b) the annual cultivation schedule used in the N reaction module, including timing of planting, fertilizer loading, irrigation application, harvest, and plowing.NH 4 fertilizer has a split loading, with 40 % of the loading occurring 2 weeks before planting, and the remainder applied 6 weeks after planting.

Figure 4 .
Figure 4. (a) The spatial distribution of aquifer thickness (m) of the alluvium in the study region, (b) the finite-difference grid used in the calibrated and tested MODFLOW-UZF1 groundwater flow model, using 250 m by 250 m grid cells, (c) spatial distribution of soil water content simulated by the MODFLOW-UZF1 model, for June 2006, and (d) average-simulated water table elevation for the 1999-2009 time period.

Figure 5 .
Figure 5. Values of (a) fertilizer loading F NH 4 (kg ha −1 ) for corn, and (b) first-order rate constant of autotrophic denitrification λ auto NO 3

Figure 6 .
Figure 6.Summary of typical UZF-RT3D model results for the study region, showing spatial distribution of (a) C O 2 and (b) C NO 3 in shallow groundwater, and spatial distribution of mass loadings of nitrate to the Arkansas River system (main stem and tributaries) for (c) 2 December 2006, and (d) 10 August 2008, showing the contrast between the winter and summer seasons.

Figure 7 .
Figure 7. Spatiotemporal average value of C NO 3 in groundwater during the 2006-2009 simulation period for each canal command area for each of the 280 UZF-RT3D model simulations.The spatiotemporal average for the non-cultivated areas also is shown (small black crosses).

Figure 8 .
Figure 8. Global sensitivity plots (σ vs. µ * ) showing influence of the nine targeted model input factors on (a) C NO 3 in layer 1 of the model (top 0.5 m of the root zone), (b) C NO 3 in layer 4 of the model (shallow saturated zone of the aquifer), (c) NO 3 mass leaching from layer 3 to layer 4 (unsaturated zone to saturated zone), and (d) total mass loading of NO 3 from the aquifer to the Arkansas River.

Figure 9 .
Figure 9. Cell-by-cell (250 m by 250 m) plots of Savage scores for (a) F NH 4 , (b) Canal NO 3 , (c) λ nit , (d) N up , (e) λ het NO 3 , (f) λ auto NO 3 , (g) λ H , and (h) λ L , indicating the ranking of influence of that parameter on C NO 3 in groundwater for each of the 7776 cells in the study region.

Figure 10 .
Figure 10.Cell-by-cell (250 m by 250 m) plots of Savage scores for (a) Canal NO 3 , (b) F NH 4 , (c) λ H , (d) λ L , (e) N up , and (f) λ auto O 2 , indicating the ranking of influence of that parameter on C O 2 in groundwater for each of the 7776 cells in the study region.

Table 1 .
Baseline agricultural management and crop parameter values for the model simulations.

Table 2 .
Parameters and values for chemical reactions involving organic matter decomposition, dissolved oxygen, and nitrogen species for the baseline simulation model.
ter values to values found in the literature and from field data in the study area.The values of F NH 4 , λ auto NO 3 , and Canal NO 3 for each of the 280 simulations are shown in Fig.

Table 3 .
Sensitivity index (µ * ) for each of the model input factors investigated, indicating the degree of parameter influence on C NO 3 in the shallow saturated zone of the aquifer (in layer 4 of the grid) for the grid cells associated with each crop type and command area, with the values of the top three influential parameters for each crop type and command area in bold.

Table 4 .
Sensitivity index (µ * ) for each of the model input factors investigated, indicating the degree of parameter influence on NO 3 mass leaching from the shallow soil zone for the grid cells associated with each crop type and command area, with the values of the top three influential parameters for each crop type and command area in bold.