**3. A Simple Hydrologic Model Using Simulink**

During the initial pilot test of the course, students also investigated the use of Simulink modeling for hydrologic simulation. As part of this work, one student created a Simulink version of the Soil Conservation Service (SCS) runoff model that makes up part of the TR55 model [35]. The method is based upon the unit hydrograph approach that convolutes a unit hydrograph with an incremental runoff time history. The student found that creating a Simulink version of the SCS runoff model was challenging, as separate models were needed to generate the unit hydrograph and calculate the runoff hydrograph (Figure 5). The level of Simulink knowledge needed to create the model was likely beyond what reasonably could be expected in the introductory modeling course. To accomplish the runoff calculation, the unit hydrograph generation model was first called, executed, and its data passed to the MATLAB workspace (Figure 5a) before the runoff generation model was executed. The runoff model then read the unit hydrograph data from the workspace and generated the incremental runoff based upon precipitation and land cover (Figure 5b). Creation of the runoff hydrograph utilized two specialized data processing blocks (memory, digital filter) from the Simulink library that had not previously been used in the course.

**Figure 5.** Soil Conservation Service (SCS) Method Runoff Calculation Model Using Simulink. Unit hydrograph generation model (**a**) is run first with output passed to workspace. Runoff calculation (**b**) reads unit hydrograph from workspace and convolutes signal with incremental runoff time history.

**4. A Multi-Constituent Water Quality Model Created with Simulink** 

As a further test of Simulink's capabilities, a Simulink model suitable for modeling water quality conditions in lakes and ponds was developed. The model uses the same finite volume formulation described earlier and the kinetic formulations of CE-QUAL-W2 [7]. This Simulink based version of CE-QUAL-W2 is currently able to simulate and output all but five of the twenty-nine possible CE-QUAL-W2 constituents that are available in version 3.6 of the model. Total inorganic carbon, alkalinity, iron, epiphyton, and zooplankton are not currently modeled, and only one phytoplankton state variable can be currently modeled in Simulink. For this initial trial version, the Simulink implementation is also limited to a one-dimensional volume balance based hydrodynamic model with a spillway outlet. Even with these limitations, the model is useful for simulating temporal and vertical eutrophication dynamics in lakes and ponds. Its utility was investigated for both teaching and research purposes.

The overall organization of the water quality model that was created for this purpose (called W2\_SL) relies on a hierarchical arrangement of subsystems that is an essential part of Simulink models. The uppermost system (Figure 6) contains subsystems for input of model kinetic parameters, initial conditions, and physical characteristics of the system. The model has a single horizontal segment, with three vertical layers. A display subsystem allows for visualization of the concentration time histories for each layer of the segment (Figure 6).

**Figure 6.** Uppermost system in the one-segment 3-layer water quality (W2SL3) model executed in Simulink. Each box is a separate subsystem for setting inputs and parameters (left) or solving volume, heat, and mass balances (right). Lines and arrows indicate data transfers between subsystems.

For this application, the water column of the segment is vertically separated into three layers (Figure 7). As implemented in CE-QUAL-W2 [7], layers can be of variable width and the top layer can have a time-varying elevation. A volume balance for the segment is performed to determine the time-varying elevation of the upper layer. Within the three layer subsystems are heat balance and mass balance subsystems. Each layer can exchange heat and mass with adjoining layers. The bottom layer exchanges heat and mass with the sediment layer. The surface layer exchanges heat with the air above (Figure 7).

**Figure 7.** W2\_SL model of segment 1. Separate subsystems (shown with boxes) calculate the heat and mass balances for each layer or the volume balance for all three layers. Lines and arrows indicate data transfers between subsystems.

Subsystems for heat and mass balances are found in each of the vertical layers within a segment. The set of constituents follows that of CE-QUAL-W2 [7], although at present only a single algal group is considered, and some other state variables are not included as described earlier. The mass balance subsystem (Figure 8) is further divided into subsystems that contain organic matter mass balances (labile and refractory dissolved and particulate matter plus algal organic matter), nutrient mass balances (nitrogen, phosphorus, silica), dissolved oxygen, and assorted other constituents (salinity, total dissolved solids, residence time, fecal coliform). Input data signals (e.g., meteorological data, water temperature) are used to calculate temperature-varying organic matter process rates (Figure 8).

