A consistent finite difference local inertial model for shallow water simulation

: Hydrological-hydraulic modeling is a core technique in assessing surface water dynamics of rivers, lakes, and floodplains. The local inertial model (LIM) as a physically simplified model of the shallow water equations is essential for efficient numerical simulator of surface water dynam‐ ics. In this paper, we point out that the conventional semi-implicit finite difference scheme for the friction slope terms, despite being convenient, is not consistent in the sense that it may lead to incorrect numerical solutions if the temporal resolution is not high. We propose an alternative discretization to resolve this issue, which is more accurate and stable, and has comparable computational efficiency. The new numerical scheme is implemented into a modern hydrological-hydraulic model, demonstrating reasonable accuracy. The new scheme is also compared with a recently-proposed implicit scheme, demonstrating compa‐ rable theoretical and computational performances. The results indicate that the proposed scheme potentially serves as a new central core for numerical simulation with the LIM.


INTRODUCTION
The local inertial model (LIM) (Bates et al., 2010) as a simplified counterpart of the conventional shallow water equations has recently been a core of efficient numerical solvers of shallow surface water dynamics (Dottori et al., 2016;Yamazaki et al., 2013;Tanaka et al., 2018a). On the horizontally 2-D x -y plane, the LIM is given by the system: where t is the time, h = h(t, x, y) ≥ 0 is the water depth, p = p(t, x, y) and q = q (t, x, y) are the line-discharges in the xand y-directions, z = z(x, y) is the bed elevation, s is the source term, and g is the gravitational acceleration. The conventional Manning formula with the roughness coefficient n is set as the friction slope terms. A staggered finite difference scheme (Bates et al., 2010) has been widely used for discretization of the LIM. It has been stated that the friction slope terms in (2) and (3) should be handled in a semi-implicit manner for numerical stability (Bates et al., 2010;de Almeida and Bates, 2013). However, recently, Tanaka and Yoshioka (2017) found that this discretization becomes unstable even when the classical Courant-Friedrichs-Lewy (CFL) condition is satisfied. An implicit scheme was developed as an alternative to the semi-implicit one (Tanaka and Yoshioka, 2017). The semiimplicit discretization has been found to be physically inconsistent when the temporal resolution is not high, while the implicit one does not encounter this issue (Tanaka et al., 2018b). Similar implicit schemes have been proposed and examined for the shallow water equations (Hou et al., 2018;Xia and Liang, 2018). The terms "implicit" and/or "semi-implicit" apply not to the whole system of equations but only their momentum equation. Several methodologies for discretization of the friction slope terms have been proposed so far as reviewed above. The governing equations of the fluid flows are essentially based on continuum physical laws, namely some differential equations. Therefore, we consider that utilizing local exact solutions to discretization may improve performance of numerical schemes. Here, the term "local exact solutions" means exact solutions to some exactly-solvable reduced problems like spatially-discretized problems. We apply this idea to the friction slope terms. Then, is it possible to develop a comparably accurate, stable, efficient, and physically more reasonable scheme? This is the motivation of our paper, and the issue is tackled with the help of the exact solution-based discretization method (Li et al., 2010) in which nonlinear source terms of partial differential equations are discretized based on the knowledge of locally exact solutions. Since it is based on (local) exact solutions, the numerical solutions are expected to inherit several properties of true solutions in some way.
The objective of this paper is to present a new discretization scheme for the friction slope terms, which is based on a local exact solution of the LIM. The implicitness of our scheme comes from an integration of the local physical law (Equation 7). We show that the new scheme does not encounter the above-mentioned consistency issue that the semi-implicit encounters. Theoretically, the new scheme is comparably efficient with the semi-implicit and fullyimplicit schemes as it does not require iterative solvers despite the fact that the nonlinear equation has to be solved at each grid point and time step. The results of numerical stability analysis demonstrate that the present scheme can be operated under a less restrictive condition than the CFL condition. Its application to a test problem demonstrates higher accuracy than the semi-implicit scheme. Finally, the new scheme is implemented in the integrated hydrologicalhydraulic model (Tanaka et al., 2018a). The implemented model is then applied to a numerical simulation of seasonal flooding and drying dynamics of Tonle Sap Lake, Cambodia, to demonstrate its versatility in practice in terms of accuracy and computational efficiency.

