A decision support system for indirect potable reuse based on integrated modeling and futurecasting

Reliable and optimal operation of a water reclamation facility (WRF) is a fundamental risk management component for an indirect potable reuse (IPR) system. It requires timely and informed decision-making in response to fluctuating operational conditions, such as weather patterns, plant performance, and water demand. Futurecasting of IPR systems consists of modeling different near-term future scenarios supplemented by medium- to long-range weather forecasts, which provides useful feedback to decision-makers. Such predictive analyses can aid in alleviating future risks associated with water availability and quality without the need to install and maintain large standby capacities.

Over the years, the interest in integrated water resource management (IWRM) has increased significantly. IWRM takes a comprehensive approach to water management by viewing the water supply, drainage, and sanitation systems holistically. To develop design and management tools for a better understanding of these integrated systems and their complex behaviors, an integrated modeling approach is required. According to Rauch et al. (2002), the concept of integrated modeling was proposed as far back as 1970. The first integrated urban drainage model was applied in 1980, though the concept did not become widely adopted until the 1990s (Mitchell et al. 2007). Since then, there have been many applications of IWRM modeling focused on the simulation of entire urban water systems (Schütze et al. 1996; Coombes & Kuczera 2002). Similar approaches have also been used in integrated watershed management. Wu et al. (2006) implemented a watershed runoff model HSPF (Hydrological Simulation Program – Fortran) coupled with a reservoir model CE-QUAL-W2 for the Swift Creek watershed in Virginia. A more complex model of seven HSPF and two CE-QUAL-W2 applications has been used for the Occoquan Watershed in Northern Virginia (Xu et al. 2007). Integration of different models to simulate an IPR system, however, has not been reported in the literature to the best of the author's knowledge.

There are several approaches to integrated modeling. Some software applications provide built-in integration to simulate a coupled system, where the models for different sub-systems are either contained within the software or are linked automatically. Coupled models are based on integrating different standalone models that simulate various components of the target system individually (e.g. watershed, reservoir and treatment plants). Integration of such models is carried out using different levels of sophistication: from manual coupling to fully automated systems. Coupled models are commonly operated sequentially (Schütze et al. 1996; Rauch et al. 2002; Xu et al. 2007). This loosely coupled modeling approach is unidirectional and flow paths are configured as a tree-like structure. Although this approach seems relatively simple, there is voluminous data input, output, and transfer between the models (Azmi & Heidarzadeh 2013). In addition, some of these models are not very user-friendly, take a significant amount of time to run, and require a significant level of familiarity with the model structure to operate. A more advanced coupling approach called iterative coupling, on the other hand, relies on the tight coupling of sub-systems to create model synchronization based on back-and-forth data transfer for each timestep using either a standard or custom protocol. An example of one such standard framework is OpenMI (Open Modeling Interface) (Moore & Tindall 2005).

Futurecasting involves strategic planning by evaluating system dynamics, predictive analysis, and a variety of other strategies to help develop an insightful vision of the future. Shah et al. (2017) examined short- to medium-range forecast data (7–45 days) to study water and agriculture resource management in India using the GEFS (Global Ensemble Forecast System) and CFS v2 (Climate Forecast System). By using forecast data as an input to simulate forecasted runoff and soil moisture, they were able to show that their methodology could provide timely information in decision-making for farmers and water managers. A near real-time drought monitor was also developed to estimate the severity and extent of agricultural and hydrologic droughts (Shah & Mishra 2015).

The overarching goal of this research is to develop a futurecasting application based on the integration of various in-house and weather forecast models and historical data sources. The outputs can provide useful feedback for plant operators to identify and analyze strategies to manage the WRF performance dynamically in response to future reservoir conditions. This paper will discuss the development and implementation of the integrated modeling application and will demonstrate its effectiveness as a decision support system (DSS) to inform and improve the reliable operation of the IPR system.

The Occoquan Reservoir is located in Northern Virginia and is one of the two major water supply sources (the other being the Potomac River) for several municipalities in the region. The reservoir spreads over an area of 6.9 km² and its drainage basin, the Occoquan Watershed (Figure 1), and spans across 1,484 km² (573 square miles). The reservoir has a full pool volume of 3.1 × 10⁷ m³, an average water depth of 5.1 m, and a maximum water depth of approximately 19 m close to the dam. The reservoir has a full pool elevation of around 37 meters above mean sea level (MAMSL), and mean hydraulic residence time is approximately 20 days (Xu et al. 2007). The mainstream of the Occoquan Reservoir is formed by the confluence of the Occoquan River and Bull Run tributaries, 14 km upstream of the Occoquan dam. The headwaters of the reservoir extend well up into these two arms of the reservoir which are called the Occoquan River Arm and the Bull Run Arm, respectively. The Upper Occoquan Service Authority (UOSA) WRF discharges into the free-flowing stream above the Bull Run Arm, upstream of stream station ST45. The WRF has a total design capacity of 204,412 m³/day (54 mgd) and currently operates at an average of 121,133 m³/day (32 mgd). The contribution of the WRF to the reservoir's inflow routinely varies throughout the year, usually going up to 50% of the total inflow during the late summers and early fall and, in some extreme cases, comprises up to 80% of the reservoir inflow. Fairfax Water's Griffith Water Treatment Plant (WTP) withdraws water from the tail end of the reservoir close to the dam, resulting in IPR of the reclaimed water.

