Proposed Addition to SWMM5 Capability to Model Long Term Continuous Rainfall Dependent Infiltration

Joseph Pang

doi:https://doi.org/10.14796/JWMM.C369

Proposed Addition to SWMM5 Capability to Model Long Term Continuous Rainfall Dependent Infiltration

Joseph Pang (2014)

Seattle Public Utilities, USA

DOI: https://doi.org/10.14796/JWMM.C369

Comments

Abstract

The current SWMM5 engine contains three components that can be used for modeling rainfall dependent infiltration and inflow (RDII) into an underground conveyance system: the Surface Runoff, Unit Hydrograph (RTK method), and Groundwater modules. When using SWMM5 to model RDII, the inflow portion of RDII is commonly modeled by using Surface Runoff. The infiltration portion is usually modeled by using either the RTK method or Groundwater. The RTK method is designed to model infiltration which exhibits seasonal repeatable patterns. The Groundwater module is designed for modeling interflow groundwater infiltration into natural conveyance systems such as river channels and ditches, and conveyance systems with regular drainage inlets (e.g. a French drain). Currently, the SWMM5 engine does not have a method or module that can be used to model continuously, in a long term simulation, the infiltration portion of RDII that has the following characteristics: does not have repeatable seasonal patterns; is composed not only of interflow and baseflow but also of preferential and matrix flow; and flows into underground conveyance systems made with artificial materials (e.g. concrete, PVC) through irregular defects along the systems which may be partially or fully submerged by a local groundwater table. This paper proposes adding an additional method into the current SWMM5 engine that can be used to simulate such long term RDI. This additional method will be referenced as the SPU-LTCRDI method (Seattle Public Utilities—long term continuous rainfall dependent infiltration method).

1 Theory

The development of the SPU-LTCRDI (Seattle Public Utilities—long term continuous rainfall dependent infiltration) method is based upon:

the current SWMM5 Groundwater module;
a model described in section 2.5 of a paper titled, “Simplified Conceptual Structures and Analytical Solutions for Groundwater Discharge Using Reservoir Equations” (Rimmer and Hartmann 2012) which will be referenced as the Section 2.5 model hereinafter;
Darcy’s equation for groundwater seepage; and
observations from RDII flow monitoring data.

The Section 2.5 model was chosen because of its similarity to the existing SWMM5 Groundwater module which simplified the additional coding necessary to implement the SPU-LTCRDI method into the current SWMM5 engine. During the development of the SPU-LTCRDI method, the aquifer configuration and discharge equations of the current SWMM5 Groundwater module and the Section 2.5 model were combined and subsequently modified for the SPU-LTCRDI method. The following sections describe the aquifer configuration and the groundwater flow components of the SPU-LTCRDI method. As the SWMM5 engine allows for both metric and U.S. traditional units, all variables presented in this paper adopt the same units as their equivalent in the current SWMM5 engine.

1.1 The aquifer Configuration

The current SWMM5 Groundwater module is a single reservoir system. The single reservoir represents an aquifer with two groundwater zones: the unsaturated (vadose) zone and the saturated zone. It contains two groundwater discharge outlets: one outlet to a user specified receiving conveyance system, and the other loses water from the model. .

Figure 1 Configuration of the current SWMM5 aquifer (Rossman, 2010).

Figure 1 shows a schematic of the current SWMM5 aquifer configuration, where (Rossman, 2010):

d_U	=	upper zone depth,
d_L	=	lower zone depth,
d_TOT	=	total depth of aquifer,
f_I	=	infiltration from the surface,
f_EU	=	evapotranspiration from the upper zone which is a fixed fraction of the unused surface evaporation,
f_U	=	percolation from the upper to lower zone which depends on the upper zone moisture content θ and d_U,
f_EL	=	evapotranspiration from the lower zone, which is a function of d_U,
f_L	=	percolation from the lower zone to deep groundwater which depends on d_L and d_TOT, and
f_G	=	lateral groundwater interflow to the drainage system, which depends on the lower zone depth d_L as well as the depth in the receiving channel or node.

Like SWMM5’s aquifer, the aquifer of the Section 2.5 model also has one reservoir but it has two outlets to the receiving conveyance system with none to lose groundwater from the model. Figure 2 shows a schematic of the aquifer in the Section 2.5 model.

Figure 2 Configuration of the Section 2.5 model (Rimmer and Hartmann, 2012).

Before the development of the SPU-LTCRDI method, the author had the opportunity to be involved in about 115 RDII calibrations using the current SWMM5 Surface Runoff and Groundwater modules over the past couple years. Furthermore, in addition to using SWMM5 for RDII modeling, the author also had the opportunity to use other hydrologic and hydraulic modeling software (e.g. InfoWorks and DHI MOUSE) to model RDII since 2002. Experience with RDII calibration using different modeling software and comparing the architectures of the two groundwater models suggest to the author that for the SWMM5 aquifer, adding a second outlet from the aquifer to the receiving conveyance system would enhance SWMM5’s ability to simulate the infiltration component of RDII for long term continuous simulations. Thus, for the SPU-LTCRDI method, its aquifer configuration combines features in the current SWMM5 aquifer with those of the Section 2.5 model aquifer.

As in SWMM5’s current aquifer configuration, the SPU-LTCRDI aquifer is also a single reservoir system. It contains the same two groundwater zones (vadose and saturated) as the current SWMM5 aquifer but has two outlets to a receiving conveyance system rather than one.

Figure 3 Configuration of new aquifer.

Figure 3 shows the configurations of the SPU-LTCRDI aquifer,
where:

H₀	=	aquifer bottom elevation, a user defined elevation which will be discussed in detail later in this paper,
H_sw	=	water depth above aquifer bottom in the receiving conveyance system,
H^*	=	water depth above aquifer bottom above which the 2nd outlet becomes active, and
Z_aq	=	depth offset between aquifer bottom and invert elevation of the receiving conveyance system; this is also the depth of water in the saturated zone above which infiltration to the receiving conveyance system begins (this parameter is not a user defined parameter: it is internally calculated based upon aquifer bottom elevation and invert elevation of the receiving node).

