Haar wavelet-based adaptive finite volume shallow water solver

Within computational hydraulics, the finite volume (FV) Godunov-type method has been extensively used for modelling shallow water flows due to its desirable properties, such as locality, numerical mass conservation and ability to inherently incorporate discontinuous flow transitions within the numerical solution (Toro 2001; Bouchut 2007; Toro & Garcia-Navarro 2007).

The original Godunov method for the gas dynamic equations assumes that the representation of the numerical solution over each discrete control volume is piecewise-constant. Connecting these approximate solutions via fluxes obtained from the solution of Riemann problems between cells leads to make the scheme well suited for solving nonlinear equations containing discontinuities (Harten et al. 1983). Further theoretical and numerical considerations have been given to apply Godunov-type FV methods to shallow water flow problems in order to properly incorporate source terms (mainly the bed slope term) and model evolving wet/dry fronts (LeVeque & George 2008; Liang & Borthwick 2009; Liang & Marche 2009; Kesserwani 2013). Following decades of research, FV Godunov-type methods have become widely applied to simulate real-scale flooding and have been adopted into commercial hydraulic modelling software packages such as TUFLOW-FV and RiverFlow2D PLUS. Nevertheless, real large-scale shallow flows have complex flow features such as shocks, contact discontinuities and a wide range of spatial scales. Typically, the computational domain is discretized uniformly using a large number of cells, given that the position of flow features is usually unknown and capturing certain small scales within a coarse mesh simulation may be difficult without causing computational cost trade-off. Therefore, automated mesh adaptation comes in handy to improve modelling efficiency and capture the various physical scales involved in shallow water flows.

Various adaptive techniques have been developed within the FV framework applied to solve shallow water equations (SWE). These include moving mesh methods (Skoula et al. 2006) or static grids with local refinement methods (Nikolos & Delis 2009; Caviedes-Voullième et al. 2012). However, most of the present techniques to date are achieved over patches of grid refinements, in a decoupled manner, which controversially gives rise to many problematic effects (Nemec & Aftosmis 2007; Liang et al. 2008; Kesserwani & Liang 2012a, 2015; Mocz et al. 2013). For instance, they require error sensors and multiple user-chosen parameters (e.g., for setting up grid resolution coarsening vs. refinement), which introduce sensitivity (e.g., can lead to inadequate or excessive resolution), inflexibility and problem dependency (e.g., due to the need to tune many parameters for each simulation problem). They also lack a rigorous strategy to accommodate flow data transfer and recovery between various inner resolution scales (given the changing nature of the mesh). Therefore, the design of a more integrated adaptation approach, which can address such problems, is still a crucial challenge (Hovhannisyan et al. 2014) and is the aim of this work.

The classical wavelet theory for local decomposition/reconstruction of (self-similar) functions (by translation and dilatation of a single mother function, or wavelet, at finer resolution) offers a natural way for adaptive compression of data. Mathematicians and engineers have widely used the multiresolution capabilities of wavelets in many applications (e.g., signal processing and denoising) (Gargour et al. 2009) including the solution of partial differential equations. The true power of such a multiresolution approach arises from the fact that the local piecewise-constant data can be described as a set of scaling coefficients encapsulating higher-resolution details thereby allowing (Harten 1995: Keinert 2003): (a) genuine information exchange between those (heterogeneous) elements with matching resolution (promotion and demotion of the details via application of high- and low-pass filters); (b) achievement of grid resolution adaptivity by selecting certain details from the local compression dataset; and (c) quantitatively control the variation of the adaptive mesh solution with reference to the underlying uniform mesh solution relevant to the finest resolution.