The Occoquan Watershed Monitoring Laboratory (OWML) is responsible for monitoring the watershed water quality. Over the years, OWML has established monitoring stations throughout the watershed to measure the flow and water quality in different streams (currently operational: ST00/01, ST10, ST25, ST30, ST45, ST50, ST60, and ST70). Continuous flow measurements are carried out at all the stations, and data are logged and transmitted to the OWML's database servers on an hourly basis. To monitor the water quality, grab samples are taken from each station at different frequencies ranging from once a month to four times a month. During storm events, flow-weighted composite samples are collected using autosamplers to evaluate the average concentration of different constituents. Grab samples are also collected from seven locations (RE02, RE05, RE10, RE15, RE20, RE30, and RE35) on the Occoquan Reservoir for quality monitoring (Figure A1, available with the online version of this paper). Temperature, dissolved oxygen (DO), nitrate, pH, oxidation–reduction potential, and conductivity data are collected in situ at different water depths (usually at 1.5 m intervals, starting at 0.3 m from the surface) at frequencies ranging from once a month to three times a month. Many other constituents (such as nitrogen and phosphorus forms and total organic carbon) are measured, via sample retrieval and transport to the laboratory for analysis, at the top and bottom of the reservoir stations. OWML also operates a weather station located close to the WRF (called OWML weather station at UOSA), which measures different meteorological parameters including rainfall, air temperature, solar radiation, wind speed and direction, and humidity, which are transmitted on an hourly basis to OWML's database servers.

Regulation of nitrate in reclaimed water throughout the year is a unique operational challenge at UOSA. Its operating permit limits the total annual nitrogen load to 6.0 × 10⁵ kg (1.316 × 10⁶ lb). More than 90% of this load is discharged in the form of oxidized nitrogen, mainly as nitrate (NO₃). Although the MCL (maximum contaminant level) for nitrate in drinking water is 10 mg/L as nitrogen, the Occoquan Policy (VSWCB 1971) specifies a Water Quality Objective (WQO) of 5 mg/L at the dam. The WRF is, therefore, required by permit to reduce its nitrate discharge concentration when the nitrate concentrations at the drinking water intakes rise above 5 mg/L.

It has been observed that during thermal stratification and extremely low hypolimnetic oxygen concentrations in the summer, the supply of nitrate from UOSA actually benefits the reservoir. Under these conditions, microbes in the sediment utilize nitrate as an electron acceptor in the absence of oxygen, advancing the denitrification process. As a result, the release of phosphorus (P) and other less preferential electron acceptors, Mn²⁺ and Fe²⁺, into the water column is inhibited. The reservoir overturns at the beginning of fall and remains well-mixed throughout the winter and early spring. The mixing redistributes oxygen across the entire water column and restores the aerobic conditions, thereby diminishing the denitrification capacity of the reservoir.

A nitrate discharge optimization study (Bartlett 2013) concluded that the total nitrogen load delivered into the reservoir by the WRF can be more effectively distributed each month based on the reservoir's denitrification capacity. Based on these recommendations, UOSA changed its operational strategy by reducing the degree of denitrification in summer to ensure a higher concentration of nitrate in the effluent and then transitioning back to increased denitrification just before reservoir turnover, generally in early fall. Using the optimized load distribution provided in that study, average monthly concentrations were calculated using the average monthly flows from the last 5 years (2013–2017) compared to the observed average monthly concentrations during the same time period, as shown in Table 1. The results indicated that the WRF is unable to meet the low nitrate loads recommended during winters mostly due to various operational constraints, including limited organic carbon and low water temperatures which hinders denitrification. During at least one winter drought in the past, UOSA reportedly had to add methanol under emergency conditions to improve the denitrification process when exceedance of the WQO was imminent. On the other hand, in the summer months, there is insufficient total Kjeldahl nitrogen (TKN) in the raw influent to generate desired nitrate loads. Although there is underutilized Occoquan Reservoir denitrification capacity in the summer, several in-house studies have shown that certain storm events at low pool elevations may still briefly increase nitrate concentrations to levels above the WQO due to high nitrate loads from the WRF.

As also noted by Bartlett (2013), these optimized loads cannot be used as a fixed allocation schedule for the WRF due to the dynamic nature of the system. In winter months, the optimized nitrate load distribution is quite conservative, as it is calculated based on the worst-case scenario using extremely low natural streamflows (winter drought). The effect of any change in the discharged load from the WRF on the nitrate concentrations at the dam is based on the hydraulic retention time of the reservoir, which may vary from a few days to months, depending on the inflows. Therefore, many factors including the pool elevation, volume and quality (temperature and background nitrate concentration) of stream inflows, withdrawal by Griffith WTP, and weather forecast can all affect the future conditions of the reservoir and consequently its denitrifying capacity. Hence, selecting the desired monthly nitrate concentration for the effluent is still a trial-and-error process, largely based on the operator's judgment.

The growing water demand caused by urbanization complicates this process further. The safe yield of drinking water from the Occoquan Reservoir, including the reclaimed water, is 3.0 × 10⁶ m³/day (79 mgd). UOSA's contribution to the inflows has been steadily increasing over the years and, based on the current projections, will exceed 50% of the safe yield after 2025. In reference to Table 1, this corresponds to an optimized effluent nitrate concentration of around 3.5 mg/L for winters, which is unachievable based on the existing operational capability. Moreover, WRF plants are always subject to perturbation due to maintenance issues and/or system outages which can temporarily compromise plant performance. UOSA's operation is regulated by the Occoquan Policy which requires increasing the plant's standby capacity and engineered storage to address worst-case scenarios which results in significant financial obligations.

