One of the approaches to flood modelling is numerical simulation of the diffusive wave approximation of the shallow water equations. Improving these models in various aspects is still an open area of research. In this study, a new diffusive wave model with explicit time integration was developed which includes some novel features: (1) time steps are determined using a novel stability criterion which resulted in more dynamic time steps (i.e., broader range) compared to the conventional Courant–Friedrichs–Lewy stability condition; (2) stability constraints are reduced, considering the flow processes within surface ponds; (3) besides Manning's formula, which is the common equation for computing velocities in diffusive wave models, the free fall velocity and a new equation for wave-front velocity are employed; and (4) the influence of upstream surface ponds on downstream flow is considered. This paper introduces the enhanced diffusive wave model, the so-called Overland Flow Simulator Cellular Automata (OFS-CA), and its results for five test cases. Available analytical solutions and an experimental study were used for verification. Two other shallow water models were used for comparison and benchmarking. Overall, good agreements were observed and OFS-CA was computationally less expensive compared to the other two shallow water models.

## ABBREVIATIONS AND NOTATION

- CA
cellular automata (singular: automaton)

- CFD
computational fluid dynamics

- CFL
Courant–Friedrichs–Lewy

*Cr*Courant number (dimensionless)

- LTS
local time stepping

- SWE
shallow water equations

- 2D
two-dimensional

*d*water depth [L]

effective water depth [L]

*g*gravitational acceleration [LT

^{−2}]*i*index of central cell

*j*index of basin cell (integer values between 1 and 4)

*n*Manning's roughness coefficient [TL

^{−1/3}]*Nb*number of basin cells assigned to a central cell

*outV*outflowing water volume [L

^{3}]total discharge (total outflow) [L

^{3}T^{−1}]total inflow [L

^{3}T^{−1}]lateral discharge (lateral outflow): from cell

*i*to basin cell*j*[L^{3}T^{−1}]*x*-direction discharge [L^{3}T^{−1}]*y*-direction discharge [L^{3}T^{−1}]*Ref*index of reference cell

*S*water surface gradient (dimensionless)

ground slope (dimensionless)

cross-sectional averaged flow velocity [LT

^{−1}]*wl*water level: sum of ground elevation and water depth [L]

*z*ground elevation [L]

water depth variation [L]

size of time step [T]

difference of water levels [L]

cell size: distance between cell centroids [L]

threshold on water surface gradient (dimensionless)

## INTRODUCTION

Pluvial flooding is an unavoidable natural hazard that may cause a threat to life and damage to properties and infrastructures. In recent decades, assessment of this phenomenon has been facilitated by enhancements in computational fluid dynamics (CFD). In the framework of CFD, various models have been developed to simulate flooding two dimensionally by *numerical simulation* of the shallow water equations (SWE). These models are often referred to as 2D flood models.

Besides those 2D flood models which solve the *fully dynamic* SWE (see Hunter *et al.* 2008; Hinkelmann *et al.* 2015), *diffusive wave* models (also called *zero-inertia*) are common in which the SWE are simplified by neglecting the inertia terms in the momentum equation (see e.g., Bradbrook *et al.* 2004; Ettrich *et al.* 2005; Wang *et al.* 2011; Costabile *et al.* 2012; López-Barrera *et al.* 2012). This diffusive wave approximation simplifies the computation of flow velocity, leading to computational gains. Nevertheless, some recent studies (Hunter *et al.* 2008; Cea *et al.* 2010) reported higher computational costs for diffusive wave models compared to fully dynamic models. Regarding these reports, it should be noted that most likely it is not the diffusive wave approximation which has been responsible for the longer run times, but the different stability criteria of the models: in both mentioned studies the fully dynamic model applied the Courant–Friedrichs–Lewy (CFL) stability condition, whereas the diffusive wave model used the von Neumann stability analysis introduced by Hunter *et al.* (2005), which is reported by Wang *et al.* (2011) and Mendicino *et al.* (2015) as being stricter than the CFL condition. It should be mentioned that the methodology in Hunter *et al.* (2005) was later enhanced by Bates *et al.* (2010) and they changed their approach from von Neumann to CFL condition.

As pointed out above, the efficiency of flood models is not only affected by the level of simplification of the SWE, but also by the considered stability criterion. Stability is a condition on the numerical solution, namely that all errors must remain bounded as computation advances in time (Hirsch 2007). Various criteria have been implemented in flood models to fulfil the stability by controlling errors. In this regard, *explicit* schemes commonly restrict the size of time steps using the CFL stability condition, which requires a dimensionless Courant number (*Cr*) to be specified.

The level of restriction of the CFL stability condition is influenced by the choice of *Cr*. In some studies, CFL condition has been combined with correction parameters, resulting in computational gains: Ghimire *et al.* (2013) applied a large Courant number of *Cr* = 1 and ensured the stability further via a relaxation parameter that has to be calibrated, i.e., the relaxation parameter damps the flow oscillations which occur due to the size of the time step not being sufficiently restricted. A similar approach is presented in Parsons & Fonstad (2007), where *Cr* = 1 is used and the flow oscillations are smoothed out with calibrating the unknown physical parameters.

The other factor which influences the efficiency of flood models is their time integration method (implicit, explicit). Implicit schemes are not the concern of the present research. However, for the sake of completeness, we refer to other studies where implicit diffusive wave models are introduced and compared with explicit ones (Lal 1998; Lal & Toth 2013; Fernández-Pato & García-Navarro 2016).

