Marine Science and Engineering

: A stormwater runoff model based on the Soil Conservation Service (SCS) method and a finite-volume based water quality model have been developed to investigate the use of Simulink for use in teaching and research. Simulink, a MATLAB extension, is a graphically based model development environment for system modeling and simulation. Widely used for mechanical and electrical systems, Simulink has had less use for modeling of hydrologic systems. The watershed model is being considered for use in teaching graduate-level courses in hydrology and/or stormwater modeling. Simulink’s block (data process) and arrow (data transfer) object model, the copy and paste user interface, the large number of existing blocks, and the absence of computer code allows students to become model developers almost immediately. The visual depiction of systems, their component subsystems, and the flow of data through the systems are ideal attributes for hands-on teaching of hydrologic and mass balance processes to today’s computer-savvy visual learners. Model development with Simulink for research purposes is also investigated. A finite volume, multi-layer pond model using the water quality kinetics present in CE-QUAL-W2


Introduction
Mechanistic models of water quality have long been used for predicting the impact of wastewaters on natural water bodies.While the first use of models relied on analytical solutions to problems of point or distributed sources of conventional pollutants (e.g., biochemical oxygen demand (BOD), suspended solids) into rivers, most of the water quality modeling applications these days rely on very large computer codes written with a high-level programming language such as FORTRAN or C++.
Many of the most used water quality models were originally developed in the 1980s and 1990s using FORTRAN as the high-level programming language.FORTRAN is still widely used; for instance, in recent a review of fifteen lake ecosystem models [1], five were written in FORTRAN and another five were written in C++, with the remainder written in other high level languages.These models are similar in concept and execution in that they are capable of modeling the dynamics in multi-dimensional physical systems of water quality properties, nutrients, dissolved oxygen, organic matter, and biota by simultaneously solving conservation equations for heat, momentum, and constituent mass in numerous water volumes.Examples of these sorts of models include the Water Analysis Simulation Program (WASP) [2], the Environmental Fluid Dynamics Code (EFDC) [3] with the water quality formulations described in HEM-3D [4], Qual2E [5], CE-QUAL-ICM [6], and CE-QUAL-W2 [7].These models are very commonly used in regulatory water quality management (e.g., [8,9]) or to investigate the dynamics of water quality indicators as influenced by dynamic anthropogenic pollutant loading (e.g., [10]).
All of the above referenced models are essentially very large, complex computer codes that are maintained, modified, and manipulated by only a very small group of expert model developers.Manipulation of the code for site-specific or research purposes is not practical for the vast majority of users.There are circumstances, however, when model development is warranted, either when new model state variables are needed, when a special set of water quality processes are appropriate, or when prediction accuracy could be improved by assimilating monitoring data rather than modeling it.Water quality model development is also an important part of education in the field.Making these powerful models more flexible in their formulation through user development could greatly improve their overall utility.
MATLAB/Simulink has been adopted in other fields of engineering as a means for modeling and simulating complex systems without the need to write thousands of lines of computer code during model development.Simulink systems are drawn using a graphical user interface as "block and arrow" diagrams.The Simulink block library is very large and includes math and logic functions, signal generation and processing, visualization, plus many specialized toolboxes (e.g., fuzzy logic, DSP, control systems, neural networks).Lines connecting the blocks represent data signals.Signals can be bundled into buses and selectively unbundled as input into "downstream" blocks.Simulink systems are often operated in a time stepping mode to model the dynamic behavior of the system.As an extension of the MATLAB environment, Simulink models can be easily integrated to read and write data from the workspace and to interact with scripts written with its own high-level programming language.
Simulink has been utilized previously to model the dynamics of engineered systems in a variety of disciplines.Examples of its use in non-environmental systems analysis include process modeling in the sugar industry [11] and building systems modeling [12].In the area of environmental engineering, Simulink has been used for simulating wastewater treatment plants, either as individual unit processes [13,14], as benchmark simulations of an entire wastewater treatment plant [15,16], for simulation of stormwater systems [17,18], or for integrated models that include both engineered and natural surface water systems [19,20].Simulink has also been used to model the hydrodynamics and water quality of a few surface water natural systems [21][22][23].
In this article we explore the benefits of Simulink for water quality model development in teaching and research.Of particular interest is how the graphical modeling environment of Simulink can assist new model developers without programming experience as they are introduced to mass-balanced based water quality simulation in a graduate-level course that intentionally includes both instruction and collaborative model development by students.Simulink has been frequently used in engineering education for situations where construction of powerful computational tools for system simulation is an integral part of the pedagogy.For example, Simulink has been used as a teaching aid in courses on digital and analog modulation [24], multiphase electric machines [25], river water quality [26], and wastewater treatment plant operation and control [27].It has been recognized in the engineering education community that use of tools such as Simulink can provide a "learner-centered" environment whereby students are able to create simulation tools that allow personal experimentation of system behavior (e.g., [28]).Courses structured to include such features can positively affect student interactivity, participation, and course satisfaction [29].
As a case study of water quality modeling instruction and model development with Simulink, we present three examples that were developed either for use in the classroom and/or in research.The first model was developed as an alternative to a spreadsheet solution of conservative flushing and transport through a pond and channel.The second application is a Simulink version of a widely used stormwater runoff hydrology model.The third example describes the finite volume implementation in Simulink of the water quality model CE-QUAL-W2.This application represents one of the first uses of Simulink for modeling dynamic, spatially varying water quality in natural systems.