Various models have been used by stakeholders to understand and improve the water quality of the Occoquan Watershed.

UOSA is currently using a wastewater and reuse process simulation software IViewOps (Intelligent View of Operations) as a day-to-day tool for simulating changes to operation and effect of assets out of service for maintenance (Sen et al. 2015). IViewOps is a multi-layer model that analyzes and optimizes the plant on three levels: (1) biochemical process modeling, (2) asset's condition and capacity, and (3) hydraulic configuration and bottlenecks in each possible operating configuration. It reads data directly from SCADA (Supervisory Control and Data Acquisition system), LIMS (Laboratory Information Management System), and CMMS (Computerized Maintenance Management System) to run the most recent operational configuration. Although IViewOps is a very helpful tool for plant operators, a more comprehensive methodology is needed to understand the effect of various nitrate concentrations on the current and future reservoir conditions in order to reevaluate the operational strategy if required. For instance, it would be helpful to know: (1) the target effluent nitrate concentration, (2) the effect of higher WRF contribution during droughts, (3) the length of time available before the WQO is exceeded, if the WRF effluent nitrate concentration increases due to any operational constraint, and (4) the effect of an expected rain event, especially after a long dry period.

The National Oceanic and Atmospheric Administration (NOAA) Climate Forecast System (CFS) is a climate model which provides medium- to long-range numerical weather predictions. The atmospheric component of the model has a horizontal spectral triangular truncation of 126 waves (T126), equivalent to ∼100 km grid resolution. The long-range (∼9 months) hourly forecast consists of four cycles per day at 00:00, 06:00, 12:00 and 18:00 UTC (Saha et al. 2014). Forecasts for a specific location can be extracted from the model output by either simply using the value from the nearest grid point or using some form of interpolation between the surrounding grid points (WMO 2012).

The NOAA National Water Model (NWM) is a hydrological model that forecasts streamflows across the continental United States and has been operational since August 16, 2016. A long-range 16-member ensemble forecast is produced every day going out to 30 days in the future. There are four ensemble members in each cycle, each forced with a downscaled and biased-corrected CFS model. The core component of NWM is the Weather Research and Forecasting Hydrologic (WRF-Hydro) model by the National Center for Atmospheric Research (NCAR). WRF-Hydro is configured to use the Noah-MP land surface model (LSM) to simulate land surface processes in NWM (Niu et al. 2011). At present, NWM provides reforecast data and there is no archive available for old forecasts (only 2 days of forecast data are maintained on the FTP server at any point in time). Therefore, the raw forecast is used without any calibration in this model. It must, however, be noted that the current version of the NWM does not include the operation and management of lakes and reservoirs within the model.

The Occoquan Reservoir is modeled using the U.S. Army Corps of Engineers software CE-QUAL-W2, a two-dimensional longitudinal/vertical hydrodynamic and water quantity model (Cole & Wells 2011). It is an open-source model written in Fortran which simulates the water quality in a waterbody, using longitudinal segments and vertical layers while assuming complete lateral mixing (appropriate for long, narrow, riverine reservoirs such as the Occoquan). The reservoir is modeled using four branches and 69 longitudinal segments (Figure A2, available with the online version of this paper). Branch 1 is the mainstream of the reservoir, starting from the Occoquan River and extending up to the dam. Branch 2 is the Bull Run Arm, discharging from the northern part of the watershed and the WRF to Branch 1. Branch 3 and Branch 4 are Sandy and Hooes Run, respectively, which are relatively small tributaries.

The major boundary conditions for CE-QUAL-W2 are (1) flow and water quality parameters for each branch and the reservoir periphery (i.e. distributed flows), (2) meteorological data, and (3) any water withdrawals (e.g. Griffith WTP). Time series for these parameters are entered for a simulated time period in various fixed-width text files using Julian days. Each branch and distributed flow in the model requires a separate input for flows (m³/s), constituent concentrations (TSS, OP, NH₃-N, NO_X-N, and DO) (mg/L), and water temperatures (°C). Additional precipitation data are required for each branch including rainfall intensity (m/s), constituent concentrations (mg/L), and temperature (°C). A separate text file is used to input the meteorological variables for the entire waterbody, including air temperature (°C), dew point temperature (°C), wind speed (m/s), wind direction (radian), cloud cover (0–10), and radiation (W/m²). The discharge from the WRF is modeled as a point source to Branch 2. Water withdrawal from the Griffith WTP is inputted as Branch 1 outlet flows using a separate file. At the start of each model run, initial water level and constituent concentrations are also required. All the major input parameters required for the CE-QUAL-W2 model run are summarized in Figure 2.

A Hypolimnetic Oxygenation System (HOS) has been operational in the Occoquan Reservoir since 2012. It primarily consists of a line diffuser (Figure A3, available online) that runs from the dam to the deepest section of the reservoir along the original stream channel, almost equally divided between segment 47 and 48 of the model. The system is generally started at the end of spring, gradually increases the oxygen flow rate to its maximum during August and September, throttles back to its minimal flow, and then is stopped sometime during fall.

Since CE-QUAL-W2 cannot simulate complex bubble-plume dynamics, oxygen mass rate (kg/day) must be individually inputted for each layer and mixing must be simulated using a vertical mixing coefficient. This method is an oversimplification of complex bubble-plume dynamics and therefore is usually unable to predict the correct oxygen concentrations and temperatures caused by mixing in the water column. A dedicated line diffuser model was proposed by Singleton (2008) to simulate a linear two-dimensional or planar bubble plume governed by the differential equations for mass, momentum, and heat. This linear plume model (LPM) uses the discrete bubble approach to simulate the oxygen distribution and mixing in the vertical layers of the water column. It requires the boundary conditions, such as temperature and oxygen distribution, along the water column to compute oxygen transfer rates and mixing flows.