In applications, the new SPU-LTCRDI aquifer is parameterized primarily (with a couple of additions) by the same set of parameters as the current SWMM5 aquifer with some of the parameters having different definitions when they are used in this new aquifer. The parameters with different definitions between the two aquifers are the moisture content parameters, the lower groundwater (gw) loss rate, and aquifer bottom elevation.

Taking the approach of using the same set of parameters does not only simplify the modifications needed to implement this new aquifer configuration into the existing SWMM5 code, but also ensures that the effort needed to use this aquifer would be approximately the same as using the existing aquifer configuration. However, unlike the current SWMM5 Groundwater module, an SPU-LTCRDI aquifer can only be shared among RDI generating subcatchments that have the same surface elevation, aquifer bottom elevation and discharge location. Subcatchments that have different surface and aquifer bottom elevations, or discharge locations must have their own aquifers.

The following describes in detail the differences between how the aforementioned parameters are used in the SPU-LTCRDI aquifer and the existing SWMM5 aquifer.

The Moisture Content Parameters

In the current SWMM5 Groundwater module, an aquifer of a subcatchment takes the shape of a prism below the surface of the subcatchment. The surface area of the prismatic aquifer is the same as its associated subcatchment. However, observations from RDII flow monitoring data and general soil seepage knowledge revolving around flow nets suggest that the potential volume of soil contributing to groundwater infiltration into an underground conveyance system would most likely be limited to the soil around the system, and take an elliptical trough shape similar to that of a flow net. Figure 4 shows a sketch of a flow net around two underground conveyance systems, illustrating the potential volume of soil influencing groundwater infiltration.

Figure 4 Flow net around underground conveyance system.

In order to limit the potential volume and surface area of soil contributing to groundwater infiltration without doing a major modification to the mass balance component of the existing SWMM5 engine, a SPU-LTCRDI aquifer also uses the same prismatic aquifer with the same surface area as its associated subcatchment. However, the following current SWMM5 aquifer parameters will be defined differently in the SPU-LTCRDI aquifer.

Porosity;
Field Capacity;
Wilting Point;
Unsaturated Zone Moisture;
Lower GW Loss Rate; and
Bottom Elevation

The following describes how these parameters differ and should be adjusted and used. The definitions of the rest of the aquifer parameters between the two aquifer configurations remain the same.

Porosity, Field Capacity, Wilting Point

The porosity of soil is calculated by using the following equation:

italic eta equal V subscript v o i d end subscript over V subscript t o t a l end subscript

(1)

where:

η	=	porosity,
V_void	=	volume of void, and
V_total	=	total volume of soil.

For continuity calculations, V_void is the variable of focus as it defines the actual volume in soil where water is stored and drained. As the porosity relationship holds true in a soil column of any size within a subcatchment, it can be used to express the volume of void in soil around an underground conveyance system in terms of the total volume of soil under an entire subcatchment. Such a relationship is critical if the new aquifer configuration is to use the same surface area as its associated subcatchment but still be able to limit the amount of soil contributing to groundwater infiltration. The volume of void in soil localized around an underground conveyance system can be calculated from the following relationship:

italic eta equal open parentheses V subscript v o i d end subscript over V subscript t o t a l end subscript close parentheses subscript l o c a l i z e d end subscript

(2)

Expressing (V_total)_localized as a fraction (r) of (V_total)_subcatchment:

open parentheses V subscript t o t a l end subscript close parentheses subscript l o c a l i s e d end subscript equal r asterisk times open parentheses V subscript t o t a l end subscript close parentheses subscript s u b c a t c h m e n t end subscript space

(3)

(V_void)_localized can be calculated by the following:

open parentheses V subscript v o i d end subscript close parentheses subscript l o c a l i s e d end subscript equal italic eta asterisk times r asterisk times open parentheses V subscript t o t a l end subscript close parentheses subscript s u b c a t c h m e n t end subscript space space

(4)

Then the porosity value to be used in a SPU-LTCRDI aquifer can be calculated by the following equation:

open parentheses V blank presubscript v o i d end presubscript close parentheses subscript l o c a l i z e d end subscript over open parentheses V subscript t o t a l end subscript close parentheses subscript s u b c a t c h m e n t end subscript equal italic eta asterisk times r equal italic eta subscript S P U minus L T C R D I end subscript

(5)

Thus, η * r instead of η should be used as the porosity value when the SPU-LTCRDI aquifer is used so that the correct volume of void is obtained when the subcatchment area is multiplied by the unit area volume and the appropriate groundwater level is simulated.

The same logic follows for Wilting Point and Field Capacity and an r-value should also be applied to Wilting Point and Field Capacity. As soil underneath a subcatchment is most likely non-homogeneous (e.g. the presence of backfills, soil pores filled with organic matter, voids created by worms and other animals), the r-value for Wilting Point and the r-value for Field Capacity need not to be the same as the r-value for Porosity. In fact, each of these moisture content parameters can have its own r-value. The following summarizes the equations that can be used to calculate Porosity, Wilting Point, and Field Capacity when they are used in a SPU-LTCRDI aquifer:

italic eta subscript S P U minus L T C R D I end subscript equal italic eta asterisk times r subscript italic eta

(6)

W P subscript S P U minus L T C R D I end subscript equal W P asterisk times r subscript W P end subscript

(7)

F C subscript S P U minus L T C R D I end subscript equal F C asterisk times r subscript F C end subscript

(8)

where:

η_SPU-LTCRDI	=	Porosity for SPU-LTCRDI aquifer,
WP_SPU-LTCRDI	=	Wilting Point for SPU-LTCRDI aquifer,
FC_SPU-LTCRDI	=	Field Capacity for SPU-LTCRDI aquifer,
WP	=	Wilting Point,
FC	=	Field Capacity,
r_η	=	r-value for Porosity,
r_WP	=	r-value for Wilting Point, and
r_FC	=	r-value for Field Capacity.