In recent years some flood models have been named cellular automata (CA). With respect to this designation the bridges between the definition of CA and shallow water models should be noticed (sometimes they are equivalent). A cellular automaton (pl: automata) is an array of cells, where each cell has a state and at each computation step new states are determined according to some transition rules. The cell states, the transition rules, and other properties of CA such as initial and boundary conditions are defined depending on the problem to be simulated (Weimar 1997), keeping in mind that CA are not limited to rectangular grids and the transition rules do not necessarily have to define an exact outcome (Wolfram 2002; Shiffman 2012). Surface flow can be simulated by setting up CA considering the water depths as the cell states. The transition rules might be defined without solving SWE and by integrating physical parameters which require calibration (Thomas & Nicholas 2002; Rinaldi *et al.* 2012). They also might be defined according to the SWE (whether fully dynamic or diffusive wave) and according to the stability criteria (Parsons & Fonstad 2007; Dottori & Todini 2011; Ghimire *et al.* 2013; Mendicino *et al.* 2015). In the latter case, CA are equivalent to shallow water models with explicit time integration, which was also the case in the present research.

In this study, a diffusive wave model with explicit time integration and of first order accuracy was developed. The model is called Overland Flow Simulator Cellular Automata (OFS-CA) and has the following novel features:

Feature

*MJ-t*: a new stability criterion is developed which restricts the size of global time step to meet a particular local criterion, while avoiding negative water depths elsewhere. This criterion takes advantage of the Bernoulli equation and does not rely on calibration or correction parameters.*Reducing stability constraints*: to deal with the limitation of global time stepping, certain cells are dismissed from computations and stability checks. In spite of local time stepping (LTS) techniques (see Crossley*et al.*2003), which categorize the cells based on size of time step, in this approach cells are grouped according to topographical criteria (i.e., surface ponds).Feature

*MJ-v*: in some conditions Manning's formula, which is the common equation to compute velocities in diffusive wave models, is replaced with either the free fall velocity or a new equation. This new equation was developed to estimate the velocity of dam-break induced wave-front and, in spite of the conventional equations (e.g., Ritter 1892; Lauber & Hager 1998), is independent of time and distance from the dam.Feature

*MJ-p*: with the purpose of considering the influence of upstream surface ponds on downstream flow, the so-called*MJ-p*algorithm is developed to recognize water volumes which are trapped within surface ponds as simulation proceeds.

The present paper introduces the enhanced model and the results of subjecting it to a number of analytical, experimental and benchmark tests for the purpose of verification. The model's results and performance are compared with those of two shallow water models: a second order accurate fully dynamic model and a first order accurate diffusive wave model, both applying CFL stability condition and explicit time integration. Additionally, a direct comparison between the enhanced MJ-t criterion and the CFL stability condition is conducted in terms of the resulting size of time steps. The results are then discussed and conclusions drawn.

## METHOD

*basin*is formed for each cell consisting of those adjacent neighbours with a lower water level. Very small water level differences are ignored by applying a user-specified threshold on water surface gradient . This was set to for the results presented in this paper. Figure 1 illustrates a neighbourhood of cells. The basin of the central cell consists of those two neighbour cells with lower levels. The importance of defining basins is that further computations (including stability checks) take place only for the cells which have at least one basin cell.

### Temporal depth variations

^{2}T

^{−1}], and v is the velocity vector [LT

^{−1}].

*d*and

*z*are water depth [L] and ground elevation [L], respectively, and g stands for the gravitational acceleration [LT

^{−2}]. In a

*diffusive wave approximation*the momentum equation is simplified to by neglecting local acceleration and advection terms . By applying resistance equations the external forces are often formulated as . Hence, the momentum equation is reduced to a simple algebraic representation of the 2D velocity vector:

Considering that term represents water surface gradient, Equation (2) indeed represents the cross-sectional average velocity at interfaces according to Manning's formula where the hydraulic radius is approximated by the depth of water.

^{3}T

^{−1}] and total outflow [L

^{3}T

^{−1}] (also called total discharge) including sink and source terms, and is the size of cell. In further descriptions we leave out the superscript

*Tout*and indicate the total discharge by . The total discharge at cell

*i*is computed by summing up the lateral outflows: where stands for the number of basin cells. is the outflow from cell

*i*to basin cell

*j*, and is the cross-sectional averaged velocity. stands for the ‘

*effective*’ depth of water including sink and source terms and is described shortly later. Considering that the computed outflows between neighbour cells (Equation (4)) contribute to the two cells with opposite signs, the total inflow to cell

*i*, required for Equation (3), is computed by summing up the corresponding outflows of its neighbour cells.

*d*) if ground elevation of the cell is lower than the ground elevation of its basin cell (Figure 2). In this case, only the surplus is taken into account in discharge computations. An example is displayed in Figure 2, where and show the effective water depths for computing outflows from cell

*i*to basin cells 1 and 2, and a volume of water with a depth of remains immobile at cell