Models of Pollutant Transport for Teaching Water Quality Modeling
Many environmental engineering curricula include a course on mass-balance based surface water quality modeling.These courses build and expand upon the steady-state, one-dimensional dissolved oxygen/BOD modeling that is usually taught at the undergraduate level.The courses rely on process-based descriptions of fate and transport of constituents in surface water systems that are analyzed by solving mass balance equations for one or more of the constituents.The assignments in the courses typically are computationally intensive in keeping with the subject matter.A number of textbooks are available to teach such a course [30][31][32][33][34]; here we examine the use of Simulink in a graduate-level course that has been taught several times using the Chapra text.
An example problem from the Chapra text is used here to illustrate alternative means of model development for quantitative problem solving.In this particular problem (Problem 10.3), conservative dye is being flushed from a completely mixed pond into a channel by a river diversion (Figure 1).The dye exits the pond as an exponentially decaying source of mass to a downstream channel.No mixing occurs in the channel.After a travel period downstream in the channel, the diverted river water containing the dye is conservatively mixed back into the river.The student is asked to calculate the concentration time history at the pond exit, the channel exit, and after mixing into the river.Students have typically solved this problem with a spreadsheet solution created using Microsoft Excel (Figure 2).Spreadsheet equations are used to calculate the three needed time histories; Excel charts are created to visualize the results.While not an onerous problem for the students, use of these spreadsheet solutions in the class can present several problems.Students generally know Excel, but they vary widely in their expertise with it.Complex problems can overwhelm their capabilities and their patience with cut and paste solutions, even when problem specific parameter names are used.Solutions often involve long complex equations that are difficult to impossible to debug.Solutions typically have no common look and feel and are therefore time consuming to grade.Finally, there is often little opportunity to reuse spreadsheet solutions from problem to problem.
The use of Simulink as a model development tool for problem solving in the graduate-level water quality modeling course was first piloted in an independent study version of the course.Two students who had previously taken the modeling course were offered an independent study to learn Simulink and demonstrate its use through problem-based model development.After a successful pilot program, Simulink was used as the primary means for student problem solving when the water quality modeling course was offered the following year.The course included one week of training in Simulink followed by another week in the computer lab solving some simple example problems from the Chapra text.Over the remainder of the semester, approximately two Simulink problems were assigned each week over the remainder of the fourteen-week semester.While not all textbook problems were amenable to Simulink solutions, there were sufficient problems available in all sections of the course that could be solved with a Simulink model.As an example, the corresponding Simulink solution created by one of the students for the problem described earlier is shown in Figure 3.The Simulink solution couples constituent mass-balance based solutions for the pond, channel, and river.The model also includes a separate subsystem for the volume balance in the pond (Figure 3).Simulink "scope" blocks are used to visualize the concentration time histories at the three locations.The system parameters and concentration time histories are bundled into data buses and passed between the subsystems that need the particular data.where C and Cin (g/m 3 ) represent the constituent concentrations in the reactor and its inflow, W (g/day) is a mass load to the reactor, X (g/day) is the reaction term in the reactor, and Qin and Qout (m 3 /day) are the inflow and outflow of the reactor, which has a volume V.The pond constituent mass balance uses a solution to Equation (1) that has as its basis a solution for an earlier problem having only a single completely stirred reactor (Figure 4).As compared with previous years where students used either handwritten analytical solutions or Excel spreadsheets, student solutions created in Simulink were found to be much easier to grade.Similarity between student solutions was aided by including with the problem assignment example figures showing a recommended overall model structure and some example x-y plots.Students liked the ability to use components from previous problems in their solution and in general caught on quickly to creating solutions in Simulink.Reuse of the mass balance solution for later problems was found to be straightforward.The solution served as the basis of several problems later in the course.

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.

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).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).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).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.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.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.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.

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.

Figure 1 .
Figure 1.Diagram of a representative problem (Chapra, Problem 10.3) where mass passes from a pond to a channel and is conservatively mixed into a river.

Figure 4 .
Figure 4. Simulink model of pond mass balance for a constituent.

Figure 5 .
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.

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.

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.

Figure 8 .
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 10 .
Figure 10.Labile particulate organic matter (LPOM) constituent found with the nutrient mass balance subsystem within each segment and layer as executed in Simulink.

Figure 11 .
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.

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.