Spatial and Temporal Variations in Phosphorus Loads in the Illinois River Basin, Illinois USA

,


INTRODUCTION
Phosphorus (P) enrichment of surface waters has been a long-standing problem (Schindler 1971) contributing to eutrophication and harmful algal blooms. Billions of dollars have been invested to reduce point and nonpoint sources of P, and yet the problem persists in much of the United States (U.S.). Reasons for this persistence vary with the specific circumstances within watersheds, such as changes in precipitation, river flow, point-source discharges, legacy P, agricultural practices, and river chemistry (Muenich et al. 2016;Jarvie et al. 2017). Stackpoole, Sabo, Falcone, and Sprague (2021a) and Stackpoole, Sabo, and Falcone (2021b) reported flow normalized total P (TP) loads in the lower Mississippi River increased between 1995 and 2005, despite reductions in the net P balance of the watershed. They hypothesized that factors that caused the increased P load included implementation of no-till, changes in watershed buffering capacity, the effects of tile drainage, or increased precipitation.
Increased yields of TP and soluble reactive P in three rivers flowing into western Lake Erie were associated with increased rainfall amounts, intensity, discharge, tile drainage, and adoption of conservation tillage systems (Jarvie et al. 2017;Williams and King 2020). In Iowa, average TP yields from 46 agricultural watersheds were highly correlated with land slope and storm flows, with the average yields increasing from 1.2 kilograms P per hectare-year (kg P/ha/year) during 2000-2005 to 1.5 kg P/ha/year during 2014-2017 (Schilling et al. 2020). TP loads leaving the Illinois River Basin (at Valley City) have increased 30% from the periods 1980-1996(IEPA et al. 2021) despite large reductions in point source discharges and riverine TP concentrations in upper tributaries draining the Chicago metropolitan area (Pluth et al. 2021).
In addition to changes in TP concentrations and loads, the upper Illinois River is affected by high chloride concentrations from a combination of wastewater discharge and road salt (Kelly et al. 2012). Increased chloride concentrations increase ionic strength, which is thought to shift the phosphate (PO 4 ) equilibrium toward adsorption (McDaniel et al. 2009). On the other hand, chloride and other ions reduce the solubility of oxygen (Sherwood et al. 1991), and chloride solutions are used in laboratory analyses to extract PO 4 from soil and sediment samples, where their capacity to desorb depends on several factors such as ionic strength, pH, and temperature (Wang and Li 2010). Haq et al. (2018) incubated river sediments with varying concentrations of sodium chloride and reported a significant increase in PO 4 release from the sediments at 7 of 12 sites. Causes for the differences among sites were not determined and the authors called for more research on this topic.
Sulfate concentrations in Illinois rivers have decreased as atmospheric deposition of sulfate declined following reductions in sulfur dioxide (SO 2 ) emissions from coal-burning power plants . Decreased sulfate concentrations reduce the ionic strength, which may offset the increased ionic strength from increased chloride concentrations. Reduced nitrate concentrations in the Illinois River have been reported (McIsaac et al. 2016), which would also reduce the ionic strength of the river water. Additionally, sulfate and nitrate can serve as electron acceptors in microbial respiration, thereby influencing redox potential. As nitrate becomes scarce, the redox potential in sediments can decline to a level that promotes PO 4 desorption from sediments (Olila and Reddy 1997). Sulfate reduction occurs at redox potentials less than those that promote PO 4 desorption (Pan et al. 2019). Reductions in sulfate concentrations and loads in rivers are thought to be primarily due to reduced deposition from the atmosphere. But if reduced sulfate concentrations are also a consequence of lower redox potential in sediments reducing sulfate to hydrogen sulfide, then we should also expect these conditions to release PO 4 from sediments.
Furthermore, because PO 4 is the limiting nutrient in most freshwater ecosystems (Schindler et al. 2008), any increase in PO 4 may increase primary productivity. Increases in primary productivity, in turn, lead to higher rates of respiration, decomposition, and depletion of dissolved oxygen (DO). The anoxic and or hypoxic conditions that result are harmful to fish and other animals and release P bound to bottom sediments, thereby forming a feedback loop (Correll 1999). A small increase in PO 4 may trigger larger and larger releases of PO 4 until the sediment reservoir is depleted of P or until oxygen is restored by an increase in water-column mixing.
The objectives of this study were to identify and quantify factors contributing to the increased TP loads at Valley City, Illinois, as well as variations in loads throughout the Illinois River Basin. Extensive and long-term monitoring in the basin was used to quantify the spatial and temporal variations in loads in subwatersheds and their relations to watershed characteristics.

