A kinetic flux vector splitting scheme for shallow water equations incorporating variable bottom topography and horizontal temperature gradients

This paper is concerned with the derivation of a well-balanced kinetic scheme to approximate a shallow flow model incorporating non-flat bottom topography and horizontal temperature gradients. The considered model equations, also called as Ripa system, are the non-homogeneous shallow water equations considering temperature gradients and non-uniform bottom topography. Due to the presence of temperature gradient terms, the steady state at rest is of primary interest from the physical point of view. However, capturing of this steady state is a challenging task for the applied numerical methods. The proposed well-balanced kinetic flux vector splitting (KFVS) scheme is non-oscillatory and second order accurate. The second order accuracy of the scheme is obtained by considering a MUSCL-type initial reconstruction and Runge-Kutta time stepping method. The scheme is applied to solve the model equations in one and two space dimensions. Several numerical case studies are carried out to validate the proposed numerical algorithm. The numerical results obtained are compared with those of staggered central NT scheme. The results obtained are also in good agreement with the recently published results in the literature, verifying the potential, efficiency, accuracy and robustness of the suggested numerical scheme.


Introduction
In the recent years, the study of hyperbolic systems with source terms has attracted much attention in the field of computational fluid dynamics (CFD) due to their wide rang physical and engineering applications. One example of such systems is the set of Saint-Venant equations which are widely used in the ocean currents and hydraulic engineering to describe the hydraulic jumps and shallow water flows [1][2][3]. In the early decades many numerical methods have been developed to investigate such flows.
The present work deals with the numerical approximation of shallow water equations which have temperature fluctuations of prime interest [4]. This set of shallow water equations, a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 containing horizontal temperature gradient terms, is called as Ripa system and is used to investigate ocean currents [5,6]. The main difference between the well known shallow water model and the Ripa system is the definition of steady states at rest. Indeed, the Ripa system involves the steady states at rest governed by the system of partial differential equations. This system is a non-homogenous hyperbolic system with non-zero source terms which appear after incorporating horizontal temperature gradients in the shallow water equations [7][8][9][10][11][12].
The classical nonlinear shallow water equations were derived in 1871 by de Saint-Venant [13]. Currently these equations are widely used in practice and one can find thousands of publications devoted to the applications, validations and numerical solutions of these equations [14]. In the pioneering work of Savage and Hutter [15] a shallow water type model in onespace dimension has been proposed to study aerial avalanches, followed by an extension of their approach to two-space dimensions [16].
Several well-balanced and higher order accurate computational methods have been designed for solving Saint-Venant equations and Ripa system in one and two space dimensions. Some of the well known finite volume schemes for solving systems of conservation laws include the finite volume weighted essentially non-oscillatory (WENO) scheme, conservation element solution element (CESE) scheme, central upwind (CUP) scheme, central Nessyahu Tadmor (NT) scheme, as well as Bhatnagar Gross Krook (BGK) and kinetic flux vector splitting (KFVS) schemes [17][18][19][20][21][22][23]. In this work, the KFVS scheme is extended to solve the current shallow flow model equations.
In the literature, many positivity preserving and well-balanced schemes have been suggested for solving single and two phase shallow flow models [24][25][26][27][28][29][30][31]. An efficient computational scheme should preserve the positivity of h and θ and should capture the steady states and small perturbations of the solutions. In addition to the above requirements, a good numerical method should avoid developing spurious oscillations of pressure in the neighborhood of temperature jumps. Such type of typical oscillations may appear when a Godunov-type finite volume scheme is applied [32].
Zhou et al. have developed a well-balanced and robust numerical scheme based on the Harten Lax and van Leer (HLL) Riemann solver and surface gradient method (SGM) for interpolating the initial data [33]. The main advantage of their approach is the consideration of h + B for the initial reconstruction instead of h that guarantees the achievement of steady state. In this manuscript, the same h + B is used as a conservative variable. Audusse et al. have developed a first-order accurate KFVS scheme based on the collisionless Boltzmann transport equation for the system of shallow water equations [34]. That scheme is well-balanced based on the appropriate choice of flux splitting but the accuracy of that scheme was not sufficient for some test problems. Furthermore, for the shallow water equations the source terms effects the evolution around the cell interface and has tendency to change the flux terms. Thus, special attention is needed to the discretization of source terms. Our proposed KFVS scheme, which is second order accurate and discretizes the source terms appropriately, has the capability to overcome all such difficulties.
The Godunov upwind finite volume schemes are based on the hyperbolic structure of the underlying conservation laws [4,35]. However, for some systems, it is difficult to find the explicit expressions of eigenvalues. Thus, the development of Riemann solvers are not possible for such systems. In such scenarios, the kinetic theory based numerical schemes can be very interesting if the Maxwellian distributions are available for such systems. One of the fascinating aspects of the kinetic schemes is that when applied to Euler equations of gas dynamics, they preserve the positivity of mass, density and pressure. As a result, the kinetic schemes are unconditionally stable in the L 1 -norm. Moreover, they also possess the entropy property as a consequence of the celebrated Boltzmann H-theorem [35]. Since the form of the numerical dissipation in the FVS scheme is consistent with the Navier-Stokes viscous terms, the robustness of the KFVS schemes can be easily understood, i.e., the absence of numerical shock instability. These advantages of the kinetic schemes are also very attractive when solving the current shallow flow model [31].
There are two types of gas kinetic schemes. One of the well-balanced gas kinetic schemes is KFVS scheme which is based on the collionless Boltzmann equation, while the other scheme is based on the BGK model [36,37]. On combining the dynamical effects of both gas evolution and projection stages, the governing equations are physically same for both KFVS and BGK schemes. The major difference between both the two schemes is that the particle collision time τ in BGK scheme is replaced by the Courant Friedrichs Lewy (CFL) time step Δt of the KFVS scheme [38]. As the structures of both Ripa system and shallow water model are very similar, we can expect the same performance of this well balanced scheme for the Ripa model. Furthermore, as the shallow water equations can be recovered by providing a constant temperature θ, the extended scheme for the Ripa system will also preserve the lake at rest steady states exactly for the constant θ.
The KFVS scheme is based on the Boltzmann equation of kinetic theory of fluids. Due to their intrinsic upwinding, multidimensional and free of Riemann solver properties, gas kinetic schemes have got increasing attention of the researchers in this field [39][40][41][42]. The current KFVS scheme has capability to tackle all difficulties associated with the non-conservative part of the equations. The basic idea behind this scheme is that it uses information of local propagation speeds and estimates solutions in the form of cell averages. Moreover, this scheme has an upwind nature, as it takes care about the direction of propagation of wave by measuring the one-sided local speeds. Apart from its advantages, the KFVS has some limitations as well. For example, the scheme strictly depends on the Boltzmann equation and Maxwellian distribution function. In the case of complications in the Maxwellian distribution function, the KFVS method becomes complicated and requires more computational effort. The application of KFVS is almost not possible in the absence of Maxwellian distribution function. This scheme has been successfully applied to approximate several models in gas dynamics. Mandal and Deshpande [43] have used this scheme to investigate bump in a channel problem and it was found that the explicit flux function of KFVS scheme was similar to the flux function of Van Leer [44]. For the simulation of two dimensional problems for structured and unstructured meshes, the high order KFVS schemes have been used [22,45]. In the current KFVS scheme, we start with the cell average initial data of conserved variables and calculate the cell averaged values in the same cells at updated time step. In order to achieve second order accuracy of the scheme, MUSCL-type initial reconstruction procedure and Runge-Kutta time stepping method are used.
As particles have the ability to move in forward and backward directions with the fluid motion, the fluxes of conservative variables can be divided into forward and backward fluxes at the cell interface, i.e., where, W i represents the vector of conserved variables within the cell [x i−1/2 , x i+1/2 ] and F(W i ) denotes the vector of flux functions. For validation, the results of our proposed scheme are compared with those obtained from the staggered central NT scheme [10,12]. The main reason for choosing this scheme is its simplicity, efficiency, robustness, Riemann-solver free nature and second-order accuracy. The current KFVS scheme is developed in such a way that it preserves the positivity of physical variables like height, temperature and pressure. Thus, the scheme gives stable results.
The present article is organized as follows. Sections 2 and 3 are devoted to the introduction of one-dimensional Ripa system and to the derivation of KFVS scheme, respectively. The twodimensional (2D) Ripa system and the corresponding 2D KFVS scheme are presented in the Sections 4 and 5, respectively. Numerous numerical case studies are carried out in Section 6. Finally, concluding remarks are given in Section 7.