The labile particulate organic matter (LPOM) constituent is used as an example as to how conceptual process descriptions are translated into a mass balance subsystem in Simulink. The process description block diagram (Figure 9) as specified in the CE-QUAL-W2 user's manual includes inputs to the LPOM from mortality of algal organic matter, ephiphytes, and macrophytes. Sinks of LPOM include sinking and decay that produce either RPOM or inorganic C, N, and P. Production of inorganic C, N, and P also consumes oxygen (Figure 9).

**Figure 8.** W2\_SL model of the mass balances for the surface layer of segment 1 as implemented in Simulink. Separate subsystems (shown with boxes) calculate the organic matter, nutrient, dissolved oxygen and other mass balances for each the layer. Lines and arrows indicate data transfers between subsystems.

**Figure 9.** CE-QUAL-W2 mass balance processes internal to a volume for labile particulate organic matter (LPOM).

The Simulink model's LPOM subsystem (Figure 10) includes additional transport processes not present in the box and arrow diagram of the CE-QUAL-W2 manual (Figure 9, [7]). As with the tracer mass balance shown earlier (Figure 4), a continuous time integration block takes input from a summation block. The summation block inputs represent each mass input or mass output term in the mass balance. For LPOM in the surface layer there are seven terms in the summation block, three of which are sources (AOM to LPOM via mortality, loading, diffusion from the layer below), while four are mass sinks (settling to the lower layer, LPOM decay, LPOM to RPOM, LPOM outflow). Labeling of input and output signals was found to make it relatively easy to use the block and arrow diagram as a process descriptor. The overall visual layout of the subsystem is similar to the conservative tracer subsystem shown earlier (Figure 4), which is not surprising since this subsystem, and in fact all the mass balance subsystems were created by cutting and pasting of other subsystems, using the tracer mass balance as the starting point.

**Figure 10.** Labile particulate organic matter (LPOM) constituent found with the nutrient mass balance subsystem within each segment and layer as executed in Simulink.

As part of model verification, several model tests were performed to confirm that the model was correctly implementing the mass balances and water quality kinetics described by CE-QUAL-W2. One such test was a comparison of corresponding predictions for the two models for a case with a steady input (both inflow and concentrations held steady) of water and nutrients (N, P, Si) into the surface layer of a pond with a spillway outlet. Identical forcings (inflows, inflow concentrations, meteorological forcings) were specified in the two model applications. All kinetic parameters were set to the same set of default values for CE-QUAL-W2. There were some unavoidable differences in the physical setup, as the minimum number of horizontal segments for CE-QUAL-W2 is two, while the Simulink version had only one horizontal segment. Total pond length and width were the same, and each model application had three layers of identical thickness. The pond temperature was set via setting inflow temperatures and meteorological forcings so that the temperature did not limit algal growth. Initial conditions for nutrients were specified such that there was a brief phytoplankton bloom that occurred when growth was not limited by either light or nutrient conditions. Light conditions were set such that growth limitation due to self-shading occurred only early in the simulation once algal abundances increased to sufficient levels but before nutrient depletion became the growth limiting factor. Eventually the nutrient limited algal growth rate (PO4 was the limiting nutrient) exactly balanced biomass losses through flushing, sinking, mortality, and respiration. Recycling of algal organic matter serves as a source of particulate and dissolved organic matter. As expected, an equilibrium condition was finally established where each water quality constituent asymptotically approached a constant concentration.

**Figure 11.** Comparison of model predictions in the surface layer for orthophosphate (PO4, top panel), labile particulate organic matter (LPOM, middle panel), and algal organic matter (AOM, bottom panel) for a test case of a pond modeled as one (Simulink) or two (CE-QUAL-W2) horizontal segments. For the CE-QUAL-W2 case, the most downstream of the two horizontal segments is shown.

