Skip to content

Pollutant Module: get and set pollutants during simulation

Jump to bottom Edit New page
Abhiram edited this page Aug 20, 2021 · 23 revisions

1.0 Summary

To address the limitations of SWMM's water quality module, SWMM's C source code was modified to allow users to access and modify pollutant concentrations during any routing time step of a simulation through external tools like PySWMM.

2.0 Current Pollutant Process

The current pollutant process is quite complex. When swmm_run() is called, it initializes three main processes:

  1. it parses the relevant pollutant and treatment information from the SWMM input file [Section 2.1];

  2. it initializes all the relevant pollutant information [Section 2.2]; and

  3. it runs pollutant related functions during the simulation [Section 2.3]. Each of these three processes are explained below. We also direct the reader to a detailed flow chart here.

2.1 Input File

The following functions are called to parse the input file to get the treatment equation.

  1. swmm_run() calls swmm_open()

  2. swmm_open() calls project_readInput()

  3. project_readInput() calls input_readData()

  4. input_readData() calls parseLine() which parses the treatment equation

2.2 Initialize

The following functions are called to initialize water quality states.

  1. swmm_run() calls swmm_start()

  2. swmm_start() calls proj_init() and routing_open()

  3. proj_init() initializes the internal states of all objects

    1. Calls link_initState() which initializes link's water quality states Link[j].oldQual[p] = 0.0 and Link[j].newQual[p] = 0.0

    2. Calls node_initState() which initializes node's water quality states Node[j].oldQual[p] = 0.0 and Node[j].newQual[p] = 0.0

  4. routing_open() calls treatmnt_open() and qualrout_init()

    1. treatmnt_open() allocates memory for commuting pollutant removals by treatment

    2. qualrout_init() initializes water quality concentrations in all nodes (Node[i].oldQual[p] = c; Node[i].newQual[p] = c) and links (Link[i].oldQual[p] = c; Link[i].newQual[p] = c)

2.3 Simulation

The following functions are called during a simulation to handlε pollutants and pollutant treatment.

  1. swmm_run() calls swmm_step()

  2. swmm_step() calls execRouting()

  3. execRouting() calls routing_exec()

  4. routing_exec() calls the following functions:

    1. link_setOldQualState() which replaces link's old water quality state values with new ones using the following code: Link[j].oldQual[p]=Link[j].newQual[p]; Link[j].newQual[p] = 0.0

    2. node_setOldQualState() which replaces nodes's old water quality state values with new ones using the following code: Node[j].oldQual[p]=Node[j].newQual[p]; Node[j].newQual[p] = 0.0

    3. qualrout_execute() routes water quality constituents through the network over current time step by calling several functions.

  5. qualrout_execute() calls the following functions:

    1. findLinkQual() which finds new quality in a link at the end of the current time step. Link quality is that of the upstream node when the link is not a conduit or is a dummy link. Otherwise, it calls several mass balance functions to calculate the link concentration which is set to Link[i]/newQual[p] = c2

    2. findNodeQual() which finds new quality in a node with no storage volume. If there is no flow into the node, then concentration is equal to mass inflow divided by node flow. Otherwise, the concentration is zero.

    3. findStorageQual() calls several mass balance functions to calculate the concentration in a storage node which is set to Node[j]/newQual[p] = c2

    4. treatmnt_treat() which updates the pollutant concentrations at a node after treatment.

  6. findLinkQual() and findNodeQual() calls the following mass balance functions:

    1. massbal_addSeepageLoss() which computes the mass balance accounting for seepage loss

    2. getReactedQual() which applies 1st order reaction to a pollutant over a given timestep

    3. getMixedQual() which finds pollutant concentration within a completely mixed reaction

    4. massbal_addToFinalStorage() which adds mass remaining on dry surface to routing totals

  7. treatmnt_treat() does the following:

    1. Removal is zero if there is no treatment equation in the input file

    2. Removal is zero if the inflow is zero

    3. Otherwise, it evaluates the treatment expression to find R[p] by calling getRemoval(p)

      1. getRemoval(p) which finds new quality in a node with no storage volume

    4. It then updates the nodal concentrations and mass balances by calling massbal_addReactedMass() and sets Node[j].newQual[p] = cOut

      1. massbal_addReactedMass() which adds mass reacted during the current time step to routing totals

