Meteorologic Data

Meteorologic data used in the hydrodynamic model included air-temperature, dew-point temperature, cloud-cover, wind-speed, and wind-direction time series that were obtained from the Automated Surface Observing Site station C09 at Morris Municipal Airport, Morris, Ill., (fig. 1; Iowa State University, 2023) and were available at 20-minute intervals throughout the 2020–22 analysis period. The qualitative Automated Surface Observing Site cloud-cover parameters of “clear,” “scattered,” “broken,” “overcast,” and “vertical visibility” were converted to numerical values for use in CE-QUAL-W2 using values of 0 for clear, 3 for scattered, 6 for broken, 8 for overcast, and 10 for vertical visibility. Hourly short wave solar radiation data were computed by CE-QUAL-W2 using the numerical values for ASOS cloud cover parameters. Daily air-temperature data used in determining water temperature time series for the Fox River at Dayton were obtained from the Global Historical Climatology Network daily data for station Ottawa 4SW (National Centers for Environmental Information, 2023a). Projections of future changes in monthly mean temperatures at Ottawa, Ill., used in the model sensitivity analyses were obtained from the U.S. Climate Resilience Toolkit Climate Explorer (U.S. Federal Government, 2023). These projections used the representative concentration pathway with 4.5 megawatts per square meter (RCP 4.5; vanVuuren and others, 2011) for the year 2075.

The monthly mean air temperatures during the January 2020 to December 2022 study period were similar to the long-term (1980–2010) monthly mean temperatures at the Morris, Ill., monitoring location (National Centers for Environmental Information, 2023b; fig. 3). The greatest deviations from the long-term normals were in the months of January (2020 and 2022) and February (2021) and outside of the expected HAB development period. Monthly precipitation was more variable during the study period compared to the long-term monthly mean precipitation (1980–2010), with several spring months of each year having substantially greater precipitation than the long-term monthly average (fig. 3).

Water-Temperature Data

Continuous and discrete water-temperature data were used for model development, calibration, and validation. Continuous data were collected in 15-minute intervals at Seneca (U.S. Geological Survey, 2022c) and Starved Rock (U.S. Geological Survey, 2022d) stations and every 30 minutes at the Marseilles Lock and Dam pool (John Brammier, USACE, written commun., 2022; fig. 1; table 1). The Seneca time-series data were used as the upstream model boundary condition, whereas the Starved Rock time series were used for model calibration and validation at a fixed near-surface (0.5 m) depth. Additional discrete temperature profile data were collected by USGS personnel near Seneca and Ottawa, Ill., and used for model calibration with depth. The discrete profile data were collected at multiple locations within the study reach between May 22, 2022, and July 15, 2022. Water-temperature time-series data were unavailable for the Fox River and, therefore, water temperatures were developed using a regression model with observed air temperature, daily mean streamflow, water temperature, and day of year normalized to 0–2π (fig. 4) as variables. This relation was developed using 97 observed air and water-temperature data pairs collected at the Fox River at Dayton, Ill. (USGS station 05552500) streamgage between 1986 and 1997, daily mean flow at the Fox River at Dayton streamgage, and daily air-temperature time-series data from the Ottawa 4SW weather station (National Centers for Environmental Information, 2023a). The resulting predicted input water-temperature time series for the Fox River and the observed time series for the Illinois River at Seneca were similar (fig. 5).

Hydrodynamic Model Parameters

In addition to the required model geometry files, many required initial conditions, boundary condition time series, and model parameters were specified to simulate water-surface elevation, streamflow, velocity, and temperatures in each model segment. Initial estimates for water-surface elevation, slope, and water temperature at each segment were specified. Boundary condition time series at the upstream model segments of the Fox River and Illinois River included streamflow and water temperature. Time series for simulating system wide analyses included air temperature, solar radiation, wind speed, wind direction, and cloud cover. Additional specified model parameters used in determining heat exchange and hydraulic characteristics included longitudinal and vertical eddy viscosity, longitudinal eddy diffusivity, and bottom friction. Recommended or default values were used for most starting parameters. Any parameters modified for this specific application are discussed in the “Hydrodynamic Model Calibration” section of this report.