*i*. In other studies, similar approaches can be found where ‘only the depth through which water can flow’ has been considered to compute discharges at cell interfaces (e.g., Hunter

*et al.*2005; Bates

*et al.*2010; Dottori & Todini 2011; Wang

*et al.*2011; Ghimire

*et al.*2013). In the present work, the concept of effective water depth is not restricted to the cell interfaces, but is expanded to the surface ponds, i.e., as depicted in Figure 3 (top), the accumulated water within a surface pond is considered being trapped. In doing so, the algorithm

*MJ-p*was developed to specify the trapped water volumes as the flow simulation proceeds. The imaginary sequence in Figure 3 (bottom) describes the principle of this algorithm. Notice the flow from cell

*b*to

*c*in step 2, and from cell

*c*to

*d*in step 5. This shows that in the case of an obstacle (i.e., ground surface) against flow, water remains immobile behind it (displayed darker). The immobile water in cell

*c*(step 5) blocks the flow further and water remains trapped in cell

*b*as displayed in step 6. The reduction in the effective depth of water at cell

*b*, as simulation proceeds from step 5 to 6, will be observed as a drop in the discharge hydrograph of cell

*b*located within the pond.

### Velocities

In computing the cross-sectional averaged velocity , required for discharges in Equation (4), the following groups of cells are handled differently:

(a) water surface gradient ≤ ground slope,

(b) water surface gradient > ground slope,

*vg*[LT

^{−1}], is also considered in order to improve the accuracy of estimated flow velocities:

*vg*is a 1D dam-break, where Manning's formula overestimates the wave-front velocity significantly. In order to solve this issue, in an initial effort the speed of a mass freely falling from distance , i.e.,, was considered as a threshold for velocity, with the explanation that if the water depth at cell

*i*belonging to group (b) exceeds a certain threshold, the particles at the top of the water column have a dominant tendency to free fall rather than a flow driven by the bed roughness. This approach significantly improved the temporal change of flow depth and the shape of the wave-front in a 1D dam-break test case, compared to when only Manning's formula is used to estimate velocities. However, the speed of propagation of the wave-front was underestimated. Therefore, the formulation of

*vg*was further enhanced as below:

is the depth of a column of water at rest and is specified for simulating dam-break, where the stored water in a reservoir collapses after removing the gate. stands for the depth integration of the free fall velocity which is a function of the fall distance (*h*).

*term*1 and *term*2 were developed to represent the *dynamic* and *initial* waves, respectively. Lauber & Hager (1998) described the ‘dynamic wave’ as the wave ‘originating from the surface particle at the dam section’, and the ‘initial wave’ as the wave which ‘propagates from the bottom of the dam section, much in analogy to gate or orifice flow’. According to Lauber & Hager (1998), after an extremely small time period the flow depth at the dam section will reach . Therefore, in the present research, the free fall velocity was integrated from 0 to to derive *term*2, as shown in Equation (7). The orifice velocity and Manning's formula are included as thresholds on velocity (also through formulation of ) to account for bed friction, resulting in an agreeable formation of tip region and a good estimate of flow velocity not only at locations near the dam section but also far downstream (see Results and discussion). For ease of following the general numerical scheme we skip further details about the development of Equation (7) and only point out, that in spite of the available analytical solutions for the wave-front (e.g., Ritter 1892; Lauber & Hager 1998), this new equation is independent of time and distance from the dam. We refer to the above described approach as feature *MJ-v*.

### Time step

*reference cell*and to avoid negative water depths elsewhere. This method, which we refer to as

*MJ-t*, is summarized in Figure 4. As can be seen, first the cell with the maximum total discharge (

*Q*) (Equation (4)) is identified and referred to as the reference cell. In case of maximum total discharge (

*Q*) occurring at multiple cells with unequal effective water depths , the one with the shallower is chosen. Note that is the depth of the water volume which is available to be transferred and in a case like Figure 2 is given by:

^{3}T

^{−1}] at reference cell and is specified according to a local rule as explained later. Note that the variables at the reference cell are shown with subscript

*Ref*while subscript

*i*is used for any other cell.

*artificial*’ total discharge specified as described later. In summary, two conditions should be fulfilled at each cell:

#### Determining (*outV*_{Ref})^{l} at reference cell

Within the neighbourhood of the reference cell, the basin cell with the maximum water level (*wl*) is identified, referred to as *max basin cell*. The available water at the reference cell is distributed into all of its basin cells (in portions relative to lateral outflows) as long as the water level of the reference cell does not drop below the new water level of *max basin cell*. The volume of water which has left the reference cell by the time this threshold is reached is referred to as , where superscript *l* stands for local.

#### Determining for each cell

*i.*By replacing the lateral effective depths with , the right hand side of Equation (12) can be expanded further: where is the sum of magnitude of lateral velocities . Using the specific energy (refer to Bernoulli equation) a relation is formulated between the sum of magnitude of velocities at reference cell and at cell : Finally, the artificial total discharge is given by placing Equation (14) into Equation (13):

#### Reducing stability constraints

Inspired by LTS techniques (Crossley *et al.* 2003), in OFS-CA some computations (including stability checks) are skipped for a certain group of cells. However, in spite of LTS methods which categorize cells based on the size of time steps, the classification in OFS-CA is simply based on surface topography and overland flow processes: the water flowing into a surface pond accumulates there until the storage capacity of the pond is reached and water overflows. Accepting some oscillations in the flow within ponds during the process of filling up can lead to computational gains, keeping in mind that the influence of these oscillations on the overland flow outside the pond should be (and can be) kept to a minimum.