In the context of the FV Godunov-type modelling, a couple of papers have successfully integrated the wavelet theory for adaptive solution of homogenous conservation law (scalar and system) (Harten 1995; Müller 2003). However, the implementation and implication of this idea in addressing practical aspects of shallow water flow simulation is unexplored as yet. To fulfil this gap, the Haar wavelet basis is used to reformulate a FV Godunov-type method to obtain a new adaptive multiresolution scheme, which will be referred to as Haar wavelets finite volume (HWFV). The technical development of the HWFV solving the one-dimensional (1D) SWE with source term and wetting and drying is described. Particular focus is mainly placed on how wavelet-based adaptivity is achieved within the HWFV SWE numerical solver. Selected test cases are used to systematically verify the performance of the HWFV scheme addressing issues of adaptivity parameterizations, modelling of irregular topography with/without wetting and drying, accuracy preservation and mesh convergence.

The rest of the paper is organized as follows: a brief overview of 1D SWE is first presented, followed by a section describing the baseline FV method that will be used. Second, the multiresolution analysis (MRA) and its mathematical properties are introduced, with particular focus on how the Haar wavelets basis and its scaling basis are incorporated into the FV method. Finally, the performance of the HWFV model is tested, analysed and discussed, and conclusions are drawn.

The techniques adopted for well-balanced discretization of the bed slope source term with wetting and drying are well-established and verified. For more details see Liang & Marche (2009), Kesserwani & Liang (2012b) and Kesserwani (2013). Herein, the threshold for dry cell definition is fixed to ⁠. The stability condition on the time step size (Cockburn & Shu 2001) is controlled by a Courant–Friedrich–Lewy condition, which is chosen to be equal to 0.98 for simulations involving fully wet domains and 0.5 when wet/dry fronts are present.

In the expansion (Equation (11)), the representation of the function f at higher resolution (n > 0) is enabled by the detail coefficients when they are non-zero.

By recursive application of the two-scale transform equations, the expansion of scaling coefficients describing a function f can be demoted or promoted across any different resolutions (Figure 3). Consequently, promoting and demoting the coefficients across different resolutions, i.e., by using Equations (15) and (16), is the keystone of the HWFV scheme.

This section shows how to exploit the multiresolution description of function f (recall the section ‘MRA’) to decompose the local piecewise-constant solution (recall the section ‘Overview of the FV Godunov type-framework’) as a compressed data set of the solution's information, which drives adaptivity. The strategy is to first generate a FV discretization at the highest resolution (herein presented for n = 2 without loss of generality) to produce a fine uniform reference mesh. Then, over this mesh, MRA is performed to select up to which resolution the local numerical solution needs to be described (i.e., form a heterogeneous grid for actual FV calculations). In doing so, three main steps are needed: prediction for careful mesh refinement, multiresolution FV update and thresholding step for mesh coarsening through solution decompression.

Here, the translating index j is related to the resolution level of the local cell n.

Equations (17) and (18) need to be only applied at a local level resolution (i.e., on the adaptive mesh). This mesh can be selected by, first, casting each local solution as in Equation (11). This requires demoting high-resolution data in order to store detail coefficients and produce the scaling coefficient at level (0). Then, a mesh prediction step is needed to locally promote the solution and refine the grid. This will yield the adaptive mesh on which calculation is performed. Finally, hard thresholding is performed to reduce the density of the mesh by discarding all non-significant details that fall below a certain threshold value (details in the sections below).

The components in are associated with the components of the conserved variables in (i.e., h and q). In real computations, it is impossible to know the optimal threshold value but a range of options is feasible (Gerhard & Müller 2013; Hovhannisyan et al. 2014) as will be investigated in the following section. By default, in this work, the threshold value (at the coarsest level) is set to and it is normalized according to the resolution levels. In addition, the magnitudes of the detail coefficients are scaled according to the maximum value of the numerical solution.

The topography is projected into the highest resolution (n = 2) in the same way as the conserved variables. Thus the compressed data set of bed information over each cell is obtained. Finally, the adaptive bed mesh is obtained by applying the thresholding step. However, the topography adaptation is performed once (i.e., initially when ⁠) and it remains constant throughout the simulation time. To ensure the preservation of an accurate water surface elevation and mass conservation across levels, the bed mesh should be considered as a reference; this means that the HWFV scheme does not allow demoting the local numerical solution to coarsen the mesh to a resolution level lower than the topography refinement, even if the local numerical solution allows for a lower refinement level.

