Matania Ben-Artzi

Joseph Falcovitz

Philippe G Lefloch
email: pglefloch@gmail.com

Hyperbolic conservation laws on the sphere. A geometry-compatible finite volume scheme

Keywords: hyperbolic conservation law, sphere, entropy solution, finite volume scheme, geometry-compatible flux PACS: 35L65, 76L05

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 onedimensional 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 timedependent solutions supported in an arbitrarily specified subdomain of the sphere. Finally, representative numerical examples and test-cases are presented.

Introduction

In this paper, building on our earlier analysis in [START_REF] Ben-Artzi | The well-posedness theory for geometry-compatible hyperbolic conservation laws on manifolds[END_REF][START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF] we study in detail the class of scalar hyperbolic conservation laws posed on the two-dimensional unit sphere

S 2 = (x, y, z) ∈ R 3 , x 2 + y 2 + z 2 = 1 .
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 [START_REF] Ben-Artzi | Generalized Riemann problems in computational fluid dynamics[END_REF].

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 [START_REF] Haltiner | Numerical weather prediction[END_REF]. 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 [START_REF] Ben-Artzi | Wave interactions and numerical approximation for twodimensional scalar conservation laws[END_REF] 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 [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF], 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 S 2 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 S 2 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 [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF], 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 [START_REF] Ben-Artzi | The well-posedness theory for geometry-compatible hyperbolic conservation laws on manifolds[END_REF] together with collaborators [START_REF] Amorim | Finite volume schemes on Lorentzian manifolds[END_REF][START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF][START_REF] Lefloch | Hyperbolic conservation laws on manifolds. Error estimate for finite volume schemes[END_REF][START_REF] Lefloch | Hyperbolic conservation laws on spacetimes. A finite volume scheme based on differential forms[END_REF]. An important condition arising in the theory is the "zerodivergence" 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 [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF]. 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 R 3 . We denote by x a variable point on the sphere S 2 , which can be represented in terms of its longitude λ and its latitude φ. Following the conventional notation in the geophysical literature we assume that 0

≤ λ ≤ 2π, - π 2 ≤ φ ≤ π 2 ,
so that the "North pole" (resp. "South pole") is at φ = π 2 (resp. -φ = π 2 ) and the equator is φ = 0, 0 ≤ λ ≤ 2π . (See Figure 1.) The coordinates in R 3 are denoted by (x 1 , x 2 , x 3 ) ∈ R 3 and the corresponding unit vectors are i 1 , i 2 , i 3 . Thus, at each point x = (λ, φ) ∈ S 2 , the unit tangent vectors (in the λ, φ directions) are given by

i λ = -sin λ i 1 + cos λ i 2 , i φ = -sin φ cos λ i 1 -sin φ sin λ i 2 + cos φ i 3 .
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 n x , to S 2 at some point x by

n x = cos φ cos λ i 1 + cos φ sin λ i 2 + sin φ i 3 .
Then, any tangent vector field F to S 2 is represented by

F = F λ i λ + F φ i φ
and the tangential gradient operator is

∇ T = 1 cos φ ∂ ∂λ , ∂ ∂φ .
Thus, the (tangential) gradient of a scalar function h(λ, φ) is given by

∇ T h = 1 cos φ ∂h ∂λ i λ + ∂h ∂φ i φ , (1.1)
and the divergence of a vector field F is

∇ T • F = 1 cos φ ∂ ∂φ (F φ cos φ) + ∂ ∂λ F λ . (1.2)
Given now a vector field F = F(x, u) depending on a real parameter u, the associated hyperbolic conservation law under consideration is

∂u ∂t + ∇ T • F(x, u) = 0, (x, t) ∈ S 2 × [0, ∞), (1.3)
where u = u(x, t) is a scalar unknown function, subject to the initial condition

u(x, 0) = u 0 (x), x ∈ S 2 (1.4)
for some prescribed data u 0 on the sphere. As mentioned above, we will impose on the vector field F(x, u) an additional "geometry compatibility" condition. An outline of this paper is as follows. In Section 2, we consider the construction of geometry-compatible flux vectors, while Section 3 is devoted to a description of several families of special solutions associated with the constructed flux vectors. In Section 4 we discuss our (first-order) finite volume scheme, which can be regarded as a Godunovtype scheme. We prove that it satisfies all of the assumptions imposed on general finite volume schemes in [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF], and we conclude that it converges to the exact (entropy) solution.

In Section 5 we describe the (second-order) GRP extension of the scheme. Finally, in Section 6 we present a variety of numerical test cases.

Families of geometry-compatible flux vectors

As pointed out in [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF], every smooth vector field F(x, u) on S 2 can be represented in the form