To allow some minor oscillations within partially filled ponds, the threshold on water surface gradient is increased there. is set back to the smaller value as soon as the filled level is reached. A larger value of reduces the number of basin cells (Figure 1). This improves the performance not only due to skipping some discharge computations, but more importantly, due to a reduction in the stability constraints which allows larger size of time steps: dismissing some basin cells of a cell *i* is equivalent to ignoring the corresponding lateral outflows and accordingly reducing the total discharge . If cell *i* was to be a constraint on stability due to Equation (10) not being satisfied, a smaller value of might now satisfy Equation (10) for cell *i* (note that if this equation is not fulfilled, the size of time step is reduced, as described previously).

This approach might lead to a wavy flow (i.e., oscillations) within ponds during the process of filling up, but the oscillations are expected to damp after is set back to the smaller value, unless the choice of is too large (then the oscillations will be transmitted downstream). For the results presented in this paper an increase factor of 10 was applied (i.e., from 5 × 10^{−4} m to 5 × 10^{−3} m).

### Initial conditions

The initial condition is constructed by assigning a water depth of 5 × 10^{−4} m to those cells of the domain which are imposed to rainfall.

### Boundary conditions

Those cells at the edge of the domain are referred to as *boundary cells*.

At

*closed boundaries*, sufficiently high ground elevations are assigned to boundary cells to ensure that no water flows into them.*Open inflow boundaries*are not considered at this stage of the work, i.e., the domain receives no inflow from the boundary cells.- At
*open outflow boundaries*, the volume of water received by boundary cells is stored as the domain's outflow. At the end of each computation step, artificial water levels are assigned to each boundary cell, so that: with*B*and*i*as the indices of boundary cell and its upstream cell, respectively. and stand for water surface gradient and ground slope, described in Equation (5).

## RESULTS AND DISCUSSION

The following sections describe five test cases presented in order of complexity. The first two are 1D cases. The simplest one (Case 1) is an inclined open channel imposed to rainfall. Results were compared with those of Hydroinformatics Modelling System (HMS) intended for benchmarking. HMS is a robust fully dynamic shallow water model based on finite volume method with second order MUSCL reconstruction (see Simons *et al.* 2014). In Case 2, a 1D dam-break flow in the absence of a bed slope was modelled. For this case, experimental data were available to verify the results of OFS-CA.

Case 3 is a V-shaped schematic catchment consisting of an inclined channel receiving inflow from the two hillsides. Results were verified using an analytical solution and were compared with those of HMS.

Cases 4 and 5 are 2D real-case watersheds with different sizes, which allowed assessment of the capability of OFS-CA not only in 2D flow simulation but also in handling the complexities of terrain. For the purpose of benchmarking, Cases 4 and 5 were modelled with both HMS as well as to HE2D (itwh). The latter simulates the diffusive wave SWE using the first order numerical scheme described in Ettrich *et al.* (2005) (Jahanbazi & Egger 2014). In both HMS and HE2D the size of time steps are determined according to CFL stability condition.

For Cases 4 and 5 the overland flow was also simulated with OFS-CA** _{CFL}** which differs from OFS-CA solely in its stability condition, in order to particularly assess the influence of MJ-t method on model performance. A summary of simulation run times of all test cases is also presented.

### Case 1: inclined open channel

As the first test case, a simple open channel (length = 100 m, width = 10 m) with a forward slope of 1% and a Manning's roughness coefficient of 0.02 sm^{−1/3} was modelled with a cell resolution of 1 m. All boundaries are closed except for the channel outlet. The initially dry channel was imposed to a constant rainfall with intensity of 10^{−5} m/s and duration of 1,800 s. Simulations ran for 1 hour. The channel was modelled also with HMS with the same grid and parameters.

This simple example showed that OFS-CA simulates variation of flow depth and discharge in open channel as good as a second order accurate fully shallow water model (HMS), while being considerably faster (see Table 1 in section Performance comparison).

OFS-CA | HMS | HE2D | OFS-CA_{CFL} | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Courant number (Cr): | 0.1 | 0.3 | 0.05 | 0.08 | 0.3 | 0.02 | 0.05 | 0.1 | 0.3 | |

Velocity in CFL condition: | ||||||||||

Case 1: Open channel; Run: 1 h | 0.2 s | 3.7 s | 1.4 s^{a} | |||||||

Case 2: Dam-break (10 cm, 1 cm); Run: 18 s | 41 s, 195 s | |||||||||

Case 3: V-catchment; Run: 3.3 h | 16 s | 295 s | 134 s^{a} | |||||||

Case 4: Small real-case; Run: 800 s | 0.6 s | 20 s | 6.3 s | 0.8 s | 0.6 s^{b} | |||||

Case 5: Medium real-case; Run: 3 h | 1.8 h | 16.9 h | Instable | 4 h | 2.5 h^{c} | Instable | 2.7 h | 1.1 h^{c} | Instable |

OFS-CA | HMS | HE2D | OFS-CA_{CFL} | |||||||
---|---|---|---|---|---|---|---|---|---|---|

Courant number (Cr): | 0.1 | 0.3 | 0.05 | 0.08 | 0.3 | 0.02 | 0.05 | 0.1 | 0.3 | |

Velocity in CFL condition: | ||||||||||