An iterative coupling approach is required to apply LPM to the reservoir model. By default, CE-QUAL-W2 is compiled into an executable file which cannot be accessed by any external software. It was, therefore, reprogrammed into a Dynamic Link Library (W2Model.dll). A number of new subroutines were added to the existing code that can be accessed externally to (1) initialize the model, (2) run a single timestep, and (3) read and write different configurational and simulated data. A new software, called W2Coupler, was developed in C#.NET to iteratively solve the plume model and transfer data back and forth between the reservoir model and LPM at each timestep. A new graphical user interface (GUI) was developed in W2Coupler, which replaced the default CE-QUAL-W2 Form (screen showing the run progress) and incorporated additional information including the depth of maximum plume rise (DMPR) and oxygen mass rates and entrainment flows for each layer. The GUI contains options for entering LPM parameters and setting properties of the oxygenated segments. The time series for the oxygen flow rate (N-m/min) of each segment is inputted into W2Coupler using a separate text file, oin_x.npt, where postfix ‘x’ is the segment number.

The coupler first runs the reservoir model for each user-defined timestep and reads the simulated boundary conditions, including the temperature and oxygen depth profiles, for each oxygenated segment in the reservoir. It then solves the LPM's nonlinear differential equations using the fourth-order Runge–Kutta method to calculate the oxygen transfer and entrainment flow for each layer, as well as the DMPR. To simulate mixing, W2Coupler automatically adds dynamic pumps (a feature in CE-QUAL-W2 to model pumps in the reservoir) to the reservoir model equivalent to the number of layers in each of the oxygenated segments. Each pump is set to withdraw water from its associated layer, equal to the entrainment flow, and discharge it at the DMPR layer. The calculated parameters are then inputted back into the reservoir model for the next timestep. This iterative process is carried out throughout the simulation until the end of the run. The source code for W2Coupler and the modification in CE-QUAL-W2 are provided in the Supplement Material of this paper for free distribution and reuse.

A new software application, called URUNME, is developed in C# language for model integration and automation (Lodhi et al. 2018). This application works as a thin layer between the users, models, and data sources, practically hiding the mechanics of the modeling process. It provides a user-friendly GUI based on drag-and-drop functionality and uses simple flowchart-type diagrams for project configuration (saved as a .urm file). A single diagram is called a process, which consists of different functions and/or sub-processes interconnected through links. Each function is used to accomplish different tasks including running external models, reading and writing data from different formats (database, text files, excel, netCDF, and WDM), uploading and downloading from web servers, statistics (Built-in R support), time series transformation (aggregation and missing data infilling), querying database (SQL), and expression evaluator (Excel-style built-in formulas). URUNME contains an embedded SQLite database for storing data (saved as a .db file). Data generated by functions are saved in the embedded database and are represented as vector and scalar variables, shown as nodes in the project explorer. These variables support various engineering units and conversions, and can be dragged and dropped onto other functions as an input (e.g. writing to a file) or onto other GUI items for data visualization. URUNME employs a powerful graphics engine for creating insightful and information-rich dashboards for data visualization. Multiple views can be used in the application to create different screens for users. Creating a dashboard is as simple as adding a GUI element (chart, pivot table, gauge, map, grid, image, text, and event manager) on the screen and dragging and dropping variables onto these panels for data binding. These GUI elements are automatically updated in case of any change to the underlying data (e.g. after a model run). One of the key features of URUNME is the capability of creating and running multiple scenarios of the model using different inputs. Any function parameter (e.g. file name) can be made local (specific to that scenario) to run a different condition (e.g. input different meteorological data). Furthermore, different scenarios and processes can be run in a batch mode manually or by using the Scheduler feature, which allows runs to be scheduled at different times or on a periodic basis.

The DSS was set up using two separate projects in URUNME. The first project (model.urm) is responsible for running the integrated model and stores all the simulated and observed data in its embedded database (model.db). The database file is automatically uploaded to a dedicated FTP server after every model run. The second project (dashboard.urm) is created as an interactive dashboard, containing various menus and visual elements (charts, tables, and diagrams) to show the most up-to-date futurecasting results and historical trends (Appendix B, available online). The dashboard project updates its elements directly from the model.db file. URUNME, installed on the computer of different users, uses the dashboard project to download the most recent data from the server to update the outputs on a periodic basis. This structure allows multiple users to run the DSS dashboard application independent of the model implementation.

URUNME reads from various data sources dynamically (data being continuously updated in a different database) and performs extensive data analysis and manipulation before feeding it to the model as inputs or storing it as historical trends. All the operations carried out during the entire simulation process are fully automated and use hundreds of functions in various process flow diagrams. The overall data flow path to operate the model is shown in Figure 3.

