Purpose and Scope

The purpose of this report is to describe the methods and results of a study to evaluate salinity and nutrients in the Heart River Basin in western North Dakota. Water-quality and streamflow data used in the study were compiled from 1970 to 2020 using the National Water Quality Monitoring Council Water Quality Portal (WQP; National Water Quality Monitoring Council, 2020) and National Water Information System (NWIS) database (U.S. Geological Survey, 2020). Streamflow characteristics were analyzed from 1970 to 2020 to examine hydroclimatic patterns that may explain changes in water quality in the Heart River Basin over time. Statistics were computed for selected constituents and sites (fig. 1) to describe spatial patterns in concentration across the basin. Water-quality trend analysis was completed for chosen sites with sufficient data to evaluate how selected water-quality constituents have changed over time. Water-quality trends were analyzed for a recent trend period (1999–2019) at two tributary sites (Green River near New Hradec, North Dakota [site 1] and Antelope Creek near Carson, N. Dak. [site 10]; fig. 1) and two main-stem sites (Heart River near Richardton, N. Dak. [site 5] and Heart River near Mandan, N. Dak. [site 22]) and for a historical trend period (1974–2019) for one tributary site (site 1) and two main-stem sites (sites 5 and 22). Geochemical modeling was used to evaluate chemical reactions relating to salinity in two reaches of the Heart River for two different periods (1974–99 and 1999–2019). Constituent loads and yields were estimated for selected sites in the Heart River Basin from 2013 to 2020 for total dissolved solids (TDS), sulfate, sodium, chloride, and total phosphorus. Monthly, annual, and total loads were computed for each site and constituent, and mean annual yields were computed to normalize loads to the drainage area upstream from the selected sites. Total loads for the period of 2013–20 were used to estimate a simplified mass balance in the lower Heart River Basin to identify areas in the basin that have the largest contributions of constituent mass in the Heart River.

Description of Study Area

The Heart River Basin, located in western North Dakota, is approximately 120 miles long and about 30 miles wide with a contributing drainage area of about 3,350 square miles (Maderak, 1966; North Dakota Department of Environmental Quality, 2021a; fig. 1). The elevation ranges from about 1,650 feet above North American Vertical Datum of 1988 (NAVD 88)at its mouth near Mandan, N. Dak., to about 2,900 feet above NAVD 88 in Billings County, N. Dak (not shown). The Heart River is a dendritic system (Ritter and others, 2011) located between the Knife River to the north and the Cannonball River to the south (not shown). The Heart River is impounded by two dams, the Dickinson Dam (E.A. Patterson Lake) in the upper basin near Dickinson, N. Dak., and the Heart Butte Dam (Lake Tschida) in the central portion of the basin (fig. 1). In this report, the “upper” Heart River Basin refers to the area above Heart Butte Dam (Lake Tschida), and the “lower” basin is the area below Heart Butte Dam (Lake Tschida). Dickinson Dam was completed in 1950 (Linenberger, 1996), and Heart Butte Dam was completed in 1949 (Simonds, 1996). Both reservoirs are operated by the Bureau of Reclamation as water supplies for irrigation, recreation, and flood control as part of the Pick-Sloan Missouri Basin Program. Additionally, E.A. Patterson Lake serves as a water supply for the City of Dickinson, N. Dak.

Western North Dakota geology generally consists of shales, sandstones, and alluvial materials. Specific formations underlying the Heart River Basin are the Ludlow, Cannonball, Bullion Creek, Sentinel Butte, Golden Valley, Chadron, and Brule Formations (figs. 2A and 3; Murphy and others, 2009).

The dominant formations in the basin were the Ludlow, Cannonball, and Bullion Creek Formations. The Paleocene Ludlow Formation is considered the basal unit of the Fort Union Group (fig. 3; Carlson, 1983). The Ludlow Formation is characterized by brown or gray with yellowish-brown or gray claystone, siltstone, and sandstone beds alternating with lignite or lignitic shale beds (Carlson, 1982, 19838; Murphy and others, 2009). The Paleocene Cannonball Formation is a marine deposit that overlayed the Ludlow Formation in western North Dakota (fig. 3; Murphy and others, 2009). This formation is characterized by alternating beds of sandstone, siltstone, and mudstones, with mudstones as the dominate lithology (Carlson, 1982, 19838; Murphy and others, 2009). The Paleocene Bullion Creek Formation (formerly Tongue River Formation pre-1977) is a nonmarine deposit that overlies the Cannonball Formation in western North Dakota (fig. 3; Murphy and others, 2009). The Bullion Creek Formation consists of interbedded layers of clay, silt, sand, and lignite (Clayton and others, 1977). The descriptions of these units were used to determine mineral types for geochemical modeling.

Land use in the basin is dominated by agriculture, with 45 percent of the land used for grasslands, 39 percent used for cultivated crops, and 9 percent used for hay or pastureland (Homer and others, 2015; fig. 2B). Grasslands are defined as areas with 80 percent or more of the total vegetation being graminoid or herbaceous (Homer and others, 2015). Additionally, grassland is not subjected to intensive management practices but may be used for grazing (Homer and others, 2015). Cultivated crops are areas used to produce annual crops and are managed intensively (Homer and others, 2015). Hay/pastureland are areas of grasses, legumes, or grass-legume mixes planted for grazing or production of seed (Homer and others, 2015). In the upper basin, there is a higher percentage of cultivated crops compared to the lower basin (43 and 35 percent, respectively), whereas in the lower basin there is a higher percentage of grassland compared to the upper basin (51 and 39 percent, respectively) (Homer and others, 2015). To support the agriculture in the subhumid climate of the basin, irrigation began in the late 1940s (Bureau of Reclamation, 2021). Conversion from flood to pivot irrigation practices has been occurring within the basin since the late 1990s (Chad Skretteberg, Western Heart River Irrigation District, written commun., 2021; James Weigel, Bureau of Reclamation, written commun., 2021). Only about 1 percent of the land use is urban (Homer and others, 2015), with two major cities in the basin, Dickinson and Mandan, N. Dak., having populations of 17,787 and 18,331, respectively (U.S. Census Bureau, 2010).

Soil characteristics in the basin have slight variations but are dominated by hydrologic soil groups C and D that cover about 71 percent of the basin (Soil Survey Staff, 2020; fig. 2C). Hydrologic soil groups are based on the saturated hydraulic conductivity of the least transmissive layer, depth to impermeable layer, and depth to a water table (U.S. Department of Agriculture, 1972). The hydrologic soil group generally is defined by how well the soil transmits water when fully saturated in each unit (U.S. Department of Agriculture, 1972). The soils are grouped by A, B, C, and D soils where the A and B soils transmit water through the soil column well and the C and D soils do not transmit water well and generally are not well drained naturally (U.S. Department of Agriculture, 1972). Within the basin, the hydrologic group A and B soils generally are present along the streams and rivers in the basin, which was expected because these consist of alluvium and terrace deposits (fig. 2C). The group C and D soils were generally formed on top of the local bedrock, which consists of mostly shales (fig. 2A, B).

Streamflow Analysis

Streamflow changes were analyzed in the Heart River Basin using the flow history component of Exploration and Graphics for RivEr Trends (EGRET) program (Hirsch and De Cicco, 2015). The flow history component provides tabular and graphic output for several streamflow statistics that can be used to assess how streamflow may be changing over time. Sites were selected for use in the EGRET analysis that had long term (50 years minimum) flow records with little to no missing data for the period of 1970–2020 and were also used for water-quality trend analysis. Sites 1, 5, and 22 (table 1 and fig. 1) were selected for the EGRET analysis. The flow history in EGRET provides a statistical description of long-term variability in streamflow using the daily streamflow record for a given streamgage.

Four streamflow statistics were selected for streamflow analysis and plotted to analyze changes in streamflow from 1970 to 2020. Using the flow history component of EGRET, plots of the smoothed trend in annual maximum daily, annual mean daily, and annual 7-day minimum were used to identify long-term changes in streamflow. The smoothing methods used by EGRET are described in detail in Hirsch and De Cicco (2015). The results of the EGRET analyses aided in the interpretation of the historical and recent water-quality trends.

A complete record of daily streamflow is required for the analysis of water-quality trends and the estimation of constituent loads. Sites on Antelope, Big Muddy, and Sweetbriar Creeks (site 10, Big Muddy Creek near Almont, N. Dak. [site 18], and Sweetbriar Creek near Judson, N. Dak. [site 21], respectively; fig. 1 and table 1) were operated as seasonal streamgages, where streamflow data were generally only recorded from April through September of each year, requiring an estimation of daily streamflow for the remaining months to use in the various analyses. Streamflow was estimated between the months of October and March using nearby streamgages with recorded data for those months and similar contributing drainage areas. The USGS streamgage on the Cannonball River at Regent, N. Dak (USGS site number 06350000) was used to estimate streamflow at site 10; Spring Creek at Zap, N. Dak. (USGS site number 0634000) was used to estimate streamflow at site 18; and Square Butte Creek below Center, N. Dak. (USGS site number 06342260) was used to estimate streamflow at site 21 (fig. 1). Streamflow from 1950 to 2020 from paired sites was compared using a regression analysis of daily streamflow. The Leaps package in R (Lumley, 2020) was used to evaluate the possible subsets of models, and the coefficient of determination (R²) was used to evaluate the possible models.