Geographic Setting: Illinois River Basin
The Illinois River Basin drains approximately 44% of the land in Illinois, and smaller portions of Indiana, Michigan, and Wisconsin ( Figure 1). For this study, we define the upper basin as the drainage area upstream from the monitoring location at Marseilles on the Illinois River. The upper basin comprises 30% of the Illinois River Basin drainage area upstream of Valley City, including the Chicago, Des Plaines, and Kankakee River Basins. The Chicago and Des Plaines River Basins include extensive urban areas and high population density associated with the city of Chicago, Cook County, and the surrounding suburbs (Table S1). Approximately 5.2 million people reside in Cook County (2,100 people/km 2 ), which has remained relatively constant from 1980 to 2019, while the population density of its neighboring counties increased from 239 to 410 people/km 2 during that time (U.S. Census Bureau 2021). Land cover in 2019 in the upper basin was 25% developed and 60% planted/cultivated, with much of the cultivated land planted to corn and soybeans in the Kankakee River Basin (Table S1).
We define the lower basin as the drainage area downstream of the monitoring station at Marseilles, where land cover was 7.7% developed and 76% planted/cultivated (Dewitz and USGS 2021). The population density was relatively low and steady over the study period averaging 36 people/km 2 , but with large spatial variation ranging from Peoria County with 115 people/km 2 to a few counties with less than 10 people/km 2 . Much of the cropland is underlain by subsurface ("tile") drainage systems, especially in the Vermilion, Mackinaw, and Sangamon River Basins (McIsaac and Hu 2004). Irrigation is uncommon except for an area with sandy soils adjacent to the Illinois River near the confluences with the Mackinaw and Sangamon Rivers (Illinois State Water Survey 2015).
The Illinois River downstream from Marseilles has a very low gradient (about 0.007%), and extensive sediment deposition in the lower river has been reported by Demissie et al. (2016). Seven tributaries to the Illinois River downstream from Marseilles are monitored, which we refer to collectively as the lower FIGURE 1. Illinois River Basin, major subwatersheds and monitoring locations.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
tributaries. Combined, the drainage areas of the upper basin and lower tributaries cover 80% of the Illinois River Basin drainage area upstream from Valley City, Illinois. We refer to the remaining 20%, which includes the segment of the Illinois River between the monitoring locations at Marseilles and Valley City, as the lower mainstem subwatershed, which also includes 13,700 km 2 of drainage area. Macoupin Creek joins the Illinois River downstream from Valley City and thus does not contribute to the TP load at Valley City but was included in the study to more fully characterize the southern portions of the Illinois River Basin.

DATA AND METHODS
Monitoring locations (n = 41) used in this study were selected based on the availability of daily discharge and at least 60 concentration observations collected in proximity during most years from 1978 through 2019 (Table 1). Sites 12 and 19 were used in our initial study, but do not appear in Table 1 because of insufficient data at these sites during 1989-1996 that was necessary for the analysis presented here. An additional U.S. Geological Survey (USGS) streamgage on the Illinois River at Henry,  (Table S1) of these watersheds are based on National Land Cover Data 2019 (Dewitz and USGS 2021).
Concentrations of constituents examined in this study (Table 2) were mostly obtained from the State of Illinois Ambient Water Quality Monitoring Network (AWQMN) coordinated by the Illinois Environmental Protection Agency (IEPA). Samples were typically collected every 40 days (9 samples per year) with some variation. Very few samples were collected in 2008, and loads were not calculated for that year for most sites. Most AWQMN data are available from the National Water Quality Monitoring Council Water Quality Portal (WQP; http://www. waterqualitydata.us), but at the time of this study, some historical data were still being migrated from the State's internal database to WQP. These missing data are limited to the 1998 water year and were available from the STORET Legacy Data Center (https://www3.epa.gov/storet/legacy/gateway.htm). For a few locations, notably the Illinois River at Valley City and Marseilles, the USGS had also collected samples and provided concentration data through the WQP. For most sites and constituents approximately 300 concentration values were obtained, but the number varied due to differences in sampling frequency and the first year of concentration measurement, which ranged from 1977 to 1985. The number of samples used for each annual load estimate is reported in Hodson (2021).
Riverine constituent loads were estimated using Weighted Regressions on Time, Discharge, and Season (WRTDS) with Kalman filtering (WRTDS_K; Zhang and Hirsch 2019). The basic WRTDS method accounts for the seasonal and long-term changes in relations between river nutrient concentrations and river discharge (Hirsch et al. 2010(Hirsch et al. , 2015. The most recent version of this approach included a residual adjustment inspired by Kalman filtering (Kalman 1960) and was adopted for this study because of its superior accuracy compared to other methods (Lee et al. 2019).
Using WRTDS_K, daily mean discharge and constituent concentrations from periodic water-quality samples were used to estimate the daily loads of each constituent. Ten versions of the WRTDS_K model were generated by resampling the data using a block bootstrap with a block span of 200 days (Hirsch et al. 2015). For each of these re-estimated models, 20 Kalman-filtered time series of daily values were computed, each using all the sample values in the original dataset, generating a total of 200 simulated time series. The mean of these simulations was used to estimate annual loads. The lag-1 autocorrelation coefficients for all models were set to 0.9. The resulting annual load estimates and associated uncertainty estimate statistics are available as a separate data release (Hodson 2021).
Dissolved P (DP) and particulate P (PP) loads were modeled separately using WRTDS_K. Because DP and PP may have different relations with time, discharge, and season in some settings, we calculated TP loads as the sum of DP and PP loads, rather than TP loads calculated from TP concentrations and discharge. Constituent yields were calculated by dividing loads by drainage area. For watersheds with multiple monitoring locations, incremental loads and yields were calculated for each watershed segment between monitoring locations by calculating the change in load between the downstream and upstream locations and dividing by the difference in drainage areas between the downstream and upstream monitoring locations. Temporal changes in incremental yields were evaluated by comparing the average incremental yields for 1989-1996. The 1989-1996 period is the earliest of the 1980-1996 baseline period adopted by the Gulf of Mexico Hypoxia Task Force and IEPA (IEPA et al. 2021) with the most complete discharge and concentration data across the IRB. The 2015-2019 period was the most recent five-year period for which consistent concentration data were available at the initiation of this study. Precipitation and river flows varied across the basin in these years, and in general 1989-1996 were highly variable with relatively dry (1989) and wet (1993) years, whereas 2015-2019 ranged from average to above average (2019) years. Hydrologic variation was treated as a potential explanatory variable. Annual flow-weighted concentrations were calculated by dividing the annual load by the annual discharge and using unit conversion factors to obtain discharge-weighted concentrations in milligrams per liter (mg/L). Correlations among water discharge and water quality constituents were calculated.

Point Source Data
Point source location and discharge data were obtained from several sources. The Metropolitan Water Reclamation District of Greater Chicago (MWRD) operates seven wastewater treatment plants that serve approximately 90% of the homes and businesses of Cook County. Daily to monthly discharge and concentration data from these facilities from 1983 are available from the MWRD website (https:// mwrd.org/final-effluents).
The Sanitary District of Decatur (SDD) provided discharge and P concentration data for 1989-2019 (Keith Richard, Sanitary District of Decatur, written commun., 2018, 2019, and 2020). SDD treats wastewater from two large grain-processing facilities as well as domestic wastewater from municipalities totaling about 90,000 residents and associated businesses. During 2015-2018, TP discharge from SDD averaged about 750 megagrams P per year (Mg P/year), which was among the highest in Illinois for an individual facility. Water year TP discharge values of MWRD and SDD facilities are in Table S2.
Locations of major dischargers (greater than 1 million gallons per day) and TP loads after about 2010 from other point-source facilities were obtained from the U.S. Environmental Protection Agency (USEPA) ECHO system (Figures S1 and S2). Starting around 2010, IEPA has been increasingly requiring larger facilities to monitor and report P discharges as part of their National Pollutant Discharge Elimination System (NPDES) permit. For those facilities required to monitor and report, discharge monthly reports (DMRs) are made available from the USEPA ECHO system.
We also considered point source discharge estimates compiled by Ivahnenko (2017), but these were often substantially lower than the values from USEPA and MWRD ( Figures S3 and S4). Data from MWRD and IEPA were generally in close agreement. The Ivahnenko (2017) point source estimates are part of a widely cited national dataset used in several studies on factors affecting riverine nutrient loads. Accurate representation of point source discharges over large areas has been an ongoing challenge.

Incremental TP Yields
During 1989-1996, the greatest incremental TP yields tended to occur in and downstream from the more densely populated areas (Figure 2; Table S3). Incremental TP yields ranged from 6.68 kg P ha/year for the Thorn Creek at Thornton (site 13 in Figure 2) to À1.30 kg P/ha/year for the lower mainstem (site 42). A yield less than zero indicates greater incoming than outgoing TP or net deposition of TP in this drainage area. The Illinois River downstream from Marseilles has a very low gradient, and extensive sediment deposition in the lower river has been reported by Demissie et al. (2016). Because much riverine P is attached to sediment, sediment deposition in this stretch of the river may account for much of the net deposition of P.
During 2015-2019, the greatest incremental TP yields still tended to be associated with the more densely populated Chicago region. Annual incremental TP yields ranged from 6.50 kg P/ha/year for the West Branch of the DuPage River at West Chicago (site 14) to 0.24 kg P/ha/year for the lower mainstem (42; Table S3). The lower mainstem had been a net sink of TP during 1989-1996 and became a net exporter of TP during 2015-2019. Even though it is a net exporter, TP may still be accumulating in the lower segment of the Illinois River if the TP yield from the lower mainstem drainage area is greater than 0.24 kg P/ha/year.
The largest increase in incremental TP yield from 1989-1996 to 2015-2019 was the 2.06 kg P/ha/year increase from Kickapoo Creek at Waynesville (site 37 in Figure 3; Table S3), followed by a 1.57 kg P/ha/ year increase in the West Branch of the DuPage River between Warrenville (site 15) and West Chicago (site 14). The Sangamon River between Fisher and Monticello (site 33) had an increase of 1.41 kg P/ ha/year. All three of these subwatersheds receive point-source wastewater discharge from municipalities that have experienced population growth in recent decades. Expanded dairy operations in the upper Sangamon (sites 32 and 33) may contribute to increased P yields there. The portion of the Sangamon River between Monticello and Riverton (site 34) had an increase of 0.99 kg P/ha/year, which was largely influenced by increased point source P discharge from the SDD.
Across all 41 subwatersheds, TP yield during 2015-2019 was weakly correlated (coefficient of determination [R 2 ] = 0.23, p < 0.002) with the percentage of developed land cover in 2019, whereas change in TP yield from 1989-1996 to 2015-2019 was uncorrelated with change in water yield. For 18 subwatersheds with less than 9% developed land cover in 2019, changes in incremental TP yield were weakly correlated (R 2 = 0.37, p < 0.01) with changes in incremental water yield (Figure 4). The 9% threshold was identified through trial-and-error analysis of the data. For the 23 subwatersheds with greater than 9% developed land cover, the change in incremental TP yield from 1989-1996 to 2015-2019 was uncorrelated with changes in the incremental water yield. In these more developed watersheds, changes in TP yields were more influenced by factors other than a precipitation change, such as changes in population and/or wastewater treatment. The largest increases in water yield occurred in watersheds with greater than 9% developed land cover, and developed land cover and change in water yield were correlated (R 2 = 0.55). Water yield decreased in 8 of the 41 subwatersheds.
For the subwatersheds with less than 9% developed land, three outliers were in the lower range of change  Table 1 for monitoring location names).

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
in water yield. When the smallest of these outlier watersheds is removed (site 29, Indian Creek, 162 km 2 ), the R 2 for the relation increases from 0.37 to 0.50. Accurate monitoring and load estimation for small watersheds is challenging because high flow and high load events are relatively short-lived and are sampled infrequently. Additionally, in small watersheds, changes in land cover or agricultural management over a large percentage of the drainage area are more likely than in large watersheds. The reason for the decline in TP yield in Indian Creek remains unknown but worthy of additional investigation.

Upper Basin
Phosphorus loads from the upper basin are dominated by point source discharges from Chicago and the surrounding suburban watersheds ( Figure 5). The more agricultural Kankakee River watershed upstream from Wilmington occupies 62% of the upper basin drainage area but contributed only 15% to 25% of the TP load at Marseilles. Although loads had considerable year-to-year variation from both the Kankakee River and the upper basin, 2015-2019 average TP loads from the Kankakee increased slightly from 1989-1996 to 2015-2019, while loads from the upper basin decreased slightly. Reductions in point source discharges from the MWRD, Thorn Creek Basin Sanitary District, and other facilities in the upper basin appear to have been offset by a combination of (1) increased point source discharge from suburban wastewater facilities, (2) legacy P in the upper tributaries, and/or (3) increased TP load from agricultural areas in the upper basin (e.g., Mazon River, site 17). The lack of an increased TP load at Marseilles indicates that the  Table 1 for monitoring location names.
increased TP load at Valley City originates downstream from Marseilles.

Lower Tributaries and Lower Mainstem Subbasin
The TP load in the Illinois River at Valley City increased 39% from 1989-1996 to 2015-2019 (or 6,892 to 9,599 Mg P/year;  water yield for watersheds with more than 9% developed land cover (n = 18) and those with less than 9% developed land cover (n = 23). Numbers associated with the more extreme points refer to subwatershed outlet ID numbers in Table 1 and Figures 2 and 3.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
limited to remobilization of legacy P that had been deposited in the lower mainstem in previous years. During 1989-1996, the combined TP loads from the upper basin and lower tributaries were greater than the load at Valley City, indicating a net deposition of TP along the lower mainstem. Conversely, the average annual load at Valley City was greater than the sum of the upstream loads during 2015-2019, indicating a net export of P from the lower mainstem drainage area (Table 3).
TP loads in six of the seven lower tributaries increased from 1989-1996 to 2015-2019, with the exception being the Spoon River (site 30), which decreased by 69 Mg P/year (9.5%). The largest increase in a tributary load was 597 Mg P/year from the Sangamon River at Oakford (location 39; Table 3). This increase in load represents about 22% of the 2,707 Mg P/year increase at Valley City and is similar in magnitude to the increase in P discharge from the SDD (590 Mg P/year) during the same time interval (Table S2). The sum of the changes in TP loads from the lower tributaries excluding Spoon River at Oakford (64 Mg P/year) offsets a 62 Mg P/year decrease from the upper basin. After accounting for these contributions, the increase in TP from the lower mainstem (2,108 Mg P/year) represents about 78% of the increase in TP load at Valley City.
Dissolved P (DP) loads increased in all lower tributaries, with the largest increase coming from the Sangamon River (Table 4; location 39). The DP load at Valley City increased 52%, from 3,612 to 5,486 Mg P/year between 1989-1996. During 1989-1996, the average DP load at Valley City was less than the sum of the upstream loads, indicating that a portion of the incoming DP was either deposited, transformed to particulate P load, or lost to groundwater along the lower mainstem. Conversely, the DP load at Valley City was slightly greater than the sum of the upstream loads during 2015-2019.
The combined increase in DP load from the upper basin and lower tributaries (1,023 Mg P/year) was 55% of the increase observed at Valley City. The single largest tributary increase (462 Mg P/year) was from the Sangamon River. The largest percentage increase in DP load (146%) was from Big Bureau Creek, which also had the highest DP yield (1.32 kg Values in bold font represent 1) the sum of loads from upstream (Upper basin + lower tributaries), 2) the load at Valley City, and 3) the difference between 1) and 2). Values in bold font represent 1) the sum of loads from upstream (Upper basin + lower tributaries), 2) the load at Valley City, and 3) the difference between 1) and 2). P/ha/year) among the lower tributaries during 2015-2019, and was only slightly less than the yield from the upper basin (1.44 kg P/ha/year). DP yields from the other lower tributaries were all less than 0.92 kg P/ha/year during 2015-2019. Particulate P (PP) loads at Valley City increased 25% from 3,281 Mg P/year during 1989-1996 to 4,113 Mg P/year during 2015-2019, despite reduced PP loads from the upper basin and every lower tributary except the Sangamon River (Table 5). The increased PP load at Valley City may have been a consequence of mobilization of PP from sediments and/or transformation of DP to PP during the transit along the Illinois River. The PP load from the Sangamon River increased 135 Mg P/year, some of which may have been the result of the transformation of DP discharged as wastewater into PP.
Total suspended solids (TSS) loads decreased at all lower tributaries, in contrast to the upper basin, where TSS increased 8% from 1989-1996 to 2015-2019 (Table 6). The largest reductions in TSS load were from the Spoon River (232 gigagrams per year [Gg/ year]) and the Vermilion River (108 Gg/year), although the TSS yield from the Spoon River remained the highest among all monitoring locations in both time periods: 1.75 Mg/ha/year in 1989-1996 and 1.21 Mg/ha/ year in 2015-2019. TSS yields from the other tributaries were less than 0.85 Mg/ha/year. The largest percentage reductions in TSS loads were in the Vermilion River (49%) and Big Bureau Creek (49%).
The TSS load at Valley City decreased by 14% (423 Gg/year), but during both time periods, the annual average TSS load at Valley City was less than the combined loads from upstream, indicating TSS deposition in the lower mainstem but with less deposition occurring in 2015-2019. Other long-term sediment budgets have also indicated considerable sediment deposition in the mainstem of the Illinois River (Demissie et al. 2016). This accumulated sediment can be a source of legacy P through remobilization of previously deposited sediments and desorption of P from sediments.
Loads of VSS increased at Marseilles and Valley City but decreased in all the lower tributaries (Table 7). Because much of VSS is organic matter on eroded sediment, the increase in VSS load at Marseilles may be due to increased TSS load at that site Values in bold font represent 1) the sum of loads from upstream (Upper basin + lower tributaries), 2) the load at Valley City, and 3) the difference between 1) and 2). Values in bold font represent 1) the sum of loads from upstream (Upper basin + lower tributaries), 2) the load at Valley City, and 3) the difference between 1) and 2). (Table 6). Of the sites in Tables 6 and 7, Valley City is the only location where VSS increased while TSS decreased, which may be the result of increased algal productivity in the lower mainstem. Reduced sediment loads would permit greater light penetration into the upper water column and potentially increase algal production, respiration, and reduced nighttime DO. In both time periods, more VSS came into the lower mainstem than left at Valley City. Whether from increased algal production or organic matter with eroded sediment, the accumulation of organic matter in the lower mainstem can contribute to oxygen depletion, low redox potential, and release of PO 4 from sediments in the river, backwater lakes, and wetlands.

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
The incremental sulfate load in the lower mainstem increased 6% from 1989-1996 to 2015-2019 while sulfate loads from the upper basin and the lower tributaries decreased from 7% to 44% (Table S4). A decrease in sulfate loads in most of the basin is largely due to decreased atmospheric deposition following reduced sulfate emissions from coalfired power plants . The increased incremental sulfate load in the lower Illinois River may be due to a variety of sources such as increased loads from unmonitored point sources and/or increased drainage from livestock wastes and/or legacy sulfate emerging from groundwater.

Major Point Sources to the Lower Mainstem Subbasin
The point-source dischargers required to monitor and report their P discharges to the IEPA in the lower mainstem reported a total of about 259 Mg P/ year discharged on average during 2011-2014 and about 184 Mg P/year during 2015-2019 (Table S5). Reductions in P discharge may be due to reduced population and/or changes in wastewater treatment methods in response to more restrictive NPDES permits. In addition to declining point source P discharge in recent years, these totals were a small percentage of the 2,108 Mg P/year increase in TP load from the lower mainstem. The TP load increase Values in bold font represent 1) the sum of loads from upstream (Upper basin + lower tributaries), 2) the load at Valley City, and 3) the difference between 1) and 2).

Correlation Analysis
Correlation analysis with annual loads from 1989 to 2019 (excluding 2008 because of a lack of concentration data) indicated that the incremental TP load between Marseilles and Valley City was highly correlated with the DP:TP load ratio of the combined loads from the upper basin and lower tributaries (r = 0.66), incremental TSS load (r = 0.65), and annual flow weighted chloride concentration at Valley City (r = 0.61) ( Table 8). These three variables are also correlated with each other. The upstream DP:TP ratio is highly correlated with chloride (r = 0.82) and nitrate (r = À0.73) concentrations at Valley City. Of note, incremental TP and TSS loads were negatively correlated with the average annual discharge at Valley City (r = À0.39 and À0.49, respectively). This may be a consequence of greater deposition in flood plains and backwaters during high flows, whereas moderate and low flows were more constrained to the main channel where deposition is less common.
The high degree of correlation among variables renders it difficult or impossible to identify causation or quantify how much influence each variable may have had. Increased DP from upstream is likely less susceptible to deposition. Increased chloride concentrations may have had a role in desorbing PO 4 from sediments, while reduced nitrate and sulfate concentrations may have allowed lower redox potentials in sediments, thereby releasing PO 4 from sediments. Disentangling these interactions is probably best achieved by controlled experiments and more intensive and dedicated measurements in the lower Illinois River and elsewhere.

DISCUSSION
The increase in DP and reduction in TSS loads in the lower tributaries may be a result of the adoption of conservation practices such as conservation tillage aimed at reducing soil erosion. One potential consequence of many conservation tillage systems, particularly no-till, is that P fertilizer applied to the soil surface is not deeply incorporated into cropland soils. This contributes to greater concentrations of P in the top layers of the soil and in eroded sediment (Zanon et al. 2020). Conservation tillage often reduces sediment and particulate P in surface runoff, but it may increase DP in runoff compared to inversion tillage (moldboard plowing) (McIsaac et al. 1995;Jarvie et al. 2017). The effect of alternative tillage systems on TP loads depends on the relative contributions of DP and PP under conventional and conservation tillage. Additionally, increased installation of subsurface drainage tiles may be a factor in reducing PP and increasing DP (Gentry et al. 2007).
Another factor that may contribute to increased TP loads from agricultural areas is the expansion of concentrated livestock feeding operations (Waller et al. 2021). Although these facilities are required to have a nutrient management plan, the rules governing manure application allow for the accumulation of soil test P in the soil up to 330 kg P/ha. At elevated concentrations, P is more vulnerable to transport in surface runoff and subsurface drainage (Heckrath et al. 1995).
The overall TP yield from the Illinois River at Valley City increased from 1.0 to 1.39 kg P/ha/year from 1989-1996 to 2015-2019. In northwest Ohio, TP yield from the Maumee River watershed remained relatively steady at 1.2 kg P/ha/year, while PO 4 loads increased and were highly correlated with increased river discharge (Williams and King 2020). DP and PO 4 are not identical, but the similarity with the Maumee River is notable, and the differences in TP load warrant further study.
The range of incremental TP yields from the 41 Illinois River Basin subwatersheds (À1.30 to 6.68 kg P/ha/year) is similar to the range Schilling et al. (2020) reported for 46 agricultural watersheds in Iowa: 0.44 to 7.71 kg P/ha/year. The high-yielding watersheds in Iowa, however, were not urban watersheds, but steeply sloping agricultural watersheds, with high rates of upland and river channel erosion. Schilling et al. (2020) did not calculate incremental TP yields and therefore could not report any TP yields <0. Based on their published results, we calculated incremental yields for three watersheds, but none had incremental TP yields <0. These watersheds had average land slopes ranging from 1.8% to 3.0%,

JOURNAL OF THE AMERICAN WATER RESOURCES ASSOCIATION
which is much greater than the 0.007% of the Illinois River in the lower mainstem.
The drainage area of the Illinois River Basin is 2.1% of the Mississippi River Basin. The 2015-2019 average TP load at Valley City (Table 3) averaged about 5% of the Mississippi River loads flowing to the Gulf of Mexico (USGS, n.d.). Annual TP loads at Valley City 1979-2019 are significantly correlated with TP loads to the Gulf of Mexico (r = 0.65). Similar processes may have contributed to increased P loads at both scales.
Based on our results, we recommend further research to quantify the P legacy in Illinois River sediments and to quantify the roles of algal production and river chemistry (e.g., chloride, nitrate, sulfate, DO concentrations) on P transport and desorption of PO 4 from sediments in different settings (see also Haq et al. 2018). Additionally, the role of increasingly concentrated animal production and associated manure P management deserves further study. The potential to reduce riverine P loads by using practices that incorporate manure and fertilizer P into the soil without increasing soil erosion warrants further exploration. Finally, detailed investigations of subwatersheds with unusually large increases or decreases in TP yields could help identify causes for these changes.

SUMMARY AND CONCLUSIONS
Variations in TP loads in the Illinois River Basin were investigated, with a focus on identifying causes of increased loads at Valley City after 1996. Net deposition of P in the lower mainstem subwatershed during 1989-1996 reduced loads at Valley City. This portion of the river system is very flat and known to accumulate sediment. During 2015-2019, the lower mainstem subwatershed had shifted from being a net sink to being a net source of riverine P load, and this shift represented 78% of the increased TP load at Valley City. We did not conclusively identify causes for this shift but found correlations that additional research could investigate for causality: increased DP:TP and chloride concentrations and reduced nitrate and sulfate concentrations from upstream tributaries.
An increased TP load from the Sangamon River at Oakford represented 22% of the increased load at Valley City, and this increase was similar in magnitude to an increase in point source P discharge from the SDD. The combined changes in TP loads from the other monitored tributaries downstream from Marseilles (including a reduction in load from the Spoon River) represented only 2.4% of the increased load at Valley City, which offset a similarly sized load reduction at Marseilles.
In the urbanized portion of the upper basin, reductions in riverine TP loads occurred downstream from major point source facilities operated by MWRD, and Lake County WRD and Thorn Creek water reclamation districts. But these and possibly other load reductions were offset, in part, by increased loads from suburban areas where the population had increased (DuPage River) and from agricultural areas (e.g., Mazon River). Consequently, TP loads from the upper basin changed very little from 1989-1996 to 2015-2019. Changes in incremental TP yields from 1989-1996 to 2015-2019 from 23 subwatersheds with more than 9% developed land cover in 2019 were unrelated to changes in water yield, as these watersheds were more likely to be influenced by point sources. In subwatersheds with less than 9% developed land cover and less likely to be affected by point sources, changes in TP yield were weakly correlated with changes in water yield.

DATA AVAILABILITY STATEMENT
Annual river loads of phosphorus and other constituents estimated using WRTDS-K at each Illinois EPA monitoring site are available in an accompanying data release (Hodson 2021) at https://doi.org/10. 5066/P9JX8YVL.

SUPPORTING INFORMATION
Additional supporting information may be found online under the Supporting Information tab for this article: Land cover characteristics for incremental watersheds, point source discharge information, TP loads and yields for the incremental watersheds, and sulfate loads for the upper basin and lower tributaries.

ACKNOWLEDGMENTS
This work was supported with funding from the Illinois Nutrient Research and Education Council and U.S. Geological Survey, as well as data provided by the U.S. Geological Survey and Illinois Environmental Protection Agency. The work also benefitted from point source discharge data provided by the Sanitary District of Decatur and the Metropolitan Water Reclamation District of Greater Chicago. We also gratefully acknowledge editorial assistance from Lisa Sheppard of the Illinois State Water Survey, and valuable comments provided by the USGS reviewers Janet Carter, Kymm Barnes, and Galen Hoogestraat.