This section shows the ability of the new HWFV scheme to solve 1D shallow water problems with fewer cells when comparing to the reference FV scheme, while retaining accuracy, total mass in the system, well-balancing and positivity of water depth. Five well-known benchmark problems are considered for validation of the HWFV scheme. Furthermore, the root mean square error (RMSE) and the relative mass error (RME) are computed. Two of the test cases are idealized dam-break cases over frictionless beds (regular and irregular), the third and fourth test cases consider steady flow over a hump (transcritical flow with shock and supercritical flow) and the fifth test case considers a deeper analysis on the test case of the oscillatory flow in a parabolic bowl.

The purpose of this case is to test the capability of the HWFV scheme to efficiently and accurately solve the homogenous SWE. The solution is compared to the exact solution, using the RMSE and the maximum error metrics. A 1D channel with a horizontal frictionless bed is considered. The length of the channel is 2,000 m and an imaginary dam is located at 1,000 m from the upstream end. Initially, the upstream water level is 20 m whereas at the downstream end, the water level is 5.0 m. The initial discharge was set to zero in every cell. Boundary conditions, although set to be numerically transmissive, are effectively irrelevant in this case as the propagating wave does not reach the boundary. At the instant of dam failure, water is released producing a shock wave travelling downstream meanwhile a rarefaction wave is formed propagating upstream. The computational domain at the coarse level (n = 0) is discretized with 71 uniform cells and the computational model was run up to 40 s after the dam-break. Since the HWFV model allows for up to three levels of mesh refinement, the size of the reference fine uniform mesh is 568 cells. The numerical results are illustrated in Figure 6. The highest level of resolution is noted to be reached at the shock wave and the kink at the tail of the rarefaction wave. The other zones of the rarefaction wave are achieved with the intermediate resolution levels 1 and 2. Meanwhile, at the rest of the domain, where the solution is smooth, the HWFV has retained the baseline coarse level. In this test, the adaptive solution required a maximum of 124 cells (22% of the fine reference mesh).

In the following test cases, the highest level will be fixed to n = 2 and the focus will be placed on exploring applied aspects that are relevant to hydraulic modelling.

The initial water surface and flow rate per unit width are set to h + z = 0.33 m and q = 0.18 m²/s, respectively. Physical boundary conditions consisted of the steady discharge at the inflow and the initial water level at the outflow. The coarse baseline mesh comprised 50 uniform cells and simulations are noted to converge at around t = 170 s. For this test, the adaptive HWFV model has been further implemented and assessed along with the upwind discretization of the source term (Vázquez-Cendón 1999). Figure 9 presents the numerical results of the adaptive scheme (with n = 2) as compared to the analytical solutions. The numerical water surface profile shows a very good agreement with the analytical solution for both the upwind and the standard cell-centred source term discretization. In terms of prediction to the flow rate, it seems to be typically improved by the upwind source term discretization apart from the peak caused by the presence of the water jump. Arguably, this is a known deficiency in standard FV schemes and, thus, does not relate to the proposed wavelet-based adaptivity (Garcia-Navarro & Vazquez-Cendon 2000; Toro & Garcia-Navarro 2007)). In this case, the majority of the domain features required coarse resolution level except at the hump, which dictated local level 1 of refinement from the onset and at the discontinuities (i.e., starting kink of the transcritical transition and shock) where the level was refined to highest.