Constituent and Site Selection

All available water-quality data in the Heart River Basin were initially compiled for the period of 1970–2020 from the WQP (National Water Quality Monitoring Council, 2020) to evaluate their use in the various analyses described in this report. The WQP contains water-quality data collected by multiple agencies for the area of interest including data from the U.S. Environmental Protection Agency, the USGS NWIS, and North Dakota Department of Environmental Quality databases. USGS streamflow data associated with water-quality sites were compiled from the USGS NWIS (U.S. Geological Survey, 2020).

Constituents were selected based on their relevance to salinity and nutrient dynamics in the Heart River Basin. For evaluating salinity, TDS, dissolved ions (sulfate, sodium, chloride, potassium, calcium, magnesium, bicarbonate), silica, sodium adsorption ratio (SAR), and physical properties (specific conductance and pH) were selected. For nutrient dynamics, ammonia, nitrate plus nitrite, filtered (dissolved) phosphorus and unfiltered (total) phosphorus were selected (table 2).

Sites were initially selected based on available data for the selected water-quality constituents. Sites that had fewer than 10 samples collected between 1970 and 2020 were removed. Some sites were colocated because data were collected by different agencies or groups at the same location and had different site identification numbers depending on the agency. Where appropriate, the data from colocated sites were combined. Six sites in the basin had colocated sites where data were combined (site 5; site 10; Heart River near Carson, N. Dak. [site 11]; site 18; Heart River at Stark Bridge near Judson, N. Dak. [site 20]; and site 22; table 1).

Many of the selected constituents had censored values or values that represent concentrations too low to be accurately quantified (Oblinger Childress and others, 1999). The censoring level can change through time and among laboratories because of changes in analytical methods, sensitivity of laboratory equipment, or dilutions during laboratory analysis. Censored values are an important consideration in statistical analysis, especially for trends and loads analysis (Helsel and others, 2020). For comparability of data and consistent statistical analysis, constituents with multiple censoring levels were recensored to a common value, and the newly censored dataset was used for all analyses. When choosing the common censoring level, the censoring level that was most consistently observed, or most recently observed, and was considered the most likely to represent an actual measured concentration was selected. When values were censored at a value greater than the common censoring level, they were left unchanged because R–QWTREND (Vecchia and Nustad, 2020) and R–LOADEST (Runkel and others, 2004; Runkel, 2013) can handle multiple censoring levels. The remaining censored values were recoded to the common censoring level, and uncensored values less than the common censoring level were recoded as censored values at the common censoring level. For most constituents, a higher censoring level was chosen because the lower censor levels were in the earlier period of the data and likely did not represent an actual measured concentration. For example, chloride data had censoring levels ranging from 0.1 to 30 milligrams per liter (mg/L), but 3 mg/L was the most frequently and recently observed censoring level. The laboratory analysis diluted the samples, causing the samples to have a censored level greater than 3 mg/L, and therefore a higher censored level was not used to recensor the data (table 2).

Selected constituents and sites varied depending on the analysis because of different data requirements for the various models and analyses. R–QWTREND was used for trend analysis, which requires at least 60 water-quality samples distributed among seasons, a complete daily streamflow record, and preferably no more than 25 percent of the dataset containing censored values (Vecchia and Nustad, 2020). Two main-stem sites (sites 1 and 5; fig. 1) and two tributary sites (sites 10 and 22; fig. 1 and table 1) had sufficient TDS, sulfate, sodium, chloride, potassium, calcium, magnesium, and SAR data for trend analysis related to salinity. Two trend periods were evaluated based on data availability: a recent period (1999–2019) and a historical period (1974–2019). All four sites were used for the recent period and only sites 1, 5, and 22 had sufficient data for the historical period. For example, site 2 appears to have met the requirements for water-quality data availability, but a complete streamflow record was not available, so site 2 was not selected for trend analysis (fig. 1 and table 1). Only two main-stem sites (sites 5 and 22) had sufficient data for trend analysis of nitrate plus nitrite and total phosphorus (fig. 1 and table 1). Based on data availability, nitrate plus nitrite and total phosphorus could only be analyzed for the recent period (1999–2019).

R–LOADEST (Runkel and others, 2004; Runkel, 2013) was used for computing loads and yields for a constituent mass balance in the lower Heart River Basin below Lake Tschida. R–LOADEST has different data requirements than R–QWTREND and a different period (2013–20) was used for the mass balance; therefore, different sites were selected for analysis. Like R–QWTREND, R–LOADEST requires a complete daily streamflow record, so site selection was dependent on availability of water-quality data and streamflow data. Sites selected for load estimation and mass balance included site 5, Heart River above Lake Tschida near Glen Ullin, N. Dak. [site 7], and site 20 on the main-stem Heart River and sites 10, 18, and 21 on Antelope, Big Muddy, and Sweetbriar Creeks, respectively, that are tributaries to the Heart River (fig. 1 and table 1). Loads were estimated for TDS, sulfate, sodium, chloride, and total phosphorus at the selected sites. High variability and prevalence of censored data prevented the estimation of loads for nitrate plus nitrite.

PHREEQC was used for inverse geochemical modeling on two selected reaches of the Heart River based on locations of sites with available data. One reach (model zone 1) included site 5 and Heart River above Lake Tschida near Glen Ullin, N. Dak. [site 6] and another reach (model zone 2) included sites 20, 21, and 22 (fig. 1 and table 1). Constituents selected for use in the PHREEQC model included sulfate, sodium, chloride, potassium, calcium, magnesium, bicarbonate, silica, pH, and water temperature to model geochemical reactions. Although inclusion of iron and aluminum data would have resulted in more accurate geochemical models, too few data were available for those constituents at the selected sites.

Descriptive Statistics

Descriptive statistics were computed for all sites in the Heart River Basin with at least 10 samples collected between 1970 and 2020 to describe the spatial variability of concentrations in the basin. Although 10 or more samples were collected at the sites, the number of values for specific constituents varied by site because samples were collected by various agencies or groups for different purposes. Statistics included minimum, maximum, and the 10th, 25th, 50th (median), 75th, and 90th percentiles of values for individual constituents at each site. Median concentrations for selected constituents were plotted on a map of the Heart River Basin to show spatial patterns in concentration across the basin.

Trend Analysis

R–QWTREND, a publicly available software package developed by USGS for analyzing trends in stream water quality, was used to evaluate water-quality trends (Vecchia and Nustad, 2020). R–QWTREND is a parametric time-series model for streamflow and concentration data (Vecchia and Nustad, 2020). The theory and parameter estimation for the time-series model are described in detail in several publications (Vecchia, 2003, 2005; Vecchia and Nustad, 2020). R–QWTREND can run the time-series model, determine model fit, verify trend models, and produce an output for interpreting and evaluating the results (Vecchia and Nustad, 2020). One of the model’s variables, FRVAR, was designed to capture as much natural flow-related variability in logarithmically transformed concentrations as possible. FRVAR is a function of variables called flow anomalies that depend on concurrent and antecedent streamflow (Vecchia and Nustad, 2020). Flow anomalies address the relation between a constituent concentration and concurrent and lagged streamflow at annual, seasonal, and daily time scales. In addition, the periodic functions of sine and cosine are used to capture the seasonal variation in concentration that is not captured by the flow anomalies (Vecchia and Nustad, 2020). Because of the variable hydroclimate in the basin, characterizing flow-related variability in concentrations at multiple time scales is important because concentrations of water-quality constituents interact with the streamflow in complex ways that cannot be captured using a simple regression between concentration and concurrent streamflow. By accounting for the natural flow-related variability, the ability to detect trends in concentration that are independent of streamflow and climate variability is greatly increased (Vecchia, 2003).

The trends detected by R–QWTREND indicate long-term (10 or more years) changes in annual “flow-averaged” geometric mean concentrations (FAGMCs) that are unrelated to streamflow changes from year to year (Vecchia and Nustad, 2020). To determine the significance of the trend results, the generalized likelihood ratio test statistic was used and is described in Vecchia and Nustad (2020). Three levels of significance were used for this study: a probability (p) value of less than or equal to 0.01 was considered significant; a p-value between 0.01 and 0.05 was considered mildly significant; and a p-value greater than 0.05 was considered nonsignificant. A small p-value indicates that a real trend was detected by R–QWTREND. For example, for a trend with a p-value less than 0.05, the chance that this trend could have occurred given the null (no trend) model that the flow-adjusted concentrations were trend free would be less than 5 percent. A nonsignificant trend indicates that a trend was inconclusive given the available data (Helsel and others, 2020), which does not mean the data have no trend, but that the data may have a trend that is too small to be detected in relation to the natural variability of the data.