F(x, u) = n(x) × Φ(x, u), (2.1)
where Φ(x, u) is a restriction to S 2 of a vector field (in R 3 ) defined in some neighborhood (i.e., a "spherical shell") of S 2 and for all values of the parameter u. The basic requirement imposed now on the flux vector F(x, u) is the following divergence free or geometric compatibility condition: For any fixed value of the parameter v ∈ R,

∇ T • F(x, v) = 0. ( 2

.2)

A flux vector F(x, u) satisfying (2.2) is called a geometry-compatible flux [START_REF] Ben-Artzi | The well-posedness theory for geometry-compatible hyperbolic conservation laws on manifolds[END_REF]. 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 S 2 the condition (2.2) can be recast in terms of a condition on the vector field Φ(x, u) appearing in (2.1). See [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF]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 Φ(x, u) = Φ(u) is independent of x (in a neighborhood of S 2 ), then the corresponding flux vector F(x, u) given by (2.1) is geometry-compatible.

Proof. The following decomposition applies to any vector Φ(u) ∈ R 3 in the form

Φ(u) = f 1 (u) i 1 + f 2 (u) i 2 + f 3 (u) i 3 , (2.3) so that F(x, u) = F λ (λ, φ, u) i λ + F φ (λ, φ, u) i φ , with F λ (λ, φ, u) = f 1 (u) sin φ cos λ + f 2 (u) sin φ sin λ -f 3 (u) cos φ, F φ (λ, φ, u) = -f 1 (u) sin λ + f 2 (u) cos λ. (2.4)
We can directly apply the divergence operator (1.2) to F(x, u) and the desired claim follows. 2

Claim 2.2 (Gradient flux vectors.) Let h = h(x, u) be a smooth function of the variables x (in a neighborhood of S 2 ) and u ∈ R, and consider the associated threedimensional flux Φ(x, u) = ∇h(x, u) (restricted to x ∈ S 2 ). Then, the flux vector F(x, u) given by (2.1) is geometry-compatible.

Proof. We use the divergence theorem in an arbitrary domain D ⊆ S 2 with smooth boundary ∂D:

D ∇ T • F(x, v) dσ = ∂D F(x, v) • ν(x) ds = ∂D n(x) × ∇h(x, v) • ν(x) ds,
where ν(x) is the unit normal (at x) along ∂D ⊂ S 2 , dσ is the surface measure on S 2 , and ds is the arc length along ∂D.

In particular, n(x) × ν(x) = t(x) coincides with the (unit) tangent vector to ∂D at x.

It follows that the triple product n(x) × ∇h(x, u) • ν(x) = ∇h(x, u) • t(x) is nothing but the directional derivative ∇ ∂D of h along ∂D. Since ∂D ∇ ∂D h ds = 0, we thus find D ∇ T • F(x, u) dσ = 0,
and since this holds for any smooth domain D, we conclude that

∇ T • F(x, v) = 0 for all v ∈ R. 2 Remark 2.3 1. Claim 2.1 is a special case of Claim 2.2. Indeed, by taking in the latter h(x, u) = x 1 f 1 (u) + x 2 f 2 (u) + x 3 f 3 (u)
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).

2. The steps in the construction of the gradient flux vector in Claim 2.2 are "linear in nature", namely if h(x, u) = h 1 (x, u) + h 2 (x, u) then the corresponding (geometrycompatible) flux vectors satisfy F(x, u) = F 1 (x, u)+F 2 (x, u). However, it is clear that the corresponding solutions to (1.3) do not add up linearly, due to the nonlinear dependence in u.

Special solutions of interest

Periodic equatorial solutions

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 [START_REF] Li | The two-dimensional Riemann problem in gas dynamics[END_REF][START_REF] Ben-Artzi | Wave interactions and numerical approximation for twodimensional scalar conservation laws[END_REF] 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 F(x, u) on S 2 . This is accomplished by making special choices of Φ(x, u) in the general representation (see (2.1)) F(x, u) = n(x) × Φ(x, u), where Φ(x, u) is a restriction to S 2 of a vector field (in R 3 ) defined in some neighborhood (i.e., "spherical shell") of S 2 and for all values of the parameter u.

We begin our discussion with the case of periodic equatorial solutions, defined as follows. Taking f 1 (u) = f 2 (u) ≡ 0 in the general decomposition (2.3) so that, by (2.4),

F λ (λ, φ, u) = -f 3 (u) cos φ, F φ (λ, φ, u) = 0, the conservation law (1.3) takes the particularly simple form ∂u ∂t - ∂ ∂λ f 3 (u) = 0, (x, t) ∈ S 2 × [0, ∞). (3.1)
In particular, obtain the following important conclusion.