To ensure that the disturbances observed in the discharge in the previous test case (Figure 9) are not induced by adaptivity over the local cells (see the section ‘MRA’), the supercritical flow case over a hump is also considered. Channel geometry and bed topography are identical to the previous test case, but the unit inflow rate and elevation water surface at the upstream of the channel are set to h + z = 2 m and q = 25.0567 m²/s, respectively. Herein, both of these physical values are used as steady state inflow boundary conditions, whereas a free outlet is numerically set. The simulation time is 10 s and 50 cells are chosen for the baseline coarse mesh. The results are shown in Figure 10. The constant flow rate and surface water profile compared with the analytical solution are well captured. Again, the mesh refinement was obtained only in the regions where both bed elevation and flow are varying. However, no artefacts are noted in the prediction of the discharge for both the upwind and cell-centred discretization. Consequently, the two test cases demonstrate the good performance of the adaptive shallow water flow model in reproducing well-balanced steady flows over topography, and in performing selection of resolution levels in relevance with flow and topographic regions.

This test case is employed to verify the capability of HWFV to preserve the total mass in the system in order to ensure that the adaptivity process does not introduce or lose mass even in the presence of shocks and wet/dry fronts. A hypothetical dam-break over a triangular hump test is set up (see Figure 11).

In Figure 13, the results show that the RME profiles for both the adaptive mesh levels and their uniform mesh counterparts are the same. This means the HWFV is not introducing any additional mass error beyond the capability of the discretization relative to the fine reference uniform scheme.

This case is well known and recognized as a challenging test case for numerical models because it involves both moving wet/dry interfaces and it has an uneven topography. It is selected here to study accuracy and mesh convergence abilities of the adaptive HWFV scheme and to explore the sensitivity of the HWFV model, i.e., in performing the adaptivity process, considering various choices for the threshold value parameter and for the baseline coarsest mesh.

To understand the effect of the threshold value parameter on the adaptivity process, a baseline mesh with = 40 cells is fixed, while considering the following threshold values: = 0.0, 0.001, 0.01, 0.1 and 1.0. In Figure 14, the numerical water surface profiles at t = 275 s (T/5) and t = 2,020 s (1.5 T) are compared with the analytical solution. As seen Figure 14, the HWFV scheme refines in the region where the wet/dry interface is moving through. Meanwhile, other parts of the domain stay at the coarsest and the intermediate levels of refinement. Furthermore, it can be seen that varying the threshold value leads to different refinement patterns. In particular, it is shown that, as increases fewer cells are refined to higher levels during the simulation, since smaller detail coefficients are selectively omitted. In Figure 14(a), it is clear that the adaptive HWFV scheme activates all the detail coefficients when = 0.0 and all the computational cells go to the highest level, thus resulting in a uniform mesh. Figure 14(b) shows an almost uniform mesh prediction when = 0.001. Hence, it is not reasonable, in terms of efficiency, to use a threshold value < 0.001. A sensible refinement is obtained with = 0.01 (Figure 14(c)). For = 0.1 and = 1.0 (shown in Figure 14(d) and 14(e), respectively) poor predictions at the wet/dry interfaces appear, which may be caused by the fact that a large threshold value leads to omitting some small detail coefficients, relevant to the wet/dry front, during the promotion process.

The choice of a baseline mesh with = 40 is rather arbitrary, and it might be already too fine, or on the contrary, never allow for a fine enough mesh for this test case. This may affect the performance of the adaptivity process. Therefore, the influence of the baseline mesh should be studied. Several baseline meshes at coarse level = 20, 40, 80, 160 and 320 are introduced to address this. The same settings of the threshold value as reported in Figure 14 are used. The evolution of the number of active cells is presented in Figure 15. The results confirm that all considered combinations of and are able to perform the adaptive solutions. Moreover, when refining the baseline mesh, the HWFV requires a reduction of the threshold value to better perform adaptivity process. This is due to the fact that most of the flow region results in rather smooth solutions; therefore the value of the detail coefficients is small. In Figure 15, reduces as the baseline mesh is refined, regardless of the varying threshold. However, in Figure 15(d) and 15(e) for = 0.1 and = 1.0, the magnitude of is relatively the same and with values bounded between 1.0 and around 2.5. These values are less when compared to other threshold values, but they strongly influence the quality of the numerical solution. Notably, with = 0.01, regardless of ⁠, optimal results are obtained, in comparison to other threshold values. The value of is bounded between around 1 to 4. This case shows a particular trend of how the pattern of active cells varies with ⁠. This trend indicates that = 0.01 is the most sensitive to ⁠. Owing to this sensitivity, this threshold value allows for a wide, automatic response of the adaptive process (contrary to, for example = 0.1) and therefore is likely to be the best to perform adaptivity in a prompt way. It is clear that an inefficient adaptive process is obtained with = 0.001 except in Figure 15(e), when a very fine baseline mesh is used.