Numerical scheme
The proposed scheme is presented following the notation of Tanaka et al. (2018a). The quantity Q evaluated at the time t = kΔt at the location (x, y) = (iΔx, jΔy) is denoted as Q i, j k , where i, j, and k are integers or half integers numbering the locations and the time steps. Equation (1) is discretized as Set the positivity-preserving water depths at time k at (i + 1/2, j) and (i, j + 1/2) as The idea of new scheme is simple: rewrite Equation (2) at each ((i + 1/2)Δx, jΔy) as the semi-discrete form and then directly solve it to find p i + 1/2, j k + 1/2 employing p i + 1/2, j k − 1/2 as the initial condition. Note that p i+1/2 is now continuous along t ((k -1/2)Δt ≤ t ≤ (k + 1/2)Δt).
A new idea proposed in this paper is to exactly solve the nonlinear differential equation in time. In fact, our previous research (Tanaka et al., 2018b) used a scheme where the quadratic nonlinearity of the friction slopes is linearized. What they used was an exact solution to a linearized spatially semi-discretized equation, but not a nonlinear one as used in this paper. This point is the difference between their and our discretization methods. Set Equation (7) is non-linear, but is exactly solvable and evaluated at t = k + 1/2 Δt as follows. If A i+1/2, j, k = 0, where the sub-scripts have been omitted for the sake of simplicity of description. Similarly, if A i+1/2, j, k > 0 and Δt . (14) If A i+1/2, j, k > 0 and p i + 1/2, j k − 1/2 > 0 and B i+1/2, j, k Δt >E i+1/2, j, k , If A i+1/2, j, k < 0, then apply the transformations (3) is discretized in the same way. The semiimplicit and fully-implicit schemes are not presented here,