Single-phase one-dimensional Ripa model
Firstly, we consider the one-dimensional single-phase shallow flow model with variable bottom topography [5][6][7] @ t h þ @ x ðhuÞ ¼ 0; ð2Þ where h is the flow height, g is the gravitational acceleration constant, u is the flow velocity along the x direction, B ¼ BðxÞ; x 2 R denotes the bottom topography from a given level and θ represents the temperature fluctuation function. The stationary steady states for this model are expressed as Eqs (2)-(4) can be written in the compact form as where W is the vector of conserved variables, F represents the vector of fluxes and τ denotes the source term. The conserved variables are give as and the conservative fluxes are given as Moreover, the non-conservative source terms are expressed as Above equations can be expressed in quasi-linear form as where W = (w 1 , w 2 , w 3 ) T , R(W) = [0, −ghθ@ x B, 0] T and A(W) is the Jacobian matrix whose element at the k-th row and l-th column is @f k @w l for k, l = 1, 2, 3 and the functions f k (w 1 , w 2 ) are defined in Eq (8). The matrix A(W) is given as The eigenvalues of the Jacobian matrix A(W) are λ 1 = u, l 2 ¼ u þ ffiffiffiffiffiffiffi ghy p and l 3 ¼ u À ffiffiffiffiffiffiffi ghy p . Since these eigenvalues are all real, so the one-dimensional Ripa system under consideration is hyperbolic.