None of the r-values are currently part of the SWMM5 parameter set. Therefore, η_SPU-LTCRDI, WP_SPU-LTCRDI and FC_SPU-LTCRDI would need to be calculated outside of the SWMM5 environment and then entered into the SPU-LTCRDI aquifer as Porosity, Wilting Point and Field Capacity respectively. These r-values are the only parameters used by the SPU-LTCRDI method that are not existing SWMM5 parameters.

The η_SPU-LTCRDI, WP_SPU-LTCRDI and FC_SPU-LTCRDI and their associated r-values can be determined through a variety of methods. Initially, η_SPU-LTCRDI, WP_SPU-LTCRDI and FC_SPU-LTCRDI can be set equal to their equivalents in a current SWMM5 aquifer. Typical values of these parameters for a current SWMM5 aquifer can be found in the current SWMM5 user manual. In such a case, the r-values would be 1. Then the r-values can be adjusted through calibration. If the area of trenches along an underground conveyance system is known, the r-value for porosity can be initially estimated by using the ratio between the sum of the area of trenches and the subcatchment area. In that case, the r-values of the corresponding wilting point and field capacity parameters could also be set to be equal to the r-value of porosity and then adjusted during calibration. If there are different types of soil underneath a subcatchment, or one type of soil is underneath several subcatchments, then the r-values would be determined by the soil types rather than subcatchments. In cases when one soil type is under multiple subcatchments, the initial r-value estimates can be determined by:

fraction numerator sum open parentheses V subscript t o t a l end subscript close parentheses subscript l o c a l i z e d end subscript over denominator sum open parentheses V subscript t o t a l end subscript close parentheses subscript s u b c a t c h m e n t end subscript end fraction equal r

(9)

In cases when there are multiple soil types under one subcatchment, new moisture content values rather than the r-values would need to be determined for the subcatchment. For example, the SPU-LTCRDI Porosity value can be determined by the following equation. The equations for Wilting Point and Field Capacity are analogous to those of Porosity.

italic eta subscript S P U minus L T C R D I end subscript equal fraction numerator begin display style sum from i minus 1 to k of end style italic eta subscript i r subscript i left parenthesis V subscript t o t a l end subscript right parenthesis subscript i over denominator begin display style sum from i equal 1 to k of end style left parenthesis V subscript t o t a l end subscript right parenthesis subscript i end fraction

(10)

where:

i	=	ith soil type in the subcatchment, and
k	=	total number of soil types in the subcatchment.

Unsaturated Zone Moisture

The Unsaturated Zone Moisture parameter is used to set the initial moisture content of the unsaturated zone. Its value is set after determining η_SPU-LTCRDI, WP_SPU-LTCRDI and FC_SPU-LTCRDI as it must lie between or be equal to Wilting Point and Porosity. For simulations with a short duration between the start of simulation and when simulation results are compared with flow monitoring data (e.g. from a couple of months to 2 y), this parameter can be used as one of the calibration parameters as it would affect how well the simulated results match with flow monitoring data for the initial wet weather events. For simulations with a longer duration between the start of simulation and when simulation results are compared with flow monitoring data (e.g. ≥2 y), this parameter can be set to be equal to Wilting Point.

Because of their effects on water depth in the saturated zone, experience has shown that it is fairly straight forward to determine the values of the moisture content parameters iteratively during calibration. The key is to ensure that the flow monitoring data used for calibration contains hydrographs from rainfall events of various intensities and durations, and then observe during calibration how these parameters change the shape of a simulated hydrograph together with other parameters in the SPU-LTCRDI method.

Lower GW Loss Rate

A new algorithm for calculating the Lower GW Loss Rate has been developed for the SPU-LTCRDI method. The current SWMM5 aquifer calculates the lower groundwater loss rate to deeper aquifer (f_L) by the following equation:

f subscript L equal D P asterisk times d subscript L over d subscript T O T end subscript

(11)

where:

user provided Lower GW Loss Rate parameter.

The equation suggests that groundwater percolation to deeper aquifer and lost from the model (f_L) increases as the depth of the saturated zone (d_L) increases. However, for RDII application, observation from flow monitoring data suggests that the deep groundwater loss rate decreases as the volume of percolated water increases. Such observation is consistent with the Green–Ampt method (as shown by Equation 12) that infiltration capacity (f_P) decreases and converges to K_sat as cumulative infiltration F(t) increases.

f subscript p open parentheses t close parentheses equal K subscript s a t end subscript asterisk times open parentheses fraction numerator italic psi increment italic theta over denominator F open parentheses t close parentheses end fraction plus 1 close parentheses

(12)

where:

ψ	=	suction head, and
θ	=	moisture content.

As a result, the SPU-LTCRDI aquifer also assumes the percolation rate to deeper groundwater aquifer to decrease as the volume of percolated water increases. Staying within the constraint of using only parameters and variables which are used by the current SWMM5 Lower GW Loss Rate equation, yet reflecting the inverse relationship between Deep Groundwater Loss Rate and Total Percolated Water as discussed above, the following empirical equations and algorithm are developed:

f subscript L equal D P asterisk times e to the power of open parentheses 1 minus d subscript L over d subscript T O T end subscript close parentheses end exponent

(13)

If f_L * ?t > d, then:

f subscript L equal fraction numerator d subscript L over denominator increment t end fraction

(14)

The exponential term allows f_L to be initially higher than DP. But as percolation volume increases (assuming the ratio of −dL/−dTOT to deeper aquifer is equal to dL /dTOT in aquifer), f_L decreases and converges to DP, analogous to how f_p(t) decreases and converges to K_sat as F(t) increases. Equation 14 sets the upper bound of the loss rate to maintain continuity.

Aquifer Bottom Elevation

As discussed earlier in this paper, the potential volume of soil contributing to infiltration of groundwater into an underground conveyance system is most likely limited to the soil around the system and takes an elliptical trough shape. However, after water is infiltrated into the soil, not all of the water will be infiltrated into the receiving conveyance system due to the configuration and the location of defects along the system. The water would need to be in the path of flow channels of flow nets formed between the soil and the conveyance system defects in order for the water to be infiltrated into the conveyance system. For the portion of water in the soil that does not get into a flow channel, it may have some effect on the overall hydraulics of infiltration but it continues to infiltrate deeper into the soil and is lost from the conveyance system.

