Artificial viscosity technique for direct runoff calculation

. In this study, an Artificial Viscosity (AV) technique was developed for direct runoff calculation, instead of using the conventional Synthetic Unit Hydrograph (SUH) methods. We solved the Shallow Water Equations (SWE) with second-order accurate Godunov finite volume model and fourth-order Runge-Kutta explicit scheme. The AV technique was devised with a Laplacian and a biharmonic operator, and employed to solve the convective terms of the SWE. This technique was applied to rainfall-runoff laboratory cases with the measured rainfall and observed direct runoff values. For comparison purpose, several SUH methods commonly used such as SCS, Snyder, GAMA-1, ITB-1, ITB-2, and Nakayasu, were also used to compute the direct runoff values. The results showed that the AV technique could predict three parameters (i.e., peak discharge, time-to-peak, and shape of the direct runoff hydrograph) accurately. Meanwhile, significant discrepancies were shown by the SUH methods in estimating such parameters, thus indicating that the SWE modeling with an AV technique is significantly more accurate than the SUH methods in predicting the direct runoff hydrograph. This study shows an interesting example of how modern numerical computations can improve the runoff prediction and may be included as a standard technique for runoff computation in the future.


Introduction
Hydrology computation is required for flood discharge prediction.In Indonesia, especially for ungauged basins, hydrology computation should be carried out according to SNI 2415:2016 [1], in which three synthetic unit hydrograph (SUH) methods are defined as standard.First, it is the Soil Conservation Service (SCS) method, which was originally developed in 1954 by the United States Department of Agriculture (USDA) for agricultural areas in the United States [2][3].This method was devised based on the curve number (CN) empirical equation in relation with a comparison between the discharge value and the peak discharge as well as a comparison between the time value and the time-to-peak that can be estimated using an empirical formula such as the Kirpich method [4].
The second method is the Snyder method, which was developed in 1938 based on three parameters, namely effective rainfall duration, peak discharge, and lag time [5][6].The third one is the Gadjah Mada 1 (GAMA-1) method, which was developed in 1993 in Indonesia by Prof. Sri Harto [7][8].The time-to-peak of this method is computed based on a comparison of the sum of the river/channel length in a basin, while the peak discharge is calculated based on an empirical factor of the number of river confluence.
In addition to these three methods, several SUHs are also mentioned in SNI 2415:2016 [1] such as Institut Teknologi Bandung (ITB), Snyder-Alexeyev, Nakayasu, and Limantara.It was stated that these SUH methods still need further investigation.According to *Corresponding author: bobbyminola.g@unpar.ac.idGrimaldi et al. [9] and Singh et al. [10], the SUH methods can be used for midscale basins (7-10 km 2 ) but in some cases are also applicable to large-scale ones (> 10 km 2 ).However, in our opinion, this must be investigated as the SUH methods include empirical formulas, and therefore, they may not readily applicable even to other similar cases.
Snyder [5] stated that the Snyder method is applicable only to the Appalachian highlands.Hence, this method needs adjustments for applications in the other areas.According to Sahu et al. [11], there are no clear guidelines in determining the antecedent moisture conditions (AMC) for the SCS method.Therefore, this can be misinterpreted for some conditions.In some studies, the GAMA-1 method produced significant discrepancies for the observed discharge of up to 80% [12][13].This was expressed by Limantara [14] that ideally, each basin should have its own unit hydrograph, and thus, it can be used for the other basins with similar characteristics.This becomes the foundation in developing the Limantara method that was tested for some locations in Indonesia such as Java, Bali, Lombok, and Kalimantan.However, Andiese [12] and Kristianto et al. [13] showed that the Limantara method still produced significant discrepancies to estimate the observed discharge.
In this study, instead of using the SUH methods, we applied the shallow water modeling to the prediction of flood discharge.The physical process of the rainfallrunoff phenomena is calculated directly by solving the shallow water equations (SWE).The solution of the SWE is obtained by means of numerical modeling with second-order accurate Godunov cell-centered finite volume (CCFV) model and fourth-order Runge-Kutta explicit scheme.Among the Riemann solvers, e.g., the Roe, HLL, and HLLC schemes, which are popularly used for shallow water modeling, we particularly interested in developing a Riemann-solver-free method for computing the convective terms of the SWE.This is an artificial viscosity (AV) technique, devised with a combination of a Laplacian and a biharmonic operator.The spectral radius of the Jacobian matrix is used as the variable scaling factor.The rainfall-runoff computations were performed in this research in the framework of an in-house code NUFSAW2D (Numerical Simulation of Free Surface Shallow Water 2D) developed by the author of this paper.This code has been tested and proven accurate in several works [15][16][17][18][19][20][21][22][23][24].
For this paper, we selected a case study carried out by Kirstetter et al. [25], which investigated the rainfallrunoff process experimentally.We simulated this case with the AV technique.For comparison purpose, six SUH methods (SCS, Snyder, GAMA-1, ITB-1, ITB-2, and Nakayasu) were employed.Choosing this case study is due to two reasons.First, there are uniform rainfall data simultaneously measured with the discharge ones, without any obstacles in the basin, thus allowing us to make a fair comparison between the shallow water model and the SUH methods in directly transforming the rainfall to the runoff values.Secondly, the basin is made of impermeable material, allowing for no infiltration, thus the measured rainfall can directly be considered as an effective rainfall value and the measured discharge as a direct runoff.
This study is expected to be an interesting example of how modern numerical computations can improve the runoff prediction.It is also expected that the shallow water modeling may be considered as a standard technique for runoff computation in the future.

