On the Efﬁciency of Staggered C-Grid Discretization for the Inviscid Shallow Water Equations from the Perspective of Nonstandard Calculus

: This paper provides a rationale for the commonly observed numerical efﬁciency of staggered C-grid discretizations for solving the inviscid shallow water equations. In particular, using the key concepts of nonstandard calculus, we aim to show that the grid staggering of the primitive variables (surface elevation and normal velocity components) is capable of dealing with ﬂow discontinuities. After a brief introduction of hyperreals through the notion of inﬁnitesimal increments, a nonstandard rendition of the governing equations is derived that essentially turns into a ﬁnite procedure and permits a convenient way of modeling the hydraulic jumps in open channel ﬂow. A central result of this paper is that the discrete formulations thus obtained are distinguished by the topological structures of the solution ﬁelds and subsequently provide a natural framework for the staggered discretization of the governing equations. Another key of the present study is to demon-strate that the discretization naturally regularizes the solution of the inviscid ﬂow passing through the hydraulic jump without the need of non-physical dissipation. The underlying justiﬁcation is provided by analytically studying the distributions of the ﬂow variables across an inﬁnitesimal thin hydraulic jump along with the use of hyperreal Heaviside step functions. This main ﬁnding is shown to be useful to comprehend the importance of the application of staggered ﬁnite difference schemes to accurately predict rapidly varying free-surface ﬂows. A numerical experiment is provided to conﬁrm this result.