Censored values, or values less than the method detection limit for which an exact value is not known (Oblinger Childress and others, 1999), need to be considered during trend analysis (Helsel and other, 2020). R–QWTREND handles censored values and will estimate a value, but it is recommended that no more than 25 percent of the data set be censored values. For the current (2022) study, to include as many sites and constituents as possible in the analysis, sites with constituents having as much as 55 percent censored data were analyzed for trends. Results for sites and constituents analyzed with more censored data than recommended should be interpreted with caution. Chloride, nitrate plus nitrite, and total phosphorus all had censored values with multiple censoring levels and were recensored to a common censoring level (table 2). After recensoring, some sites had more than 50-percent data censored and were not analyzed. Censoring levels for chloride ranged from 0.1 to 30 mg/L (table 2). A higher censoring level of 3 mg/L was used to recensor chloride data with the lower censoring level of 0.1 mg/L to a common censoring level of 3 mg/L. Values with censoring limits higher than 3 mg/L were left as is for R–QWTREND. After recensoring the chloride data, site 10 had more than 50 percent censored data, and therefore chloride trends were not analyzed for this site. For nitrate plus nitrite, censoring levels ranged from 0.005 to 0.1 mg/L (table 2). Nitrate plus nitrite concentrations were recensored to 0.03 mg/L for trend analysis. After recensoring the nitrate plus nitrite data, site 22 had about 50 percent of the data censored and site 5 had about 55 percent of the data censored, but trends were analyzed because of the lack of sites with available data for nitrate plus nitrite trend analysis. For total phosphorus, censoring levels ranged from 0.004 to 0.031 mg/L. Total phosphorus concentrations were recensored to 0.02 mg/L for trend analysis. After recensoring the total phosphorus data, sites 5 and 22 had less than 50 percent of the data censored, but trends were analyzed because of the lack of sites with available data for total phosphorus trend analysis.

Two trend periods were evaluated in this report: a historical period (1974–2019) and a recent period (1999–2019). Streamflow and water-quality data from 1995 to 2020 were used for the recent period and data from 1970 to 2020 were used for the historical period. Data were available for evaluating historical trends at sites 1, 5, and 22 and for recent trends at sites 1, 5, 10, and 22. Historical and recent trends were analyzed for selected constituents including TDS, sodium, sulfate, chloride, potassium, calcium, magnesium, and SAR. Only recent trends were analyzed for nitrate plus nitrite and total phosphorus owing to data availability.

Geochemical Modeling

Inverse geochemical modeling was performed using PHREEQC to simulate geochemical reactions related to increasing salinity (for example, dissolution of calcite [CaCO₃] or gypsum [CaSO₄•2H₂O]; Parkhurst, 1995; Parkhurst and Appelo, 1999). PHREEQC inverse modeling is a mole-balance formulation that solves a set of linear equalities that account for the changes in the number of moles for each ion based on the expected minerals present in the system (Parkhurst, 1995; Parkhurst and Appelo, 1999, 2013). The calculations of the mole-balance formulas are based on the stoichiometry of the potential chemical reactions that have been identified as reacting along the selected flow path (Parkhurst and Appelo, 1999). PHREEQC uses several equations and inequality constraints that allow for uncertainty in the data to calculate the mole transfers (Parkhurst and Appelo, 1999). Specific equations and inequalities used within the PHREEQC program are detailed in several publications (Parkhurst, 1995; Parkhurst and Appelo, 1999). The inverse modeling approach determined the chemical reactions that had occurred in the basin by working backwards from current conditions to historical conditions. Parkhurst and Appelo (1999) contains additional information regarding the PHREEQC inverse model.

Uncertainty terms were included in the model for each aqueous solution and could also be used for each element in the suite of aqueous solutions. These uncertainty terms were related to the charge balance of a solution and adjust the analytical data by the defined percentage to create a charge balance or were related to the analytical uncertainty of a particular constituent (Parkhurst and Appelo, 1999). For example, if an aqueous solution had a charge balance of plus or minus 20 percent, the uncertainty term for this solution cannot be less than 0.2, meaning the analytical data are adjusted within 20 percent for each constituent to reach a suitable charge balance. These uncertainty terms manifest in the model results as residuals. The PHREEQC inverse model produced a residual value that indicated how many parameters and how much the model adjusted the input concentrations to create the balance (Parkhurst and Appelo, 1999). For example, a lower sum of residuals indicated that the data were adjusted less than a model with a higher sum of residuals.

PHREEQC requires a thermodynamic and mineral database for comparison to chemical reactions. The phreeqc.dat database (Parkhurst and Appelo, 1999) was used for the Heart River Basin. Many of the minerals previously identified to be present in the Heart River Basin were already included in the phreeqc.dat database, but the minerals of mirabilite, thenardite, konyaite, biotite, and sodium montmorillonite (Na-montmorillonite) were added (table 3). Inverse modeling only used the stoichiometry of the minerals dissolution reaction, but thermodynamic calculations were performed on thenardite and mirabilite to obtain accurate saturation index (SI) values for these phases. Balanced dissolution reactions were required for all the added minerals.

A list of minerals used for the models were identified from natural sources such as soil and geologic formations, and from anthropogenic sources such as fertilizers used in agricultural practices (table 3). Although cation exchange processes were likely occurring owing to the high cation exchange capacity of clays in the local geology (Langmuir, 1997), these processes were not simulated in the PHREEQC model because the models would not converge when cation exchange processes were included. Minerals present in the Ludlow, Cannonball, and the Bullion Creek Formations, located in the lower Heart River Basin, were identified through previous mineralogical and chemical analyses (Fenner, 1974; Jacob, 1975; Brekke, 1979). Previous soil chemistry studies outlined the possible minerals that could be present in soils across North Dakota (Keller and others, 1986a). A necessary mineral input was to assume mineral precipitation or dissolution, and this assumption was made based on other geochemistry information at each site. To determine whether individual minerals were most likely precipitating or dissolving at a site, a dissolution plot was used in which paired constituent concentrations in millimole per liter were plotted against one another and compared against a dissolution line. The dissolution line was based on the stoichiometry of the dissolution reaction for a given mineral, which is the ideal dissolution of the given mineral. For example, in the dissolution of calcite (CaCO₃ + 2H⁺ = Ca²⁺ + 2HCO₃^-), the ratio of moles [Ca²⁺] to [HCO₃^-] is 1:2. The ratio of 1:2 represents the expected slope of 0.5 for the dissolution of calcite. If dissolution is occurring, when bicarbonate (HCO₃^-) concentrations in millimole per liter are plotted against calcium (Ca²⁺) concentrations in millimole per liter, the points should fall along the dissolution line. Also, SI values were calculated using PHREEQC to test whether the mineral would be dissolving or precipitating for a given mineral in each model. SI values greater than zero indicated precipitation because the solution was supersaturated with respect to that mineral and SI values less than zero indicated dissolution owing to the solution being undersaturated with respect to the given mineral.

Two model zones were delineated to evaluate geochemical changes in salinity in the Heart River Basin. Each model zone was delineated along bedrock geology to understand the control geology has on geochemical changes in salinity. Two model periods for each model zone were selected based on the trend analysis: period 1 included 1974–99, and period 2 included 1999–2019. Model zone 1 was delineated as the reach between site 5 and site 6 upstream from Lake Tschida, and model zone 2 was delineated for the reach between site 20 and site 22 that includes Sweetbriar Creek (site 21; table 1). Model zone 1 was the farthest upstream reach and is underlain by the Bullion Creek Formation (fig. 2A), whereas model zone 2 is underlain by the Ludlow and Cannonball Formations (fig. 2A).

For each model zone and model period, a single sample from each site that was considered most representative of the water-quality changes from the beginning to the ending of the model period was needed to simulate the geochemistry along the reach. Several considerations were made when selecting a representative sample for each site and model period. The most important consideration for selecting a representative sample was the suite of constituents collected for a given sample. For a sample to be included it needed to have results for sulfate, sodium, calcium, magnesium, chloride, potassium, and bicarbonate. Ideally, each sample would also have results for dissolved silica. Another consideration was selecting two samples that closely resemble the results from the trend analysis. Trend analysis results from sites 5 and 22 were used to select samples for model zones 1 and 2, respectively, that showed similar trends for each constituent used in trend analysis. For example, if the trend results at site 22 showed a 45-percent increase in sulfate and a 46-percent increase in chloride between 1999 and 2019, selected samples from sites 21 and 22 might have a 60-percent increase in sulfate and a 50-percent increase in chloride. Although the percentages of each individual constituent did not exactly match the trend result, the sample was overall considered to be representative of the trends observed in the basin. To eliminate seasonal effects, samples for the two sites were selected from the same season. The model periods determined the range of samples from which to select, but the selected samples did not need to come from the specified years at the beginning and end of the model periods but needed to have changes similar to the trend analysis results. Groundwater water-quality data were not considered for the PHREEQC modeling because there were not enough data for the aquifers of interest. This selection method limited the bias in the sample selection process and did not induce an artificial trend when providing representative samples.

Results of the PHREEQC inverse modeling produces multiple models for each model zone and period. Because PHREEQC inverse modeling produces multiple models, a criterion was needed to determine which models were reasonable and unreasonable. Comparison of mole transfers among models was used to ensure that the model was converging on a reliable solution. Models with mole transfer results that were one or more order of magnitude greater than the others were deemed unreasonable. These unreasonable models fit the solution but do not provide accurate information about the model zone and period.

Load Estimation and Mass Balance Analysis