Case 1: Open channel; Run: 1 h | 0.2 s | 3.7 s | 1.4 s^{a} | |||||||

Case 2: Dam-break (10 cm, 1 cm); Run: 18 s | 41 s, 195 s | |||||||||

Case 3: V-catchment; Run: 3.3 h | 16 s | 295 s | 134 s^{a} | |||||||

Case 4: Small real-case; Run: 800 s | 0.6 s | 20 s | 6.3 s | 0.8 s | 0.6 s^{b} | |||||

Case 5: Medium real-case; Run: 3 h | 1.8 h | 16.9 h | Instable | 4 h | 2.5 h^{c} | Instable | 2.7 h | 1.1 h^{c} | Instable |

^{a}Minor overshooting at the moment of reaching the maximum discharge.

^{b}Some oscillations at control point 3.

^{c}Some oscillations at control points 1 and 2.

### Case 2: dam-break flume

*et al.*(2009, 2011). They set up a flume with bed slope = 0, length = 6 m, width = 0.25 m and height = 0.6 m. The upstream reservoir has an initial water depth of

*h*

_{0}= 0.325 m, while downstream is dry with an open outlet. They captured the images of flow after gate removal. Figure 6 (top) illustrates the images of flow at two times, t = 0.9 s and t = 0.18 s. The noticeable particles were added to water for the purpose of extracting velocity fields as described in Aleixo

*et al.*(2011). Additional to image capturing, Baumer

^{®}ultrasonic probes were placed at certain locations above the channel to measure water levels during the course of time. This allowed capture of the temporal change of water depth, as displayed in Figure 7, with non-dimensional coordinates: , where

*x*is the distance from gate.

The flume was modelled with OFS-CA on two different grids:

3 by 62 square grid, resolution: 10 cm,

3 by 602 square grid, resolution: 1 cm.

Considering the flume's bed materials, Manning's roughness coefficient was set to 0.012 sm^{−1/3} for the reservoir (polished wood) and 0.01 sm^{−1/3} for downstream (glass). Simulations ran for 18 seconds, i.e., up to T = 100.

Flow profiles at several instances, simulated by OFS-CA, were compared with the experimental data. Two of the profiles are illustrated in Figure 6, where x = 0 m is the location of gate. At t = 0.18 s (T = 1) the water levels obtained by OFS-CA show promising agreement with the flow profile extracted from raw images. The tip region is reasonably formed, specially in the 1 cm grid, showing that the introduced formulation in Equation (7) accounted for bed friction to a good extent. This could be considered as an advantage over the classical Ritter's solution: (Lauber & Hager 1998), which is also displayed in Figure 6.

At t = 0.9 s, OFS-CA overestimated to some extent the water depths in the reservoir and near the gate compared to experimental data and Ritter's solution: around 3 to 4 cm depth difference compared to the experimental data, and around 5 to 6 cm depth difference compared to Ritter's solution. The overestimation is observable not only in Figure 6 (right), but also in Figure 7: at gauge X = −0.15 (i.e., x = −0.05 m) located right before the gate, water levels computed by OFS-CA for 1 < T < 15 are higher than the levels obtained from the experimental study, although they agree for T ≤ 1. However, at downstream gauges, depth hydrographs matched very well with the experimental data. Results of the 1 cm grid and 10 cm grid converged well, especially near the gate and at far downstream (X = 3.69, i.e., x = 1.2 m).

This test case demonstrated that Equation (7) is suitable to simulate the wave-front in the absence of a bed slope. The transition from wave-front to normal flow occurred smoothly, since the magnitude of and Manning's velocity converged around T = 50.

### Case 3: V-shaped schematic catchment

*et al.*(1996) produced this solution for a catchment with the geometry illustrated in Figure 8. Later on, Simons

*et al.*(2014) reproduced this analytical solution, which we used to verify the results of OFS-CA, while employing HMS as the benchmark model. Manning's roughness coefficient was set to 0.015 sm

^{−1/3}for the hillsides and 0.15 sm

^{−1/3}for the channel. A constant rainfall with intensity of 10.8 mmh

^{−1}and duration of 1.5 h was applied to the initially dry hillsides. Simulations ran for 3.3 h. All boundaries are closed except for the channel outlet. The grid consists of 164 by 101 cells with a resolution of 10 m. In this case, compared to the simple open channel in Case 1, some additional conditions are considered. The dry channel receives rainwater in the form of inflow from the two hillsides. Flow direction changes when run-off flows from hillsides into the channel. There is a considerable drop between hillsides and channel, which varies along the channel length: 1 m upstream and around 20 m downstream.

Figure 9 (bottom) demonstrates the temporal change of water depths. At channel outlet, OFS-CA and HMS produce identical depth hydrographs, but at hillside edge, in spite of agreement in discharge computations, they differ in respect of flow depths: in OFS-CA, water depths at the hillside edge vary along the channel length, whereas in HMS, they are constant while being larger than those of OFS-CA (Figure 9 (bottom)). The water depths obtained by OFS-CA were considered as the reliable ones since they reasonably demonstrated the effect of varying ground slopes and also it is already known that the positivity preserving hydrostatic reconstruction implemented in HMS introduces some numerical errors for cases with large discontinuities (Simons *et al.* 2014).