Introduction
The inviscid shallow water equations describe the behavior of a shallow incompressible and inviscid fluid layer and are suitable to model hydrodynamics in coastal seas, estuaries, lakes and rivers. They are derived from the depth-integrated Euler or Navier-Stokes equations under the hydrostatic pressure assumption. Since both these equations and the Euler equations of gas dynamics have a similar mathematical structure, reflecting the resemblance between the hydraulic jumps and shock waves in compressible flow, conventional high resolution schemes for treating shocks are likewise a commonplace tool in the numerical modeling of rapidly varying shallow water flows.
The theory of weak solutions is the traditional approach in the framework of numerical analysis of nonlinear hyperbolic partial differential equations [1,2]. This theory necessitates these equations to be expressed in conservation form, and has thus created a natural ground for the development of finite volume methods [3,4]. However, such methods often cannot prevent the obtained solutions from being unphysical when subjected to discontinuities (entropy-violating solutions), and heuristic attempts are then made to adequately recover the discontinuous behavior. Shock-capturing techniques, usually with an entropy fix [5,6], and well-balanced methods are particularly appealing in addressing such needs. They routinely regularize the governing equations by some degree of non-physical dissipation.
For example, applications of approximate Riemann solvers to shallow water equations are found in, e.g., [7][8][9]. They aim to exploit the hyperbolicity property of the conservation laws with source terms, which is, however, substantially problematic to master the development of numerical codes. The lack of well-balancing in the resulting numerical schemes is perhaps the most prominent one [10]. This imbalance typically leads to artificial flow transitions over discontinuous bed topography in water initially at rest. Although the design of well-balanced Riemann solvers has been frequently reported in the literature with varying degrees of success, one may still encounter the risk of approximating a non-physical solution, even though by refining the grid.
Alternative methods have been proposed that avoid the necessity to solve a Riemann problem, such as the central-upwind methods [11,12], the artificial viscosity approach [13][14][15][16], and the energy-stable schemes using numerical diffusion operators based on entropy variables [17], just to mention a few. Though this abundance of methods may have proven to be a reasonable means of predicting the evolution of shock waves, their mathematical theory does in general not constitute a conclusive approach, as it provides no clear-cut measures to regularize the inviscid flow equations in a physically consistent manner. It remains to be seen whether any of such discretizations will produce the physically realistic solution, especially when it is highly nonlinear and rapidly varying.
An example involves the introduction of artificial viscosity into a space-centered discretization. Although artificial viscosity has a direct analogy to physical dissipation, it poses a major problem for overall physical accuracy, as its generic forms have an explicit choice regarding shock layer thickness, while one has control over the amount of dissipation through the associated tuning parameters. It is not a trivial task to find suitable values for these parameters because they are typically problem and mesh dependent and, as such, their underlying methods are likely to fail to converge to physically based solutions [15,16].
Over the years, staggered C-grid methods have been developed for the modeling of rapidly varied flows, and have shown remarkable potential in predicting the formation of hydraulic jumps, bores and breaking waves; see, e.g., the papers of Stelling and Duinmeijer [18] and others [19][20][21][22]. The staggered spatial discretization schemes of the type proposed by Arakawa [23] are primarily designed for accurate and efficient computation of incompressible (often gradually varying) flows. The computational efficiency along with the physical accuracy is also the central theme in the above mentioned papers. Furthermore, staggered C-grid schemes exhibit distinguished properties such as the absence of spurious pressure modes and the local mass conservation. These attractive numerical properties are due to the staggering of the primary unknowns in a grid, with pressure at the cell center and the normal velocity components at the center of the cell faces. Yet, to the best of author's knowledge, such methods are the only class of discretization methods that are capable of accurately resolving the hydrodynamic shocks under various flow regimes without recourse to any explicit solution regularization. See [22] for details.
The main contribution of this paper is to provide a rationale for this finding by using the language of nonstandard calculus. Here, we first give a brief background to this branch of modern mathematics and in particular its natural discrete representation of physical phenomena, which is then followed by an outline of the current work.
Historically, the Newtonian mechanics was closely connected with the early development of calculus that proceeded for centuries on heuristic arguments of infinitesimals [24,25]. The paradoxical difficulties raised by infinitesimals were thoroughly bypassed by employing the ε − δ definitions of limit, as introduced in the nineteenth century. In the 1960s, a rigorous framework for dealing with infinitely small and infinitely large quantities as numbers has entered the area of pure mathematics. This new application of infinitesimals and infinite numbers (not to be confused with infinity ∞), termed nonstandard analysis, was invented as an alternative to the ε − δ limit to put the means of calculus on a sound footing. Even more so, the use of infinitesimals is, to a greater extent, effective and insight-ful than the traditional analysis in terms of ε and δ. Today, nonstandard analysis helps to reconcile the way infinitesimals are used in, e.g., applied mathematics, physics and engineering with some formalism.
Nonstandard analysis basically extends the real numbers to the hyperreal numbers that involve infinite and infinitesimal numbers [26,27]. Because of the existence of infinitesimals in the hyperreal number system, the process of limit is no longer an issue. The implication is that the interpretation of the derivative is viewed as a finite process rather than an infinite limiting procedure [27]. Specifically, either a derivative is utilized for finding the limit of a difference quotient of vanishing quantities, that is, measuring the instantaneous rate of change, or a derivative is considered as the ultimate limit of the ratio of small but finite quantities. While the former statement is related to the ε − δ concept of calculus, the last assertion has a physical meaning, in the sense that no physical quantity can be measured to infinite precision. Yet, once the limit is performed, a continuum model is obtained that does no longer resolve features on the various length scales; they actually contracted to (zero-dimensional) points. Moreover, all physically relevant quantities must then be infinitely differentiable.
The extension to the hyperreals comes equipped with a hierarchical lattice-like structure on the Euclidean space. Such a structure shares many of the properties of topological spaces and is essentially simpler than the continuum. In particular, the presence of a non-vanishing length scale is distinctly defined. It thus provides a means for developing a discrete model to describe a physical system with very large degrees of freedom. As it will be demonstrated in the present study, this becomes particularly relevant for problems with discontinuities. In fact, as we will see, this approach inherently deals with the regularization of the discontinuity.
Viewed in this way, the discrete model is more rich in physics than the continuum model. This may seem counterintuitive, as it is deeply ingrained in our conceptual understanding of the differential equation that it provides a truer description of a physical problem than any discrete formulation that approximates the original equation in a way that it produces solution sufficiently close to the "exact" solution. Note that we have put the word exact in quotes to emphasize the difference with the actual, that is, physical-based solution to the problem in question.
That said, the crux of nonstandard calculus lies in the development of mathematical models with a degree of computability, which also have a meaningful physical sense. Even though the significance of this view has long been recognized in the field of nonstandard analysis, see, e.g., [28,29], it is still relatively unknown to applied mathematicians, physicists and engineers, given the fact that the large majority of mathematical models in physics and engineering are based on continuous partial differential equations.
The organization of the paper is as follows. In Section 2, first a basic understanding of the concepts of nonstandard analysis is presented. We address the natural extension of differentiation to the hyperreal numbers by means of infinitesimal increments. An important feature of this extension is that the derivative turns into a discrete process while handling the discontinuities adequately. Inspired by this key observation, a nonstandard rendition of the inviscid shallow water equations is formulated that expresses exactly the discrete conservation of mass and momentum across a hydraulic jump in an open channel with uniform bed. Such discrete formulations additionally advocate the use of the staggered C-grid arrangement and are amenable to be numerically solved. We then highlight that the continuous form of the traditional shallow water equations is the standard part of the obtained finite difference equations, while they have lost properties typically associated with hydrodynamic shocks. To gain some insight on the solution regularization, Section 3 analyses the distribution of the flow field across the hydraulic jump, where the nonstandard Heaviside functions are employed to model smooth jump functions for the flow variables such as water depth, velocity and mechanical energy. A key result of this nonstandard analysis is that, unlike the primitive variables, the energy flux displays an undershoot in the hydraulic jump, suggesting that the viscosity does not play a crucial role, as similarly observed by Salas and Iollo [30], who demonstrated a local maximum of the entropy profile inside an inviscid gaseous shock. Numerical computations to study the formation of a hydraulic jump in a rectangular channel are performed as detailed in Section 4. Finally, conclusions are stated in Section 5.

