Strong bounded variation estimates for the multi-dimensional finite volume approximation of scalar conservation laws

We prove a uniform bounded variation estimate for finite volume approximations of the nonlinear scalar conservation law $\partial_t \alpha + \mathrm{div}(\boldsymbol{u}f(\alpha)) = 0$ in two dimensions with an initial data of bounded variation. We assume that the divergence of the velocity $\mathrm{div}(\boldsymbol{u})$ is of bounded variation instead of the classical assumption that $\mathrm{div}(\boldsymbol{u})$ is zero. A uniform bounded variation estimate provides compactness for finite volume approximations in $L^p$ spaces, which is essential to prove the existence of a solution for a partial differential equation (p.d.e.) with nonlinear terms in $\alpha$, when uniqueness of the solution is not available. The finite volume schemes analysed in this article are set on nonuniform Cartesian grids. A uniform bounded variation estimate for finite volume solutions of the conservation law $\partial_t \alpha + \mathrm{div}(\boldsymbol{F}(t,\boldsymbol{x},\alpha)) = 0$, where $\mathrm{div}_{\boldsymbol{x}}\boldsymbol{F} \not=0$ on nonuniform Cartesian grids is also proved. Results from numerical tests are provided, and they complement theoretical findings.


Introduction
Consider the following scalar hyperbolic conservation law in R 2 with a homogeneous source term and an initial data of bounded variation (BV ): ∂ t α + div(uf (α)) = 0 in Ω T and α(0, ·) = α 0 in Ω, (1.1) where α is the unknown, u = (u, v) is the advecting velocity, Ω T := (0, T ) × Ω, Ω := I × J, I := (a, b) ⊂ R and J := (c, d) ⊂ R are intervals. For technical simplicity assume that u = 0 on ∂Ω. The flux f is Lipschitz continuous with Lipschitz constant Lip(f ). Finite volume methods are extensively used to discretise and compute numerical solutions to (1.1) since such schemes respect the conservation property associated with the underlying p.d.e. Usually, conservation laws of the form (1.1) model density or concentration of conserved quantity in a coupled system, wherein the velocity u is governed by a p.d.e. that contains nonlinear terms in α. For instance, consider the Peaceman model [3] used to study oil recovery using a pressurised solvent in a porous medium. Here, the Darcy velocity is described by u = −K ∇p/µ(α), (1.2) where µ(α) is the viscosity, K(x) is the permeability tensor, and α follows the conservation law Here, a + := max(a, 0), a − := − min(a, 0) for a ∈ R, φ(x) is the porosity of the medium, and D(x, u) is the diffusion tensor. Also, consider the two-phase tumour growth problem [9], where the cell phase velocity u is modelled by −div µα(∇u + (∇u) T ) + λαdiv(u) where µ and λ are the viscosity coefficients, k is the traction coefficient, α * is a positive parameter that controls intra-cellular attraction, p is the pressure, I d is the 2 × 2 identity tensor, and α is modelled by (1.1) with a nonlinear source function in α. To show that a possible limit of discrete solutions obtained from a finite volume scheme applied to (1.1) or (1.3) satisfies (1.2) or (1.4), respectively and hence to prove the existence of a solution, we need that the discrete solutions converge to the limit in strong L p -norm, where p ≥ 1. Otherwise, pass to the limit arguments on the nonlinearities (in α) appearing in (1.2) and (1.4) can be challenging. A feasible way to obtain strong L p -norm convergence is to show that the discrete solutions have uniform BV and invoke Helly's selection theorem [10, p. 176] to extract a strongly converging subsequence. However, a uniform BV estimate on finite volume solutions on unstructured grids is yet an unsolved problem. Moreover, the velocity vector field u is not necessarily divergence-free of which (1.2) is a direct example. Hence, while attempting to obtain a uniform BV estimate on discrete solutions of (1.1), we need to account for divergent velocity vector fields also.
Total variation properties of weak and entropy solutions of (1.1) are rather classical results. E. Conway and J. Smoller [5] studied conservation laws of the form with BV initial data and (f j ) j=1,...,d are assumed to be in C 1 (R; R). They studied a finite difference scheme on a uniform Cartesian grid (see Definition 2.2) and showed that discrete solutions have uniform BV . The limit solution obtained from a strongly convergent subsequence is then showed to be a weak solution and is of BV . N. Kuznetsov [14] provided early results on BV properties of entropy solutions of (1.5). This paper [14] establishes that the BV seminorm of the entropy solution to (1.5) at any time t is bounded by the BV seminorm of the initial data. M. G. Crandall and A. Majda [6] considered monotone finite difference approximations of (1.5) with BV initial data on uniform Cartesian meshes and established uniform BV estimate for discrete solutions. This estimate is used to prove the convergence of the discrete solutions to the unique entropy solution in strong L 1 -norm and to prove that the entropy solution also inherits the BV property of the discrete solutions. Later, this work was extended to nonuniform Cartesian meshes by R. Sanders [17]. B. Merlet and J. Vovelle [15,16] considered linear advection equations of the form (1.1) with f (α) = α, u ∈ W 1,∞ (R + × R d ; R d ), and div(u(t, ·)) = 0. The BV seminorm of the unique weak solution of this problem constructed using the characteristic method is bounded and the bound depends on the BV seminorm of the initial data. However, discrete solutions corresponding to this problem obtained by using finite volume schemes on general polygonal meshes are not proved to satisfy a uniform BV estimate (see the Remark 1.5 in [16, p. 7]). In fact, to show that the finite volume solutions converge to the entropy solution, whose existence is known a priori, it is enough to have a weak BV estimate [2, p. 143] [11, p. 161] of the following form N n=0 δ e |f (α p e ) − f (α n e )| ˆe u(t n , ·) · n e ds ≤ Ch −1/2 , (1. 6) where δ is the temporal discretisation factor, h is the spatial discretisation factor, e is an edge of a polygon K in the mesh, n e is the outward unit normal to e with respect to K, α are the values of a discrete solution on the neighbouring polygons of e. The weak BV estimate ensures convergence in nonlinear weak- * sense (see Definition 6.3 in [11, p. 100]) to a Young measure, called a process solution. We can establish that the process solution is indeed a function by proving the uniqueness of the process solution. In this scenario, the nonlinear weak- * convergence actually becomes strong L p convergence (see Theorem 6.4 and 6.5 in [11, p. 187-188]). Uniqueness of the process solution is crucial in this technique and hence, it is difficult to use it in the case of coupled systems like (1.2) and (1.4). The relationship between process solution and function solution is not very clear in this case and an a priori compactness result like a uniform BV estimate is necessary to obtain strong L p convergence.
A recent uniform BV estimate on finite volume solutions of conservation laws of the form (1.5) on uniform Cartesian grid is obtained by K. H. Karlsen and J. D. Towers [13]. They consider (1.5) with an auxiliary boundary condition f · n Ω = 0, where n Ω is the outward unit normal to ∂Ω. Chainais-Hillairet [1] also provided a uniform BV estimate on finite volume solutions of fully nonlinear conservations laws on uniform square Cartesian grids (see subsection 3.1 for details).
In all of the works reviewed above, either the advecting velocity vector is component-wise constant (see (1.5)) or the advecting velocity is assumed to be divergence-free. However, this may not be the physical case as we have already shown by the examples of (1.2) and (1.4). While discretising physical models, it is imperative to refine the regions where discontinuities of the solution are expected and to retain other regions relatively coarse so that the scheme remains economical. A uniform BV estimate is crucial in enabling the nonlinear terms to converge and hence to prove existence of a solution. In this article, we analyse a class of finite volume schemes that serves these three requirements. The main contributions are (C.1) In the conservation law (1.1), the assumption that div(u) = 0 is relaxed.
(C.2) A finite volume scheme on nonuniform Cartesian grids is considered and the analysis is kept general by considering the class of monotone numerical fluxes. The nonuniformity of Cartesian grid can be used to refine the mesh adaptively and economically.
(C.3) Finite volume solutions satisfy a uniform BV estimate in space and time and this result is extended to the case of fully nonlinear conservation laws analysed by Chainais-Hillairet [1].
The main idea used to achieve (C.3) is to compute the variation of the discrete solution along x and y axes separately and it involves two major difficulties. Firstly, the term α div(u) serves as an additional source function, thanks to (C.1). The difference of α div(u) at time step t n+1 across neighbouring control volumes is estimated in terms of the BV seminorm of div(u) and L ∞ bound of α at time step t n . Secondly, while estimating the difference of the discrete solution across two control volumes in x direction, we obtain terms that contain differences of numerical fluxes across the y direction and vice-versa. This is a potential obstacle for the standard technique of writing the variation of the discrete solution at t n+1 across two control volumes as a convex linear combination of variations of the discrete solutions across neighbouring control volumes at t n . We introduce the idea of an intermediate nodal flux, which is the numerical flux across control volumes sharing only a single vertex, to transform the vertical differences into horizontal ones and vice-versa. This helps to obtain a relation of the form where BV (n) is the BV seminorm of the discrete solution at t n , and A(t) depends on BV seminorm of div(u) and ∇u(t, ·) L ∞ (Ω) . Finally, an application of induction on (1.7) yields the BV estimate on the discrete solution. This article is organised in the following fashion. In Section 2, we present the necessary notations, assumptions, function spaces, and the finite volume scheme. The main theorem is also presented in Section 2. The uniform strong BV estimates are presented in Section 3 and the results are extended to the fully nonlinear case in subsection 3.1. The numerical results and discussion are presented in Section 4. The conclusions are provided in Section 5.

Main theorem
In this section, we present the finite volume scheme and state the main theorem that establishes the BV property of the discrete solution.
Also, define the following norms for a function v : For a function β : (a, b) → R, define the total variation by where P := {a = x 0 , . . . , x I+1 = b} is a partition of (a, b). It is a classical result that |β| BVx(a,b) = T.V.(β) (see Appendix A of [12]). Next, we define an admissible grid.

Presentation of the numerical scheme
Define the spatial discretisation factor h max by h max = max i,j {k i , h j }, which quantifies the size of the Cartesian grid X k × Y h . Let T δ defined by T δ := {t 0 , . . . , T N } be a discretisation of (0, T ), where t 0 = 0 and t N = T . Define the spatial discretisation factor by δ = max n δ n , where δ n = t n+1 − t n . For technical simplicity a uniform temporal discretisation is taken, wherein δ n = δ ∀ n. However, note that the results in this article hold with a nonuniform temporal discretisation also. Integrate (1.1) on the time-space control volume (t n+1 , t n )×K ij , where K ij = (x i−1/2 , x i+1/2 )× (y j−1/2 , y j+1/2 ) and apply the divergence theorem to obtain where n ij is the outward unit normal to ∂K ij and u = (u, v). Term I 2 in (2.1) is approximated by the numerical flux g : R 2 → R. Replace I 1 by the difference formula k i h j (α n+1 i,j − α n i,j ) and Locations of the discrete unknowns α n i,j , velocities u i−1/2,j and v i,j−1/2 in an admissible grid is shown in Figure 1. The discrete scheme is given by where µ i = δ/k i and λ j = δ/h j . We set the discrete initial data as follows The terms F i,j and G i,j can be expressed as where (2.3c) Observe that D i,j (α n i−1,j , α n i,j ), D i,j (α n i,j , α n i−1,j ), D i,j (α n i,j , α n i,j−1 ) and D i,j (α n i,j−1 , α n i,j ) (hence, M x i−1/2,j and M y i,j−1/2 ) are nonnegative due to the monotonicity of g and to g(α n i,j , α n i,j ) = f (α n i,j ). Use (2.3a) and (2.3b) to transform the right hand side of (2.2a) into a convex linear combination of the terms α n l,m , where (l, m) ∈ {(i, j), (i − 1, j), (i + 1, j), (i, j + 1), (i, j − 1)}, and this yields The Lipschitz continuity of the function g yields . . , I and j = 0, . . . , J. We show that under a proper Courant-Friedrichs-Lewy (CFL) condition the iterates α n i,j can be bounded independently of the choice of discretisation factors. The piecewise time-reconstruct is defined next.
The assumptions used in this article are presented here.
(AS.1) The flux f and the numerical flux g are Lipschitz continuous with Lipschitz constants Lip(f ) and Lip(g), respectively. (AS. 2) The numerical flux g is monotonically nondecreasing in the first variable and nonincreasing in the second variable, and satisfies g(a, a) = f (a) ∀a ∈ R.
(AS.3) There exists a generic constant C ≥ 0 such that The main result in this article is as follows.
Theorem 2.4 (bounded variation). Let X k ×Y h be an admissible grid, and assumptions (AS.1)-(AS.3) and the Courant-Friedrichs-Lewy (CFL) condition  [1,11], which yields (a) on compact subsets of R d × R + . Though (b) apparently seems to be restrictive, it is pivotal in bounding the difference of div(u) between two control volumes (see (3.18)). Indeed, we can relax this assumption to div(u) ∈ L 1 t L ∞ (Ω T ), which is the formally correct choice and is used in the seminal work [8] by DiPerna and Lions. However, with this less restrictive assumption, Proposition II.1 in DiPerna and Lions [8] only guarantees the existence of a weak solution α ∈ L ∞ t L 1 (Ω T ). Therefore, (b) is justified for the reason that a stronger convergence of the finite volume solutions is obtained and the limiting solution also has the higher BV -regularity.

Variation in spatial and temporal coordinates
We let the hypotheses of Theorem 2.4 to hold throughout the sequel of this article and recall that α h,δ is the time-reconstruct in the sense of Definition 2.3. The proof of Theorem 2.4 is accomplished through three steps: establish the • boundeness of α h,δ in Proposition 3.1, • spatial BV estimate of α h,δ in Proposition 3.2, and The proof of Proposition 3.1 is based on writing α n+1 i,j as convex linear combination of values of α h,δ at the previous time step.
where B u := exp C ∇u L 1 t L ∞ (Ω T ) and C has the same dependencies as in Theorem 2.4. The proof of Proposition 3.2 is achieved in five intermediate steps, which are as follows Step 2 Use the intermediate nodal fluxes (see Figure 4) to transform the vertical differences in J i,j into horizontal differences.
Step 4 Estimate variations of ∂ x u and ∂ y v in terms of the BV seminorm of div(u).

Proof.
Step 1: Consider the difference between the scheme (2.2a) written for α n+1 i,j and α n+1 The term H i,j gathers the variation in the x-direction and use (2.3a) to write Step 2: The goal of this step is to transform the horizontal difference of variations between the vertical levels (i − r, j + s) and Figure 2: Spatial locations of the numerical fluxes g + i,j+1/2 and g − i,j+1/2 .
in J ij of (3.4) so that the resulting terms can be combined to form a convex linear combination of differences of α h,δ (t n , ·) between neighbouring rectangles. Split J i,j into two components . The numerical fluxes involved in J + i,j and J − i,j can be assigned with spatial locations as in Figures 2(a) and 2(b). Observe that J i,j /λ j is the horizontal variation between differences (of a l,j+1/2 and a l,j−1/2 , l ∈ {i, i − 1}) across two vertical levels as in Figure 3(a). However, this form does not yield any terms like α n i,r − α n p,r , where p ∈ {i + 1, i − 1} and r ∈ {j + 1, j, j − 1}, and thereby annihilates any chance of expressing α n+1 i,j − α n+1 i−1,j as a linear combination of such terms, which is crucial in controlling the growth of spatial variation over time. This problem can be fixed by considering the terms J + i,j and J − i,j as vertical variations between differences across two horizontal levels as in Figure 3 Figure 3: Differences between horizontal and vertical levels.
The term T j+1/2 1 in (3.7) can be combined with T j−1/2 1 coming from the T + 2 of (3.6) to obtain a bound using the total variation of div(u). The difference between the numerical fluxes on the edges in T j+1/2 2 needs to be carefully analysed. We cannot transform the difference g(α n i,j , α n i,j+1 ) − g(α n i−1,j , α n i−1,j+1 ) into a difference quotient as we did in (2.3c), since g(α n i,j , α n i,j+1 ) and g(α n i−1,j , α n i−1,j+1 ) have different first and second arguments. As a result, the Lipschitz continuity of g cannot be employed to bound any resulting difference quotients. To overcome this problem, we introduce an artificial nodal flux g(α n i−1,j , α n i,j+1 ) arising from two diagonally opposite control volumes as in Figure 4. The nodal flux helps to split the term J i,j as follows: where the difference quotients E m : R 3 → R are defined by (3.8b) Note that the sums (1 ± (±)) used in (3.8b) are understood as (1 ± (±1)). Use the identity a + = a + a − to transform the differences (v n + i,j+1/2 − v n + i−1,j+1/2 ) and (v n + i,j−1/2 − v n + i−1,j−1/2 ) in J + i,j and combine the resulting negative parts with the corresponding negative parts in J − i,j . This yields Step 3: Combine (3.4), (3.5a), (3.8a), and (3.9) to obtain Note that in (3.10b) the terms E − are nonpositive and E + are nonnegative. This fact along with the CFL condition ensures that 1 − c i,j is nonnegative. Take absolute value on both sides of (3.10a), multiply by h j , sum on i = 1, . . . , I and j = 0, . . . , J, and use the condition that u = 0 on ∂Ω to change the indices appropriately to obtain . (3.11) The term 1 − c i,j and coefficients of |α n i,j − α n i−1,j | in the second to the sixth sum on the right hand side of (3.11) adds up to one, and this yields Use the Lipschitz continuity of x → x − (with constant 1) and g, Lipschitz continuity of v in the x-direction, and grid regularity condition of Definition 2.2 to obtain Step 4: Apply sum-product identity (A.1) on K g i,j (see (3.9)) to obtain Write the term K f i,j (see (3.5b)) as Use the Lipschitz continuity of g, Lipschitz continuity of v in the x-direction, and the grid regularity condition of Definition 2.2 to bound K g,1 i,j . This yields (3.14) A use of sum-product identity (A.1) on K f,1 i,j yields Therefore, |K f,1 i,j | can be bounded by The sum K g,2 i,j + K f,2 i,j can be written as The Lipschitz continuity of g and f and the fact that g(a, a) = f (a) for a ∈ R yield | − 2f (α n i,j ) + g + i,j+1/2 + g + i,j−1/2 | ≤ Lip(g)|α n i,j − α n i,j−1 | + Lip(g)|α n i,j − α n i,j+1 | and (3.17a) |2f (α n i−1,j ) + g + i,j+1/2 + g + i,j−1/2 | ≤ 2Lip(f )|α n i,j − α n i−1,j | + Lip(g)|α n i,j − α n i,j−1 | + Lip(g)|α n i,j − α n i,j+1 |.

Proposition 3.3 (temporal variation). The function α h,δ satisfies
Take absolute value on both sides of (3.22), multiply by h j k i , sum over n = 0, . . . , N , i = 0, . . . , I and j = 0, . . . , J, and use the homogeneous boundary condition on u to obtain (AS.4) S ∈ L 1 t L ∞ (Ω T ) and S(t, x, z) is Lipschitz continuous with respect to z (with constant Lip z (S)), uniformly with respect to t and x, and is Lipschitz continuous with respect to x (with constant Lip x (S)), uniformly with respect to t and z.
In this case, we obtain the following corollary to Theorem 2.4.
where the constant C additionally depends on Lip x (S), Lip z (S), and |S| L 1 t BVx,y . Proof. It is enough to estimate variation of the source term in the x direction, which can be written as Add and subtract´t n+1 tn ffl K i,j S(t, x, α n i−1,j ) dt dx to (3.25) and group the terms appropriately to obtain Use the Lipschitz continuity of S with respect to the third argument to bound V 1 by δ Lip z (S)|α n i,j − α n i−1,j |. Sum (3.26) for i = 1, . . . , I to obtain Rest of the proof follows by adding the terms in the right hand side of (3.27) to the right hand side of (3.19) and by following the steps from there on.

BV estimate with fully nonlinear flux
Theorem 2.4 can be extended to the case with fully nonlinear flux such as ∂ t α + div(F (t, x, α)) = 0 in Ω T and α(0, ·) = α 0 in Ω. (3.28) The strong BV estimate on finite volume schemes for (3.28) on square Cartesian grids is obtained by Chainais-Hilairet [1] under the assumption that div x (F ) = 0. This is equivalent to div(u) = 0 in the case, where F (t, x, α) = uf (α). We can relax this condition to obtain the following set of assumptions: and is Lipschitz continuous with respect to z (with constant Lip(F )), uniformly with respect to (t, x), and ∂ z F is Lipschitz continuous with respect to x (with constant Lip(∂ s F )), uniformly with respect to t and z, (AS.6) |div x (F )| L 1 t BVx,y < ∞ and div x (F ) is Lipschitz continuous with respect to z (with constant constant Lip(div x (F ))), uniformly with respect to t and x.
Use (AS.5) to write the flux F as F := (F 1 , F 2 ), F 1 = a + b, and F 2 = c + d, where a and c are monotonically nondecreasing and b and d are monotonically nonincreasing in z, uniformly with respect to t and x. In this case, we can set the following finite volume scheme on an admissible grid X h × Y k : with the initial condition 2.2b, where the numerical fluxes are defined, for γ ∈ {a, b}, and ∈ {c, d}, by x i−1/2 (t, x, y j+1/2 , s)dx dt.
Theorem 3.4 (bounded variation for fully nonlinear flux). Let the assumptions (AS.5)-(AS.6) and the following CFL condition hold: Then the piecewise time-reconstruct α h,δ : Ω T → R re-constructed from the values α n i,j obtained from the scheme (3.29) satisfies where C depends on T , α 0 , |div x (F )| L 1 t BVx,y , and Lip(div x (F )). The proof of Theorem 3.4 is based on two main ideas. Firstly, the terms in the scheme (3.29) are re-arranged and grouped appropriately so that the term´t n tn ffl K i,j div x (F )(t, x, α n i,j ) dxdt can be separately estimated (see (3.32)). Secondly, we employ the Lipschitz continuity of div x (F ) to bound difference of the terms {´t n+1 tn ffl K l,j div x (F )(t, x, α n+1 l,j ) dxdt : l = i, i + 1} by the BV seminorms´t n+1 tn |α h,δ (t, ·)| BVx dt and´t n+1 tn |div x (F )(t, ·, ·)| BVx dt.
Proof. Note that the scheme (3.29) can be expressed as It is enough to estimate |α h,δ | L 1 y BVx as we did in the proof of Proposition 3.2. Take the difference between the scheme (3.32) written for α n+1 i+1,j and α n+1 i,j . The difference T 1,i+1 − T 1,i can be estimated exactly as in the proof of Lemma 8 of [1]. Consider the difference |T 2,i+1 − T 2,i |: The term Q 2 can be estimated as Estimates from the proof of Lemma 8 in [1], (3.33) and (3.34) yield Apply induction on 3.1 and use similar arguments as in the proof of Proposition 3.3 to obtain (3.31).