2.4 Motivation for SWMM Water Quality Modifications

SWMM has some pollutant modeling limitations that must be addressed. The water quality module is limited by the range of treatment measures that can be modeled, specifically, limited nutrient treatment capabilities inside storage nodes (e.g., basins, wetlands). SWMM cannot simulate pollutant treatment inside links (e.g., conduits, channels) or pollutant generation processes (e.g., resuspension, erosion) inside any stormwater asset. Pollutant treatment cannot be turned on or off based on site conditions or other parameters, requiring treatment to run for the entire simulation. All of these constraints limit a user's ability to model complex pollutant transformations, necessitating a more generalizable and scalable approach.

3.0 Proposed Pollutant Process

Our idea is to extend the water quality modeling capabilities of SWMM by providing an API that enables users to inject and modify the pollutant in SWMM. We believe that this approach might be sustainable in the long run rather than augmenting the pollutant treatment models provided in SWMM. Through this API, users can read the inflow pollutant concentration into a node or link, model the pollutant transformation in a high-level programming language like Python, and set the computed values back into the SWMM. This approach relies upon the existing quality routing in SWMM's engine for transporting the pollutants through the stormwater network.

Modeling pollutant transformations in Python enables us to leverage the scientific computing stack to develop complex pollutant models. We plan on releasing a Python package (StormReactor) that builds on the newly added getters and setters to provide an easy to use toolkit for modeling pollutant transformations. StormReactor will provide implementations of CSTR, erosion, and other established pollutant transformation models for users to build on. We also plan on submitting a pull request to PySWMM with the pollutant getters and setters as well. StormReactor would then have PySWMM as a dependency.

In terms of specific use cases, we have updated the pollutant implementation in SWMM such that users can introduce pollutants in links to model flow driven erosion in channels. In nodes, as alluded above, users can model pollutant removal as a combination of reactors; this would vital in pollutant based control of stormwater systems. More details on the use case are provided in Section 4.0.

In terms of implementation details, we introduced extPollutFlag in nodes and links. This flag is pollutant and object (i.e., link or node) specific. It enables users to control which pollutant they want to model externally; the flag controls when SWMM skips its pollutant removal in favor of external pollutant values. Externally set pollutant values are accounted for in SWMM's mass balance functionality. In this implementation, all the pollutants can be considered as internal, which can be modified externally. More details on the code modifications are provided in Section 3.1.

3.1 Code Modifications

To address the limitations of SWMM's water quality module, we modified SWMM's C source code. An overview of the changes include:

  • Added new pollutant-based variables to objects.h to enable external pollutant handling [Section 3.1.1]

  • Allocate memory for the new pollutant-based variables in project.c

  • Added new getters to toolkitenums.h, toolkit.h, and toolkit.c to access pollutant-based variables [Section 3.1.2]

  • Added new setters to toolkit.h and toolkit.c to set pollutant concentrations [Section 3.1.3]

  • Added extPollutFlag as a global variable in global.h and initialize it by setting extPollutFlag = 0 in swmm5.c [Section 3.1.1]

  • Modified the functions qualrout_execute(), findLinkQual(), findStorageQual() in qualrout.c for external pollutant handling [Section 3.1.4]

  • Modified the function treatmnt_treat() in treatmnt.c to enable external pollutant handling [Section 3.1.5]

3.1.1 Variables

The following summarizes the variables and their type and description added in the SWMM source code.

  • extQual (int): External Quality State: Container for storing pollutant value that user wants to set in the node or link.
  • extPollutFlag(int): External Pollutant Flag: Flag that controls when the value in extQual is set in a node/link. If it is 0, pollutant is computed by native SWMM. When it is 1, pollutant in the node/link is set to the value in extQual.
  • hrt(double): Hydraulic Retention Time: HRT (i.e., water age) in a node is computed by SWMM. This exposes the computed HRT to users.
  • inQual (double): Inflow Quality State: Container for storing inflow quality state used in pollutant mass balance calculations.
  • reactorQual(double): Concentration in the Mixed Reactor: Container for storing reactor concentration used in pollutant mass balance calculations.