To evaluate what portions of the basin contribute to the mass of selected constituents in the Heart River downstream from Lake Tschida (fig. 1), a mass-balance was computed from estimated loads (mass per time) at selected sites. Constituent load (L) is a function of the volumetric rate of water passing a point in the stream, known as streamflow (Q), and the constituent concentration within the water (C). Constituent yield is a function of constituent load and the drainage area contributing to flow at the site.

The program R–LOADEST (Runkel and others, 2004; Runkel, 2013) was used to estimate constituent loads using regression methods. These methods use natural logarithm (ln)-transformed relations between Q and C to estimate daily C (or L) for a particular constituent at a site (Cohn and others, 1989; Cohn and others, 1992; Cohn, 1995). The regression method can account for non-normal data distributions, seasonal and long-term cycles, censored data, biases associated with using logarithmic transformations, and serial correlations of the residuals (Cohn, 1995). Regression methods use discrete water-quality samples often collected during several years and a daily streamflow hydrograph. A regression model for estimating constituent load can be expressed as the following equation:

In this model, relations between streamflow and load are identified by the β₁ coefficient, temporal trends are identified by β₂, and seasonal effects are identified by β₃ and β₄. Transforming the results of the model from logarithmic space to linear space was accomplished using an adjusted maximum likelihood estimator (AMLE) (Cohn and others, 1992). The AMLE method can also handle censored values. For this application, the variables included in the load model (eq. 1) were selected by using various diagnostic statistics provided by R–LOADEST to evaluate the significance of the variables to include in the model for each site and constituent combination. The explanatory variables were considered statistically significant if the p-value was less than 0.05 for the t-statistic. In using the AMLE method, normal distribution (normality) of the dataset is assumed. The validity of the normality assumption for the residuals was examined using the Turnbull-Weiss likelihood ratio normality statistic (Turnbull and Weiss, 1978). If the p-value from the Turnbull-Weiss statistic was less than 0.05, the residual plots were examined for homoscedasticity (equal statistical variances) and normality. There were some cases where the p-value was less than 0.05, but the AMLE method was used because the dataset contained censored data. As a measure of how much variability in the dependent variable is explained by the independent variable and the R–LOADEST regression equation, R² values were computed and expressed as a percentage. The R² value is a number, 0–1, that when multiplied by 100 is interpreted as the percentage of the variability in the dependent variable explained by the independent variable(s) and the regression equation (Helsel and others, 2020). Generally, a larger R² value indicates a better relation. For example, an R² of 100 percent indicates that all the variability in the dependent variable is explained by the independent variable(s). However, a large R² value does not guarantee the relation is useful (Neter and others, 1996). For example, if estimates require extrapolation outside of the observed independent variables, the estimates may not be accurate. Unless constituent concentrations were highly variable, R² values were expected to be large for the R–LOADEST models because the dependent variable in the R–LOADEST models (constituent load) is a function of one of the independent variables (streamflow).

As a measure of uncertainty in the load estimates, the standard error of prediction (SEP) was provided in R–LOADEST output (Runkel and others, 2004). To compare uncertainty among sites with large differences in loads, the SEP was expressed as a percentage of the total estimated load during the 7-year period for each site and constituent (table 4).

Loads were estimated for selected dissolved ion constituents related to salinity concerns including TDS, sulfate, sodium, and chloride. Loads also were estimated for total phosphorus to evaluate the mass balance of nutrients. High variability and prevalence of censored data prevented the estimation of loads for nitrate plus nitrite. The selected load models, calibration periods, R², and SEP values for each site and constituent are given in table 4.

Loads were computed for sites selected by location and data availability. Site 5 was selected to represent loads in the Heart River Basin upstream from Lake Tschida. For the mass balance analysis, loads were computed for sites 7, 10, 18, 20, 21, and 22. Site 7, located just downstream from Lake Tschida, represents the upper end of the reach for mass balance analysis. Site 22 represents the lower end of the reach for analysis. Sites 10, 18, and 21 represent tributary inputs to the Heart River (fig. 4). Loads for site 20 were computed to evaluate the change in constituent mass in the reach of the Heart River between tributaries for TDS, sulfate, sodium, and chloride. Loads were not computed at site 20 for total phosphorus because no data were available for this constituent.

To avoid long-term changes in streamflow, a short period of 2013–20 was used for developing load models for the selected constituents and sites, and estimation of loads was based on concentration data availability. Data from 2004 to 2020 were used for developing load models for site 7 because there were fewer sample data compared to other sites, although loads were only computed for the 2013–20 period. Constituent concentration data were not collected at site 7, so concentration data collected in Lake Tschida (fig. 1) were used with outflow discharge data from Heart Butte Dam (not shown) to build the load models for the site. Water-quality samples at Lake Tschida near Glen Ullin, N. Dak. [site 23] were collected near the outlet structure in the lake at the midpoint of the total reservoir depth. It was assumed that the water quality near the outlet structure in Lake Tschida would be similar to the outflow location at site 7, although no data were available to evaluate the comparison. For total phosphorus, data from 2004 to 2020 were also used for building the load models for all the selected sites because of the high number of censored data. Estimated loads were only computed for the period of 2013–20.

Monthly, annual, and total loads were estimated for selected sites and constituents for the period of 2013–20. Yields in tons per year per square mile were also calculated for each site by dividing the annual loads by the drainage area (in square miles) contributing to flow at the sampling site for each constituent. Total phosphorus yields were computed as pounds per year per square mile.

A simplistic mass balance was computed for the lower Heart River Basin using the reach on the Heart River between site 7 and site 22 (fig. 4) and the total load for the period of 2013–20. The total load at site 22 was assumed to be a sum of the load from sites 7, 10, 18, and 21, and any intervening flow. Intervening flow includes any input of flow not accounted for by measured tributary inflow and could include groundwater inflow, irrigation return flow, local runoff, and smaller unmeasured tributaries. Not all processes that can affect the mass of constituents in the Heart River were accounted for in the mass balance analysis. Although most of the dissolved ions are generally conservative (the mass of a constituent at an upstream site will be transported to a downstream site with minimal loss of mass), nutrients generally are not considered to be conservative because many processes can affect nutrients such as biological activity and atmospheric exchange (Hem, 1985).

Total Dissolved Solids, Dissolved Ions, Sodium Adsorption Ratio, and Silica

In the Heart River Basin, TDS is more likely to reflect patterns in sulfate, sodium, and bicarbonate because these generally represent a majority portion of the dissolved constituents in a water sample. In contrast, changes in chloride and potassium, which are present in smaller concentrations, have less effect on TDS concentrations. Total dissolved solids, a measure of the sum of all dissolved constituents such as sulfate, sodium, chloride, calcium, magnesium, potassium, bicarbonate, and many other constituents present in small amounts, can be used as an indication of salinity. Most sites had median sulfate and sodium concentrations greater than about 400 and 200 mg/L respectively, compared with most sites having median concentrations of chloride and potassium less than about 30 and about 12 mg/L, respectively (figs. 8 and 9; table 1.1).

The spatial patterns of median concentrations in the basin tended to be similar among TDS, sulfate, sodium, chloride, and bicarbonate with the lowest median concentrations in the South Branch Heart River near South Heart, N. Dak. (site 4) and the highest median concentrations at one of three sites, Heart River near South Heart, N. Dak. (site 2), Hailstone Creek (site 17), or site 18 (table 1.1; figs. 8 and 9). Median TDS concentrations ranged from 510 mg/L at site 4 to 2,150 mg/L at site 17, sulfate concentrations ranged from 170 mg/L at site 4 to 1,160 mg/L at site 17, median sodium concentrations ranged from 120 mg/L at site 4 to 485 mg/L at site 2; median chloride concentrations ranged from 4.8 mg/L at site 4 to 30 mg/L at site 17; and median bicarbonate concentrations ranged from 199 mg/L at site 4 to 606 mg/L at site 18. Spatial patterns among calcium, magnesium, and potassium concentrations were similar; some of the highest concentrations were also observed at site 17 and some of the lowest concentrations at site 4 (fig. 9). At many of the sites, the median concentrations are less representative because fewer observations or years of data are available compared to other sites. Hailstone Creek (sites 12 and 13), Sims Creek (site 14), Big Muddy Creek (site 15), and Wilson Creek (site 16) are located on smaller tributaries and had too few data for a statistical summary.

For one-third of the selected sites in the Heart River Basin with SAR data, median values were relatively high in relation to usability for irrigation water. The SAR is a relation between concentrations of sodium, calcium, and magnesium in the sample and is indicative of the suitability of irrigation water (Richards, 1954; Hem, 1985). In North Dakota, it is recommended that water for continuous irrigation should have SAR values less than 6, and some sites in the Heart River Basin had SAR values greater than 6 (Scherer and others, 2017). Five sites out of 15 had median values for SAR that were greater than 6 (table 1.1), meaning that for these five sites, about one-half of the measured SAR values were less than 6 and the other one-half were greater than 6. Of the five sites, two were in the lower basin in the Big Muddy Creek Basin (sites 17 and 18) and the other three sites were in the upper basin, with one on the Heart River (site 2) and two on tributaries (sites 3 and 4) upstream from E.A. Patterson Lake. The SAR values exceeded 6 in 25 percent of the measured values (the 75th percentile) from 1970 to 2020 at seven sites (sites 2, 3, 4, 17, 18, 21, and 22; table 1.1) in the basin, of which three were in the upper basin and four in the lower basin. At 11 sites (sites 1, 2, 3, 4, 5, 6, 17, 18, 20, 21, and 22), 10 percent of the SAR values (the 90th percentile) exceeded 6, with six in the upper basin and five in the lower basin.

