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 [L3]

  • total discharge (total outflow) [L3T−1]

  • total inflow [L3T−1]

  • lateral discharge (lateral outflow): from cell i to basin cell j [L3T−1]

  • x-direction discharge [L3T−1]

  • y-direction discharge [L3T−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

The space is discretized with a 2D regular square grid. At each computation step, a 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.
Figure 1

A central cell and its two basin cells.

Figure 1

A central cell and its two basin cells.

Temporal depth variations

The governing equations in diffusive wave models are the simplified form of the SWE. The fully SWE can be displayed as: 
formula
1
According to the conservation of mass, the temporal change of water depth within a cell is given by the difference of inflows and outflows and by including sink and/or source terms . is the vector of discharge per unit of width [L2T−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: 
formula
2

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.

In a first order accurate finite volume scheme, the conservation of mass (see Equation (1)) can be formulated for each cell as below to compute the temporal variation of water depth: 
formula
3
where shows the variation of water depth and is the size of time step. and stand for total inflow [L3T−1] and total outflow [L3T−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: 
formula
4
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 celli, required for Equation (3), is computed by summing up the corresponding outflows of its neighbour cells.
Effective depth of water differs from the actual water depth (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.
Figure 2

Effective depth of water for computing outflows.

Figure 2

Effective depth of water for computing outflows.

Figure 3

Accumulation of surface water within ponds (top). An imaginary sequence showing the principle of algorithm MJ-p (bottom).

Figure 3

Accumulation of surface water within ponds (top). An imaginary sequence showing the principle of algorithm MJ-p (bottom).

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,

where water surface gradient and ground slope are computed using the respective ground elevations (z) and water levels (wl) of the central cell i and the basin cell j: 
formula
5
For group (a) of cells, velocities are computed according to Manning's formula (refer to the diffusive wave approximation of SWE displayed in Equation (2)), while for group (b) of cells, an additional formulation, so-called vg [LT−1], is also considered in order to improve the accuracy of estimated flow velocities: 
formula
6
A good example to show the necessity of including 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: 
formula
 
formula
7

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).

term1 and term2 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 term2, 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

At each computation step, the size of global time step is restricted to meet a certain condition at a 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: 
formula
8
Figure 4

The schematic algorithm of method MJ-t to determine at each computation step.

Figure 4

The schematic algorithm of method MJ-t to determine at each computation step.

As demonstrated in Figure 4 (block 1), the global time step is given by the time required to distribute a water volume of from reference cell into its basin: 
formula
9
where is the total discharge [L3T−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.
The global time step is restricted further so that no negative water depth occurs elsewhere. That means the outflowing water volume at other cells of the domain should be smaller than the available water volume there . At each cell, the outflowing water volume is once computed according to the total discharge specified in Equation (4) and once according to an ‘artificial’ total discharge specified as described later. In summary, two conditions should be fulfilled at each cell: 
formula
10
As displayed in Figure 4 (block 2), in the case of these conditions not being satisfied, the global time step is reduced to fulfil them: 
formula
11

Determining (outVRef)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

To specify the artificial total discharge for any cell , first Equation (4) is rewritten as below since is a constant value in the summation: 
formula
12
As pointed out before, is considered as the depth of available water at cell i. By replacing the lateral effective depths with , the right hand side of Equation (12) can be expanded further: 
formula
13
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 :  
formula
14
Finally, the artificial total discharge is given by placing Equation (14) into Equation (13): 
formula
15

Although the artificial total discharge seems to be a rough approximation, it provided a suitable additional restriction for . It should be noted that application of is restricted to MJ-t criterion (Equations (10) and (11)) and it does not contribute to solving Equation (3).

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: 
    formula
    16
    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-CACFL 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.

Results of OFS-CA and HMS were compared at three control points: x = 50, 80 and 100 m. The two models produced identical depth and discharge hydrographs at all three control points. Figure 5 shows that discharge and water depth increase consistently and reach constant values at a certain time (i.e., time of concentration). By the end of rainfall, flow depth and velocity starts dropping, which is also simulated consistently, agreeing perfectly with the results of HMS.
Figure 5

Temporal variation of discharge (left) and water depth (right) at a control point in the open channel.

Figure 5

Temporal variation of discharge (left) and water depth (right) at a control point in the open channel.

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).

Table 1

Model run times

  OFS-CA HMS HE2D OFS-CACFL 
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 sa     
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 sa 
Case 4: Small real-case; Run: 800 s 0.6 s   20 s   6.3 s   0.8 s 0.6 sb 
Case 5: Medium real-case; Run: 3 h 1.8 h 16.9 h Instable 4 h 2.5 hc Instable 2.7 h 1.1 hc Instable 
  OFS-CA HMS HE2D OFS-CACFL 
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 sa     
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 sa 
Case 4: Small real-case; Run: 800 s 0.6 s   20 s   6.3 s   0.8 s 0.6 sb 
Case 5: Medium real-case; Run: 3 h 1.8 h 16.9 h Instable 4 h 2.5 hc Instable 2.7 h 1.1 hc Instable 

aMinor overshooting at the moment of reaching the maximum discharge.

bSome oscillations at control point 3.

cSome oscillations at control points 1 and 2.

Case 2: dam-break flume

One of the innovations of the present study was formulation of Equation (7) (feature MJ-v). This was verified using the information gained by the experimental study introduced in Aleixo 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 h0 = 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.
Figure 6

Two instances of the 1D dam-break flow. Wave front, simulated by OFS-CA, is compared with the camera images and water levels obtained from an experimental study (labelled as Exp).

Figure 6

Two instances of the 1D dam-break flow. Wave front, simulated by OFS-CA, is compared with the camera images and water levels obtained from an experimental study (labelled as Exp).

Figure 7

Temporal change of water depth at three locations of the 1D dam-break flume. Experimental data are illustrated in the background.

Figure 7

Temporal change of water depth at three locations of the 1D dam-break flume. Experimental data are illustrated in the background.

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

Overton & Brakensiek (1970) presented a kinematic solution of surface runoff in general dimensionless form for a V-shaped watershed. Di Giammarco 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 8

Geometry of the V-shaped catchment. Control cells at channel outlet and on the hillside edge are marked.

Figure 8

Geometry of the V-shaped catchment. Control cells at channel outlet and on the hillside edge are marked.

Figure 9 (top) displays the temporal change of discharge at hillside edge and at the channel outlet. OFS-CA produced smooth discharge hydrographs with excellent agreement with the available analytical solution and those obtained by HMS. The gentle early curve in the rising limb of hydrograph at channel outlet shows the gradual increase of discharge in the channel due to rapidly increasing inflow from the hillsides. As can be seen this is captured very well by OFS-CA. The two models also agreed ideally in capturing the subsidence of flow discharge, conforming to the analytical solution for the hillside (note: for the channel outlet, the analytical solution does not cover the receding limb of the hydrograph).
Figure 9

Temporal variation of discharge (top) and water depth (bottom) at the edge of hillside and at the outlet of channel in the V-shaped catchment.

Figure 9

Temporal variation of discharge (top) and water depth (bottom) at the edge of hillside and at the outlet of channel in the V-shaped catchment.

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

A small watershed with an area of 874 m2 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.
Figure 10

Small real-case watershed: (a) ponds and potential flow paths and (b) water depths 50 s after rain stops simulated by OFS-CA.

Figure 10

Small real-case watershed: (a) ponds and potential flow paths and (b) water depths 50 s after rain stops simulated by OFS-CA.

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.

For the purpose of benchmarking, the domain was modelled also with HMS using the same grid as in OFS-CA. Results were inspected at four control points marked in Figure 10. Each control point represents a special case (see their ground elevations in Figure 11), allowing a detailed comparison between OFS-CA and HMS:
  • 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.

Figure 11

Hydrographs at the control points of small real-case watershed (note: in HE2D the grid is transformed).

Figure 11

Hydrographs at the control points of small real-case watershed (note: in HE2D the grid is transformed).

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 m2, varying in the range of 0.1–1.5 m2. 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

A larger area was modelled to assess the performance and results of OFS-CA further. This domain, with an area of about 7 ha, was modelled on a 1 m resolution grid with 216 by 322 cells. For the purpose of comparison, the domain was modelled also with HMS and HE2D. HMS is compatible with both unstructured triangular and regular square grids. This allowed the same grid as in OFS-CA to be used. However, HE2D accepts only unstructured triangular grids. Therefore, the regular grid was transformed into an unstructured triangular grid by keeping the smoothing of elevation data to minimum. This consists of almost the same number of cells as in the regular grid, with an average cell area of 0.95 m2, varying in the range of 0.25–1.95 m2. 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).
Figure 12

Medium real-case watershed: (a) ponds and potential flow paths and (b) water depths 30 minutes after rain stops, simulated by OFS-CA.

Figure 12

Medium real-case watershed: (a) ponds and potential flow paths and (b) water depths 30 minutes after rain stops, simulated by OFS-CA.

The initially dry domain was imposed to a constant rainfall with intensity of 1.3 × 10−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

Hydrographs at the control points of the medium real-case watershed.

Figure 13

Hydrographs at the control points of the medium real-case watershed.

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

This section provides a comparison between the performance of HMS and HE2D with OFS-CA according to the previous test cases. For ease of following the arguments, we briefly refer to the determination of time steps using the CFL stability condition (applied in HMS and HE2D): 
formula
where and stand for the magnitude of velocity and the distance between cell centroids, respectively. Some shallow water models like HE2D consider only velocity , while others like HMS include celerity as well: . Although the non-dimensional Courant number (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-CACFL, in which CFL condition with restricts the size of time steps. OFS-CA and OFS-CACFL 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-CACFL 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-CACFL 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-CACFL. 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.

Table 2

Comparison of the size of time steps

  Case 4: Small real-case watershed; Run: 800 s
 
Case 5: Medium real-case watershed; Run: 3 h
 
OFS-CA OFS-CACFL; Cr = 0.1 OFS-CA OFS-CACFL; 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-CACFL; Cr = 0.1 OFS-CA OFS-CACFL; 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.

REFERENCES

REFERENCES
Aleixo
R.
Spinewine
B.
Soares-Frazão
S.
Zech
Y.
2009
Nonintrusive measurements of water surface and velocity profiles in a dam-break flow
. In:
Proceedings of the 33rd IAHR Conference
,
Vancouver
, pp.
6898
6905
.
Aleixo
R.
Soares-Frazão
S.
Zech
Y.
2011
Velocity-field measurements in a dam-break flow using a PTV Voronoï imaging technique
.
Experiments in Fluids
50
,
1633
1649
.
Bradbrook
K. F.
Lane
S. N.
Waller
S. G.
Bates
P. D.
2004
Two-dimensional diffusion wave modelling of flood inundation using a simplified channel representation
.
International Journal of River Basin Management
2
(
3
),
211
223
.
Costabile
P.
Costanzo
C.
Macchione
F.
2012
Comparative analysis of overland flow models using finite volume schemes
.
Journal of Hydroinformatics
14
(
1
),
122
135
.
Crossley
A. J.
Wright
N. G.
Whitlow
C. D.
2003
Local time stepping for modeling open channel flows
.
Journal of Hydraulic Engineering – ASCE
129
,
455
462
.
Di Giammarco
P.
Todini
E.
Lamberti
P.
1996
A conservative finite elements approach to overland flow: the control volume finite element formulation
.
Journal of Hydrology
175
(
1–4
),
267
291
.
Ettrich
N.
Steiner
K.
Thomas
M.
Rothe
R.
2005
Surface models for coupled modeling of runoff and sewer flow in urban areas
.
Water Science and Technology
52
(
5
),
25
33
.
Fernández-Pato
J.
García-Navarro
P.
2016
2D Zero-inertia model for solution of overland flow problems in flexible meshes
.
Journal of Hydrologic Engineering
21
(
11
),
04016038
.
Ghimire
B.
Chen
A. S.
Guidolin
M.
Keedwell
E. C.
Djordjević
S.
Savić
D. A.
2013
Formulation of fast 2D urban pluvial flood model using cellular automata approach
.
Journal of Hydroinformatics
15
(
3
),
676
686
.
Hinkelmann
R.
Liang
Q.
Aizinger
V.
Dawson
C.
2015
Robust shallow water models
.
Environmental Earth Sciences
74
,
7273
7274
.
Hirsch
C.
2007
Numerical Computation of Internal and External Flows
,
2nd edn
.
Butterworth-Heinemann
,
Oxford
.
Hunter
N. M.
Horritt
M. S.
Bates
P. D.
Wilson
M. D.
Werner
M. G. F.
2005
An adaptive time step solution for raster-based storage cell modelling of floodplain inundation
.
Advances in Water Resources
28
(
9
),
975
991
.
Hunter
N. M.
Bates
P. D.
Neelz
S.
Pender
G.
Villanueva
I.
Wright
N. G.
Liang
D.
Falconer
R. A.
Lin
B.
Waller
S.
Crossley
A. J.
Mason
D. C.
2008
Benchmarking 2D hydraulic models for urban flooding
.
Proceedings of the Institution of Civil Engineers – Water Management
161
,
13
30
.
Lal
A. M. W.
1998
Weighted implicit finite-volume model for overland flow
.
Journal of Hydraulic Engineering
124
(
4
),
342
349
.
Lal
A. M. W.
Toth
G.
2013
Implicit TVDLF methods for diffusion and kinematic flows
.
Journal of Hydraulic Engineering
139
(
9
),
974
983
.
Lauber
G.
Hager
W. H.
1998
Experiments to dambreak wave: horizontal channel
.
Journal of Hydraulic Research
36
(
3
),
291
307
.
López-Barrera
D.
García-Navarro
P.
Brufau
P.
Burguete
J.
2012
Diffusive-wave based hydrologic-hydraulic model with sediment transport. I: model development
.
Journal of Hydrologic Engineering
17
(
10
),
1093
1104
.
Mendicino
G.
Pedace
J.
Senatore
A.
2015
Stability of an overland flow scheme in the framework of a fully coupled eco-hydrological model based on the macroscopic cellular automata approach
.
Communications in Nonlinear Science and Numerical Simulations
21
,
128
146
.
Overton
D.
Brakensiek
D.
1970
A kinematic model of surface runoff response
. In:
Proceedings Wellington Symposium
.
UNESCO/IAHS
,
Paris
, pp.
100
112
.
Parsons
J. A.
Fonstad
M. A.
2007
A cellular automata model of surface water flow
.
Hydrological Processes
21
(
16
),
2189
2195
.
Rinaldi
P.
Dalponte
D.
Vénere
M.
Clausse
A.
2012
Graph-based cellular automata for simulation of surface flows in large plains
.
Asian Journal of Applied Science
5
(
4
),
224
231
.
Ritter
A.
1892
Die Fortpflanzung der Wasserwellen (The propagation of water waves)
.
Zeitschrift Verein Deutscher Ingenieure
36
(
33
),
947
954
.
Shiffman
D.
2012
The Nature of Code – Chapter 7. Cellular Automata [online]. Available from
: http://natureofcode.com/book/chapter-7-cellular-automata
(accessed 15 February 2014)
.
Simons
F.
Busse
T.
Hou
J.
Özgen
I.
Hinkelmann
R.
2014
A model for overland flow and associated processes within the Hydroinformatics Modelling System
.
Journal of Hydroinformatics
16
(
2
),
375
391
.
Thomas
R.
Nicholas
A. P.
2002
Simulation of braided flow using a new cellular routing scheme
.
Geomorphology
43
,
179
195
.
Wang
Y.
Liang
Q.
Kesserwani
G.
Hall
J. W.
2011
A positivity-preserving zero-inertia model for flood simulation
.
Computers &. Fluids
46
(
1
),
505
511
.
Weimar
J. R.
1997
Simulation with Cellular Automata
.
Logos-Verl
,
Berlin
.
Wolfram
S.
2002
A New Kind of Science
.
Wolfram Media
,
Champaign, IL
.
Available from
: http://www.wolframscience.com
(accessed 15 February 2014)
.