KFVS scheme for one-dimensional single-phase shallow flow model
The flux of conserved variables is related to the particle motion across cell interfaces. For the one-dimensional flow, the particle motion in this direction determines the flux function F(W).
In statistical mechanics, the distribution of moving particles in the x-direction can be considered by local Maxwellian distribution function. The Maxwellian distribution function in normal direction n 2 {x} is given as [22,46], where u n is the average fluid velocity in the n-direction, v n is the individual particle velocity in the same direction, and λ is the normalization factor of the distribution of random velocity. The transport of any flow quantity is due to the movement of particles. Considering the distribution function f M in Eq (12), particles can be split into two groups. One group move to the right side with positive velocity (u n > 0) and the other group move to the left side with negative velocity (u n < 0). Before splitting the fluxes, let us define The above two moments are sufficient to split all the fluxes. In order to simplify the notation, one can define [22] hv 0 i þn ¼ and hv 1 i þn ¼ Here, the complementary error function is defined as In order to apply KFVS scheme, first we divide the domain of interest into N sub-domains. Let 1 2 represents the uniform width of cell, the points x i = iΔx denote to the cells center and the points x iAE 1 2 ¼ x i AE Dx=2 represent the cells faces. We start the process with a cell averaged initial data W n i at time step t n and compute the cell average updated solution W nþ1 i over the same cell at the next time step t n+1 .
Let us consider the one dimensional single-phase shallow flow model. Using the above relations, the flux functions in Eq (8) can be splitted as Following [22], the flux vectors at the left and right interfaces of the cell C i are given as where, F AE i ≔ F ± (W i ). Thus we have, Integration of Eq (6) over the cell gives the following semi-discrete kinetic upwind scheme where, W i denotes the cell averaged values of the conserved variables. Moreover, following [32] a modified Godunov type approach is taken into account for the discretization of the non-conservative term. In this approach, non-conservative terms appear on the right hand side are approximated as cell averaged values and the differential terms are approximated by using flux splitting technique, Here, the cell averaged values W i are defined as and F iþ 1 2 and F iÀ 1 2 are given by (21). The above scheme is only first order accurate in space. To achieve high order accuracy, the initial reconstruction strategy must be applied for the interpolation of cell averaged variables W i . Here, a second order accurate MUSCL-type initial reconstruction procedure is employed. Starting with a piecewise constant solution W i , one can reconstructs a piecewise linear (MUSCL-type) approximation by selecting slope vector (differences) W x . Then, boundary extrapolated values are given as A possible computation of these slopes is given by the family of discrete derivatives parameterized with 1 ϑ 2, for example where, Δ denotes central differencing and 1 ϑ 2 is a parameter. For ϑ = 1 the limiter is most dissipative and is least dissipative for ϑ = 2. Moreover, and MM denotes the min-mod nonlinear limiter MMfx 1 ; x 2 ; :::g ¼ On the basis of above reconstruction, a semi discrete high resolution kinetic solver is formulated as In the above equation, To obtain second order accuracy in time, a second order total variation diminishing (TVD) Runge-Kutta (R-K) scheme is used to solve Eq (30). Let us denote the right-hand side of Eq (30) as L(W). A second order TVD R-K scheme updates W through the following two stages [22] W ð1Þ ¼ W n þ DtLðW n Þ; ð32Þ where, Δt = t n+1 − t n represents the time step. In all our 1D computations, we have chosen the dynamic time steps according to the following CFL condition: where σ i is the maximum eigenvalue of the Jacobin matrix A(W) given by Eq (11).

