Automatic Calibration of the U.S. EPA SWMM Model for a Large Urban Catchment

Major metropolitan areas are characterized by continuous increases in imperviousness due to urban development. Increasing imperviousness increases runoff volume, and maximum rates of runoff, with possible negative consequences for natural systems. To avoid environmental degradation, new development standards often prohibit increases in total runoff volume and may limit maximum flow rates. Methods to reduce runoff volume and maximum runoff rate are required, and solutions to the problems may benefit from the use of advanced models. Developing management strategies for large urban watersheds is especially difficult when the flow network is large or complex, such as the Ballona Creek Watershed located in west Los Angeles, and is the subject of the present study. Designers and managers are now using mathematical models to address the problems of such a large and complex watershed.

Various models are available to manage urban runoff, including HEC-1 (U.S. Army Corps of Engineers 1985), TR-20 and TR-55 (Soil Conservation Service 1983, 1986), MOUSE (Danish Hydraulic Institute 1995), HydroWorks (HR Wallingford Ltd. 1997), and storm water management model (SWMM) (Huber and Dickinson 1988). The Environmental Protection Agency’s (EPA) SWMM is a dynamic rainfall-runoff model for simulation of quantity and quality problems associated with runoff from urban areas (Huber and Dickinson 1988). The model simulates aspects of urban hydrologic and water quality cycles, including rainfall, snow melt, surface runoff, transport through the drainage network, storage and treatment, and receiving water effects. Both single event and continuous simulation can be performed on catchments having storm drains, or combined sewers and natural drainage.

In this study SWMM was modified and adapted for the upper Ballona Creek Watershed. The drainage system and subcatchments were defined using records from the City and County of Los Angeles, and a ground survey of the undefined areas. Imperviousness was obtained from land-use data and the spatial distribution of precipitation was developed using the geographic information system (GIS) and isohyetal map. Calibration was performed using a constrained optimization procedure, and a total of ten storms were used for calibration and validation. The calibrated model was able to predict the observed outputs with reasonable accuracy.

SWMM, developed for the EPA, has been applied on numerous watersheds in United States cities and other part of the world (Selvalingam et al. 1987; Warwick and Tadepalli 1991; Bhaduri et al. 2001). It has been applied to all types of storm water management from urban drainage (Zaghloul 1998; Campbell and Sullivan 2002) to flood routing (Hsu et al. 2000). An alternative approach to commercial packages with graphical user interfaces is to couple SWMM or a suitable model to a GIS. Martin et al. (2005) presented a state-of-the-art critical review of current trends in interfacing GIS with predictive water resources models.

Physically based deterministic rainfall-runoff models require initial estimates of one or more parameters. Manual calibration is labor intensive, especially when the catchment is large and complex. Automatic parameter estimation and calibration methods have been implemented to overcome this difficulty. Maalel and Huber (1984) proposed a calibration procedure for SWMM using multiple, single events, and continuous simulation. Baffaut and Delleur (1989) described an automated parameter estimation and calibration procedure that used expert system technology. Liong et al. (1991) described a knowledge-based system for automating the calibration of the SWMM’s runoff block. Methodologies for estimating global optimum sets of calibration parameters have been presented by Ibrahim and Liong (1992, 1993), and Liong and Ibrahim (1994) and Liong et al. (1995). The proposed methods consist of two operations. The first derives functional relationships using the response surface method and the second estimates the optimal set of parameters using a probabilistic approach. Liong et al. (1995) and Balascio et al. (1998) have described genetic algorithms for calibrating SWMM. Zaghloul and Abu Kiefa (2001) used an artificial neural network for sensitivity calibration of the SWMM model.

Box (1965) developed a pattern-search optimization method called the complex method, which has been applied to many problems. The method uses a set of parameter values, called vertexes, which are successively created, evaluated, and eliminated to locate the optimum solution. Advantages of the complex method are that it can be applied to any cost function, does not require the calculation of any derivatives, both explicit constraints (upper and lower bounds of the parameters) and implicit constraints (limits on properties of the resulting solution) can be used, and the number of vertexes can be increased to better manage problems with local optima. The complex method has been used by many researchers for different problems and with stable calibration results. Yuan et al. (1993) applied this method on the high-purity oxygen-activated sludge process and found the algorithm to be very robust, converging to the same optimal values for a wide range of starting points and constrains, and concluding that if the model has only random errors (no systematic errors), the algorithm is capable of estimating the parameters to any acceptable precision. Haque (1996) used the method for structural optimization, and Subramanian et al. (2005) used it for warpage analysis for an optimal housing of compact disks. The complex method was used in this problem because the Ballona Creek Watershed is far too complicated for manual calibration, and anticipated future uses of the calibrated model need to incorporate implicit constraints.

