Flood Routing with the Soil and Water Assessment Tool : A Review and Verification

The Soil and Water Assessment Tool (SWAT) is one of the most widely used hydrologic models. SWAT has been undergoing constant changes since its development. However, compartment review and testing of SWAT are comparably limited, especially the flood routing functions. In this study, the daily flood routing subroutines of different SWAT versions were reviewed and tested using a well observed segment of the Weser river located in Germany. Results show several problems with the flood routing subroutines of SWAT. The variable storage subroutine of SWAT (revision 664) does not transform the flood wave. Unphysical results could be obtained with the variable storage routing of SWAT (revision 528). The Muskingum subroutine of SWAT (revisions 664 and 528) overestimates channel evaporation (resulting in a bias of 14% to 19% in streamflow) and underestimates transmission losses. Simulated results show that the timing and shape of the flood wave could be improved with a corrected Muskingum subroutine. Based on the results of this study, we suggest the SWAT user community to review their existing SWAT models to see how the aforementioned issues will affect their methods, findings and conclusions.


Introduction
Flood routing is predicting the timing and magnitude of a flood wave at a downstream point of a reach from the known data at an upstream point [1].It plays an important role in flood forecasting, reservoir design, and flood control.It is also one of the key components of hydrologic models and serves as a basic function for sediment and nutrient routing [2].In addition, hydrologic models are often calibrated based on streamflow, sediment, and/or nutrient yields.Therefore, having a robust and well-tested flood routing function in hydrologic models is important.
The Soil and Water Assessment Tools (SWAT) is a conceptual semi-distributed hydrological model used to predict the effect of land use management practices on water, sediment, and nutrient yield [2,3].SWAT has been undergoing constant changes since its development [4].SWAT have been tested and applied worldwide [4][5][6].However, compartment testing and review of SWAT are limited, especially regarding its flood routing functions.
In SWAT, users can select either the variable storage method [7] or the Muskingum method [8,9] for flood routing [3].There have been some modifications to the flood routing functions in SWAT [e.g., 10,11].Kim and Lee [10] suggested to implement a nonlinear storage equation obtained by coupling continuity and Manning's equations, to avoid underestimation of peak flows and false signals during the recession period with the Muskingum used in SWAT.Pati et al. [11] replaced the Muskingum used in SWAT with the variable parameter McCarthy-Muskingum to enhance channel routing.In the official SWAT revisions (https://swat.tamu.edu/)however, the model uses the two aforementioned methods [7][8][9] for flood routing.Despite that, our preliminary studies show that there is a significant difference in terms of simulated streamflow between the two methods and between the same method in different SWAT versions (as shown in section 4).
The main objective of this study is to explain the differences in simulated streamflow and motivate the users of the model to draw more attention on the flood routing processes.This study provides (1) an overview into the technical aspects of the two flood routing methods used in SWAT and (2) an insight to the code implementation of these concepts in different SWAT revisions, and (3) an improvement of the SWAT flood routing functions.A case study for the Weser river segment, Germany is used as a verification example.

Variable Storage Method
The variable storage method is based on the continuity equation [7]: where I (m 3 /s) and O (m 3 /s) are the inflow and outflow rate for a river reach, respectively, t (s) is the time, S (m 3 ) is the storage.In discrete form, Eq. 1 becomes: where the subscripts "1" and "2" refer to the start and end of the routing time interval ∆t (s), respectively.Eq. 2 can be rearranged to have the following form: where I a = 0.5(I 1 + I 2 ) is the average inflow rate during the time interval.The travel time, T (s), is calculated as follows: Eq. 3 can be rewritten using Eq. 4 to obtain the relation between the storage coefficient and the travel time: where C is the storage coefficient: The condition C ≤ 1 must be satisfied to avoid unphysical results.

Muskingum Method
The Muskingum method is based on the continuity equation (Eq. 1 and 2) and the empirical linear storage equation [9,12]: where S (m 3 ) is the total storage in channel, K (s) is the storage constant, X (-) is a weighting factor, ranging from 0 to 0.5, I (m 3 /s) and O (m 3 /s) are inflow and outflow rate.From Eqs. 2 and 9, the following relation can be obtained: where: where C 1 , C 2 , C 3 are coefficients.To avoid numerical instability, the following condition must be satisfied:

Flood Routing with SWAT
In this section, the daily variable storage and Muskingum method in SWAT2000, SWAT2005, SWAT2009 (revision 528) , and SWAT2012 (revision 664) were reviewed.The source code of these SWAT versions was downloaded from the official SWAT website (https://swat.tamu.edu/).The daily variable storage and Muskingum routing subroutines in these SWAT versions are named "rtday" and "rtmusk", respectively.

Variable Storage Routing Subroutine
The flood routing in SWAT was originally performed only with the variable storage method [2,13].Although the four aforementioned SWAT versions were reported to incorporate the variable storage routing method as described in previous sections [3,[14][15][16], our study showed that the rtday subroutine of SWAT2012 did not perform a transformation of the outflow compared to the inflow.In this subroutine, the outflow from a reach was calculated as follows: where rtwtr (m 3 ) is the outflow volume during a day, rcharea (m 2 ) is the cross-sectional area of flow, 86400 (s) is the number of seconds in a day, vc (m/s) is the average flow velocity, calculated as using the equation: where sdti (m 3 /s) is the average flow rate, calculated as follows: where vol (m 3 ) is the volume of water in reach at the beginning of a day and the inflow volume.Eq. 18 shows that the rtday subroutine of SWAT2012 does not perform flood routing.The outflow volume is simply assigned as the volume of water in reach at the beginning of a day, i.e. the inflow volume.
In the rtday subroutine of other SWAT versions (SWAT2000, SWAT2005, SWAT2009), flood routing is performed using Eq.7 .In these SWAT versions, C is assigned as 1, if C > 1 to avoid unphysical results.It means that the outflow volume is assigned the total inflow volume and storage volume.

Muskingum Routing Subroutine
The Muskingum routing method in SWAT was reported to follow the concept described in the previous section [3,[14][15][16].However, revision of the "rtmusk" subroutine in four SWAT versions shows that SWAT2000 and SWAT2005 do not check the numerical stability condition (Eq.14).Within SWAT2009 and SWAT2012, if the numerical stability condition is not satisfied, the daily time step is divided into smaller time steps, 12 hours or 6 hours or 1 hour.However, it does not always guarantee the numerical stability condition.In addition, (1) the calculated evaporation in a reach for each sub-daily time step in SWAT2009 and SWAT2012 was taken as daily evaporation and (2) the calculated transmission losses were not summed up during each time step to have the total amount of transmission losses during a day.These result in an overestimation of reach evaporation and underestimation of transmission losses, which, in combination, can effect the calibration of SWAT models in different ways depending on the characteristics of the case study area.In this study, the rtmusk routing subroutine of SWAT2012 was modified to overcome the aforementioned problems.The source code of this modified subroutine is available at https://github.com/tamnva/updated_SWAT2012.The storage constant K and weighting factor X in four aforementioned SWAT versions are user defined.

Study Area and Data
The study area is a part of the Weser river located in Lower Saxony and North Rhine-Westphalia, Germany (Fig. 1).This Weser river segment has a total length of 142 km, an average slope of 0.01 % and an average width of 108 m.Daily inflow data at the Karlshafen gauging station and outflow data at the Vlotho gauging station were obtained from the the Wasser-und Schifffahrtsdirektion Mitte.Daily streamflow data at three tributary channels were taken from the Niedersächsische Landesbetrieb für Wasserwirtschaft, Küsten-und Naturschutz.The drainage areas corresponding to the Karlshafen and Vlotho gauging stations are 14790 km 2 and 17620 km 2 , respectively.Weather data in this area were taken from the Deutscher Wetterdienst.The Digital Elevation Model (DEM) obtained from the Niedersächsische Landesbetrieb für Wasserwirtschaft, Küsten-und Naturschutz (NLWKN) shows that the elevation of the study area varies between 42 m to 522 m above mean sea level (a.m.s.l).Land use/Land cover and soil maps were obtained from the CORINE Land Cover project and the Bundesanstalt für Geowissenschaften und Rohstoffe (BGR), respectively.

Simulation Scenarios
Due to the similarities (1) in the variable storage subroutine between SWAT2000, SWAT2005, and SWAT2009 and (2) in the Muskingum routing subroutine between SWAT2009 and SWAT2012, therefore, simulations in this study were only performed with 4 different routing subroutines (Tab.1).The numerical instability due to neglecting the condition in Eq. 14 in the Muskingum routing subroutine of SWAT2000 (or SWAT2005) was discussed in detail by Kim and Le [10].In order to compare the effects of different routing subroutines on simulated streamflow, different flood routing subroutines (Tab. 1) were called by SWAT2012.To minimize the error of simulated flow at the outlet of the study area, (1) observed flow at the Uchdorf, Afferde, and Welsede gauging stations was used as model inflow condition, and (2) contributing flow from other ungauged tributary channels and intermediate subbasins was simulated by using parameters obtained when calibrating for streamflow at the Welsede gauging station.The Manning's value was selected as 0.025 for all channels in all simulations because it gives the best results in case of simulation with the Muskingum method.In case of simulations with the Variable storage method, the Manning's value has very minor effect on simulated streamflow.

Results and Discussions
Simulated streamflow at the Vlotho gauging station shows that V2009 and V2012 failed to perform flood routing (Fig. 2a and c) despite having good Nash-Sutcliffe efficiency [NSE,17] and percentage bias (PBIAS) due to the high influence of inflow taken from observed values (Tab.2).For example, (1) the timing of simulated peaks is earlier than the observed and (2) the simulated peaks are higher than the observed.The simulated outflow from V2009 and V2012 is almost identical to inflow.The variable storage method does not show a transformation of the flood wave as it moves downstream.V2012 simply assigns inflow as outflow as previously discussed.The similarities between simulated outflow from V2009 and V2012 are due to the fact that V2009 assigns C = 1 most of the time (to avoid C > 1).We found that the unphysical oscillation of simulated flow with V2009 (as shown in the red circle in Fig. 2a and c) occurs when flow in the channel flood plain occurs (or disappears).This is due to the assumption about the flood plain geometry and the use of the Manning's equation.The flood plain width in SWAT is assumed five times the top width of the channel.When flow starts to change from non flood-plain flow to flood-plain flow (or from flood-plain flow to non flood-plain flow), there is a sudden increase (or decrease) in the flow velocity, the travel time T, the storage coefficient C, and ultimately the simulated outflow.We also found that, in case C > 1, changing the simulation time step during a simulation in order to have C ≤ 1 could also cause the same problem.The simulation time step should be the same for the entire simulation time.A solution to this problem could be selecting a sufficiently small time step (small ∆t) and a sufficiently long river section (high T), however, this is event-based and reach specific.Therefore, the V2009 and V2012 are not practical for long-term simulation at basin-scale.Fig. 2b and d shows simulated outflow at the Vlotho gauging station using M2012 and M2012*.It is seen that simulated outflow from M2012 and M2012* is better than that from V2009 and V2012 (Fig. 2a and c and Tab. 2).The timing of simulated peaks and the shape of the flood wave matched well with observed outflow.The effect of the overestimation of channel evaporation with M2012 is negligible in this case (Fig. 2b and d) because the flow rate is high compared to the evaporation rate.The underestimation of flow from M2012 and M2012* as shown by the PBIAS (Tab.2) is also due to the underestimation of flow from tributary reach and intermediate subbasins.
The effect of overestimation of channel evaporation with M2012 on the streamflow is more pronounced for head water catchments.For example, Fig. 3 shows the simulated streamflow at Wesede using M2012 and M2012*.It is clearly seen that simulated streamflow using M2012 is lower than that from M2012*.When streamflow is smaller than a certain thereshold (<3 m 3 /s), it could be completely lost (Fig. 3).For a long-term simulation (from 1990 to 2000), the results show that about 18% of streamflow is lost due to overestimation of channel evaporation.Similar results were also found for simulated flow at the Uchtdorf and Afferde gauging stations.About 19% and 14% of streamflow at the Uchtdorf and Afferde gauging stations, respectively, is lost due to the overestimation of channel evaporation.Furthermore, in semi-arid areas, the effect of overestimation of channel evaporation with M2012 is expected to be higher compared to humid regions as shown here.

Conclusions and Recommendations
In this study, we reviewed and tested the flood routing subroutines of different SWAT revisions.Results show that there are major issues in both flood routing subroutines of SWAT, the variable storage and Muskingum routing subroutines.In case of simulation with the variable storage method, (1) SWAT2012 (revision 664) does not perform flood routing, and (2) unphysical oscillation of simulated flow was observed with SWAT2009 due to the assumption of the flood plain geometry and the use of the Manning's equation.In case of simulation with the Muskingum subroutine (SWAT2009 revision 528 and SWAT2012 revision 664), (1) about 14% to 19% of streamflow volume could be lost due to overestimation of channel evaporation , (2) unphysical results could be obtained due to violating the numerical stability condition, and (3) channel transmission loss is underestimated if it is activated.A corrected Muskingum subroutine from SWAT2012 (revision 664) shows a good estimation of the flood wave.
Based on the results in this study, we suggest the SWAT user community to check and upgrade their existing SWAT models, which used the aforementioned affected SWAT versions to see how our improvement will affect their methods, findings and conclusions.This is especially important for studies focused on losses are significant.A bias contribution of streamflow in an order of more than 14% is considered significant for the calibration of other components of the model.Although SWAT+ is announced as the next generation of SWAT, there are still many SWAT users working with existing SWAT models.Therefore, review and testing other functions of SWAT are also suggested.

Figure 1 .
Figure 1.Location of the Weser river segment and its tributary rivers.

Figure 2 .
Figure 2. Simulated outflow at the Vlotho gauging station for two flood events using (a, c) V2009 and V2012 and (b, d) M2012 and M2012*.

Table 1 .
List of the flood routing subroutines used for simulation.

Table 2 .
Model performance statistics for two flood events (Fig.2).