The calibration scenario, called ‘base’, is used to evaluate the current model calibration using the last 3 years of observed data. Note that for the Occoquan case described in this paper, these choices on data and simulation length are those taken for the Occoquan model application and are not limitations imposed by URUNME. As long as there is sufficient storage space and memory available, URUNME has no documented limits on the size of the problem it is called upon to run. For this Occoquan run, inflows for all the branches and distributed areas were calculated from observed flows at ST10 and ST45. These flows are converted to model branch inputs by scale factors (Table A2, available online) obtained by dividing the drainage area of each branch by the drainage area of the associated stream station. Water quality data for these stations are read separately for base and storm flows from the laboratory database. Base flow concentrations, measured intermittently at these stations, are interpolated to create a time series for each input constituent and water temperature. The constituent values are then replaced by storm flow concentrations for the duration of each storm event. The meteorological data from the UOSA weather station are read from a text file updated hourly. HOS oxygen flow rates and water withdrawal by the WTP are provided by Fairfax Water through email and read into the embedded database. All the data are formatted and written to different model input text files.

The first step after the base run is to validate whether the simulated pool elevations, temperature, and different constituents are consistent with the values observed during the calibration period. In case of unsatisfactory calibration performance for any water quality parameter, recalibration has to be done manually before running the remaining forecast scenarios. It must, however, be noted that after successful initial calibration, frequent changes in model performance are not expected. Both graphical and statistical techniques are employed to keep track of the model's calibration performance. After every model run, time series plots of observed vs. simulated data are created by the DSS to provide an overview of the calibration and to identify the differences in time and magnitude of the data visually. Generally, an agreement between observed and simulated frequencies indicates an adequate model performance (Moriasi et al. 2007). Four different statistical indices are used for evaluating the model calibration including (1) Nash–Sutcliffe efficiency (NSE), (2) percent bias (PBIAS), (3) root mean square error, and (4) coefficient of determination (R²). NSE is a dimensionless normalized statistic that determines the relative magnitude of residuals compared to the measured data variance (Nash & Sutcliff, 1970). It ranges between –∞ and 1, with NSE = 1 as optimal. Moriasi et al. (2007) have recommended some guidelines for acceptable NSE values for water quality calibration where NSE > (0.5, 0.65) is considered ‘satisfactory’, NSE > (0.65, 0.75) as ‘good’, and NSE >0.75 as ‘very good’. Although these recommendations are based on monthly means, the thresholds can still provide a reasonable basis for comparison for daily data. PBIAS is an error index, which measures the average tendency of the simulated data in comparison with the observed values (Gupta et al. 1999). The optimal value of PBIAS is 0. It can be either positive or negative, indicating under- or overestimation bias, respectively. Similar to NSE, Moriasi et al. (2007) has recommended guidelines for acceptable PBIAS for water quality based on a monthly timestep, where PBIAS < ±(70%, 40%) is considered ‘satisfactory’, PBIAS < ±(40%, 25%) as ‘good’, and PBIAS <25% as ‘very good’. RMSE is also an error index which is the average magnitude of calibration error. Moriasi et al. (2007) have provided guidelines only for RMSE normalized by standard deviation but not for standalone RMSE values. R² describes the proportion of variation in observed data explained by the model. It ranges from 0 to 1, with 1 being optimal. It is very commonly used for model evaluation. Similar to the RMSE, no accepted guidelines are available for R².

To run the forecast scenarios, the latest CFS operational forecast GRIB2 files are automatically downloaded from the NOAA website for the four selected grid points closest to the UOSA weather station (Table A1, available online). The data are first interpolated for the UOSA weather station using the Inverse Distance Method (Flitter et al. 2016) and then bias-corrected before storing into the embedded database. To estimate streamflows, the four most recent ensembles are downloaded from the NWM website. There are a total of 120 netCDF files in every ensemble, each containing data for one timestep of 6 h. Forecasted flows are read from these files as a ‘streamflow’ variable (m³/s) using a feature_id key to identify the required stream's reach. Files for each ensemble are combined, and flow rates for ST10 (feature_id: 22339083) and ST45 (feature_id: 22339527) are extracted from each ensemble using a scale factor of 0.01 (as indicated in the metadata of the netCDF files).

A total of 24 forecast scenarios were created using various combinations of boundary conditions including (1) natural streamflows, (2) WRF discharge, (3) WRF effluent nitrate concentrations, and (4) WTP withdrawal. These user-defined time series include two possible natural streamflows (Q₁ and Q₂ [m³/day]), two possible UOSA discharge flow rates (D₁ and D₂ [m³/day]), three possible UOSA effluent nitrate concentrations (N₁, N₂, and N₃ [mg/L Nitrate-N]), and two possible WTP's withdrawal flow rates (W₁ and W₂ [m³/day]). Future stream temperatures and water quality constituents are calculated using historical trends. Monthly averages for the base (C_b) and storm flow (C_s) concentrations are calculated from the last 5 years for ST10 and ST45. In addition, average monthly base flows (Q_b) are calculated for the same period to identify any storm events in the forecasted flows. Time series for stream concentrations are created using C_b, unless forecasted flow for a given day is greater than Q_b, in which case C_s is used.

Two different forecasts are produced by the DSS. A 30-day forecast is based on NWM streamflows going out 30 days in the future. Q₁ and Q₂ are calculated using the average and either of the maximum or minimum flows for each day in the four ensembles for ST10 and ST45. Long-range 90-day forecast, on the other hand, is used as a What-If analysis to simulate any worst-case scenario. For winter months, this forecast is generally used to simulate an extended drought period and its effect on the reservoir. Such droughts have been observed in the past, during both winter and summer, the worst of which lasted for around 9 months in 1977. Such a scenario potentially increases the WRF's contribution to the inflows of the reservoir, which may lead to high nitrate concentrations, especially during winters. The 90-day forecast is simulated after the first 30-day forecast for all the 24 scenarios. The expected streamflows during a drought condition for these 90 days are calculated using low-flow frequency analysis (discussed later).