Median silica concentrations were low compared to other constituents, ranging from 3.2 mg/L at site 11 to 9.1 mg/L at site 4 (table 1.1). Silica data were not available for sites 12–17. Generally, main-stem sites had lower concentrations of silica than tributaries.

Spatial patterns in TDS and dissolved ions in the basin may be related to evaporation, groundwater effects from geology, and lake processes in Lake Tschida. The process of evaporation may contribute to higher concentrations within the basin, particularly in tributaries and upper basin sites that have lower volumes of water. The concentration of a constituent is determined by dividing its mass by the volume of water. Evaporation reduces the volume of water and, therefore, produces higher concentrations. Differing bedrock across the basin may contribute to the spatial water-quality patterns observed in the Heart River Basin. The Ludlow and Cannonball, Bullion Creek, and Sentinel Butte Formations have aquifers that contribute flow to streams in the Heart River Basin. The Ludlow and Cannonball Formations and Bullion Creek aquifer systems have been defined as sodium-bicarbonate dominated water (Naplin and Shaver, 1978; Randich, 1979; Ackerman, 1980), whereas the Sentinel Butte aquifer system has been defined as a sodium-sulfate type water (Trapp and Croft, 1975). Comparing these aquifer systems, the Sentinel Butte system has higher dissolved ion concentrations than the Ludlow and Cannonball and Bullion Creek systems (Naplin and Shaver, 1978; Randich, 1979; Ackerman 1980). Groundwater can have a major effect on a stream’s water quality during low flow or base-flow conditions (Smakhtin, 2001; Brutsaert, 2008). Site 7, downstream from Lake Tschida, generally had lower dissolved ion concentrations than site 5, upstream from Lake Tschida, which indicates that processes within the lake are reducing dissolved ion concentrations (figs. 8 and 9). Using the average storage volume for Lake Tschida (60,737 acre-feet) and average discharge (12 cubic feet per second) between 1999 and 2019 (Bureau of Reclamation, 2020), residence time of Lake Tschida was estimated to be 253 days. This residence time is sufficient to allow the lake processes to have some effect on the concentration disparity in dissolved ion concentrations between upstream and downstream sites. A potential process could be the settling of negatively charged clay particles, which occurs because the water does not turn over as frequently, and this allows the adsorption of positively charged cations to these clay particles (Langmuir, 1997). Reservoirs have been observed elsewhere in the United States to reduce concentrations in dissolved ion concentrations between upstream and downstream sections and reduce the variability in concentrations in the downstream section (Hem, 1985).

Nutrients

Spatial patterns for median concentrations of nitrate plus nitrite differed from the other nutrient constituents. The lowest median nitrate plus nitrite concentrations were below the recensored reporting level of 0.03 mg/L at multiple sites, and the highest median concentration was 0.8 mg/L at Antelope Creek (site 8; fig. 10 and table 1.2). Spatial patterns for median concentrations in total ammonia, dissolved phosphorus, and total phosphorus were similar with the lowest median concentrations at site 11 and the highest median concentrations at site 4 and site 15. Median concentrations of total ammonia ranged from less than 0.03 mg/L at many sites to 0.09 mg/L at site 4; median dissolved phosphorus concentrations ranged from less than 0.02 mg/L at four sites (sites 5, 10, 11, and 22) to 0.1 mg/L at site 4; and median total phosphorus concentrations ranged from less than 0.02 mg/L at sites 11 and 20 to 0.51 mg/L at site 15 (fig. 10 and table 1.2). At many of the sites, the median concentrations were less representative because fewer observations or years of data were available to compare to other sites, such as North Creek near South Heart, N. Dak. (site 3) or site 11. Sites 6, 7, and 20 had too few data for a statistical summary.

A large percentage of censored values in the nutrient data, in many cases 50 percent or more, made it difficult to compare concentrations spatially. Many of the median values in the statistical summaries for the four nutrient constituents of ammonia, nitrate plus nitrite, dissolved phosphorus, and total phosphorus were censored values. Dissolved phosphorus had the least amount of data of the selected nutrient constituents.

Nutrient patterns in the Heart River Basin may be related to land use differences, energy production, and lake processes in Lake Tschida. Land use differences in the basin can contribute to the higher median nutrient concentrations because of different practices on the landscape that can contribute nutrients to streams. Some sources of nutrients include runoff from agricultural areas, where fertilizers are applied or livestock production occurs, and runoff from urban areas, where fertilizers are applied to lawns, shrubs, and trees (Hem, 1985). As described above, the upper basin has higher land use of cultivated crops and hay/pasture than the lower basin (fig. 2B). The Antelope Creek Basin and Big Muddy Creek Basin (not shown) both had sites with the highest median concentrations for some of the nutrients (fig. 10), and these basins have 63 and 36 percent land use for cultivated crops and hay/pasture, respectively (fig. 2B). Energy production in North Dakota has increased dramatically since about 2008 (U.S. Energy Information Administration, 2013). Energy production in the Heart River Basin is generally only in the upper basin near Dickinson, N. Dak. The burning of fossil fuels is a source of nitrogen in the atmosphere (Hem, 1985), which could be an atmospheric source of nitrogen in the basin.

Like dissolved ions, processes in Lake Tschida can reduce the concentrations of nutrients in the basin. Lake processes such as eutrophication and adsorption of nutrients can decrease nutrient concentrations in the basin. Algal blooms, a result of eutrophication, have been observed in Lake Tschida. Algae use the nitrogen and phosphorus as a nutrient source, which can result in reduced the concentrations dissolved in the water column (North Dakota Department of Environmental Quality, 2021b). The adsorption of positively charged nutrients to negatively charged clay particles can occur because the particles settle out of suspension when they reach the lake (Langmuir, 1997). Reservoirs elsewhere in the United States were observed to reduce nutrient concentrations between upstream and downstream sections and reduce the variability in concentrations in the downstream section (Hem, 1985). Because many of the sites had censored data, it was difficult to observe the effects Lake Tschida had on the spatial pattern of nutrient concentrations in the Heart River Basin.

Physical Properties

The spatial pattern of specific conductance was similar to TDS, and the median specific conductance values ranged in the Heart River Basin from 641 microsiemens per centimeter at 25 degrees Celsius (µS/cm) at site 4 to 2,320 µS/cm at site 2 on the Heart River (table 1.3). Specific conductance is an indirect measure of the collective concentration of dissolved ions in solution (U.S. Geological Survey, variously dated), and can be used to evaluate the salinity hazard of irrigation water (Wilcox, 1955). Many of the same processes that can affect TDS and dissolved ion concentrations will also affect the specific conductance of a sample. Basinwide median specific conductance values were greater than 1,000 µS/cm at all sites except site 4.

No spatial pattern for pH was evident in the basin, and the median pH value in the basin ranged from 8.1 at site 4 to 8.4 at three sites (sites 11, 18, and 20; table 1.3). The pH of water is defined as the negative base-10 logarithm of the hydrogen ion activity and is used as an indicator of the health of any given aquatic system (U.S. Geological Survey, variously dated). River waters generally have a pH range between 6.5 and 8.5 and median values for pH in the Heart River Basin fall in this range (Hem, 1985).

Total Dissolved Solids, Dissolved Ions, and Sodium Adsorption Ratio

During the historical trend period, TDS concentrations have increased since the mid-1970s through 2019 (fig. 11, table 5). For all sites, a one-trend model was used. For site 1, a significant 22-percent increase in FAGMC was detected between 1974 and 2019 where the FAGMC increased from 590 to 723 mg/L (table 5). The significant trend detected at site 5 was a 23-percent increase where the FAGMC increased from 1,100 to 1,350 mg/L (table 5). Finally, at site 22 a significant 44-percent increase was detected in which the FAGMC increased from 830 to 1,200 mg/L (table 5).

During the recent trend period, increasing concentrations in TDS were observed across the Heart River Basin, and the magnitude of the increases was smaller at tributary sites compared to main-stem sites (table 5). All sites had similar percentage changes in FAGMC of TDS except for site 10, which had a nonsignificant 0.5-percent increase with a FAGMC increase in TDS of 3 mg/L. Percentage changes at the other sites ranged between 15 (site 5) and 26 percent (site 1), but the FAGMC increases were larger in the main-stem sites (sites 5 and 22). Both increases in FAGMC on the main-stem sites were greater than 150 mg/L, whereas the tributary sites had increases less than 150 mg/L (table 5).