Single-phase two-dimensional Ripa model
Now, we consider the two-dimensional Ripa system in which the water temperature fluctuations are taken into account. The two dimensional Ripa model has the following form @ t ðhyÞ þ @ x ðuhyÞ þ @ y ðvhyÞ ¼ 0: Here, h represents the water height, u and v are respectively the fluid velocities in the x and y directions, B = B(x, y) denotes the bottom topography and the gravitational constant is represented by g. The variable θ is referred to as the potential temperature field and hu and hv are the momentums along x and y directions, respectively. The stationary steady states for this system are In compact form, the above system can be rewritten as where W is the vector of conserved variable, F(W) and G(W) represent fluxes along x-axis and y-axis respectively, and τ denotes the vector of source terms. The conserved variables are given as and the conservative fluxes are expressed as Moreover, the non-conservative source terms are defined as The quasi-linear form of the above system is as follows Eigenvalues of the Jacobian matrix H(W) are α 1 = α 2 = u, a 3 ¼ u þ ffiffiffiffiffiffiffi ghy p and Similarly, the eigenvalues of the Jacobian matrix J(W) are . As the eigenvalues of two-dimensional Ripa system are all real, so it is a hyperbolic system having independent eigenvectors. To calculate the time step of the numerical scheme, eigenvalues are needed in the time step relation for ensuring the stability of the scheme. If we take u ¼ v ¼ AE ffiffiffiffiffiffiffi ghy p , the Jacobian of the Ripa system will not provide a complete set of eigenvalues, i.e. the Ripa system produces a resonance phenomenon. For that reason, the solution of the Riemann problems becomes very difficult in such situations. The current KFVS for the Ripa system avoids such complexities because it calculates the cell interface fluxes from the Maxwellian distribution function and it does not require a Riemann solver.

