Stéphane Popinet

Graham Rickard

A tree-based solver for adaptive ocean modelling

The development of an adaptive (in space and time) ocean model from an existing adaptive finite-volume Navier-Stokes model is described. A flexible and efficient quadtree spatial discretisation is used which requires collocation of all variables (i.e. an A-grid discretisation). We demonstrate that the use of an approximate projection method allows for implicit damping of instabilities generally associated with the A-grid, at the expense of a relatively small amount of numerical energy dissipation, while accurately preserving dispersive properties and geostrophic balance. The finite-volume formulation also maintains second-order spatial accuracy at all solid boundaries. Test cases demonstrate the efficacy of the adaptive ocean model, and the advantages it has in terms of efficient representation of multi-scale behaviour within a single model. The model is freely available as open-source code.

Introduction

This paper details the extension of an adaptive, finite-volume, three-dimensional, incompressible, Navier-Stokes fluid solver (the Gerris Flow Solver, [START_REF] Popinet | The Gerris Flow Solver[END_REF][START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]) into a dynamical core able to model geophysical fluid flows. Here the focus will be on ocean modelling, but the core can be extended to the atmosphere as well [START_REF] Marshall | Atmosphere-ocean modeling exploiting fluid isomorphisms[END_REF]. The ability of the model to adapt in space and time constrains the structure of the solver. Nevertheless, the tests described here show that many of the fundamental requirements of a geophysical solver can be more than adequately met by this novel adaptive solver.

Why is having adaptivity relevant? In nature, the global dynamics of oceanic or atmospheric flows are controlled by processes occurring on widely different spatial scales. In the ocean, narrow straits control the exchanges between large ocean basins, energetic but thin western boundary currents define latitudinal heat transport and a few geographically small areas control much of the deep ocean circulation through deep water formation. Getting to even smaller spatial scales, the dynamics of near-shore mixing, internal tides, sediments transport and freshwater input is of crucial importance to ecosystems and human activities but depends directly on larger-scale offshore processes (and conversely).

The regular Cartesian spatial discretisations traditionally used in geophysical models are easy to analyse mathematically and easy to implement but lead to a constant spatial resolution (or quasiconstant depending on the type of spherical-to-planar projection) independent from any physicsrelated scale dependence. While large increases in computational resources have allowed ever finer spatial resolutions to be reached using traditional models, the computational efficiency of the numerical schemes is still a very important issue for many applications. A logical way to decrease the computational cost is to discretise the solution using a more optimal spatial distribution of the discretisation elements. Simple applications of this idea to Cartesian meshes lead to grid nesting or stretched grids using curvilinear coordinates. In the wider field of computational fluid dynamics, finite-volume or more rarely finite-element discretisations on unstructured meshes provide almost full flexibility in the choice of local spatial resolution.

It is important to make the distinction between methods where the spatial resolution can vary but is fixed in time (and usually defined a priori by the user) and adaptive methods where the resolution is continuously adjusted according to the properties of the solution. Adaptive methods do not require a priori assumptions about the solution but instead use error estimates to ensure an optimal distribution of discretisation elements. In the context of geophysical fluid flows, adaptive methods have the potential to greatly improve the description of spatially restricted but time-varying processes such as: strong density gradients associated with oceanic or atmospheric fronts, eddying of western boundary currents, tidally-driven small-scale mixing/transport near coastlines etc... A number of recent initiatives have sought to apply the more flexible techniques developed for computational fluid dynamics to geophysical fluid modelling, in recognition of the shortcomings of regular Cartesian discretisation (as well as ad hoc grid nesting). The focus has been mainly on using finite-volume or finite-element methods on statically-refined unstructured grids for barotropic coastal ocean modelling [START_REF] Chen | An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[END_REF][START_REF] Walters | A finite element model for tidal and residual circulation[END_REF][START_REF] Walters | Ocean tides around New Zealand[END_REF]. Unstructured grids allow an accurate description of complex coastlines and the coastal ocean is a natural extension of the river flows for which these models were initially developed [START_REF] Walters | Robust, finite element model for hydrostatic surface water flows[END_REF]. The extension of finite-element or finite-volume models to large-scale oceanic flows is non-trivial because of the need to guarantee specific properties of the large-scale flows. Geostrophic balance must be represented accurately at the discrete level and -more generally -energy must be better conserved because of the low energy dissipation in the open ocean. Cartesian methods have benefited from decades of development and describe the large-scale ocean dynamics accurately. Optimal schemes for finite-volume or finite-element formulations on unstructured grids remain a subject of active research [START_REF] Chen | An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[END_REF][START_REF] Le Roux | Finite elements for shallow-water equation ocean models[END_REF][START_REF] Le Roux | A semi-implicit semi-lagrangian finite-element shallowwater ocean model[END_REF][START_REF] Walters | Comparison of H and P finite element approximations of the shallow water equations[END_REF].

Different adaptive methods have been applied to atmospheric and ocean modelling. For instance [START_REF] Blayo | Adaptive mesh refinement for finite difference ocean models: first experiments[END_REF], [START_REF] Bacon | A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA)[END_REF] and [START_REF] Giraldo | The Lagrange-Galerkin method for the two-dimensional shallow water equations on adaptive grids[END_REF] have studied adaptive techniques, and indeed [START_REF] Bacon | A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA)[END_REF] use their system OMEGA (The Operational Multiscale Environment Model with Grid Adaptivity) for real time hazard prediction. More recent adaptive developments include, e.g., [START_REF] Power | Adjoint goal-based error norms for adaptive mesh ocean modelling[END_REF] (and references therein to their adaptive finite element ocean model), [START_REF] Herrnstein | Structured adaptive mesh refinement using leapfrog time integration on a staggered grid for ocean models[END_REF], [START_REF] Barros | A Global Semi-Implicit Semi-Lagrangian Shallow-Water Model on Locally Refined Grids[END_REF], indicating the on-going level of activity and interest for the application of adaptive methods in geophysical modelling.