Mathematical formulation
The full version of NUFSAW2D employs the 2D-3D non-hydrostatic shallow water equations (NH-SWEs) [22].However, only the 2D hydrostatic SWEs are employed for this study, and thus, the 2D SWEs can be written as Equation 1 where Q, C x , C y , S b x , S b y , S f , and S o are the vectors denoting the conservative variables, convective fluxes, bed-slope terms, bed friction term, and rainfall/infiltration term, respectively.Each vector consists of several variables, written as Equations 2-9.
The variables h (m), u (m/s), and v (m/s) are the depth and velocities in x and y directions, respectively, g (m/s 2 ) is the gravity acceleration, z b (m) is the bed contour, n m (s/m 1/3 ) is the Manning coefficient, and R (m) and I (m) are the rainfall and infiltration, respectively.
For the spatial discretization, the second-order cellcentered finite volume (CCFV) method with the Gauss divergence theorem is applied to Equation 1 over a cell domain Ω e , resulting in Equation 10.
The symbols  x and  y are the unit outward vectors normal to edge i, computed in counter-clockwise direction.
In terms of the temporal discretization, Equation 10 is solved using the Runge-Kutta fourth order method.For the sake of brevity, only the convective fluxes of Equation 10 are explained in this paper.The convective fluxes (second term of the left-hand side) of Equation 10is simply calculated as Equation 11.
where C X  and C x  are the left and right states of the cell fluxes C x for its edge i.A similar way is applicable to C y  and C y  .Unless using the Riemann solvers (e.g., the Roe, HLL, and HLLC methods) or the other Riemannsolver-free method such as the central-upwind scheme, Equation 11 is unstable when discontinuity appears.To deal with this problem, the AV technique is employed so that Equation 11 now becomes Equation 12.
where D x and D y are the artificial diffusive terms for edge i.These terms are constructed using the Laplacian (d 2 ) and the biharmonic (d 4 ) operators of the vector Q at edge i.For example, the term D x for edge i is written as Equation 13.
where ∆Q i and ∆ 4 Q i are the differential operators, and  i (2) and  i (4) are the Laplacian and biharmonic coefficients, respectively.Based on Jameson et al. [26], the artificial diffusive terms must be devised in a flexible, uncomplicated manner.This can be done in such a way that fits the dimensionality of the convective fluxes.Following Ginting and Ginting [22], the variable Λ i max is therefore formed using the spectral radius of the Jacobian matrix computed as Equation 14 where the subscripts j and k denote two adjacent cells of edge i.
The boundary conditions are treated based on the flux computation.To this regard, the characteristic method is applied to the outflow boundary as supercritical flow, for which the depth and velocities of the right states are set equal to the left ones.For solid boundary, we apply the ghost cell technique ensuring that no flow emerges through walls.Another boundary condition used is the rainfall intensity values that are given to all the computational cells.Meanwhile, the infiltration is set to zero.