KFVS scheme for two-dimensional single-phase shallow flow model
Now consider the two-dimensional single-phase shallow flow model in Eqs (35)- (38). In this case, the flux is associated with the particles motion across the boundaries of cell interfaces. In the two-dimensional case, the flow of particles along each characteristic direction x and y is considered using the macroscopic flux functions F(W) and G(W) through boundaries of the mesh cells. For example, for flow in the x-direction, the flux function is determined by the particle motion in that direction. The remaining quantities, e.g. velocity in y-direction and water height may be treated as passive scalars which are moved with the x-direction particle's velocity. The same Maxwellian distribution f M in Eq (12) is considered with n 2 {x, y} to split the flux functions in the x and y-directions. In the x-direction, we have where, Similarly, in the y-direction where, Let N x and N y be large integers in the x and y-directions, respectively. We assume a Cartesian grid with a rectangular domain [x 0 , x max ] × [y 0 , y max ] which is covered by cells C ij ≔ [x i−1/2 , x i+1/2 ] × [y j−1/2 , y j+1/2 ] for 1 i N x and 1 j N y . Here, (x N x +1/2 , y N y +1/2 ) = (x max , y max ) The representative coordinates in cell C ij are denoted by (x i , y j ). Therefore, and Thus, fluxes at the cell interfaces are described as Integration of Eq (17) over C i,j gives the following semi-discrete kinetic upwind scheme where, the cell averaged values W i,j are defined as The cell interface fluxes F iþ 1 2 ;j and G i;jþ 1 2 are defined by Eq (54). Moreover, T iþ 1 2 ;j − T iÀ 1 2 ;j and T i;jþ 1 2 − T i;jÀ 1 2 are expressed as The same reconstruction procedure can be applied in direction by direction manner which was discussed in Eqs (26)-(29) of the previous section. Therefore, the final form of the 2D-KFVS scheme is given as The second order TVD RK scheme (Eqs (32) and (33)) is used to solve the system of ordinary differential equations in Eq (59). In all our 2D computations, we have chosen the dynamic time steps according to the following CFL condition: where s x i;j and s y i;j are the maximum eigenvalues of the Jacobian matrices H(W) and J(W) given by Eq (46). The proposed KFVS schemes for 1D and 2D flows can be analogously extended to the 3D flows. In that case, the same Maxwellian distribution can be used to split the flux functions along the characteristic directions x, y and z. For example, for flow in the x-direction, the flux function is determined by the particles motion in that direction. The remaining quantities, such as velocities in the y and z directions and water height may be treated as passive scalars moving with x-direction velocity of particles. Moreover, the current KFVS scheme can be easily extended to solve multidimensional incompressible and compressible multiphase flow models, see for example [39,46].
There are so many techniques and procedures which can be used to further enhance the accuracy and reduce the computational cost of the above KFVS schemes. For example one can increase the order of accuracy by using WENO-limiters and can apply adaptive meshing techniques to reduce the computational cost [47]. One can also adopt Singh's numerical scheme for further improvements [48]. Moreover, one can also try to establish a variational formulation for the discussed systems [49].

Numerical case studies
In this section, several case studies are carried out in one and two-space dimensions. The initial data of these Riemann problems contain a single discontinuity. There are conditions, called as "shock conditions", which decide the occurrence of shocks. These conditions can be used to find the state on the other side of the shock if a state on the one side of the shock is already known [4,50]. The results obtained from KFVS scheme are compared with the results of central NT scheme. The main reason for choosing NT scheme is its simplicity, robustness and efficiency. Moreover, the scheme does not requires any exact or approximate Riemann solver. If we want to apply a Godunov upwind scheme then we have to derive the exact Riemann solver and a complete set of eigenvalues which are very difficult or almost impossible. The KFVS scheme is applied for the first time to the Ripa system incorporating variable bottom topography and horizontal temperature gradients. In all the test problems outflow boundary conditions are used.