During the historical trend period, sulfate and sodium concentrations have increased since the mid-1970s at all sites, but the increase was greater on main-stem sites from 1999 to 2019 than from 1974 to 1999 (table 5, figs. 12–13). For site 1, a one-trend model was used for sulfate and sodium and significant 46- and 30-percent increases, respectively, were detected from 1974 to 2019 (table 5). For site 5, a two-trend model was used for sulfate and sodium and a nonsignificant 4-percent increase for sulfate and a 16-percent mildly significant increase for sodium was detected from 1974 to 1999, but a larger significant increase of 37 percent for sulfate and 30 percent for sodium was detected from 1999 to 2019 (table 5). Sulfate and sodium at site 22 had a mildly significant 27- and significant 34-percent increase detected from 1974 to 1999, respectively (table 5). A larger significant increase in sulfate and sodium of 44 and 36 percent, respectively, was detected during 1999–2019 (figs. 12 and 13, table 5). For sulfate and sodium, the percentage increases on the main-stem sites for the 20-year period between 1999 and 2019 were nearly equivalent to the increase for site 1 during the 45-year period. Taking the 45-year period into account, during 1974–2019, the FAGMC nearly doubled at site 22 for sulfate (351 to 641 mg/L) and sodium (161 to 293 mg/L) (table 5).

During the recent trend period, increases in sulfate and sodium concentrations were detected basinwide, and sulfate trends had larger increases in concentration compared to sodium (table 5). For sulfate and sodium, the largest increases in magnitude and percentage change occurred on the main-stem sites (sites 5 and 22). Site 10 had the smallest increases in magnitude and percentage change for sulfate and sodium, but site 1 had a similar percentage change to the main-stem sites (table 5). However, the FAGMC increase in sulfate and sodium at site 1 was smaller than the changes on the main-stem sites. The FAGMC increases for sulfate and sodium at sites 5 and 22 were nearly identical, but site 5 had a higher starting and ending concentration compared to site 22. The large disparity in concentrations between site 1 and site 5 in the upper basin indicated that increases in both constituents were likely occurring in other Heart River tributaries upstream from site 5.

Land use, agricultural practices, and hydroclimatic changes in the Heart River Basin, combined with naturally occurring and readily available sulfate and sodium in geologic formations, soils, and aquifers, likely contributed to the moderate to large increases of sulfate and sodium concentrations. A primary land-use change in the basin during the recent trend period was the conversion of “considerable acreage” from flood to pivot irrigation (Chad Skretteburg, Western Heart River Irrigation District, written commun., 2021). Pivot irrigation can cause increases in soil salinity because there is not a large enough volume of water applied to leach the salts through the soil column (Al-Ghobari, 2011; Shahid, 2013). The salts that form in the Heart River Basin are sodium-sulfate evaporites (Keller and others, 1986a), and during runoff events these evaporites dissolve and flow into the Heart River owing to the high solubility of the sodium-sulfate evaporites (Keller and others, 1986b). Hydroclimate changes from a dry to wet climatic period around 1995 in the basin likely contributed to the increasing concentrations of sulfate and sodium because increases in base flows were observed basinwide (figs. 5–7). Although R–QWTREND does remove streamflow variability when computing the trends, it may not fully capture these streamflow changes in the basin (Nustad and Vecchia, 2020). Increases in base flow correspond to increases in groundwater discharge, and sulfate and sodium are naturally occurring in the geology of the basin (Trapp and Croft, 1975; Naplin and Shaver, 1978; Randich, 1979; Ackerman, 1980). Concentrations of sulfate and sodium likely increase in the streams owing to the increasing discharge from the aquifers. With the change to a wet climatic period, increasing groundwater levels will have more contact with sodium-sulfate evaporites that were in the previously unsaturated zone, increasing the concentrations in the groundwater. The NRCS has implemented conservation practices on about 300 acres for salinity and sodicity management in the basin for saline or sodic soils since 2006 (Rita Sveen, NRCS, written commun., 2021). Two common methods are used to treat saline and sodic soils: one is the application of elemental sulfur to form gypsum and the other is the direct application of gypsum (Franzen and others, 2014; Diaz and Presley, 2017). In both methods, the dissolution of the gypsum is used to replace the sodium ions that have adsorbed to the soil particles and flush the system of sodium ions (Franzen and others, 2014; Diaz and Presley, 2017).

Chloride concentrations showed large, consistent, and significant increases starting in the mid-1970s through 2019, whereas potassium concentrations remained mostly constant with some small fluctuations over the same period (figs. 14 and 15, respectively; table 5). From the one-trend chloride model at site 1, a significant 84 percent increase was detected from 1974 to 2019, where the FAGMC increased from 5.6 to 10 mg/L. From the two-trend potassium model at site 1, a significant 32 percent increase was detected from 1974 to 1999 and a mildly significant 15 percent decrease was detected from 1999 to 2019. Considering both trend periods for potassium at site 1, the FAGMC remained largely unchanged (5.6 to 6.3 mg/L; table 5). At site 5, the two-trend chloride model detected a significant 30 percent increase in chloride concentrations from 1974 to 1999 and a significant 87 percent increase from 1999 to 2019. Considering both trend periods for chloride at site 5, the FAGMC more than doubled from 13 to 32 mg/L during 1974–2019 (fig. 14; table 5). At site 5, the three-trend potassium model detected fluctuating concentrations with a significant 40 percent increase from 1974 to 1984, a mildly significant 11 percent decrease from 1984 to 1999, followed by a significant 12 percent increase from 1999 to 2019 (table 5). Considering the three trend periods for potassium at site 5, overall, the FAGMC slightly increased from 7.9 to 11 mg/L. At site 22, the same significant increase of 41 percent was detected for chloride for both trend periods, and the FAGMC of chloride nearly doubled between 1974 and 2019 from 10 to 19 mg/L (table 5).

During the recent period, large increases (greater than 40 percent) in chloride concentrations were detected, and potassium concentrations were mostly constant, although small (0.9 mg/L or less) decreases on tributaries and small (1.3 mg/L) increases on the main-stem sites were detected (table 5). Although percentage increases in chloride were large, actual changes in concentration seem small because FAGMCs of chloride were less than about 30 mg/L (table 5). For site 10, recent trends for chloride were not analyzed because of too many censored values. Site 1 and 22 had similar percentage changes for chloride (46 and 44 percent, respectively) and FAGMC increases (2.9 and 5 mg/L, respectively, table 5). The increase in chloride for the main-stem site upstream from Lake Tschida (site 5) was substantially larger (14 mg/L) than the increase at the most downstream main-stem site (5 mg/L; site 22), which indicates larger increases in chloride in the upper basin tributaries compared to the lower basin tributaries (table 5). The FAGMC of potassium generally hovered around 10 mg/L at all sites, and the tributary sites had small nonsignificant decreasing concentrations (15 percent or smaller) compared to the main-stem sites that had small significant (15 percent or smaller) increasing concentrations (table 5). Overall, the FAGMC changes in potassium at all sites ranged from −1.2 to 1.3 mg/L and the percentage changes ranged from −12 to 14 percent. The small FAGMC and percentage changes were indicative that the potassium concentrations were remaining mostly constant basinwide.

Chloride concentrations had larger increases between 1999 and 2019 compared to 1974–99, which were likely a result of sources such as deicing and dust-control chemicals, agricultural management and practices, and energy production (Granato and others, 2015). Many roadways are located within the basin, notably Interstate 94 (fig. 1), that may have had deicing chemicals applied to road surfaces during the winter, which may have contributed to increased trends in chloride (Granato and others, 2015). Dust-control chemicals, such as calcium and magnesium chloride, used for gravel and oil field roads likely affected the chloride concentration increases in the basin. Agricultural fertilizers with potash (potassium chloride) are used in North Dakota, and an average of 40 percent of wheat, corn, and soybean fields in the State have potash applied (U.S. Department of Agriculture, 2019). Runoff from hay/pasture lands in the Heart River Basin also likely contributed to the increasing chloride concentrations from animal waste (fig. 2B; Granato and others, 2015). Documentation from the Western Heart River Irrigation District stated that conversion from flood to pivot irrigation has been occurring in the Heart River Basin (Chad Skretteburg, Western Heart River Irrigation District, written commun., 2021). These practices may produce more salts on the soil surfaces, which can include chloride salts (Al-Ghobari, 2011; Shahid, 2013). Many of the salts that would be precipitated on the surface are highly soluble when in contact with water and would likely be dissolved and transported during runoff events. Energy production in North Dakota has increased dramatically since about 2008 and could be a possible source for increasing chloride concentrations (U.S. Energy Information Administration, 2013). Energy production in Heart River Basin is generally only in the upper portion of the basin, mostly in Stark and Dunn Counties. Saline waters or brines are produced as part of the extraction process. Historically, production water was disposed of using evaporation pits until 2011 when it was mandated that all production waters either be disposed of in deep injection wells or be recycled (North Dakota Administrative Code, 2020). Releases of brines into the environment can occur from the transport of the brines to disposal sites, pipeline breaks, vehicle accidents, and leakage from historic evaporation pits.

Potassium concentrations did not follow the trends observed with other constituents and remained mostly constant from 1974 to 2019, likely because of stable processes during the trend periods. These processes include, but are not limited to, agricultural fertilizer application, weathering of biotite and potassium feldspars into clays, and possible biological processes (Langmuir, 1997; Franzen and Bu, 2018). The Heart River Basin is a predominantly agricultural basin, and 40 percent of wheat, soybean, and corn fields in North Dakota have potash applied (U.S. Department of Agriculture, 2019). Naturally occurring potassium is present in the geologic formations but is likely a minor contributor to the potassium concentrations (Langmuir, 1997).

