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
. 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.A central cell and its two basin cells.
Temporal depth variations
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: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.
Effective depth of water for computing outflows.
Accumulation of surface water within ponds (top). An imaginary sequence showing the principle of algorithm MJ-p (bottom).
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,
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:
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
The schematic algorithm of method MJ-t to determine at each computation step.
is given by the time required to distribute a water volume of from reference cell into its basin: 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.
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:
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
, first Equation (4) is rewritten as below since is a constant value in the summation:
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.16
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.
Temporal variation of discharge (left) and water depth (right) at a control point in the open channel.
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).
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
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).
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).
Temporal change of water depth at three locations of the 1D dam-break flume. Experimental data are illustrated in the background.
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
Geometry of the V-shaped catchment. Control cells at channel outlet and on the hillside edge are marked.
Geometry of the V-shaped catchment. Control cells at channel outlet and on the hillside edge are marked.
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.
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
Small real-case watershed: (a) ponds and potential flow paths and (b) water depths 50 s after rain stops simulated by OFS-CA.
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.
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.
Hydrographs at the control points of small real-case watershed (note: in HE2D the grid is transformed).
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
Medium real-case watershed: (a) ponds and potential flow paths and (b) water depths 30 minutes after rain stops, simulated by OFS-CA.
Medium real-case watershed: (a) ponds and potential flow paths and (b) water depths 30 minutes after rain stops, simulated by OFS-CA.
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
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.
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.
%20cropped.png?versionId=7142)