Corollary 3.1 (Solutions with one-dimensional structure.) Let u = u(λ, t) be a solution to the following one-dimensional conservation law with periodic boundary condition

∂ u ∂t - ∂ ∂λ f 3 ( u) = 0, 0 < λ ≤ 2π, u(0, t) = u(2π, t),
and let u = u(φ) be an arbitrary function. Then, the function u(λ, φ, t) = u(λ, t) u(φ) 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 twodimensional, 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.

Steady states

Let F = F(x, u) be a flux vector and u 0 : S 2 → R be an initial function such that ∇ T • F(x, u 0 (x)) ≡ 0. Then, clearly u 0 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 h = h(x, u) be a smooth function defined for all x in a neighborhood of S 2 , and consider the associated gradient flux vector Φ = ∇h (as in Claim 2.2). Suppose the function u 0 : S 2 → R satisfies the condition

∇ y h(y, u 0 (x))| y=x = ∇ x H(x), x ∈ S 2 , (3.2)
where H = H(x) be a smooth function defined in a neighborhood of S 2 . Then, u 0 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 D ⊆ S 2 with smooth boundary ∂D, we obtain

D ∇ T • F(x, u 0 (x)) dσ = ∂D F(x, u 0 (x)) • ν ds = ∂D n(x) × ∇ x H(x) • ν(x) ds.
where, as before, ν(x) is the unit normal, dσ the surface measure, and ds the arc length.

In particular, n(x) × ν(x) = t(x), the (unit) tangent vector to ∂D at x. It follows that the triple product (n(x

) × ∇ x H(x)) • ν(x) = (∇ x H(x)) • t(x) is the directional derivative of H along ∂D. Thus, D ∇ T • F(x, u 0 (x)) dσ = 0,
and since this holds for any smooth domain D, it follows that ∇ T • F(x, u 0 (x)) ≡ 0, which concludes the proof. 2

The above claim yields readily a large family of non-trivial stationary solutions, as expressed in the following corollary.

Corollary 3.3 Consider the flux vector F = F(x, u) given by

F(x, u) = n(x) × f 1 (u) i 1 ,
for an arbitrary choice of function f 1 = f 1 (u). Then, any function u 0 = u 0 (x 1 ) depending only on the first coordinate x 1 is a stationary solution to the conservation law (associated with this flux). In particular, in polar coordinates (λ, φ) any function of the form u 0 (λ, φ) = g(cos φ cos λ) is a stationary solution.

Proof. According to Claim 2.1 this flux vector is associated with the scalar function

h(x, u) = x 1 f 1 (u). So we can invoke Claim 3.2 with H(x) = H(x 1 ) such that H (x 1 ) = f 1 (u 0 (x 1 )). 2
Remark 3.4 This corollary enables us to construct stationary solutions supported in "bands" on the sphere. This is accomplished by taking u 0 = u 0 (x 1 ) to be supported in 0 < α < x 1 < β < 1. Observe that this band is not parallel neither to the latitude curves (φ = const) nor to the longitude curves (λ = const).

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 R 3 . Corollary 3.5 Consider the flux vector F = F(x, u) be given by

F(x, u) = n(x) × (f 1 (u) i 1 + f 2 (u) i 2 + f 3 (u) i 3 ) = f (u) n(x) × (i 1 + i 2 + i 3 ),
in which all three components coincide:

f 1 (u) = f 2 (u) = f 3 (u) = f (u).
Then, any function of the form u 0 (x) = u 0 (x 1 + x 2 + x 3 ), where u 0 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

H(x) = H 0 (x 1 +x 2 +x 3 ), where H 0 (ξ) = f ( u 0 (ξ)). 2
Remark 3.6 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.

Confined solutions

If in the conservation law (1.3) we have F(x, u) ≡ 0 for x in the exterior of some domain D ⊆ S 2 , identically in u ∈ R, and if the initial function u 0 (x) vanishes outside of D, then clearly the solutions satisfy u(x, t) = 0 for x / ∈ D and all t ≥ 0. We label such solutions as confined (to D) solutions. In view of equation (2.1) a sufficient condition for the vanishing of F(x, u) outside of D is obtained by Φ(x, u) = 0 for x / ∈ D, identically in u ∈ R. In view of Claim 2.2, this will follow if we choose h(x, u) such that h(x, u) = 0 for x only in D. In particular, let ψ = ψ(ξ) be a twice continuously differentiable function on R supported in the interval (α, β) ⊆ (0, 1) and such that 3β 2 > 1 and 3α 2 < 1. With an eye to computable test cases, we can use this function to generate solutions which are confined within the intersection of S 2 with the (three-dimensional) cube [α, β] 3 . Claim 3.7 (A family of confined solutions.) Let ψ be as above and let f = f (u) be any (smooth

) function of u ∈ R. Define h = h(x, u) by h(x, u) = ψ(x 1 )ψ(x 2 ) ψ(x 3 ) f (u),
and let F(x, u) be the gradient flux vector determined in terms of h(x, u) as in Claim 2.2. Let D ⊆ S 2 be the spherical patch cut out from S 2 by the inequalities α < x i < β, i = 1, 2, 3. Then, if the initial data u 0 (x) is supported in D, the solution u = u(x, t) of the conservation law (1.3) associated with F(x, u) is supported in D for all t ≥ 0.