Recently, Popinet [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Popinet | The Gerris Flow Solver[END_REF] showed that a hierarchical tree-based spatial discretisation can efficiently provide the flexibility of fully unstructured grids while preserving the simplicity and conceptual framework of Cartesian grids. This work also showed that adapting the discretisation dynamically at each timestep was practical and only represented a small fraction of the total computational cost. The adaptive solver led to several order of magnitude savings in computational time for complex three-dimensional problems.

In this article, we propose to extend the work of [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] to obtain an adaptive method for large-scale oceanic flows with good mass and energy conservation properties and accurate geostrophic balance representation. The overall numerical scheme will be described but we refer the reader interested in the details of the quadtree/octree discretisation and multilevel solver implementation to [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. We will also try to reproduce previously published test cases designed to assess several properties of the numerical scheme: representation of inertia-gravity waves, geostrophic balance and bathymetry. An application illustrating the potential of the adaptive method for practical problems will also be presented.

In the following section the system under consideration is defined. In sections 3 and 4 the spatial and temporal discretisations, respectively, consistent with the Gerris flow solver are described, followed by section 5 detailing our application of the adaptive aspects implied by Gerris for oceanic flows. Section 6 presents a series of applications to two-and three-dimensional test problems, and finishes with an example of adaptivity in the context of multi-scale flow evolution (in this case that of tidally driven flow in Cook Strait, New Zealand). Conclusions are presented in section 7.

Problem definition

Our initial goal is to solve the Boussinesq, incompressible, hydrostatic equations of oceanic motion in a Cartesian reference frame. We will also make a "shallow water" approximation and neglect the coupling terms arising from the vertical motion (Marshall, Hill, Perelman, and Adcroft, 1997b). The equations can then be written as

∂ t U + U • ∇U + w∂ z U = BU + S -∇p -g∇η for the horizontal momentum, ∂ t H + ∇ • HU = 0 for continuity, ∂ t T + U • ∇T + w∂ z T = S T for temperature, ∂ t S + U • ∇S + w∂ z S = S S
for salinity, with U ≡ (u, v) the horizontal components of the velocity,

B ≡ 0 f -f 0
where f is the Coriolis parameter, H ≡ h + η where h is the depth of the ocean and η the surface elevation. S, S T and S S contain the forcing and dissipation terms. U is the depth-averaged fluid velocity. The Boussinesq hydrostatic pressure p is obtained as

p = z 0 -g δρ ρ 0 ,
where δρ is the deviation from the reference density ρ 0 obtained using the equation of state ρ ≡ ρ(T, S, p). The vertical component w of the velocity is obtained by vertical integration of the continuity equation as

w = - z 0 ∇ • U (1)
If we assume that η ≪ h the continuity equation can be linearised as

∂ t η + ∇ • hU = 0
These approximations are representative of current ocean models designed for large-scale flows (as noted by Marshall, Hill, Perelman, and Adcroft (1997b) they are often incorrectly referred to as the "primitive equations" of oceanic motion). The hydrostatic assumption breaks down at scales corresponding to convective overturning in the ocean (≈1 km) (Marshall, Hill, Perelman, and Adcroft, 1997b). The linearised free-surface approximation can cause problems for global tracer conservation [START_REF] Roullet | Salt conservation, free surface, and varying levels: A new formulation for ocean general circulation models[END_REF][START_REF] Campin | Conservation of properties in a free-surface model[END_REF] and becomes inconsistent when the amplitude of free-surface motion is comparable to the depth (e.g. tidal perturbations near coastlines). None of these approximations is a necessary prerequisite for the approach described in this article but they make the initial implementation simpler while preserving all the important physical properties of the oceanic system. In the future we expect to be able to follow the work of Marshall et al. to relax both the hydrostatic and linear free-surface assumptions (Marshall, Adcroft, Hill, Perelman, and Heisey, 1997a). Our starting point being a full 3D Navier-Stokes solver should make this easier.

Spatial discretisation

In this article, we would like to explore how well existing geophysical discretisation techniques fit within the tree-based approach presented in [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The tree-based discretisation was shown to be a good compromise in term of flexibility and computational efficiency. By tree-based we mean either a quadtree structure in two dimensions or an octree structure in three dimensions. Here, the discussion focuses on the quadtree structure, but all properties map equally well into the third dimension.

[Figure 1 about here.]

In the quadtree discretisation, individual cells are classical Cartesian square finite volumes but they are organised hierarchically rather than as an array. Each cell can be the parent of up to four children cells; this is the process at the heart of adaptivity. The resulting recursive structure is a tree (Figure 1) and is easily represented using memory pointers. While this discretisation is not as flexible as a fully-unstructured mesh, it does not require additional geometric information and as importantly preserves the simplicity of the discretisation given by orthogonal Cartesian discretisation volumes.

[Figure 2 about here.] A finite-volume approach (where the solution is averaged on the control volume defined by each cell) imposes an important restriction on the type of discretisation permitted on a quadtree. Let us consider a C-grid discretisation where the quadtree hierarchy is centred on the pressure (Figure 2). When refining the grid, the new control volumes for each velocity component do not necessarily delimit a full subset of the coarser control volumes (the single-hashed area in Figure 2.b is covered by the coarser control volume but not by any of the finer control volumes). It is thus not clear how to consistently and conservatively define the refined values of the momentum components as well as fluxes at coarse-fine boundaries. If one relaxes the finite-volume constraint and uses finite differences instead this problem does not occur but the important conservation properties of finite-volume schemes are lost. A solution would be to maintain independent staggered quadtree hierarchies but this greatly complicates the computation of terms coupling variables defined on different hierarchies. A collocated discretisation (Arakawa A-grid) eliminates this specific problem and is used in the present work in order to allow the use of adaptivity.

Collocated A-grid discretisations are usually avoided in ocean models because of the strong computational modes induced by pressure-velocity decoupling as well as the associated degradation of the dispersive properties of gravity waves. Ocean models traditionally use a staggered Arakawa C-grid or B-grid. C-grids are preferred for finer-scale ocean modelling due to their superior discretisation of the dispersion properties of gravity waves. B-grids show superior discretisation of large-scale inertial waves because they use collocated definitions of the velocity components which allows the Coriolis terms to be computed without any need for interpolation. It is important to note that both B-grids and C-grids suffer from spatial computational modes. In the case of C-grids the computational mode arise from the need to interpolate the values of each component of the velocity to the location of the other component in order to compute the Coriolis term. The resulting computational mode can be problematic for relatively coarse-resolution global ocean modelling [START_REF] Adcroft | A new treatment of the Coriolis terms in C-grid models at both high and low resolutions[END_REF]. B-grids suffer from a pressure computational mode caused by the collocation of the velocity components. Neither discretisation is thus entirely mode-free and filters are often used to keep the computational modes in check. More complex discretisations have been proposed which attempt to solve this problem (C-D grid of [START_REF] Adcroft | A new treatment of the Coriolis terms in C-grid models at both high and low resolutions[END_REF]).

Several notable exceptions that use an A-grid have relied on higher-order discretisations [START_REF] Dietrich | Application of a modified Arakawa 'A' grid ocean model having reduced numerical dispersion to the Gulf of Mexico circulation[END_REF] or other forms of built-in filtering [START_REF] Kar | Stable centered-difference schemes, based on an unstaggered a grid, that eliminate two-grid interval noise[END_REF] in order to control computational modes. More generally, the pressure-velocity decoupling problem of collocated schemes had also been recognised early in the history of incompressible flow solvers [START_REF] Harlow | Numerical calculation of time-dependent incompressible flow of fluid with free surface[END_REF] and more recently solutions have been developed [START_REF] Lai | A projection method for combustion in the zero Mach number limit[END_REF][START_REF] Rider | Approximate projection methods for incompressible flow: Implementation, variants and robustness[END_REF][START_REF] Almgren | Approximate projection methods: Part I. Inviscid analysis[END_REF]. Building on this previous work, we propose in section 5 a semi-implicit barotropic solver on a collocated A-grid.

Vertical discretisation has been the topic of a large number of recent studies [START_REF] Bleck | An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates[END_REF][START_REF] Adcroft | Rescaled height coordinates for accurate representation of freesurface flows in ocean circulation models[END_REF][START_REF] Song | Parametric vertical coordinate formulation for multiscale, Boussinesq, and non-Boussinesq ocean modeling[END_REF]. While our work relies on horizontal collocation, it does not impose restrictions on the vertical discretisation. In section 5.4 we describe a z-level vertical discretisation with an accurate finite-volume description of solid boundaries.

Temporal discretisation

The time discretisation of the equations of motion can theoretically be chosen independently from the spatial discretisation although in practice we will see that some combinations of time-discretisations are much easier to implement using specific spatial schemes.

Advancing the equations of motion in time requires the solution of three independent sub-problems:

1. Pressure update through the barotropic equation.

2. Computation of advection terms in the momentum equation.

3. Computation of Coriolis terms in the momentum equation. Time-explicit discretisation of the barotropic equation is simple but stability requires the resolution of the fast external gravity waves which imposes a strong restriction on the maximum timestep. This problem has been addressed in two ways. The barotropic equation can be advanced independently from the momentum equation using a stable, shorter timestep (split-explicit schemes [START_REF] Shchepetkin | The regional oceanic modeling system (ROMS): A splitexplicit, free-surface, topography-following-coordinate oceanic model[END_REF]). Changing the discretisation to a fully-or semi-implicit scheme removes the constraint on the timestep altogether [START_REF] Dukowicz | Implicit free-surface method for the Bryan-Cox-Semtner ocean model[END_REF] and allows the use of a single timestep for the whole scheme.

Both schemes have advantages and drawbacks. An implicit discretisation requires the resolution of a spatially-coupled linear system at every timestep. This can be expensive if naive algorithms are used, however modern techniques such as multigrid or pre-conditioned conjugate gradients can be very efficient for this type of problem (arguably more efficient than the multiple timesteps of a splitexplicit solution method). Implicit schemes are stable but do not guarantee a physically-meaningful aliasing of the unresolved fast gravity waves. While a split-explicit technique should describe correctly all the external gravity wavelengths, the synchronisation step necessary to reconcile the fast pressure solution and the slow momentum solution adds complexity [START_REF] Shchepetkin | The regional oceanic modeling system (ROMS): A splitexplicit, free-surface, topography-following-coordinate oceanic model[END_REF]. As a result, for both schemes the impact of the coupling of the fast gravity waves with the slow momentum equations is problem-dependent and is unclear in most cases. In this work we use a semi-implicit discretisation with a multigrid solver for the linear system.

A common way to compute the advection term of the momentum equation is to use a leap-frog scheme. This scheme is simple to implement and has the theoretical advantage of providing a purely non-diffusive solution. This solution is not guaranteed to be oscillations-free (monotonic) however and the leap-frog discretisation introduces a temporal computational mode. Some form of explicit filtering (added diffusion) or implicit filtering (e.g. damping provided by the implicit discretisation of the barotropic equation [START_REF] Dukowicz | Implicit free-surface method for the Bryan-Cox-Semtner ocean model[END_REF]) is used in practice, thus cancelling the non-diffusive property of the initial scheme. More recently, ocean modellers have started to use more sophisticated advection schemes initially developed for the solution of compressible flows. These schemes are stable and non-oscillatory even in the presence of strong discontinuities. These properties are obtained through a judicious construction of built-in diffusivity. They are also more computationally expensive than simpler schemes. In our previous work [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] we have successfully used the incompressible version of a two-time-levels, predictor-corrector, unsplit scheme developed by Bell, Collela and Glaz (BCG) [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF] and this is what we make use of in this work.

The Coriolis term can also be discretised using a leap-frog scheme. An alternative is to use an Adams-Bashforth (AB) discretisation which can be adapted to any underlying time-discretisation (leap-frog, predictor-corrector (Marshall, Adcroft, Hill, Perelman, and Heisey, 1997a) or semi-implicit). The AB scheme uses the previous two time-levels to extrapolate the value of the Coriolis term at the time-level required to ensure second-order accuracy in time. In the case of an AB discretisation of the Coriolis terms coupled with the barotropic equation, the stability criterion is non-trivial. We present in Appendix A a detailed analysis of the stability properties of an AB discretisation of the Coriolis terms coupled with a semi-implicit discretisation of the barotropic equation. The analysis is difficult due to the third time-level introduced by the AB scheme. We show that the complex behaviour of inertia-gravity waves gives rise to a stability criterion which can be quite restrictive for this combination.

A simpler and more easily analysed alternative is to use a semi-implicit discretisation of the Coriolis terms, which fits well with a semi-implicit treatment of the barotropic equation. The resulting twotime-levels scheme is easy to analyse and is unconditionally stable. If the velocity components are collocated (e.g. A-or B-grid) the semi-implicit discretisation leads to a compact linear system which can be inverted analytically (see section 5). In the case of a staggered discretisation of the velocity components (C-grid), the semi-implicit discretisation of the Coriolis terms leads to a globally-coupled sparse linear system which needs to be solved numerically at each timestep. This greatly complicates the solution procedure and is the main reason why semi-implicit discretisations of Coriolis terms are avoided on C-grids.

A quadtree, adaptive, semi-implicit scheme

Taking into account the previous discussion and keeping in mind that our starting point is the adaptive incompressible solver presented in [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] we chose to use a two-time-levels, semiimplicit discretisation of both the barotropic and Coriolis terms. If the components of the velocity are collocated, the discrete time evolution of the barotropic momentum can be written

U n+1 -U n ∆t = B[(1 -α)U n + αU n+1 ] -g∇[(1 -α)η n + αη n+1 ] + S n+ 1 2 , ( 2
)
where ∆t is the timestep, n is the time level and α controls the degree of implicitness. All the other source terms (advection terms, viscous diffusion, forcing) are contained in S n+1/2 . Similarly the continuity equation can be written

η n+1 -η n ∆t + ∇ • h[(1 -θ)U n + θU n+1 ] = 0, ( 3
)
where θ is another free parameter. If α and θ are set to 1/2, the discretisation is a second-order accurate Crank-Nicolson scheme. It is also unconditionally stable and conserves the total energy. When expressed in this form, the semi-implicit schemes for the momentum and the continuity equations are coupled through the Coriolis terms. A simpler uncoupled second-order accurate approximation can be obtained through a splitting technique. If we pose

Û ≡ U n+1 + ∆tg∇[α ′ η n+1 -(1 -α ′ )η n ],
where α ′ is another control parameter, we can rewrite (2) as

(I -α∆tB) Û + α∆t 2 gB∇[α ′ η n+1 -(1 -α ′ )η n ] = [I + (1 -α)∆tB]U n -2(1 -α ′ )∆tg∇η n + ∆tS n+ 1 2 ,
where I is the identity matrix. The second term on the left-hand side is clearly second-order in time, moreover if α ′ is set to 1/2 the prefactor is only the (small) relative change in surface-height gradient. This leads to the second-order accurate, explicit approximation for Û

(I -α∆tB) Û = [I + (1 -α)∆tB]U n -2(1 -α ′ )∆tg∇η n + ∆tS n+ 1 2 + O(∆t 2 ). (4)
The resulting 2 × 2 system is easily inverted analytically. If we now pose

U ⋆ ≡ Û + (1 -α ′ )∆tg∇η n , (5)
equation ( 3) can be rewritten

η n+1 -θα ′ ∆t 2 g∇ • h∇η n+1 = η n -∆t∇ • h[(1 -θ)U n + θU ⋆ ]. (6)
Equations ( 4), ( 5) and ( 6) are now entirely decoupled. The full solution procedure can be summarised as follows:

1. Compute Û using U n , η n , S n+ 1 2 and (4)

2. Compute U ⋆ using Û , η n and (5)

3. Compute η n+1 using U n , U ⋆ , η n and ( 6)

4. Compute U n+1 = U ⋆ -α ′ ∆t∇η n+1
The overall scheme is very similar to the projection method (also known as "pressure-correction" method) often used to solve incompressible flow problems [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Peyret | Computational methods for fluid flow[END_REF] or equivalently oceans with a rigid-lid approximation (Marshall, Adcroft, Hill, Perelman, and Heisey, 1997a). In a first phase a provisional velocity field U ⋆ is obtained by adding to U n all the (explicit) source terms of the momentum equation with the exception of the pressure-gradient term. In a second phase the pressure compensating the eventual divergence of U ⋆ is obtained through the solution of a Poisson-like equation and the corresponding pressure gradient term is added to U ⋆ in order to obtain a non-divergent U n+1 . The semi-implicit scheme presented here differs slightly from this outline because a linear freesurface approximation is made rather than a stricter rigid-lid approximation. The free-surface equation (6) (equivalent to the pressure equation for incompressible flows) is then of Helmholtz-type (x + a∇ 2 x = . . .) rather than Poisson-type (∇ 2 x = . . .). In both cases, the spatial coupling induced by the Laplacian operator results in a sparse linear system. In the case of incompressible flows, the solution of this system largely dominates the computational cost of the overall scheme. Much effort has thus gone into efficient solution procedures for this type of systems. Multigrid techniques in particular can achieve an optimal O(n) asymptotic computational cost (where n is the number of unknowns) [START_REF] Brandt | Guide to multigrid development[END_REF][START_REF] Bramble | Multigrid Methods[END_REF].

As noted by ocean modellers [START_REF] Dukowicz | Implicit free-surface method for the Bryan-Cox-Semtner ocean model[END_REF]), the linear system described by equation ( 6) is easier to solve than the Poisson problem of a rigid-lid approximation. This is due to the η n+1 term on the left-hand side of (6) which increases the weight of the diagonal when ∆t is small enough and thus leads to a better-conditioned linear system. If ∆t is large, the added diagonal term becomes small compared to the Laplacian and the solution converges toward the incompressible Poisson solution, with a corresponding stiffness of the linear system. We expect that the multilevel quadtree Poisson solver we developed in [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] will perform well in both limits.

Exact and approximate projections

Up to this point we have only assumed that the velocity components are collocated. The discrete gradient and divergence operators appearing in (4), ( 5) and (6) will depend on the details of the spatial discretisation. To simplify the argument we will consider a one-dimensional problem with g = h = 1. The discrete divergence operator can then be expressed as

∇ i • u = u i+1/2 -u i-1/2 ∆x ,
which is also a consistent finite-volume discretisation (when u i+1/2 and u i-1/2 are interpreted as fluxes through the boundaries of the control volume centred on u i ). For a collocated scheme, u is defined only at integer locations. The half-integer values can be approximated to second-order using a linear interpolation from cell centres. This gives the second-order accurate discretisation

∇ i • u = u i+1 -u i-1 2∆x . ( 7
)
Similarly the discrete gradient operator for the collocated surface-height can be expressed as

∇ i η = η i+1 -η i-1 2∆x . ( 8
)
The discrete Laplacian appearing in (6) results from the combination of these operators as

(∇ • ∇) i η ≡ ∇ i • (∇ i η) = η i+2 -2η i + η i-2 4∆x 2 . ( 9
)
The corresponding stencil is twice as wide as the initial discretisation and results in two discrete solutions for odd and even η i s only coupled through the boundary conditions. This is another way to express the existence of a computational mode for a collocated discretisation. Aside from the inherent instability of a scheme containing computational modes [START_REF] Lai | A projection method for combustion in the zero Mach number limit[END_REF][START_REF] Kar | Stable centered-difference schemes, based on an unstaggered a grid, that eliminate two-grid interval noise[END_REF], the wide stencil of the Laplacian greatly complicates the solution of ( 6) when a non-uniform resolution is used [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF]. The standard compact Laplacian

∇ 2 i η ≡ η i+1 -2η i + η i-1 2∆x 2
is also second-order and could be used instead in (6), however this would mean that the discrete Laplacian used is not anymore the convolution of the discrete divergence operator (7) and discrete gradient operator (8). The difference is O(∆x 2 ) and can be written

F i (η) ≡ (∇ • ∇) i η -∇ 2 i η = η i+2 -2η i+1 + 2η i -2η i-1 + η i-2 4∆x 2 .
Equation ( 6) then becomes

η n+1 i -θα ′ ∆t 2 ∇ 2 i η n+1 -θα ′ ∆t 2 F i (η n+1 ) = η n i -∆t∇ i • [(1 -θ)u n + θu ⋆ ].
Replacing the wide-stencil Laplacian with the compact-stencil Laplacian would thus be equivalent to adding a θα ′ ∆t 2 F i (η n+1 ) term to the left-hand-side of equation ( 6). This term can be seen as a spatial filter on η.

The initial wide-stencil Laplacian is constructed so that the discretisation of the continuity equation (6) guarantees volume conservation i.e. the discrete gradient of η n+1 applied to U ⋆ guarantees the discrete non-divergence (in the finite volume sense described in section 5.2) of U n+1 . This is the definition of an exact projection [START_REF] Peyret | Computational methods for fluid flow[END_REF]. Replacing the wide-stencil Laplacian with a compact stencil guarantees volume conservation only in an asymptotic manner (the error on the discrete divergence is of order ∆x 2 ). This defines an approximate projection [START_REF] Rider | Approximate projection methods for incompressible flow: Implementation, variants and robustness[END_REF][START_REF] Almgren | Approximate projection methods: Part I. Inviscid analysis[END_REF]). An important consequence is that the total energy is not exactly preserved anymore. This is clearly apparent in the form of the F i error which acts as a dissipative term on the free-surface height (and thus on the potential energy). In the case of incompressible flows, this dissipative term has been shown to have the beneficial side-effect of almost eliminating the pressure oscillations associated with a collocated discretisation [START_REF] Rider | Approximate projection methods for incompressible flow: Implementation, variants and robustness[END_REF]. Moreover, the relaxation of strict energy conservation has only a small effect on the quality of the solutions obtained and the transport of momentum and tracers can still be formulated as to guarantee exact conservation [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF]. In section 6.1 we study whether these conclusions still apply in the case of a linear free-surface ocean formulation.

Advection and conservation properties

The advection terms appearing in the momentum and tracer equations (contained in S n+1/2 in (2)) are computed using a conservative finite-volume predictor-corrector scheme. If we consider a tracer c contained in a cubic cell C of boundary ∂C the integrated advection term can be written

C A n+ 1 2 ≡ C [U • ∇c] n+ 1 2 + w n+ 1 2 ∂ z c n+ 1 2 = C ∇ • (U c) n+ 1 2 + ∂ z (wc) n+ 1 2 = (10) ∂C c n+ 1 2 (U n+ 1 2 • n + wn z )
where n = (n x , n y , n z ) is the unit normal to the boundary and we have used the Stokes theorem and the non-divergence (incompressibility) condition ∇ • U + ∂ z w = 0. A discretised form can be written as

C A n+ 1 2 i,j,k ≡ ∆ 2 [ (uc) i+ 1 2 ,j,k -(uc) i-1 2 ,j,k + (vc) i,j+ 1 2 ,k -(vc) i,j-1 2 ,k + (wc) i,j,k+ 1 2 -(wc) i,j,k-1 2 n+ 1 2 (11)
where ∆ is the size of the cell. Values of the tracer and normal components of the velocity are needed on cell faces at time level n + 1/2. We use a Godunov procedure where face values are extrapolated in space and time from time-and cell-centred values using a second-order Taylor series expansion. This scheme as been described in detail in [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] and [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF]) and we will just summarise its properties here. Upwinding is performed in order to reconcile values extrapolated on both sides of the faces and standard limiters can be used for fields with sharp discontinuities. The resulting scheme is the second-order, non-oscillatory, unsplit upwind scheme of Bell, Collela and Glaz [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF]. Stability is ensured independently of an explicit dissipation term, provided the CFL number is smaller than one.

It is clear from (11) that local conservation of c is guaranteed only if

u i+ 1 2 ,j,k -u i-1 2 ,j,k + v i,j+ 1 2 ,k -v i,j-1 2 ,k + w i,j
,k+ 1 2 -w i,j,k-1 2 = 0 which is the discrete expression of the non-divergence condition. The face velocities constructed using the Godunov scheme do not necessarily verify this condition. In our original scheme for incompressible flows an exact projection is applied to the face-based normal velocity field at time n + 1/2 to restore this condition. The resulting advection scheme is thus exactly conservative.

In the case of our hydrostatic, linear free-surface ocean model, the vertical component of the velocity w is obtained diagnostically through the vertical continuity condition (1) which is just another way to express the non-divergence condition. Noting that w appears only in the advection terms of the equations of motion and is thus needed only at cell faces, we can write a discretised finite-volume expression of (1) as

w i,j,k+ 1 2 = w i,j, 1 2 + k ′ =k k ′ =1 u i-1 2 ,j,k ′ -u i+ 1 2 ,j,k ′ + v i,j-1 2 ,k ′ -v i,j+ 1 2 ,k ′ (12)
where w i,j,1/2 is the boundary condition for w at the bottom of the ocean. The discrete non-divergence condition is thus trivially satisfied at all levels k and the advection scheme is locally conservative.

As noted in other studies [START_REF] Roullet | Salt conservation, free surface, and varying levels: A new formulation for ocean general circulation models[END_REF][START_REF] Griffies | Tracer conservation with an explicit free surface method for z-coordinate ocean models[END_REF], this does not guarantee global conservation however. This is due to the linear free-surface approximation which allows (small) fluctuations in the volume of each discretised vertical column of fluid. As a result, the vertical velocity at the top of the column obtained using ( 12) is non-zero (it is equal to the speed of the vertical free-surface motion). This generates either an outgoing tracer flux at the top of the column which is lost to the system (because the free-surface layer is not represented) or an extrapolated incoming flux which may not match the required boundary conditions. While these fluxes should be small when the linear approximation is valid, nothing guarantees that they balance over the whole system and strict global conservation is not ensured.

Adaptivity

The original incompressible flow solver Gerris uses a fully three-dimensional adaptive octree discretisation. Each dimension is treated in exactly the same way and refinement or coarsening occur simultaneously on all three dimensions (i.e. adaptivity is isotropic). A direct extension of this approach to oceanic flows is possible but poses several problems. Isotropic adaptivity means that the aspect ratio of cells is independent of the resolution. Taking a typical coarse oceanic discretisation with an horizontal resolution of 9.6 km and vertical resolution of 16 m, refining six times would lead to a cell with an horizontal resolution of 150 m and a vertical resolution of 25 cm. This is probably too high and would lead to a very large number of "layers" in the vertical. Non-isotropic refinement would allow the vertical resolution to be changed independently from the horizontal resolution. Such a capability would require substantial changes in the Gerris code base, however. For this first implementation we chose not to adapt the vertical resolution.

The model can then be seen as a fixed number of quadtree-discretised layers of constant thickness stacked vertically. The solution of the barotropic equation, hydrostatic pressure and diagnostic vertical velocity require local vertical integration over the whole depth of the water column. This is very simple when the discretised control volumes are aligned vertically i.e. when the horizontal resolution is independent of depth. We also add this constraint for our initial implementation. This is probably a more important restriction than the constant-thickness-layers approach because surface and deep currents are expected to have significantly different structures and associated spatial scales. Relaxing this constraint is certainly possible however and this will be the subject of future work.

When refining or coarsening the mesh, care must be taken to conservatively interpolate values from coarse to fine or fine to coarse cells. Simple averaging guarantees a conservative fine to coarse interpolation. For coarse to fine interpolation, values are interpolated linearly using the cell-centred gradient computed on the coarse cell. This also guarantees exact mass and momentum conservation.

Gerris provides a flexible framework for specifying adaptive refinement criteria. Several criteria can be specified simultaneously (such as vorticity, tracer gradients) and the threshold for refinement/coarsening as well as minimum and maximum number of refinement levels can themselves be functions of space and time. A limit on the total size of the simulation can also be set. When reached, the code automatically redistributes the cells in order to globally minimise the refinement cost given by the user-defined criteria.

We also note that Gerris does not use time-adaptivity in the sense that a single global timestep is used irrespective of the spatial resolution. Subcycling in time depending on the level of refinement is possible in principle but leads to complex synchronisation problems to ensure the discrete nondivergence of the global velocity field [START_REF] Almgren | A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations[END_REF]. In practice the finest grid cells often represent as much as half the total number of cells. Thus, an ideal subcycling implementation would only halve the computational cost compared to simple global time-stepping.

Representation of boundaries

Traditional height-coordinate ocean models on Cartesian grids usually represent the bottom shape of ocean basins using a "staircase" topography. The coastlines are also often restricted to following the edges of the square discretisation elements. This is a simple but rather crude approach which has been shown to substantially affect the accuracy of the solution [START_REF] Adcroft | Representation of topography by shaved cells in a height coordinate ocean model[END_REF][START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF].

A terrain-following vertical coordinate can be used instead of the height-coordinate [START_REF] Haidvogel | A semi-spectral primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear horizontal coordinates[END_REF][START_REF] Song | A semi-implicit ocean circulation model using a generalized topographyfollowing coordinate system[END_REF]. This theoretically ensures an accurate representation of the bottom topography (but does not improve the coastline representation). In practice the σ-coordinate transformation causes the surfaces of constant pressure (which are nearly horizontal) to intersect the vertical discretisation at large angles near steep topography. This causes large errors in the discretisation of the horizontal pressure gradients which lead to strong spurious currents near basin boundaries [START_REF] Mellor | The pressure gradient conundrum of sigma coordinate ocean models[END_REF][START_REF] Mellor | Sigma coordinate pressure gradient errors and the seamount problem[END_REF]. Substantial research has been done to minimise this problem but it is still a fundamental limitation of terrain-following models [START_REF] Shchepetkin | A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate[END_REF][START_REF] Ezer | Developments in terrain-following ocean models: intercomparisons of numerical aspects[END_REF][START_REF] Griffies | Developments in ocean climate modelling[END_REF].

Fully-unstructured meshes combined with a finite-element or finite-volume approach provide full flexibility for the discretisation of complex domains. This is the main reason for their widespread use in engineering problems. This flexibility comes at the expense of complexity and increased computational cost however. The sub-problem of mesh generation is itself far from trivial, especially in three dimensions [START_REF] Bern | Mesh Generation and Optimal Triangulation[END_REF][START_REF] Gorman | Optimisation based bathymetry approximation through constrained unstructured mesh adaptivity[END_REF]. We also note that most of the unstructured finite-volume or finite-element techniques recently developed for 3D geophysical fluid modelling [START_REF] Casulli | An unstructured grid, three-dimensional model based on the shallow water equations[END_REF][START_REF] Chen | An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[END_REF][START_REF] Bacon | A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA)[END_REF][START_REF] Fringer | An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator[END_REF]) can be seen as layered 2D discretisations. This solves the problem of coastline representation but does not improve the representation of the bathymetry (a standard terrain-following coordinate is used for the vertical discretisation). A notable exception is the fully three-dimensional model currently developed by Ford, Pain, Piggott, Goddard, de Oliveira, and Umpleby (2004) which uses unstructured tetrahedral meshes.

The quadtree discretisation does not solve the boundary representation problem, although it could be used to increase the resolution near a staircase representation of the topography and coastline. A staircase representation is not the only choice however. A variety of approaches have been developed which can combine the simplicity of a Cartesian discretisation with an accurate representation of complex boundaries: "embedded boundaries" [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Ye | An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries[END_REF], "cut-cell techniques" [START_REF] Quirk | An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies[END_REF][START_REF] Berger | A rotated difference scheme for Cartesian grids in complex geometries[END_REF][START_REF] Almgren | A Cartesian grid projection method for the incompressible Euler equations in complex geometries[END_REF] and for ocean models the "shaved cell" and "partial cell" approach of [START_REF] Adcroft | A new treatment of the Coriolis terms in C-grid models at both high and low resolutions[END_REF]. In our previous work [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] we have demonstrated how a "cut-cell technique" can be applied to an octree discretisation. The same technique is applied here and we will just give a general overview.

The technique relies on the finite-volume formulation of the system of equations. In a manner analogous to the derivation in (10) all the conservation laws can be recast as flux integrals over the boundaries of the discretisation elements. If the (cubic) cell is cut by a solid boundary the advection term in (11) can be written as the more general form

C A n+ 1 2 i,j,k ≡ a∆ 3 A n+ 1 2 i,j,k = ∆ 2 [ (suc) i+ 1 2 ,j,k -(suc) i-1 2 ,j,k + (svc) i,j+ 1 2 ,k -(svc) i,j-1 2 ,k + (swc) i,j,k+ 1 2 -(swc) i,j,k-1 2 n+ 1 2 (13)
where s i+1/2,j,k and a are the surface (resp. volume) fractions of the face (resp. cell) occupied by the fluid. All the discrete gradient, divergence and Laplacian operators can be generalised in a similar manner. This formulation is also used to compute consistent second-order-accurate fluxes between cells at different levels of refinement (see [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] for a detailed description). Second-order accuracy in space is maintained for the barotropic equation in complex domains. A "merged cell" technique is used to avoid the CFL restriction caused by very-small cut cells [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Quirk | An alternative to unstructured grids for computing gas dynamic flows around arbitrarily complex two-dimensional bodies[END_REF].

Numerical results

The semi-implicit scheme (θ = α = α ′ = 1/2) has been implemented as part of the Gerris flow solver. Most of the routines are common to both the ocean solver and the incompressible flow solver. In particular, the 2D multilevel Poisson algorithm is used with minor variations to solve the Helmholtz equation for the linear free-surface. The BCG advection scheme is also unchanged. The only routines strictly specific to the ocean model are the equation of state and the vertical integrations required for the hydrostatic pressure and vertical velocity.

In the following we present both two-dimensional test cases designed to highlight important properties and/or shortcomings of the model in the context of ocean modelling, as well as three-dimensional examples of realistic applications. The 2D simulations use the 3D model with a single layer in the vertical. We also try to give quantitative estimates rather than just a qualitative description of the model results.

Geostrophic adjustment problem

The approximate projection method has been shown to give good results for incompressible flows, however it has not been tested in the context of linearised inertia-gravity waves we consider here. The two main issues regard the stability of the computational mode associated with a collocated A-grid discretisation and the strength of the energy damping induced by the approximate projection.

As a first test case we consider the geostrophic adjustment problem studied by [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF]Le Roux, Staniforth, and[START_REF] Le Roux | Finite elements for shallow-water equation ocean models[END_REF]

. A Gaussian bump η(x, y) = η 0 e - x 2 +y 2 R 2
is initialised in a 1000×1000 km, 1000 m deep square basin. A reduced gravity g = 0.01 m/s is used to approximate a 10 m-thick stratified surface layer. The corresponding geostrophic velocities on an f -plane are initialised as

u(x, y) = 2gη 0 y f 0 R 2 e -x 2 +y 2 R 2 , v(x, y) = - 2gη 0 x f 0 R 2 e -x 2 +y 2 R 2
, where f 0 is the Coriolis parameter. Following Dupont we set f 0 = 1.0285 × 10 -4 s -1 , R = 100 km, η 0 = 599.5 m which gives a maximum geostrophic velocity of 0.5 m/s. In the context of the linearised shallow-water equations, the geostrophic balance is an exact solution which should be preserved by the numerical method. In practice, this would require an exact numerical balance between terms computed very differently: the pressure gradient and the Coriolis terms in the momentum equation. If this numerical balance is not exact, the numerical solution will adjust toward numerical equilibrium through the emission of gravity-wave noise which should not affect the stability of the solution. This problem is thus a good test of both the overall accuracy of the numerical scheme and its stability properties when dealing with inertia-gravity waves. We note in particular that a standard A-grid discretisation would develop a strong computational-mode instability in this case. Also, as studied by Le Roux, [START_REF] Le Roux | A semi-implicit semi-lagrangian finite-element shallowwater ocean model[END_REF], an inappropriate choice of finite-element basis functions will result in growing gravity-wave noise.

[Figure 3 about here.]

[Figure 4 about here.] Figures 3 and4 summarise the results obtained when running the geostrophic adjustment problem on a 64×64 uniform grid with a timestep ∆t = 1000 s. The maximum error on the height field (Figure 3) is small even after 18 days. After a strong initial transient corresponding to the emission of gravity waves, the error reaches a minimum at day 3 and then slowly grows with time with modulations due to the reflections of the initial gravity waves on the domain boundaries. As illustrated on figure 4, this growth is not due to any instability of the solution but to the slow decrease of the maximum amplitude of the Gaussian bump due to numerical energy dissipation.

Numerical energy dissipation is itself a problem and it is important to quantify it. As shown by [START_REF] Sadourny | The dynamics of finite-difference models of the shallow-water equations[END_REF] it is possible to formulate a strictly enstrophy-and energy-preserving finite-difference model on a C-grid. A finite-volume A-grid formulation using a semi-implicit time-discretisation and an exact projection would also be exactly energy-preserving [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF], however it would have a problematic computational mode as discussed previously. [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF] quantifies the energy conservation properties of various schemes using a non-linear variant of the geostrophic adjustment problem. A β-plane, f = f 0 + βy is used and the advection terms are included in the momentum equation. No explicit dissipation is added. As in [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF] we chose β = 1.607 × 10 -11 m -1 s -1 . The geostrophic eddy moves slowly westward through the emission of Rossby waves and southward due to the non-linear advection terms. The resulting evolution of the total energy is shown on figure 5. For our method, the slow decrease in the total energy is due both to the dissipation of potential energy induced by the approximate projection operator and to the dissipative properties of the BCG upwind advection scheme. Another run with the advection terms switched off (figure 5, square symbols) confirms that the dissipation induced by the approximate projection operator dominates the total dissipation. The results however compare favourably with the finite-element formulations tested by Dupont which all show significantly larger energy dissipation.

[Figure 5 about here.]

Coastally-trapped wave

To test the ability of the method to represent coastlines accurately, we reproduce the coastally-trapped wave test case of [START_REF] Curchitser | On the transient adjustment of a mid-latitude abyssal ocean basin with realistic geometry and bathymetry[END_REF]. The simulation is initialised with the exact solution of [START_REF] Lamb | Hydrodynamics[END_REF]. A constant-depth ocean on an f -plane is bounded by a circular coastline. The parameters are the same as used in IMCS ( 2006): h = 1500 m, g = 3.92 × 10 -2 ms -2 , f 0 = 8.34 × 10 -5 s -1 , radius of the ocean r = 600 km. The simulation runs for three wave periods (5.2465 days) with a timestep of 1199 seconds for all simulations. The linearised shallow-water equations (i.e. using a linearised free-surface and neglecting advection of momentum) are solved with no added explicit dissipation. The error on the final height field is estimated using the correlation technique of [START_REF] Curchitser | On the transient adjustment of a mid-latitude abyssal ocean basin with realistic geometry and bathymetry[END_REF]. The correlation function C is defined as:

C(θ 0 ) = M L(r, θ -θ 0 ) dx dy L 2 dx dy ,
where L(r, θ) is the analytical solution of Lamb and M (x, y) is the numerical solution. The phase error and the correlation error are then computed by finding the value θ max of θ 0 maximising C. Several runs were performed for different spatial resolutions (by varying the number of quadtree levels from 5 to 8) and using a constant or adaptive refinement strategy. The final solutions obtained with 7 levels of refinement for both a constant and an adaptive resolution are illustrated on Figure 6. The vorticity adaption criterion [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] was used to follow the waves as they travelled around the coastline.

Table 1 summarises the results for Gerris and reproduces the results obtained using FVCOM [START_REF] Chen | An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries[END_REF] and ROMS [START_REF] Shchepetkin | The regional oceanic modeling system (ROMS): A splitexplicit, free-surface, topography-following-coordinate oceanic model[END_REF]. FVCOM is a finite-volume solver using an unstructured triangular discretisation of the circular ocean. The constant resolution version of Gerris shows a convergence and accuracy which is better than FVCOM. The adaptive version shows a similar convergence and accuracy with a slight decrease of accuracy for highresolution due to accumulating interpolation noise when refining/coarsening the adaptive grid (this leads to slightly larger amplitudes close to the coastline which explains the correlation coefficients larger than one). ROMS is using a staircase representation of the coastline and consequently suffers from much larger phase errors than both Gerris (using the cut-cell technique) and FVCOM. The ROMS solution is also very noisy compared to the solutions illustrated in Figure 6. We also note that because both FVCOM and ROMS rely on an explicit barotropic solver, they used a timestep of around 250 seconds on the coarsest mesh (proportionally less on finer meshes).

An often-heard argument against the use of variable-resolution meshes for ocean modelling is the occurrence of spurious inertial and gravity waves reflections/scattering at refinement boundaries [START_REF] Griffies | Developments in ocean climate modelling[END_REF]. We believe one of the causes of spurious reflections is the use of discrete spatial operators which do not preserve second-order spatial accuracy when the resolution changes. The operators used in Gerris are designed to be second-order accurate everywhere [START_REF] Popinet | Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF] and we do not expect spurious reflections to occur at refinement boundaries. To test whether this is indeed the case, we repeated the coastallytrapped wave test case on a mesh where half the circular domain is resolved with double the resolution (9.375 km) of the other half (18.75 km). The results are given at the end of Table 1 and are close to those obtained on a mesh with a constant resolution equal to that of the coarser half of the domain (18.75 km), which confirms that no significant spurious reflections occur at refinement boundaries.

[Figure 6 about here.]

[Table 1 about here.]

Wind-driven circular ocean

Dupont presents another interesting test case for the representation of a circular coastline [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF]: the wind-driven circular ocean circulation. For a constant-depth ocean on an f -plane with a linear bottom friction, the linearised shallow-water equation admits the following stationary solution for the surface-height:

η(r, θ) = W f RgHκ R 2 8 + r 2 4 κ f sin 2θ -1 ,
where R is the radius of the ocean, g the acceleration of gravity, f the Coriolis parameter, H the ocean depth, κ the linear friction coefficient and W the gradient of the wind forcing:

τ x (r, θ) = W r sin θ R .
We follow Dupont and study the convergence with spatial resolution of the normalised error in η defined as:

E(η) = |η -η| dx dy dx dy dx dy η 2 dx dy ,
where η is the model solution. The model is run to convergence with f = 10 -4 s -1 , W = 10 -4 m 2 s -2 and κ = 10 -3 s -1 . Figure 7 illustrates the solution obtained with 5 levels of refinement.

The convergence of the normalised error with spatial resolution is plotted on Figure 8 together with the results obtained by Dupont for a streamfunction-vorticity C-grid model and the Hua-Thomasset finite-element model (HT) [START_REF] Hua | A noise-free finite element scheme for the two-layer shallow water equations[END_REF]. The HT model gives the best result of all the finite-element models tested by Dupont for this particular test case. An unstructured triangular mesh is used for the finite-element model.

[Figure 7 about here.]

[Figure 8 about here.]

The C-grid model shows much larger errors and only first-order convergence due to the stairstep representation of the circular coastline. Both the HT model and Gerris show second-order convergence. The errors obtained with Gerris are also substantially smaller than that obtained with the HT model.

Flow over a Gaussian bump

The test cases presented this far are only two-dimensional. An interesting three-dimensional configuration was used by Adcroft et al. [START_REF] Adcroft | A new treatment of the Coriolis terms in C-grid models at both high and low resolutions[END_REF][START_REF] Haidvogel | A semi-spectral primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear horizontal coordinates[END_REF] to validate their shaved cells and partial step representations of topography. A Gaussian bump is placed in a periodic channel bounded by solid walls with a slip boundary condition. The width of the channel is 320 km and its length 640 km. The 3D "primitive equations" of section 2 are solved with the parameters listed in Table 2.

[Table 2

about here.]

No explicit diffusion is added either on the density or momentum. The shape of the Gaussian bump is given by h

(x, y) = H -H 0 exp -(x -x 0 ) 2 + (y -y 0 ) 2 /L 2 ,
with x 0 = 200 km and y 0 = 150 km. The constant barotropic inflow u 0 is imposed by adding a spatially-constant body force in the x-direction. After 10 days the flow pattern illustrated on Figure 9 is obtained. The same contour intervals as in [START_REF] Adcroft | A new treatment of the Coriolis terms in C-grid models at both high and low resolutions[END_REF] are used for the representation. The result compares very well with Figure 14 and 18 of Adcroft et al. In particular the solution is as free of high-frequency noise as that obtained by [START_REF] Haidvogel | A semi-spectral primitive equation ocean circulation model using vertical sigma and orthogonal curvilinear horizontal coordinates[END_REF] with the spectral SPEM model (reproduced in Figure 14 of Adcroft et al.). The vorticity maximum in the detached vortex is well captured and does not show the excessive damping of the "full step" representation of topography (Figure 18(c) of Adcroft et al.). The Gerris solution is also cleaner than that obtained by Adcroft et al. with either the "piecewise linear" or "partial step" methods.

[Figure 9 about here.]

Tidally-driven flow in Cook Strait

We conclude with an example of adaptive barotropic flow with a complex coastline and bathymetry and a large ratio between the smallest and largest scales. Cook Strait, the stretch of water between the North and South Islands of New Zealand is a complex, highly energetic marine environment. The Marlborough Sounds -a drowned mountain river system at the northern end of the South Island -have a very convoluted coastline and provide a good test of the robustness of the treatment of embedded solid boundaries.

The simulation domain is illustrated on Figure 10. The eastern part of the strait is relatively shallow while the western part drops off steeply. Because we only solve for the barotropic flow in this example, we limit the maximum depth to 400 meters as an approximation to a fixed-depth thermocline. A "Flather" boundary condition is used on all four boundaries [START_REF] Flather | Tidal computations for Morecambe bay[END_REF] where only the M 2 component of the surface elevation is imposed, using results from the Walters and Goring tidal model [START_REF] Walters | Ocean tides around New Zealand[END_REF]. The main feature of the M 2 tides around New Zealand is a single Kelvin-like wave propagating anti-clockwise around both islands. As a result, the tides on the eastern and western sides of Cook Strait are in opposition of phase with an amphidrome close to the narrowest point of the passage. The resulting tidal pressure gradient then forces strong currents both in the main passage and in some branches of the Marlborough Sounds.

[Figure 10 about here.]

The full 3D model is used to solve for the barotropic part of the solution by using only a single vertical layer. The bathymetry is described using the full 3D volume-of-fluid technique described in section 5.4. A small amount (4×10 -6 s -1 ) of linear bottom friction is added to dampen large velocities in very shallow areas. The resolution of the cells defining the coastline is fixed at 12 quadtree-levels i.e. 500 km / 2 12 ≈ 122 meters. The minimum resolution is set to 6 levels i.e.≈ 7.8 km. Away from the coastline the resolution is adapted dynamically at every timestep according to the local vorticity. The timestep is also adapted dynamically in order to verify the CFL condition required for the stability of the advection scheme. The average timestep is approximately 40 seconds. In established regime the total number of cells is approximately constant at ≈ 140, 000 (compared to (2 12 ) 2 ≈ 17, 000, 000 for a regular Cartesian grid of equivalent resolution). Each timestep takes 10 seconds of CPU time on a Pentium 4 2.6 GHz processor, of which less than 5% is used to adapt the grid. The simulation required approximately 50 MBytes of memory.

Figure 11 and 12 shows part of the depth-averaged vorticity field for t = 6 days. The combination of strong currents and sharp headlands creates vorticity "dipoles" close to the coastline which are then advected away from the coast by the mean tidal circulation. The weak dissipation due to bottom friction in deeper water as well as the good energy conservation properties of the numerical scheme lead to long-lived eddies which interact with newer eddies produced at each tidal cycle. This leads to the complex vorticity field illustrated in the figures. This process is more clearly illustrated by the accompanying animation [START_REF] Popinet | Animation of the vorticity field in Cook Strait[END_REF].

[Figure 11 about here.] Figure 12 displays a detail of the quadtree mesh for t = 6 days. High-vorticity areas are tracked with a maximum resolution of 12 levels. Using a statically refined mesh with a high resolution near the coastline decreasing offshore leads to a very different picture. Eddies are strongly damped by numerical diffusion as soon as they leave the coastline and the process of inter-tidal interactions between eddies does not occur.

While the code has been validated in the Cook Strait area for the prediction of tidal amplitudes and currents [START_REF] Msadek | Hydrodynamic tidal model of Cook Strait[END_REF], the high-resolution results presented here still need to be compared with available observational measurements. However, satellite imagery provide convincing first evidence that the processes of eddy generation off headlands, subsequent vorticity transport offshore and intertidal eddy interactions illustrated here are important in reality. The International Space Station pictures in the accompanying material (NASA, 2002a,b) show structures very similar to the small offshore eddies and vorticity layers off headlands illustrated in Figure 11.

One may argue that trying to resolve individual eddies at the scales illustrated here is not necessary as these scales are meant to be adequately represented by sub-grid turbulence models, however most turbulence models rely on assumptions (e.g. spatial isotropy) which are not verified in strongly sheared flows such as boundary layers. We believe that the high resolution made possible by adaptive techniques will be a key requirement in accurately modelling the on-shore/off-shore transport processes controlling sediment dynamics, nutrient fluxes in the marine ecosystem and bio-geochemical cycles generally.

[Figure 12

Conclusion

The development of an adaptive ocean model has been described. The dynamic adaptivity -allowing for multi-scale representation within a single model that can evolve in time -is achieved through a conservative finite-volume discretisation which requires collocation of all the geophysical variables (A-grid discretisation). It is demonstrated, however, that the extension to a free-surface ocean model of the approximate projection method overcomes the A-grid limitations at the expense of a small amount of numerical energy dissipation. Test cases clearly underline the accuracy and efficiency of this adaptive ocean model; accurate geostrophic balance, robust and accurate representation of bathymetry and coastlines are demonstrated.

As a first demonstration of the functional adaptivity, the present three-dimensional ocean model has a linear free surface, and is incompressible, hydrostatic, and Boussinesq. Here, the adaptivity is also limited to the horizontal plane only, and is projected down through a series of z-level, stacked layers in the vertical. Nevertheless, the boundaries are represented with a "cut-cell" technique, so the familiar "staircase" pitfalls of typical finite difference, z-level models are avoided [START_REF] Adcroft | Representation of topography by shaved cells in a height coordinate ocean model[END_REF]. Furthermore, the range of relevant geophysical solutions for such an ocean model remains large, and therefore permit the exploration of a number of testing problems.

For the case study of coupling between a complex coastal and offshore barotropic flow, we demonstrate how the adaptive mesh results in a reduction of two orders of magnitude in the number of grid points. Moreover, the dynamic adaptivity gives the model an advantage over other models that have an unstructured spatial grid that is fixed in time, as it has the ability to track evolving fine-scale features (e.g. eddies shed from headlands, meandering frontal boundaries etc...). This is particularly important for the accurate simulation of processes which are not adequately described by current subgrid-scale models (e.g. turbulence near complex coastlines). Compared to more classical unstructured grid approaches, the "semi-structured" quadtree grid we employ is very easy to generate and -as a result -the cost of adapting the grid at every timestep is a small fraction of the total computational cost. While time-adaptivity is certainly possible on unstructured grids, it is expensive and cannot usually be applied at every timestep [START_REF] Bacon | A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA)[END_REF]. In the context of geophysical fluid flows, we are only aware of two studies demonstrating time-adaptive unstructured grids [START_REF] Bacon | A dynamically adapting weather and dispersion model: The operational multiscale environment model with grid adaptivity (OMEGA)[END_REF][START_REF] Ford | A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows. part I: Model formulation[END_REF].

From a user perspective, the fact that adaptive mesh generation is part of the solution procedure allows flexible control of the level of refinement, so that issues of convergence -for example -can be addressed by simply changing an input file, thereby avoiding having to do a model re-build each time the resolution is altered. The hierarchical quadtree structure is also naturally suited to multigrid algorithms which leads to robust and efficient solvers for the implicit free-surface discretisation even with arbitrarily complex coastlines.

We are currently working on generalising this approach to a non-hydrostatic, nonlinear free-surface formulation with adaptivity in all dimensions. The underlying adaptive Navier-Stokes fluid solver Gerris naturally deals with non-hydrostatic, nonlinear free surface problems, in which adaptivity in all three spatial dimensions is routinely performed. In the oceanic case, decoupling of the hydrostatic and non-hydrostatic pressures can be employed to overcome the stiffness of the pressure equation due to the large aspect ratio between horizontal and vertical scales (Marshall, Adcroft, Hill, Perelman, and Heisey, 1997a). The extension of adaptivity to the vertical coordinate could be more problematic and has not yet been studied in any significant detail by us or others. The main issue will be ensuring the accurate discretisation of the horizontal hydrostatic pressure gradients across levels of refinement. Another important general issue for adaptive solvers which has received little attention is the development of criteria for refinement to control errors in subgrid-scale parametrisation.

We are confident that these issues can be overcome and that the resulting dynamic adaptivity will be an essential component of the next generation of geophysical fluid flows models.

A Stability of Adams-Bashforth discretisation

Model simulations of the geostrophic adjustment problem tested an Adams-Bashforth scheme to discretise the Coriolis terms in time. The simulations seemed to require a timestep much shorter than expected as the simulation proceeded. This appendix confirms the need for such a short timestep, and motivates the use of the semi-implicit scheme for the Coriolis terms which is possible within the adaptive model because of the collocated velocity components.

A simple, one-dimensional shallow water system is used to explore the stability properties of an inertia-gravity wave system using Adams-Bashforth discretisation for the Coriolis terms, viz,

∂u ∂t = f v -g ∂η ∂x , (14)
∂v ∂t = -f u, (15)
∂η ∂t = -H ∂u ∂x . ( 16
) Taking u = ûcoskx, v = v cos kx, η = ηsinkx, analytic solutions are, u = u o cosωt cos kx, v = -u o f ω sinωt coskx, η = Hku o f ω sinωt sinkx, (17)
for each wavenumber k with frequency ω, where,

ω 2 = f 2 (1 + k 2 R 2 ), R 2 = gH f 2 , (18)
and u o is a constant. Following Marshall et al. (1997), an Adams-Bashforth method (AB2) is used to discretise the Coriolis terms, while the remaining terms are allowed to range from fully explicit to fully implicit. The discrete (in time) system can then be written,

u n+1 -u n = a(v n + α(v n -v n-1 )) -b(dη n+1 + eη n ), (19)
v n+1 -v n = -a(u n + α(u n -u n-1 )), (20)
η n+1 -η n = c(du n+1 + eu n ), (21)
where a = f ∆t, b = gk∆t, c = Hk∆t, α = 1/2 + χ, d = θ, and e = 1 -θ, where θ is the implicitness parameter such that 0 ≤ θ ≤ 1. As in Marshall et al. (1997), the offset in the AB2 scheme χ is set to be 0.1, and is consistent with the stability of this AB2 scheme alone, i.e., a ≤ 1/2 when χ > a 2 /4 + a 4 /2. In terms of an amplification factor λ such that η n+1 = λ η n , the discrete system can be reduced to the following fourth order polynomial;

λ 4 (1 + cbd 2 ) + λ 3 (2bcde -2) + λ 2 (1 + bce 2 + a 2 (1 + α) 2 ) + λ(-2a 2 α(1 + α)) + a 2 α 2 = 0. ( 22
)
Stable numerical solutions arise for |λ| < 1. The fourth order system reduces to two complex conjugate pairs, one representing the computational mode which is found to be always stable for the range of parameters explored here, the other the physical mode. Solutions for the physical mode are plotted in Figure 14 as a function of the dimensionless numbers f ∆t and Rk. The region of stability is that below (or to the right of) each curve, and the different curves are labelled with the value of the implicitness parameter θ. The horizontal dashed line is the AB2 stability boundary in the absence of gravity waves. For pure gravity waves with θ ≥ 0.5 the system is always stable.

The important result is that the stability region for the full inertia-gravity wave system falls well below the AB2 line at f ∆t = 1/2 for the centred-in-time θ = 1/2 discrete system, particularly as Rk increases. Such Rk changes represent either an increase in the depth H of the domain, or an increase in the wavenumber k (or equivalently a reduction in the resolved wavelengths). This theoretical solution explains the need for significant reductions in the timestep as the generation of resolved, fine-scale structure in the numerical model solution occurs. Even though the separate gravity and inertia waves are theoretically stably resolved by f ∆t = 1/2, their interaction as an inertia-gravity wave system significantly brings down the available region of stability.

Allowing for more implicitness in the gravity wave terms via increasing θ expands the stability range, and with θ = 0.6 a permissible timestep around f ∆t = 1/2 is recovered for all values of Rk. Further increases in θ allow for significant increases in f ∆t at larger values of Rk, with smaller gains for smaller Rk. However this comes at the cost of increased dissipation of gravity waves.

[Figure 14 [START_REF] Sadourny | The dynamics of finite-difference models of the shallow-water equations[END_REF] and implemented by [START_REF] Dupont | Comparison of numerical methods for modelling ocean circulation in basins with irregular coasts[END_REF]. The finite-element formulations are those studied by Dupont. LW: [START_REF] Lynch | Three-dimensional hydrodynamics on finite elements. Part I: Linearized harmonic model[END_REF], LLS: Le Roux, [START_REF] Le Roux | A semi-implicit semi-lagrangian finite-element shallowwater ocean model[END_REF], PZM: [START_REF] Peraire | Shallow water problems: A general explicit formulation[END_REF].

  about here.] [Figure 13 about here.]

Figure 1 :Figure 2 :

 12 Figure 1: Example of quadtree discretisation and corresponding tree representation. The dotted lines on the right-hand-side define the level of the cells.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Evolution of the maximum error on the surface height for the geostrophic adjustment problem.

Figure 6 :

 6 Figure 6: Surface-height isolines for the coastally-trapped wave solution at t = 5.2465 days (3 wave periods). (a) Constant resolution 9.375 km. (b) Adaptive resolution, white 9.375 km, black 150 km.

Figure 7 :

 7 Figure 7: Velocity and height field for the wind-driven circular ocean test case.

Figure 8 :

 8 Figure 8: Convergence of the normalised error with spatial resolution for different models. Both the HT and C-grid results are reproduced from the study of Dupont (2001).

Figure 9 :

 9 Figure 9: Non-dimensional depth-integrated vorticity ζ/f at t = 10 days. Contour interval = 100.

Figure 10 :

 10 Figure 10: Coastline and isobaths for the Cook Strait tidal model. The domain extent is 500 × 500 km. Both the coastline and the isobaths are drawn using their exact representation on the quadtree mesh. The light grey area indicates depths deeper than 400 metres. The black box indicates the area zoomed in on Figure 11.

Figure 11 :

 11 Figure 11: Detail of the vorticity field at t = 6 days. The domain extent is 83 × 73 km. The black boxes indicate the areas represented on Figure 12 and 13.

Figure 12 :Figure 14 :

 1214 Figure 12: Detail of the adaptive quadtree mesh. The coastline is statically refined at the finest level while individual eddies are tracked dynamically. The domain extent is 36 × 29 km. The minimum spatial resolution is 122 metres.

Table 1 :

 1 Results for the coastally-trapped waves in a circular basin. The FVCOM and ROMS results are reproduced from IMCS (2006).

	• )

Table 2 :

 2 Parameters for the Gaussian bump test case.

List of Figures

List of Tables