The overall execution flow path of the model is shown in Figure 4. First, the base scenario is run for the last 3 years. Using the restart feature in CE-QUAL-W2 v4.1, the boundary conditions are saved in a file at the end of the run for the last day of simulation. Data for each of the forecast scenarios are created and inputted to the model by URUNME before it is run, from the present day to 120 days in the future. All 24 scenarios are run sequentially using the present-day boundary conditions from the saved file, and the output data are saved into the embedded database.

The results presented in this section are for RE02, extracted from the DSS output generated on October 10, 2018. RE02 is located at around 0.5 km from the dam and is the closest sampling station to the WTP intakes (Figure 2).

The base scenario was run for the last 3 years starting from January 1, 2015 to October 10, 2018. Observed and simulated time series plots of pool elevation, temperature, dissolved oxygen, and nitrate concentrations are shown in Figure 5. The temperature and concentration plots contain the actual data collected from the reservoir bottom (shown as dots) along with the simulated model output (shown as solid lines) from the bottom-most layer. The calibration indices for all the reservoir stations, as calculated by the DSS are shown in Table 2.

Due to inadequate hypolimnetic mixing, initial runs of the model underpredicted the bottom water temperatures at RE02. The poorer mixing was caused by a lower value of plume rise (DMPR) computed by the LPM than the observed value of plume rise caused by hypolimnetic aeration.

The DMPR is most influenced by the initial bubble radius and entrainment coefficient (α) in addition to the oxygen flow rate (Singleton 2008). The value of α was decreased from the recommended value of 0.11–0.05. An initial bubble radius of 0.75 mm was used (which is the average of the manufacturer specified bubble radius of 0.5–1 mm). This fine-tuning increased the average DMPR and consequently the bottom water temperatures at RE02. The final calibration results for water temperature show NSE values between 0.73 and 0.97 for the bottom layer at all the reservoir stations. The PBIAS for all the reservoir stations remained within ±10%. RMSE values remained between 1.0 and 2.8 °C, while R² ranged from 0.85 to 0.98.

The nitrate calibration has shown mostly ‘satisfactory’ NSE values, with the exception of RE20 and RE35 which are ‘unsatisfactory’. Bartlett (2013) observed similar results due to high nitrate spikes at RE35 during low flows caused by the backflow from Branch 2 (Bull Run Arm) to the Branch 1 (Occoquan River Arm) of the model. This phenomenon is known to occur in the reservoir during summer and is poorly predicted by the model due to its inability to simulate an extended period of continuous bottom-layer anoxia upstream of RE35. This condition also affected the calibration of RE20 where NSE is slightly below 0.5. The PBIAS, however, has ‘good’ to ‘very good’ for all the stations. RMSE values remained within 0.4–1.26, where higher values are for stations close to Branch 2 (i.e. RE20 and RE30) due to higher inflow nitrate-nitrogen concentration from WRF. R² ranged from 0.51 to 0.78 excluding RE35.

To determine the effectiveness of the bias correction in improving the forecasting skill, the method was first trained for 10 years, between 2006 and 2015, and then applied to the raw forecast data for 2016 and 2017 for validation. Figure 6(a) indicates a good correlation between observed and reforecasted air temperatures in the training period (R² = 0.73, NSE = 0.67, and PBIAS = 11.5%). However, the CDF of these datasets plotted together in Figure 6(b) shows a significant bias, where the reforecast air temperatures are underpredicted in the entire range in general and at lower values in particular. After applying the bias correction to the validation period, the R², NSE, and PBIAS between observed and reforecast data improved from 0.72, 0.63, and 15.5% to 0.73, 0.71, and 3.6%, respectively. Other parameters were also tested using the same method (not shown), although the improvements in the forecasting skill were not significant for some parameters.

As explained earlier, different combinations of boundary conditions (Q, N, D, and W) are used to create forecasting scenarios. Most of these scenarios are selected in consultation with the stakeholders and are reviewed periodically based on different factors including time of the year, plant performance, and water demand.

For the 30-day forecast, the natural streamflows, Q₁ and Q₂ at ST10 and ST45, were obtained using the average and minimum of the four ensembles in the NWM forecast, respectively. The plant's effluent nitrate-N concentration was set to 13.6 mg/L as N (N₂), the last 7-day average from the date of model run. Upper and lower values for the effluent nitrate concentrations, i.e. N₁ and N₃, were set to 0.75 and 1.25 of N₂, respectively. The WRF effluent discharge D₁ was taken from the projections for September and October (Table 1). To simulate higher than projected discharge conditions, D₂ was set to 149,524 m³/day (39.5 mgd) which is approximately 15% higher than average D₁ for the year. Due to unavailability of the projected withdrawals from the WTP for 2018, W₁ for the forecast was taken from the average monthly withdrawals in the last 5 years, i.e. 2013–2017 (Table A3, available with the online version of this paper). Finally, W₂ was set to the safe yield of the reservoir, i.e. 299,048 m³/day (79 mgd), to simulate the highest possible withdrawals from the reservoir.