Nonstandard Calculus
A complete introduction to nonstandard analysis can be found in the textbooks of Goldblatt [26] and Keisler [27]. In the following, the building blocks of nonstandard calculus, in particular, those used in this work, are briefly reviewed.
Basically, nonstandard analysis deals with an extension of the real number line R including nonzero infinitesimals and infinitely large numbers, along with the real numbers. This extended set is denoted by * R, and its elements are called hyperreals (or nonstandard reals). An infinitesimal is, in absolute value, a number that is smaller than every positive real number. Zero is the only real number that is an infinitesimal. The reciprocal of a nonzero infinitesimal is a hyperlarge number (or infinitely large or simply infinite), which is larger than any real number, though strictly smaller than infinity. If a number is not hyperlarge, it is called finite (or limited). Accordingly, every real number is finite. Finally, a number which is limited and not infinitesimal is called appreciable.
Neither nonzero infinitesimals nor infinitely large numbers exist on the real line. Yet, they can be treated in much the same way as standard numbers. This property is known as the transfer principle. As an example, the reals are an ordered field, and so are the hyperreals; infinitesimals differ in magnitude from other infinitesimals, and hyperlarge numbers from other hyperlarge numbers. For details on the construction of the nonstandard reals, we refer to Goldblatt [26].
Let ε be an infinitesimal and x and y be finite hyperreals. We define the relation x y as x − y is infinitesimal or zero. It immediately follows that ε 0. This corollary essentially captures the idea of "arbitrary close" that lies at the heart of the well-known ε − δ limit in calculus.
Furthermore, let a be a real number. The standard part of x, denoted by st(x), is a unique real a infinitely close to x, thus a = st(x). Consequently, st(ε) = 0. A hyperlarge number cannot have a standard part. The standard part function provides an essential link between the finite hyperreals and the real numbers.
Next to the hyperreal numbers, one also deals with hyperfunctions. Given a function f from R to R, there is a unique function * f : * R → * R, called the * -transform of f , such that * f (x) = f (x) if and only if x ∈ R. Thus, * f is a nonstandard extension of f .
Let f : R → R, then the derivative of f at point x, denoted f (x), is defined by for every nonzero infinitesimal δ, provided that f (x) exists independent of δ. Function f is said to be differentiable at point x. Notice that the derivative of f is not defined by a vanishing limit δ → 0, but as the standard part of a difference quotient; f (x) is infinitesimally close to that quotient. As a basic example, we consider f (x) = x 2 . Its derivative is calculated as follows: This example clearly demonstrates the effectiveness of using nonzero infinitesimals to determine the derivative without the need of an endless ε − δ limiting process and, in turn, provides a rigorous basis for the calculus. For further illustrative examples, we refer to Keisler [27]. The transfer principle guarantees that all familiar rules involving derivatives, e.g., the product, quotient, and chain rules, also hold for hyperreal-valued functions.
Next, suppose that f (x) exists at point x and let ∆x be a nonzero infinitesimal. Then, the increment of f , denoted ∆ f , due to change ∆x is given by Thus, the standard part of the quotient ∆ f /∆x equals the derivative f (x). This result is known as the increment theorem. By virtue of the transfer principle, this property can be generalized to functions of multiple variables. Note that the infinitesimal ε can be viewed as an error term, which depends on ∆x.
The use of infinitesimals in nonstandard analysis enables us to study topological spaces of reals induced by hyperfinite lattices. (Hyperfinite sets may be treated as finite sets, though they are generally infinite.) A halo (or a monad) of a point x on the hyperreal line is the interval of hyperreals that are infinitesimally close to x. Let µ(x) denote the halo of x, then µ(x) = {y ∈ * R | y x}. While every point on the real line has zero dimensions, the size of each halo on the hyperreal line is infinitesimal. A halo is said to possess a microstructure. There exist infinitely large halos containing the point x. Furthermore, if x is finite, then st(x) ∈ µ(x), while not including other real numbers. Additionally, the halos of different real numbers do not overlap.