Mesh convergence is also studied in terms of the -norm as defined in Equation (28). Simulations were performed up to 2,020 s and for analysis the t = 270 s (T/5) and t = 2,020 s (1.5 T) are selected. Several baseline meshes at coarse level = 10, 20, 30, 40, 80, 160 and 320 are introduced with using the same threshold values as defined previously. It can be seen in Figure 16 that the behaviour of the -norm is asymptotic regardless of and ⁠. The uniform mesh (⁠ = 0.0) is taken as a reference curve for comparison. Each point within a curve is associated with each baseline mesh. The difference between the convergence curve for = 0.001 and the reference are relatively small. For = 0.01 the convergence is slightly better than all as shown in Figure 16(a) and 16(b). The large error obtained with = 1.0 and = 20, 30 cells and this is due to too few cells being activated during the adaptivity process (over-filtering). Thus the quality of the numerical solution is affected. Nevertheless, the magnitude of the error becomes smaller as increases but still remains larger than the error of the reference curve. Furthermore, the better results are obtained for = 0.1 but with sensible differences compared to = 1.0. The magnitude of the -norm is increased from Figure 16(a) to 16(b). This is merely because of numerical diffusion which is an anticipated issue for FV first-order schemes. In Figure 17, the same analysis as the one reported in Figure 16 is considered to illustrate the relative performance of CPU time (RPCPU). Figure 17(a) shows the ratio of the CPU times of the adaptive schemes to the CPU time obtained from the associated fine uniform reference schemes (i.e., = 0.0); while Figure 17(b) further describes the normalized CPU times, which are obtained by dividing all CPU time by the maximum one (i.e., N₀ = 320, = 0.0). The results show clearly that when > 0.0, less time is required for achieving a simulation as the baseline mesh N₀ density increases. They also show that the efficiency of the adaptive HWFV schemes is near their equivalent uniform mesh FV schemes when the baseline mesh has a size N₀ ≤ 40 and despite the choice of ⁠.

Adaptive mesh refinement schemes are useful tools to efficiently model various scales of shallow water flow. A new adaptive formulation is proposed that combines the Haar wavelets with the FV method (HWFV). The appeal of the formulation can be easily exploited to drive spatial resolution adaptation from the solution itself according to a single threshold value ɛ. A series of numerical tests has been performed, considering issues of well-balanced property, convergence of scheme, sensitivity relating to different choice of the thresholds values and ability to treat wet/dry fronts.

Numerical evidence confirms that the proposed HWFV scheme can accurately solve the SWE with source term. Adaptive solutions are shown to be mass conservative. Notably, the new model is proven to be able to self-decide appropriate resolution levels following the dynamics of the flow, including shocks, wet/dry fronts and the presence of non-flat topography. Future research will consist of integrating the friction source term in a well-balanced manner, exploring sensitivity issues in relation to increasing the resolution levels with varying threshold values and extending the approach to 2D.

The first author acknowledges the support of the Kurdistan Regional Government, Ministry of High Education and Scientific Research. The research is further supported by the UK Engineering and Physical Sciences Research Council (grant ID: EP/K031023/1) and by the Pennine Water Group platform grant (grant ID: EP/I029346/1).