Low-flow frequency analysis (Gustard & Demuth 2009) was carried out to develop the streamflows during the drought conditions simulated in the 90-day forecast. A 10-year return period (T) is used for this analysis, which is considered as an acceptable standard to calculate drought conditions. For instance, a commonly used term ‘7Q10’ in low-flow analysis represents the minimum 7-day flow expected to occur every 10 years. Historical daily flow data from USGS (US Geological Survey) and the OWML databases were obtained for ST10 and ST45 from 1968 and 1973 onwards, respectively. First, an n-day (where n = duration) centered moving average is passed through the daily data to calculate the time series for the duration of 7, 15, 30, 60, 90, and 180 days. Annual minima were then derived for each duration by selecting the lowest flow every year. A CDF was developed using the yearly minimum flows and low flows were then calculated by selecting the exceedance probability (p) of 0.1 (T = 1/p) on the CDF (Figure 7). Low flows for 90Q10 were selected corresponding to the forecast duration of 90 days. Using these flows, UOSA contribution for 90-day forecast would remain close to 62.5% of the total reservoir inflows. For effluent nitrates, N₂ was set to 10 mg/L based on the average effluent nitrate concentration from October to January (Table 1). N₂ and N₃ were set at 7.5 mg/L and 12.5 mg/L, respectively. Discharge and withdrawal scenarios were set similar to 30-day forecast.

The forecasted bottom-layer water temperature profiles are shown in Figure A4 (available online). The predicted values are consistent with the observed water temperatures during the forecasted months and therefore indicate an acceptable model performance. The nitrate profiles were outputted for the middle layer (Figure 8), which is around the same level from where most of the water is withdrawn by the WTP. These results are for the streamflows in Q₁ scenario followed by 90Q10 drought flows (there was no significant difference between the output from Q₁ and Q₂). The model run showed no significant change in the nitrate during the 30-day forecast (shown with green background) and concentrations remained well below 1 mg/L. NWM predicted relatively higher streamflow for the 30 days which reduced the UOSA contribution to the total inflows to the reservoir to around 11% on average. The dilution caused by the higher streamflows kept the nitrate levels below 3 mg/L at ST45 during all the N scenarios. Additionally, the reservoir denitrification capacity remained significantly high during September due to relatively higher water temperatures causing low nitrate concentrations.

For D₁W₁ (Figure 8(a)), the nitrate concentrations at RE02 reached a maximum of 1.7, 2.3, and 2.9 mg/L for N₁, N₂, and N₃, respectively, at the end the forecast. A drawdown of 0.63 m (2.1 ft) in the pool elevation was simulated by the model due to higher WTP withdrawals compared to the total inflows to the reservoir on average. Although D₂W₁ (Figure 8(b)) scenarios were simulated with higher effluent discharge corresponding to 15% additional nitrate-N load to the reservoir, the results did not show any significant changes in nitrate concentrations compared to D₁W₁. This can be attributed to the lower drawdown in this scenario (0.4 m) due to higher WRF contribution. The relatively higher hydraulic residence time (HRT) in this scenario is counteracting the effect of higher nitrate loads. This signifies that in severe drought conditions, the higher WRF effluent flows would actually be beneficial for the reservoir, delaying the increase in nitrate levels at RE02 due to higher HRT. In both D₁W₂ and D₂W₂ scenarios (Figure 8(c) and 8(d), respectively), the nitrates at RE02 increased faster when the withdrawal was increased from W₁ to W₂ (299,048 m³/day). A higher effluent nitrate-N from UOSA (N₃ = 12.5 mg/L Nitrate-N) resulted in 3 mg/L nitrates at RE02 after approximately 79 days. Again, the main reason was higher drawdowns, 1.5 m and 1.25 m for D₁W₂ and D₂W₂, respectively, which essentially lowered the capacity of the reservoir to resist change.

The nitrate levels shown in the above results remained under the WQO limit in all the 12 scenarios. This would indicate that an extended winter drought based on 90Q10 is not going to adversely affect the reservoir. However, it must be noted that these results are based on a single model output where the boundary conditions were highly conducive for such favorable results. Due to higher than average precipitation in summer and early fall of 2018, the reservoir was already at full pool at the start of the 90-day forecast. The starting HRT of the reservoir was calculated to be approximately 120 days based on the withdrawal at the time. In addition, the background concentrations were significantly low (∼0.5 mg/L). Even at the lower denitrifying capacities of the reservoir during the winter months, the higher HRT provided better resiliency against any abrupt increase in nitrate at RE02. The relatively higher simulated nitrate concentrations for scenarios involving W₂ are only because of the reduction in HRT due to higher drawdown caused by more withdrawals.

As elaborated in the above discussion, the initial pool elevation and, hence, the HRT of the reservoir play a crucial role in the outcome of the forecast. To understand the effect of how lower pool elevations can affect the nitrate concentrations in the reservoir, a hypothetical condition was considered where an already existing drought condition in the fall was extended into winter (the simulation was only done for this study and is not a part of regular DSS forecast). All the 12 scenarios were rerun by creating an artificial drawdown of 2 m just before the start of the 90-day forecast, which corresponds to a reduction of 33% in the reservoir volume (Figure 9). It can be seen that in all the scenarios, the higher effluent nitrate concentration of 12.5 mg/L (N₃) will eventually push the nitrates to the WQO limit in the forecasted period. The worst-case scenario remains D1W2 where the nitrate reaches 5 mg/L at RE02 in around 77 days.

The above results presented how the dynamic nature of the system can affect the reservoir based on different boundary conditions. Different meteorological parameters can also have significant effects on the outcome. For instance, colder than normal air temperature in fall or winter can also increase the nitrate concentration at RE02 by reducing the denitrification capacity of the reservoir. On the contrary, any significant rain event can washout the entire reservoir diluting the nitrate concentrations to much lower levels.