This test case showed the capability of OFS-CA in capturing the dynamics of build-up and subsidence of flow depth and discharge in the considered special conditions. Flow was simulated consistently (i.e., smooth hydrographs) and in excellent agreement with the available analytical solution. At the channel outlet, results of OFS-CA and HMS conformed perfectly in terms of both flow discharges as well as flow depths. At the hillside edge, in spite of excellent agreement in flow discharges, the models differed in computation of flow depths: as mentioned previously, those obtained by OFS-CA were considered as the reliable ones. Furthermore, as listed in Table 1 (see section Performance comparison), run time of OFS-CA was considerably shorter than that of HMS.

### Case 4: small real-case watershed

^{2}was modelled with a grid of 23 by 38 cells and resolution of 1 m. In an initial step, the terrain topography was analysed by means of ArcGIS Spatial Analyst toolbox. The generated potential flow paths (where flow may occur or not, depending on the rainfall event), surface ponds, and the respective digital elevation model (DEM) are illustrated in Figure 10(a). All boundaries are open and Manning's roughness coefficient was set to 0.02 sm

^{−1/3}.

The initially dry domain was imposed to a constant rainfall with duration of 300 s and intensity of 3 × 10^{−5} m/s, and overland flow was simulated for 800 s. Illustrating water depths at different time intervals allowed inspection of the flow propagation. An example is Figure 10(b), showing water depths at t = 350 s.

Control points 1 and 2 allowed inspection of the models in dealing with negative bed slopes: OFS-CA considers only the effective water depth in discharge and velocity computations (Figure 2), while this is not the case in the second order scheme of HMS. At control point 1, HMS computes small values of , whereas OFS-CA computes no due to water level not having reached the bed elevation of left- and right-hand neighbours. At control point 2, OFS-CA computes only when water begins to flow over the right-hand neighbour (after reaching a sufficiently high level), whereas before this point of time HMS computes greater than zero. This point of time is exactly where a drop is observed in depth hydrograph of HMS.

Control point 3 is a special case, as it is located in a pond. Here, the water depth hydrographs produced by OFS-CA and HMS match perfectly and the discharge hydrographs have similar forms in both

*x*- and*y*-directions. However, the maximum discharge values differ considerably, due to the approach of OFS-CA to ignore the trapped water volumes in discharge computations (feature MJ-p), leading to lower discharge hydrographs compared to HMS.At control point 4 it was observed that the two models vary to some extent in the ratio of lateral discharges even for positive bed slopes. At this control point, where the ground slopes in both directions are zero, the flow simulated by HMS occurs dominantly in

*y*-direction, whereas OFS-CA computes fairly similar discharges in both directions. Therefore, in spite of good agreement in the total discharges (*Q*), there is a considerable difference in the maximum , as displayed in Figure 11. This could be explained by the fact that in the second order scheme of HMS the state variables are extrapolated to the cell interfaces using a three-point stencil-based MUSCL reconstruction (Simons*et al.*2014), which means the variables in one cell are more significantly influenced by the neighbouring cells than in OFS-CA.

A further comparison was conducted using HE2D, introduced in the preface of Results and discussion. As HE2D accepts only an unstructured triangular grid, the regular grid was transformed by keeping the smoothing of elevation data to minimum. The transformed grid consists of 821 triangular cells with an average cell area of 0.5 m^{2}, varying in the range of 0.1–1.5 m^{2}. The general topography of the unstructured grid in HE2D model agrees with that of the regular grid, accordingly, similar flow paths and flow ponding as in Figure 10(b) were formed. However, the special topographic complexities of the control points (described previously) could not be fully reproduced by HE2D: these complexities depend greatly on the ground elevation of adjacent cells, which do not remain identical through the transformation from the regular square grid to the unstructured triangular one. The influence of the mentioned difference in grid geometry is observed in the shape of the rising limb of hydrographs and also in the peak values (Figure 11). However, one general similarity is observed: a good agreement between the forms of receding limbs of the hydrographs obtained by the three models (Figure 11).

This small real-case watershed allowed a detailed assessment of OFS-CA, showing its capability for consistent flow simulation (i.e., smooth hydrographs) in spite of topographic complexities. The hydrographs obtained by OFS-CA concurred very well with those of HMS in terms of their form. The total discharges (*Q*) agreed well except for within ponds, since OFS-CA considers the trapped water volumes as immobile. The total discharges agreed better than the lateral discharges, since as mentioned above, in HMS the influence of state of neighbouring cells on the flow is greater than in OFS-CA. In terms of lateral discharges, results of HMS were more reasonable (second order accuracy), whereas in dealing with negative bed slopes OFS-CA showed an advantage over HMS. A detailed comparison of HE2D and OFS-CA regarding dealing with complexities was not possible, since different grid geometries did not allow generation of identical complexities. However, it was valuable to observe their general agreement in simulating the flood extent and also their moderate agreements in the form of hydrographs. Both OFS-CA and HE2D (diffusive wave models, first order of accuracy) outperformed HMS (fully dynamic model, second order of accuracy), as summarized in Table 1 (see section Performance comparison), while OFS-CA was less expensive than HE2D.

### Case 5: medium real-case watershed

^{2}, varying in the range of 0.25–1.95 m

^{2}. All boundaries are open and Manning's roughness coefficient was set to 0.02 sm