Calcium and magnesium concentrations increased since the mid-1970s at all sites (table 5, figs. 16–17, respectively), except for a decrease at site 5 between 1974 and 1999. At site 1, one-trend models were used for calcium and magnesium and a nonsignificant 4-percent increase in calcium and a significant 29-percent increase in magnesium was detected from 1974 to 2019. A two-trend model was used for site 5, and nonsignificant decreases in calcium and magnesium concentrations of 11 and 5 percent, respectively, were detected from 1974 to 1999 (table 5). At site 5 during 1999–2019, a significant 34-percent increase in calcium and a significant 54-percent increase in magnesium was detected (figs. 16 and 17, respectively, table 5). At site 22, a two-trend model was used for calcium and magnesium and increases in concentrations were detected for both trend periods. For calcium, a nonsignificant 8-percent increase was detected from 1974 to 1999, and a significant 24-percent increase was detected from 1999 to 2019 (table 5, fig. 16). For magnesium, a mildly significant 22-percent increase was detected from 1974 to 1999, and a significant 49-percent increase was detected from 1999 to 2019 was detected (table 5, fig. 17).

During the recent period, calcium and magnesium concentrations increased basinwide (table 5). Sites 1, 10, and 22 all had similar percentage increases (13, 17, and 14 percent, respectively) for calcium concentrations, whereas the percentage change at site 5 was nearly double (32 percent) compared to the other three sites (table 5). FAGMC from the first to last year of calcium started and ended at similar concentrations at sites 5, 10, and 22. The FAGMC at site 1 was one-half that of the other three sites. Site 5 had the largest percentage change and FAGMC increase in calcium in the basin, which indicates that the upper basin tributaries had larger increases compared to the lower basin tributaries. At the Heart River sites (sites 5 and 22), the magnesium percentage increases (52 and 46 percent, respectively) were larger than on the tributary sites (sites 1 and 10), with the increases of 39 and 13 percent, respectively (table 5). The FAGMC magnesium increases at sites 1 and 10 were similar (6 and 7 mg/L, respectively; table 5), however site 10 has starting and ending concentrations that were more than double the concentrations at site 1 (table 5). Main-stem sites (sites 5 and 22) had similar starting and ending magnesium concentrations and similar FAGMC increases (26 and 21 mg/L, respectively; table 5).

Overall, increasing concentrations of calcium and magnesium were consistent between the upper and lower basin and are likely related to anthropogenic effects and naturally occurring geologic effects. Anthropogenic effects in the basin include dust-control measures on gravel roads or agricultural practice changes. Calcium and magnesium chloride chemicals are used to control dust on gravel and oil field roads and can be a source of increased concentrations in streams (Granato and others, 2015). The conversion from flood to pivot irrigation is the major agricultural change in the basin and can cause the accumulation of salts on the surface and in the shallow subsurface (Al-Ghobari, 2011; Shahid, 2013). Changing climatic conditions in the basin include the increase in base flow that was observed across the basin (figs. 5–7). The increasing base flow likely affects the increasing concentrations in calcium and magnesium because both are present in the aquifer materials (Trapp and Croft, 1975; Naplin and Shaver, 1978; Randich, 1979; Ackerman, 1980). The most common soluble salts in soils in the basin were sodium-, calcium-, and magnesium-sulfates (Keller and others, 1986a). Dissolution of these calcium- and magnesium-sulfates owing to wetter conditions was likely contributing to the increasing concentrations.

Results of the historical trend analysis of SAR values detected significant increases basinwide since the mid-1970s (fig. 18, table 5). Site 1 had a significant 19-percent increase where the FAGMC increased from 4.1 to 4.9 from 1974 to 2019 (table 5). Site 5 had a significant 25-percent increase where the FAGMC increased from 4.6 to 5.8 from 1974 to 2019 (table 5). Site 22, the farthest downstream site on the Heart River, had a significant 35-percent increase in SAR, where the FAGMC increased from 4.3 to 5.8 from 1974 to 2019 (table 5). The magnitude of the increases in SAR was greater in the lower basin at site 22 (35 percent) compared to the upper basin at site 1 (19 percent) (table 5).

Increasing values of SAR were observed basinwide during the recent trend periods (table 5). Both tributary sites (sites 1 and 10) had nonsignificant increases and both Heart River sites (site 5 and 22) had significant increases. The percentage change of the increases was similar at sites 1, 5, and 22 (13, 10, and 17 percent, respectively) and site 10 had a smaller percentage increase (4 percent). The smallest increase occurred at site 10 on Antelope Creek, but overall by the end of the trend period the tributary sites (sites 1 and 10) had lower SAR values than main-stem sites (sites 5 and 22) and Heart River sites had values approaching the recommended SAR value threshold of 6 (Scherer and others, 2017).

Because SAR is computed as the ratio of sodium to calcium and magnesium in the soil water, increasing concentrations in these constituents resulted in increasing SAR values. SAR is an indicator used in the management of irrigation water along with TDS and specific conductance. In North Dakota, the recommendation is that for continuous irrigation SAR should remain below 6 mg/L (Scherer and others, 2017). Poorer quality water may be used for irrigation on land that is irrigated sporadically, which is defined as 1 year out of every 3 or more years (Scherer and others, 2017). The small to moderate trends observed in SAR were a function of sodium, magnesium, and calcium each having positive trends. SAR is related to all of those constituents and with each having a different positive trend, the effect on the overall SAR was reduced. If the changes occurred in only one or two of the constituents, it would have a greater effect on the SAR trend. A difference in the SAR trends compared to the other constituents was that site 22 had the largest magnitude increase. Site 5 consistently had the largest increases with the other constituents, but not the largest increase in SAR. Results of the trend analysis at site 5 and 22 were similar for sodium, calcium, and magnesium (table 5). Sodium concentrations were similar at sites 5 and 22, but site 22 had lower concentrations for calcium and magnesium, which corresponded to larger SAR values.

Nutrients

Unlike TDS and dissolved ion concentrations, significant increases in nutrient concentrations were not detected from 1999 to 2019, and nitrate plus nitrite concentrations most likely decreased upstream from Lake Tschida (table 6). Because of limited data availability, only nitrate plus nitrite and total phosphorus at two sites for the recent trend period could be analyzed. Additionally, because 55 percent of the nitrate plus nitrite data for site 5 were censored and 41 percent of the nitrate plus nitrite data were censored for site 22, results should be interpreted with caution (table 1.2). No significant changes in nutrient concentrations were detected in the Heart River Basin during the recent trend period (table 6). Site 5, upstream from Lake Tschida, had a 47-percent decrease in nitrate plus nitrite concentrations (table 6). From visual inspection of the nitrate plus nitrite data for site 5 (not shown), most of the censored values occurred toward the end of the trend period, supporting the trend result of decreasing concentrations. Nonsignificant changes do not necessarily mean that a trend is not present, but variability in the data combined with a small trend can make a trend undetectable. Results from site 22, the most downstream site, indicate that the nutrient concentrations did not have significant changes in the lower basin.

Trend results did not indicate increasing concentrations of nitrate plus nitrite and total phosphorus concentrations on the two Heart River sites (sites 5 and 22), but anthropogenic changes, conservation practices, and land use can affect nutrients in the basin. When nutrients become overabundant in lakes and reservoirs, harmful algal blooms can form (Paerl and others, 2001). Anthropogenic changes in nutrient concentrations can be related to crop management, livestock, fertilizers, land-use changes, geology, and industrial or municipal effluents (Hem, 1985; Tornes and Brigham, 1993; Dubrovsky and others, 2010). Crop and livestock management practices can also affect nutrient concentrations through the application of fertilizer and runoff from feed lots (Hem, 1985; Tornes and Brigham, 1993). Conservation practices are implemented to limit the amount of nutrients being transported to streams and lakes. Between 2006 and 2020, conservation practices to reduce nutrients were applied to about 128,000 acres throughout the Heart River Basin (Rita Sveen, NRCS, written commun., 2021). The upper basin has a higher percentage of agricultural land use than the lower basin does (fig. 2B), which likely corresponds to more acres in the upper basin with conservation practices being implemented. Site 5, located in the upper basin, had decreasing nitrate plus nitrite concentrations, which may be related to conservation practices implemented in the upper basin. Nutrient dynamics are complex and may vary over time and location. Increased energy production in the western portion of the basin may contribute to the nitrate plus nitrite concentrations in the basin (U.S. Energy Information Administration, 2013) because the burning of fossil fuels, such as coal or petroleum products, release nitrogen into the atmosphere (Hem, 1985).

Model Zone 1

