<emphasis role="strong">Title: Nitrate loads from land to stream are
balanced by in-stream nitrate uptake across seasons in a dryland
stream network</emphasis>
Amalia M. Handler, Ashley M. Helton, and Nancy B. Grimm
2 Materials and Methods
2.1 Site Description
Oak Creek is a tributary of the Verde River in the transition zone
between the Basin and Range Province and the Colorado Plateau in
Arizona, U.S. The creek originates at an elevation of 2,300 m in
ponderosa pine forest and descends through pinyon-juniper and high
desert ecosystems to its confluence with the Verde River at 950 m
elevation. Central Arizona has two rainy seasons separated by two
dry seasons. The winter-spring rainy season (January-April) is
characterized by frontal systems dominating precipitation in the
form of rain and snow. The dry summer season (May-July) can feature
daytime high temperatures that regularly exceed 36°C. The summer
rainy season (August-September) is characterized by convective
monsoon storms. The post-summer rainy season season
(October-December) is a dry period characterized by lower
temperature and less rainfall. Rainfall amounts and timing in both
rainy seasons exhibit high interannual variability.
Oak Creek drains a 1200 km2 watershed.
The main stem of Oak Creek is perennial, as are short sections of
three spring-fed tributaries. The remainder of the network is
temporary (including intermittent and ephemeral streams), with some
channels having seasonal surface water and many others supporting
surface flow only during large storms. The Oak Creek main stem and
perennial tributaries are bordered by a riparian gallery forest of
cottonwood and willow; however, the stream channel’s canopy is open
except for the smallest, high-elevation tributaries (most of which
are not perennial).
Land cover in the watershed is mostly undeveloped (>95%), with
some green space, residential, and commercial developed area (4%)
around the cities of Sedona and Cornville and a small proportion of
agricultural cover (<1%). There are two fish hatcheries in the
watershed that rely on springs for water supply and discharge water
to the main channel of Oak Creek (Oak Creek Watershed Council,
2012). Residential and agricultural areas are irrigated through a
combination of diverted stream water and groundwater wells.
To test our hypothesis that seasonal differences in precipitation
increases hydrologic connectivity that would change the balance of
lateral NO3
–
loading and in-stream
NO3
– uptake at the
network scale, we conducted synoptic surveys coupled with
NO3
– uptake
experiments. We then integrated these data in a mass-balance model
to scale in-stream
NO3
– uptake to the
entire stream network.
2.2 Synoptic Sampling
We conducted four seasonal sampling campaigns of stream water
chemistry and discharge across Oak Creek. Campaigns took place in
2017 during the winter wet season (February; minimum and maximum air
temperature for sample dates: 0 and 11°C), summer dry season (June;
20 and 36°C), summer wet season (September; 15 and 29°C), and winter
dry season (November; 10 and 23°C). During each campaign, we sampled
25–29 sites distributed across the watershed. Sites along the main
stem were selected to target locations above and below confluences
and irrigation divisions, and otherwise were evenly spaced. Both
perennial and intermittent (when flowing) tributaries as well as
irrigation ditches were sampled. Often, only one site along an
irrigation ditch was accessible. We identified and sampled two sites
that drain fish hatcheries, and noted a trout farm adjacent to the
main stem of Oak Creek in the upper portion of the watershed. The
main stem sites included locations below these potential point
sources. Most sites were accessed via public lands with a small
subset accessed via private land with permission. The campaign
during the winter wet season was the first completed and took place
during a period of high discharge relative to the latter campaigns.
As a result, many sampling locations were difficult or impossible to
access during this season. Subsequent campaigns were adjusted and
thus there are fewer sampling sites in the winter wet season in
common with those collected in other seasons.
To estimate stream discharge and surface area, we measured stream
wetted width and recorded 10-15 depth measurements at regular width
intervals at each site. We calculated cross-sectional area by
summing the area bounded by each depth measurement. For the main
channel, tributaries, and larger irrigation ditches, we measured
water velocity at each depth measurement. We calculated discharge as
the cross-sectional area times the mean water velocity. For smaller
irrigation ditches, we used pulse injection of sodium chloride
(NaCl) to measure discharge via dilution gauging (Day & Day,
1977).
At each site, we collected duplicate, field-filtered water samples
(0.7 µm GFF syringe filter, Fisher Scientific) for analysis of
NO3
–, ammonium
(NH4
+), and
chloride (Cl–). Chloride is a
non-biologically reactive ion that can indicate a change in water
source. Samples were stored in a cooler on dry ice until returned to
the lab where samples were stored at -20℃ until analysis. Nitrate
(NO3
– plus
NO2
–; hereafter
NO3
– in units of
µg N/L), NH4
+ (µg
N/L), and Cl– (mg/L) concentrations were
analyzed via Lachat QC 8000 flow injection analyzer (Lachat
Instruments, Loveland, CO). Ammonium was uniformly low across
seasons with approximately 20% of samples below the detection limit
(10 µg N/L) and 80% of samples below 20 µg N/L. As a result,
NH4
+ will not be
discussed further.
2.3 Nitrate Uptake Experiments
To determine the
NO3
– uptake rate
in Oak Creek and tributaries, we conducted seasonal
nutrient-spiraling experiments to pair with the sampling campaigns
(Newbold, Elwood, O'Neill, & Winkle, 1981). We conducted a total
of nine experiments by pulse injection (Covino, McGlynn, &
Baker, 2010). The experiments were conducted once per season on the
main stem of Oak Creek in 2017 in the spring (April), summer dry
season (June), summer rainy season (September), and winter dry
season (November). Discharge in the main stem was prohibitively high
earlier in the winter wet season (February and March) and only
dropped to a level that that allowed the experiment in April.
Experiments also took place once per season in Spring Creek, a
perennial tributary, in the winter wet, summer dry, summer wet, and
winter dry seasons. An additional experiment was conducted in Dry
Creek, an intermittent tributary, in the winter wet season. The
winter wet season was the only period for which Dry Creek had
flowing surface water during 2017.
Each experiment was conducted by dissolving NaCl (conservative
tracer) and sodium nitrate (NaNO3, reactive
tracer) into 10-50 L of stream water collected from the site. The
mass of NaCl and NaNO3 were adjusted such
that the solution would raise the
NO3
– concentration
by 50 µg N/L and chloride by 6.5 mg/L at the downstream collection
location when added in a single pulse. The reach lengths from the
injection point to the collection point were 710, 400, and 370 m for
Oak, Spring, and Dry Creeks, respectively. Reach lengths were
determined during pre-experiment tests to ensure that both tracers
had fully mixed across the stream width and depth. Following
collection of background water samples, the tracer-enriched water
was added to the stream instantaneously. Monitoring at the
collection location took place with an electrical conductivity meter
(either YSI 556 MPS, YSI 600 XLM, YSI DSM Pro, or Eureka Manta) and
a SUNA NO3
– sensor
(Sea-Bird Scientific, Bellevue, WA). The SUNA is factory calibrated
by Sea-Bird Scientific. For the Dry Creek experiment, due to SUNA
instrument failure, we instead collected 23 grab samples distributed
across the breakthrough curve at the collection location for
analysis of NO3
–.
An injection was considered complete when electrical conductivity
had returned to within 1 µS/cm of the background measurement or more
than three hours had elapsed since the injection time. In cases
where experiments were concluded after the full three hours,
electrical conductivity was 2-3 µS/cm above the pre-experiment
background. All grab samples were collected in duplicate in 500-mL
HPDE Nalgene® bottles; subsamples were poured into 50-mL centrifuge
tubes for transport to the lab on ice and stored at -20℃ until
analysis as described above.
Data from the experiments were processed using the TASCC method
described by Covino et al. (2010). Briefly, electrical conductivity
together with NO3
–
data were background-corrected, and interpolated to one-minute
intervals over the breakthrough curve. The breakthrough curve was
adjusted to have an equal number of observations for the rising and
falling limbs of the curve. The first-order uptake-rate coefficient
was calculated by taking the log of the background-corrected ratio
of
NO3
–:Cl–
of each observation divided by the same ratio in the injection
solution all divided by the distance between the injection and
collection points. Uptake length (Sw) was
calculated for each observation by taking the inverse of the uptake
coefficient. A regression between Sw and the
associated NO3
–
concentration was used to estimate Sw for the
stream at the background
NO3
– concentration
by extrapolation. Either a first-order (linear) or loss of
efficiency (log-transformed) regression was used according to
O'Brien, Dodds, Wilson, Murdock, and Eichmiller (2007), including
estimation of 95% confidence intervals for the uptake parameters at
the background
NO3
–
concentration. The vertical uptake velocity
(vf) was calculated by multiplying the
inverse of Sw by the reach discharge and
dividing by the mean reach width. We used a travel-time correction
to control for the differing amounts of time each sample had in the
stream prior to collection (Covino et al., 2010).
2.4 Stream Network Model
We used a steady-state mass-balance model to combine the field
measured in-stream
NO3
–
concentrations and
NO3
– uptake rates
at the network scale in order to evaluate how much in-stream uptake
attenuates lateral
NO3
– loads (Helton
et al., 2011; Mulholland & et al., 2008). Following Helton et
al. (2011), we implemented the model using an inverse approach to
estimate patterns of lateral
NO3
– loading from
land to stream reaches. The model estimates the lateral
NO3
– load
necessary to reproduce the observed spatial patterns in stream
NO3
– load, given
the field-based stream
NO3
– uptake rate,
discharge, and stream
NO3
–
concentration. The model accounts for both serial processing of
NO3
– along the
network and water diversions.
Discharge (Q, in volume/time) for each reach was calculated by
subtracting outgoing water from incoming water according to the
following equations
Qp = (∑Qp-1i +
QL) – (QW +
Qp+1i) Eq 1
where
QL = Ap ∙
Yp Eq 2
where Qp is the discharge in-stream reach
p, ∑Qp-1i is the sum of
the discharge of upstream reaches, p-1i,
contributing discharge to stream reach p,
QL is the lateral discharge from the adjacent
drainage area, QW is water loss from reach p,
Qp+1i is the discharge to the next downstream
reaches p+1i, Ap is
the area of the catchment draining directly to stream reach
p, and Yp is the per
unit subcatchment area lateral water yield to stream reach
p. We modified the existing model to include
the term QW, which represents the sum of
water loss to diversion or groundwater recharge.
For flowing reaches, we parameterized the lateral water yield
(Yp in Eq 2) for each season and stream
segment based on the discharge measured during the respective
synoptic sampling campaign. The change in discharge for each reach
was calculated as the discharge at the reach outlet minus the sum of
discharge from the upstream reach and any tributary confluences and
irrigation returns. Generally, only one site was accessible along
irrigation ditches and we assumed discharge to be uniform along
their length for the model. For network reaches containing an
irrigation diversion, the discharge in the ditch was subtracted from
the main channel. Each season included both net gaining and net
losing stream reaches. We parameterized lateral water yield as a
gross process where each stream reach has both lateral water input
and loss. Lateral water yield was calculated by dividing the change
in discharge by the lateral drainage area to the reach (Eq 2). We
implemented this by applying a lateral water yield to each
subcatchment that would produce the maximum discharge measured in
each season. From this, we calculated a stream water loss term
necessary to produce the field-measured discharge
QW (in Eq 1). Water loss encompasses stream
water transfer to the hyporheic zone and groundwater (Lange, 2005;
Valett, Fisher, Grimm, & Camill, 1994), mechanical withdrawal to
supply potable or irrigation water (Alger, Lane, & Neilson,
2021), and evapotranspiration (Dahm & Molles, 1992). By
implementing the lateral water flows as a gross rather than net
process, we allow for lateral inputs of
NO3
– to stream
reaches that have net discharge loss.
Stream NO3
– load
(N, in mass/time) is modeled similarly by
subtracting the outgoing
NO3
– load from
incoming load
Np = (∑Np-1i +
NL) – (NR +
NW + Np+1i) Eq 3
NL = Ap ∙
Lp Eq 4
where Np is the
NO3
– load in
stream reach p, ∑Np-1i
is the sum the of
NO3
– load from all
upstream reaches contributing
NO3
– to stream
reach p, NL is the
lateral NO3
– load
from the adjacent drainage area, NL is the
stream NO3
– load
taken up from reach p,
NW is
NO3
– loss due to
water loss from the reach, Np+1i is the
stream NO3
– load
to the next downstream reach p+1,
Ap is the area of the catchment draining
directly to stream reach p, and
Lp is the lateral
NO3
– load per unit
drainage area to stream reach p. Nitrate loss
stemming from water loss (Nw) is calculated
by multiplying the water loss by the
NO3
– concentration
in stream reach p. Lateral
NO3
– load per unit
drainage area is the primary model term of interest in this
analysis. Negative lateral
NO3
− loading
(Lp) signifies that there was not sufficient
uptake capacity in a stream reach or decrease from water loss to
achieve the measured
NO3
− concentration
(i.e., that NO3
−
concentration was lower than would be predicted based on upstream
concentration and uptake rate). For each stream reach, the mass of
NO3
– removed
(NR, in mass/time) is equal to the total
NO3
– load in the
stream reach times the fractional removal factor
NR = R ∙ Np Eq 5
The fractional removal factor is determined according to the
following equation from Wollheim, Vörösmarty, Peterson, Seitzinger,
and Hopkinson (2006)
R = 1 – e(-vf/HL) Eq 6
HL = Qp /
SAp Eq 7
and vf is the experimentally determined
vertical uptake rate for
NO3
– (in m/d).
We parameterized in-stream
NO3
– uptake in the
model based on field data. The vertical uptake rate
(vf in Eq 6) was determined based on the
median of all stream
NO3
– uptake rates
measured in the main channel and tributaries (Sec 2.3) because
vf did not vary significantly based on season
or stream type. Initial model testing with the minimum and maximum
NO3
– uptake values
measured experimentally changed only the magnitude of the lateral
NO3
– load, but
patterns and overall stream network
NO3
– attenuation
remained high. The vf is normalized by
hydraulic load (HL in length/time), which is
a measure of the rate of water passage through the stream relative
to the benthic surface area. Hydraulic load is calculated by
dividing the discharge (Qp in volume/time) by
the surface area (SAp, calculated as stream
length times average width). Average stream width
(w) is estimated as
w = aQb Eq 8
where a and b are the
width coefficient, which controls the scaling, and width exponent,
which controls the rate of increase, respectively (Leopold &
Maddock, 1953). Both a and
b were determined from field survey data from
each season.
We implemented four steady-state models, one for each synoptic
sampling campaign that represents the dynamics of each season:
Winter wet, summer dry, summer wet, and winter dry. We derived
network structure from the National Hydrography Dataset Plus Version
2 (NHDPlusV2) flowlines, flow direction, and drainage area. This
included the two largest irrigation ditches in the watershed that
have both diversions from and returns to the main channel of Oak
Creek. We derived subcatchments for each synoptic sampling location
for each season by delineating the drainage area for each reach
between sampling locations. Reaches were further subdivided into
segments between network junctions (e.g., tributary confluences,
irrigation diversions, irrigation returns) or 1000 m, whichever was
shorter.
The extent of the flowing portion of the drainage network varied
between campaigns. For the summer dry, summer wet, and winter dry
seasons, the extent of the flowing network followed the perennial
flow designations in the NHDPlusV2. For the winter wet season, one
large tributary and two headwater streams had flowing surface water
that were coded as intermittent or ephemeral in the NHDPlusV2.
Across surveys, flowlines that did not have flowing surface water
were parameterized to have negligible lateral water yield. We
estimated lateral
NO3
– loading
(Lp) by using a model-independent parameter
estimator (PEST 16.0, Model-Independent Parameter Estimation &
Uncertainty Analysis). PEST adjusts model parameters, based on their
values and model outputs, compared to observed data. In this
instance, PEST was used to adjust the lateral
NO3
– loads to
reproduce the measured in-stream
NO3
– load observed
at each sampling location. The lateral
NO3
– load lack
error estimates because these would be based largely on error
associated with the in-stream
NO3
– uptake.
Initial model testing revealed that varying the stream
NO3
– uptake
changed only the magnitude of the lateral
NO3
– load, but
patterns in lateral
NO3
– load and
overall stream network
NO3
– attenuation
remained similar.
2.5 Land-cover data
To test whether lateral
NO3
– loads are
associated with land cover, we used the 2011 National Land Cover
Database (NLCD) (Dewitz, 2014). We calculated zonal statistics for
land-cover types at two scales: (1) The subcatchment area and (2) a
200-m buffer area extending perpendicular to the direction of flow
(100 m on both sides of the stream). For each subcatchment and
buffer area, we calculated the proportions of total agricultural
land (sum of NLCD cropland and hay/pasture), developed land (sum of
NLCD developed open spaces and low-, medium-, and high-intensity
development), and wetland area (sum of NLCD herbaceous and woody
wetlands). The 2011 NLCD was the best available data at the time of
the analysis; however, we reviewed changes in the percent cover of
the classes included in this study between the 2011 and 2016 NLCD
and found that most changes were negligible to small (<2.5%).
2.6 Data analysis
We evaluated seasonal differences and spatial patterns within Oak
Creek watershed using the measured discharge, stream
NO3
– and
Cl– concentrations,
NO3
– uptake
experiments, and hydrologic and lateral
NO3
– loading model
outputs. We evaluated the fluvial network patterns in stream
discharge, NO3
–
concentration, and Cl– concentration by
plotting the data against distance from the watershed outlet, and
plotted the lateral
NO3
– loads in map
format. We tested for seasonal differences in discharge, chemistry,
NO3
– uptake, and
lateral NO3
– loads
using a one-way analysis of variance (ANOVA). In addition, we tested
for differences in-stream chemistry among site types (i.e., main
channel, tributary, irrigation ditch) using ANOVA. We varied the
base group in the ANOVA to evaluate pairwise differences between
groups. We evaluated the ratio of
NO3
– to
Cl− to identify outliers that may
indicate a difference in source water. We evaluated differences in
spatial patterns of discharge and
NO3
– and
Cl– concentrations from the synoptic
campaigns by calculating Spearman’s correlation coefficient between
data for the same sampling locations in different seasons, as a
measure of synchrony. We expected that seasons would either (1) have
a similar spatial pattern in sources and sinks for discharge,
NO3
–
concentration, and Cl– concentration and
therefore exhibit synchrony (i.e., positive pairwise correlations
between seasons) or (2) would have a different pattern and thus be
asynchronous (i.e. negative or no pairwise correlation between
seasons), potentially indicating a different set of sources and
sinks. An assumption of this analysis is that the
NO3
– and
Cl– concentrations in the sources do not
change much across seasons, which is supported based on findings
from another desert stream showing consistency of source chemistry
across seasons (Dent & Grimm, 1999) and years (Dong, Ruhí, &
Grimm, 2017). We completed a similar correlation analysis for the
lateral NO3
– loads
from the model, but applied an additional constraint that the
subcatchment area needed to cover identical areas among seasons for
pairwise comparison. In both correlation analyses, the number of
sites included differed in each pairwise comparison because the
sample locations differed between the campaigns. In particular, the
winter wet season was the first and took place during a period with
high discharge, making sampling difficult or impossible in many
locations. As a result, the winter wet season had comparatively
fewer sites and subcatchments that matched the summer dry, summer
wet, and winter dry seasons. Finally, we evaluated the relationship
between lateral
NO3
– loading and
land cover variables as well as stream characteristics (e.g.,
discharge, width, depth, and chemistry) using correlation analysis.
Correlations with land-cover variables were evaluated at the
subcatchment scale and for a 200-m buffer area centered on the
stream channel. We evaluated the correlations at these two scales
because the subcatchments could be very large, and land cover
several kilometers from the stream may be hydrologically
disconnected from the stream channel.
In all cases, correlations were calculated using Spearman’s rank
correlation coefficient (rho) because of non-linear relationships
among variables. For correlations between
NO3
– loading and
land cover variables, we excluded sites that we suspected were
influenced by point sources of nutrients, including the two fish
hatcheries and one trout farm. At these locations, we reasoned that
NO3
– loading would
be influenced more by the point source than the land cover variables
included in the analysis. All statistical analyses were performed in
R version 4.3.2 (2023-10-31 ucrt) (R Core Team, 2023).
References
Alger, M., Lane, B. A., & Neilson, B. T. (2021). Combined
influences of irrigation diversions and associated subsurface return
flows on river temperature in a semi-arid region.
Hydrological Processes, 35(8), e14283-e14283.
doi:https://doi.org/10.1002/hyp.14283
Covino, T. P., McGlynn, B. L., & Baker, M. A. (2010). Separating
physical and biological nutrient retention and quantifying uptake
kinetics from ambient to saturation in successive mountain stream
reaches. Journal of Geophysical Research: Biogeosciences,
115(4), 1-17. doi:10.1029/2009JG001263
Dahm, C. N., & Molles, M. C. (1992). Streams in Semiarid Regions
as Sensitive Indicators of Global Climate Change. In S. G. Fisher
& P. Firth (Eds.), (pp. 250-260). New York, NY, USA:
Springer-Verlag.
Day, T. J., & Day, T. T. (1977). Field procedures and evaluation
of a slug dilution gauging method in mountain streams.
Journal of Hydrology (New Zealand), 16(2),
113-133. Retrieved from http://www.jstor.org/stable/43944411
Dent, C. L., & Grimm, N. B. (1999). Spatial heterogeneity of
stream water nutrient concentrations over successional time.
Ecology, 80(7), 2283-2298.
doi:https://doi.org/10.1890/0012-9658(1999)080[2283:SHOSWN]2.0.CO;2
Dewitz, J. (2014). National Land Cover Database (NLCD)
2011 Land Cover Conterminous United States.
Dong, X., Ruhí, A., & Grimm, N. B. (2017). Evidence for
self-organization in determining spatial patterns of stream
nutrients, despite primacy of the geomorphic template.
114(24), E4744-E4752.
doi:doi:10.1073/pnas.1617571114
Helton, A. M., Poole, G. C., Meyer, J. L., Wollheim, W. M.,
Peterson, B. J., Mulholland, P. J., . . . Zeglin, L. H. (2011).
Thinking outside the channel: Modeling nitrogen cycling in networked
river ecosystems. Frontiers in Ecology and the
Environment, 9(4), 229-238. doi:10.1890/080211
Lange, J. (2005). Dynamics of transmission losses in a large arid
stream channel. Journal of Hydrology, 306(1),
112-126. doi:https://doi.org/10.1016/j.jhydrol.2004.09.016
Leopold, L. B., & Maddock, T. (1953). The hydraulic
geometry of stream channels and some physiograhpic
implications. Retrieved from
Mulholland, P. J., & et al. (2008). Supplementary information:
Stream denitrification across biomes and its response to
anthropogenic nitrate loading. Nature, 254.
doi:10.1038/nature06686
Newbold, J. D., Elwood, J. W., O'Neill, R. V., & Winkle, W. V.
(1981). Measuring nutrient spiralling in streams. Canadian
Journal of Fisheries and Aquatic Sciences, 38, 860-863.
O'Brien, J. M., Dodds, W. K., Wilson, K. C., Murdock, J. N., &
Eichmiller, J. (2007). The saturation of N cycling in Central Plains
streams: 15N experiments across a broad gradient of nitrate
concentrations. Biogeochemistry, 84(1), 31-49.
doi:10.1007/s10533-007-9073-7
Oak Creek Watershed Council. (2012). Improvement plan for
the Oak Creek Watershed, Arizona. Retrieved from
R Core Team. (2023). R: A Language and Environment for Statistical
Computing. Vienna, Austria: R Foundation for Statistical Computing.
Retrieved from https://www.R-project.org/
Valett, H. M., Fisher, S. G., Grimm, N. B., & Camill, P. (1994).
Vertical Hydrologic Exchange and Ecological Stability of a Desert
Stream Ecosystem. Ecology, 75(2), 548-560.
Wollheim, W. M., Vörösmarty, C. J., Peterson, B. J., Seitzinger, S.
P., & Hopkinson, C. S. (2006). Relationship between river size
and nutrient removal. Geophysical Research Letters,
33(6), 2-5. doi:10.1029/2006GL025845