NUMERICALLY CONSISTENT LOCAL INERTIAL MODEL
but can be found in the literature (Tanaka et al., 2018a(Tanaka et al., , 2018b. The only difference between the present and existing schemes are the discretization of Equations (2) and (3).

Consistency and stability tests
Among the present, semi-implicit, and implicit schemes, the present one is the most physically reasonable since its discretization is directly based on the exact solution to the differential equation. In addition, it is based on the time integration of the momentum equations, which avoids the issue of implicitness that the implicit scheme encounters. Its other theoretical advantage is the consistency (Xia and Liang et al., 2018) in the sense that, for large Δt, we asymptotically have the uniform flow limit where the sign of p i + 1/2, j k + 1/2 is the same as that of the water surface gradient. In the present scheme, the same applies to q i, j + 1/2 k + 1/2 . This consistency property, which has been said to be asymptotic-preserving in computational physics (Jin, 1999), implies that a scheme equipped with the property has smaller risk of generating physically unrealistic water flows than schemes without it. In addition, this property guarantees that the present scheme gives a physically consistent numerical solution to the conventional diffusion wave model (Bates et al., 2010) when the temporal partial differential terms of Equations (2) and (3) are small. It is remarkable to note that the implicit scheme satisfies the same property, while the widely-used semi-implicit schemes do not (as shown in Equation (29) in Tanaka et al. (2018b)). Stability analysis against a perturbed uniform flow is carried out to quantify stability conditions of the three schemes. The flow considered here is uni-directional, and the 1-D counterpart of the LIM (q = 0 and ∂h ∂ y = ∂z ∂ y = 0) is employed (Tanaka and Yoshioka, 2017). The initial condition is a uniform flow on a slope having a constant gradient. The Manning's roughness coefficient is set as 0.03 m -1/3 s. The perturbation has a magnitude of 0.01% of the initial water depth and line discharge to both the lower and upper boundaries. A scheme is judged to be stable if it recovers the original uniform flow as time elapses with the relative errors of both discharge and water depth less than 0.1% from those of the uniform flow. Otherwise (i.e. if water depth becomes negative or diverges), the scheme is judged to be unstable. The time step that serves as the threshold of stability is referred to as the maximum allowable time step (MATS). A larger MATS means higher stability. The MATS is computed against the range of spatial resolution from 100 m to 10,000 m. The initial water depth is set as 0.1 m. Figure 1 shows the MATS for the semiimplicit, fully-implicit, and proposed schemes, demonstrat-ing that the semi-implicit scheme has narrower stability than that designated by the CFL condition, while the others have higher stability than the CFL condition. The fullyimplicit scheme has higher stability than the proposed scheme, and their difference becomes larger as the spatial resolution decreases, i.e. spatial step length becomes longer. The stability condition of the proposed scheme approaches the CFL condition for low spatial resolution. The performance of the proposed and implicit schemes, which have wider stability conditions than that of the CFL, is not trivial because they are not truly fully-implicit as the continuity equation is discretized explicitly. We consider that the performance of the implicit and proposed schemes with the numerical stability even wider than the CFL condition is due to the non fully-explicit treatment of the momentum equation. Theoretical stability analysis of this phenomenon would be mathematically complicated and is beyond of the scope of this paper. In summary, the stability of the proposed scheme lies between those of the implicit and semi-implicit schemes and is wider than the CFL condition.

Accuracy test
The accuracy of the three schemes is compared with a benchmark test of constant wave propagation (Bates et al., 2010). The computational domain is the 1-D horizontal channel (x = 0-5000 m). The analytical solution was derived in Bates et al. (2010) and is used as a reference solution. Figure 2 shows the simulated results for each scheme, setting spatial and temporal resolution as Δx = 10 m and Δt = 1 s, respectively. All the schemes have the largest magnitude of errors at the wave front as in Bates et al. (2010) and de Almeida and , while its magnitude is far smaller than the scale of the waves. The time series of the l 2 error (the square root of the sum of square error at each point) is plotted in Figure 3. The error of the fully-implicit scheme is the smallest until 6,000 s, but the largest afterward. The proposed scheme shows comparable accuracy to the semi-implicit scheme. Robustness of the schemes is also examined for different resolution (Δx, Δt): (10 m, 1 s), (25 m, 5 s), and (100 m, 20 s). The l 2 error of each scheme at time 6,000 s is summarized in Table SI. The error magnitude increases with the increase in space and time steps for all the schemes.

Application to Tonle Sap Lake
The three schemes are finally applied to simulating shallow surface water dynamics of the Tonle Sap Lake and its floodplains, Cambodia (Tanaka et al., 2018a). The simulation domain is shown in Figure 4. Tonle Sap Lake is located in the southern part of Cambodia and flows into the Mekong River through the Tonle Sap River. Surface water dynamics in and around Tonle Sap Lake flow are driven by upstream river flow and backwater flow from the Mekong River in the rainy season. The 2-D LIM with the semiimplicit, fully-implicit, and the proposed schemes is applied to lake flow simulation from 1998 to 2003. The upstream and downstream boundary conditions are given by tributary river discharges and the water stage at the Prek Kdam station on the Tonle Sap River, respectively, following Tanaka et al. (2018a). The modelling and meteorological data, i.e. rainfall and evapotranspiration, during the tar-get period was also presented in Tanaka et al. (2018a). The spatial resolution is set to 500 (m), which has been found to be sufficiently fine (Tanaka et al., 2018a). The time step is adaptively chosen by multiplying 0.7 by the MATS of the CFL condition so that the semi-implicit scheme can be applied to the simulation. The simulated water stage at the Kampong Luong station (see Figure 4) is shown in Figure  5. All the schemes reasonably reproduced the observed data over the period and demonstrated small differences with each other. The time series of the difference of the simulated level from the observed one for each scheme is shown in Figure 6. The semi-implicit scheme had the lowest accuracy when the water level was low. On the other hand, the proposed scheme showed comparable accuracy with the fully-implicit scheme. Although the simulated water level with the fully-implicit scheme is slightly closer than the proposed scheme, it should be noted that the computation , and 6,000 s (green). Spatial and temporal resolutions were set as Δx = 10 m and Δt = 1 s, respectively. The cross, dots, plus plots show the semiimplicit, fully-implicit, and the proposed schemes, respectively. The analytical solution is shown as the solid line here is based on a set of input observation data, which essentially has observation errors. Another reason may be using a spatially-constant roughness coefficient and limited rainfall and evapotranspiration (Tanaka et al., 2018a).
For practical applications, computational efficiency is as important as its accuracy. The computational time of the semi-implicit, fully-implicit and proposed schemes were 4.017, 4.0667, and 4.6333 hours on a Linux machine with 34 cores (68 threads with Hyper-Threading) of Intel (R) Xeon (R) E5-2695 v4 (2.10 GHz). We implemented OpenMP parallelization to divide the calculation of the momentum and continuity equations among the 2dimensional grids into the 34 cores. The computational time of the proposed scheme was 15% larger than that of the semi-implicit scheme due to "if" sentences for Equations (13) to (16), but the difference is not significant, and computational efficiency of the three schemes are consid-  ered to be comparable. Overall, the computational results obtained indicate that the proposed scheme can serve as an effective alternative to the semi-implicit scheme.

CONCLUSIONS
This paper proposed a new scheme for the momentum equation of the LIM with the exact solution-based discretization, as an alternative to the conventional semi-implicit scheme. It was found that the new scheme can be implemented without any iterative solvers avoiding high computational costs. In addition, the scheme is physically more consistent and has higher stability than the semi-implicit scheme, demonstrating theoretical advantages over the semi-implicit one. The proposed scheme was also compared with the recently-developed implicit scheme, showing comparable theoretical properties. Numerical computation of the three schemes against a test case revealed that the proposed and implicit schemes have almost the same accuracy that was higher accuracy than that of the semiimplicit scheme. Therefore, this paper clearly demonstrated that the semi-implicit scheme was less stable and less accurate than the remaining two schemes. This is a remarkable Figure 5. Water level at the Kampong Luong station on the Tonle Sap Lake, Cambodia. The dots show the observation; the red, blue, and green lines show the semi-implicit, fullyimplicit, and proposed schemes, respectively Figure 6. Differences in the simulated water level at the Kampong Luong station from the observed one for the semi-implicit (red), fully-implicit (blue), and proposed (green) schemes finding since the semi-implicit scheme is currently used in various hydrological simulators. Application of the three schemes to surface water dynamics of Tonle Sap Lake demonstrated that they are applicable to the flow with backwaters and wetting and drying interfaces. It was found that their computational results are not critically different, but the semi-implicit scheme had the lowest accuracy when the water level was low.
In summary, this paper demonstrated that the proposed scheme is theoretically and computationally advantageous to the semi-implicit scheme and has almost the same performance as the implicit scheme. What is important in this paper is that the proposed scheme can serve as an accurate, stable, and consistent alternative to the conventional semiimplicit scheme. A future research topic is comparison of the three schemes against the shallow water equations. It will be interesting to see whether the present approach can be applied to the discretization of the shallow water equations.