Inputs for the PHREEQC inverse model for model zone 1 (Heart River reach from site 5 to site 6; fig. 1) included minerals, dissolution or precipitation assumptions, and uncertainty terms. The 11 minerals included in model zone 1 for both model periods were mirabilite, thenardite, konyaite, gypsum, sylvite, calcite, dolomite, quartz, K-feldspar, Na-montmorillonite, and illite (Keller and others, 1986a; Fenner, 1974; Jacob, 1975; Brekke, 1979; Granato and others, 2015). Within the shales and sandstones of Bullion Creek Formation that underlies model zone 1, gypsum, calcite, dolomite, quartz, K-feldspar, Na-montmorillonite, and illite are present (Brekke, 1979). From the dissolution plots at site 5, it was determined that mirabilite, thenardite, konyaite, sylvite, calcite, and dolomite were most likely dissolving (fig. 19A–F). The SI values for mirabilite and thenardite from PHREEQC of approximately −5 and −7, respectively, indicated dissolution of these minerals. The SI values for sylvite were less than zero, indicating dissolution, and the dissolution plot for sylvite also indicates dissolution as the data plots close to the line (fig. 19D). The sample selection for model zone 1 was determined by the trend analysis completed at site 5. Two samples were selected for sites 5 and 6 for each model period (table 7). The uncertainty term specified for period 1 was 2.5 percent and for period 2 was 10 percent. The uncertainty terms for each model were low, reflecting the quality of the methods used to collect the data.

From PHREEQC inverse modeling for period 1 (1974–99) and model zone 1, eight reasonable models indicated that the clay mineral-water interactions and dissolution of evaporites control the geochemistry (table 8); however, some models have lower residuals (for example, model 5 has a model residual of 2.6), indicating less adjustment of selected water-quality samples, so all eight models were taken into consideration when interpreting geochemical changes. Hydrolysis of K-feldspar to clays in the Bullion Creek Formation was determined to be a control on potassium concentrations as seen in the positive mole transfers of K-feldspar in five models and the negative mole transfers of illite, which occurred in all eight models (table 8). The dissolution of sylvite, as seen in every model with 0.18 millimole per kilogram of water (mmol/kg H₂O) of sylvite dissolving, indicates that agricultural fertilizers were another source of potassium and likely the main source of chloride in the system. The amount of konyaite, mirabilite, and thenardite dissolving during this period ranged from 0.50 to 1.16 mmol/kg H₂O. Dissolution of konyaite, mirabilite, and thenardite was the major source of sodium, sulfate, and magnesium in the system because one of these minerals was dissolving in 7 of 8 models. Dolomite and calcite were sources of magnesium, calcium, and dissolution of these minerals, which occurred in 4 out of 5 models, accounting for the increases in concentrations of magnesium and calcium in the Heart River. The amount of gypsum precipitating ranged from 0.81 to 2.2 mmol/kg H₂O, and precipitation occurred in 7 of 8 models (table 8). From trend analysis, the sulfate concentration was unchanging and calcium concentrations were slightly decreasing during this period, which was most likely caused by the precipitation of gypsum.

Results of the inverse modeling for period 2 (1999–2019) for model zone 1 had eight reasonable models that indicated that the dissolution of evaporites was the major geochemical control (table 8). Three major differences were observed in period 2 compared to period 1: (1) gypsum was dissolving in five models, (2) dissolution of sylvite was the dominant source of potassium over hydrolysis of K-feldspar (reduction of approximately 98 percent in K-feldspar mole transfers), and (3) evaporite (konyaite, mirabilite, and thenardite) dissolution was 50 percent greater compared to period 1. Dissolution of gypsum, konyaite, mirabilite, and thenardite increased from an average of about 0.9 mmol/kg H₂O in period 1 to about 1.45 mmol/kg H₂O in period 2 and was the major source of calcium, sodium, sulfate, and magnesium (table 8). Dissolution of dolomite, ranging from 0.67 to 1.72 mmol/kg H₂O, was another major source of magnesium in the system and occurred in five of the models (table 8). Precipitation of calcite occurred in four of the models and was in response to the overabundance of bicarbonate and calcium from the dissolution of dolomite and gypsum. Sylvite dissolution increased from 0.18 to 0.20 mmol/kg H₂O and was the major source of potassium and chloride (table 8). Hydrolysis of K-feldspar was no longer a dominant source of potassium, and this explains why the potassium trends discussed in the “Water-Quality Trends for Selected Sites” section decreased during the later periods in the three-trend model at site 5 (fig 15).

Model Zone 2

Inputs to the inverse geochemical model for model zone 2 required inputs for minerals, precipitation or dissolution assumptions, and uncertainty terms. The 11 minerals included in model zone 2 for both periods were mirabilite, thenardite, konyaite, gypsum, sylvite, calcite, dolomite, quartz, biotite, Ca-montmorillonite, and illite (Keller and others, 1986a, 1986b; Fenner, 1974). The Ludlow and Cannonball Formations underlies model zone 2, and calcite, dolomite, quartz, biotite, Ca-montmorillonite, and illite are present in these formations (Fenner, 1974). From the dissolution plots at site 22, mirabilite, thenardite, konyaite, and sylvite were most likely dissolving (fig. 20). Additional information supported that thenardite and mirabilite were assumed to be dissolving in these models owing to these minerals having a high solubility and each having a SI value less than zero (table 7). No dissolution or precipitation assumption was made for konyaite owing to the mineral having the ability to be stable under multiple conditions that could be present in this system (Keller and others, 1986b). Further investigation determined that calcite and dolomite were likely precipitating owing to SI values greater than zero for the selected samples (table 7). The PHREEQC SI value for sylvite was less than zero and supports the dissolution plot for sylvite (table 7). Sample selection for this model zone was determined from the trend analysis at site 22; two samples were selected from sites 20, 21, and 22 and were used to develop the PHREEQC inverse models (table 7). The sample selected at site 5 was collected in 1998 but was determined to be the most representative sample for the starting water-quality conditions for this model. Uncertainty terms used for both model periods were 7.5 percent to account for the charge balance of the samples.

Results of the geochemical modeling for period 1 (1974–99) produced seven reasonable models, and the geochemical control of the system was the dissolution of the sulfate evaporite minerals of mirabilite, thenardite, konyaite, and gypsum (table 9). Dissolution of either mirabilite, thenardite, or konyaite occurred in each model with amounts of gypsum dissolution less than the other three evaporites. These sulfate evaporite minerals were the likely sources of sulfate, sodium, calcium, and magnesium in the system and likely contributed to the increasing concentrations. Minor components in this system included the conversion of biotite (positive mole transfers in all models) to illite (negative mole transfers in all models), the dissolution of sylvite (all models), the weathering of Ca-montmorillonite (all models), and precipitation of calcite or dolomite (2 of 7 models; table 9).

Geochemical modeling results for period 2 (1999–2019) produced 11 reasonable models that were also controlled by the dissolution of sulfate evaporite minerals (table 9). Mole transfers increased between period 1 and period 2 for the sulfate evaporites of mirabilite, thenardite, and konyaite. Mole transfers for these three sulfate evaporites ranged from 1.1 to 1.7 mmol/kg H₂O during period 1 with a combined average of about 1.40 mmol/kg H₂O and ranged from 0.79 to 2.9 mmol/kg H₂O during period 2 with a combined average of about 2.00 mmol/kg H₂O (table 9). Mole transfers of gypsum decreased from about 0.80 to about 0.55 mmol/kg H₂O between period 1and period 2. Sulfate evaporites likely contributed to increasing concentrations of sulfate, sodium, calcium, and magnesium. Another difference of note is that the dissolution of sylvite decreased from period 1 to period 2 from about 0.32 mmol/kg H₂O to about 0.08 mmol/kg H₂O (table 9). Minor controls on the system for model zone 2 during period 2 included agricultural fertilizers (dissolution of sylvite), conversion of biotite (positive mole transfers in 10 models) to illite and montmorillonite (negative mole transfers in 10 models), and the physical weathering of clay minerals (four models).

Model Comparison

Differences between models in the upper and lower basin (model zones 1 and 2, respectively) indicated that geology controls some of the water-quality changes in the Heart River Basin. Carbonate dissolution was occurring in the upper basin, whereas these minerals generally precipitated out in the lower basin. The dissolution of carbonates in the upper basin likely oversaturated the river in these minerals, causing the precipitation of carbonates in the lower basin. The varying composition of sandstones produced different chemical reactions between zones in the basin. Hydrolysis of K-feldspar had a larger effect in the upper basin compared to the chemical weathering of biotite in the lower basin. Groundwater was not considered because of too few data; therefore, geochemical evolution through the aquifers was not modeled.

Sulfate evaporite minerals in soils in the Heart River had the most geochemical control over the system. In the upper and lower basin, large mole transfers of sulfate evaporites likely were a result of natural climate variations and irrigation practices. An increase in base flow beginning in the 1990s was likely caused by a water table increase, which changes the geochemical evolution of the water (figs. 5–7). In the unsaturated zone, evaporites and other minerals are precipitated and left behind during the infiltration of water. With a groundwater table rise, these evaporites and other minerals are in contact with water, which can begin to dissolve the minerals, increasing the dissolved ion content of the groundwater (fig. 21). The increase in base flow corresponds to an increase in groundwater discharge into the Heart River Basin, which can increase the loads of the dissolved solids in the streams (fig. 21). The conversion from flood to pivot irrigation can create soil salinity issues at the surface owing to less water infiltrating the soil column (Al-Ghobari, 2011; Shahid, 2013). To manage the development of salts on the surface, more water can be used in irrigation water application to flush the salts below the root zone where accumulation generally happens (Shahid, 2013).