Possible choices for a function ψ : [α, β] → R as in the claim are ψ(ξ) = sin 2 (kξ) for some integer k such that kα and kβ are multiples of π, or else ψ(ξ) = (ξα) 2 (ξβ) 2 .

Design of the scheme

Computational grid

The general structure of our grid is shown in Figure 1, and its essential feature is the following. Every cell R is bounded by sides which lie either along a fixed latitude circle (φ = const.) or a fixed longitude circle (λ = const.). We have

R := λ 1 ≤ λ ≤ λ 2 , φ 1 ≤ φ ≤ φ 2 , (4.1)
as represented in Figure 2. In most cases, ∂R consists of the four sides of R. However, across special latitude circles we reduce the number of cells, so that the situation (for a reduction by ratio of 2) is as in Figure 3. In this case the boundary ∂R consists of five sides, (so that the intermediate point (λ 3 , φ 2 ) is regarded as an additional vertex), and even in this five-sided cell R every side satisfies the above requirement. The length of a side e = λ 1 ≤ λ ≤ λ 2 , φ = const. equals (λ 2λ 1 ) cos φ, while the length of a side e

= φ 1 ≤ φ ≤ φ 2 , λ = const. is φ 2 -φ 1 . Consequently, the area A R of the cell R is A R = λ2 λ1 dλ φ2 φ1 cosφ dφ = (λ 2 -λ 1 )(sin φ 2 -sin φ 1 ).

Geometry-compatible discretization of the divergence operator

Given any rectangular domain R of the form (4.1), the approximate flux divergence is now derived as an approximation of the integral of the flux along the boundary ∂R, divided by its area, as follows:

∇ T • F(x, u) approx = I R A R , I R = ∂R F(x, u) • ν ds approx , (4.2)
where ds is the arc length along ∂R and ν is the outward-pointing unit normal to ∂R ⊂ S 2 . In the limit λ 2 , φ 2 → λ 1 , φ 1 the approximation (4.2) to the divergence term approaches the exact value (1.2). 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 I R defined in (4.2) and consider the cell shown in Figure 2, under the assumption that u = u(λ, φ, t) is smooth on R. We propose to approximate the flux integral along each edge of R in the following way. As in Section 2, let us decompose the flux into its (λ, φ) components:

F(x, u) = F λ (λ, φ, u)i λ + F φ (λ, φ, u)i φ .
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 e as λ e,m = (λ 1 + λ 2 )/2 and φ e,m = φ 1 (see Figure 2), and likewise for the edge e .

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 u as constant u = u e,m along the side e ∈ ∂R, the total approximate flux is given by e

F(x, u) • ν ds approx = -h(e 2 , u e,m ) -h(e 1 , u e,m ) , (4.3)
where e 1 , e 2 are, respectively, the initial and final endpoints of e (with respect to the sense of the integration).

Summing up over all edges we obtain:

Claim 4.1 (Discrete geometry-compatibility condition.) Consider the gradient flux vector constructed in Claim 2.2. Then, if u ≡ const., I R = 0, so that

∇ T • F(x, u) approx = 0,
and thus a discrete version of the divergence-free condition (2.2) holds.

Remark 4.2

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 geometrycompatible flux F(x, u), such a result can be obtained only if the dependence on x is integrated exactly along each side, a requirement that must be imposed on the scheme.

Godunov-type approach to the numerical flux

We continue to deal with the gradient flux given in Claim 2.2. We assume different (constant) values of u = u(λ, φ, t) in grid cells and evaluate the numerical flux values at each edge from the solution to a Riemann problem with data comprising these values u(λ, φ, t) in the cells on either side of that edge. At the midpoint (λ e,m , φ e,m ) of each side e we solve the Riemann problem in a direction perpendicular to e, and denote the resulting solution u e,m . The corresponding fluxes are then evaluated as F(λ e,m , φ e,m , u e,m ).

We can split Eq. ( 1.3) by invoking the explicit form of the divergence (1.2), getting

