Hyperbolic conservation laws on the sphere.
A geometry-compatible finite volume scheme
We consider entropy solutions to the initial value problem associated with scalar nonlinear hyperbolic conservation laws posed on the two-dimensional sphere. We propose a finite volume scheme which relies on a web-like mesh made of segments of longitude and latitude lines. The structure of the mesh allows for a discrete version of a natural geometric compatibility condition, which arose earlier in the well-posedness theory established by Ben-Artzi and LeFloch. We study here several classes of flux vectors which define the conservation law under consideration. They are based on prescribing a suitable vector field in the Euclidean three-dimensional space and then suitably projecting it on the sphere’s tangent plane; even when the flux vector in the ambient space is constant, the corresponding flux vector is a non-trivial vector field on the sphere. In particular, we construct here “equatorial periodic solutions”, analogous to one-dimensional periodic solutions to one-dimensional conservation laws, as well as a wide variety of stationary (steady state) solutions. We also construct “confined solutions”, which are time-dependent solutions supported in an arbitrarily specified subdomain of the sphere. Finally, representative numerical examples and test-cases are presented.
keywords:hyperbolic conservation law, sphere, entropy solution, finite volume scheme, geometry-compatible flux.
We propose a Godunov-type finite volume scheme that satisfies certain important consistency and convergence properties. We then present a second-order extension based on the generalized Riemann problem (GRP) methodology .
It should be stated at the outset that an important motivation for this paper is the need to provide accurate numerical tools for the so-called shallow water system on the sphere. This system is widely used in geophysics as a model for global air flows on the rotating Earth . In its mathematical classification it is a system of nonlinear hyperbolic PDE’s posed on the sphere. Its physical nature dictates that it can be described “invariantly”, namely in a way which is independent of any particular coordinate system. Locally, it has the (mathematical) character of a two-dimensional isentropic compressible flow, whereas globally the spherical geometry plays a crucial role in shaping the nature of solutions –which, as expected for nonlinear hyperbolic equations, may contain propagating discontinuities such as shock fronts or contact curves. Thus, the relation of the present study to the shallow water system is analogous to the connection between Burgers’ equation and the system of compressible fluid flow (say, in the plane). In fact, in light of this analogy it is somewhat surprising that in the existing literature so far, virtually all treatments, theoretical as well as numerical, were confined to the Cartesian setting. In particular, to the best of our knowledge, there have been no systematic numerical studies of scalar conservation laws on the sphere.
Having introduced the scalar conservation law as a simple model for more complex physical systems, we should emphasize here also the intrinsic mathematical interest of the model under consideration. It is already known (see  and references there) that even in the Cartesian setting, the two-dimensional scalar conservation law displays a wealth of wave interactions typical of the physical phenomena (such as triple points, sonic shocks, interplay of rarefactions and shocks coming from different directions and more). As we show here, “geometric effects”, superposed on the (necessarily) two-dimensional framework, carry the scalar model still further. For example, the concept of “self-similar” solutions makes no sense here. In particular, one loses the Riemann solutions, a fundamental building block in many schemes (of the so-called “Godunov-type”). On the other hand, it allows for large classes of non-trivial steady states, periodic solutions, and solutions supported in specified subdomains. All these have natural consequences in developing numerical schemes; they offer us a variety of test-cases amenable to detailed analysis, to be compared with the computational results.
In practical applications a finite volume scheme requires a specification of a coordinate system, where the symmetry-preserving latitude–longitude coordinates are the “natural coordinates” of preferred choice. The proposed finite volume scheme in this paper is based on these natural coordinates, but should pay attention to the artificial singularities at the poles.
In , a general convergence theorem was proved for a class of finite volume schemes for the computation of entropy solutions to conservation laws posed on a manifold. As a particular example, the case of the sphere was discussed, both from the points of view of an “invariant” formalism and that of an “embedded” coordinate-dependent formulation. In the present study we focus on the sphere and we actually construct, in a fully explicit and implementable way, a finite volume scheme which is geometrically natural and can be viewed as an extension of the basic Godunov scheme for one-dimensional conservation laws. Furthermore, we prove that our scheme fulfills all of the assumptions required in , which ensures its strong convergence toward the unique entropy solution to the initial value problem under consideration. We then describe the GRP extension of the scheme, whose convergence proof is still a challenging open problem.
The theoretical background about the well-posedness theory for hyperbolic conservation laws on manifolds was established recently by Ben-Artzi and LeFloch  together with collaborators [1, 2, 8]. An important condition arising in the theory is the “zero-divergence” or geometric-compatibility property of the flux vector; a basic requirement in our construction of a finite volume scheme is to formulate and ensure a suitable discrete version of this condition.
We conclude this introduction with some notation and remarks connecting the present paper to the general finite volume framework presented in . Following the terminology therein, we use an “embedded” approach to the spherical geometry, namely, we view the sphere as embedded in the three-dimensional Euclidean space . We denote by a variable point on the sphere , which can be represented in terms of its longitude and its latitude . Following the conventional notation in the geophysical literature we assume that
so that the “North pole” (resp. “South pole”) is at (resp. ) and the equator is (See Figure 1.) The coordinates in are denoted by and the corresponding unit vectors are . Thus, at each point , the unit tangent vectors (in the directions) are given by
It should be observed that while a choice of a coordinate system is necessary in practice, it always introduces singularities and the unit vectors given above are not well-defined at the poles and, therefore, in the neighborhood of these points it cannot be used for a representation of smooth vector fields (such as the flux vectors of our conservation laws). We also emphasize that the status of these two poles is equivalent to the one of any other pair of opposite points on the sphere. When such local coordinates are introduced, special care is needed to handle these points in practice, and this is precisely why we advocate a different approach.
Continuing with the description of our “embedded” approach, we define the unit normal , to at some point by
Then, any tangent vector field to is represented by
and the tangential gradient operator is
Thus, the (tangential) gradient of a scalar function is given by
and the divergence of a vector field is
Given now a vector field depending on a real parameter , the associated hyperbolic conservation law under consideration is
where is a scalar unknown function, subject to the initial condition
for some prescribed data on the sphere. As mentioned above, we will impose on the vector field an additional “geometry compatibility” condition.
An outline of this paper is as follows. In Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme, we consider the construction of geometry-compatible flux vectors, while Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme is devoted to a description of several families of special solutions associated with the constructed flux vectors. In Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme we discuss our (first-order) finite volume scheme, which can be regarded as a Godunov-type scheme. We prove that it satisfies all of the assumptions imposed on general finite volume schemes in , and we conclude that it converges to the exact (entropy) solution. In Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme we describe the (second-order) GRP extension of the scheme. Finally, in Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme we present a variety of numerical test cases.
As pointed out in , every smooth vector field on can be represented in the form
where is a restriction to of a vector field (in ) defined in some neighborhood (i.e., a “spherical shell”) of and for all values of the parameter . The basic requirement imposed now on the flux vector is the following divergence free or geometric compatibility condition: For any fixed value of the parameter ,
A flux vector satisfying (2.2) is called a geometry-compatible flux . Note that this condition is equivalent, in terms of the nonlinear conservation law (1.3), to the following requirement: constant initial data are (trivial) solutions to the conservation law. In the case of the sphere the condition (2.2) can be recast in terms of a condition on the vector field appearing in (2.1). See [2, Proposition 3.3].
Our main aim in the present section is singling out two (quite general) families of geometry-compatible fluxes of particular interest, which are amenable to detailed analytical and numerical investigation.
The flux-vectors of interest are introduced by way of the following two claims.
Claim 2.1 (Homogeneous flux vectors.)
If the three-dimensional flux is independent of (in a neighborhood of ), then the corresponding flux vector given by (2.1) is geometry-compatible.
Proof. The following decomposition applies to any vector in the form
so that , with
We can directly apply the divergence operator (1.2) to and the desired claim follows. ∎
Claim 2.2 (Gradient flux vectors.)
Let be a smooth function of the variables (in a neighborhood of ) and , and consider the associated three-dimensional flux (restricted to ). Then, the flux vector given by (2.1) is geometry-compatible.
Proof. We use the divergence theorem in an arbitrary domain with smooth boundary :
where is the unit normal (at ) along , is the surface measure on , and is the arc length along
In particular, coincides with the (unit) tangent vector to at . It follows that the triple product is nothing but the directional derivative of along . Since
we thus find
and since this holds for any smooth domain , we conclude that for all . ∎
1. Claim 2.1 is a special case of Claim 2.2. Indeed, by taking in the latter we obtain the conclusion of the former. However, we chose to single out Claim 2.1 as a special case since it will serve in obtaining special solutions (Section 3) and in dealing with numerical examples (Section 6).
The scalar conservation laws discussed in this paper have two basic features:
The problem is necessarily two-dimensional (in spatial coordinates).
The geometry plays a significant role, inasmuch as the flux vectors are subject to geometric constraints.
It should be noted that even within the framework of Euclidean two dimensional conservation laws there is a great wealth of special solutions, displaying complex wave interactions, such as triple points, sonic shocks and more. We refer to [9, 5] for detailed treatments of the theoretical and numerical aspects.
In the situation under consideration in the present paper, geometric effects yield a large variety of non-trivial steady states, solutions supported in arbitrary subdomains, etc. In this section we consider such solutions by selecting some special flux vectors on . This is accomplished by making special choices of in the general representation (see (2.1)) , where is a restriction to of a vector field (in ) defined in some neighborhood (i.e., “spherical shell”) of and for all values of the parameter .
the conservation law (1.3) takes the particularly simple form
In particular, obtain the following important conclusion.
Corollary 3.1 (Solutions with one-dimensional structure.)
Let be a solution to the following one-dimensional conservation law with periodic boundary condition
and let be an arbitrary function. Then, the function is a solution to the conservation law (3.1).
It follows that all periodic solutions from the one-dimensional case can be recovered here as special cases. However, in numerical experiments the computational grid is two-dimensional, so it is not obvious that the accuracy achieved in the computation of the former can indeed be achieved in the numerical scheme implemented on the sphere. This issue will be further discussed below, in Section 6.
Let be a flux vector and be an initial function such that . Then, clearly is a stationary solution (or steady state) to the conservation law. In fact, we can show that there exist many (analytically computable) non-trivial steady state solutions, as follows.
Claim 3.2 (A family of steady-state solutions.)
Let be a smooth function defined for all in a neighborhood of , and consider the associated gradient flux vector (as in Claim 2.2). Suppose the function satisfies the condition
where be a smooth function defined in a neighborhood of . Then, is a stationary solution to the conservation law (1.3).
Proof. We follow the proof of Claim 2.2 and the notation therein. Using the divergence theorem in an arbitrary domain with smooth boundary , we obtain
where, as before, is the unit normal, the surface measure, and the arc length. In particular, the (unit) tangent vector to at It follows that the triple product is the directional derivative of along Thus,
and since this holds for any smooth domain , it follows that which concludes the proof. ∎
The above claim yields readily a large family of non-trivial stationary solutions, as expressed in the following corollary.
Consider the flux vector given by
for an arbitrary choice of function . Then, any function depending only on the first coordinate is a stationary solution to the conservation law (associated with this flux). In particular, in polar coordinates any function of the form is a stationary solution.
This corollary enables us to construct stationary solutions supported in “bands” on the sphere. This is accomplished by taking to be supported in Observe that this band is not parallel neither to the latitude curves () nor to the longitude curves ().
There is yet another possibility of obtaining stationary solutions, where all three coordinates are involved, as stated now. This example can also be derived from the previous one by applying a rotation in .
Consider the flux vector be given by
in which all three components coincide: . Then, any function of the form , where depends on one real variable, only, is a stationary solution to the conservation law associated with the above flux.
Proof. Following the proof of the previous corollary, we now take where . ∎
In analogy with Remark 3.4, this result allows us to construct stationary solutions in a spherical “cap” (a piece of the sphere cut out by a plane). In Section 6 below, we will provide numerical test cases for such stationary solutions.
If in the conservation law (1.3) we have for in the exterior of some domain identically in and if the initial function vanishes outside of , then clearly the solutions satisfy for and all We label such solutions as confined (to ) solutions. In view of equation (2.1) a sufficient condition for the vanishing of outside of is obtained by for identically in In view of Claim 2.2, this will follow if we choose such that for only in In particular, let be a twice continuously differentiable function on supported in the interval and such that and With an eye to computable test cases, we can use this function to generate solutions which are confined within the intersection of with the (three-dimensional) cube .
Claim 3.7 (A family of confined solutions.)
Let be as above and let be any (smooth) function of . Define by
and let be the gradient flux vector determined in terms of as in Claim 2.2. Let be the spherical patch cut out from by the inequalities Then, if the initial data is supported in the solution of the conservation law (1.3) associated with is supported in for all
Possible choices for a function as in the claim are for some integer such that and are multiples of , or else .
The general structure of our grid is shown in Figure 1, and its essential feature is the following. Every cell is bounded by sides which lie either along a fixed latitude circle () or a fixed longitude circle (). We have
as represented in Figure 2. In most cases, consists of the four sides of . However, across special latitude circles we reduce the number of cells, so that the situation (for a reduction by ratio of ) is as in Figure 3. In this case the boundary consists of five sides, (so that the intermediate point is regarded as an additional vertex), and even in this five-sided cell every side satisfies the above requirement.
The length of a side equals , while the length of a side is . Consequently, the area of the cell is
Given any rectangular domain of the form (4.1), the approximate flux divergence is now derived as an approximation of the integral of the flux along the boundary , divided by its area, as follows:
We need to check that the geometric compatibility condition (2.2) is satisfied for the approximate flux divergence. This requirement will be taken into account in formulating our finite volume scheme for (1.3).
Consider now the actual evaluation of the term defined in (4.2) and consider the cell shown in Figure 2, under the assumption that is smooth on . We propose to approximate the flux integral along each edge of in the following way. As in Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme, let us decompose the flux into its components:
On each side the integration is carried out by (i) taking midpoint values of the appropriate flux component, and (ii) using the correct arc-length of the side. We designate the midpoints of the edge as and (see Figure 2), and likewise for the edge .
Throughout the rest of this section we restrict attention to the gradient flux vector constructed in Claim 2.2. In particular, it comprises the class of homogeneous flux vectors, given by (2.3)–(2.4).
Taking as constant along the side , the total approximate flux is given by
where are, respectively, the initial and final endpoints of (with respect to the sense of the integration).
Summing up over all edges we obtain:
Claim 4.1 (Discrete geometry-compatibility condition.)
The claim above applies to gradient flux vectors in Claim 2.2, and, in particular, to homogeneous flux (2.3)–(2.4). On the other hand, for a more general geometry-compatible flux , such a result can be obtained only if the dependence on is integrated exactly along each side, a requirement that must be imposed on the scheme.
We continue to deal with the gradient flux given in Claim 2.2. We assume different (constant) values of in grid cells and evaluate the numerical flux values at each edge from the solution to a Riemann problem with data comprising these values in the cells on either side of that edge. At the midpoint of each side we solve the Riemann problem in a direction perpendicular to , and denote the resulting solution . The corresponding fluxes are then evaluated as .
Consider two adjacent cells, as in Figure 4 or in Figure 5. By fixing (resp. ) in (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)(resp. (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)) we can evaluate as a one-dimensional solution at (resp. ).
We include here some remarks that will be useful in the implementation of the scheme.
Consider an homogeneous flux vector as in Claim 2.1 so that its components are given by (2.4). Suppose that (resp. ) in the cell (resp. ), as in Figure 4. At the point Eq. (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)takes the form
we see that equation (4.5) is the scalar one-dimensional conservation law
subject to the initial data (resp. ) for (resp. ).
Likewise, we repeat the former analysis for -adjacent cells by taking the constant states , in cells , , as depicted in Figure 5. At the point , the equation (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)then takes the form
We then set the -flux function
so that equation (4.8) is the scalar one-dimensional conservation law
subject to the initial data (resp. ) for (resp. ).
The solution at the discontinuity at the initial time is given by the Riemann solution to (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme). For simplicity of the presentation we specialize here to the flux (4.7). Since the dependence of on is smooth, this solution is obtained by fixing , thus solving the classical conservation law
subject to the initial jump discontinuity of .
We denote this solution by . Observe that the flux in (4.11) is in general non-convex. The Riemann solution may therefore consist of several waves. It is a self-similar solution depending only on the slope . The value is the value along the line . It therefore corresponds either to a sonic wave, namely , or to an “upwind value” (resp. ) in the case where all waves propagate to the right (resp. left).
Actually, the procedure for solving the Riemann problem in the case of a nonconvex flux function is well-known and goes back to classical works by Oleinik and others. We recall it here briefly. Assume first that . Consider the convex envelope of , namely, the largest convex continuous function , over the interval , such that at all points. Clearly, in “convex sections” of the graph of , while it consists of linear segments when . It is easy to see that the “convex segments”, where , represent rarefaction waves (in the full Riemann solution) while the linear segments represent jumps (i.e., shock waves). In particular, the solution is given by the following formula:
There are in fact three possibilities for this solution:
, which implies that (a sonic point).
, the whole wave pattern moves to the right.
, the whole wave pattern moves to the left.
Similarly, in the case , we construct the “concave envelope” of , namely, the smallest concave continuous function such that . Again the linear segments correspond to jump discontinuities while the concave segments () correspond to rarefaction waves. The solution to the Riemann problem is now given by , where , . As above, there are three possibilities for the solution (sonic, left-upwind, or right-upwind).
Replacing in the foregoing analysis the -flux function by the -flux function , the equation (4.10) reads
We get the Riemann solution to (4.13) in the three cases a), b), c) as above.
The computational elements (“grid cells”) are denoted in  by . Their sides are denoted by and the flux function across is given by where is the (constant) value in and is the value in the neighboring cell (sharing the same side ) In our grid of the sphere, some cells are actually pentagons; these are the cells whose lower-latitude side (along a latitude ) borders the two higher-latitude sides of the two lower-latitude neighbor cells, as shown in Figure 3 for the southern hemisphere grid. For such cells, the lower-latitude side consists of two faces, each one of them common with one of the lower-latitude neighboring cells.
With this construction of the grid, we can check the conditions in  imposed on the numerical flux. It is important to keep in mind that we are dealing with the gradient flux vectors given by Claim 2.2.
Claim 4.3 (Convergence of the proposed scheme.)
Consider the first-order finite volume scheme described above. Assume that the flux vector has the gradient form in Claim 2.2. Let be the numerical flux calculated on the side of the computational cell using (4.3), where the midpoint value of is obtained from the Riemann solution. Then satisfies the assumptions (5.5)-(5.7) of , and the numerical solution converges to the exact solution as the maximal size of the grid cells shrinks to zero.
Proof. Consider the flux across a longitude side , which is given by in the equation (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme). The procedure for integrating the flux across is described by (4.3), while in Subsection Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme the calculation of is described. It can be summarized as follows.
First, the solution to the Riemann problem associated with equation (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)is found, assuming to be the values on the two sides. However, note that depends explicitly on , and to be precise we need to replace in (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)the mean value by Thus, we find
Clearly, in the case we get identically and so the exact flux satisfies
and its integration will give exactly the approximate value
Clearly, the conservation property (5.6) is satisfied even with the approximate definition.
Also, the flux as defined in (4.3) makes it easy to check (5.7), as the flux is independent of and the monotonicity is thus a result of general properties of the Riemann solver (even for nonconvex fluxes). For example, if , one considers the convex envelope of , as defined in (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)(with ) and then considers as the minimal value on this envelope (over ). Clearly changing upward will either change upward or leave it unchanged. This completes the proof. ∎
To improve the order of accuracy, we consider again the cell and assume that is linearly distributed there. We use (resp. ) to denote the slopes in the cell to the left (resp. right) of the side . We also denote by (resp. ) the limiting value (linearly distributed) of at (resp. ). Clearly, the solution to the Riemann problem across the discontinuity is a function of , and we denote it by , which conforms to our notation in Subsecion Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme above (where was constant on either side of the discontinuity). The value of is obtained by solving the Riemann problem associated with Eq. (Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme)with replaced by , subject to the initial data . Restricting to the middle point , the solution (at ) is in one of the three categories listed above (i.e., sonic, left-upwind, right-upwind). By continuity, the solution will still be in the same category for sufficiently small. The solution at varies in time and the GRP method deals with the determination of its time-derivative at that point.
Accounting for the variation of the solution over a time interval enables us to modify the Godunov approach to the determination of edge fluxes , as presented in Section Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme. We assume that the flux vector depends explicitly on as in (2.1). In what follows we use for simplicity the “imbedded” notation for a point on the sphere (see the Introduction), along with the corresponding spherical coordinates We further assume that the vector field is given by the following extension of (2.3)
The zero-divergence identity is obtained as a result of expressing as a gradient in the sense of Claim 2.2.
For our choice of such a representation of as gradient of is obtained when is taken as
The numerical approximation to this equation requires an operator splitting approach,
where the derivatives with respect to and are considered separately.
We note that such a splitting has already been implemented in the Godunov case,
(Hyperbolic conservation laws on the sphere.
A geometry-compatible finite volume scheme), in the most general case. In that case, no use has been
made of the geometry-compatibility property. Indeed, this has no bearing on the
first-order scheme since the solution to the Riemann problem is obtained by
“freezing” the explicit dependence on (and, in particular, ignoring
the terms involving the derivatives with respect to this explicit dependence).
In the present (second-order) situation we proceed as follows.
The “-split” equation obtained from (5.3), is
Note that the coefficients are retained as functions of and are not “frozen” at . This is of course due to the fact that in employing the GRP scheme we consider -derivatives on either side of the edge, so as in any limiting analysis, we must first let , then substitute .
The -edge flux function (compare (4.6)), is now extended to as
and the scalar one-dimensional conservation law under consideration is now rewritten as an equation with a source term (a balance law)
subject to the initial data (for and its slope) , (resp. , ) for (resp. ). Observe that the equation is written in a “quasi-conservative form”, which offers more convenience in the GRP treatment [3, Chap. 5]. The right-hand side term is just the result of the differentiation of the flux Obviously, the geometry-compatibility condition implies that this source term should cancel out with the corresponding source term in the “-split” equation. The solution to the Riemann problem is obtained by freezing the coordinate at its edge value, so that, in particular, the source term in (5.5) can be taken as zero at this stage.
In the framework of the GRP analysis, the source term is added to terms arising from the piecewise-linear initial data, in producing the time-derivative of the solution . As explained above, , the solution to the associated Riemann problem, is obtained by using the “edge values” . It remains, therefore, to determine the instantaneous time-derivative , as is outlined below.
The time-derivative of is given by