As mentioned earlier, the shape of a SPU-LTCRDI aquifer adopts the same vertical prism shape as the current SWMM5 aquifer rather than a general elliptical trough shape similar to a flow net. Also, in a model, all water that is infiltrated into an aquifer has an effect on the rate and volume of water infiltrated into a conveyance system. To simulate the effect of water in the soil around a conveyance system that does not get infiltrated into the system (as discussed above) and the effect of simulating elliptical trough soil by a prismatic soil aquifer, the volume of soil between the bottom of a SPU-LTCRDI aquifer and the invert of its discharging node is used. This part of the aquifer is used because, as discussed in the following sections, the equation algorithms used to simulate infiltration from an aquifer into a conveyance system are functions of saturated zone depth in the aquifer and water depth inside the receiving conveyance system measured from the same datum. Water below the invert of a discharging node has some effect on the amount of water infiltrated into a conveyance system but the water itself does not get infiltrated into the system. It is lost from the system at a rate controlled by the Lower GW Loss Rate and Lower Evap. Depth. This water bears some similar properties and characteristics to the water around a conveyance system that does not get infiltrated into the conveyance system as described above.

As a result, the Bottom Elevation of an aquifer is a calibration parameter. It should always be lower than the invert of its discharging node. Other field information such as the field elevation of an impervious layer below the conveyance system or the elevation of a long term water table should also be considered if known when setting the bottom elevation of a SPU-LTCRDI aquifer. A good starting value can be the invert of the discharge node minus half of the average of the maximum depths of nodes in the associated subcatchment. However, ultimately, its final value is determined from calibration.

Effects of Aquifer Parameters on Groundwater Discharge

It should be noted that aquifer parameters have a direct impact on the groundwater discharge rate from an aquifer. As mentioned above and shown in the following sections, this is because the groundwater discharge equation algorithm is a function of the depth of water in an aquifer and in the receiving conveyance system. Thus, during calibration, time should be spent on sizing an aquifer appropriately, through choosing appropriate aquifer parameter values such that an aquifer can work together with the discharge equation algorithm to generate simulation results that replicate flow monitoring data in a long term continuous simulation with various hydrologic conditions (e.g. varying evaporation rates and rainfall events of various intensities and durations). If an aquifer is not sized appropriately, experience shows that it is very difficult, if not impossible, to continuously match simulated results to flow monitoring data regardless of what parameter values are used in the discharge equation algorithm. Such behavior of the SPU-LTCRDI method helps to prevent simulated results from being force-fitted onto flow monitoring data. To help size an aquifer properly, it is recommended that the duration between the start of a simulation and when the simulated result would be compared with flow monitoring data to be as long as possible with the use of an appropriate value for the Unsaturated Zone Moisture parameter.

Figure 5 Preferential flow (Cornell University, 2013).

Figure 6 Matrix flow (Cornell University, 2013).

Figure 7 Types of groundwater infiltration.

1.2 The Groundwater Flow Component

RDI is mainly generated by two sources: infiltration from surface water percolation and infiltration from the local groundwater table. After rain has fallen onto a pervious area, if the hydraulic gradient for water infiltration into the soil is >0, rain water infiltrates into the vadose zone of the soil. Depending on the moisture content and matrix framework of the soil, rain water can either stay in the soil matrix or continue to percolate downward and infiltrate into a conveyance system in the form of the quicker preferential flow and slower Green–Ampt matrix flow.

Concurrently, if the local groundwater table is high enough, groundwater from the local groundwater table can also infiltrate into a conveyance system in the form of baseflow infiltration. Figure 7 illustrates all of the aforementioned infiltration processes into a conveyance system. The water can also be lost from the system as discussed above. Figure 5 and Figure 6 illustrate the percolation of preferential flow and matrix flow. Interflow refers mainly to the lateral movement of matrix flow in the vadose zone. Thus, to simplify terminology, matrix flow and interflow will be referenced together as matrix–interflow hereinafter unless their distinction is required.

Current SWMM5 Groundwater Discharge Equation Algorithm

If preferential flow, matrix–interflow, and baseflow infiltration processes into a conveyance system are to be simulated individually, a 2D application of Darcy’s Equation would need to be solved for each type of infiltration in their respective vadose or saturated zones around an underground conveyance system. In the current SWMM5 Groundwater module, only the vertical 1D infiltration and percolation of matrix flow through the soil is simulated, not the quicker preferential flow infiltration. Under the current scheme in SWMM5, after water has infiltrated through the surface, if an aquifer is present, the infiltrated water is collected in the vadose zone of the aquifer. After the field capacity of the vadose zone is exceeded, the infiltrated water will percolate into the saturated zone of the aquifer. Water from the saturated zone then either flows into a user specified receiving conveyance system in the form of interflow infiltration or loses from the model as percolation to deep groundwater.

The SWMM5 Groundwater module uses a six-parameter groundwater discharge equation to simulate the infiltration of groundwater into a conveyance system. Equation 15 shows the current SWMM5 groundwater discharge equation. Details of this equation can be found in the current SWMM5 user manual (Rossman, 2010).

f subscript G equal A 1 open parentheses d subscript L minus H asterisk times close parentheses to the power of B 1 end exponent minus A 2 open parentheses H subscript S W end subscript minus H asterisk times close parentheses to the power of B 2 end exponent plus A 3 d subscript L H subscript S W end subscript

(15)

where:

A1, B1, A2, B2, A3	=	user defined SWMM5 Groundwater discharge equation calibration parameters, and
H^*	=	user defined threshold groundwater depth above aquifer bottom, above which groundwater discharge into a user specified conveyance system begins.

Figure 1 above shows a schematic of this infiltration process. Baseflow infiltration in the current SWMM5 engine is simulated in the form of inflow at a junction and is not part of the Groundwater module, although the Groundwater module can be configured to simulate baseflow infiltration, if needed, through configuring the aquifer Bottom Elevation, Water Table Elevation and Lower GW Loss Rate parameters.

SPU-LTCRDI Groundwater Discharge Equation Algorithm