Another technical problem in developing an SWMM model for a large watershed is defining precipitation data. Rainfall data drive the model and produce runoff, which means that an accurate estimation of rainfall data determines the success of the modeling effort. Singh and Chowdhury (1986) reviewed 13 different methods for computing mean, areal rainfall, including the inverse distance-squared method (IDSM). The result of their investigation showed that all 13 methods yielded comparable estimates.

In its standard form, SWMM version 4.3 has ten processes and is able to simulate 500 subcatchments and 500 channels/pipes. For the current study, seven processes (or blocks) were removed; Process 2—Statistics, Process 3—Graph, Process 4—Combine, Process 5—Rain, Process 6—Temperature, Process 9—Extran, Process 10—Storage/Treatment. A new process was added which uses an optimization technique for calibration. The executive block (Process 1) was replaced by a new main block to accommodate the addition of the calibration process. The runoff and transport processes were used for routing in lieu of Extran, and were modified to handle the larger watershed, with its greater number of subcatchments and channels/pipes. The runoff process was also modified to incorporate the rainfall weighting factor for each subcatchment, which was based upon an isohyetal map, described later. The modified SWMM code had three processes and the maximum number of subcatchments and channels/pipes was increased to 1,580 and 4,300, respectively. The modified version GIS/SWMM is illustrated in Fig. 1.

The GIS used in this study was Environmental System Research Institute ARC/INFO (ESRI, Redlands, Calif.) and Compaq Visual FORTRAN 6.5 was used to compile the SWMM code.

Ballona Creek, which has become a concrete lined channel to facilitate flood routing, drains the entire watershed and enters Santa Monica Bay at Marina del Rey. There is a USGS gauging station and a Los Angeles County Department of Public Works (Alhambra, Calif.) monitoring station near the intersection of the Creek and the IS 405 Freeway. The area upstream the gauging station, called the upper Ballona Creek Watershed, was modeled in this study because of the availability of monitoring data. The total contributing area of the catchment above the gauging station is to $217{\mathrm{km}}^2$ , and 65% of the total area is residential land use and over 19% is either commercial or public/industrial land use. A total of 1,579 subcatchments are delineated within the watershed. Rainfall in each subcatchment becomes overland surface runoff and flows are assumed to exit each subwatershed through a single outlet (i.e., catch basin). Surface runoff flow entering the inlet is then routed through the connected channel or pipe (open or closed, depending on the location). Finally, all surface flow is summed up at the watershed outlet after routing through 2,648 channels and pipes, and the total length of all channels and pipes is over $283\mathrm{km}$ . The slope of the channels/pipes averages 2.6% and ranges from 0.005 and 47%. Fig. 2 shows the watershed with subcatchments, channels/pipes, and the flow gauge at the watershed outlet. Land-use data for the watershed were obtained from the Southern California Association of Governments (SCAG) Database, which contains land-use records based on Anderson level II (Anderson et al. 1976) that were reduced to eight types based upon environmentally similar characteristics (Park and Stenstrom 2006).

The precipitation data were obtained from the 1997 Edition of EarthInfo, Inc. (Boulder, Colo. 80301-2846) NCDC hourly precipitation CDROM to obtain the average rainfall for each storm. Ten rain gauges were available. Precipitation in each subcatchment was calculated from the average rainfall using weighting factors from an isohyetal map based upon 50-year record rainfall events (Stenstrom and Strecker 1993). The GIS was used to create the rainfall weighting factor corresponding to the rainfall in each subcatchment, which was read by the SWMM runoff process. Hourly flow data were obtained from the Los Angeles County Department of Public Works (Alhambra, Calif.). Ten representative storms from the 1994 to 1996 wet seasons were selected and Table 1 shows the dates, rainfall, and duration.

Each subcatchment and all channel/pipe in Ballona Watershed were manually digitized and converted into GIS databases based on the information provided by the drainage maps from the Los Angeles County Department of Public Works and the Bureau of Engineering, City of Los Angeles (LADWP) (Los Angeles). Each subcatchment was delineated based on the location of inlets (catch basins) and the flow direction of each street provided by the drainage maps. Fig. 2 shows the imperviousness which are based on the ground surveys by LADWP. The Horton (1942) equation was used for estimating infiltration in the pervious areas all soil was assumed to be highly permeable (ASCE 1996). The USGS $7.5\min$ digital elevation model (DEM) spatial data transfer system (SDTS) combined with the GRID module of ARCINFO GIS were used to compute the subcatchment and channels/pipes slopes. The slopes were then used to compute the impervious depression storage coefficients based on methods by Kidd (1978) and Viessman et al. (1989). These parameter values were used as the initial estimates or starting values for the optimization, which is described later. More information on the methodology is available (Wong et al. 1997).

