D Sarmany
email: d.sarmany@leeds.ac.uk

M Hubbard
email: m.e.hubbard@leeds.ac.uk

M Ricchiuto
email: mario.ricchiuto@inria.fr

Unconditionally stable space-time discontinuous residual distribution for shallow-water flows

Keywords: hyperbolic conservation laws, shallow-water equations, space-time discontinuous representation, residual distribution hyperbolic conservation laws, shallow-water equations, space-time discontinuous representation, residual distribution

This article describes a discontinuous implementation of residual distribution for shallow-water flows. The emphasis is put on the space-time implementation of residual distribution for the time-dependent system of equations with discontinuity in time only. This lifts the timestep restriction that even implicit continuous residual distribution schemes invariably suffer from, and thus leads to an unconditionally stable discretisation. The distributions are the space-time variants of the upwind distributions for the steady-state system of equations and are designed to satisfy the most important properties of the original mathematical equations: positivity, linearity preservation, conservation and hydrostatic balance. The purpose of the several numerical examples presented in this article is twofold. First, to show that the discontinuous numerical discretisation does indeed exhibit all the desired properties when applied to the shallow-water equations. Second, to investigate how much the time step can be increased without adversely affecting the accuracy of the scheme and whether this translates into gains in computational efficiency. Comparison to other existing residual distribution schemes is also provided to demonstrate the improved performance of the scheme.

Introduction

The framework of residual distribution (RD) has a nearly thirty-year-old history, having first been introduced in [START_REF] Roe | Fluctuations and signals[END_REF] as an alternative to finite volume schemes for the numerical discretisation of hyperbolic conservation laws. They can more naturally represent the most important underlying physical properties of steady-flow problems -especially avoiding spurious oscillations around discontinuities -while still providing second-order accurate approximations. For time-dependent problems, however, it initially proved to be rather more difficult to construct an efficient, second-order accurate, RD scheme that is also free of spurious oscillations. Much of the recent research in the field, therefore, has been aimed at devising truly time-dependent discretisations that retain all the desired properties of the underlying physical equations. Recent reviews of the field can be found in [START_REF] Abgrall | Residual distribution schemes: current status and future trends[END_REF][START_REF] Karman | Institute of Fluid Dynamics[END_REF][START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF].

Almost all of these schemes, however, use continuous discretisations of the physical variables. This leads to numerical methods that are both implicit, because of the presence of a global mass matrix, and require a time-step restriction to obtain stability, positivity or both. For many practical applications in physics or engineering, these methods can prove to be computationally inefficient and uncompetitive with more established discretisation methods, such as finite volumes [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Eleuterio | Riemann solvers and numerical methods for fluid dynamics[END_REF], finite elements [START_REF] Zienkiewicz | The finite element method[END_REF] or finite differences [START_REF] Bertil Gustafsson | Time dependent problems and difference methods[END_REF]. To remedy this shortcoming, two different approaches have been proposed in the literature for RD schemes that improve the efficiency of the time-step algorithm as well as keep the discrete representation continuous. One is to make the discretisation genuinely explicit and thus remove the need to solve an algebraic system at each time step [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: second order case[END_REF]. The other is to introduce double layers in a space-time implicit scheme, which removes the restrictions on the time step in the second layer [START_REF] Csík | Space-time residual distribution schemes for hyperbolic conservation laws on unstructured linear finite elements[END_REF].

In contrast, it is also possible to allow the discrete representation to be discontinuous across cell/prism interfaces. It is only relatively recently that this concept has been introduced -at first for steady-state equations only [START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF][START_REF] Hubbard | A framework for discontinuous fluctuation distribution[END_REF]. The potential benefits of the discontinuous RD (DRD) framework in space include: a) the possibility to implement h-adaptivity more easily; b) a mass-matrix structure that may render the inversion of the mass matrix computationally less expensive. As a further development, the work in [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF] took the first steps towards unconditionally stable space-time RD schemes by introducing discontinuous representation in time. One can also view this construction as a special, and slightly simpler, case of the space-time double-layer scheme [START_REF] Csík | Space-time residual distribution schemes for hyperbolic conservation laws on unstructured linear finite elements[END_REF][START_REF] Ricchiuto | Residual distribution for general timedependent conservation laws[END_REF], where the first layer is of zero height, while there is still no formal condition on the height of the second layer.

In this article, we extend the concept of space-time DRD to the shallow-water equations with possibly non-flat bottom topography. Mathematically, the non-flat bottom topography manifests itself as a source term in the governing equations. This means that apart from the usual requirements -positivity, linearity preservation, conservation and upwinding -a successful numerical scheme should also preserve the hydrostatic balance. Various RD schemes have been developed for shallow-water flows in the past one and a half decades. Both wave decomposition schemes [START_REF] Garcia-Navarro | Genuinely multidimensional upwinding for the 2D shallow water equations[END_REF][START_REF] Hubbard | Conservative multidimensional upwinding for the steady twodimensional shallow water equations[END_REF] and matrix distribution schemes [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] proved to be successful in representing the underlying physical properties. However, they are often computationally rather expensive owing partly to a time-step restriction mentioned above, and partly to poor iterative convergence. As a result and more recently, the explicit RD approach has been applied to the shallow-water system [START_REF] Ricchiuto | Explicit residual discretizations for shallow water flows[END_REF].

We take a different approach in this work and instead of replacing an intrinsically implicit scheme with an explicit one, we aim to eliminate the time-step restriction and construct an unconditionally stable space-time RD scheme. The focus of this article is therefore on discontinuous representation in time, which in this way promises to improve the efficiency of the implicit space-time discretisation. Our spacetime DRD discretisation is also designed to be upwind, conservative and hydrostatically well-balanced. Through a number of numerical experiments, we investigate whether these properties are satisfied, as well as positivity and linearity preservation -two other properties that are formally lost but experimentally often observed. Crucially, we also address the question of how far the time step can be increased before the discretisation adversely affects accuracy and positivity. Comparison with other implicit RD schemes is used to demonstrate the improved performance of this approach.

The remaining part of the article is organised as follows. In Section 2, we describe the framework of the space-time DRD for general hyperbolic conservation laws in the case of both scalar and system variables. Section 3 applies the scheme to the shallow-water equations, while Section 4 outlines details about the numerical implementation of the discretisation. An extensive study of the scheme via numerous computational examples is carried out in Section 5. Section 6 compares the space-time DRD scheme developed in this work with other existing RD schemes. Finally, we draw conclusions and provide outlook in Section 7.

Space-time discontinuous residual distribution

In this section, we give a general description for space-time DRD schemes for hyperbolic conservation laws with homogeneous right-hand side. Assume a two-dimensional spatial domain Ω ⊂ R 2 , its triangular tessellation Ω h , and the corresponding prismatic space-time mesh Ω t h . Let, furthermore, E denote a given triangle, E t denote the space-time prism defined over this triangle, D i = ∪ i∈E E, dΩ h = dx dy, dΩ t h = dx dy dt, d(∂Ω) = ds(x, y) and d(∂Ω t ) = ds t (x, y, t). In the last two definitions, s and s t denote the mappings s : Ω → ∂Ω and Ω t → ∂Ω t , respectively.

Scalar equations

Consider the scalar conservation law

∂ t u + ∇ • f = 0 or ∂ t u + a(u) • ∇u = 0 (1)
with appropriate initial conditions and Dirichlet boundary conditions at the inflow part of the domain.

Here f represents the conservative flux vector and a(u) = ∂f /∂u is the wave speed. The associated residual over each space-time prism E t is given by

φ Et = t n+1 t n E (∂ t u + ∇ • f ) dx dy dt. (2)
For the original (non-discretised) equation ( 1), it therefore holds that

φ = t n+1 t n Ω h (∂ t u + ∇ • f ) dx dy dt = E∈Ω h φ Et ,
provided that the solution u is bounded and piecewise differentiable. We carry out the RD numerical discretisation in the spirit of [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] and it consists of the following general steps.

1. In every space-time prism, replace the unknown u with an approximation u h that is linear in space and linear in time.

Transform the prism residual into a space-time boundary integral

φ Et = t n+1 t n E (∂ t u h + ∇ • f ) dx dy dt = Et ∇ t • f t dΩ t = ∂Et f t • n t d(∂Ω t ), (3)
where ∇ t = (∂t, ∂x, ∂y), f t = (u, f ), and n t is the outward-pointing unit vector normal to the space-time prism. 3. Let L denote the space of linear one-dimensional Lagrangian basis functions on the interval t n , t n+1 and W denote the space of linear two-dimensional basis functions associated with a triangle. Assume furthermore that u h ∈ L × W . Then the mid-point rule in time results in

φ Et = ∂Et f t • n t d(∂Ω t ) ≈ E u n+1 h -u n h dΩ + ∆t 2   ∂E f (u n ) • n d(∂Ω) + ∂E f (u n+1 ) • n d(∂Ω)   = E u n+1 h -u n h dΩ + ∆t 2 (φ n E + φ n+1 E ), (4)
where ∆t = t n , t n+1 is the 'height' of the prism. This is a second-order accurate approximation of the prism residual, as long as the spatial residuals, φ n E and φ n+1 E , are computed to be at least one order more accurate than the discretisation itself. 4. Distribute the prism fluctuation φ Et (4) to the six vertices of the prism in a conservative manner.

That is, the fractions of the residual sent to vertex i at time levels n and n + 1 are defined as

φ E i,n = β E i,n φ Et and φ E i,n+1 = β E i,n+1 φ Et , (5)
with i∈E β E i,n + i∈E β E i,n+1 = 1.

Solve the algebraic system

E∈Di

φ E i,n = 0, E∈Di φ E i,n+1 = 0, (6)
∀i ∈ Ω h at each time step. How to precisely distribute the residual does greatly depend on the direction of the flow, which is typically indicated by 'inflow' parameters. For the space-time scheme, these are defined as

k i,n = - ∆t 4 a • n i - |E| 3 , k + i,n = max(0, k i,n ), k - i,n = min(0, k i,n ), k i,n+1 = - ∆t 4 a • n i + |E| 3 , k + i,n+1 = max(0, k i,n+1 ), k - i,n+1 = min(0, k i,n+1
),

N t = 1 i∈E k + i,n + i∈E k + i,n+1 , (7)
where n i is the outward pointing normal vector opposite node i with length that of the edge opposite node i. The quantity a is a prism-averaged state of the values of a at the vertices of the prism.

If the discrete representation of the solution u h is continuous over the whole space-time domain, the scheme is only consistent if there is no residual sent back to the previous time level t n . This condition results in a CFL-type restriction on the height of the space-time prism ∆t,

k i,n = - ∆t 4 a • n i - |E| 3 ≤ 0. (8)
Although this time step is generally larger than a more traditional CFL condition in an upwind finitevolume method -as we show it in the Appendix -it is an undesirable condition for a numerical discretisation that is inherently implicit. The properties of a given distribution depend to a great degree on the precise definition of the coefficients. The typical requirements for a broadly successful RD scheme include [START_REF] Karman | Institute of Fluid Dynamics[END_REF] • positivity, which warrants that the numerical approximations are free of spurious oscillations;

• linearity preservation, which ensures that a (k -1)th-order polynomial representation leads to a kth-order accurate scheme; • conservation, which guarantees that discontinuities are captured correctly; • compactness, which is primarily for computational efficiency and requires that the cell/prism residual be distributed to its own vertices only; • continuous dependence of the coefficients, which enhances the iterative convergence of the algebraic solver; and • upwinding, which dictates that the discretised model propagates information in the same direction and at the same velocity as its non-discretised counterpart. By Godunov's theorem [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Wesseling | Principles of computational fluid dynamics[END_REF], only nonlinear schemes can satisfy the conditions for both positivity and linearity preservation. Very often, though, linear schemes are used as the basis for constructing these nonlinear schemes. If a is a conservative linearisation of a(u) = ∂f /∂u over the space-time prism, then we can compute the discrete residual as

φ Et = i∈E k i,n u n i + i∈E k i,n+1 u n+1 i , (9)
and the resulting scheme will be conservative as long as the constraint following (5) is satisfied [START_REF] Deconinck | A multidimensional generalization of Roe's flux difference splitter for the Euler equations[END_REF].

Otherwise, the prism residual needs to be computed by a sufficiently accurate quadrature rule in order to obtain conservation, even when (5) holds. This latter version is also known as the conservative RD (CRD) formulation [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF].

There are a large number of numerical schemes that can be (re)cast in the RD framework. We now describe only the space-time upwind schemes that are the focus of our investigation in this article. For this, the parameters defined in [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] are used.

• The space-time N (STN) scheme [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF] is here defined so that it is suitable for the CRD formulation [START_REF] Ricchiuto | Residual distribution for general timedependent conservation laws[END_REF],

u in t = N t   j∈E k + j,n u n j + j∈E k + j,n+1 u n+1 j -φ Et   , φ E i,n N = k + i,n u n i -u in t , φ E i,n+1 N = k + i,n+1 u n+1 i -u in t . (10)
This is a linear scheme that has all the desired properties except linearity preservation. We note, however, that positivity is only formally guaranteed if we use [START_REF] Deconinck | Status of multidimensional upwind residual distribution schemes and applications in aeronautics[END_REF] for the computation of φ Et ; see [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF] for space-time DRD schemes. Nevertheless, oscillation-free behaviour is also often observed in computational experiments [START_REF] Deconinck | Residual distribution schemes: foundation and analysis[END_REF] for the CRD scheme, when a quadrature rule is used to evaluate φ Et . • The space-time LDA (STLDA) [START_REF] Roe | Linear advection schemes on triangular meshes[END_REF][START_REF] Ricchiuto | Residual distribution for general timedependent conservation laws[END_REF] scheme is defined as

φ E i,n LDA = k + i,n N t φ Et , φ E i,n+1 LDA = k + i,n+1 N t φ Et , (11)
which is also a linear scheme and has all the desired properties except positivity. • The space-time blended (STB) scheme combines two linear schemes, typically the STN and STLDA schemes, through a nonlinear blending coefficient,

φ E i,n B = θ φ E i,n N + (1 -θ) φ E i,n LDA , φ E i,n+1 B = θ φ E i,n+1 N + (1 -θ) φ E i,n+1 LDA . (12)
The blending coefficient determines how 'well' the required properties, especially positivity, are satisfied. The definition adopted in this article [START_REF] Deconinck | Status of multidimensional upwind residual distribution schemes and applications in aeronautics[END_REF] is

θ = |φ Et | i∈E φ E i,n N + i∈E φ E i,n+1 N . ( 13
)

Discontinuity in time for the scalar equation

When the representation of the discrete solution is allowed to be discontinuous in time, one needs to introduce additional residuals. These can be viewed as prism residuals with ∆t → 0, see [START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF][START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF], and thus read

ψ E = lim ∆t→0 φ Et = lim ∆t→0 ∂Et f t • n t d(∂Ω t ) = ∂E [u n h ] dΩ = |E| 3 i∈E [u n i ] , (14)
where [•] represents the jump through the face. The simple vertex-centred distribution,

ψ E i,n = |E| 3 [u n i ] , (15)
was shown in [START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF][START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF] to be positive and linearity preserving, and this is the formulation we use throughout this article. As a result, in the case when only discontinuities in time are introduced, the algebraic system takes the form

E∈Di

φ E i,n + ψ E i,n = 0, E∈Di φ E i,n+1 = 0, (16)
∀i ∈ Ω h .

Systems of equations

Much of the RD framework described for the scalar equation can be applied to the nonlinear hyperbolic system of conservation laws,

∂ t U + ∇ • F = 0 or ∂ t U + A(U ) • ∇U = 0, (17)
where A(U ) = [A x , A y ] = [∂F x /∂U, ∂F y /∂U ] = ∂F /∂U is the wave-speed tensor. The two key elements where the extension is not so straightforward are positivity and upwinding. Positivity of the general system ( 17) is not clearly defined; instead, an algorithmic property that preserves the nonnegativity of certain physical quantities -water height, density, pressure, concentration, etc. -are required.

The definitions of the upwind directions for [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF] are not immediately apparent, either. Different approaches have been proposed, the most popular of which is matrix distribution [START_REF] Van Der Weide | Euler and Navier-Stokes solvers using multi-dimensional upwind schemes and multigrid acceleration[END_REF], which we also adopt in the current work. For space-time prisms, the inflow parameters used for the prism distribution are then defined as

K i,n = - ∆t 4 A • n i - |E| 3 I, K i,n+1 = - ∆t 4 A • n i + |E| 3 I, ( 18
)
where I is the identity matrix and A represents a prism-averaged state of the flux Jacobian A.

Conservative linearisation for systems of equations

As in the scalar case, A is preferably derived from a conservative linearisation. For systems of equations, however, the exact form of the conservative linearisation is rarely straightforward. It often assumes linear variation in quantities other than the conservative variables U . The set of linearly varying quantities is traditionally called the Roe-parameter vector. The (vector-valued) discrete residuals of the system [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF] can then be computed as

Φ Et = i∈E K i,n U n i + i∈E K i,n+1 U n+1 i with U n i = ∂U ∂Z Z n i , U n+1 i = ∂U ∂Z Z n+1 i , ( 19
)
where Z is the vector of the Roe-parameter variables [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Roe | Fluctuations and signals[END_REF] and the derivative matrix ∂U ∂Z is evaluated at the arithmetic mean states of the parameter-vector variables over the space-time prism.

If the conservative linearisation is not known or its implementation is not practical, it is also possible to achieve conservation by using a sufficiently accurate quadrature rule to compute the prism residual [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF].

Upwind RD schemes for systems of equations

Assuming that the inflow matrices defined in [START_REF] Hubbard | Discontinuous fluctuation distribution[END_REF] are diagonalisable, as is the case in almost all physical applications, we have K i,n = RDR -1 , where D is the diagonal matrix of eigenvalues, R -1 is the matrix of the left and R of the right eigenvectors. Defining D ± = 1 2 (D ± |D|) with |D| denoting the absolute values of the entries, we can introduce the other inflow matrix parameters (cf. (7)),

K + i,n = RD + R -1 i,n , K + i,n+1 = RD + R -1 i,n+1 , K - i,n = RD -R -1 i,n , K - i,n+1 = RD -R -1 i,n+1 , N t = i∈E K + i,n + i∈E K + i,n+1 -1 . (20)
These, in turn, are used to define the upwind RD schemes along the lines of the scalar case ( 10)-( 12).

• The space-time N (STN) scheme for systems is now defined as

U in t = N t   j∈E K + j,n U n j + j∈E K + j,n+1 U n+1 j -Φ Et   , Φ E i,n N = K + i,n U n i -U in t , Φ E i,n+1 N = K + i,n+1 U n+1 i -U in t , (21)
where U is defined as in [START_REF] Klaij | Pseudo-time stepping methods for spacetime discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[END_REF] if a conservative linearisation exists and thus Φ Et is computed as in [START_REF] Klaij | Pseudo-time stepping methods for spacetime discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[END_REF]. On the other hand, U = U if the CRD formulation is used and Φ Et is computed by means of a quadrature rule. In the system case, positivity should still be understood as positivity of the underlying linear advection. For nonlinear systems, positivity in the strict sense is not and should not be satisfied. Instead, the more generally property that the solution is free of spurious oscillations is required. • The space-time LDA (STLDA) scheme for systems is defined as

Φ E i,n LDA = K + i,n N t Φ Et , Φ E i,n+1 LDA = K + i,n+1 N t Φ Et . (22)
• The space-time blended (STB) scheme for systems is still defined as a combination of STN and the STLDA schemes,

Φ E i,n B = Θ Φ E i,n N + (I -Θ) Φ E i,n LDA , Φ E i,n+1 B = Θ Φ E i,n+1 N + (I -Θ) Φ E i,n+1 LDA . (23)
The matrix blending parameter Θ can now be computed in a number of different ways [START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Sermeus | An entropy fix for multi-dimensional upwind residual distribution schemes[END_REF]. The direct analogue of ( 13) is

Θ 1 = diag   |Φ Et | i∈E Φ E i,n N + i∈E Φ E i,n+1 N   , (24)
where the division should be understood as an elementwise operation. Another possibility, proposed in [START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF], is to choose a particular direction ξ = (ξ x , ξ y ) and compute the decomposition

A • ξ = R ξ D ξ R -1 ξ .
The blending ( 23) is then carried out on the 'characteristic' residuals

Φ N i,n = R -1 ξ Φ N i,n , Φ N i,n+1 = R -1 ξ Φ N i,n+1 , Φ LDA i,n = R -1 ξ Φ LDA i,n , Φ LDA i,n+1 = R -1 ξ Φ LDA i,n+1
, with the blending parameter computed as

Θ 2 = diag i∈E Φ N i,n + i∈E Φ N i,n+1 i∈E Φ N i,n + i∈E Φ N i,n+1 , (25)
where we drop the superscript 'E' to avoid clutter. Finally, we calculate the blended residuals based on the original variables by

Φ B i,n = R ξ Φ B i,n and Φ B i,n+1 = R ξ Φ B i,n+1 .

Discontinuity in time for systems of equations

The treatment of the discontinuity in time is directly analogous to [START_REF] Hubbard | A framework for discontinuous fluctuation distribution[END_REF], and is defined as [START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF]]

Ψ E i = |E| 3 U n i,n , (26)
which leads to the algebraic equations

E∈Di

Φ E i,n + Ψ E i,n = 0, E∈Di Φ E i,n+1 = 0, (27)
∀i ∈ Ω h , which needs to be solved at each time step.

Application to the shallow-water system

In this article, we apply the space-time RD framework for systems of equations to the frictionless shallowwater equations with non-flat bottom topography, i.e. we seek solution to the system

∂ t U + ∇ • F (U ) + S(U ) = 0 on Ω t , (28)
where

U =   d du dv   , F = F x F y =    du dv du 2 + gd 2 2 duv duv dv 2 + gd 2 2    , S = -    0 gd ∂b(x,y) ∂x gd ∂b(x,y) ∂y    . ( 29
)
Here d is the water height, u = (u, v) is the flow velocity and b is the height of the bottom topography.

The level of the free surface is defined as η = d + b. In the case of b ≡ 0, the nonlinear system [START_REF] Roe | Fluctuations and signals[END_REF], with the variables (29), expresses the conservation of water height d and discharge du -it can, therefore, admit discontinuous solutions (hydraulic jumps). In order to capture these discontinuities correctly, the discretisation of ( 28) should also be done in a conservative manner.

Conservative schemes

There are essentially two different ways to achieve a conservative RD scheme for [START_REF] Roe | Fluctuations and signals[END_REF]. One is to derive a conservative linearisation for the flux Jacobian by assuming linear variation of the Roe-parameter vector [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Roe | Fluctuations and signals[END_REF]. The other is to use the CRD formulation [START_REF] Csík | A conservative formulation of the multidimensional upwind residual distribution schemes for general nonlinear conservation laws[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF], where a nonconservative linearisation of the conservative variables is used and conservation is achieved by computing the prism residuals via quadrature rules.

A conservative linearisation

In the steady state, conservative linearisation of (28) can be obtained by assuming the linear variation of the Roe-parameter vector

Z = 1 √ d U . Then ∂U ∂Z =    2 √ d 0 0 √ du √ d 0 √ dv 0 √ d   
so that the evaluation of A at the following cell-averaged states [26], together with [START_REF] Klaij | Pseudo-time stepping methods for spacetime discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[END_REF], result in a conservative discretisation,

u = √ du √ d , v = √ dv √ d , d = 0.9d + 0.1 √ d 1 d 2 d 3 √ d . ( 30
)
Here • signifies the arithmetic mean of the values at the vertices (represented by the subscripts 1, 2, 3) of the triangle.

For the space-time prism, one can use the linearisation (30) at the bottom and top of the prism, and apply the trapezium rule in time to achieve conservation.

The CRD approach

As an alternative to conservative linearisation, one can instead use the arithmetic means d, u, v (over either triangles or prisms), and apply the CRD formulation. In this case, however, the positivity of the N and STN schemes is formally lost, even though it is often observed in numerical experiments.

Hydrostatic balance: the C-property

The shallow-water system (28)-( 29) is, by construction, in hydrostatic balance as this is one of the main assumptions in its derivation from the Navier-Stokes equation [START_REF] Stoker | The mathematical theory with applications[END_REF]. Almost all numerical discretisations preserve this property exactly over a flat bed. Over a general (even smooth) non-flat bed, however, the exact preservation of the hydrostatic balance (also called the C-property) is not always straightforward.

The CRD approach

As it was proved in [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF], as long as the water-height d is assumed to be linearly varying and so is the bottom topography b, linearity preserving CRD schemes also satisfy the C-property [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. This simply means that the source term's contribution to the prism residual needs to be computed as

Φ b Et = ∆t gd 2 i∈E 0 b i n i . ( 31
)
The conventional definition of the STN scheme ( 21), however, does not satisfy the hydrostatic balance, which means that the corresponding B and STB schemes will also lack this property. So instead of ( 21) we need to apply a slightly modified version of the scheme. For the CRD formulation, it means that U = [d, du, dv] T is replaced by V = [η, du, dv] T so that we have

V in t = N t   j∈E K + j,n V n j + j∈E K + j,n+1 V n+1 j -Φ Et   , Φ E i,n N = K + i,n V n i -V in t , Φ E i,n+1 N = K + i,n+1 V n+1 i -V in t . (32)
Proposition. The STN scheme (32) exactly satisfies the hydrostatic balance (C-property). As a result, the STB scheme defined in (12) also exactly satisfies the same property.

Proof. In each space-time prism, with dropping the superscript Es, the STN scheme (32) can be rewritten as [START_REF] Ricchiuto | Residual distribution for general timedependent conservation laws[END_REF] Φ

N i,n = Φ LDA i,n + K + i,n N t j∈E K + j,n V n i -V n j + K + i,n N t j∈E K + j,n+1 V n i -V n+1 j , Φ N i,n+1 = Φ LDA i,n+1 + K + i,n+1 N t j∈E K + j,n V n+1 i -V n j + K + i,n N t j∈E K + j,n+1 V n+1 i -V n+1 j .
The fluctuations Φ LDA i,n and Φ LDA i,n+1 are exactly zero for V = [η, du, dv] T = [const, 0, 0] T as long as ( 31) is used to compute the source term's contribution to the prism residual. Each of the other terms is exactly zero because the assumed linear variation of U = [d, du, dv] T and b implies the linear variation of

V = [η, du, dv] T , i.e. V = [η, du, dv] T = [const, 0, 0] T .

Conservative linearisation

We are not aware of a conservative linearisation for the shallow-water equations that preserves the hydrostatic balance exactly as well as conserves water height and discharge.

Since the linear variation of √ d is assumed in the conservative linearisation proposed in this work, it is natural to assume the linear variation of √ b, too, and compute the source term's contribution in the prism residual as

Φ b Et = ∆tg d i∈E 0 √ b √ b i n i . ( 33
)
This formulation introduces a small error in the hydrostatic balance because in general

η i = d i + b i = const √ d √ d i + √ b √ b i = const.
Conversely, it is possible to achieve exact hydrostatic balance by assuming the linear variation of √ d and √ η (instead of √ d and √ b) and compute the source term's contribution in the prism residual as

Φ b Et = ∆t g d 2 i∈E   0 √ η η n i + η n+1 i - √ d d n i + d n+1 i n i   . ( 34
)
This formulation preserves the lake at rest, V = [η, du, dv] T = [const, 0, 0] T , exactly but in general introduces an error in the conservation of water height and discharge because the term

i∈E   0 √ η η n i + η n+1 i - √ d d n i + d n+1 i n i  
does not stay constant in time even when b does. We emphasise, however, that in both cases the errors are rather small and in most computational experiments they do not affect the results substantially.

Implementation details

In the numerical examples presented in this article, we assume the discrete representation to be discontinuous in time only. In this case, we need to solve the algebraic system [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF], which is the space-time DRD discretisation of (28), at each time step. For this a simple pseudo-time-stepping algorithm is used

U n i U n+1 i m+1 = U n i U n+1 i m - τ s t i E∈Di Φ E i,n + Ψ E i,n Φ E i,n+1 m , (35)
where s t i = ∆ts i is the volume of the dual space-time cell (with s i being the volume of the spatial dual cell). The pseudo-time step τ is given by τ = 0.9 min

i s t i E∈Di ρ(K + i,n+1 ) , ρ(K + i,n+1 ) = max diag D + i,n+1 , (36)
with ρ(M) denoting the spectral radius of a given matrix M. In all the computations, we use the stopping criterion

rel_tol = Ψ n , Ψ n+1 T m 1 (Ψ n , Ψ n+1 ) T 0 1 < 10 -3 .
Often, though, we also apply an intermediate criterion to freeze the blending parameter, i.e. to stop its recomputation and continue the iteration with a constant value. In our numerical experiments, this happens once rel_tol = 10 -1.5 is reached.

Only results with the CRD formulation are reported. We note, however, that for flat bottom topography, the scheme with the conservative linearisation provides identical results. Even when the bottom topography is not flat, we only observe significant difference between the results provided by the two different formulations when the exact satisfaction of the hydrostatic balance is investigated (cf. Section 5.3).

Numerical results

In this section, we present a number of time-dependent test cases -most of which are nonlinearto validate the discontinuous schemes presented in the previous sections. Throughout the section, the representation is discontinuous in time, which results in an unconditionally stable and globally positive scheme.

Circular dam break over wet bed

Dam-break examples are used to experimentally assess whether the numerical solution is free of spurious oscillations. The first dam-break example we consider is that of a circular dam with radius r = x 2 + y 2 = 60 separating water levels d = 10 and d = 0.5 in a square basin. The computational domain is the top right quarter, Ω = [0, 100] 2 , of the entire basin with solid-wall boundary conditions. The blending parameter Θ 2 is used in this example.

The velocity field tends to be more sensitive to spurious oscillations than the water height so we plot Froude-number contours for the STB scheme.

Figure 1 shows that the STB scheme captures the discontinuity well without introducing spurious oscillations for CFL = 1, CFL = 2 and CFL = 4. Interestingly, however, the schemes with CFL = 2 and CFL = 4 seem to be less diffusive than the scheme with CFL = 1. Three-dimensional plots of the water height are given in Figure 2 and they reveal no discernible oscillations for CFL = 1, 2, 4, either.

Improving the iterative convergence

Any increase in the time step for any implicit space-time scheme will only translate into gains in computational work if the schemes with different CFL numbers perform similarly in terms of iterative convergence at each time step. We use this test to investigate this property when the STB scheme is applied. At each time step, we integrate until a relative tolerance rel_tol = 10 -3 is reached. This is typically more than enough to achieve the accuracy required by a space-time RD discretisation because we can use the result from the previous time step as an initial guess. We plot the convergence histories in Figure 3 for different CFL numbers for two slightly different blending strategies. In both cases the blending parameter Θ 2 is applied, but in the second instance we freeze the value after reaching the relative tolerance rel_tol = 10 -1.5 and carry on with the iteration using this constant value. It is clear from the results that the step from CFL = 4 to CFL = 8 does not overall lead to any reduction in computational work. Doubling the time step is obviously beneficial, and moving to CFL = 4 can still provide efficiency gains as long as we freeze the blending parameter (at rel_tol = 10 -1.5 in this example).

Lake at rest: the preservation of the hydrostatic balance

This time-dependent example is to experimentally verify that the modified STN scheme [START_REF] Stoker | The mathematical theory with applications[END_REF] does indeed preserve the C-property, i.e. it is hydrostatically well-balanced [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF][START_REF] Xing | High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. The blended scheme should inherit this property since it is a linear interpolation between the STN and the STLDA schemes. Also, the actual form of the blending parameter does not influence the lake-at-rest property so for the sake of brevity we only show results for the STB scheme with Θ 1 , which we freeze after rel_tol = 10 -1.5 . In the numerical tests we assume a smooth bathymetric function, b(x, y) = 0.8e -5(x-0.9) 2 -50(y-0.5) 2 ,

and 'still-water' initial conditions, [η 0 , u 0 , v 0 ] = [1, 0, 0], over the domain Ω = [0, 2] × [0, 1], using a mesh with a typical edge resolution of h = 0.01. Table 1 shows the errors after integrating until t = 0.5.

The errors remain around machine precision until the end of the time marching, which indicates a well-balanced scheme. The numbers are for test runs with solid-wall boundary conditions but weak characteristic boundary conditions provide near identical results. As a second variant of this test case, we put a perturbation on the initial, 'still-water' condition,

η 0 = 1.01 if 0.05 < x < 0.15 1 otherwise . ( 38
)
For this example, the boundary conditions are weakly enforced everywhere: symmetric (i.e. solid wall) for the bottom and top boundaries, and characteristic freestream for the left and right boundaries including all corners. Figures 4-7 plot 50 contours of the free surface η for CFL = 1, 2, 4, 8, respectively, at times t = 0.24 and t = 0.48, while Figure 8 shows the same along the line y = 0.5. They reveal that qualitatively the scheme with CFL = 1 and CFL = 2 appear to be the most accurate while increasing the CFL number further results in the solution being more diffusive. This phenomenon is partly because of the loss of accuracy associated with a large time step relative to the mesh size, and partly because of the lack of local positivity of space-time implicit schemes.

Nevertheless, some of the most important qualitative properties are satisfied for all CFL numbers: a) they preserve the C-property in front of the perturbation; and b) they capture the interaction between the gravitational wave and the non-flat bottom; c) they settle back to lake at rest after the wave has passed.

Travelling vortex

To evaluate the accuracy and (grid) convergence properties of the STLDA and STB schemes, we include the example of a travelling vortex with known exact solution [START_REF] Dobeš | Second order blended multidimensional upwind residual distribution scheme for steady and unsteady computations[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. Given a flat bottom topography, the exact velocity field is expressed as u ∞ + u , with

u = Γ (1 + cos(ωr c )) (y c -y, x -x c ) if ωr c < π (0, 0) otherwise ,
and u ∞ being constant. The constant Γ is the vortex intensity parameter, (x c , y c ) are the coordinates of the centre of the vortex, r c is the distance from the centre of the vortex, and ω is the angular wave frequency associated with the diameter of the vortex. The water height is then given as

d(r c ) = d ∞ + 1 g Γ ω 2 (κ(ωr c ) -κ(π)) if ωr c < π 0 otherwise , with κ(x) = 2 cos(x) + 2x sin(x) + 1 8 cos(2x) + x 4 sin(2x) + 3 4 x 2 and d ∞ = 1.
For the grid-convergence study, we set u ∞ = (6, 0), Γ = 15, ω = 4π, g = 9.80665 and use a sequence of five unstructured triangulations of the domain Ω = [0, 2] × [0, 1] with characteristic mesh sizes h = 1/10, 1/20, 1/40, 1/80, 1/160, respectively. At the initial state the centre of the vortex is at (x c , y c ) = (0.5, 0.5) and the time marching stops at t = 1/6, when (x c , y c ) = (1.5, 0.5). Freestream characteristic boundary conditions are used everywhere.

Figure 9 shows grid convergence of the STLDA and STB schemes with CFL = 1, 2, 4, 8. Second-order accuracy is observed for the STLDA scheme, although in the case of CFL = 8 this is only reached at the finest mesh. Also, it is clearly between CFL = 4 and CFL = 8 that the larger time step has a significant effect on the accuracy of the scheme.

The convergence rate for the STB scheme is slightly suboptimal -at around 1.8 -but still better than existing results of nonlinear space-time RD schemes [START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Hubbard | Discontinuous upwind residual distribution: A route to unconditional positivity and high order accuracy[END_REF]. The blending in this example is applied to the residuals of the characteristic variables Θ 2 .

Partial dam break

This example is similar to the one in Section 5.1 but has a more complex geometry. We consider the domain An unstructured mesh with characteristic mesh size of h ≈ 2 is used and solid-wall boundary conditions are imposed everywhere. The blending parameter is defined on the characteristic values, Θ 2 , and we also freeze the parameter once rel_tol = 10 -1.5 is reached. Figure 10 shows water-height contours at the end of the time integration t = 7.2, while Figure 11 shows slice plots along the line y = 135. The schemes capture both the rarefaction wave (left of the break) and the shock wave (right of the break) with acceptable accuracy and compare favourably to published results in the literature [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF]. As in the previous examples, however, the scheme with CFL = 8 is markedly less accurate than the schemes with lower CFL numbers.

Circular dam break over nonsmooth bed

This test case considers a two-dimensional variant of the Riemann problem over discontinuous bottom topography, proposed in [START_REF] Bernetti | Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry[END_REF]. The computational domain is now Ω = [0, 30] = max Θ 2 , is applied to all variables of the residual in order to achieve an additional stabilising effect. The simulation follows the wave hitting the underwater wall, then partially reflecting from it and partially moving forward and exiting the domain. There is also a stationary shock wave along the discontinuity of the bed.

30 contours of the free surface η are depicted in Figure 12 for CFL = 4. The figure shows four snapshots of the solution at intervals of exactly 2.5 in time. All three waves -the outgoing, the reflected the stationary -are well captured. The contour plots for the STB schemes with other CFL numbers are omitted for this example because they show very similar behaviour to what we observe in the previous dam-break problems: the schemes with CFL = 2 and CFL = 4 provide the best qualitative results while the one with CFL = 1 is more diffusive and the one with CFL = 8 is both more diffusive and more oscillatory. This general pattern is also apparent from Figure 13, which shows slice plots for all CFL numbers along the diagonal, x = y, of the domain.

Comparison with existing schemes

This section compares the space-time scheme developed in this work with two other RD discretisations that exist in the literature. Both of these schemes are implicit but also require a time-step restriction because the representation in continuous in time as well as in space. The development of the explicit RKRD scheme [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: second order case[END_REF] for the shallow-water equations is ongoing work [START_REF] Ricchiuto | Explicit residual discretizations for shallow water flows[END_REF]. That scheme and its comparison to some implicit formulations will be reported elsewhere. The two other schemes we consider from the literature are the following.

1. The blended LDA-N scheme where the LDA scheme is defined as in [START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF]. This is a similar upwind blended scheme to what is developed in this article but it is not a space-time formulation. 2. The stabilised LLF scheme in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. This is currently one of the most robust RD schemes for timedependent shallow-water simulations. It is essentially a central scheme and it is computationally cheaper than an upwind schemes that is continuous in both space and time. We choose three test cases to compare the space-time scheme investigated here with the two other schemes. The first is the travelling vortex in Section 5.4 to compare the computational work against accuracy. The second is the partial dam break in Section 5.5 and the third is the circular dam break over non-smooth bed topography in Section 5.6. We measure computational work in the total number of pseudo-time iterations over the entire time integration. This, however, does not reflect the fact that the upwind blended schemes are computationally more expensive per space-time prism than the LLFs scheme. The blended scheme requires about four times as much computational work per prism as the LLFs schemes, while the discontinuous STB schemes requires about eight times as much. To reflect this, we define one work unit as being the amount of computational work the LLFs scheme needs per space-time prism.

Travelling vortex

This test case is used to compare different RD schemes based on the computational work needed to achieve a given accuracy. The set-up of the test case is the same as in Section 5.4. Tables 2, 3 and 4 show the computational performance of the blended scheme, the LLFs scheme and the discontinuous STB scheme with CFL = 4. The computational work is measured both in the total number of pseudo-time iterations and in work units defined above. Figure 14 plots accuracy achieved for computational work. The results indicate the discontinuous STB scheme (with CFL = 4 in this case) is the most computationally efficient overall.

Partial dam break

We use this example, which is the same as the one in Section 5.5, to compare the computational performance of the space-time scheme with the two other schemes when discontinuities in the solution have to be captured. Figure 15 shows the same water-height contours for these two schemes as Figure 10 for the space-time blended scheme. The blended LDA-N scheme of Ferrante appears to give the better result of the two. The STB scheme with CFL = 1, 2, 4 are of comparable quality to these but the one with CFL = 8 is clearly inferior. Figure 16 shows slice plots to directly compare two of the discontinuous STB schemes with these methods.

The computational performance of the two schemes and that of the STB scheme with CFL = 1, 2, 4, 8 is listed in Table 5. The computational work is measured in three ways until the final time of the simulation, t = 7.2, is reached: as the number of physical-time steps; as the number of total pseudo-time iterations; and as the number of work units defined at the beginning of this section. By the first two measures, the discontinuous space-time scheme outperforms all other schemes even with CFL = 1. This is in part because the past-shield condition is often (but not always) larger than the CFL-type restrictions used in other schemes. But it is also because the iterative convergence of the STB scheme is relatively good. In terms of work units, the LLFs scheme is still less expensive but it is also less accurate, though accuracy can only be assessed qualitatively for this example.

In particular, moving from CFL = 1 to CFL = 2 provides the most gains by more than halving the total computational cost. The benefits of increasing the time step from CFL = 2 to CFL = 4 are less obvious if the mild deterioration in the quality of the numerical simulation is also taken into account. The performance of the STB scheme with CFL = 8 is clearly poor. In this test case, it seems to generate spurious modes that spoil both the iterative performance of the scheme and quality of the results.

Circular dam break over nonsmooth bed

We use this relatively difficult dam break problem to compare the performances of the different schemes on both a quasi-uniform and a locally refined mesh. First, we compare the schemes on the same quasi-uniform mesh as used in Section 5.6. Table 6 lists the computational work associated with the continuous-in-time blended scheme, the LLFs scheme, and the space-time schemes with CFL = 1, 2, 4, 8. It shows that the computational cost, measured as the total number of pseudo-time iterations, is reduced to about a third of the LLFs scheme and to about half of blended scheme with the LDA part defined as in [START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF]. Second, since there is a stationary shock wave over the discontinuous bed topography, it is natural to use local refinement there. The typical edge length in the region where the refinement takes place is 0.1, which is one-third of the typical edge length elsewhere in the domain. We show results for the same test case but with a mesh that is locally refined where the stationary shock occurs. The space-time scheme is now run with CFL = 9, which would approximately correspond to CFL = 3 in the region where the mesh is not refined. show 30 free-surface contour plots for the blended scheme of Ferrante [START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF], for the LLFs scheme [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF], and for the discontinuous space-time scheme developed here with CFL = 9, respectively. Slice plots of the free surface along the line x = y are shown in Figure 21, while 3D free-surface plots at t = 7.5 are shown in Figure 20.

The space-time scheme seems to be less oscillatory than the two other schemes studied in this test case. In particular, the LLFs scheme develops rather large spurious spikes along the stationary shock in the second half of the integration time interval (see the bottom two subfigures in Figure 18 and also the 3D plots in Figure 20). By contrast, the space-time DRD schemes do not suffer from these artifacts, albeit at the price of a somewhat larger numerical diffusion. Overall, the space-time DRD schemes appear to represent the physical processes most accurately. Table 7 lists the computational work for five schemes, including three space-time DRD schemes. The space-time schemes need a relatively large number of iterations to converge in each physical-time step .95 for the blended scheme and the LLFs scheme at the end of the integration t = 7.2. This figure is directly comparable to Figure 10 for the space-time DRD scheme. so the benefits to computational cost are less pronounced in this example. They are still about twice as computationally efficient as the other schemes when measured in total number of pseudo-time iterations.

When measured in number of physical-time steps, the gain is about 7-10 times. This opens up the possibility of further efficiency gains if the performance of the pseudo-iterative algorithm is improved. Also, this example shows a bigger difference in accuracy between the space-time scheme and the two other schemes (in favour of the space-time schemes) than other dam break problems.

Concluding remarks and outlook

This article applies the framework of discontinuous residual distribution (DRD) to the shallow-water equations with non-flat bottom topography. The focus is on the space-time representation that is discontinuous in time only. This choice is motivated by the fact that discontinuity in time lifts the time-step restriction on the size of the space-time prism and thus results in an unconditionally stable discretisation.

As the numerical experiments demonstrate, we can indeed increase the time step of discontinuous space-time schemes and still retain the most important properties of the shallow-water system: conservation, linearity preservation, upwinding, hydrostatic balance and local positivity. In particular, we emphasise that our interest here is restricted to 'pure' upwinding and therefore no stabilisation term [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] is included in these discretisations.

Comparison to two other implicit RD schemes show that the space-time DRD algorithm developed in this work provides more accurate results (measured either quantitatively or qualitatively) than the currently available RD schemes for the shallow-water equations. The comparisons were carried out on both quasi-uniform and locally refined meshes. Depending on the particular test case and implementation, the best-performing space-time DRD schemes require around 7-10 times fewer physical-time steps than the other schemes. This, on current implementation, translates into 2-4 times fewer pseudo-time iterations over the entire integration. In terms of estimated total computational work, the discontinuous space-time schemes are still less expensive for a given level of accuracy.

Although more accurate than other currently available RD schemes, it has also emerged from the study that the benefits from increasing the time step in the space-time DRD scheme has practical limitations. One of these is the loss of accuracy associated with larger time steps. Another is the relatively large number of pseudo-time iterations needed per physical-time step. It is possible to mitigate these drawbacks by focusing on improving the iterative convergence of the schemes with medium-sized time steps, say between CFL = 3 and CFL = 5 (or higher in case of local mesh refinement) for these genuinely timevarying cases. Improving the iterative convergence by applying a more advanced pseudo-time-stepping algorithm [START_REF] Klaij | Pseudo-time stepping methods for spacetime discontinuous Galerkin discretizations of the compressible Navier-Stokes equations[END_REF] is among the most attractive options. In particular, a further study should compare the computational performance of the space-time DRD scheme with that of the explicit RKRD scheme [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: second order case[END_REF][START_REF] Ricchiuto | Explicit residual discretizations for shallow water flows[END_REF] when applied to the same or similar test cases presented in this work.

Overall, the likeliest areas of application where the proposed scheme could prove advantageous will probably include a relatively large degree of stiffness that comes from the physical equations when the viscous term is included. Figure 17: Circular dam break over discontinuous bed using a locally refined mesh. 30 free-surface contours for the blended scheme where the LDA part is computed as [START_REF] Ferrante | Solution of the unsteady Euler equations using residual distribution and flux corrected transport[END_REF]. Figure 18: Circular dam break over discontinuous bed using a locally refined mesh. 30 free-surface contours for the LLFs scheme taken from [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF].

Case 1: advection along the axis

The upwind FV scheme gets contributions from edges 1, 2, and 3 (red, blue, and yellow in Figure 22). For a vertical speed, the angles between the speed and normals are

θ 1 = π 4 , θ 2 = π 8 , θ 3 = π 2 - π 8 .
One can similarly show that the edge lengths are Case 2: advection along the diagonal

The upwind FV scheme gets contributions from edges 1, 2, and 6 (red, blue, and green in Figure 22). In this case, the angles between the speed and normals are

θ 1 = 0, θ 2 = π 8 , θ 6 = π 2 - π 8
As a consequence, the limiting explicit Euler time step from positivity analysis reads

8 Figure 1 :

 81 Figure 1: Circular dam break. 30 Froude-number contour plots of the discontinuous STB scheme for CFL = 1 (top left), CFL = 2 (top right), CFL = 4 (bottom left) and CFL = 8 (bottom right).

Figure 2 :

 2 Figure 2: Circular dam break. 3D water-height plots of the discontinuous STB scheme for CFL = 1 (top left), CFL = 2 (top right), CFL = 4 (bottom left) and CFL = 8 (bottom right).

Figure 3 :

 3 Figure3: Circular dam break. Convergence of pseudo-time iteration for the space-time DRD scheme with CFL = 1, 2, 4, 8 when the blending parameter is computed at every pseudo-time step (left) and when it is frozen once rel_tol = 10 -1.5 is reached (right).

Figure 5 :

 5 Figure 5: Perturbation to the lake at rest when solved with the discontinuous STB scheme with CFL = 2. 50 water-height contours between values 0.992 and 1.012 at t = 0.24 (left) and t = 0.48 (right) are plotted.

Figure 7 :

 7 Figure 7: Perturbation to the lake at rest when solved with the discontinuous STB scheme with CFL = 8. 50 water-height contours between values 0.992 and 1.012 at t = 0.24 (left) and t = 0.48 (right) are plotted.

  0

Figure 8 : 8 Figure 9 :

 889 Figure 8: Perturbation to the lake at rest when solved with the discontinuous STB scheme with CFL = 1, CFL = 2, CFL = 4 and CFL = 8. Slice plots are shown along the line y = 0.5 at t = 0.24 and t = 0.48.

  [0, 200] 2 with a dam that separates water levels d = 10 and d = 5. The dam is situated in the region [95, 105] × [0, 200] and it breaks between y = 95 and y = 170 at initial time t = 0. The computational domain is thus Ω = [0, 200] \ Ω dam , where Ω dam = ((95, 105) × (0, 95)) ∪ ((95, 105) × (170, 200)).

8 Figure 10 :

 810 Figure 10: Partial dam break. 30 water-height contours between 4 and 9.95 for the discontinuous STB scheme with CFL = 1, 2, 4 and 8.

Figure 11 :

 11 Figure 11: Partial dam break. Slice plots of water height at final time t = 7.2 along the line y = 135 for the discontinuous STB scheme with CFL = 1, 2, 4 and 8.

Figure 12 :

 12 Figure 12: Circular dam break over discontinuous bed. 30 free-surface contours for the discontinuous STB scheme with CFL = 4.

Figure 13 :

 13 Figure 13: Plots of the free surface along the line x = y. The slices are taken at t = 2.5, t = 5.0, t = 7.5 and t = 10.0. In each figure the results for CFL = 1, CFL = 2, CFL = 4 and CFL = 8 are shown.

Figure 14 :

 14 Figure 14: Travelling vortex. Computational work needed to achieve a given accuracy for the LLFs scheme, for the blended non-space-time scheme and for the discontinuous STB scheme with CFL = 4.

Figure 15 :

 15 Figure15: Partial dam break. 30 water-height contours between 4 and 9.95 for the blended scheme and the LLFs scheme at the end of the integration t = 7.2. This figure is directly comparable to Figure10for the space-time DRD scheme.

Figure 16 :

 16 Figure 16: Partial dam break. Slice plots of water height at final time t = 7.2 along the line y = 135 for the blended non-space-time scheme, for the LLFs scheme and for the discontinuous STB scheme with CFL = 2, 4.

l 1 = l 4 2 = l 3 = l 5 = l 6 .

 142356 As a consequence, the limiting explicit Euler time step from positivity analysis reads ∆t(al 1 cos(θ 1 ) + al 2 cos(θ 2 ) + al 3 cos(θ 3 )) = |C| = In the same case, one easily shows that the past-shield condition gives

Figure 19 :

 19 Figure19: Circular dam break over discontinuous bed using a locally refined mesh. 30 free-surface contours for the space-time blended scheme with CFL = 9.

∆t(al 1

 1 cos(θ 1 ) + al 2 cos(θ 2 ) + al 6 cos(θ 6 )) = |C| =In the same case, one again shows that the past-shield condition still gives

Figure 20 :

 20 Figure20: Circular dam break over discontinuous bed using a locally refined mesh. 3D free-surface plots at t = 7.5 for the blended LDAN scheme, the LLFs scheme and the discontinuous STB scheme (CFL = 9).

Figure 21 :Figure 22 :

 2122 Figure21: Circular dam break over discontinuous bed using a locally refined mesh. Plots of the free surface along the line x = y. The slices are taken at t = 2.5, t = 5.0, t = 7.5 and t = 10.0. In each figure the results for the blended LDAN scheme, the LLFs scheme and the discontinuous STB scheme (CFL = 9) are shown.

Table 1 :

 1 Lake at rest. Errors at time t = 0.5 for the discontinuous STB scheme with different CFL numbers.

	CFL = 1.0	η -1	u	v
	L 1	5.405934e -19	1.255336e -16	5.870068e -17
	L 2	8.000922e -18	5.165644e -16	2.647004e -16
	L ∞	2.220446e -16	1.034175e -14	6.178231e -15
	CFL = 2.0
	L 1	1.534435e -18	9.505484e -17	4.305477e -17
	L 2	1.309467e -17	4.054003e -16	1.890289e -16
	L ∞	2.220446e -16	8.781880e -15	4.230023e -15
	CFL = 4.0
	L 1	3.744022e -18	1.078974e -16	5.057157e -17
	L 2	2.054486e -17	4.354319e -16	2.091099e -16
	L ∞	2.220446e -16	8.131857e -15	4.450119e -15
	CFL = 8.0
	L 1	6.180234e -18	1.177755e -16	5.675335e -17
	L 2	2.657800e -17	4.540011e -16	2.285096e -16
	L ∞	2.220446e -16	7.692615e -15	5.239349e -15

  Perturbation to the lake at rest when solved with the discontinuous STB scheme with CFL = 4. 50 water-height contours between values 0.992 and 1.012 at t = 0.24 (left) and t = 0.48 (right) are plotted.

					CFL = 4, t = 0.24										CFL = 4, t = 0.48
	1											1
	0.8											0.8
	0.6											0.6
	0.4											0.4
	0.2											0.2
	0											0
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2		0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	0 0.2 0.4 0.6 0.8 1 Figure 6: 0	0.2	0.4	0.6	0.8 CFL = 8, t = 0.24 1 1.2	1.4	1.6	1.8	2	0 0.2 0.4 0.6 0.8 1	0	0.2	0.4	0.6	0.8 CFL = 8, t = 0.48 1 1.2	1.4	1.6	1.8	2

  + y 2 ) is the radius of the dam. Solid-wall boundary conditions are used at the left and bottom boundaries while homogeneous Neumann at the right and top ones. The characteristic mesh size is h ≈ 0.3 and the time integration stops at time t = 10. The largest value of the blending parameter, θ max

					2 with bathymetric function,
	b(x, y) =	0 0.2	if x + y < 30 otherwise	,
	and initial condition,
	η(x, y) =	1.461837 0.308732	if r < 15 otherwise	,
	where r = (x 2 2

Table 2 :

 2 Travelling vortex. Error and computational work for the blended LDAN scheme with Ferrante's LDA.

	Mesh size	d ex -d h L 2	No. time steps	Total no. iterations	Work units
	1/10	3.9667e-01	41	1061	2.01e+06
	1/20	2.7576e-01	77	1857	1.38e+07
	1/40	9.2256e-02	163	3973	1.18e+08
	1/80	2.2355e-02	330	3957	4.69e+08
	1/160	5.7876e-03	649	6485	3.07e+09

Table 3 :

 3 Travelling vortex. Error and computational work for the LLFs scheme.

	Mesh size	d ex -d h L 2	No. time steps	Total no. iterations	Work units
	1/10	4.1190e-01	40	1223	5.80e+05
	1/20	3.1413e-01	77	2297	4.26e+06
	1/40	1.0290e-01	163	3797	2.82e+07
	1/80	2.6203e-02	330	4340	1.29e+08
	1/160	6.6451e-03	649	7788	9.23e+08

Table 4 :

 4 Travelling vortex. Error and computational work for the discontinuous STB with CFL = 4.

	Mesh size	d ex -d h L 2	No. time steps	Total no. iterations	Work units
	1/10	4.2018e-02	10	376	1.43e+06
	1/20	3.5085e-02	18	781	1.16e+07
	1/40	1.2608e-02	38	1784	1.06e+08
	1/80	3.3537e-03	73	3578	8.49e+08
	1/160	8.6824e-04	155	7414	7.03e+09

Table 5 :

 5 Partial dam break on a quasi-uniform mesh. Computational work is shown for two existing implicit RD schemes and the space-time scheme with CFL = 1, 2, 4, 8.

	Scheme	No. time steps	Average no. iterations	Total no. iterations	Work units
	B-Ferrante	134	29.44	3945	3.59e+08
	LLFs	127	66.39	8432	1.92e+08
	STDRD1	105	34.39	3611	6.57e+08
	STDRD2	54	30.80	1663	3.02e+08
	STDRD4	28	52.32	1465	2.66e+08
	STDRD8	14	132.64	1857	3.38e+08

RR n°7958

Sarmany & Hubbard & Ricchiuto

Acknowledgement

We thank Philip Roe for his contribution to the conservative linearisation of the shallow-water equations.

Appendix

Based on Figure 22, in this Appendix we give an estimate for the CFL-condition of the vertex-centred upwind FV scheme with explicit Euler time discretisation on structured triangulation. This is to be compared with the past-shield condition (8) of the space-time scheme.

Sarmany & Hubbard & Ricchiuto

Table 6: Circular dam break over non-smooth bed using a quasi-uniform mesh. Computational work is shown for two existing implicit RD schemes and the space-time scheme with CFL = 1, 2, 4, 8.

Scheme

No


Limited to one research pioneer and instrument pioneer at the beginning. Then extended to other researchers (in the same and in different fields), to the instrument manufacturer(s), etc. Selection of participants to use the exploratory device. Steering committee limited to pioneers.

Platform

Industry platform = professional rules Industry platform = Production process Industry platform = validation bench for Research field platform = new services for

Table 2 :

 2 Main Features of the Design Process

		Biomaterials				ITRS				Cockpit	Biotech
	Goals	Build a new industry			Support industry growth	Industry growth through innovative	Build new research fields
										cockpits
	Inputs	"Building with Hemp" + a list of	Several issues regarding concepts and	"Future cockpits". A single player at the	New	instrument	concept	(actionable
		potential partners + some isolated	processes i n t h e n o n -competitive	beginning: a cockpit supplier.	phenomenon) + an original research question
		experiments				domain + m a i n i n d u s t r y p l a y e r s
						(suppliers, manufacturers, users)
	New	Yes: characterization of materials,	Yes				Yes: on aircraft integrator utility and	Yes: on the research question (animals and
	capabilities	processes, usages, the construction					constraints, on aircraft companies, on	pathologies, etc.), o n t h e i n s t r u m e n t s
		business								pilots, on technologies, etc.	underlying the phenomena, etc.
	Phases	Sub-sets of one or a few partners	Each partner conducts (or sub-	In-house exploration of conceptual	An original experiment with new instrument
		conducting or financing experiments,	contracts) individual investigations on	alternatives for cockpits.	facility (+publication)
		research projects, validation projects,	clear,	well-identified	technological	Design of cockpit simulators enabling	Other experiments with other researchers,
		communication projects			issues				collective explorations with aircraft	based on the same instrument principle.
										integrator, companies and pilots	Improvement and innovation on the instrument
										principles
	Coordination Regular steering committee sharing	Three conferences per year, with all	Two levels:	Two levels:
		results of the explorations, listing	ITRS group members. Working groups	-inside company: gathering results,	-device level: steering committee (instrument
		open questions to be explored,	responsible for one process issue.	orienting explorations and choosing	manager + pioneer researchers) choosing
		discussing	possible	ways	and	Each group updates its related	partners	relevant experiment to be accepted on the
		collective facility for each related	roadmap, identifies emerging solutions	-simulator level: steering committee	experimental device
		exploration				and pending questions	sharing results with main partners	-experiment level: regular meetings between
										instrument manager and researcher
	Resources	Classic association funding; partners	Members	pay	for	their	own	Cockpit supplier provides the simulation	Initial fund-raising through research contracts.
		pay for their own explorations.		explorations.	Collective	expenses	bench	Experimental devices then financed through
						(meeting) shared		public funding or foundations
	Leading	Leading builders, architects, hemp	Main players of semiconductor industry	Cockpit designer alone at the
	partners	transformers, lime suppliers, grouped	(process suppliers, manufacturers,	beginning. Then component suppliers,
		in an association "Building with	users	such	as	cell	phone	aircraft integrator, aircraft companies,
		Hemp". Growing number of members	manufacturers, PC integrators, etc.).	pilots	also	participate	(slow
		in the association. Limited and	Open participation.		enlargement).
		varying n u m b e r o f m e m b e r s taking	ITRS groups are led by designers	Simulation level steering committee
		part in the association's steering	working		for		semiconductor	limited to aircraft integrator and cockpit
		committee (decreasing first and then	manufacturers.		supplier.
		increasing)