NOTE: hrt, inQual, and reactorQual variables were created to expose SWMM's internal states. The other two variables were created to switch between SWMM's pollutant routing and external pollutant handling.

NOTE: The pollutant-based variables added in SWMM's source code to enable external pollutant handling. All variables are defined in object.h except extPollutFlag, which is defined in global.h.

3.1.2 Getters

Three getters were added to toolkit.h, and toolkit.c in SWMM's source code. Getters enable users to access pollutant-based variables which are defined in toolkitenums.h in SWMM's source code. They can be obtained by calling the following functions:

  • swmm_getNodePollut(int index, SMNodePollut type, double **pollutArray, int *length)

    • index refers to the index of the node in the SWMM input file

    • type refers to the SMNodePollut property type code

    • **pollutArray and *length are not input parameters, they are the function output, which is an array of the pollutant data

  • swmm_getNodeResult(int index, SMNodeResult type, double *result)

    • index refers to the index of the node in the SWMM input file

    • type refers to the SMNodeResult property type code

    • *result is not an input parameter, it is the function output

  • swmm_getLinkPollut(int index, SMLinkPollut type, double **pollutArray, int *length)

    • index refers to the index of the node in the SWMM input file

    • type refers to the SMLinkPollut property type code

    • **pollutArray and *length are not input parameters, they are the function output, which is an array of the pollutant data

The following summarizes the getter functionality, specifically result property and description:

  • SMNODEQUAL: Get the current concentration in a node
  • SMNODECIN: Get the inflow concentration in a node
  • SMNODEREACTORC: Get the reactor concentration in a node
  • SMHRT: the hydraulic residence time (hours) in a node
  • SMLINKQUAL: Get the current concentration in a link
  • SMTOTALLOAD: Get the total quality mass loading in a link
  • SMLINKREACTORC: Get the reactor concentration in a link

NOTE: The pollutant-based variables that can be accessed by using the getters described above. These pollutant-based variables were added to toolkitenums.h in SWMM's source code.

3.1.3 Setters

Two pollutant concentration setters, swmm_setNodePollut and swmm_setLinkPollut, were added to toolkit.h and toolkit.c in SWMM's source code which enable a user to set the pollutant concentration in a node or link at any routing time step, respectively.

  • swmm_setNodePollut(int index, int pollutant_index, double pollutant_value)

    • index refers to the index of the node in the SWMM input file

    • pollutantindex refers to the index of desired pollutant in the SWMM input file

    • pollutantvalue refers to the value of the pollutant concentration you want to set in SWMM

  • swmm_setLinkPollut(int index, int type, int pollutant_index, double pollutant_value)

    • index refers to the index of the link in the SWMM input file

    • type needs to be set to SMLINKQUAL which sets the link's qual and allows accounting for loss and mixing calculation

    • pollutantindex refers to the index of desired pollutant in the SWMM input file

    • pollutantvalue refers to the value of the pollutant concentration you want to set in SWMM

When a setter is called, SWMM does the following:

  1. changes the global variable extPollutFlag from 0 to 1. extPollutFlag = 0 by default, meaning SWMM uses its traditional pollutant functionality. When extPollutFlag = 1, SWMM follows the external pollutant handling functionality as outlined in Sections 3.1.4 and 3.1.5.

  2. enacts the following depending on if it is a node or link setter:

    1. node: sets Node[index].extQual[pollutantindex] = pollutantvalue and Node[index].extPollutFlag[pollutantindex] = 1, turning on the extPollutFlag for the specific pollutant in the specific node

    2. link: sets Link[index].extQual[pollutantindex] = pollutantvalue and Link[index].extPollutFlag[pollutantindex] = 1, turning on the extPollutFlag for the specific pollutant in the specific link

NOTE: The setters added to toolkit.h and toolkit.c in SWMM's source code. The setters enables a user to change the pollutant concentration in a link or node at any time step.

3.1.4 qualrout.c

First, we modified findStorageQual() by adding Node[j].reactorQual[p] = c2 so that users can access the node reactor concentration using the appropriate getter [Section 3.1.2].