The complex method was used for calibration and was incorporated into the modified GIS/SWMM as illustrated in Fig. 1. Calibration begins by generating at least $m+1$ sets of parameters, called vertexes, where $m$ =number of parameters being identified. In this case the vertexes were calculated from the upper and lower constraints using random numbers as followswhere $X_{i,j}$ =value of the $i$ th parameter in the $j$ th set of parameters; $U_i$ =upper bound of the $i$ th parameter; $L_i$ =lower bound of the $i$ th parameter; $r_{ij}$ =random real number ranging between 0 and 1; $i$ is $i$ th parameter; and $j$ = $j$ th vertex. It is also possible to manually select one or more of the initial vertexes, and the first vertex was manually entered in this study. The parameters for the first vertex were based on their most probable values. It is also possible to test for local optima by selecting different initial vertexes, different constraint intervals, or a different series of random numbers used in calculating the parameters for the remaining vertexes. The objective function, $F$ , for each vertex is evaluated and the $X_k$ vertex having the greatest function value $F$ is selected (rejected). The centroid of the remaining points (in the case of Cartesian coordinates, the average value of the coordinates of the remaining vertexes), is then computed. The set $X_k$ is projected toward the centroid using a projection factor $\gamma$ (Box recommended the use of $\gamma$ = 1.3). This creates a new vertex $X_g$ that replaces the rejected vertex if and only if the error of the new vertex is less than the rejected vertex. Conversely, if the new set has an error greater than or equal to the error of the rejected set, the projection distance, $\gamma$ , is reduced by 50%, and another vertex is projected. The projection factor is decreased until a new vertex with a lower error than the rejected vertex is obtained. The technique is repeated by selecting the next vertex with the largest error for replacement. The whole process continues until the termination criterion assigned by the user is satisfied. The termination criteria can be the maximum number of iterations, the minimum size of $\gamma$ , or a relative improvement criterion of $F$ . Constraints are handled by reducing the value of $\gamma$ . If a vertex is projected outside of the feasible region, defined by the upper and lower bounds on the parameters, the value of $\gamma$ is reduced until a feasible vertex is obtained. A violation of an implicit constraint is handled in the same way.

The procedure is not guaranteed to obtain a global optimum and can fail if the initial set of vertexes surrounds a nonconvex area, but is robust and rarely failed to find the global optimum, even when repeating an optimization many times using different initial vertexes, and a wide range of explicit constraints. Increasing the number of vertexes relative to the number of parameters generally improves the likelihood of obtaining a global optimum, but requires more iterations. To ensure that a global optimum is obtained, the problem can be solved multiple times with different starting points. The most likely difficulty a user might have in implementing this method for complicated problems with many constraints is the creation of a suitable, initial set of vertexes.

The maximum number of iterations depends upon the number of parameters and the complexity of the problem. The average number of iterations for this problem was 60 and rarely was the required number of iterations greater than 200 iterations. Yuan et al. (1993), when simulating the activated sludge process using 3–7 parameters, found that 2,000 iterations were always sufficient. The large number of required iterations means that the complex method will be useful for problems which have easily evaluated objective functions and implicit or explicit constraints. It is not necessary to calculate derivatives of the model with respect to the parameters, which can sometimes be difficult. For problems with more computationally intensive objective functions, other approaches may be justified (Becker and Yeh 1972; Tzeng et al. 2003).

The objective function used in this study iswhere $Q$ =total flow volume; $P$ =peak flow rate; $f$ =instantaneous flow rate at the outlet hydrograph; the superscript ${}^{*}$ denotes the $i$ th predicted value, the subscript $i$ denotes the $i$ th observation of the outlet hydrograph; $w_1$ , $w_2$ , and $w_3$ =weights; and $n$ =total number of observed data points in the outlet hydrograph. The division by the predicted value normalizes the objective function and the three terms have equal weights for equal values of $w_i$ . Other objective functions and methods of weighting are possible, such as a function to locate the peak flow as a function of time. For the initial calibrations, weights of 1.0 were used for $w_1$ and $w_2$ and zero for $w_3$ .

Four parameters were chosen for calibration: subcatchment imperviousness, width, impervious depression storage coefficient, and channel Manning’s roughness coefficient. The parameters were varied uniformly across the watershed and the constraints on the four parameters were selected based on the physical meaning of the parameters (i.e., physically plausible values, see Table 2), and trial runs of the calibration algorithm. Only explicit constraints were needed in this problem, and were generally set to 50 and 150% of most probable parameter values. Reducing the range of constraints will increase the speed of convergence, but may exclude appropriate parameter values, if the initial vertex is poorly estimated.

The US EPA SWMM model and a raster-based GIS were applied to a large urban catchment in west Los Angeles. The GIS was used to manage and process input and output data from the model. The model included more than 2,500 nodes. The complex method of Box was used to calibrate the model for ten storms using several strategies and weighting functions. The following specific conclusions were made: