Salt Flushing, Salt Storage, and Controls on Selenium and Uranium: A 31‐Year Mass‐Balance Analysis of an Irrigated, Semiarid Valley

Salinity, selenium, and uranium pose water‐quality challenges for the Arkansas River in southeastern Colorado and other rivers that support irrigation in semiarid regions. This study used 31 years of continuous discharge and specific conductance (SC) monitoring data to assess interannual patterns in water quality using mass balance on a 120‐km reach of river. Discrete sampling data were used to link the SC records to salinity, selenium, and uranium. Several important patterns emerged. Consumptive use reduced discharge by a median value of 33% and drove corresponding increases in salinity and uranium concentrations. Increased water availability for irrigation from rainfall and upstream snowpack in 1995–1999 flushed additional salinity and uranium into the river in 1996–2000; average annual total dissolved solids (salinity) concentrations increased 25%, and loads increased 131%. Smaller flushing events have occurred, sometimes lagging an increase in water availability by about one year. The pattern indicates flushing of salts temporarily stored, evaporatively concentrated, or of geologic origin. Mobilization of selenium from the reach was minor compared to salinity and uranium, and net selenium removal from the river was suggested in some years. Several processes related to irrigation could be removing selenium. The results provide context for efforts to improve water quality in the Arkansas River and rivers in other semiarid regions.


INTRODUCTION
Rivers in arid and semiarid settings worldwide are harnessed for irrigation to boost crop productivity from fertile soils where rainfall is limiting, but water quality is a common challenge in such settings as evaporation and transpiration drive less water to hold generally the same quantity of both natural and anthropogenic salts (Hillel 2000;Duncan et al. 2008;Szynkiewicz et al. 2011;Neissi et al. 2019).
Additionally, concentrations of problematic trace elements such as Se, As, V, Mo, and U can become elevated by the same processes (Bradford et al. 1990;Seiler et al. 1999;Xie et al. 2015). Managing salinity and problematic constituents is crucial for water users both within the irrigated landscape and downstream. Such management is best achieved when the controls on salinity are understood, and there have been successes in improving water quality in irrigated regions (Duncan et al. 2008;Herrero and Castañeda 2018). Controls on salinity plus Se and U for a reach of the Arkansas River in southeastern Colorado are the focus of the present paper.
Water quality has been a challenge in the Arkansas River Valley since the earliest studies (Gilbert 1896), and salinity is the issue that has been longest recognized (Miles 1977). Reductions in crop yields caused by salinity in Arkansas River water are estimated to range 6-17% along different sections in Colorado (Burkhalter and Gates 2005;Morway and Gates 2012) and are presumably greater downstream in Kansas (Whittemore 2000). More recently, concentrations of Se and U have gained attention as problematic (Zielinski et al. 1995;Gates et al. 2009;Miller et al. 2010). Selenium in surface water exceeds the Colorado chronic wildlife toxicity standard of 4.6 µg/L in many locations (CDPHE 2016), as well as the 3.1 µg/L lotic water standard recently recommended (USEPA 2016). Uranium exceeds the 30 µg/L standard for public drinking-water supplies (USEPA 2018) particularly in downstream regions. Not just surface water users are affected. As pumping has lowered water levels in the High Plains aquifer in Kansas, Arkansas River waters now recharge and mix with the higher quality waters of the aquifer (Whittemore 2002).
The causes of water-quality challenges in the Arkansas River Valley have been generally attributed to interactions between irrigated agriculture, a semiarid climate, and Cretaceous marine shales, limestones, and other sedimentary rocks that weather to produce salts, Se and U (Zielinski et al. 1995;Sutherland 2008;Gates et al. 2009;Miller et al. 2010;Bern and Stogner 2017). Contributions of U from fertilizers are minimal compared to natural sources (Zielinski et al. 1997). Large quantities of water are diverted from the Arkansas River for irrigation. Evaporation and transpiration consume water and thus concentrate salts, Se, and U already dissolved in the water. Water leaking from irrigation structures or applied in excess of crop demand becomes deep percolation. In the subsurface, deep percolation mobilizes additional salts, Se and U derived from marine sedimentary rocks and weathered materials within canal command areas before rejoining the river as subsurface return flow. The process is repeated by subsequent diversions farther downstream. Efforts to improve water quality face both the constraints of the agricultural-hydrologic system, along with socioeconomic and legal challenges (Sutherland 2008;Triana et al. 2010;Bailey et al. 2015;Sharp et al. 2016). In particular, reducing deep percolation through improved irrigation efficiency has been difficult to achieve because downstream users have rights to the return flows generated by inefficiency (Triana et al. 2010;Sharp et al. 2016).
The conceptual model above generally explains the spatial pattern of substantial increases in salinity, Se, and U concentrations with distance downstream, along with general decreases in streamflow (Cain 1987;Miller et al. 2010). Detailed understanding of the controls on water quality has generally been pursued by multiple efforts to model interactions between surface water and groundwater in combination with water quality. Earlier and simpler models focused primarily on evapotranspiration and return-flow effects on salinity (Konikow and Bredehoeft 1974;Goff et al. 1998;Gates et al. 2002;Burkhalter and Gates 2005). Conclusions drawn from some models based upon short-term data have been found incorrect (Konikow and Person 1985). Subsequent efforts incorporated reactive transport of carbonate and sulfate species (Lin and Garcia 2008). Most recently, oxidation-reduction reactions have been incorporated to account for mobilization and immobilization of Se along with sulfate, the dominant component of salinity (Bailey et al. 2014;Tavakoli-Kivi and Bailey 2017). This recent work has emphasized the influence of nitrate, dominantly sourced from fertilizer, as an influence on oxidation-reduction reactions (Bailey et al. 2012). Generally, all such modeling efforts have used four years or less of field data to calibrate or evaluate key parameters.
In this paper, a mass-balance analysis of salts (as total dissolved solids, TDS), Se, and U in surface water is conducted for a reach of the Arkansas River in southeastern Colorado that supports substantial irrigated agriculture. Volumes of water (as discharge), concentrations, and loads of the above constituents are calculated on an annual basis over 31 years from 1988 to 2018. One emphasis is therefore on interannual variability and covers a longer period than previous studies, but at a coarser resolution. Another emphasis is consistency of patterns. Changes in factors like surface-water or groundwater storage can be substantial seasonally and from year to year, but such changes cancel out over the 31 years. In this complex and highly managed reach of river, where water is diverted from the river and some fraction is returned multiple times, straightforward calculations on long time steps can reveal the broadest controls on water quality. By examining the quantities of water and dissolved constituents entering and exiting the study reach, along with their timing, new insights into the controls on salinity, Se, and U in the Arkansas River are achieved, and a pattern of salt mobilization is recognized that may exist in other settings.