One-dimensional case studies
Here, we present the one-dimensional test problems. In the first problem, a variable bottom topography and smooth initial data are considered to calculate L 1 -errors and     In the literature, the accuracy and efficiency of the KFVS scheme has been analyzed and verified for several nonlinear hyperbolic systems, see for example [21-23, 31, 38, 39, 41, 46]. Problem 2: Dam break over the flat bottom. This dam break problem over the flat bottom was considered in [7]. Here, we have solved it for the Ripa The numerical solution is calculated at t = 0.2 in the computational domain [−1, 1] using the KFVS scheme and its results are compared with those obtained from the central NT scheme. The computational domain is divided into 200 grid points. We noticed that when waterbed is flat, the source term appearing in the Ripa system vanishes and the system can be easily solved by using any finite volume scheme. Fig 5 shows the plots of h, u, θ, p, hu and hθ. A good agreement can be observed between the results of KFVS and central NT schemes. The variable bottom topography function is defined as Both the schemes are found well-balanced and preserve the positivity of w, h and θ. Problem 5: Small perturbation of a steady state solution. This one-dimensional experiment describes a small perturbation of a steady-state solution [7]. The initial conditions are given as ðw; u; yÞ ¼ Our last one-dimensional problem of variable bottom topography is a discontinuous waterbed problem which is taken from [6] for the case of Ripa system. The waterbed for the rectangular bump is defined as In this experiment, we compare the results of KFVS and central schemes for flat bottom topography. The 3D plots of KFVS method and the 1D plots for comparison of the two schemes are shown in Fig 10 at t = 0.15 using the uniform mesh Δx = Δy = 0.02. When the dam breaks, a shock wave is generated which moves radially outwards, a rarefaction wave travels radially inward and a contact wave moves between the shock and rarefaction waves. One can easily see that a good agreement is found between both schemes. It is observed that some small oscillations are generated in the contact area, which can be seen in the contour plots of KFVS scheme presented in Fig 11. Problem 9: Perturbation of steady states on an irregular waterbed. This 2D steady state problem was considered for Ripa system in [6,7]. It contains a bottom topography function having two Gaussian shaped humps. The initial conditions, featuring the steady states of a lake at rest are given below ðw; u; v; yÞ ¼  The numerical solutions are obtained at t = 0.15 and perturbation propagation can be seen in the plots given in Figs 14 and 15. The 3D plots of KFVS scheme and the 1D plots for comparisons of both schemes are given in Fig 14. The solution obtained from central scheme develops radial shaped small oscillations in the pressure solution, which interact with the perturbation. Due to the interaction of circular shape pressure oscillations and perturbation, parasitic waves appear in the contour of pressure that can be seen easily in Fig 15. Moreover, KFVS scheme shows its capability to capture perturbation more accurately.

Conclusion
In this paper, we have developed a high resolution KFVS scheme to approximate the one and two-dimensional Ripa systems. The KFVS method was applied for the first time to solve the Ripa systems for both flat and non-flat bottom topographies. The proposed scheme has the ability to capture sharp discontinuities, efficiently handles dry bed regions, and is extendable for multi-phase flow problems. The scheme preserves the non-negativity of the flow height and temperature. To satisfy the well-balanced nature of the scheme for the Ripa system, we have discretized the source term analogously to the discretization of flux divergence terms. The scheme was validated by taking into account some classical test problems already available in the recent literature. The results of KFVS scheme were compared with those obtained from the staggered central NT scheme and the reference solutions. The numerical results of both schemes were found in good agreement with each other and those available in the literature. However, it was observed that KFVS gives better resolution of sharp discontinuities as compared to the central scheme. The scheme has preserved the well-balanced property and nonnegativity in the problems of steady sate perturbed variable waterbed. This confirms the efficiency of the KFVS method to accurately solve the one and two dimensional Ripa systems. The considered numerical simulations include the dam-break problems, small perturbation problems of variable bottom topography, and varying temperature gradients.