The SPU-LTCRDI method uses the same vertical 1D infiltration and percolation processes as the current SWMM5 Groundwater module to simulate the groundwater infiltration and percolation processes. But for the discharge equation algorithm, it uses a different algorithm derived from the current SWMM5 groundwater discharge equation algorithm, the discharge algorithm of the Section 2.5 model, Darcy’s Equation of groundwater seepage, the current SWMM5 unsaturated hydraulic conductivity equation K(θ), and observation from flow monitoring data. The SPU-LTCRDI groundwater discharge equation algorithm can be summarized by the follow equation:

Q subscript R D I end subscript equal Q subscript P M I end subscript plus Q subscript B

(16)

where:

Q_RDI	=	RDI,
Q_PMI	=	preferential flow and matrix–interflow (PMI) infiltration, and
Q_B	=	baseflow infiltration.

Dividing flow by area normal to flow, Equation 16 can be expressed as:

q subscript R D I end subscript equal q subscript P M I end subscript plus q subscript B

(17)

where:

q equal Q over A

(18)

The SPU-LTCRDI groundwater discharge equation algorithm uses the same six parameters as the current SWMM5 groundwater discharge equation, but each with a different application. Using the same infiltration and percolation processes and the same number of parameters minimizes the modification of coding needed to implement the discharge equation algorithm into the existing SWMM5 code. It also helps to keep the amount of effort needed to use the SPU-LTCRDI method to be approximately the same as using the existing SWMM5 Groundwater module. The following sections describe how each term of Equation 17 is developed.

Preferential Flow and Matrix–Interflow Infiltration (q_PMI)

Although the SPU-LTCRDI method uses the same vertical 1D infiltration and matrix flow percolation processes as the current SWMM5 Groundwater module, the SPU-LTCRDI method attempts to model all of the preferential flow, matrix–interflow and baseflow infiltration into a conveyance system. During development, each component of Equation 17 was developed starting with Darcy’s Equation:

q equal K subscript s a t end subscript asterisk times i with rightwards arrow on top

(19)

where:

q	=	seepage flow per unit area,
K_sat	=	saturated hydraulic conductivity, and
$i with rightwards arrow on top$	=	hydraulic gradient across two points along a streamline in soil.

In order to avoid solving a 2D seepage analysis, as mentioned earlier, the hydraulic gradient term of the Darcy’s Equation (a vector quantity in x– and z– directions) was split into two scalar quantities for simulating q_PMI in the SPU-LTCRDI method. These two scalar quantities were: a new coefficient called coefficient of preferential flow and matrix–interflow (PMI(θ)); and G as shown in Equations 20 and 21.

i equal P M I open parentheses italic theta close parentheses asterisk times G

(20)

Equation 21 below shows Darcy’s Equation expressed in terms of PMI(θ) and G:

q equal K subscript s a t end subscript asterisk times P M I open parentheses italic theta close parentheses asterisk times G

(21)

PMI(θ) was used to imitate the faster preferential flow infiltration within the constraint of a matrix flow percolation process. The equation of PMI(θ) was empirically developed based upon the exponential term of the current SWMM5 unsaturated hydraulic conductivity equation K(θ). The difference between the exponential term in that equation and what was in PMI(θ) was at the exponent of the exponential term. The exponential term of the SWMM5 K(θ) equation is shown below:

e to the power of open parentheses italic theta minus italic eta close parentheses end exponent asterisk times H C O

(22)

where:

HCO

current SWMM5 hydraulic conductivity slope.

This exponential term says that the value of the exponential increases and approaches 1 as moisture content increases and approaches porosity. Observations from Figures 5 and 6 above suggest that preferential flow travels faster than matrix flow. Thus infiltration of preferential flow would cause q to start sooner and have larger magnitude than what otherwise matrix flow percolation alone would predict. As the moisture content of the unsaturated zone increases, matrix flow increases until eventually the wetting front of the matrix flow catches up with that of preferential flow and the overall wetting front in the unsaturated zone reaches the saturated zone. At that point, the soil becomes fully saturated. By modifying the exponent of Equation 22 by multiplying the exponent by −1, a term can be obtained and used in the q_PMI equation, together with the aquifer field capacity parameter, to imitate the process described above. Equation 23 shows the equation for PMI(θ).

P M I open parentheses italic theta close parentheses equal e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent

(23)

where:

HCO

≥

When implemented into Darcy’s Equation as shown by Equation 24, the exponential term of the equation allows q to be large at the beginning when the θ of the soil matrix is small (preferential flow), and q to decrease to matrix flow when soil is saturated (η = θ).

q equal K subscript s a t end subscript asterisk times e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times G

(24)

The equation for G was empirically developed based upon observing the shapes of field RDII hydrographs—especially the shapes of the falling and recession limbs and their relationships with the rest of the hydrographs—under evaporation and rainfall events of various intensities and durations. By plotting flow on a logarithmic scale and noting the locations of where changes of slopes occurred along the falling and recession limbs, as shown in Figure 8, observation from field data suggested that the shapes and magnitudes of preferential and matrix–interflow infiltration components of a RDII hydrograph could be simulated by a combination of exponential and power functions.

Figure 8 Components of RDII.

Other infiltration equations also supported the use of exponential form (e.g. Horton) and power form (e.g. the current SWMM5 groundwater discharge equation) to describe infiltration flow. Thus, borrowing from the idea of the exponential form of Horton equation, the power function form of the current SWMM5 groundwater discharge equation, and field data observation, the equation for G was developed with an exponential term and a power term as shown below (in matrix notation for clarity):

when d_L ≤ H^*,

G equal open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row 0 row cell h to the power of B subscript 2 end exponent end cell end table close square brackets

(25a)

else when d_L > H^* and H_sw ≥ H^*,

G equal open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of open parentheses B 1 asterisk times h close parentheses end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets

(25b)

else when d_L > H^* and H_sw < H^*,

G equal open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of open parentheses B 1 asterisk times h close parentheses end exponent minus h subscript L asterisk times e to the power of open parentheses B 1 asterisk times h subscript L close parentheses end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets

(25c)

where:

h	=	d_L − H_sw,
h_L	=	H^* − H_sw,
A1, B1, A2, B2	=	user defined calibration parameters, which are the same as those in the current SWMM5 groundwater discharge equation except that they are utilized differently in the SPU-LTCRDI method as shown above, and
H^*	=	user defined parameter which is the same parameter as the current SWMM5 groundwater discharge equation except it is being redefined as the depth above aquifer bottom above which the exponential term and the second outlet in the SPU-LTCRDI aquifer become active.

In addition to simplifying a 2D seepage problem to a 1D problem, Equations 25a to 25c with all of their parameters also attempted to capture in a simplistic manner the heuristic aspects of RDI. Different from natural channel infiltration, the process of infiltration into an underground conveyance system made with artificial materials (e.g. concrete, PVC) through irregular defects along the system contains two components: physical and heuristic. The physical components are infiltration from the surface and percolation through soil, both of which can be described well by seepage physics. The heuristic component is how the hydraulic gradient needed for groundwater to infiltrate from soil into an artificially constructed conveyance system is affected by the material (e.g. concrete, PVC) and condition (e.g. new, old) of the system. For example, a new PVC conveyance system under a high groundwater table would most likely experience less groundwater infiltration than an old concrete conveyance system under a lower groundwater table. This is because the PVC conveyance system would most likely contain fewer defects along the system than the old concrete conveyance system. Such material and condition dependent phenomenon cannot be easily described if at all possible by using soil seepage physics principles. This results in the use of a more empirical approach for developing the equation for G. In summary, the equation for q_PMI is shown below:

when d_L ≤ H^*,

q subscript P M I end subscript equal K subscript s a t end subscript asterisk times e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row 0 row cell h to the power of B 2 end exponent end cell end table close square brackets

(26a)

else when d_L > H^* and H_sw ≥ H^*

^{$q subscript P M I end subscript equal K subscript s a t asterisk times end subscript e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of open parentheses B 1 asterisk times h close parentheses end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets$}

(26b)

else when d_L > H^* and H_sw < H^*,

q subscript P M I end subscript equal K subscript s a t end subscript asterisk times e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of open parentheses B 1 asterisk times h close parentheses end exponent minus h subscript L asterisk times e to the power of open parentheses B 1 asterisk times h subscript L close parentheses end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets

(26c)

where:

HCO

≥

Baseflow Infiltration (q_B)

As soil is saturated under a groundwater table, a baseflow infiltration equation can be developed by starting with the flow net equation:

q equal K subscript s a t end subscript open parentheses N subscript f over N subscript d close parentheses increment H

(27)

where:

N_f	=	number of flow channels,
N_d	=	number of equipotential drops, and
?H	=	difference in total head between the first and last equipotentials.

Based upon the flow net equation, the baseflow infiltration component of the SPU-LTCRDI discharge equation algorithm is shown as:

q subscript B equal K subscript s a t end subscript asterisk times A 3 asterisk times h

(28)

where:

a user defined calibration parameter, which is the same as that in the current SWMM5 groundwater discharge equation except that it is utilized differently in the SPU-LTCRDI method, as shown above.

Since we generally do not know N_f?/N_d, this ratio is represented by the calibration parameter A3. In addition, A3 is also used as a lump parameter for incorporating the heuristic effect of the condition of sewers on baseflow infiltration.

Final Form of the Groundwater Discharge Equation Algorithm

The following summarizes altogether the different components of the SPU-LTCRDI discharge equation algorithm:

when H_sw ≥ d_L,

q equal presubscript R D I end presubscript 0

(29a)

else when d_L ≤ H^*,

q subscript R D I end subscript equal K subscript s a t end subscript open curly brackets e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row 0 row cell h to the power of B 2 end exponent end cell end table close square brackets plus A 3 asterisk times h close curly brackets

(29b)

else when d_L > H^* and H_sw ≥ H^*,

q subscript R D I end subscript equal K subscript s a t end subscript open curly brackets e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of left parenthesis B 1 asterisk times h right parenthesis end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets plus A 3 asterisk times h close curly brackets

(29c)

else when d_L > H^* and H_sw < H^*,

q subscript R D I end subscript equal K subscript s a t end subscript open curly brackets e to the power of open parentheses italic eta minus italic theta close parentheses asterisk times H C O end exponent asterisk times open square brackets space A 1 space space space A 2 space close square brackets open square brackets table row cell h asterisk times e to the power of open parentheses B 1 asterisk times h close parentheses end exponent minus h subscript L asterisk times e to the power of open parentheses B 1 asterisk times h subscript L close parentheses end exponent end cell row cell h to the power of B 2 end exponent end cell end table close square brackets plus A 3 asterisk times h close curly brackets

(29d)

2 Sample Results

The ability of the SPU-LTCRDI method to simulate RDI was tested in two long term continuous simulations of 3.04 y and 2.45 y. The inflow component of RDII was simulated using the current surface runoff component of SWMM5 as discussed at the beginning of this paper, and the infiltration component of RDII was simulated using both the current SWMM5 Groundwater module and the SPU-LTCRDI method. Results were generated independently for comparison. These two long term continuous simulations were conducted each with its own model basin: model basin 1 for the 3.04 y simulation, and model basin 2 for the 2.45 y simulation. Seattle Public Utilities uses four goodness of fit (GOF) measures to evaluate the quality of a calibration (MGS, 2010). They are:

GOF peak flow: the closer the value to 100%, the better the match between measured and simulated peak flow. The equation of GOF peak flow is as follow:

G O F subscript P e a k end subscript equal 100 asterisk times open parentheses fraction numerator P e a k D i s c h a r g e subscript s i m u l a t e d end subscript over denominator P e a k D i s c h a r g e subscript r e c o r d e d end subscript end fraction close parentheses

(30)

GOF flow volume: the closer the value to 100%, the better the match between measured and simulated volume. The equation of GOF flow volume is as follow:

G O F subscript V o l u m e end subscript equal 100 asterisk times open parentheses fraction numerator D i s c h a r g e V o l u m e subscript s i m u a t e d end subscript over denominator D i s c h a r g e V o l u m e subscript r e c o r d e d end subscript end fraction close parentheses

(31)

Standardized bias (SB): this is a standardized measure of the average amount of under or over simulation, expressed as a percentage. The closer the value to 0%, the better the match between measured and simulated flow. The equation of standardized bias is:

S B equal 100 asterisk times open square brackets fraction numerator 1 over n begin display style sum from i equal 1 to n of end style open parentheses r e s i d u a l subscript i close parentheses over denominator stack Q r e c with bar on top end fraction close square brackets

(32)

where:

residual_i	=	Qsim_i − Qrec_i,
Qsim_i	=	simulated flow at ith time step,
Qrec_i	=	recorded flow at ith time step, and
$top enclose Q r e c end enclose$	=	average recorded flow.

Standardized root mean square error (SRMSE): this is a standardized measure of the variability of the difference between simulated and recorded flows expressed as a percentage. The closer the value to 0%, the better the overall match between simulated and recorded flow. The equation of SRMSE is:

S R M S E equal 100 asterisk times fraction numerator square root of begin display style sum from i equal 1 to n of end style open parentheses begin display style fraction numerator r e s i d u a l subscript i to the power of 2 over denominator n end fraction end style close parentheses end root over denominator stack Q r e c with bar on top end fraction

(33)

In addition to the four GOF measures discussed above, a fifth GOF is added for this paper: correlation coefficient (R) as defined by the CORREL function in MS Excel. For a correlation coefficient, the closer the value to 1, the better the match between simulated and measured flow.

These GOF measures are used to evaluate the quality of the results from the two sets of simulations against flow monitoring data time series for various wet weather events as well as the whole time series.

2.1 Model Basin 1

Table 1 summarizes the characteristics of model basin 1 and the total duration of simulation and observed data time series. Figure 9 shows a plan view of the basin model. Table 2 summarizes the characteristics of the overall observed data and wet weather events used for comparison.

Figures 10 to 14 summarize the aforementioned five statistics of goodness of fit measures for two independent calibrations of the 3.04 y simulation for model basin 1. One calibration was completed by a consultant by using the current SWMM5 Groundwater module and the current SWMM5 Runoff component before the SPU-LTCRDI method was developed. The consultant used an automated calibration process that utilized 500 Monte Carlo simulations for initial calibration. Result from the automated calibration was further refined by the consultant through manual calibration. The other calibration was a manual calibration completed by the author using the SPU-LTCRDI method together with the current SWMM5 Runoff component.

Table 1 Characteristics of model basin 1.

Characteristics of Model	Value
Area, acre (ha)	221.4 (89.6)
Number of subcatchments	365
Number of subcatchments with aquifer	126
Number of aquifer	126
Number of nodes	140
Total length of conduits, miles (km)	5.9 (9.5)
Average % impervious	40
Total duration of simulation (y)	3.04
Total duration of observed data (y)	1.51

Figure 9 Plan view of model basin 1.

Table 2 Characteristics of flow monitoring data.

Events	Observed Peak Flow mgd (L/s)	Observed Flow Volume MGal (ML)	Duration of Event (d)
Nov 2008	2.94 (128.81)	6.05 (22.90)	26
Jan 2009	5.06 (221.69)	18.01 (68.18)	31
Nov 2009	3.46 (151.59)	15.51 (58.71)	19
Jan 2010	3.20 (140.20)	19.98 (75.63)	22
Over entire duration of observed data	6.31 (276.46)	166.82 (631.48)	552.73

Figure 10 Model basin 1 GOF peak flow.

Figure 11 Model basin 1 GOF flow volume.

Figure 12 Model basin 1 standardized bias.

Figure 13 Model basin 1 SRMSE.

Figure 14 Model basin 1 correlation coefficient.

When compared with results from the current SWMM5 Groundwater module, the GOF measure results indicate that the SPU-LTCRDI method is able to enhance the agreement between simulation results and flow monitoring data over the selected wet weather events as well as the whole 1.51 y of observed data in a continuous long term simulation of 3.04 y.

2.2 Model Basin 2

Table 3 summarizes the characteristics of model basin 2 and the total duration of the simulation. Figure 15 shows a plan view of the basin model. Table 4 summarizes the characteristics of the overall observed data and wet weather events used for comparison.

Figures 16 to 20 summarize the aforementioned five statistics of goodness of fit measures for two independent calibrations of the 2.45 y simulation for model basin 2. One calibration was completed by a consultant using the current SWMM5 Groundwater module and the current SWMM5 Runoff component before the SPU-LTCRDI method was developed. The consultant used an automated calibration process that utilized 500 Monte Carlo simulations for initial calibration. Result from the automated calibration was further refined by the consultant through manual calibration. The other calibration was a manual calibration completed by the author using the SPU-LTCRDI method together with the current SWMM5 Runoff component.

Table 3 Characteristics of model basin 2.

Characteristics of Model	Value
Area, acre (ha)	76.3 (30.9)
Number of subcatchments	146
Number of subcatchments with aquifer	143
Number of aquifer	50
Number of nodes	63
Total length of conduits, miles (km)	2.5 (4.0)
Average % impervious	48
Total duration of simulation (y)	2.45
Total duration of observed data (y)	1.41

Figure 15 Plan view of model basin 2.

Table 4 Characteristics of flow monitoring data.

Events	Observed Peak Flow mgd (L/s)	Observed Flow Volume MGal (ML)	Duration of Event (d)
Nov 2008	2.83 (123.99)	6.19 (23.43)	26
Jan 2009	2.93 (128.37)	8.76 (33.16)	31
Nov 2009	3.49 (152.91)	6.32 (23.92)	19
Jan 2010	2.27 (99.45)	8.59 (32.52)	22
Over entire duration of observed data	3.67 (160.79)	78.59 (297.49)	515.19

Figure 16 Model basin 2 GOF peak flow.

Figure 17 Model basin 2 GOF flow volume.

Figure 18 Model basin 2 standardized bias.

Figure 19 Model basin 2 SRMSE.