∂u ∂t + 1 cos φ ∂ ∂λ F λ (λ, φ, u) = 0 for the side e : λ = λ 2 , (4.4) λ ∂u ∂t - 1 cos φ ∂ ∂φ F φ (λ, φ, u) cos φ = 0 for the side e : φ = φ 2 , (4.4) φ
Consider two adjacent cells, as in Figure 4 or in Figure 5. By fixing φ = φ e,m (resp. λ = λ e,m ) in (4.4) λ (resp. (4.4) φ ) we can evaluate u = u e,m as a one-dimensional solution at λ = λ e,m (resp. φ = φ e,m ). 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 u(λ, φ,

λ 1 λ 2 λ 3 φ 1 φ 2 φ e ,m u = u L u = u R u = u e ,m

M e

t n ) = u L (resp. u(λ, φ, t n ) = u R ) in the cell λ 1 < λ< λ 2 , φ 1 < φ < φ 2 (resp. λ 2 < λ < λ 3 , φ 1 < φ < φ 2 ), as in
∂u ∂t + tan φ e ,m ∂ ∂λ f 1 (u) cos λ + f 2 (u) sin λ - ∂ ∂λ f 3 (u) = 0. (4.5) Setting g(λ, u) = tan φ e ,m f 1 (u) cos λ + f 2 (u) sin λ -f 3 (u), (4.6)
we see that equation (4.5) is the scalar one-dimensional conservation law

∂u ∂t + ∂ ∂λ g(λ, u) = 0, t≥ t n (4.7)
subject to the initial data u = u L (resp. u = u R ) for λ < λ 2 (resp. λ > λ 2 ). Likewise, we repeat the former analysis for φ-adjacent cells by taking the constant states u(λ, φ, t n )=u L , u(λ, φ, t n )=u R in cells λ 1 <λ<λ 2 , φ 1 <φ < φ 2 , λ 1 <λ<λ 2 , φ 2 <φ<φ 3 , as depicted in Figure 5. At the point M (λ = λ e,m , φ = φ 2 ), the equation (4.4) φ then takes the form

∂u ∂t + 1 cos φ ∂ ∂φ -sin λ e,m cos φf 1 (u) + cos λ e,m cos φf 2 (u) = 0. (4.8)
We then set the φ-flux function k(φ, u) = -sin λ e,m f 1 (u) + cos λ e,m f 2 (u) cos φ (4.9) so that equation (4.8) is the scalar one-dimensional conservation law

∂u ∂t + 1 cos φ ∂ ∂φ k(φ, u) = 0, t≥ t n (4.10)
subject to the initial data u = u L (resp. u = u R ) for φ < φ 2 (resp. φ > φ 2 ).

Solution to the Riemann problem

The solution at the discontinuity λ=λ 2 at the initial time t = t n is given by the Riemann solution to (4.4) λ . For simplicity of the presentation we specialize here to the flux (4.7). Since the dependence of g(λ, u) on λ is smooth, this solution is obtained by fixing λ = λ 2 , thus solving the classical conservation law

∂u ∂t + ∂ ∂λ g(λ 2 , u) = 0, t≥ t n (4.11)
subject to the initial jump discontinuity of u.

We denote this solution by u 2,m . Observe that the flux g(λ, u) 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 (λλ 2 )/(tt n ). The value u 2,m is the value along the line λ=λ 2 . It therefore corresponds either to a sonic wave, namely g (λ 2 , u 2,m )=0, or to an "upwind value" u=u L (resp. u=u R ) 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 g(λ 2 , u) is well-known and goes back to classical works by Oleinik and others. We recall it here briefly. Assume first that u L <u R . Consider the convex envelope of g, namely, the largest convex continuous function g c , over the interval [u L , u R ], such that g c ≤g at all points. Clearly, g c =g in "convex sections" of the graph of g, while it consists of linear segments when g c <g. It is easy to see that the "convex segments", where g=g c , represent rarefaction waves (in the full Riemann solution) while the linear segments represent jumps (i.e., shock waves). In particular, the solution u 2,m is given by the following formula:

u 2,m = v min , where g(λ 2 , v min ) ≤ g(λ 2 , v) for all v ∈ [u L , u R ].
(4.12)

There are in fact three possibilities for this solution: a) u L < u 2,m < u R , which implies that g (λ 2 , u 2,m ) = 0 (a sonic point). b) u 2,m = u L , the whole wave pattern moves to the right. c) u 2,m = u R , the whole wave pattern moves to the left.