^{−1/3}. First, the terrain topography was analysed by means of ArcGIS Spatial Analyst toolbox. The generated potential flow paths (where flow may occur or not, depending on the rainfall event), ponds, and the respective DEM are illustrated in Figure 12(a).

^{−5}ms

^{−1}(48 mmh

^{−1}) and duration of 1 h. Overland flow was simulated for 3 h. Figure 12(b) illustrates flow depths 30 minutes after rain stops. Animating the flow depths at time intervals of 1 minute demonstrated the flow propagation: rainwater flows into the ponds and accumulates there until they are filled up. The surplus then overflows into downstream ponds. As can be seen in Figure 12, the domain drains mainly into control point 1. At this point, and two additional control points, discharge and water depth hydrographs were compared with those obtained by HMS and HE2D (Figure 13), showing similarity in their general form. The rising limbs of hydrographs have staircase developments. A stair rises (i.e., temporal changes increase) as soon as the surplus of an upstream pond reaches the control point.

Figure 13 demonstrates that the three models agree generally to a good extent, except for one main difference: OFS-CA computes smaller discharges over the filled ponds, as it considers some water volumes as being trapped within surface ponds, which we previously referred to as feature MJ-p (Figure 3). This explains the considerable difference between maximum discharges at control point 2 located within a pond, whereas water depths agree well. As a result of smaller discharges over ponds, in OFS-CA the surplus from upstream ponds reaches the control points with a delay compared to the two other models. This is observable in the hydrographs as a delayed rise of stair (Figure 13), which can be especially recognized at control points 1 and 3. This also explains the steeper receding limbs of the hydrographs obtained by OFS-CA. In HMS and HE2D models, the surplus from upstream ponds still feeds the control points after rain stops, leading to slower temporal change of velocity and discharge. It was investigated whether feature MJ-v (see Equations (6) and (7)) is responsible for the steep receding limbs and not MJ-p, which proved not to be the case: velocities were computed once solely according to Manning's formula, resulting in very similar hydrographs as those presented here, with very slightly larger peak values.

This test case not only showed the consistency and stability of OFS-CA in simulating 2D overland flow, but also allowed inspection of the influence of upstream surface ponds on downstream flow according to feature MJ-p. OFS-CA differed from HE2D and HMS with respect to discharges over ponds and also downstream flow discharges, since it considers some part of overland flow as being trapped in surface ponds (Figure 4). This was observed as delayed rise in hydrographs of OFS-CA and steeper receding limbs. Similar to in the previous test case, OFS-CA agreed very well with HMS in terms of maximum water depth within the pond, but the flow depths at other points were lower, which is very likely due to their different schemes (first order in OFS-CA vs. second order in HMS) which influences the lateral water distribution as described in Case 4: Small real-case watershed. Both OFS-CA and HE2D (diffusive wave models, first order of accuracy) overtook HMS (fully dynamic model, second order of accuracy) with respect to run times (see Table 1), while OFS-CA was less expensive than HE2D.

### Performance comparison

*Cr*) is time and space dependent, often a constant value is specified for the whole simulation over the whole domain.

Table 1 summarizes run times of the models presented in previous sections. The reported times correspond to computations while no output file is generated. For the purpose of comparison, all simulations ran on the same desktop computer (Intel Core i7 CPU 2.8 GHz), except for the HMS model for medium real-case watershed which ran on Intel Core i7-2600 CPU 3.4 GHz.

As can be seen in Table 1, the computational expenses of those models with CFL condition depend a great deal on the choice of Courant number. Less strict *Cr* reduces the run time, but causes instability observed as flow oscillations. Note that the demonstrated hydrographs in previous sections correspond to stable results. Accordingly, run times of those simulations with stable results are compared with OFS-CA.

It was concluded that OFS-CA is certainly computationally more efficient than HMS. Several factors could contribute to this higher performance: (1) the numerical scheme in OFS-CA is of first order accuracy, whereas the scheme in HMS is of second order accuracy and therefore computationally more expensive; (2) OFS-CA simplifies the computation of flow velocity, while HMS solves the fully dynamic SWE; (3) the MJ-t criterion in OFS-CA is less strict than the CFL stability condition, as described shortly later; (4) OFS-CA reduces the stability constraints; and (5) OFS-CA reduces computational costs by considering the trapped water within ponds as immobile (feature MJ-p).

A comparison with HE2D is more challenging, since the models have different grid geometries. However, in the two cases investigated here, OFS-CA showed a higher performance. This should not only be due to their different stability conditions, i.e., point (3) mentioned above, but also due to feature MJ-p and reducing the stability constraints, i.e., points (4) and (5).

In order to provide a more meaningful comparison between the MJ-t criterion and CFL condition, a further version of OFS-CA was developed, named OFS-CA_{CFL}, in which CFL condition with restricts the size of time steps. OFS-CA and OFS-CA_{CFL} are in all other aspects identical. Table 1 reports the run times of this version for the two real-case watersheds. For the medium watershed, in order to produce hydrographs as smooth as those gained by OFS-CA, the Courant number was limited to a fairly small value (0.02), leading to longer run times compared to OFS-CA. In comparison with HE2D, OFS-CA** _{CFL}** has shorter run times, although both apply CFL condition with . This is probably due to the fact that the distance between the cell centroids are generally smaller in HE2D (unstructured triangular grid), leading to smaller time steps.