Then, we modified qualrout_execute(). In qualrout_execute() we check the value of the extPollutFlag. When extPollutFlag is set to 1, SWMM does the following for nodes:

  1. runs treatmnt_setInflow in treatmnt.c to compute the inflow concentrations to a node [Section 3.1.5]

  2. runs treatmnt_treat() in treatmnt.c to update the pollutant concentrations at a node after treatment [Section 3.1.5]

  3. at the end of the time step, sets extPollutFlag=0 so that if user does not want to have external control during the next time step, SWMM can reverts back to its traditional pollutant functionality

Finally, we modified findLinkQual(). Here is where the link concentration is set when the setter is used. We had to use this function because treatmnt.c is only for updating node concentrations. In findLinkQual(), we added lossExtQual, a new double variable which is the loss value for external water quality. When extPollutFlag is set to 1, SWMM does the following in findLinkQual():

  1. calculates lossExtQual

  2. inputs lossExtQual to the function massbaladdReactedMass()

  3. sets Link[i].newQual[p] = Link[i].extQual[p]

  4. resets the flag Link[i].extPollutFlag[p] = 0 so that if user does not want to have external control during the next time step, SWMM can reverts back to its traditional pollutant functionality

3.1.5 treatmnt.c

First off, treatmnt.c is only for nodes since SWMM only allows pollutant treatment in nodes. When Node[index].extPollutFlag[pollutantindex] = 1, SWMM follows the external pollutant handling functionality added to treatmnt.c. First, SWMM runs treatmnt_setInflow() which computes and saves array of inflow concentrations to a node. Then, SWMM runs treatmnt_treat(), doing the following:

  1. Sets R[p] = 0.0, meaning it will compute removal using a treatment equation provided in the input file for that pollutant in that node

  2. Sets cOut = Node[j].extQual[p]

  3. Calculates mass loss (as it normally would)

  4. Resets the flag to default SWMM treatment if (Node[j].extPollutFlag[p] == 1) Node[j].extPollutFlag[p] = 0

  5. Adds mass balance totals and revises nodal concentration (as it normally would)

It should be noted that we do not set newQual from an external function call directly, because we would have to re-implement mass balance. By updating cOut, we can build on SWMM's existing mass balance functionality. Furthermore, if we do not have the container variable extQual and have to set these pollutants at each step directly, we would have to interrupt the treatmnt_treat() function in qualrout.c, make an external function call to Python, read the variable, set it, and repeat this for every node/link and pollutant. By setting the pollutants in a container before the SWMM step is called, we can limit the number of external function calls we have to make, thus reducing the computational overhead.

4.0 Mass Balance

Pollutants injected using this API are accounted for in SWMM's mass-balance. When a pollutant is injected or removed in a node or link, it is accounted for in SWMM's mass balance using the same sub-routines used by SWMM for accounting pollutant interactions.

Mass Balance in links for external pollutants is handled in qualrout.c and findlinkqual function

Mass Balance in nodes for external pollutants is handled in treatmnt.c and treatmnt_treat function

5.0 Use Case

Users can access and modify pollutant concentrations in any node or link during any routing time step. Users can use new getters [Section 3.1.2] to access pollutant-based variables to model pollutant transformations and then set the new new pollutant concentration using new setters externally (e.g., Python) [Section 3.1.3].

Users can either:

  1. build their model directly in a Python script using the appropriate getters and setters; or

  2. add their new method to StormReactor's code base.

Example code snippet shows how a Python reactor model can be implemented for a node. The code uses getters to obtain both water quality and quantity variables which are then used to compute the new nitrate concentration. Then the setter is used to set the new node concentration in SWMM.

# Import required modules
from pyswmm import Simulation, Nodes, Links
from scipy.integrate import ode 
import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------#
# CSTR Equation
def CSTR_tank(t, C, Qin, Cin, Qout, V, k):
    dCdt = (Qin*Cin - Qout*C)/V - k*C
    return dCdt