Concentration time histories for three representative constituents (orthophosphate—PO4, labile particulate organic matter—LPOM, and algal organic matter—AOM) as predicted by CE-QUAL-W2 and W2\_SL (Figure 11), show nearly identical values for the asymptotic concentrations, but some transient differences in the concentrations in the first few days of the simulations. The final concentrations for the three constituents varied between 0.0 percent (PO4) and 0.87 percent (AOM). The PO**4** concentration declined more rapidly in the Simulink case, reaching 10 percent of the initial concentration in 10.5 rather than 13.1 days (Table 1). Smaller relative differences were observed between the corresponding LPOM and AOM peak concentrations (1.63 and 0.88 percent) and the times to peak AOM concentration (11.1 percent, Table 1). These differences are thought to be due to the different physical configuration of the two systems and not to model errors or limitations of the modeling approach. We are currently developing a Simulink model having two horizontal segments that will allow for a better comparison test between the two models.

**Table 1.** Comparison of representative statistics for predictions of orthophosphate (PO4), labile particulate organic matter (LPOM), and algal organic matter (AOM, bottom panel) for a test case of a pond modeled as one (Simulink) or two (CE-QUAL-W2) horizontal segments.


For teaching purposes, a working but not completely functional version of the model was given to students who were working cooperatively as part of a final project for the water quality modeling course. The project task assigned was to add the necessary functionality and then use the model to predict dissolved oxygen concentrations for a system receiving pollutant inputs of organic matter and ammonia. Small groups of students were assigned responsibility for creating submodels of individual constituents such as ammonia, dissolved oxygen, or BOD. Other students were responsible for connecting the submodels together and generating model solutions for the group.

The students successfully used the single segment three-layer water quality model based upon CE-QUAL-W2 as the basis for a final project in the water quality modeling course. During the project students were able to work together to combine their individual contributions into a single model. It was necessary to provide at the start of the project a working, but not fully functional version of the model to the students. In particular, several constituents (e.g., ammonia, nitrate, dissolved oxygen) and several process calculations (e.g., dissolved oxygen reaeration, sediment oxygen demand) were left out of the model that served as the project's starting point. Despite the fact that students differed quite significantly in their numerical problem solving ability and in their experience with high-level programming languages, they each were able to contribute to the group's efforts, and collectively they did complete the model development part of the project. Overall both the students and the instructor found the model development project to be successful.

#### **5. Discussion and Conclusions**

Our investigation of Simulink has shown that it is possible to create multi-dimensional, multi-constituent finite-volume based water quality models without the use of a high-level programming language. A version of the model with the complete set of CE-QUAL-W2 kinetics is being tested for research use. Currently the model is being tested with the data sets included with CE-QUAL-W2. Simulation times for a 400-day run are approximately 5 minutes. Significantly faster run times are expected once a fully compiled version of the model is created. Once the model has been thoroughly tested, we plan to use it as part of a study investigating toxic byproducts of algal production that have been observed in agricultural ponds in the Limousin region of France. This novel Simulink-based multi-dimensional water quality model seems like a particularly good choice for this application because additional constituents may be needed to simulate the processes that lead to toxin accumulation in the pond. The Simulink model environment has been found to be quite convenient for adding and deleting of constituents in the model depending upon project need.

The model as currently written could be applied for systems with a reasonably small number (1–5) of segments and layers. Additional development is focused on creating a model that can conveniently be applied to systems with many segments and many layers, as these models are commonly used. The model as shown in this article requires manual construction of data links between any new segments and layers, which makes expansion to many segments and layers difficult. In testing now is a procedure that uses the "for each subsystem" Simulink construct that allows for specification of a single set of subsystems (e.g., mass, heat, volume balances) that can simultaneously be applied to all segments and layers of the physical system. Use of this new construct should make the model much easier to expand to multiple horizontal segments and vertical layers.

Our experience with Simulink in the classroom is also promising. Students involved in the pilot test class and the first-time use in the regular course, report that problem solutions created with Simulink take approximately the same time as those earlier made with Excel and are much less tedious to debug and grade. Students report that the copy and paste feature for subsystems and objects makes reuse from problem to problem relatively simple. Overall they had a positive experience in using Simulink as a model development environment in the water quality model course. Based upon these experiences, Simulink will be used for model development the next time the course is taught.

#### **Acknowledgments**

This material is based on work supported by the National Science Foundation/National Institutes of Health joint program in ecology and evolution of infectious diseases under grant no. OCE-0813147.

#### **Conflicts of Interest**

The authors declare no conflict of interest.