Results
In this paper, our shallow water model is compared with an experimental case study covering an ideal configuration of overland flow by rainfall.As previously mentioned, we selected this case because a fair comparison between NUFSAW2D and several SUH methods to the observed discharge values can be made by reducing some uncertainties, if real cases are selected.For instance, real cases often include complex topography, which are not always known.In addition, rainfall data are usually available only as point values and not measured everywhere in the catchment area.Also, real cases include many uncertain physical parameters such as infiltration, erosion, deposition, etc., which need special investigations in order to fully understand them.
According to Kirstetter et al. [25], the laboratory experiment was carried out at the Rainfall Simulation Hall of the French Institute for Agricultural Research (INRA, Orléans, France).A flat, rectangular flume channel with a length of 4.04 m and a width of 11.5 cm was tilted by several slope values (here 2% and 5% are selected).A nozzle-type rainfall simulator was located above the channel and used to produce the rainfall.In this study, a constant rainfall intensity of 22 mm/hour during 600 s was chosen.The outflow discharge was recorded at the beginning and after the end of the rainfall (the whole run) by collecting the water in a bucket through a funnel.The sketch of this experimental setup is shown in Fig. 1. [25].

Fig. 1. Overland flow produced by Homogeneous Rainfall Intensity on a Tilted Flume Channel After
The cumulative weight of the water was measured using an electronic scale (a range of 30 kg with 1 g resolution at a rate of 10 Hz).The derivative of such cumulative weight was used to produce the outflow hydrograph.With this manner, the produced hydrograph was quite noisy (small weight increment of 7 g/s).Thus, a Kalman filter was applied to smoothing the hydrograph [25].
The results for both cases (the slope values of 2% and 5%) using NUFSAW2D and the six SUH methods are presented in Fig. 2 within 1,000 s timeframe.Note that we used the Manning value of 0.025 s m -1/3 for NUFSAW2D like employed by Kirstetter et al. [25].Generally, one can see that the AV technique agrees with the observed data by predicting the hydrograph accurately.All the hydrograph parameters can be wellrepresented by this technique.First, the rising limb are captured appropriately until it reaches the steady state condition.This indicates that the time-to-peak as well as the peak flow are computed properly.As the rainfall intensity becomes zero (after 600 s), the decreasing limb starts, which can be predicted accurately by the AV technique.It can also be noticed from Fig. 2 that the peak flow increases as the channel slope increases.
Fig. 2 also shows the results of SUH methods, which were computed using a single basin analysis.Interestingly, none of the SUH methods is able to predict the rising or decreasing limbs.For the slope value of 2%, the SCS method can achieve the steady state condition; however, the time-to-peak is inappropriately captured.Meanwhile, the other SUH methods fail to achieve the steady state condition.Indeed, the Snyder method can compute the peak discharge but not the time-to-peak.Similarly, for the slope value of 5%, the SCS method can still achieve the steady state but predict a smaller peak discharge value.A similar phenomenon is also experienced by the other SUH methods, which not only fail to compute the rising limb but also the steady state condition and the decreasing limb.In order to present the errors produced by the AV technique and the SUH methods, we use two indicators: (1) Root Mean Square Error () and (2) Relative Error (), which are computed respectively as Equations 15-16 where  is the number of data,  denotes the -th variable,   is the discharge obtained from the computation and   is the observed discharge.Equation 16may result in a negative value that indicates a smaller value than that of the observed data.Note that  is employed to assess the computed hydrograph shape, while  is used to assess the computed time-to-peak 0.0 0.5  It is shown from Table 1 that the AV technique can produce the smallest  values among the SUH methods.This indicates that the AV technique can give the most accurate results among the others for predicting the direct runoff shape.In Table 2 one can notice that the AV technique is again able to produce the most accurate results among the other methods for predicting the time-to-peak and the peak discharge.Meanwhile, among the SUH methods, ITB-1 and Nakayasu become the best methods.
Note that the accuracy of the numerical model is affected by the bed friction model used.In this research, we employed the friction law of Manning.According to Kirstetter et al. [25], the Poiseuille friction model should be used to obtain more accurate results for the velocity and thus the discharge.Also, the Darcy-Weisbach model is applicable to the prediction of rainfall-runoff if the flow regime is laminar.Nevertheless, the Manning model is still used for this study because it is quite simple to be integrated into the time integration scheme for Equation 1. Obviously, both Poiseuille and Darcy-Weisbach models can be combined with the AV technique.
From Tables 1-2, one may also observe that there is none of the SUH method being consistent and able to compute the direct runoff accurately, even for a very simple domain like investigated in this study.This shows a notable flaw of such conventional SUH methods to be a reliable technique for direct runoff prediction.