Similarly, in the case u L >u R , we construct the "concave envelope" of g, namely, the smallest concave continuous function g c such that g c ≥g. Again the linear segments correspond to jump discontinuities while the concave segments (g=g c ) correspond to rarefaction waves. The solution to the Riemann problem is now given by u 2,m =v max , where

g(λ 2 , v max )≥g(λ 2 , v), v ∈ [u R , u L ].
As above, there are three possibilities for the solution (sonic, left-upwind, or right-upwind).

Replacing in the foregoing analysis the λ-flux function g(λ 2 , u) by the φ-flux function k(φ 2 , u), the equation (4.10) reads ∂u ∂t

+ ∂ ∂φ -sin λ 2,m f 1 (u) + cos λ 2,m f 2 (u) = 0, t≥ t n . (4.13)
We get the Riemann solution to (4.13) in the three cases a), b), c) as above.

Convergence proof

The computational elements ("grid cells") are denoted in [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF] by K. Their sides are denoted by e and the flux function across e is given by f e,K (u, v), where u is the (constant) value in K and v is the value in the neighboring cell (sharing the same side e) K e . In our grid of the sphere, some cells are actually pentagons; these are the cells whose lowerlatitude side (along a latitude φ = const) 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 [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF] 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 f e,R (u, v) be the numerical flux calculated on the side e of the computational cell R, using (4.3), where the midpoint value of u is obtained from the Riemann solution. Then f e,R (u, v) satisfies the assumptions (5.5)-(5.7) of [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF], 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 e : λ = λ 2 , which is given by F λ in the equation (4.4) λ . The procedure for integrating the flux across e is described by (4.3), while in Subsection 4.4 the calculation of F λ (λ 2 , φ 2,m , u 2,m ) is described. It can be summarized as follows.

First, the solution u 2,m to the Riemann problem associated with equation (4.4) λ is found, assuming u, v to be the values on the two sides. However, note that F λ depends explicitly on φ, and to be precise we need to replace in (4.4) λ the mean value φ 2,m by φ. Thus, we find u 2,m = u 2,m (φ).

Clearly, in the case u = v we get identically u 2,m (φ) = u = v and so the exact flux satisfies

F λ = F λ (λ 2 , φ, u 2,m )
and its integration will give exactly the approximate value f e,K (u, v) =h(e 2 , u 2,m )h(e 1 , u 2,m ) , as in (4.3). Thus, condition (5.5) in [START_REF] Amorim | Hyperbolic conservation laws on manifolds. Total variation estimates and the finite volume method[END_REF] is satisfied.

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 u < v, one considers the convex envelope of F λ , as defined in (4.4) λ (with φ = φ 2,m ) and then considers u 2,m as the minimal value on this envelope (over [u, v]). Clearly changing u upward will either change u 2,m upward or leave it unchanged. This completes the proof. 2

Second-order extension based on the GRP solver

To improve the order of accuracy, we consider again the cell λ 1 <λ<λ 2 , φ 1 <φ<φ 2 and assume that u is linearly distributed there. We use u L,λ , u L,φ (resp. u R,λ , u R,φ ) to denote the slopes in the cell to the left (resp. right) of the side λ=λ 2 . We also denote by u L (φ) (resp. u R (φ)) the limiting value (linearly distributed) of u at λ=λ 2 -(resp. λ=λ 2 +). Clearly, the solution to the Riemann problem across the discontinuity is a function of φ, and we denote it by u 2,m (φ), which conforms to our notation in Subsecion 4.4 above (where u was constant on either side of the discontinuity). The value of u 2,m (φ) is obtained by solving the Riemann problem associated with Eq. (4.4) λ with φ 2,m replaced by φ, subject to the initial data u L (φ), u R (φ). Restricting to the middle point φ = φ 2,m , the solution u 2,m (φ 2,m ) (at λ = λ 2,m ) is in one of the three categories listed above (i.e., sonic, left-upwind, right-upwind). By continuity, the solution u 2,m (φ) will still be in the same category for φφ 2,m sufficiently small. The solution at (λ 2,m , φ 2,m ) 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 4.3. We assume that the flux vector depends explicitly on x, as in (2.1). In what follows we use for simplicity the "imbedded" notation x = (x 1 , x 2 , x 3 ) 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)

Φ(x, u) = ∇ x h(x, u) = q 1 (x 1 )f 1 (u) i 1 + q 2 (x 2 )f 2 (u) i 2 + q 3 (x 3 )f 3 (u) i 3 , (5.1)
The zero-divergence identity is obtained as a result of expressing Φ as a gradient ∇h in the sense of Claim 2.2.

For our choice of Φ such a representation of Φ as gradient of h is obtained when h is taken as