Every week, the DSS moves forward in time by predicting the reservoir water quality based on the most recent initial and future boundary conditions. When the forecasted water quality going out 30–60 days is above the 3 mg/L Nitrate-N threshold, the DSS automatically runs an optimization routine on IViewOps to determine means to reduce the nutrient discharge from the UOSA plant by simultaneously adjusting asset availability (changing maintenance schedules to bring unavailable assets such as reactors, clarifiers, and pumps back online) and operational setpoints. Sen et al. (2016) have discussed the nitrate optimization in detail for winter by using different combinations of process configuration (tanks in service), mixed liquor recycle, DO setpoints, wasting, and supplement carbon, while specifying the constraints on the effluent quality.

The DSS can help the plant operators to decide the extent of denitrification, especially during the winter months. Instead of running the plant at full throttle to maximize denitrification, the operators can target a more realistic effluent nitrate concentration based on the existing condition of the reservoir. Less denitrification means potential cost savings due to the following:

• Reduced amount of carbon required for denitrification and therefore more BOD (biochemical oxygen demand) available for the digesters to increase biogas production for heating and electricity generation.

• Less supplemental carbon (methanol and MicroC) required for denitrification.

• Less anoxic zone mixing requirements and reduce nitrate recycle pumping rates.

• Lower anoxic volumes reduce undesired luxury bio-P uptake resulting in less soluble phosphate release and consequently lower struvite formation in digesters. This reduces the potential cost of descaling of the centrifuge system and struvite control chemical feed.

This paper presented the development and implementation of a DSS for the Occoquan IPR system to forecast and regulate the nitrate concentrations in the augmented reservoir. The entire DSS application is powered by newly developed software, called URUNME, to fully automate the operation of the integrated reservoir and process models, forecasting products, and various data sources. URUNME was also used to develop information-rich and user-friendly interactive dashboards to output the updated historical data and forecasting results for the stakeholders.

The reservoir model is operated once every week, running multiple future scenarios based on different combinations of natural stream inflows, plant effluent flows and nitrate concentrations, and water withdrawal by the WTP. The future weather conditions are simulated using the forecasts obtained from the CFS. A 30-day forecast predicted the reservoir conditions based on the NWM streamflow forecast going out 30 days. A 90-day forecast, on the other hand, is used as a What-If analysis to simulate the effects of any possible future winter drought.

The forecasted nitrate-N concentrations under different scenarios can be used by the WRF to identify and analyze operational strategies dynamically in response to the reservoir future conditions. In the event that nitrates going out 30–60 days exceeds 3 mg/L, the DSS can be used to run different optimization routines on IViewOps to determine the means to reduce the nutrient discharge. Even when the target nitrates effluent concentrations are difficult to achieve due to operational constraint, the DSS output provides a fair estimate of the time it would take for the reservoir to reach unacceptable nitrate levels. This provides valuable feedback to the plant managers for contingency planning, e.g. rescheduling the maintenance of certain assets, adding supplemental carbon chemicals, and changing plant configuration, to avoid any future emergency conditions. As the DSS progresses forward in time, the optimization of the WRF is revisited every 7 days with updated outputs due to the dynamic nature of the system. For instance, any significant rain event can partially or completely wash out the entire reservoir, diluting the nitrate concentrations to much lower levels.

The file transfer system makes it easier for the WRF to run the simulations based on reservoir model runs done at OWML. Each week, URUNME is run to generate a database file that is sent to the WRF via FTP. The plant uses the URUNME interface to run this together with IViewOps and populate relevant data in a screen that is a sister screen to SCADA. The operations and maintenance staff can then make decisions about how to optimize the plant and what the future conditions would be with optimization.

Integrated modeling using the URUNME system provides a comprehensive approach to environmental management by evaluating different components of a system in a holistic manner. It can be used as an effective tool for design, planning, and risk assessments, and can facilitate informed decision-making, enhancing the transparency and collaboration between different stakeholders. Moreover, it can be used for simulating, understanding, and managing a combination of wastewater, drinking water, and groundwater interaction, such as the following:

• 1.

  IPR applications in Southern California (Orange County).

• 2.

  Wastewater treatment, discharge to a reservoir, extraction for drinking water supply, and drinking water treatment (IPR, Wichita Falls, Texas).

• 3.

  Raw water supply from a reservoir, water treatment, and sludge disposal to the wastewater treatment plant (WSSC (Washington Suburban Sanitary Commission), Maryland).

Although developing a reservoir model from scratch is a major undertaking that requires a significant amount of data for setting up and calibrating the model. However, there are certain components of this DSS which can be reused for similar applications without any major modification:

• Downloading the data, reading the binary file formats (GRIB2, NetCDF), and calibrating the forecasts from the CFS and NWM.

• Obtaining data from different data sources (stream stations, reservoir stations, and weather stations), manipulating and transforming data, and writing the boundary conditions and input parameters to the reservoir model.

• Open-source HOS model developed for CE-QUAL-W2.

• Different screens created as part of the DSS dashboard.

The use of a DSS driven by URUNME can be used in various systems in addition to water resources. These can include chemical and biological engineering systems and their models, manufacturing and supply chain management systems and their models. URUNME is available free of cost for academic use (www.urunme.com).

The authors acknowledge the support of UOSA and OWML for supporting and funding this research. The authors would also like to honor the memory of the late Dr Thomas J. Grizzard for his unwavering commitment in the establishment and operation of the Occoquan program and the OWML for over 40 years.