Numerical examples
We consider three examples to demonstrate the conclusions of Theorem 2.4 and Theorem 3.4. In Example 4.1, we manufacture a source term such that the conservation law (4.1) has a smooth solution. In Example 4.2, the source term is set to be zero and a discontinuous function is chosen as the initial data, and as a result the exact solution also becomes discontinuous. Example 4.2 helps to understand how the discontinuities in the solution, which is a common phenomenon observed in hyperbolic conservation laws, affect the growth of BV seminorm. In   (x, y) ∈ Ω, and an appropriate source term S such that the problem ∂ t α + div(uf (α)) = S in Ω 1 and α(0, x, y) = α 0 (x, y) ∀ (x, y) ∈ Ω, (4.1) has the unique smooth solution α(t, x, y) = exp(t(x + y)) ∀ (t, x, y) ∈ Ω 1 , where Ω 1 = (0, 1) × Ω.
Note that in the case of nonlinear flux, the vector u is zero on the boundary of the square (−3, 3) 2 , and as a result we can take the boundary data (αu) |∂Ω · n |∂Ω = 0, where n |∂Ω is the outward normal to ∂Ω. This homogeneous boundary condition on u is useful since the exact solution to the problem (4.1) with a nonlinear flux is not available. The source term and the initial condition remain the same as in the case of linear flux.
The source term S N is chosen such that (4.2) has the smooth solution α(t, x, y) = exp(t(x + y)). Note that the divergence of the flux div (sin(( We consider two fluxes in the tests: (i) linear flux, f (s) = s and (ii) sinusoidal flux, f (s) = sin(2πs). The numerical flux used is Godunov, which is described by The BV and L 1 rates are defined by BV rate = log α h k+1 ,δ k+1 (T, ·) BVx,y / |α h k ,δ k (T, ·)| BVx,y log(h k+1 /h k ) and Discretisation factors, CFL ratio, BV norms, and BV rates corresponding to Cartesian and perturbed Cartesian grids are presented in Tables 2-5 and Tables 7-10. The L 1 errors and L 1 rates are also included whenever an exact solution is available. Arrangement of the contents in Tables 2-10 are outlined in Table 1 for clarity. The BV rates corresponding to hexagonal, triangular, and staggered are presented in Table 6 and Table 11. The quantities L 1 error and L 1 rate are omitted for these three families of meshes since they follow a trend exactly similar to that of pertubred Cartesian grids. The BV norms of these meshes are also not included considering the space constraints. The L 1 and BV rates of the discrete solutions obtained by applying scheme 3.29 to (4.2) of Example 4.3 is provided in Table 12.  Table 2  Table 7 linear Cartesian Table 3  Table 8 sinusoidal Cartesian Table 4  Table 9 linear perturbed Cartesian Table 5  Table 10 sinusoidal perturbed Cartesian Table 1: Arrangement of contents in Tables 2-5 and Tables 7-10.
We recall three classical results from the theory of convergence analysis of finite volume schemes for conservation laws of the type (1.1).
(R.1) For a BV initial data, finite volume approximations of conservation laws of the type (1.1) on structured Cartesian meshes converge with h 1/2 rate with respect to L ∞ t L 1 norm [14], and this result is extended to nearly Cartesian meshes by B. Cockburn et al. [4]. For generic meshes the L ∞ t L 1 convergence rate is h 1/4 [11, p. 188]. (R. 2) The BV seminorm of the finite volume solution grows with a rate not greater than h −1/2 [11, p. 168]. Further details can be found in [4, p. 1777] and the references therein.

Observations
In Tables 2-9, it can be observed that the order of convergence with respect to the L 1 -norm is well above 1/4. The BV seminorm grows as h decreases indicated by the negative values of BV rate in Tables 2-9. But the growth rate is well below h −1/2 except in the case of initial coarse meshes. The trend in L 1 rate is related to the trend in BV rate. A reduced L 1 rate is attributed to the fact that finite volume solutions on generic grids lack a uniform strong BV estimate. The weak BV estimate (1.6) diminishes the L 1 rate from h 1/2 to h 1/4 in the case of non-Cartesian meshes.

Observations for Example 4.1
When the flux in linear, mesh is Cartesian (uniform or nonuniform), and (1.1) possesses a smooth solution, we obtain first order L 1 rate and the BV rate decreases in magnitude but with oscillations. In the case of sinusoidal flux, the L 1 rate shows a slight reduction for coarse meshes but readily becomes well above h 1/2 , which is the theoretical L 1 rate. Here also, BV rate decreases in magnitude but with oscillations as h decreases.
The numerical tests with perturbed Cartesian meshes also shows a similar behaviour. The linear flux exhibits first order L 1 rate and the sinusoidal flux a slightly reduced L 1 rate but well above h 1/2 . However, the BV rate shows a steady reduction in magnitude in both the linear and sinusoidal case. The BV rate for hexagonal, triangular, and staggered meshes also show a steady decrease in magnitude as provided in Table 11.

Observations for Example 4.2
In Example 4.2, we see a prominent reduction in the L 1 rate and this is due to the discontinuities in the weak solution to (1.1). The explicit finite volume scheme introduce considerable numerical diffusion in the discrete solution by smearing out the sharp fronts, and thereby reducing the convergence rate. This reduction in the L 1 rate is visible in both the Cartesian and perturbed Cartesian cases (see Tables 7 and 9). The BV rate is also decreasing in magnitude similar to Example 4.1. In the sinusoidal flux case also BV rates show the same pattern (see Tables 8  and 10). For other non-Cartesian meshes also BV rate seems to be decreasing in magnitude as presented in Table 6  It is clear from Table 12 that the BV rate is decreasing in magnitude steadily as h decreases. This justifies estimate (3.31) in Theorem 3.4. The L 1 rate is also greater than the theoretical rate of h 1/2 except for the initial coarse mesh (see result (R.3)). Table 12 also complements Lemma 8 and Theorem 4 in [1], which provides the boundedness of the BV seminorm of discrete solutions corresponding to uniform square Cartesian grids.
Remark 4.4. Choice of the functions a, b, c, and d for the scheme (3.29) is not arbitrary. It is crucial that a and c and nondecreasing, b and d are nonincreasing, and the CFL condition (3.30) holds. We use the following pairs to obtain the results provided in Table (

Final remarks
In the case of Cartesian mesh, note that the BV rate decreases in magnitude as h decreases and the BV seminorm stabilises eventually, which agrees with the conclusion of Theorem 2.4. This is also supported by the higher values of L 1 rate than the theoretically predicted ones and the fact that the reduced convergence rate stems from lack of a strong BV estimate (see result (R.1)).

Conclusions
A uniform estimate on total variation of discrete solutions obtained by applying finite volume schemes on conservation laws of the form (1.1) on nonuniform Cartesian grids is proved. We relaxed the standard assumption that the advecting velocity vector is divergence free. This enables us to apply the finite volume scheme to problems in which the advecting velocity vector is a nonlinear function of the conserved variable. Since the underlying meshes are nonuniform Cartesian it is possible to adaptively refine the mesh on regions where the solution is expected to have sharp fronts. A uniform BV estimate is also obtained for finite volume approximations of conservation laws of the type (3.28) that has a fully nonlinear flux on nonuniform Cartesian grids. Numerical experiments support the theoretical findings. The counterexample by B. Després and numerical evidence from Table 13 indicate that nonuniform Cartesian grids are the current limit on which we can get uniform BV estimates and extending Theorem 2.4 to perturbed Cartesian grids ( Figure 5(e)) might be the immediate future step.

Acknowledgement
• The author is grateful to A/Prof. Jérôme Droniou, Monash University, Australia, Prof. Neela Nataraj, IIT Bombay, India and A/Prof. Jennifer Anne Flegg, Univeristy of Melbourne, Australia for their supervision, valuable comments, and suggestions. The author also thanks Prof. Claire Chainais-Hilairet, Prof. Thierry Gallouët, and Prof. Julien Vovelle for their comments and suggestions.
• Additionally, the author expresses gratitude towards ANZIAM Student Support Scheme, IIT Bombay, IITB-Monash Research Academy, Monash University, A/Prof. Jennifer Anne Flegg, and A/Prof. Jérôme Droniou for funding the travel and hospitality expenses during the author's stay at Monash University on February-March, 2020 during which this work was carried out.