Figure 20 Model basin 2 correlation coefficient.

When compared with results from the current SWMM5 Groundwater module, the GOF measure results indicate that the SPU-LTCRDI method is able to enhance the agreement between simulation results and flow monitoring data over the selected wet weather events as well as the whole 1.41 y of observed data in a continuous long term simulation of 2.45 y.

3 Conclusions and Next Steps

Together with the current Surface Runoff component of the SWMM5 engine for modeling RD Inflow, the ability of the SPU-LTCRDI method to model RDI adds to SWMM5’s current capability in modeling long term continuous RDII that does not display seasonal repeatable patterns and where the infiltration portion of RDII enters into an underground conveyance system through irregular defects along the system. The SPU-LTCRDI method uses a single reservoir configuration as in the current SWMM5 Groundwater module. In terms of groundwater discharge, it also uses the same six parameters as in the current SWMM5 groundwater discharge equation algorithm. Both of these measures enhance SPU-LTCRDI’s compatibility with the rest of the SWMM5 modules, minimize the code modification needed for its implementation into the SWMM5 code, and ensure that the effort needed to use this method would be approximately the same as using the current SWMM5 Groundwater module.

However, despite of their similarities, there are a number of differences between the two modules. The SPU-LTCRDI method has three groundwater discharge outlets rather than two from an aquifer. It uses different definitions for the moisture content parameters, different usage of the aquifer Bottom Elevation parameter, and has a different ower groundwater loss rate equation algorithm. A SPU-LTCRDI aquifer cannot be shared among RDI generating subcatchments that do not have the same surface elevation, aquifer bottom elevation and discharge location. For the groundwater discharge component, a new groundwater discharge equation algorithm has been developed. It uses the same six parameters as the current SWMM5 Groundwater module, but each parameter is being used differently in the new equation algorithm.

The SPU-LTCRDI method has been tested with two different basin models in a 3.04 y and a 2.45 y continuous simulation. Simulation results were compared with 1.51 y and 1.41 y continuous flow monitoring data respectively. Five different goodness of fit measures were used to evaluate the quality of the agreement between simulation results and flow monitoring data. Results showed that the SPU-LTCRDI method was able to enhance the agreement between simulation results and flow monitoring data for both simulations when compared with the agreement between simulation results generated by the current SWMM5 Groundwater module.

Currently, unlike the current SWMM5 Groundwater module, the SPU-LTCRDI method does not model exfiltration from an underground conveyance system to its adjacent soil. When water level inside a conveyance system is higher than its adjacent groundwater level, the SPU-LTCRDI method stops infiltration from the adjacent soil into the underground conveyance system but does not generate flow from the conveyance system to its adjacent soil. This feature is not being currently added into the SPU-LTCRDI method because field exfiltration data is not available for use in development and calibration. This would become the next phase of development of the SPU-LTCRDI method once field exfiltration data becomes available.

Notice

The algorithm described in this chapter is the property of the City of Seattle. With the exception of USEPA implementing it into the official USEPA SWMM5 engine, the described algorithm and any part of it cannot be implemented into any software without the written approval of the City of Seattle. However, interest in testing the above algorithm is very welcome. Details on how to obtain a copy of the SWMM5.dll file with the above algorithm implemented and associated documentation for testing purposes can be obtained from Joseph Pang at joseph.pang@seattle.gov.

Acknowledgment

Special thanks to Lewis Rossman for providing comments, critiques, review of codes, and instruction on how to recompile code during the development and testing processes.

References

Cornell University. 2013. Why Preferential Flow is Important. Ithaca, NY: Cornell University Department of Biological & Environmental Engineering. http://soilandwater.bee.cornell.edu/research/pfweb/educators/intro/why.htm.
MGS Engineering Consultants, Inc. 2010. Methods for Estimating Control Volumes for CSO Reduction: Technical Guidance Manual. Seattle, WA: Seattle Public Utilities.
Rimmer, A. and A. Hartmann. 2012. “Simplified Conceptual Structures and Analytical Solutions for Groundwater Discharge Using Reservoir Equations.” Chap. 10 in Water Resources Management and Modeling, edited by P. Nayak. Rijeka, Croatia: Intech. ISBN: 978-953-51-0246-5. http://www.intechopen.com.
Rossman, L. A. 2010. Storm Water Management Model User’s Manual Version 5.0. Washington, DC: U. S. Environmental Protection Agency.

Identification
CHI ref #:	C369 587
Volume:	22
DOI:	https://doi.org/10.14796/JWMM.C369
Cite as:	JWMM 22: C369

Publication History
Received:	May 10, 2013
First decision:	N/A
Accepted:	December 01, 2013
Published:	February 15, 2014

Status
# reviewers:	2
Version:	Final published

Copyright
© 2014 CHI. Some rights reserved.

License
This work is licensed under a Creative Commons Attribution 4.0 International License.

JWMM Login

Proposed Addition to SWMM5 Capability to Model Long Term Continuous Rainfall Dependent Infiltration

Abstract

1 Theory

1.1 The aquifer Configuration

The Moisture Content Parameters

Porosity, Field Capacity, Wilting Point

Unsaturated Zone Moisture

Lower GW Loss Rate

Aquifer Bottom Elevation

Effects of Aquifer Parameters on Groundwater Discharge

1.2 The Groundwater Flow Component

Current SWMM5 Groundwater Discharge Equation Algorithm

SPU-LTCRDI Groundwater Discharge Equation Algorithm

Preferential Flow and Matrix–Interflow Infiltration (qPMI)

Baseflow Infiltration (qB)

Final Form of the Groundwater Discharge Equation Algorithm

2 Sample Results

2.1 Model Basin 1

2.2 Model Basin 2

3 Conclusions and Next Steps

Notice

Acknowledgment

References

Image

PAPER INFO

Identification

Publication History

Status

Copyright

License

AUTHORS

Joseph Pang

ADDITIONAL DATA

COMMENTS

RELATED PAPERS

TAGS

Preferential Flow and Matrix–Interflow Infiltration (q_PMI)

Baseflow Infiltration (q_B)