h(x, u) = r 1 (x 1 )f 1 (u) + r 2 (x 2 )f 2 (u) + r 3 (x 3 )f 3 (u) , (5.2)
and q j (x j ) = r j (x j ), j = 1, 2, 3. Using (1.2) together with the geometry-compatibility property, we get an explicit form of the conservation law (1.3) in our case as

∂u ∂t -sin λq 1 (x 1 ) ∂ ∂φ f 1 (u) + cos λq 2 (x 2 ) ∂ ∂φ f 2 (u) + tan φ cos λq 1 (x 1 ) ∂ ∂λ f 1 (u) + sin λq 2 (x 2 ) ∂ ∂λ f 2 (u) -q 3 (x 3 ) ∂ ∂λ f 3 (u) = 0.
(5.

3)

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, (4.4), 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

∂u ∂t + tan φ 2,m q 1 (x 1 ) cos λ ∂ ∂λ f 1 (u) + q 2 (x 2 ) sin λ ∂ ∂λ f 2 (u) -q 3 (x 3 ) ∂ ∂λ f 3 (u) = 0. (5.4)
Note that the coefficients are retained as functions of λ and are not "frozen" at λ = λ 2,m . 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 λ → λ 2,m , then substitute λ = λ 2,m . The λ-edge flux function g(λ, u) (compare (4.6)), is now extended to g(x, u) as g(x, u) = tan φ 2,m q 1 (x 1 ) cos λf 1 (u) + q 2 (x 2 ) sin λf 2 (u)q 3 (x 3 )f 3 (u), and the scalar one-dimensional conservation law under consideration is now rewritten as an equation with a source term (a balance law)

∂u ∂t + ∂ ∂λ g(x, u) = S λ , t > t S λ = tan φ 2,m f 1 (u) ∂ ∂λ q 1 (x 1 ) cos λ + f 2 (u) ∂ ∂λ q 2 (x 2 ) sin λ -f 3 (u) ∂ ∂λ q 3 (x 3 ), (5.5)
subject to the initial data (for u and its slope)

u L (φ 2,m ), u L,λ (resp. u R (φ 2,m ), u R,λ ) for λ < λ 2 (resp. λ > λ 2 ).
Observe that the equation is written in a "quasi-conservative form", which offers more convenience in the GRP treatment [START_REF] Ben-Artzi | Generalized Riemann problems in computational fluid dynamics[END_REF]Chap. 5]. The right-hand side term S λ is just the result of the λ differentiation of the flux g(x, u). Obviously, the geometry-compatibility condition implies that this source term should cancel out with the corresponding source term in the "φ-split" equation. The solution u 2,m 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 S λ is added to terms arising from the piecewise-linear initial data, in producing the time-derivative of the solution u 2,m (φ 2,m ) + ∂u ∂t (λ 2,m , φ 2,m , t n +) ∆t 2 , ∆t = t n+1t n . As explained above, u 2,m (φ 2,m ), the solution to the associated Riemann problem, is obtained by using the "edge values" u L (φ), u R (φ). It remains, therefore, to determine the instantaneous time-derivative ∂u ∂t (λ 2,m , φ 2,m , t n +), as is outlined below.

The time-derivative of u is given by

∂u ∂t (λ 2,m , φ 2,m , t n +) = -u m,λ ∂ ∂u g(x, u)| λ 2,m ,φ 2,m ,u 2,m ,
where the slope value u m,λ is obtained by "upwinding", determined by the associated Riemann problem as follows (we start with the "easy" categories b), c) above).

b) u 2,m = u L (φ 2,m ). Then, the wave moves to the right and we set u m,λ = u L,λ .

c) u 2,m = u R (φ 2,m ). Then, the wave moves to the left and we set u m,λ = u R,λ .

Finally, the first category deals with the sonic case. As noted above, it remains sonic in the neighborhood of φ 2,m , so that we have there ∂ ∂u g(x, u)| λ 2,m ,φ 2,m ,u 2,m . The timederivative of u reduces therefore to ∂ ∂t u(λ 2 , φ 2,m , t=t n +) = 0.

Finally, the "φ-split" equation obtained from (5.3), is treated in analogy with the "λ-split" procedure outlined above. As discussed in Section 3 (see the discussion of solutions to (3.1)) it is clear that the solution here (as a function of λ) is identical to the periodic solution for the Burgers equation in R 1 , with periodic boundary conditions on [0 < x < 2π]. However, we compute the numerical solution here on our spherical grid, and we need to check not only that it conforms with the one-dimensional case but that it does not "leak" beyond the band supporting the initial data. The results at the shock formation time t s = 1/2π are shown in Figure 6 for ∆λ = 2π/16, in Figure 7 for ∆λ = 2π/32 and in Figure 8 for ∆λ = 2π/64. These GRP solutions to (4.7) clearly converge to the exact solution with refinement of the λ grid, and are comparable to the corresponding solution to the scalar conservation law in R 1 with ∆x = 2π/22. We refer to Corollary 3.3 and using the notation there we take the flux vector and initial data as:

f 1 (u) = u 2 /2,
f 2 (u) = f 3 (u) = 0, u(λ, φ, 0) = cos λ cos φ.

(6.2)

Using the terminology of Corollary 3.3 we see that the initial function is the "simplest" possible function, corresponding to g(x 1 ) = x 1 .

As is shown in Figure 9, the numerical solution remains nearly unchanged in time after being subjected to integration up to t = 5 by the GRP scheme with constant time step ∆t = 0.05, the color maps of u(λ, φ, t) at the initial and final times are virtually indistinguishable. The shown grid has latitude step ∆φ = π/60, and an equatorial longitude step ∆λ = π/128. A measure u dif f to the numerical solution error is defined

Fig. 1 .

 1 Fig. 1. Web grid on a sphere

λ λ 1 λ 2 ( 1 φ 2 (Fig. 2 .

 2122 Fig. 2. Rectangular cell R as part of grid on S 2

Fig. 3 .

 3 Fig. 3. Five-sided rectangular cell R (on southern hemisphere of S 2 )

Fig. 4 . 1 λ e,m λ 2 φ 1 φ 2 φ 3 uFig. 5 .

 4121235 Fig. 4. Two λ-adjacent cells with constant states uL, uR

Figure 4 .

 4 At the point M (λ e ,m , φ e ,m ) Eq. (4.4) λ takes the form

6. Numerical tests 6 . 1 .

 61 First test case: equatorial periodic solutionsHere, the conservation law takes the form (3.1) and the flux function and initial data are given byf 1 (u) = f 2 (u) = 0, f 3 (u) = -2π (u 2 /2), u(λ, φ, 0) = sin λ, 0 < λ < 2π, 0 < φ < π/12,

Fig. 6 .- 1 Fig. 7 .- 1 Fig. 8 . 6 . 2 .

 6171862 Fig.6. Exact, GRP/SCL and GRP/SPHERE (∆λ=2π/16) solutions to the IVP (6.1) at t = 1/2π

Acknowledgments

The authors were supported by a research grant of cooperation in mathematics, sponsored by the High Council for Scientific and Technological Cooperation between France and Israel, entitled: "Theoretical and numerical study of geophysical fluid dynamics in general geometry". This research was also partially supported by the A.N.R. (Agence Nationale de la Recherche) through the grant 06-2-134423 and by the Centre National de la Recherche Scientifique (CNRS).

Fig. 9. Steady-state initial data (and solution) to the IVP (6.2) at t = 5.

Color map range scaled to (u min , umax) = (-0.998, 0.998).

as the area-weighted difference |u(λ, φ, 5)u(λ, φ, 0)|, obtained by summation over all grid cells. In this case we obtained u dif f = 0.0093, which is small relative to the full range u maxu min = 2. Hence, the GRP scheme produces an approximation to the steadystate solution u(λ, φ, t) = u(λ, φ, 0) over S 2 . This test case demonstrates that the scheme computes correctly the time-evolution for the non-constant data (6.2), by calculating an approximately zero value for the flux divergence in computational cells.

Third test case: confined solutions

We take (as in Claim 2.2) Φ(x, u) = ∇h(x, u), where h(x, u) = ψ(x 1 )x 1 f 1 (u). The function ψ(x 1 ) is defined by

The flux vector is then given by

The solution is clearly confined to the sector

2 of the sphere. Its boundary is a circle which intersects the meridian λ = 0 at φ = π 4 .

20

The flux in the subdomain x 1 ≤ 0 is given by

so if we take the initial data as ψ(x 1 )u 0 (x 1 ), where u 0 is the steady state solution of the second test case (and also the same f 1 (u)), the solution remains steady in that part, namely, in x 1 ≤ 0. Clearly, it evolves in time in the region 0 ≤ x 1 ≤ √ 2

2 , but vanishes identically (for all time) if √ 2

2 ≤ x 1 . The confined IVP was integrated in time up to t = 5 by the GRP scheme, using the same grid and time step as in the second test case (Subsection 6.2). The solution is represented by the color map in Figure 10. Comparing it to the corresponding initial map (not shown here), it seems nearly unchanged. In fact, the initial-to-final difference measure obtained is u dif f = 0.0057, which indicates a nearly steady solution in the strip 0 < x 1 < 1/2 . This test case demonstrates that the scheme computes correctly the time-evolution for the non-constant "confined" data (6.3).