Streamflows through the Marseilles and Starved Rock Dam structures were simulated in the model by specifying the dynamic weir formula option of gate simulation at the corresponding model segment locations. The dynamic weir formula used in the model is:

where where Although equation 1 calculations were done within the CE-QUAL-W2 model for both dynamic weir locations, the equation 2 calculations were done outside of the model and the resulting input time series were used in the dynamic weir equations for Marseilles and Starved Rock in subsequent model runs.

The complex lock operation at the Marseilles Lock and Dam is represented by a streamflow time series in the CE-QUAL-W2 model. The gate structure in the Marseilles Lock is depicted in the model at the downstream face of segment 30 (fig. 6, inset 2). Gate operations in the model are simplified by using a discretized daily mean computed flow time series to regulate flow between segment 30 and the downstream segment. The daily mean lock flow is calculated by the USGS using a USACE-supplied number of lock operations per day and a developed rating table between number of lockages and streamflow to compute a time series of daily mean lockage flows. In this way, the gate is continuously open in the model, but the flow is highly regulated as a result of lock operations and the representative flow time series. The resulting computed water-surface elevation at the segment downstream from the lock gate was compared to the observed water-surface time series at the USACE stage gage at this location (fig. 6, inset 2) during model calibration.

Water temperature is an important model variable and is included in all CE-QUAL-W2 hydrodynamic model applications because of its effects on water density and hydraulic characteristics. The hydrodynamic role of water temperature, while important, is secondary in this study to the objective of assessing the specific role of water temperature as a physicochemical characteristic in the development of HABs.

Hydrodynamic Model Calibration

Model calibration involves the adjustment of model parameters until simulated and observed data are within an acceptable threshold. Hydrodynamic model calibration involved adjusting model parameters until observed and simulated daily water-surface elevations, daily streamflows, and the observed and simulated streamflows, velocities, and channel areas associated with discrete streamflow measurements were similar. The period of January 1, 2021, through December 31, 2022, was used for model calibration, and the model output parameters were validated using an independent dataset for the period of January 1, 2020, through December 31, 2020. The 2021–22 period was selected for calibration because of additional temperature profiles and algal community data collected in the reach during this period. Calibration of time-series data including streamflow, water-surface elevation, and temperature was assessed using three goodness of fit metrics, including the mean absolute error (MAE), the root mean square error (RMSE), and the Nash-Sutcliffe model efficiency (NSE). The MAE is the arithmetic average of the absolute error between simulated and observed values (Moriasi and others, 2015), has a range of 0 to infinity (∞), and the units are the same as the time-series data. The MAE was determined by the following equation:

where The RMSE also has a range of 0 to ∞ and has the same units as the time series evaluated. RMSE was determined using the following equation: which is the statistical measure of the standard deviation of the residuals between observed and simulated values. Lastly, NSE (Nash and Sutcliffe, 1970), which has a range of −∞ to 1 was calculated as follows: Discrete observed streamflow, velocity, and channel area data pairs from instantaneous streamflow measurements were compared to the nearest corresponding model segment face value using the absolute percentage difference calculated as follows:

At the Marseilles Lock and Dam, the model calibration process involved matching simulated and observed dam pool water-surface elevations, simulated and observed dam outflows, and simulated and observed lock tailwater water-surface elevations. The calibration steps of the Marseilles Dam pool water-surface elevation time series and dam outflows were:

Calibrating simulated water-surface elevations to observed water-surface elevations at the USACE tailwater streamgage at the Marseilles Lock involved the modification of bed friction and branch slope to match simulated water-surface elevation values to observed values.