Discrete Formulations
A hydraulic jump forms at the transition from the upstream supercritical flow to the downstream subcritical flow. A sudden rise subsequently occurs in the free surface while the flow is rapidly varied. Let us consider the hydraulic jump that moves at a constant speed c in an open channel with a uniform bed. The rate at which the flow changes between the upstream state and the downstream state is determined by the flow continuity and the momentum balance. Following Stoker [31], they can be written, respectively, in the following form and Here, the variables h and u have their conventional meanings of water depth and flow velocity, whereas the subscripts u and d refer to the upstream and downstream states with respect to the jump; see Figure 1. Furthermore, g is the acceleration of gravity. Note that the vertical pressure distribution away from the jump is hydrostatic.
Microscopically, due to viscous forces acting in the large flow gradients within the relatively thin jump layer, the flow is irreversible. Some of the kinetic energy is converted into an increase in potential energy and the remainder is lost through turbulence to heat. Consequently, the downstream state must have a higher entropy than the upstream state.
The rate of energy loss in flow expansion can be computed as follows (see [31]): where e k = g q(h k + (u k − c) 2 /2g) is the mechanical energy (flux) per unit mass (per unit width), with k ∈ {u, d} labeling the position with respect to the jump, and is the mass flux through the jump. The physical laws of mass conservation, momentum transport (Newton's law of motion) and energy loss (the second law of thermodynamics) thus provide the macroscopic balance of flow changes across the hydraulic jump. In fact, this applies to any form of discontinuity that is intrinsically linked to an irreversible process. Yet, the resolution of Equations (1) and (2) adequately resolves the large-scale features of the flow field, while fluid motions that occur at length scales smaller than the width of the jump layer are not captured.
For the discussion below and in Section 3, and without loss of generality, we consider the channel with one dimension of space on the x−axis. The flow passes through the hydraulic jump in the positive x−direction. Accordingly, the traveling flow variables along a characteristic are a function of the similarity variable s = x − ct. By virtue of the increment theorem, we have Let the jump be located at some point S on the characteristic. Now, the extension of the real number line to include nonzero infinitesimals naturally permits to deal with the set of all points separated from S by an infinitesimal distance. This halo of S, denoted µ(S), conceptually resolves the scale of the flow variables across the jump. The size of the halo is given by ∆s, which is a nonzero infinitesimal.
Let the functions h(s) and u(s) be the water depth and flow velocity along the characteristic, respectively. These flow variables are constant on both sides of the infinitesimal jump. (Recall that the jump is uniformly moving with speed c.) Furthermore, the downstream state is distanced from the upstream one by ∆s. Then, mass balance (1) with respect to the jump can be rewritten as Using the variable transformation and a time shift of + 1 2 ∆t, this equation is cast into with h(x, t) and u(x, t) being the water depth and the flow velocity in the space-time domain, respectively. Note that point (x, t) ∈ µ(S) is finite. Furthermore, ∆t and ∆x are the nonzero infinitesimals associated with halo ∆s. Finally, ∆x a represents the infinitesimal area (per unit width) of the halo. Physically, Equation (4) describes the volume balance with which the infinitesimal rate of change in the volume h ∆x a of a halo is equated to the net mass flux q = h u (inflow less outflow) through the boundaries of the halo. This naturally leads to the staggered positioning of the variables h and q in space and in time; the water depth is located at points intermediate between the mass fluxes. It should be noted that here the water depth (or the water level) undergoes a jump in the halo ∆x a .
The partial differential equation describing the conservation of volume can be inferred from the nonstandard Equation (4). Thus, noting that ∆x a = ∆x, and owing to the increment theorem, the standard part of this quotient yields ∂h ∂t In the same way, we derive the momentum balance as follows. From Equation (2),

This equation can be reshaped into space-time form as
and its standard part is then given by ∂hu ∂t Equation (6) describes the infinitesimal rate of change of the total momentum h u ∆x l for an inviscid fluid within a halo ∆x l through which the fluid flows. The first term on the left hand side is the local rate of change in momentum. Furthermore, the second and and the third term express the advective transfer of momentum with the mass flux q as the transporting velocity. This advective acceleration causes a nonlinear spatial change in the flow velocity u.
Now, the rate of change in momentum of the fluid is due to the action of hydrostatic pressure on the fluid on both sides of the halo as given by the right hand side of Equation (6). Again, the momentum and the pressure are carried at alternate spatial locations. Moreover, the direction of the flow acceleration is aligned with that of the net streamwise force. Accordingly, ∆x l must be understood as the infinitesimal length of the path within the halo along which the momentum h u (or depth-integrated velocity) is being transported. Thus, this momentum experiences a jump inside halo ∆x l .
A key feature of the discussed approach is the close association between the primary unknowns and the topology of the halo. In principle, the unknowns are the integrals over the points, lines, surfaces, and volumes, and are referred to as discrete forms of degree 0, 1, 2, and 3, respectively. Furthermore, they constitute a discrete representation for the primitive variables on a primal-dual pair of space lattices, which naturally induces staggering of the variables. For a more extensive survey of this topic, see Tonti [32] and Ferretti [29]. In the current depth-integrated framework, the area integral of h over ∆x a (designated as the primal surface lattice of the halo) and the line integral of h u along ∆x l (i.e., the dual of the edge lattice of the halo ∆x a ) are the discrete 2-form and 1-form, respectively.

Discussion
Equations (5) and (7) are the one-dimensional nonlinear shallow water equations with the unknowns h and h u, suitable for cases with horizontal frictionless beds. (A physically consistent extension to the case of flows over varying bottom topography with bed friction is dealt with in, e.g., [18,22]) Such continuous differential equations are routinely viewed as a mathematical model for describing the dynamics of a physical system while its number of degrees of freedom is infinity; their main merit resides precisely in their uniqueness and universality, that is, they uniquely describe physics at all scales.
The mathematical model for the inviscid shallow water equations is typically presented as a hyperbolic system of conservation laws (e.g., [8,11,16]). In one dimension, this system is of the form ∂W ∂t where W is the vector of conserved variables, F is the vector function of inviscid flux, and S is a source vector. Note that such vectors are expressed as field functions at a given point (x, t) in the continuum. The mathematical structure of the hyperbolic system provides a good basis for the theory of weak solutions involving shock waves [3]. Nevertheless, this mathematical framework is general and does not take into account the physical nature of the variation of hydraulic jumps. In fact, the system may also admit non-physical discontinuous solutions. Within the framework of (usually non-staggered) finite volume methods, a stable high-resolution scheme can be adopted with the aim of constructing the numerical flux vector. Then, the residual of the resulting algebraic equations (as computed by substituting the exact solution) acts as the principal source of discretization error. Provided that the solution is continuous and smooth, this error will be infinitesimally small. However, the situation becomes different when the flow variables are discontinuous so that the solution error is no longer infinitesimal. If a hydraulic jump is located in µ(S), so that irreversible processes intervene, the solution of the difference equations is rather different from the solution of the differential equations. This divergence persists in the limit of vanishing size of µ(S). Obviously, the effect of discontinuity leaves no trace on the solution of differential equations. As explained above, this effect is basically discarded by an application of the standard part function, which in turn implies that the spatial features of the halo (here ∆x a and ∆x l ) are lost. This loss should be restored afterward by the discretization process. This generally calls for measures, often proposed heuristically, to ensure the entropy consistency of the discretization. Examples of such measures have been discussed in the introduction section of this paper.
By contrast, the nonstandard Equations (4) and (6) are principally considered to describe mathematically a physical system with a hyperlarge number of degrees of freedom [28]. This sole principle is built on the existence of nonzero infinitesimals. The resulting discrete equations establish the flow distribution within each infinitesimally small halo exhibiting a geometric structure, which may or may not contain a discontinuity. Furthermore, as it will be revealed in the next section, these equations have a built-in regularization effect on the solution.
The enrichment with a nonvanishing halo size is an essential feature because it paves the way towards grid partitioning of a flow domain. This feature provides a paradigm for constructing discretizations of the nonstandard equations that are naturally adapted for any type of grid system (Cartesian, curvilinear, and unstructured grids) and that behave correctly in a physical sense. By making the size of both ∆t and ∆x appreciable, the above system is translated into a finite but underdetermined system of equations. Its closure is commonly obtained by the addition of a number of constitutive (locally dependent) relations between the different solution fields, which involve the various interpolation (central and upwind) techniques. An example of this process of approximation can be found in Appendix A.
The development of mixed finite volume-finite difference methods that use C-grid staggering on two-dimensional structured and unstructured (triangular) meshes for the modeling of rapidly varied flows in coastal waters is documented in, e.g., [33][34][35][36][37].

Introduction
The aim of this section is to analyze one-dimensional jump distributions for an unsteady, incompressible, inviscid flow across an infinitesimal hydraulic jump on a horizontal bed by utilizing the nonstandard Heaviside functions. We adopt the methodology of Baty et al. [38], who advocate a rigorous approach based on nonstandard analysis to study the one-dimensional shock waves in a compressible, inviscid, perfect gas. However, to the best of our knowledge, no study has been reported for an inviscid hydrodynamic shock.
The significance of the present analysis is that it not only provides valuable insight into the jump solutions of water depth, flow velocity and mechanical energy, but also sheds light on the embedded regularization capability of the governing equations of motion. This is evident by the fact that the analysis deals with continuous hyperreal-valued functions mapped onto the set of hyperreal numbers instead of discontinuous real-valued functions that are mapped to real numbers. This is highly convenient, since it provides a rigorous justification for nonlinear problems that cannot be adequately addressed by the theory of distributions [39,40]. In particular, multiplications of differentiable nonstandard functions in the form of the Dirac delta function or the Heaviside step function are thus straightforward [38].

Nonstandard Heaviside Functions
The initial step of the analysis to be presented here is the introduction of nonstandard Heaviside step functions to properly represent the distribution of the flow field across the jump contained in a halo. Such functions belong to the class of infinitely differentiable functions.
A nonstandard Heaviside function is defined as follows: where ε is an infinitely small number, s = x − ct is a characteristic with c a constant jump celerity, and * h(s) ∈ C ∞ (0, ε) is a strictly increasing hyperreal-valued function representing an infinitesimal jump with thickness ε and This function defines the microstructure (or the microscopic profile) of the flow field inside a shock layer [38,40]. The standard part of * H is unique and is given by In contrast to step function H, any function * H will support all of the operations in calculus required to derive solutions of the governing equations, in particular when dealing with nonlinear terms, e.g., the products of Heaviside functions. Indeed, by virtue of the definition of the standard Heaviside function, we have, for instance, H 2 (s) = H(s), but also H 3 (s) = H(s). Taking the derivative in both cases implies 2HH = H and 3H 2 H = 3HH = H with H the Dirac delta function. Combining the last two equations yields 1 2 H = 1 3 H , which is false. By contrast, * H k (s) = * H(s) for all k > 1, thus preventing this inconsistency.
Following Baty et al. [38] and Salas and Iollo [30], each flow variable is assumed to have a different microstructure. Thus, we look for solutions to h and u of the form Recall that on both sides of the jump, the flow variables are constant. Since we are only dealing with nonstandard functions, we have omitted the superscript * for convenience. Furthermore, the Heaviside functions H(s) and U(s) are taken to have their jumps contained on the same interval (0, ε). In addition, they must satisfy the conditions at the endpoints of the halo, that is and H(ε) = U(ε) = 1 .

Jump Profiles for Flow Variables
To ease the process of deriving smooth functions of h and u from Equations (4) and (6), we adopt their standard parts expressed in nonconservative form. Note that the use of this form poses no problem for the present analysis since the Heaviside functions are differentiable [38]. (By contrast, the theories of weak solutions and distributions require the notion of conservation form.) Hence, we continue our analysis with the following nonlinear equations ∂h ∂t and ∂u ∂t It is worth noting that the governing equations are strictly model equations based on two ideal, controlled conditions: the water is inviscid and incompressible. The latter implies no thermodynamic equation of state. Accordingly, the model equations do not contain enough physical information and consequently the microscopic flow profiles inside a shock layer cannot be deduced (see also [41]).
Another key issue is the local invalidity of the hydrostatic pressure assumption; the pressure cannot be hydrostatic inside the jump. This suggests that Equations (9) and (10) cannot be satisfied by the solutions as given in Equation (8). We verify this by direct substituting and subsequently arranging the resulting equations in matrix form. We then find where the prime denotes differentiation with respect to s. A non-trivial solution [H , U ] exists if the matrix is singular. This is the case if u − c = g h which corresponds to the critical flow condition that must be valid on the whole interval (0, ε). This would imply no jump. Thus, it is plausible to assert that Equation (10) does not describe the complete physics. Along with this finding, we will pursue the analysis by considering the continuity Equation (9) for which we will show below that regularization of the jump discontinuity is accomplished through this particular equation. Since a distinct jump profile for each flow variable cannot be determined, we instead derive a relationship between the jump functions H and U. To this end, we substitute the jump solutions (8) into the continuity equation and subsequently rewrite the resulting equation by means of separation of H and U, as follows: Upon making use of boundary conditions H(0) = U(0) = 0, we obtain the integration constant We observe that the mass flux is constant inside the jump (0, ε) and is equal to q. Hence, no mass can be created nor destroyed within the hydrodynamic shock [42]. This is particularly important because it discloses the incompressibility constraint of the shallow water equations, which can be dealt with effectively by the use of staggered grids [22]. We further note that boundary conditions H(ε) = U(ε) = 1 yield the Rankine-Hugoniot relation expressing the continuity of mass flux across the jump Finally, the sought relation follows from Equation (11) This central result implies that the microstructure of the flow velocity differs from the microstructure of the water depth. Recall that these microstructures remain undetermined.
Note that besides the above jump condition of mass flux, another jump condition related to the conservation of momentum flux can be found. For this purpose, the momentum equation is written in the following form: Recalling that h = h u + [h] H and the constant mass flux q, we obtain which is a nonlinear ordinary differential equation. Integrating both sides of this equation then yields with β the constant of integration. With H = U = 0 for s = 0, the value of this constant is β = 0. However, for s = ε, we have H = U = 1 and obtain the following Rankine-Hugoniot jump condition in the case of uniform bed (cf. Equation (2)) In view of the staggered grid discretization, the non-smooth transition across a hydraulic jump, as manifested by finite changes in the primitive variables u and h on either side of the jump in line with the above jump conditions, is effectively residing in one grid cell. The cellwise continuity of both mass and momentum flux is thus enforced by the staggered spatial placement of the unknowns and remains effective regardless of the grid size. By contrast, central finite volume schemes on non-staggered (or colocated) grids necessarily spread the discontinuity over two grid cells and also allow us to provoke the decoupling between odd and even grid cells. As a result, the constancy of mass and momentum fluxes across an abrupt change cannot be guaranteed. The importance of the Rankine-Hugoniot jump relations in the suppression of odd-even grid oscillations is elaborated in [22].
The remainder of the present section provides an analysis of a direct corollary of the above key result. In particular, we examine the jump distribution for mechanical energy per unit mass, which is given by Like the other flow variables, function e(s) is assumed to have its jump on the interval (0, ε). However, since it is a composite variable, we do not adopt the following solution form: where E(s) is a nonstandard step function. Instead, we will show that function e(s) exhibits a local extremum inside the hydraulic jump, considering that H is strictly increasing on (0, ε). Let the derivative de/ds vanish at point σ p ∈ (0, ε). Then, Recall Equation (12). Differentiation yields Combining the last two equations and then using the fact that H (σ p ) > 0, we obtain which is precisely the critical flow condition, corresponding to a Froude number of unity. Hence, the energy jump e(s) is a non-decreasing function. Since energy is lost when the flow passes through the hydraulic jump, this function must take a minimum value at the critical flow.
As support for this finding, we consider the energy distribution across the jump at the microscopic scale. Within the shock layer, viscosity is then responsible for the redistribution and dissipation of mechanical energy. To show this, we first consider the following momentum equation including the viscous stress term ∂hu ∂t with ν the viscosity of water. Next, we derive the transport equation for the total energy (the sum of the kinetic energy and the potential energy but without bottom topography) by combining the momentum equation and continuity Equation (5). This equation reads ∂η ∂t where η = 1 2 hu 2 + 1 2 g h 2 is the total energy (also termed entropy, see, e.g., [17]) and G = q ( 1 2 u 2 + gh) is the entropy flux. The last term on the right hand side of Equation (14) results in a net loss of energy across the jump, while the first term redistributes the energy inside the hydraulic jump. In particular, the velocity upstream of the jump decreases rapidly, which in turn creates a sharp trough in the energy at the critical flow location, and then recovers towards the downstream side of the jump, which thus causes a rise in the energy level.
A similar phenomenon has been previously reported by Salas and Iollo [30], who have demonstrated the existence of an entropy overshoot inside an inviscid aerodynamic shock using the algebra of nonlinear generalized functions of Colombeau et al. [40]. Their finding has also been justified by the work of Baty and Margolin [41] using nonstandard analysis.
The conclusions that can be drawn from the present analysis are that despite being characterized by the complete absence of viscosity, the nonstandard Equations (4) and (6) do comply with the second law of thermodynamics, and the flow distribution inside the hydraulic jump is effectively regularized by the continuity equation. In particular, incompressibility is an essential component of the solution regularization.
The present analysis also offers support for the ability of the staggered grid approach to reproduce the macroscopic features of hydraulic jumps and bores such as the energy loss and the jump location. This will be further discussed in Section 4.
Yet, the numerical techniques for the solution of shallow water equations by which regularization is based on the addition of artificial dissipation, either explicitly (see, e.g., [14,16]) or implicitly (see, e.g., [8]), do not directly invoke the incompressibility constraint. Indeed, these methods are originally designed for the numerical solution of compressible Euler equations. (An in-depth study of the dissipative regularization of the compressible Euler (and Navier-Stokes) equations based on physical grounds is described in [43].)

An Illustrative Example
For the purpose of exploring the relation between the different microstructures of the flow field across an infinitesimal hydraulic jump, we introduce the nondimensional Heaviside function to model the microstructure for water depth h as defined by where s lies in the halo I = [0, 1] (normalized by the jump thickness ε) and m ∈ * N is a fixed hyperinteger, either appreciable or infinitely large, that determines the actual thickness. It should be noted that the proposed error function is one way to represent a jump. Other examples of Heaviside functions to model the jump can be found in [41]. Also note that function H is monotone increasing on I. The water depth profile is then given by The other variables are determined from this jump function, specifically flow velocity u(s) and energy e(s), by use of Equations (12) and (13), respectively. The nonstandard function U(s) is likewise a monotone increasing function on I.
We consider the case of a standing hydraulic jump (c = 0) in a steady state, over a frictionless horizontal bed, as illustrated in Figure 2. It can be observed that the water depth and the flow velocity are monotonically increasing and decreasing, respectively, whereas the energy flux displays a sharp trough inside the jump. Furthermore, this undershoot does not affect the overall energy loss across the jump (cf. Equation (15)). The Froude number at the location of this undershoot (s = 0.5) turns out to be exactly 1.
Overall, the present observations confirmed the results from the analysis of Section 3.3 and demonstrated that the undershoot observed in the energy profile inside the hydraulic jump is rooted in the continuity of mass flux across the shock because the microstructure associated with the water depth and the flow velocity is distinct, viz. Equation (12). Under the assumption of a strictly increasing function for the water depth in the microstructure (0, ε), this guarantees the existence of the entropy-consistent minimum value of energy in the shock.

Numerical Experiment
Gharangik and Chaudhry [14] reported both laboratory and numerical experiments to investigate the formation of a hydraulic jump in a rectangular flume with metal walls and bed. These include six test conditions with the upstream Froude number F u = u u / gh u ranging from 2.3 to 7.0, and thus provide a useful benchmark solution for validating numerical methods. The objective of this test case is to assess the numerical efficiency of the staggered C-grid discretization as outlined in Section 2.2. Its complete description and implementation is given in Appendix A.
Numerical results of the current method are presented for two conditions: F u = 2.30 (test no. 6 of [14]) and F u = 5.74 (test no. 3). The flume dimensions and flow parameters for these cases are identical to those used by Gharangik and Chaudhry [14].
The computational grid was uniform and the grid size was set to 0.15 m. At the upstream boundary, an in-going Riemann invariant, defined as u u + 2 gh u , was imposed, whereas at the downstream end the water depth h d was fixed. The time step was initially taken as 0.01 s and the Courant number was set equal to 0.5 (see Appendix A). The simulation time was set to 90 s, which was long enough to get a stabilized jump. Following Gharangik and Chaudhry [14], a Manning's bed roughness coefficient of k = 0.008 was adopted. No artificial viscosity nor flux limiting was applied in the course of simulations.
Numerical predictions and comparisons with the experimental data are depicted in Figure 4.  Black solid line: present method using Fromm's scheme; black dash line: present method using first order upwind scheme; red dots: experimental data of Gharangik and Chaudhry [14].
Here, the momentum advection is approximated by means of two upwind schemes: the first order upwind scheme (see Appendix A) and the second order Fromm's scheme [44]. The latter scheme is designed for low numerical dispersion and dissipation.
Obviously, model-predicted jump profiles agree very well with the measurements. The location of the steady-state jump is adequately captured by the numerical approach. Additionally, the steepness of the jump profile is hardly affected by the amount of numerical diffusion induced by an advection scheme. This suggests that the role of this diffusion in the formation of the jump is less relevant than commonly believed. This finding is in line with the conclusion of Section 3.3.
To verify the correctness of the simulations the well-known Bélanger formula is employed and is given by Based on the numerical outcomes, we observed a 5% lower ratio h d /h u in comparison with Equation (16), which is acceptable.
It is further observed that there are virtually no spatial oscillations, albeit a small overshoot produced by the Fromm's scheme. It should be noted that the numerical methods employed by Gharangik and Chaudhry [14] necessarily involve some amount of artificial viscosity to suppress these oscillations, and possibly numerical instability. The same conclusion is also found in [16]. The main cause of such oscillations is due to the nonstaggered arrangement of the unknowns, which induces a checkerboard pattern. In turn, this can lead to a loss of continuity of the mass flux. Lastly, the study of Ginting [16] has also demonstrated the influence of the tuning of artificial viscosity on the jump location.

Conclusions
This paper provided an explanation of the ability of staggered C-grid discretizations to produce physically realistic results for rapidly varying free-surface flows by leveraging the concepts of nonstandard calculus. It was further argued that the physically consistent nonstandard equations governing the conservation of water volume and the streamwise transport of momentum can be established as a numerical methodology for solving the inviscid shallow water equations on staggered grids. The resulting discretization is characterized by the topological structures embedded in the nonstandard equations and the primary unknowns.
Nonstandard analysis also opens the door to adequately analyze nonlinear problems in the presence of discontinuities. The distributions of the water depth, flow velocity and mechanical energy across a one-dimensional infinitesimal thin inviscid hydraulic jump were analyzed by means of nonstandard Heaviside functions. The analysis presented in this work demonstrated that the modeled microstructure of the energy profile displays an undershoot inside the hydraulic jump, which is the result of the incompressibility constraint while the viscosity, either physical or artificial, plays no role. This principal result concludes that the staggered finite difference equations derived in this paper are effective for regularizing the entropy-consistent solution to the nonlinear shallow water equations.
The numerical results reported in this paper serve to illustrate that the numerical approach based on C-grid staggering indeed provides a superior representation of hydraulic jumps to that returned by the methods using artificial viscosity.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
We present a staggered finite difference scheme in the one-dimensional case. In the following the space domain is divided into M grid cells of length ∆x with cell centers located at {x m | m = 1, . . . , M} and cell interfaces at {x m+1/2 | m = 0, . . . , M}; see Figure A1. With the staggered grid arrangement, the primitive variables are carried at alternate grid points. The water depth is located at the center of the cells, whereas the flow velocity is resided at the cell interfaces. Likewise, updates of these unknowns are staggered in time: the water depth is evaluated at each whole time step n∆t, whereas the velocity at each half time step (n + 1/2)∆t, with ∆t the time step and n = 0, . . . , N indicating the time level t n = n∆t. The associated approximations are denoted by h n m and u n+1/2 m+1/2 , respectively. Concerning the solution strategy, the nonstandard Equations (4) and (6) are solved by adopting the method described in Zijlema [22]. The underlying approach uses a mix of interpolation and upwinding techniques with the aim to design a staggered scheme which satisfies the Rankine-Hugoniot jump conditions. These conditions express the continuity of mass flux and momentum flux across a discontinuity and are necessary for preventing the odd-even decoupling problem. Below, we briefly outline the numerical procedure. For further details, which also include a detailed implementation of the varying bed topography, we refer to [22,34].
We first solve the momentum equation as follows: where h n m+1/2 = Second order accuracy can be obtained by means of a high order upwind scheme. In this paper, we employ the well-known second order Fromm's scheme, which generates a limited amount of numerical dispersion and dissipation [44].
The water depth at the cell face in Equation (A2) is approximated as follows: This approximation guarantees a non-negative water depth if the time step is limited according to max u n−1/2 m+1/2 , 0 − min u n−1/2 m−1/2 , 0 ∆t ∆x ≤ 1 .
To model the effect of bed friction on the flow motion, the following sink term is added to the left hand side of Equation (A1) k 2 g | u n−1/2 m+1/2 | (h n m+1/2 ) 1/3 u n+1/2 m+1/2 with k the Manning's roughness coefficient.
After the momentum Equation (A1) is solved, the continuity equation is then solved, which is given by With respect to the temporal order of accuracy, the present scheme is second order accurate in time in the absence of the advection term by virtue of staggering in time. To obtain second order accuracy in both space and time, a predictor-corrector technique is adopted. In the predictor step, the momentum Equation (A1) is solved to determine the provisional flow velocities u * m+1/2 . Next, in the corrector step, these values are corrected by means of the Fromm's scheme, while marching in time using the MacCormack approach to achieve second order accuracy in time. See [34] for further details. Convergence tests of the present method are demonstrated in [22].
Since the time integration is explicit, the following CFL criterion | u n−1/2 m+1/2 | + gĥ n m+1/2 ∆t ∆x ≤ 1 is employed to determine the maximum time step for stability. Note that this time step automatically satisfies condition (A3).