Conclusion
The AV technique has been presented in this research for direct runoff prediction.This technique was devised based on a combination of a Laplacian and a biharmonic operator, which was then used for the calculation of the convective fluxes in order to avoid the instability due to the averaging process of the left and right states of such fluxes.
As a case study, a laboratory experiment was selected, which was equipped with a nozzle-type rainfall simulator and a bucket to measure the discharge (through a funnel), creating overland flow on a tilted flume channel.The flume was made of impermeable bed material, thus allowing for no infiltration.Therefore, the measured rainfall can directly be considered as an effective rainfall value and the measured discharge as a direct runoff.
As a comparison to the AV technique, six conventional SUH methods, i.e., SCS, Snyder, GAMA-1, ITB-1, ITB-2, and Nakayasu, were also used to compute the direct runoff for the case study selected.This analysis was carried out using single-basin computation.Predicting three parameters became the focus of this study, namely the time-to-peak, the peak flow, and the hydrograph shape.Two error indicators, i.e., Root Mean Square Error () and Relative Error (), were used to assess the accuracy of the AV technique and the SUH methods.
The simulation results showed that the AV technique could outperform all the SUH methods by giving the smallest  and  values, thus being accurate in predicting the hydrograph shape.In addition, the AV technique could predict the time-to-peak and the peak discharge significantly more accurate than the SUH methods.In contrast, we found that none of the SUH method was consistent and able to compute the direct runoff accurately.
This study reveals a notable flaw of the SUH methods that even for a very simple domain, such methods were neither reliable nor accurate for direct runoff prediction.In contrast, the SWE model with an AV technique can be used for an accurate prediction of direct runoff.One aspect worth investigating for the future study is the effect of the bed friction model used with the AV technique because the water depth of the direct runoff is very low, and thus, the bed friction plays a significant role for the numerical computation.
Finally, we would like to point out in this study an interesting example of how modern numerical computations can improve the runoff prediction.Therefore, we expect that the shallow water modeling may be included as a standard technique for runoff computation in the future.
The author acknowledges Cleon Christopher for his help in providing the results of the SUH methods.
C x  x + C y  y ] i N i=1 ∆L i = ∬ (S b x + S b y + S f + S o ) dΩ e Ω e

Table 1 .
values for the hydrograph shape.

Table 2 .
values for the time-to-peak and peach discharge.