The 2021 through 2022 calibration results of simulated water-surface elevation and streamflow time series at the Marseilles Lock and Dam indicated that the simulated and observed values were similar (fig. 7A, C). The MAE and RMSE values from the calibration of the Marseilles pool water-surface elevations were 0.03 foot (ft) and 0.04 ft, respectively, indicating close agreement (table 3). Similarly, the simulated and observed water-surface elevations at the Marseilles Lock tailwater were similar (fig. 8B), with an MAE of 0.23 ft and an RMSE of 0.29 ft (table 3). The NSE value of 0.99 for the calibration of the simulated Marseilles Dam outflows with the observed streamflows at the USGS Marseilles streamgage also indicated close agreement. The 2020 validation period was dominated by the May 19, 2020, event for which water-surface elevations and streamflows were underpredicted by the model at the Marseilles Lock and Dam (fig. 7A, B), but this can be attributed to the Q_max value in the dynamic weir equation (eq. 1) being maintained at the maximum flow value during the calibration period. The MAE for Marseilles Lock and Dam pool and lock tailwater water-surface elevations during the validation period were 0.04 and 0.32, respectively, and similar to corresponding values during calibration, whereas the NSE for the Marseilles Dam outflows during the model validation period was 0.99 (table 3).

Calibration of pool water-surface elevations and total releases at the Starved Rock Lock and Dam only involved comparing simulated and observed combined dam and lock outflows and simulated and observed pool water-surface elevations. Calibration of the simulated and observed lock tailwater elevations was not done because the downstream model boundary is the Starved Rock Dam pool. The simulated Starved Rock Dam pool water-surface elevations were compared to observed elevation time series at the USACE streamgage at Starved Rock. The simulated total outflows were compared to the sum of the daily flows for Fox River at Dayton and Illinois River at Marseilles. The calibration steps for the Starved Rock total outflows and pool water-surface elevation were the same as those described previously for the Marseilles Lock and Dam except the $Q_{si{m}_t}$parameter in equation 1 for the Starved Rock application referred to total simulated dam and lock flows rather than just the dam outflows.

The calibration results of simulated water-surface elevation and streamflow time series at the Starved Rock Lock and Dam also indicated that the observed and simulated values were similar (fig. 7B, D). The calibrated Starved Rock pool water-surface elevations resulted in an MAE of 0.06 ft and RMSE of 0.08 ft and the NSE value of calibrated total outflows was 0.99 (table 3). Comparison of observed and simulated values during the 2020 validation period also indicated close agreement in the time series because the MAE and RMSE values of pool water-surface elevations were 0.07 ft and 0.11 ft, respectively, and the NSE metric of observed and simulated total streamflows was 0.99 (table 3).

The calibration of the water-surface elevations at remaining nonlock and dam locations including the USGS streamgage at Marseilles, the USACE streamgage Illinois River at Ottawa, and the USGS streamgage on the Fox River at Dayton (fig. 1) involved modifying the channel friction parameter (Manning’s n) and branch slopes to match the simulated to observed time series. The calibration and validation results at these locations also indicated close agreement between observed and simulated values (fig  8). Model performance metrics for the calibration of water-surface elevation indicated MAE ranged from 0.20 ft at the Fox River at Dayton to 0.35 ft at the USGS Marseilles streamgage (table 3). The RMSE values ranged from 0.29 ft at the Fox River at Dayton and Illinois River at Ottawa to 0.46 ft at the USGS Marseilles streamgage. Performance metrics during the validation period at the three sites were similar to the metrics measured during the calibration period, with a maximum mean difference in water-surface elevations of 0.58 ft (RMSE) at the USGS streamgage at Marseilles (table 3).

Calibration of Instantaneous Streamflow Measurement Data

Streamflow measurement data collected near the USGS streamgages at Marseilles and Dayton provided the opportunity to compare observed point-location streamflow, velocities, and channel areas at the measurement locations to simulated values. Seventeen streamflow measurements were made near each of the Marseilles and Dayton streamgages during the January 1, 2020, through December 31, 2022, calibration and validation periods. Corresponding simulated streamflow, velocities, and model channel area at the nearest segment face at the time of the streamflow measurements were used for comparison.

Sensitivity Analyses of Temperature Simulations

To determine the relative effects of changes in input time series and selected model parameters on simulated water temperatures, eight additional model scenarios were run and the resulting changes in water temperature were quantified. The sensitivity analysis runs were done for the period of January 1, 2021, through December 31, 2021, and results are presented for the Starved Rock Lock and Dam pool at the end of the study reach. The base model input time series was changed as follows:

Results indicated that the model was most sensitive (greatest quantitative effect on simulated water temperature) to the defined increases and decreases in air temperature followed by increases and decreases in streamflow, parameter AFW, and parameter WSC (table 8). Although the monthly change in air temperature in nonsummer months may exceed 15 percent, the temperature changes in the summer months were about 15 percent. This 15-percent change was then used as the basis for the comparable 15-percent change in the other time series and parameters included in the sensitivity analyses. The month with the greatest effects from the simulated changes in projected air temperature, wind sheltering, and surface heat exchange parameters was September, with a simulated increase of 2.65 °C under future RCP 4.5 2075 conditions, a decline of 0.30 °C with the 15-percent increase in parameter AFW, and a decline of 0.12 °C with the 12-percent increase in WSC (the model became unstable with a 15-percent target increase in WSC and a 12-percent increase was used). The months with the greatest effects from a 15-percent change in streamflows were June and December. A 15-percent increase in streamflow resulted in a decline in water temperature of 0.1 °C in June and an increase in water temperature of 0.1 °C in December. A 15-percent reduction in streamflow resulted in a 0.16 °C increase in water temperature in June and a 0.1 °C decline in December. The months that responded the least to the change in air temperature, streamflows, and the parameters WSC and AFW were February and March with resulting changes in water temperature of 0.01 to 0.88 °C, respectively.

Streamflow Characteristics

Notable streamflow characteristics leading up to and associated with the occurrence of HABs could be substantial spring runoff events that can transport nutrients into the system from agriculture and urban sources followed by a dry period and declining flows that allow nutrients to be retained in the system. In addition to the role of nutrient delivery and transport, streamflows also are a key component of heat transfer and a conveyer of phytoplankton communities. Of the 3 years of analysis, the streamflow characteristics of May to July of 2021 and May to July of 2022 were the most similar and contrasted with the higher flow conditions of May to July of 2020. Agricultural field preparation and fertilizer application in the study area generally occurs in spring (April–May generally are peak months; Illinois Department of Agriculture, 2023; fig. 11) and fall (October–December). Fertilizer application quantity and timing were similar among study years (fig. 11), but the quantity of post application precipitation was less in June 2022 than in 2020 and 2021. Although the timing and rates of fertilizer application may vary along with timing and magnitude of precipitation and resulting runoff (fig. 12), there also were consistent sources of nutrients into the study reach from upstream urban areas (Peake and Hodson, 2022). Streamflows in late June for all 3 analysis years were relatively stable and less than 10,000 ft³/s on the Illinois River in the Marseilles reach and less than 2,000 ft³/s on the lower Fox River (fig. 12).

Water Velocities

Simulated flow-weighted velocities were compared at segments representing selected geomorphic channel features of the Illinois River study reach at different longitudinal locations along the reach; the channel features were the main stem, side channel (segments fed by main stem that feed back into a main stem segment), and backwater areas (segments fed by main stem or side-channel segments but “dead end” and do not feed back into the main stem or side channel). The results of simulated water velocities were compared from each of five longitudinal groupings of segments representing the three geomorphic channel features (figs. 13–15). As expected, simulated velocities decreased substantially between the main channel, side channel, and backwater features, but there were no substantial differences between velocities in common areas between years. Mean velocities in the main channel segments in each year generally were between 1 and 2 feet per second (ft/s). Side-channel velocities generally were between 0 and 1 ft/s, whereas velocities in the backwater areas generally were within +0.1 and −0.1 ft/s and never exceeded 1 ft/s.

Residence Times

The model-computed residence times during the June 2020, June 2021, and June 2022 analysis periods at the selected comparison segments reflect the differences in velocity and segment water volume exchange characteristics (figs. 13–15). Residence times were normalized to a 500-m common segment length. Two gradients in residence times were apparent from the simulation results. There was a substantial lateral difference in residence time in segments between channel features as values ranged from 0.28 to 15.9 days per 500 m of channel. Residence times also generally increased longitudinally from upstream to downstream within common channel features; for example, residence times within the selected main stem segments increased 370 percent in 2021 (fig. 14) from upstream from the Marseilles Lock and Dam to the upper Starved Rock Lock and Dam pool. These gradients in residence time reflect the effects of differences in velocities, segment volumes, water temperature, and the operation of the lock and dams. The relatively high residence times of segment 26 (refer to fig. 6 for locations) among side-channel segments, for example, are a result of its location within the Marseilles Lock. Variation in average June residence times also was greater among channel features than among years in the study period (fig. 16). Average simulated residence times in June 2022 were 5–16 percent higher than in June 2020 and June 2021—periods of extensive HABs—depending on the channel feature and year of comparison. The results quantify the substantial differences and complexities of water movement among channel features in the study reach.