#----------------------------------------------------------------------#
# Setup simulation
with Simulation("./NO.inp") as sim:
    # Get asset information
    Wetland = Nodes(sim)["93-49759"]
    
    # Setup dt calculation        
    start_time = sim.start_time
    last_timestep = start_time

    # Setup CSTR solver
    solver = ode(CSTR_tank)
    solver.set_integrator("dopri5")

    # Step through the simulation    
    for index,step in enumerate(sim):

        # Calculate dt
        current_step = sim.current_time
        dt = (current_step - last_timestep).total_seconds()
        last_timestep = current_step

        # Get wetland parameters        
        Wt_if = Wetland.total_inflow
        Wt_d = Wetland.depth
        Wt_v = Wetland.volume
        Wt_of = Wetland.total_outflow
        k_ni = 0.000087  # rate/5 sec

        # Get parameters to calculate NO
        Cin_NO = sim._model.getNodeCin("93-49759",0)
        Wetland_Cin.append(Cin_NO)

        #Solve Nitrate ODE
        if index == 0:
            solver.set_initial_value(0.0, 0.0)
            solver.set_f_params(Wt_if,Cin_NO,Wt_of,Wt_v,k_ni)
        else:
            solver1.set_initial_value(solver1.y, solver1.t)
            solver1.set_f_params(Wt_if,Cin_NO,Wt_of,Wt_v,k_ni)
        
        # Set new concentration
        sim._model.setNodePollutant("93-49759", 0, solver1.y[0])

    sim._model.swmm_end()

Example code snippet shows how a Python erosion model can be implemented for a link. The code uses getters to obtain both water quality and quantity variables which are then used to compute the new sediment concentration. Then the setter is used to set the new link concentration in SWMM. Details on adding a new pollutant method to StormReactor's code base can be found here.

# Import required modules
from pyswmm import Simulation, Nodes, Links
from scipy.integrate import ode 
import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------#
"""
Engelund-Hansen Erosion(1967)
Engelund and Hansen (1967) developed a procedure for predicting stage-
discharge relationships and sediment transport in alluvial streams.
        
w   = channel width (SI: m, US: ft)
So  = bottom slope (SI: m/m, US: ft/ft)
Ss  = specific gravity of sediment (for soil usually between 2.65-2.80)
d50 = mean sediment particle diameter (SI/US: mm)
d   = depth (SI: m, US: ft)
qt  = sediment discharge per unit width (SI: kg/m-s, US: lb/ft-s)
Qt  = sediment discharge (SI: kg/s, US: lb/s)
"""
#----------------------------------------------------------------------#
# Setup simulation
with Simulation("./TSS.inp") as sim:
    # Get asset information
    Channel = Links(sim)["95-69044"]
    
    # Setup dt calculation        
    start_time = sim.start_time
    last_timestep = start_time

    # Step through the simulation    
    for index,step in enumerate(sim):

        # Calculate dt
        current_step = sim.current_time
        dt = (current_step - last_timestep).total_seconds()
        last_timestep = current_step

        # Get channel parameters        
        Cin = sim._model.getLinkC2("95-69044", 0)
        Q = sim._model.getLinkResult("95-69044", 0)
        d = sim._model.getLinkResult("95-69044", 1)
        v = sim._model.getConduitVelocity("95-69044")
        w = 10.0
        So = 0.001
        Ss = 2.68
        d50 = 0.7

        # Erosion calculations   
        g = 9.81            # m/s^2
        ρw = 1000           # kg/m^3
        mm_m = 0.001        # m/mm
        kg_mg = 1000000     # mg/kg
        L_m3 =  0.001       # m3/L
        if v != 0.0:
            qt = 0.1*(1/((2*g*So*d)/v**2))*((d*So/((Ss-1)*d50))\
            *(1/mm_m))**(5/2)*Ss*ρw*((Ss-1)*g*(d50*mm_m)**3)**(1/2) # kg/m-s
            Qt = w*qt       # kg/s
        else:
            Qt = 0.0
        if Q != 0.0:
            Cnew = (Qt/Q)*L_m3*kg_mg    # mg/L
            Cnew = max(Cin, Cin+Cnew)
            # Set new concentration
            sim._model.setLinkPollutant("95-69044", "SM_LINKQUAL", 0, Cnew)

    sim._model.swmm_end()
Add a custom footer
Add a custom sidebar
Clone this wiki locally