In a further step, the size of time steps determined by OFS-CA and OFS-CA** _{CFL}** were compared, showing that MJ-t criterion allows more dynamic sizes (i.e., broader range of ) compared to CFL condition with a constant Courant number. Some statistics of distributions for the two real-case watersheds are summarized in Table 2. OFS-CA resulted in smaller total number of time steps compared to OFS-CA

**. Considering that each time step corresponds to a computational loop, this leads to reduction of run times (Table 2). However, note that not only the number of time steps affects the performance, but also the total number of cells and domain's topography.**

_{CFL}Case 4: Small real-case watershed; Run: 800 s | Case 5: Medium real-case watershed; Run: 3 h | |||
---|---|---|---|---|

OFS-CA | OFS-CA_{CFL}; Cr = 0.1 | OFS-CA | OFS-CA_{CFL}; Cr = 0.02 | |

Minimum | 0.001 | 0.066 | 0.001 | 0.009 |

Maximum | 4.443 | 0.272 | 4.584 | 0.161 |

Average | 0.192 | 0.170 | 0.061 | 0.033 |

Total number of time steps | 4,066 | 4,616 | 176,181 | 323,731 |

Reduction in the total number of time steps (%) | 12% | 46% | ||

Run time | 0.6 s | 0.8 s | 1.8 h | 2.7 h |

Run time decrease (%) | 25% | 33% |

Case 4: Small real-case watershed; Run: 800 s | Case 5: Medium real-case watershed; Run: 3 h | |||
---|---|---|---|---|

OFS-CA | OFS-CA_{CFL}; Cr = 0.1 | OFS-CA | OFS-CA_{CFL}; Cr = 0.02 | |

Minimum | 0.001 | 0.066 | 0.001 | 0.009 |

Maximum | 4.443 | 0.272 | 4.584 | 0.161 |

Average | 0.192 | 0.170 | 0.061 | 0.033 |

Total number of time steps | 4,066 | 4,616 | 176,181 | 323,731 |

Reduction in the total number of time steps (%) | 12% | 46% | ||

Run time | 0.6 s | 0.8 s | 1.8 h | 2.7 h |

Run time decrease (%) | 25% | 33% |

## CONCLUSIONS

The proposed diffusive wave shallow water model, referred to as OFS-CA, showed a promising accuracy and performance in five test cases (1D and 2D problems). Water depth and discharge hydrographs for an inclined open channel matched excellently with those gained by a second order accurate fully shallow water model (HMS).

The wave-front induced by a 1D dam-break in the absence of a bed slope was simulated with good accuracy: flow profiles and water depth hydrographs matched the available experimental data very well, showing the suitability of the feature MJ-v developed in the present study.

For a V-shaped schematic catchment, discharge hydrographs agreed perfectly with an available analytical solution and the results of HMS. Depth hydrographs obtained by OFS-CA and HMS matched excellently at the channel outlets, but differed at hillside edge, where OFS-CA simulated water depths more realistically.

For two 2D real-case watersheds, results were not only compared with HMS but also with the commercial diffusive wave model HE2D (itwh). Overall, water depth and discharge hydrographs agreed reasonably, but two main explainable differences were observed: (1) OFS-CA and HMS differed in lateral water distribution, since in HMS the influence of the state of neighbouring cells on the flow is greater than in OFS-CA due to their different numerical schemes (first order in OFS-CA vs. second order in HMS); (2) in case of surface ponds being filled up with water, OFS-CA takes only the effective flow into account (feature MJ-p); this led to lower flow discharges over filled ponds and also at downstream compared to HMS and HE2D. The pond treatment in feature MJ-p could be further examined in future works using experimental data.

OFS-CA was computationally less expensive than both HMS and HE2D. When compared to HE2D, its higher performance is due to (1) its novel stability condition (called MJ-t) leading to more dynamic time step sizes (i.e., broader range) compared to those restricted to CFL condition with a constant Courant number, (2) reducing the stability constraints by dismissing some of the cells located within ponds from simulation while ponds are filling up, and (3) feature MJ-p which considers some part of the overland flow as being trapped within ponds.

Compared to HMS, the computational gains of OFS-CA are not only due to the three above-mentioned features but also due to (4) neglecting inertia terms of the momentum equation and (5) being of first order accuracy in contrary to HMS which is of second order accuracy and thus an expensive solver (these two reasons also explain the shorter run times of HE2D against HMS). OFS-CA will be enhanced for parallel computations to reduce run times further, as its structure is very suitable for that.

## ACKNOWLEDGEMENTS

This research was carried out towards a PhD at Chair of Water Resources Management and Modeling of Hydrosystems, Technische Universität Berlin. We thank the Liegenschaftsamt (LA) and Tiefbauamt (TBA) of the city of Karlsruhe and the Institute for Technical and Scientific Hydrology (itwh) for providing the elevation data used for test cases 4 and 5. The first author appreciates the financial support provided by itwh and thanks the very helpful people there, especially Przemyslaw Olejniczak and Ulrich Egger (for kindly sharing their knowledge of C + +) and Dr Lothar Fuchs. She also thanks Prof. Alejandro Clausse, Dr Albert Chen, Prof. Paul Bates, and Wasantha Lal for kindly responding to her questions about their publications, and appreciates the constructive comments of the anonymous reviewers.