Water Temperatures

Simulated daily water temperature profiles also were compared at the five selected longitudinal locations along the Illinois River for June 2020, June 2021, and June 2022, and variability among years was greatest for this physical characteristic (figs. 17–19). For all years, there was a sharp contrast in water temperatures between the beginning of June and the mid- to late-June periods with increases in water temperature from about 20 to 28 °C or higher. This contrast was most evident in June 2021 in which water temperatures of near 20 °C transitioned into an extended period of water temperatures of greater than 26 °C in under a week (fig. 18), corresponding to the formation of the HAB that year. The weakest contrast in temperatures was in the June 2022 simulated period (fig. 19). The period of water temperatures greater than 27 °C in June 2020 was relatively short compared to similar periods in 2021 and 2022 (fig. 17), but this relatively short period of warm water temperatures coincided with the HAB that year (Illinois Environmental Protection Agency, 2023). The longest period of temperatures greater than 27 °C in the selected segments in the 3-year analysis period was in the second half of June 2022, yet only a localized HAB was documented that year. There were multiple short cooling periods (near 20 °C) in the first half of June 2022, and shortly after water temperatures approached 26–28 °C; at midmonth another cooling period broke up the possible extended warm water-temperature period. This cooling period varied in duration and intensity depending on channel feature and longitudinal location but was most prevalent in the segments downstream from segments 51, 28, and 120 (fig. 19) and the Marseilles Dam (refer to fig. 6 for locations). This could account, in part, for the lack of an extensive HAB in 2022 as the cooling water temperatures were most pronounced downstream from the Marseilles Lock and Dam pool where a HAB was documented that year. It is possible that the cooling period altered the development of and composition of the algal community and limited the spread of the bloom. Warm water temperatures exceeding 27 °C, however, were again present throughout the main stem and side-channel areas from June 21, 2022, through the remainder of the month.

The timing of the water-temperature changes among the lateral transect segments generally was consistent between geomorphic channel areas and between years, but side channels and particularly backwater areas did exhibit greater variability in timing compared to corresponding main-stem segments. Examples of greater variability include the temperature increases in the side-channel segments 40 and 26 between June 6–11, 2020, and similar temperature increases in backwater segments 140, 207, and 200 from June 5 through June 10, 2020, compared to the corresponding homogenous main channel segments. Simulated temperatures generally were uniform with depth in this relatively shallow system, with exceptions at some side-channel and backwater segments (for example refer to segments 40 and 26 in mid-June 2020, 2021, and 2022 in figs. 17–19). The results also demonstrate that some consistent spatial and longitudinal patterns in temperature magnitude are likely within and among years. There generally is a consistent lateral increase in temperature magnitude from main stem to backwater areas for a specified day; however, the longitudinal gradients vary by geomorphic channel feature. The selected side channels indicate a gradient in water temperature from downstream (lower) to upstream (higher) whereas the selected backwater segments indicate the opposite with an upstream (lower) to downstream (higher) gradient. The lateral and longitudinal temperature gradients can be attributed, in part, to differences in water depths and residence times.

The earlier warm water temperatures in backwater areas compared to the main stem could provide the opportunity for these areas to serve as algal community “nurseries” for the main stem, but the residence times in these backwater areas likely exceed the temperature lag times between the backwater and main stem. The documented HABs in 2020 and 2021 at Starved Rock Dam occurred during temperatures of about 26 °C and backwater areas at some locations did attain these temperatures 2–3 days before the conditions in the main stem. Residence times in the backwater areas, however, generally exceeded 9 days, thus limiting the exchange of water with potentially higher algal biomass into the main channel. Overall, the model results indicated that CE-QUAL-W2 captures the variability in water temperatures along with the highly variable hydrodynamics that provide a complex array of habitat conditions associated with algal community dynamics and the formation of HABs.