Improvement of anisotropic porosity models with a merging technique

Anisotropic porosity shallow-water models are used to take into account detailed topographic information through porosity parameters multiplying the various terms of the shallow-water equations. A storage porosity is assigned to each cell to reflect the void fraction in the cell and a conveyance porosity is used at each edge to reproduce the impact of subgrid obstacles on the flux terms. To guaranty the numerical stability, the time step depends on the value of the porosity parameters. This may hamper severely the computational efficiency in the presence of cells with low values of storage porosity. Cartesian grids are particularly sensitive to such a case since the meshing stems directly from the choice of the grid size. In this paper, this problem is addressed by using an original merging technique consisting in merging cells with a storage porosity lower than a threshold value with neighbouring cells. The model was tested for modelling a prismatic channel with different orientations between the Cartesian computational grid and the channel direction. The results show that the standard anisotropic porosity model (without merging) improves the reproduction of the flow characteristics; but at the cost of a significantly higher computational time. In contrast, the computational time is drastically reduced and the accuracy preserved when the merging technique is used with the porosity model.


Introduction
Anisotropic porosity models (or integral porosity models) were introduced by [1] for the simulation of floods in urban areas at a coarse scale, using porosity parameters to reproduce the effect of obstacles at a finer scale on water storage and flow conveyance. Recent developments include a generalization of the anisotropic porosity model by distinguishing domain-and boundary-based flow variables (dual integral porosity model) [2], by considering the depth-dependency of the porosity parameters [3,4], by extending the definition of conveyance porosity parameters from purely geometric parameters to make them dependent on the flow and the scale of modelling and by exploring the use of distinct conveyance porosity parameters in the various fluxes of the governing equations [5].
Anisotropic porosity models are generally used with unstructured meshes which enable a gap-conforming mesh for the proper computation of the conveyance porosities [2,[6][7][8].
However, recent research has shown the benefit of such models for computations on Cartesian grids [9,10]. A key issue in the use of porosity-based models on a Cartesian grid is the handling of computational cells with low values of the storage porosity, which leads to a significant increase in the computational time to maintain the numerical stability.
Based on simulations of a uniform rectangular channel with different orientations with respect to a fixed Cartesian grid, the aim of this paper is twofold: (1) showing that anisotropic porosity models improve the reproduction of oblique boundaries on a Cartesian grid and (2) showing the benefit of adding a merging technique to the porosity-based model to handle properly the cells with low values of storage porosity.

Governing equations
For a cell j and its K edges k, the discretized formulation of the governing equations for the shallow water equations with anisotropic porosity writes [1]: with j  the total horizontal surface of cell j and k  the total horizontal length of edge k.
The variables U , fluxes F and sources terms S are given by: with h the water depth, u and v the xand y-velocity components, The storage and conveyance porosity parameters  and  are computed respectively as the void fraction (i.e. the fraction not occupied by an obstacle) within the cell and along the edge.

Numerical model
Equations (1) and (2) are solved with the hydraulic model Wolf2D [5,12] on a Cartesian grid. The stationary water level 0  in a direction xor yof a computational cell is computed as a linear combination of the water levels at the two opposite edges of this cell and orthogonal to the considered direction. The parameters of the linear combination are determined based on energy balance considerations, as detailed in [13]. The numerical stability of explicit numerical schemes for solving the governing equations is controlled by a modified Courant-Friedrichs-Lewy criterion, in which the time step depends on the value of the porosity parameters [1]: with t  the time step, cj the wave velocity in cell j and CFL the Courant number guarantying the stability of the numerical model.

Merging technique
According to Eq. (3), the computational time is expected to increase significantly in the presence of low values of the storage porosity, as a consequence of the reduction of the time step.
A simplistic solution to this problem consists in removing from the computational domain the cells whose storage porosity is lower than a threshold porosity min  . However, this induces an underestimation of both storage and conveyance capacity of the domain (Fig.  1). To address this issue, an original merging technique based on the approach developed by [14] in the context of cutcells is added to the anisotropic porosity model [5]. This merging technique preserves the real storage and conveyance properties while avoiding the explicit computation of the cells with a storage porosity lower than min  . The method is based on the four-step procedure schematized in Fig. 2:

Test case description
We consider a 15 km-long and 120 m-wide rectangular channel with a constant bed slope of 0.2% and a uniform Manning roughness coefficient of 0.025 sm -1/3 . A steady discharge of 2,000 m 3 /s is prescribed at the upstream end. At the downstream end, the normal water depth is prescribed (hu = 2.52 m) so that the computed water depth should theoretically remain the same all along the channel. This is indeed the case when the channel is oriented along the Cartesian grid. In contrast, higher water depths are computed along the channel if the channel is not aligned with the Cartesian grid, due to abrupt changes in the cross section appearing as artefacts of the numerical discretisation of the inclined lateral walls of the channel. Considering a grid size of 30 m, Fig. 3 shows that the computed water depth along the centreline of the channel increases from downstream to upstream and then fluctuates around a "numerical" normal water depth ĥ u . The numerical error on the upstream water depths due to the discretization of the channel sidewalls on a Cartesian grid is quantified by using the ration ĥ u /h.
Porosity parameters are used for the cells along the sidewalls to improve their reproduction in the numerical model. No drag loss term is used here (

Results and discussion
Simulations were performed considering angles α between the orientation of the channel and the Cartesian grid ranging between 0° and 45° and a fixed grid spacing of 30 m x 30 m. Fig.  4 compares the overestimation of the upstream water depth obtained with a classical shallowwater model (black line) with the one resulting from anisotropic porosity models with different values of the threshold porosity and with or without merging. The maximum value of the relative numerical error over all the angles α and the relative computational time are reported in respectively Fig. 5a and Fig. 5b as a function of the value of the threshold porosity.
The numerical normal depth obtained with a classical model (continuous black line) exceeds by almost 40% the theoretical value for angles α in-between 20° and 30° (Fig. 4a). This overestimation is significantly reduced when using the anisotropic porosity model without merging, especially as the value of the threshold porosity is reduced (Fig. 5a -red  color). This highlights the ability of anisotropic porosity models to improve the description of oblique boundaries on a Cartesian grid.
However, the computational efficiency is severely hampered by a reduction of the value of the threshold porosity (Fig. 5b), accordingly to Eq. (3). Using the merging technique to deal with cells with low values of the storage porosity, a good accuracy can be obtained with high values of the threshold porosity ( Fig. 4b and Fig. 5a -blue color). Therefore, adding the merging technique in a porosity-based model enables a reduction of the computational cost for a given accuracy.
For instance, the values of the threshold porosity limiting the maximum numerical error to 5% are ϕmin = 0.05 and ϕmin = 0.75 respectively without or with merging. For such values, while the anisotropic porosity model without merging is around 15 times computationally more expensive than the classical model, the anisotropic porosity model with merging is only twice more expensive.

Conclusion
Anisotropic porosity models enable improving the description of obstacles at a coarse scale by using a storage porosity parameter at the cells to represent the available volume of water within it and a conveyance porosity parameter at the edges to reflect the impact of the obstacles on the flux terms. However, the presence of these porosity parameters in the governing equations impacts the stability criterion. This may lead to a significant increase of the computational cost in the presence of cells with low values of the storage porosity.
In this paper, an anisotropic porosity model is used to improve the description of oblique boundaries of a rectangular and uniform channel on a Cartesian grid. The results show that the porosity parameters enable compensating nearly completely for the abrupt changes in cross-section resulting from the discretisation of the sidewalls. Adding a merging technique to the anisotropic porosity model reduces considerably the computational time for a given accuracy, or improves the accuracy of the result for a given threshold porosity.