STUDY AREA
The study area is an approximately 120 km reach of the Arkansas River extending from Avondale, Colorado, at the upstream end to Las Animas, Colorado downstream ( Figure 1). The within-reach contributing area is composed of five hydrologic unit code subbasins (HUC8) (Seaber et al. 1987) and is in the Great Plains physiographic region. Snowmelt, from the Rocky Mountain portion of the basin upstream to the west supplies most of the flow in the river, with lesser contributions from precipitation falling on tributary watersheds of the study reach. Flows are augmented by transmountain diversions from outside the watershed in the mountains, primarily the Fryingpan-Arkansas Project (Abbot 1985). Flows in the study reach are partially controlled by numerous upstream reservoirs, most notably Pueblo Reservoir, which was completed in 1974. Dissolved solids entering the study reach via the main-stem river were sourced 67% from the main stem of the river and 33% from Fountain Creek (Figure 1b), not far upstream from Avondale, during 2000-2006(Miller et al. 2010. Climate within the reach is semiarid, with mean monthly temperatures ranging from À1°C in winter to 25°C in summer; mean annual precipitation is around 320 mm and the majority falls during summer months (https://www.ncdc.noaa.gov/cdo-web/). Irrigated agriculture covers approximately 64,000 ha. Groundwater pumping from the alluvial valley aquifer supports some irrigation, but pumping has been substantially reduced since a 1985 lawsuit (Kansas v. Colorado 1995). Primarily, irrigation is supported by substantial diversions from the Arkansas River into several unlined canals and lateral distribution systems that parallel the river for many tens of kilometers and were originally constructed between 1874 and 1890 (Sutherland 2008). Seepage from these systems can be considerable (Dash 1994). The result is a rather dynamic shallow groundwater system that is strongly influenced by substantial quantities of canal seepage and irrigation-induced deep percolation. Groundwater travel times of 0-3 years to surface have been modeled within large portions of the canal command areas, but areas also exist where modeled travel time is >20 years (Bailey et al. 2015). Much of the groundwater pumped for irrigation within the study area can thus be understood to have been surface water from a canal diversion in the relatively recent past. Downstream, the river is dominated by canal seepage and irrigation return flows for parts of the year, and it is possible for return flow from one diversion to be diverted again downstream (Cain 1985). Six off-channel reservoirs in the study area, all north of the river, have a combined live storage of approximately 192,000 ML. Additionally, Cheraw Lake exists in a topographic low in the landscape with restricted surface drainage. Groundwater, surface water, and return flows collect in Cheraw Lake and become highly saline because of evaporation.
Furrow irrigation is the most common irrigation method in the study area, which can generate substantial deep percolation, but some areas of sprinkler and drip irrigation also exist. More widespread conversion to sprinkler irrigation has occurred downstream from and outside of the study area. Popular crops include corn, alfalfa, and melons, with lesser emphasis on vegetable crops. Cultivated soils are generally deep, well drained, and are derived from alluvium, eolian deposits, and residuum from weathering of Cretaceous marine sedimentary rocks that underlie most of the study reach (Soil Survey Staff 2019). Feedlots and wastewater treatment plants for small municipalities are present, but have minor effects on water quality. Agriculture is well-established in the study area and land-use changes between 1988 and 2018 have been relatively minor.

Synoptic Sampling
Six synoptic samplings at the sites shown in Figure 1b were completed to gain perspective on seasonal and spatial changes in water quantity and quality. Samplings were conducted in February, May, August, and November of 2017, and February and May of 2018 over 2-3 days per sampling. Rain events during sampling in August 2017 had a notable influence on results.
Synoptic samples were collected using the depthintegrated, multiple verticals compositing technique from wadable waters and multiple composited grab samples from higher flows. Sampling techniques and field measurements of specific conductance (SC) followed methods described in the United States (U.S.) Geological Survey (USGS) National Field Manual (USGS, variously dated). Water-quality samples were filtered to <0.45 lm by capsule filter and preserved by acidification to pH < 2 using nitric acid. Unfiltered splits were not collected because previous work on the Arkansas River and tributaries has shown that total concentrations (unfiltered samples) of Se and U are only slightly different (À0.2 to 0.2 and À2.2 to 0.7 lg/L, respectively) from dissolved concentrations (filtered samples) across broad concentration ranges (0.06-15.8 and 3.4-44.5 lg/L, respectively; Bern and Stogner 2017). Dissolved Se and U concentrations were determined by AGAT Laboratories (Calgary, Alberta, Canada) using inductively coupled plasmamass spectrometry (Lamothe et al. 2002). Blind standards were submitted with samples and results considered satisfactory if within 15% of expected values. Flow measurements were made by standard methods (Turnipseed and Sauer 2010). All water quality-data used in this study are archived in the National Water Information System (NWIS; https://doi.org/10.5066/ F7P55KJN) and can be accessed using the site numbers (e.g., Figure 1b).

Isotope Data and Evaporation Estimates
To quantify the contribution of evaporation to consumption of flows within the study reach, hydrogen and oxygen isotope ratios of river water were measured. Evaporation favors the lighter isotopes of oxygen and hydrogen in water and causes heavy isotope enrichment in the water left behind (Skrzypek et al. 2015). In contrast, transpiration does not cause isotopic enrichment of water. Synoptic sampling results from the Arkansas River near Avondale, Colorado (station 07109500) were used to characterize the isotopic composition of water input to the study reach ( Figure 1b). To characterize evaporated output water, monthly samples were collected at the Arkansas River at Las Animas, Colorado (station 07124000) from February 2017 through May 2018. Water isotope ratios of hydrogen and oxygen (dD and d 18 O) were measured at the USGS Reston Stable Isotope Laboratory in Reston, Virginia (R esv esz and Coplen 2008a, b).
The Hydrocalculator software (Skrzypek et al. 2015) was used to calculate the fraction of water evaporated. The nonsteady-state model was used, treating starting water as a volume that becomes progressively more fractionated as it moves downstream. Starting water was represented by a sample from the Arkansas River near Avondale, Colorado during spring snowmelt in May 2017 (dD = À106&, d 18 O = À13.91&), as spring runoff dominates discharge volumes. Hydrocalculator requires inputs of average temperature, relative humidity, and isotopic composition of precipitation. Average temperature for the Rocky Ford station of the CoAgMET network (https://coagmet.colostate.edu/) was 13°C for January-December of 2017 and relative humidity was 60%. The isotopic composition of precipitation was represented by the weighted means (dD = À74&, d 18 O = À10.6&) from northeastern Colorado, also in the Great Plains (Harvey 2005).

Continuous Monitoring Records
Continuous records of discharge and SC from streamgages at locations representing inputs and outputs from the study reach ( Figure 1b) spanning calendar years 1988 to 2018 underpin the mass-balance analysis. USGS surface-water data are hosted by the NWIS database (https://doi.org/10.5066/F7P55KJN). The two input locations to the study reach are the Arkansas River near Avondale, Colorado (station 07109500) and the Huerfano River near Boone, Colorado (07116500) (Figure 1b). The Huerfano River is the only tributary considered as an input because it is the only tributary with a substantial catchment area in the mountains that generates substantial runoff. Other tributaries in the study reach such as the Apishipa River and Timpas Creek have catchments primarily on the semiarid plains ( Figure 1). Average recharge in such nonirrigated, plains catchments is estimated to be 2-10 mm/yr (Wolock 2003) and such tributaries are generally intermittent in upstream regions outside the canal command areas, yet flow year-round inside the canal command areas. Consequently, much of the discharge from such plains tributaries is simply irrigation return flow from upstream Arkansas River diversions, or irrigation return flow from wells that tap the shallow aquifer that is influenced by canal leakage and irrigationinduced, deep percolation (Cain 1985;Dash 1994;Figure 2a). The Huerfano River may capture some return flow from the Bessemer Ditch, but that diversion is upstream from the study reach. Regardless, the Huerfano River is a small input to the study reach, averaging 3% of the annual discharge of the Arkansas River at Avondale. The three output locations are the Arkansas River at Las Animas, Colorado (07124000), the Fort Lyon Canal (Colorado State Division of Water Resources [DWR] station FLY-CANCO), and the Fort Lyon Storage Canal (DWR station FLSCANCO). The Fort Lyon Canal is included because it carries substantial flows past the Las Animas streamgage, and the Fort Lyon Storage Canal is included because it carries flows to off-channel reservoirs that can supply the Fort Lyon Canal ( Figure 1). Discharge data for the diversion point of both canals are collected and hosted by the DWR (https://www.dwr.state.co.us).
SC is not measured at either the Fort Lyon Canal or the Fort Lyon Storage Canal streamgages. To represent SC at the Fort Lyon Storage Canal streamgage, the record of continuously monitored SC from the Arkansas River at Catlin Dam near Fowler, Colorado (station 0711970) is used. The streamgage at Catlin Dam is 9 km upstream of the Fort Lyon Storage Canal diversion, with no major intervening tributaries. Two methods are employed to represent SC in the Fort Lyon Canal and the results are compared in the mass-balance analysis. One method uses continuous SC data from the Arkansas River near Rocky Ford (07120500) streamgage; this method underestimates SC in the Fort Lyon Canal. The Rocky Ford streamgage is~12 km upstream from the Fort Lyon JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION diversion point, and the confluence with Timpas Creek lies between the streamgage and the canal diversion, which could affect estimates of SC for the canal. Four measurements in the Fort Lyon Canal (380646103125701) in 2017 and 2018 showed canal SC that averaged 12% greater (range 9-17%) than at the Rocky Ford streamgage. The second method uses continuous SC data from the Las Animas streamgage; this overestimates SC in the canal. The Las Animas streamgage is~38 km downstream from the diversion, with two other tributaries in between, and using the same four measurements, canal SC averages 35% less than at the Las Animas (range 22-50%). The underestimate and overestimate provide a range for consideration, but complete calendar year continuous SC data at the Rocky Ford streamgage are only available for 2013-2018.
Discharge and SC values used from all streamgages were daily means, and some data required processing to remove discharge outliers and fill SC data gaps. Anomalous discharge values more than 29 greater than maxima in other years were removed from the Fort Lyon Canal record and were the only outliers. The SC data from all streamgages contained gaps because of equipment problems or where data had not met quality-control criteria. Gaps are unacceptable for mass-balance assessment because they compute as zero load for the missing days. Therefore, daily mean SC values were interpolated to fill the gaps. For gaps of three days or less, the single values before and after the gap were averaged to represent SC. For gaps of four days or more, or when discharge was trending across the gap, the power law relation between daily-mean discharge and daily-mean SC (e.g., Bern et al. 2015) was established using the two days preceding the gap and the two days after the gap. The power-law relation was then used to estimate daily mean SC from the daily-mean discharge values. Interpolation across these gaps adds uncertainty to the analysis that is not accounted for, but the gaps average only 3.4, 3.3, 3.5, and 4.0 days for the Avondale, Catlin Dam, Rocky Ford, and Las Animas SC records, respectively, and constituted 5%, 5%, 12%, and 11% of the respective SC records. No continuous SC data were available for the Huerfano River, but 323 discrete SC measurements were made between 1988 and 2018. The R package EGRET was used to generate daily estimates of SC using the Weighted Regressions on Time Discharge and Season methodology (Hirsch et al. 2010;Hirsch and De Cicco 2014). Daily-mean discharge values were summed to give total annual discharge for each calendar year. Daily-mean SC values were weighted by daily-mean discharge to calculate annual mean SC for the calendar year.
Archival, TDS concentrations, determined gravimetrically by evaporation at 180°C, exist for many samples collected at the USGS streamgages from which SC data were obtained, as do dissolved Se and U concentrations (https://doi.org/10.5066/F7P55KJN). The unique linear relation between SC and the constituents at each streamgage was established by regression, with standard errors as a measure of uncertainty. The data are plotted in Figures S1 and S2. Linear-regression equations were used to convert annual discharge-weighted, mean SC to annual, FIGURE 2. Conceptual model depicting (a) discharge inputs and outputs for the study reach considered by the mass-balance calculations prior to simplification and (b) load inputs and outputs as implemented by the mass-balance calculations.
discharge-weighted, mean constituent concentrations. Such concentrations were multiplied by discharge to calculate annual loads. For Se and U, data for the regressions were restricted to samples collected and analyzed beginning in 2006 up through the 2018 synoptic samples and total numbers of samples are noted in Tables 1 and 2. All regressions were conducted in R (v 3.5.1) (R Foundation for Statistical Computing, Vienna, Austria).

Mass-Balance Analysis
The mass-balance analysis not only focuses upon the area along the study reach bounded by the unlined canals but also includes the shallow aquifer outside the influence of canals and irrigation (Figure 1b). Mass-balance analysis ( Figure 2a) begins with assessment of flow volumes. Complete water mass balance for the reach would be: where Q O are flows out of the reach, Q I are flows into the reach, Q R is natural groundwater recharge from precipitation, Q CON is consumption by evapotranspiration, DS SW is change in surface-water storage, DS GW is change in groundwater storage, and Q GW-net is net inflow from groundwater outside the irrigated area of the reach and would be negative for a net gain. Data are readily available to address only some of these components, and other components are omitted after consideration of the potential error they would introduce considering the annual time step and a 31-year period of analysis. Discharge outputs (Q O ) are the sum of the Arkansas River at Las Animas (Q LA ), Fort Lyon Storage Canal (Q FLS ), and Fort Lyon Canal (Q FL ). Discharge inputs (Q I ) are the sum of the Arkansas River near Avondale (Q Av ) and the Huerfano River near Boone (Q H ). Natural groundwater recharge (Q R ) is the first component to be omitted for simplicity, due to its generally small contribution (2-10 mm/yr; Wolock 2003). Consumption (Q CON ) will be solved for later.
Three other terms (DS SW , DS GW, and Q GW-net ) are omitted from the calculations because they are both challenging to quantify with acceptable accuracy over 31 years, and because their magnitudes are relatively small compared to the nonomitted terms on an annual time step and more so for totals from the 31year period of analysis. Although these components are omitted as individual terms in the calculations they are important. However, their influence must be inferred from the outcome of the simplified mass balance analysis because by omitting them from the calculation they become implicitly included in the term to be solved for, Q CON . As detailed in the following section, however, their influence is smaller than the input term (Q I ) and they are thus judged to have little influence. Regarding DS SW , at the calendar year break on December 31, surface-water storage is not substantially changing during this low flow period. Except for the Fort Lyon Storage Canal, other canals are generally filled only during the irrigation season from March 15 to November 15 and streams and drains are no longer carrying tailwater during the nonirrigation season. Off-channel reservoir capacity is 192,000 ML, or 24% of average input flows, and thus interannual changes in storage can drive some error in individual years, but storage changes from maximum to minimum capacity are usually spread over multiple years. Additionally, only net changes in off-channel storage between 1988 and 2018 contribute to net error in the analysis and total off-channel storage capacity is equivalent to only 0.75% of the total inflow to the study area over 31 years.
Groundwater storage (DS GW ) within the study reach changes both seasonally and interannually as evidenced by fluctuations in the depth to water below land surface ( Figure S3). To assess how such changes could affect the mass-balance analysis, a 123,000 ha alluvial aquifer in the study area with a specific yield of 0.1 (Morway et al. 2013) may be considered. A 1-m change in groundwater level would reflect 123,000 ML of water, which is 15% of mean annual surface water inflow to the study area as calculated in this study, but only 0.48% of the total inflow over 31 years. Thus, large changes in DS GW could substantially influence

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
the analysis in individual years, but the net effect on the 31-year analysis is negligible, particularly as average net change in depth to water for 11 wells in the study area was 0.13 m over the study period (Figure S3). Inflow (Q GW-net ) of groundwater from the larger contributing area ( Figure 1a) outside the irrigated area of the reach (Figure 1b) is likely to be volumetrically small, because of the limited recharge, but contributions of TDS, Se, and U from such highly evaporated waters could be substantial. By omitting DS SW , DS GW , and Q GW-net , a simplified mass balance for water volumes that still provides substantial insights is possible. As detailed in the following section, the influence of these three terms on TDS, Se, and U mass balance are accounted for by a within-reach loading term. With four components omitted, the volumetric mass balance (Equation 1) is reduced to: To solve for the unknown Q CON , Equation (2) rearranges to: Note that Q CON should not be confused with actual evapotranspiration within the irrigated area, which would also include consumption of precipitation. However, Q CON is a measure of net consumption of river water that serves as an informative reference point in subsequent interpretations. All flows are considered in units of ML/yr.
Another reference point can be calculated from Q CON . Using TDS as an example constituent, the expected concentration of TDS in water flowing out of the study reach at the Arkansas River at Las Animas (C Con TDS;LA ) if there was no within-reach loading, only water consumption by evapotranspiration and conservation of TDS in remaining water, is given as follows: Here, the input load (L I ) is the combined TDS load from the Arkansas River near Avondale (L TDS,Av ) and Huerfano River near Boone (L TDS,H ) ( Figure 2b). The C Con TDS;LA reference point provides insight into how consumption alone can drive water-quality changes. Comparison of C Con TDS;LA to annual average concentrations of TDS in outflow from the study area helps illustrate how much additional change comes from within-reach loading.
The mass-balance analysis of TDS, Se, and U is based upon the load estimations from the previous section and similar simplified equations to those used for water volumes (Figure 2b). The analysis is done with downstream load output defined in two ways for comparison purposes. In both cases, L I is defined as earlier.  which means that some fraction of its diverted and gaged water can re-enter the study reach of the river from seepage (Figure 1). Similarly, because the Fort Lyon Canal diversion is about 38 km upstream from the Las Animas streamgage some fraction of its diverted and gaged water can re-enter the river from seepage and return flow (Figure 1). Both situations result in water volumes and loads being potentially counted a second time at the Las Animas streamgage. Such double counting of water would overestimate outputs from the study reach. The problem is acknowledged, but there is no easy means to address it. However, two of the effects observed in the analysis (minor salt storage and Se removal) would be underestimated by double counting of loads, yet they remain observable in the results. Observation of the other effect (salt flushing) relies on year-to-year differences to discern; and therefore a persistent issue like double counted water seems unlikely to enhance or obscure it. The load of TDS generated within the study reach (L TDS,WR ) can be determined as simple subtraction of inputs from outputs by Here, the estimates of output load using either SC from Rocky Ford (L TDS,O* ) or Las Animas (L TDS,O ) to represent SC in the Fort Lyon Canal yield underestimates and overestimates, respectively, of the load generated within the reach. The same calculations made for TDS are made for dissolved Se and U.
It should be noted that the within-reach loading terms L TDS,WR , L Se,WR , and L U,WR include the effects on dissolved constituents contributed by the components DS SW , DS GW , and Q GW-net that were omitted as (Equation 1) was simplified to (Equation 2). As was indicated for the water-volume analysis, this is not suggestion that the other terms are unimportant, but rather that their effects are all implicitly included in the within-reach loading term. Thus, the hypothetical 123,000 ML change in groundwater storage presented previously could be multiplied by a given groundwater TDS concentration to estimate how a value of DS TDS,GW might compare to L TDS,I or other components of the mass balance. However, L TDS,WR (Equation 5) gives an actual measure of the combined net annual magnitude of DS TDS,SW , DS TDS,GW , and Q TDS,GW-net along with other mobilization/immobilization mechanisms described in the discussion. To the extent that a component like change in groundwater TDS storage (DS TDS,GW ) has a large net effect on the mass balance in a given year, that large net effect is apparent through its contribution to L TDS,WR .
Effects of uncertainty were assessed for the mass-balance analysis. Uncertainty for discharge was described using AE10% bounds for the USGS streamgages for Avondale, Huerfano, and Las Animas (Novak 1985) as well as the DWR streamgages at the Fort Lyon and Fort Lyon Storage Canals. Uncertainty was also described using AE10% bounds for the SC records. Uncertainty in the linear, regressed SC-constituent relations (Figures S1 and S2) was represented by the standard error of both the slope and intercept. Uncertainties from streamgage records and regressions were assumed to be independent and propagated in quadrature for individual years in the mass-balance analysis using standard error-propagation methodology (e.g., Taylor 1997), but treated as additive across years for summing totals for 31 years. Error bars in subsequent figures represent the propagated uncertainties.

Annual Precipitation Data
Annual precipitation totals for the years 1987-2018 at two municipalities in the study reach (Rocky Ford and Las Animas) are used for perspective in interpreting mass-balance analysis results. Precipitation data were obtained from the National Oceanic and Atmospheric Administration, National Centers for Environmental Information's Climate Data Online (https://www.ncdc.noaa.gov/cdo-web/).

Synoptic Samplings
The synoptic data from the study reach illustrate seasonal differences in discharge between spring runoff and base flow (Figure 3a), though such natural variations are highly regulated by upstream structures (Figure 1). Decreases in discharge downstream reflect diversions from the river and consumption. Seasonal differences in SC also reflect differences between spring runoff, when waters are more dilute, and baseflow when solutes are more concentrated, but SC generally increases downstream (Figure 3b). Concentrations of Se in the river change seasonally, but trends through the reach are muted in several synoptics (Figure 3c). In contrast, concentrations of U (Figure 3d) largely match the patterns in SC, generally increasing with distance downstream (Figure 3b).

Evaporation Assessed by Hydrogen and Oxygen Isotopes
The fraction of water evaporated, as traced by the dD and d 18 O of water sampled monthly in the JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION Arkansas River at Las Animas is compared to discharge at the time of sampling in Figure 4. Not surprisingly, the two have an inverse relation. Effects of evaporation are most pronounced during baseflow and least pronounced during spring runoff in 2017. However, decline in the fraction evaporated from summer into the cooler weather of winter in 2017 suggests the isotopes do not simply trace runoff vs. baseflow mixing. The decline suggests relatively current evaporative conditions are also reflected.
The discharge-weighted average fraction of water evaporated for the sampling period (Figure 4) is 4.1%. The subsequent mass-balance analysis assesses the magnitude of annual consumption of water in the study reach over 31 years, but the 2017 value is 17%. Using that 2017 value, evaporation accounts for 24% of the evapotranspiration losses represented by consumption in 2017 and transpiration therefore accounts for 76% of consumption. Considerable uncertainty underlies such estimates, but the indication is that transpiration is the majority share of consumption. Prior work with water isotopes indicated that evaporation was notable in the reach (Ivahnenko et al. 2012), but these values provide numerical estimates of its effect.

TDS Mass-Balance Analysis
The regressions between SC and TDS are all significant and R 2 values of~0.99 for all indicate SC as a good proxy for TDS (Table 1). The standard errors are relatively small and contribute little uncertainty to the subsequent mass-balance analysis. Annual values for each streamgage in the mass-balance analysis, results of calculations, and associated uncertainties are in the Supporting Information. Discharge and TDS load for the Huerfano River are large enough to be notable contributions in some years but not others.
Total discharge outputs from the study reach are generally less than inputs. Because diverted water that is not evapotranspired rejoins the river as return flow, the difference represents primarily consumption of discharge by evapotranspiration (Figure 5a). Dividing the reference point Q CON by Q I gives the fraction of input discharge consumed. For 1988-2018 that  Table S1. value ranges from a consumption loss of 46% to a gain of 2% with a median loss of 33%. The gain of 2% occurred in 1999, a wet year that included storm-induced flooding along the study reach and followed other years of high water availability. The second lowest value was a loss of 6%. Annual precipitation totals ( Figure 6) provide important perspective on consumption. Consumption of discharge drops notably in some years with above average rainfall like 1995 and 1999. However, consumption also drops in drought years of low discharge input to the reach, like 2002 and 2012, though it increases as a fraction of inputs in such years. Discharge-weighted, annual TDS concentrations in inputs are dominated by the main stem near Avondale and vary relatively little in most of the record (Figure 5b), ranging 340-720 mg/L, with a median of 480 mg/L. The effect of consumption on TDS concentrations is a useful reference point, and the TDS concentrations at Las Animas that would be expected if consumption were the only influence (C Con TDS;LA ), range 420-1,200 mg/L with a median of 670 mg/L. Such values are most elevated in years of low runoff. Actual TDS concentrations at Las Animas are greater and more variable, ranging 680-2,150 mg/L with a median of 1,200 mg/L (Figure 5b). The difference between the two values is reflective of within-reach sources of TDS.
Because discharge-weighted, annual TDS input concentrations (Figure 5b), which predominantly reflect the Arkansas River at Avondale, are relatively constant they reduce variability of year to year TDS input loads (Figure 5c), which become primarily a function of input discharge (Figure 5a). In contrast, TDS output loads vary substantially across years ( Figure 5c). Although the total output load calculated using the Rocky Ford SC record to represent SC in the Fort Lyon Canal (L TDS,O* ) only covers 2013-2018, it is a useful check on the output load calculated using SC at Las Animas to represent SC for the Fort Lyon Canal (L TDS,O ). Output load L TDS,O* is acknowledged as an underestimate and L TDS,O as an overestimate, and L TDS,O exceeds L TDS,O* by 22-56% (average 22%) (Figure 5c). This is a somewhat narrower difference than might be expected based upon SC at the Rock Ford and Las Animas streamgages being 12% less than and 35% greater than that in the Fort Lyon Canal, respectively. The estimate of within-reach TDS load (Figure 5d) calculated (Equation 5) with L TDS,O is 18-36% greater (average 26%) than that calculated with L TDS,O* is. These differences illustrate the magnitude of the gap between the known overestimate and known underestimate of output loads. The next section examines how within-reach loading changes from year to year. The perspective gained is therefore not dependent on the specific magnitudes of loads, which obviously carry considerable uncertainty, but rather changes through time, for which relative differences are important.

Salt Flushing and Salt Storage
Examination of the patterns in Figure 5 reveals that, in a given year or with a possible lag of about a year, times of greater water availability flush salts from the surrounding landscape into the Arkansas River. The largest flushing event occurred in 1996-2000. It was both pronounced and sustained enough to be apparent in a locally weighted scatterplot smoothing curve for SC in the Arkansas River at Las Animas (Miller et al. 2010). Here, the event is apparent in both the output loads and within-reach loads, along with the concentrations of TDS at Las Animas (Figure 5b, 5c, 5d). The output load during 1996-2000 averaged 1,470,000 AE 154,000 Mg/yr, which is 131% greater than the 640,000 AE 63,000 Mg/yr for the 26 other years in the 1988-2018 time period. The within-reach load during 1996-2000 averaged 976,000 AE 169,000 Mg/yr, which is 222% greater than the 303,000 AE 80,000 Mg/yr for the other 26 years. Concentrations of TDS in the Arkansas River at Las Animas during 1996-2000 averaged 1,550 AE 180 mg/L which is 25% greater than the 1,240 AE 150 mg/L for the other 26 years. The difference with average TDS concentrations is less striking because the long record includes years of extreme drought like 2002 and 2012 when TDS concentrations were high because of extreme low flows. High average concentrations are expected in drought years, but not during times of water abundance.
The primary driver of the big 1996-2000 flush appears to have been the sustained, high water availability from large Rocky Mountain snowpacks in 1995-1999, apparent as above-average discharge inputs (Figure 5a). Above average water availability allowed greater diversions from the river and likely allowed farmers to irrigate fields left fallow in drier years and apply more water to all fields. The amount of deep percolation occurring within the reach would have increased in terms of both total area influenced and quantity per area. Temporary and localized increases in groundwater storage are apparent during this time ( Figure S3). The period 1995-1999 also included several years of above average rainfall, particularly in the downstream end of the study reach ( Figure 6). Such precipitation might have mobilized some salt out of uncultivated land, but would also have contributed to greater water availability in irrigated fields.
Several flushes of smaller magnitude are also apparent in Figure 5. Visible as peaks in within- In contrast to the years of salt flushing, years of likely salt storage are also apparent. It warrants consideration that within-reach loads are the net effect of multiple sources and processes. In every year, some salt enters the river from nonirrigated, noncultivated parts of the landscape via tributaries and groundwater. For within-reach load to be so low in years like 2002-2003) particularly with L TDS,O understood to be an overestimate output loads, some storage of salts diverted out of the river in irrigation waters must be occurring.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
The influence of upstream water supply and rainfall on storage and flushing of salts can be discerned from Figures 5 and 6, but scatterplots in Figure 7 further clarify the links. Regression analysis finds significant relations with within-reach TDS loading for both current-year discharge and previous-year discharge at the Avondale streamgage (Figure 7a, 7b). Multiple regression that includes those two variables and previous-year rainfall and current-year rainfall finds previous-year discharge at the Avondale streamgage (p = 0.026) and current-year discharge (p = 0.007) to be significant, but neither of the rainfall variables to be significant (p = 0.40 and 0.49, respectively). Thus, it appears that the availability of water for diversion and irrigation is a stronger influence on the flushing of salts than precipitation within the study reach.

Sources of Flushed and Stored Salts
The new understanding of a flushing effect provides important context for understanding and managing salinity in the Arkansas River Valley and perhaps other irrigated landscapes. Sources and FIGURE 6. Annual precipitation totals for two municipalities in the study reach. Mean-annual precipitation for each municipality for 1987-2018 is shown as a horizontal line in corresponding colors.
FIGURE 7. Regression analyses of within-reach TDS load as a function of (a) the current year's discharge at the Arkansas River near Avondale and (b) the previous year's discharge for Avondale.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
processes associated with variables (DS SW , DS GW , and Q GW-net ) contained in the within-reach load term seem unlikely to explain the patterns. Movement of water from surface-and groundwater storage to the river does mobilize salts, but storage in off-channel reservoirs (DS SW ) and groundwater storage (DS GW ) should increase during times of increased water availability (compare Figure 5d and Figure S3). Therefore, gains in storage would seem to be counteracting mobilization at about the same time that salts are flushed out, but another mobilization effect is larger. Increased recharge in the nonirrigated portion of the landscape might drive increased Q GW-net , but the lack of significant correlation with rainfall in the multiple regression suggests salts mobilized from the nonirrigated landscape do not play large roles in the observed flushes.
The flushing concept envisioned here is of salts in different parts of the landscape lacking only a sufficient flux of water to dissolve and mobilize them. Three categories of such salts may be considered (Figure 8). The first category is river-derived salts from irrigation water. Salts from irrigation water become concentrated by evapotranspiration in and underneath cultivated fields. They accumulate when deep percolation is insufficient to flush them completely or quickly to the water table with subsurface return flow. Such salts are stored temporarily in the landscape, precipitated as minerals like gypsum and other sulfates, and excessive accumulation is referred to as secondary salinization. During years of greater runoff and irrigation, increased deep percolation mobilizes them to groundwater and then to the river. Efflorescent salts from shallow groundwater sourced from diversions would also be temporarily stored, riverderived salts.
A second category is weathering-derived salts. Weathering-derived salts and particularly gypsum are common in soils where average precipitation is insufficient to flush them away (Herrero and Porta 2000;. Such salts can contribute substantially to loading; in an irrigation district in Spain, 82% of the exported TDS resulted from gypsum dissolution (Barros et al. 2012). In western Colorado, the Mancos Shale, which is coeval to formations in the Arkansas River Valley (Bern and Stogner 2017), has weathered over~20,000 years to generate large deposits of salts in the unsaturated zone (Tuttle et al. 2014). Under nonirrigated conditions, there is insufficient deep percolation to mobilize such salts. With irrigation, however, deep percolation mobilizes them to surface water via return flow (Mast et al. 2014). Gypsum is present in soils of the Arkansas River Valley (Cooper et al. 2006) and descriptions of soil series like the Cheraw, Hackamore, Kornman, Las Animas, Limon, Minnequa, Rocky Ford, and Timpas include gypsum in subsurface horizons, in some cases up to 5 wt.% (Soil Survey Staff 2019).
A third category is natural and evaporative salts. Such salts predominate in groundwater outside the canal command areas and would account for base flow TDS under nonirrigated conditions. Years when more water than average is diverted from the river into the surrounding landscape could perturb groundwater movement outside the immediate footprint of irrigation. A more visible source of evaporative salts is Cheraw Lake, the highly saline lake (e.g., 17,000 mg/L TDS upper layer, 60,000 mg/L lake bottom) located north of La Junta, Colorado (Figure 1) in a topographic low point with restricted surface drainage (MacDonnell 1989). Satellite imagery of Cheraw Lake shows changes in its areal extent in different years ( Figure S4). It appears exceptionally full in 1998 during the 1996-2000 flushing event.
Determining the dominant sources of salts that are flushed ( Figure 8) could shed light on interactions with best management practices (BMPs). If temporarily stored river salts predominate in flushing years, then conversion from furrow to sprinkler irrigation, deficit irrigation and reduced leaching fractions will help store salts in drier years, but leave them vulnerable to mobilization when water applications are increased. If weathering-derived salts predominate in intermittently irrigated fields, including those in lease fallow programs, they could be large contributors of mobilized salt in years of water abundance. Lining and sealing to reduce seepage from the major canals (Dash 1994), smaller laterals, and head stabilization ponds might reduce mobilization of weathering-derived salts. In Mancos Shale landscapes of western Colorado, reducing such seepage has reduced salt and Se mobilization from such sources to rivers (Butler 2001). In the Arkansas River Valley, it would be worth considering if weathering-derived salts have been depleted along certain major flow paths (Figure 8), like those from major canals operated for >120 years (Sutherland 2008). If natural and evaporative salts predominate, the locations of the salt sources will be of interest. Discrete sources, like Cheraw Lake, might be amenable to engineering solutions (Smith 2018). Widely distributed sources of natural salts, including groundwater outside canal command areas, would be difficult to address.

Selenium and Uranium Mass-Balance Analysis
The mass-balance framework applied to TDS was also applied to Se and U. Correlations with SC used to convert SC to Se or U concentrations (Figures S1 and S2; Table S2) are not as tight as those for TDS (Table 1). Correlations between SC and U are quite strong, with R 2 values ranging 0.81-0.96. However, correlations between SC and Se are weaker on the main stem of the Arkansas River (R 2 = 0.50-0.84) and nonexistent for the Huerfano River (R 2 = 0.0004). Weaker correlations for Se likely reflect partial biogeochemical immobilization or other processes that remove Se relative to other constituents. The absence of correlation between SC and Se at the Huerfano River site (Table 2) resulted in large (>2,000%) standard errors in the SC-Se regression. Propagation of that uncertainty through the mass-balance calculations resulted in large final uncertainties. Therefore, uncertainty for Huerfano River Se concentrations was set to 100% for error propagation. ). An immediately apparent pattern is that in most years Se concentrations at Las Animas are less than those expected based upon consumption (Figure 9a). This is the first piece of evidence that some fraction of Se in the river near Avondale is being retained in the study reach. In contrast, the overall patterns in U concentrations (Figure 9b) match those of TDS concentrations (Figure 5b) quite well and indicate that additional U is mobilized in the reach.  (1) river-derived salts from irrigation water temporarily stored in and below fields; (2) weathering-derived salts accumulated over geologic time; and (3) natural and evaporative salts from areas outside of canal-command areas and in low spots with limited natural drainage, like Cheraw Lake. Lengths and widths of arrows indicate relative magnitudes of water movement. Flow paths that get consistent large fluxes of water every year are represented by an "A." The analysis of loads provides the second piece of evidence that the study reach is a net sink for Se in some years (Figure 10a). Within-reach loading of Se is negative in 12 years of the 31-year analysis, in which within-reach loading is overestimated because Fort Lyon Canal loads are based upon calculations using the Las Animas streamgage SC. Positive within-reach loading occurs in years of high input discharge ( Figure 5a) and salt flushing. In the shorter record, which underestimates within-reach loading, that loading is negative in 5 of 6 years and values are more negative (Figure 10a). Across the long 1988-2018 record, Se input at the upstream end of the reach averages 6.6 AE 1.5 Mg/yr and output averages 7.6 AE 1.3 Mg/yr. On average, a comparatively small amount, 1.0 AE 2.0 Mg/yr of Se, is mobilized within the study reach. The total amount of Se that has entered from upstream over 31 years is 205 AE 46 Mg and outputs total 236 AE 39 Mg indicating a comparatively small amount of 31 AE 61 Mg has been mobilized in the reach. Over 31 years, only 13% of the Se leaving the study reach is accounted for by mobilization from within the reach. As the output loads are likely overestimates, net mobilization might be even less, as suggested by the load analysis using the Rocky Ford SC to represent the Fort Lyon Canal (Figure 10a).
In contrast to Se, the load analysis consistently indicates substantial within-reach loading of U (Figure 10b  accounted for by mobilization within the reach. Again, the shorter analysis based upon using the Rocky Ford SC record to represent the Fort Lyon Canals is similar, but shows smaller amounts of within-reach loading (Figure 10b). The greater than average within-reach loading of U in 1996-2000 and other correspondence between U and TDS loading (Figures 5d and 10b) suggests that U responds to the same flushing patterns as TDS.
Comparing the Se and U loading analyses, it is striking that Se mobilization is comparatively small in some years and net Se removal is indicated in other years. It is therefore worth considering a third line of evidence regarding Se, one found in the synoptic sampling data. Flows decrease with distance downstream in response to both diversions and consumption (Figure 3a). Concentrations of TDS (SC) increase by maximum factors of 1.6-3.3, and U concentrations increase by maximum factors of 2.0-4.3 with distance downstream (Figure 3b, 3d). In contrast, Se concentrations increase by maximum factors of 1.1-2.3 with distance downstream (Figure 3c). For Se concentrations to increase so little along the river as consumption decreases flows, Se mobilization must be minor and/or some Se must be removed. More perspective on the immobilization of Se and where it is occurring can be obtained from Se/TDS and Se/U ratios in the synoptic data ( Figure 11). Both show steady decreases with distance downstream across seasons and flow regimes. Such patterns bolster the conclusion that little Se is mobilized, or some Se is removed even as substantial additional TDS and U are mobilized. Additionally, the patterns of consistent downstream decrease in these ratios suggest that Se removal is occurring along the length of the study reach rather than one distinct area.

Selenium and Uranium Mobility
The discrepancy between mobilization of U vs. possible net removal of Se would not have been predicted based upon prior work showing strong positive correlations in their concentrations in the Arkansas River and its tributaries (Zielinski et al. 1995;Gates et al. 2009). Indeed, Se and U are positively correlated in the synoptic data collected here (R 2 = 0.44, p < 0.001). Existence of that correlation can be attributed to strong effects of both evapotranspiration and dilution driving concentrations of all solutes up or FIGURE 10. Calculated input loads and within-reach loads for selenium (a) and uranium (b). Negative loads for Se in most years suggest immobilization within the reach.
down in unison and only partial (at the river scale) removal of Se. Spurious correlations in untransformed concentration data are common (Reimann et al. 2008). Selenium could be removed from waters within the study reach in several ways. Highly mobile selenate (SeO 2À 4 ) will be the dominant species in neutral to alkaline and oxic soils and waters present throughout much of the study reach (Gates et al. 2009). Under low oxygen and nitrate concentrations, selenate can be reduced to selenite (SeO À 3 ), and subsequent adsorption to mineral surfaces and/or interactions with organic matter remove it from solution (Balistrieri and Chao 1987; Fern andes-Mart ınez and Charlet 2009). Under more strongly reducing conditions, or in the presence of a reducing phase like FeS 2 , selenite can be reduced to elemental Se or selenide (Se 2À ) and removed from solution by precipitation. Anoxic groundwater, wetlands, and riparian areas are all good settings for Se immobilization (Masscheleyn and Patrick 1993;Naftz et al. 2005). Column experiments using soil from within the study reach showed substantial Se removal from infiltrating waters, and batch experiments using shales showed complete reduction and removal of Se from solution under conditions in which nitrate was also reduced (Bailey et al. 2012). Thus, Se could be immobilized either in soil or in groundwater in the Arkansas River Valley. A survey of riverbed and bank sediment in the valley found reduced and organic forms of Se in both (Romero 2016). Volatilization is another process that could remove Se and one that has been suggested previously for the Arkansas River (Whittemore 2009). Selenium volatilization has been linked to biological processes in wetland and aquatic environments and can generate fluxes that are significant to the mass balance of such systems (Zhang and Moore 1997;Diaz et al. 2009). Another potential removal process for Se is uptake by crops that are then exported from the region (Mikkelsen et al. 1988). Data on Se concentrations in crops and the fraction of harvest exported vs. used locally are not available to estimate such Se export, but considering that corn and alfalfa production in the study area each exceed 100,000 Mg annually, such export could be a factor (USDA 2016). Although some Se is undoubtedly mobilized within the study reach, the mass-balance analysis indicates much of that mobilization is offset by removal. Both mobilization and removal processes likely operate simultaneously in a patchwork of sources and sinks.
In contrast to the retention mechanisms described for Se, U is likely to remain mobile in solution. Uranium (VI)-carbonate complexes will be the dominant dissolved species in the alkaline waters of the region and their adsorption can be negligible (Hsi and Langmuir 1985;Duff and Amrhein 1996;Ivahnenko et al. 2012). Relatively conservative behavior by U already in solution, in combination with mobilization of additional U derived from Cretaceous marine rocks (Zielinski et al. 1995;Bern and Stogner 2017) would explain the net gains of U seen here (Figure 10b).
Apparent removal of Se in the study reach is a previously unrecognized water-quality benefit. The removal processes described above can occur in soil, groundwater, wetland, and cropping systems, and much of the river's water has passed through a combination of these at least once before coming back as return flow. It is therefore likely that Se removal occurs as a byproduct of the substantial diversion of water from the river.

Comparisons to Previous Studies
Results and interpretations presented here run counter to some derived from recent modeling efforts that examined a similar reach of the Arkansas River (Bailey et al. 2014;Bailey et al. 2015;Tavakoli-Kivi and Bailey 2017). Periodic salt flushing, a major driver of water quality observed here, has not emerged from modeling efforts, which have also FIGURE 11. Ratios of (a) Se/TDS and (b) Se/U for the synoptic samples depicted in Figure 3.
underpredicted groundwater concentrations of sulfate (Tavakoli-Kivi and Bailey 2017). However, the dissolution of gypsum present in regional soils (Soil Survey Staff 2019) is acknowledged as something that warrants future consideration (Tavakoli-Kivi and Bailey 2017) and has been modeled in a way that approximates averages, but not temporal patterns (Tavakoli-Kivi et al. 2019). Net Se removal in some years, as indicated by the analysis here, has not emerged from modeling efforts (Bailey et al. 2014;Shultz et al. 2018) although immobilization is a desired outcome of management practices (Bailey et al. 2015). The different results may arise from the different approaches. This paper uses an empirical approach focused on broad patterns over long time steps using a simplified-mass balance concept. In contrast, the models have simulated complex and interacting processes like oxidation of shale, mineralization of litter and manure, and root uptake of S, N, and Se along with spatially distributed factors like crop type, fertilizer application, canal seepage, and groundwater in a three-dimensional model on a daily time step (Bailey et al. 2014;Tavakoli-Kivi et al. 2019). Such models use intensive data coverage, but generally obtained over short time frames of a few years. Modeling such complexity when there is considerable uncertainty regarding the target parameters of TDS, Se, and U loads being mobilized is a substantial challenge . Regarding salinity in particular, simpler models may perform better (Burkhalter and Gates 2005). It is desirable for any model being used to predict responses to management actions to represent the broad patterns in the system modeled (Bredehoeft 2005). This study suggests that some do not.
BMPs recommended for the region (Bailey et al. 2015) are still sensible, although there are several possibilities for types of flushed salts and U (Figure 8). Reducing deep percolation by conversion from furrow to sprinkler irrigation is a strategy that has worked in another setting with gypsum-bearing soils (Herrero and Castañeda 2018; Jim enez-Aguirre and Isidoro 2018). Piping laterals (Butler 2001), land fallowing and reduced irrigation (Bailey et al. 2015), and lining head stabilization ponds also reduce deep percolation. Lowering nitrate concentrations in groundwater and enhancing riparian buffer zones could enhance Se immobilization (Bailey et al. 2015). When conducting trend analysis (e.g., Miller et al. 2010) or assessing the effectiveness of such practices, the salt storage and flushing patterns observed here ( Figure 5) warrant consideration.
The type of annual-scale salt flushing phenomenon detailed here has apparently not been documented for other rivers, though episodic salinization occurs in urbanized rivers (Haq et al. 2018) and physical disturbance in semiarid watersheds can also mobilize salts (Bern et al. 2015). Climate variability is recognized to generate pulses of sediment transport (Kaushal et al. 2010). A general expectation is that increased water abundance dilutes salinity, a pattern observed even at the decadal scale (Tillman et al. 2019). However, most long-term studies seek trends in water quality (e.g., Hirsch et al. 2010) and the short pulses of salt flushing might not show up in trend analysis. Whether similar patterns exist in other intensively irrigated landscapes warrants investigation where robust long-term records of flows and salinity exist.

CONCLUSIONS
The analysis presented uses relatively straightforward calculations and a simplified mass-balance concept to assess an irrigated, agricultural reach of the Arkansas River on an annual time step for 31 years. Robust records of discharge and SC, and discrete sampling for TDS, Se, and U made it possible to calculate inputs and outputs for the study reach with acceptable accuracy. Other mass balance components physically located between inputs and outputs, like groundwater storage change, could not be calculated with acceptable accuracy for similar time steps and scales. The approach taken here allowed net effects of within-reach annual consumption, storage change, mobilization, and immobilization to be inferred through differences between inputs and outputs. Although absolute magnitudes of calculated values carry considerable uncertainty, relative variations in loads through time and consistency of patterns allow several conclusions to be drawn. First, consumption of water by evapotranspiration is considerable, with a median loss of 33%, and that loss drives corresponding increases in dissolved constituents. Transpiration appears to account for most of the consumption. Second, salts are flushed from the landscape following years of abundant water, with a lag time of about one year or less. Third, low to negative within-reach loading during years of low water abundance or drought suggests some salts from the river are temporarily stored in parts of the surrounding landscape. Finally, comparatively little Se is mobilized from the study reach and net Se removal appears to occur in some years. The effect is likely related to diversion of river water into the surrounding landscape for irrigation. Other evidence for Se removal comes from synoptic data that show only minor downstream Se concentration increases as consumption removes water, and steady downstream declines in Se/TDS and Se/U ratios that suggest Se is partially removed throughout the reach, rather than in one particular location, even as substantial TDS and U are mobilized. The insights presented here are only possible because of decades of continuous monitoring of discharge and SC, combined with discrete water-quality sampling. Similar patterns of salt flushing and redox-sensitive element immobilization may affect other rivers in other irrigated, semiarid landscapes.

SUPPORTING INFORMATION
Additional supporting information may be found online under the Supporting Information tab for this article: Four supplemental figures and six supplemental tables of calendar year inputs and outputs (total discharge, average concentrations, total loads) from the mass-balance calculations for salinity, Se, and U along with associated uncertainties.

ACKNOWLEDGMENTS
This work was directly supported by the Upper Arkansas, Lower Arkansas Valley, and Southeastern Water Conservancy Districts, along with Colorado Springs Utilities, Aurora Water, and the Pueblo Board of Water Works. The work was indirectly supported by the seven entities that have funded streamgage operations and continuous and discrete monitoring along the river for the last three decades. We thank employees in the USGS Colorado Water Science Center (CWSC) Pueblo office for the decades of highquality data collection. Numerous CWSC employees directly assisted with this project including William Banks, Dustin Ethredge, Abby Keith, Roderick Ortiz, and Pam Tello. We thank David Naftz, Mike Weber, and three anonymous reviewers for conversations and comments that improved the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.