## Abstract

Research into the use of unstructured mesh methods in oceanography has been growing steadily over the past decade. The advantages of this approach for domain representation and non-uniform resolution are clear. However, a number of issues remain, in particular those related to the computational cost of models produced using unstructured mesh methods compared with their structured mesh counterparts. Mesh adaptivity represents an important means to improve the competitiveness of unstructured mesh models, where high resolution is only used when and where necessary. In this paper, an optimization-based approach to mesh adaptivity is described where emphasis is placed on capturing anisotropic solution characteristics. Comparisons are made between the results obtained with uniform isotropic resolution, isotropic adaptive resolution and fully anisotropic adaptive resolution.

## 1. Introduction

Unstructured meshes have become popular in the engineering literature for problems that involve complex geometries. Even the most complex engineering problem would struggle to match the complexity of the geometry one must contend with when simulating the oceans; i.e. the fractal nature of coastlines and bathymetry. Details of both the large-scale meridional overturning circulation and important smaller scale fluid dynamical processes are heavily influenced by the geometry of ocean basins (Siedler *et al.* 2001). Hence, when simulating this system numerically, it is highly desirable to seek an optimal representation of the physical domain. Unstructured meshes represent an attractive means to achieve this. Hence, their use in oceanography represents an interesting alternative to structured mesh methods, which have undergone many decades of development and optimization; for background, see Griffies *et al.* (2000), Haidvogal & Beckmann (1999) and Kantha & Clayson (2000). For further discussion on the use of unstructured meshes to represent complex ocean domains, see Gorman *et al.* (2006, 2007), Lambrechts *et al.* (2008) and Legrand *et al.* (2000).

With an unstructured mesh, it is straightforward to use variable resolution within the domain. For example, one could place additional resolution where complex dynamics are occurring, e.g. boundary layers, or where important small-scale geometry features are present, e.g. narrow straits, if those locations are known in advance.

However, discretization methods based upon unstructured meshes are typically more costly than those based on structured meshes, measured per degree of freedom (Löhner 2001). As one is often very concerned with run times, the additional flexibility afforded by unstructured meshes must be exploited. The aim would be that unstructured mesh methods not only give better, more faithful, discrete representation of the ocean dynamics, but also that there is either no or actually improved cost implications in the resulting models.

The numerical methods employed to discretize the underlying equations on unstructured meshes typically fall into the finite-element/-volume categories. These methods can yield significant flexibility in the representation of solution variables. This flexibility can be used to impart desirable discrete characteristics and stability/accuracy properties: for example, the ability to accurately simulate geotrophically balanced states and wave propagation phenomena (Le Roux *et al.* 1998; Cotter *et al.* 2009). Examples of some of the early work into finite-element/-volume methods in oceanography can be found in Dumas *et al.* (1982), Fix (1975) and Haidvogel *et al.* (1980). For more recent examples, see Chen *et al.* (2003), Danilov *et al.* (2004, 2005), Fringer *et al.* (2006), Pain *et al.* (2005), Piggott *et al.* (2008*a*) and White *et al.* (2008).

As the dynamics of the oceans are highly unsteady, producing a static optimal variable-resolution unstructured mesh suitable for a simulation over a long time period is virtually impossible. One must therefore consider the use of dynamically adaptive methods, where the resolution is varied in time in response to the evolving solution fields. For example, when a front or eddy moves from one position to another, the mesh resolution in some sense should follow or pre-empt this. Equally, if an important process requiring enhanced resolution occurs sporadically in time, e.g. a seasonal effect, the mesh should be able to focus on the process as it forms, and then de-focus as the process dies away. For some examples of adaptive mesh methods in ocean modelling, see Bernard *et al.* (2007), Piggott *et al.* (2005, 2008*b*) and Popinet & Rickard (2007). See also Blayo & Debreu (1999) and Debreu *et al.* (2005) for the use of adaptive methods in structured mesh models.

An additional point to note is that as well as the ocean domain having a high aspect ratio, the dynamics of the ocean circulation can also exhibit high aspect ratios. Oceanic phenomena, which can possess strong anisotropies, or directional variations, include boundary layers and currents, fronts, jets and pycnoclines. There are therefore many regions of the ocean where anisotropic, in addition to simply variable, mesh resolution may be required for optimality.

The purpose of this article is to explore, through example problems, the use of an optimization approach to anisotropic unstructured mesh adaptivity. In the next section, the mesh optimization approach is presented. This includes elemental operations used to improve the shape and size profiles of elements making up the mesh, functionals used to gauge the quality of an element or a mesh and the construction of an error metric used to differentiate regions where increased/decreased resolution is required. Following this, three barotropic gyre-based examples are presented to demonstrate the advantages that anisotropic variable resolution may deliver over the use of fixed uniform isotropic resolution. The paper concludes with a summary of the findings of this work and discussions on some of the extensions required for the use of mesh adaptivity on more complex real-world problems.

## 2. Optimization-based mesh adaptivity

### (a) Mesh optimization operations

Given an unstructured mesh and information regarding the ideal shape and sizes of the elements making up the mesh, an optimization-based adaptivity algorithm can be formulated via the use of local topological operations that seeks to improve the quality of elements.

In the examples presented in this work, a two-dimensional mesh optimization algorithm (Agouzal *et al.* 1999; Vasilevskii & Lipnikov 1999) is used that employs the following local operations depicted in figure 1.

(i)

*Node insertion or edge split*. Here a node is inserted on a pre-existing edge in the mesh so that the four new elements have improved shape/size characteristics compared with the original two; while the location of this new node along the pre-existing edge can be optimized, it is common to simply split it at its midpoint.(ii)

*Node deletion or edge collapse*. Here the inverse operation is performed whereby an edge in the mesh is collapsed, and consequently a node is deleted and two elements removed from the mesh.(iii)

*Edge swap*. Here an edge between two elements is removed and replaced with the only other possible configuration in two dimensions; the number of nodes and elements is preserved through the operation, but the edge lengths and element shapes are manipulated.(iv)

*Node movement*. Here a node is moved so as to improve the quality of all the elements surrounding it; the direction in which to move the node can be calculated by considering a discrete gradient of the mesh quality functional and performing a line search in this direction.

### (b) The objective functional

In order to decide which of the operations described above should be used to improve the quality of the mesh, an objective functional measuring the quality of the mesh must be defined. There are multiple choices for such a mesh quality functional. When one considers anisotropic mesh optimization, it is common to form a functional that depends on both the shape and size of the elements making up the mesh. In the work presented here, the two-dimensional objective functional of Vasilevskii & Lipnikov (1999) is used; see also Agouzal *et al.* (1999), Buscaglia & Dari (1997), Knupp (2003) and Pain *et al.* (2001). For an element Δ in the mesh, this is defined as
2.1
where |Δ| is the area of element Δ and |∂Δ| is its perimeter. The first factor (I) is used to control the shape of element Δ. For an equilateral triangle with sides of length *l*, and |∂Δ|=3*l* and so I=1. For non-equilateral triangles, I<1. The second factor (II) controls the size of element Δ. The function *F* is smooth and defined as , which has a single maximum of unity with *x*=1 and decreases smoothly away from this with . Therefore, II=1 when the sum of the lengths of the edges of Δ is equal to 3, e.g. an equilateral triangle with sides of unit length, and II<1 otherwise. Hence, taken together, the two factors making up *Q* yield a maximum value of unity for an equilateral triangle with edges of unit length, and *Q* is less than one otherwise.

The quality functional for the entire mesh can then be taken to be the minimum value of the elemental quality functional taken over all the elements in the mesh.

In order for the quality functional, and hence the mesh optimization procedure, to yield variable-resolution anisotropic meshes, the element quantities |Δ| and |∂Δ| are actually measured with respect to the non-Euclidean metric *M*, which varies in space, hence the presence of the subscript *M* in equation (2.1). An example of an edge swap operation improving the quality of two triangles is shown in figure 2.

### (c) Error measures and metrics

In the previous two subsections, mesh optimization operations and an objective functional have been introduced. These can be used together to formulate an optimization approach to deliver optimally shaped elements within a given domain. The element functional (2.1) depends upon the area and edge lengths of the element, and these scalars were defined in terms of a metric *M*. Assuming that this is the standard Euclidean metric, then this procedure would seek to deliver uniformly sized equilateral triangles in the domain. To turn this procedure into one that delivers a variable anisotropic resolution, a non-Euclidean metric *M*(** x**) can be used. Given element Δ, one may assume, both for simplicity of explanation and also for numerical implementation, that

*M*is constant over Δ. In this case, where |Δ|

_{E}is the area of element Δ in the standard Euclidean metric,

*e*_{i}is the vector representing the

*i*th edge of element Δ and ||∂

*e*_{i}||

_{M}is its length with respect to the metric

*M*.

The choice of *M* is critical in achieving an adapted mesh that in some sense is optimal for given solution characteristics. Interpolation error estimates form an excellent starting point for the construction of a suitable metric *M*. In one dimension, the linear interpolation error of a smooth function is bounded by the product of the second derivative of the function being interpolated and the size of the region being considered. For example, for a finite-element approximation, this would be the size of the element over which the piecewise-linear interpolant is defined. In higher dimensions, a similar result holds, with the second derivative replaced by the Hessian of the function being interpolated, the tensor field of second partial derivatives (D’Azevedo 1991; George 1998; Frey & Alauzet 2005). The Hessian describes the curvature of the solution at each point in every direction; using the Hessian as a basis for the metric *M* results in enhanced resolution where the solution field has high curvature. For elliptic problems, Céa’s lemma (Brenner & Scott 2002) shows that the interpolation error bounds the discretization error in the energy norm; in practice, minimizing the interpolation error in this manner has been found to be highly effective in controlling the discretization error, even for non-elliptic problems (e.g. Frey & Alauzet 2005).

For a mesh of triangles in two dimensions, assume a sufficiently smooth function *u*≡*u*(** x**) is approximated by its piecewise-linear Lagrange interpolant Π

_{h}

*u*on element Δ. The interpolation error on element Δ satisfies 2.2 where

*c*is a constant independent of the mesh,

**is a vector in Δ,**

*v**E*

_{Δ}is the set of edges of Δ, is the infinity norm over Δ and is an element-valued Hessian (Pain

*et al.*2001; Frey & Alauzet 2005). |

*H*| denotes the positive-definite metric formed by taking the absolute value of the eigenvalues of

*H*and reflects the fact that it is the magnitude of the curvature that is of interest, rather than the sign of the curvature.

For simplicity, the constant *c* is assumed to take the value 1. Then, for a given interpolation error, *ϵ*_{u}, the metric *M* is defined such that an element edge is of unit length with respect to *M* if it has the desired interpolation error, i.e.
2.3
The user-defined interpolation error, *ϵ*_{u}, may vary spatially and temporally to yield additional control over the adapted mesh.

Construction of an adapted mesh that equidistributes the interpolation error throughout the domain can be seen from equation (2.2) to be equivalent to the construction of a uniform mesh of equilateral triangles, i.e. all element edge lengths are unity with respect to *M*. It is the functional (2.1) that defines how close each element is to this requirement of equidistribution. For other approaches to the construction of error metrics, see Micheletti & Perotto (2005), Power *et al.* (2006) and Bernard *et al.* (2007).

A symmetric positive-definite metric tensor *M* may be decomposed as
where *Λ* is the diagonal matrix consisting of the eigenvalues of *M* (the principal curvatures of the solution field for *M*=*H*) and *Q* is an ortho-normal matrix of eigenvectors *Q*^{i}. *Q* may be interpreted geometrically as a rotation and *Λ* a rescaling. The eigenvalue *λ*^{i} encodes the desired edge length in the direction *Q*^{i} as
i.e. anisotropic as well as inhomogeneous resolution results from a mesh that respects the metric *M*. Note that the transformation may be used to transform a uniform isotropic element in a computational space into an adapted anisotropic element in the physical domain, achieving the desired interpolation error everywhere (figure 3). This can be seen from the decomposition of *M* above and

### (d) Combining multiple fields and bounding aspect ratio and edge lengths

For complex problems, one will wish to simultaneously adapt the mesh for several solution fields. The approach taken here is to use a single mesh to represent all fields and hence typically a compromise needs to be made. Often this compromise will not be substantial because in many problems, regions and directions of high curvature are common between solution fields. For oceanographic problems, the base set of fields we would want to adapt to are the components of velocity, temperature and salinity. Metric fields are constructed for each with their own interpolation errors. A superposition procedure is then undertaken so that a single metric results, which respects the edge length and requirements of each individual metric. The procedure follows Pain *et al.* (2001); see also Castro-Diaz *et al.* (1997) and George (1998).

In addition, it is sometimes advantageous to modify the metric to impose maximum and minimum element sizes on the mesh: for example, to stop the optimization procedure from attempting to produce elements larger than the solution domain, or smaller than a certain scale which the user deems to be unimportant or which they wish to use subgrid-scale models to represent. Also, for problems that have known high-aspect ratio dynamics or domains, it can be necessary to impose these maximum and minimum constraints directionally. The approach taken to achieve this here is to construct extra metrics describing maximum and minimum edge lengths. For example, assuming the user wished to limit the maximum element sizes in the *x* and *y* directions to be *h*_{x,max} and *h*_{y,max}, respectively, then the metric *M*_{max}=*diag*((*h*_{x,max})^{−1/2},(*h*_{y,max})^{−1/2}) could be superimposed. For the minimum edge-length metric, the operator must be replaced with a operator. More general directions for the maximum and minimum edge lengths can be introduced through the use of the rotation matrix *Q*.

To bound the aspect ratio of elements in physical space, a limit on the ratio of eigenvalues may also be performed. For example, to limit the aspect ratio of elements to be below *a*, the eigenvalues of the metric can be modified as follows:
i.e. this operation seeks to limit the minimum eigenvalue, and hence the maximum edge length in an element. It thus achieves a maximum bound on aspect ratio without compromising the interpolation error achieved over that element.

### (e) Hessian recovery

In practice, only the discrete approximation to solution fields are available, not the continuous form of the function being interpolated. Thus, it is necessary to recover a discrete approximation to the true Hessian. Different recovery schemes such as quadratic fitting (Vallet *et al.* 2007), superconvergent patch recovery (Zienkiewicz & Zhu 1992), integration by parts (Buscaglia & Dari 1997) and double Galerkin projection (Pain *et al.* 2001) have been suggested. These schemes vary in cost, accuracy and robustness; a number of comparisons of Hessian recovery schemes have recently been published (e.g. Lipnikov & Vassilevski 2006; Vallet *et al.* 2007). Here, the double Galerkin projection approach described in Pain *et al.* (2001) is used.

## 3. Examples

### (a) The Stommel gyre

#### (i) Problem description

This test problem involves the steady-state wind-driven barotropic circulation in a rectangular domain, and compares against an analytical solution (Stommel 1948). The problem is oceanographically relevant, and numerically challenging, owing to the presence of a thin boundary current on the western side of the domain. Following Stommel (1948) and Hecht *et al.* (2000), consider the following equation for the unknown streamfunction, *Ψ*, in the domain [0,*L*]^{2}, *L*=1:
3.1
where *β*=50 is the north–south derivative of the assumed linear Coriolis parameter, *F*=0.1 is the strength of the wind forcing which takes the form and *γ*=1 is the strength of the assumed linear frictional force. Homogeneous Dirichlet conditions are applied on all boundaries. The analytical solution to equation (3.1) is given by
3.2

#### (ii) Discretization

To solve this problem numerically, the ability of the adaptive mesh ocean model employed in this work to solve general advection–diffusion equations is utilized. Consider the equation
3.3
with ** a** a transport velocity,

*κ*a tensor of diffusivities and

*f*a source term. This is consistent with equation (3.2) with the choice

*u*=

*Ψ*,

**=−(**

*a**α*,0)

^{T},

*κ*=

*I*

_{2}and .

A standard Galerkin finite-element discretization in space with continuous piecewise-linear basis function is used. Backward Euler differencing is used in time and the calculation is conducted until a steady state is reached; here integrating until *t*=1 with a time step of Δ*t*=0.0025 is sufficient. For the adaptive simulations, the mesh is adapted every 40 time steps. Between adapts, the solution field is interpolated using the simple linear interpolation implied by the underlying basis functions. For a discussion of more advanced interpolation options, see Farrell *et al.* (2009).

#### (iii) Results

Figure 4 shows a plot of the computed streamfunction for this problem along with an adapted mesh. Figure 5 shows a blow-up of the region close to the western boundary where anisotropic resolution results when no limit is placed on the aspect ratio of elements. When a limit on aspect ratio is introduced, one can see that this results in significant differences in the mesh close to the boundary and an overall loss of efficiency in the calculation.

Figure 6 plots the error between the computed and analytical solutions on different meshes. The analytical solution has been calculated on the nodes of the mesh and then quadrature over elements used to evaluate an *L*^{2} error norm. In all cases, the error can be seen to decrease linearly with the number of nodes in the mesh, which is consistent with second-order convergence for this two-dimensional problem. For uniform isotropic meshes, the error is larger than for adapted meshes, which are able to place the increased resolution near the western boundary where the streamfunction experiences rapid variation and higher numerical errors. Also shown is the impact that limiting the maximum allowed element aspect ratio has on the error–nodecount profile.

For example, in this problem, suppose one wishes to achieve an error of 10^{−6}, then from figure 6 one can see that with isotropic adaptivity (e.g. with a limit of two on aspect ratio) one would need to make use of approximately a factor of eight fewer nodes than would be required in a uniform isotropic unstructured mesh. With fully anisotropic adaptivity, a further saving of approximately a factor of five can be achieved.

Alternatively, suppose one had a machine that allowed the use of 10^{4} degrees of freedom. Then with anisotropic adaptivity, we can achieve an error approximately five times lower than isotropic adaptivity, which itself can yield an error approximately a factor of eight lower than uniform refinement for this problem.

Hence, for this problem, the advantages of adaptive resolution, and in particular anisotropic resolution, are clear. Also, these advantages increase as the thickness of the boundary layer decreases.

### (b) Advection in a Stommel gyre

#### (i) Problem description

In this problem, the advection of a passive tracer by a prescribed velocity field is considered. The problem is taken from Hanert *et al.* (2004); see also Hecht *et al.* (2000). The problem domain is an *L*×*L* box with *L*=1000 *km*. The prescribed velocity field is obtained from the Stommel streamfunction
where the magnitude of the wind stress is *F*=0.1 *N* m^{−2}, the frictional coefficient is *γ*=10^{−6} s^{−1}, the density is *ρ*=1000 *kg* m^{−3}, the fluid depth is *H*=200 m and *α*=*β*/*γ* is chosen such that *α*^{−1}=100 km.

The initial condition for the passive tracer is given by the Gaussian function 3.4

The reference solution at time *t* used to compare numerical results against is constructed as follows. The ordinary differential equation (ODE) system is integrated repeatedly from time *t* to 0, with initial condition given by the node positions on the mesh at time *t*. This yields the departure point, at time 0, for each node on the current mesh. The reference solution is found by evaluating the function (3.4) at these departure points. Integration of the system of ODEs for departure points is computed using the SciPY (Scientific Tools for Python) odeint function, which itself is a wrapper around the lsoda ODE solver (Hindmarsh 2001).

#### (ii) Discretization

Consider again the general form of the advection–diffusion equation, this time with the advection term rewritten as the difference between a conservative term and a velocity divergence term
3.5
with all notation as in equation (3.3). In order to prevent spurious under-overshoots in the advection of the tracer field in this problem, a control volume bases method is used. Equation (3.5) is multiplied through by a test function *ϕ*, which, along with the trial space, is considered to be piecewise constant within control volumes, forming a dual mesh around the nodes of the original finite-element mesh. Clearly, this space does not support derivatives, so the advection and diffusion terms are integrated by parts leaving only surface integrals around each control volume (∂*Ω*_{CV}):
3.6
and
3.7
Note that as the diffusion term involves a second derivative, an auxiliary gradient variable, ** q**, is introduced. However, as this is governed by an equation (3.7) that is local to the control volume,

*Ω*

_{CV}, it can be substituted back into equation (3.7) implicitly.

However, the values of *u* and ** q** are undefined on the boundaries of the control volumes. For the diffusion term, we follow the method of Bassi & Rebay (1997) and take the average of the two control volumes on either side of the face for both the unknown

*u*and the auxiliary variable

**. This can be handled implicitly using the theta time-stepping method discussed previously.**

*q*For the advection terms, we use the parent finite-element mesh to find the interpolated values at the control volume faces. However, this is not guaranteed to produce a bounded flux. Hence, the face value is limited using the donor downwind and projected upwind values following the method of Sweby (1984). This is a nonlinear process, and so the advection terms are pivoted on a first-order upwind flux such that the face value, *u*_{f}, is defined as (LeVeque 1992)
3.8
Hence, *u*_{upwind} may be treated implicitly in time while nonlinear iterations update the nonlinear, , terms.

In this example, equation (3.3) with *κ*=0 is solved for the evolution of the tracer field until a time of 10^{7} s. The time step is taken to be 1000 s and *θ*=0.5. For the adaptive simulations, the mesh is adapted every 20 time steps. Linear interpolation is again used to transfer the solution field between adapted meshes.

#### (iii) Results

Figure 7 shows results at four time levels for three different fixed uniform isotropic mesh resolutions. Figures 8 and 9 show the corresponding results computed with anisotropic adapted meshes. Figure 9 shows the results using a lower interpolation error in the construction of the metric and hence has higher resolution. While there are some qualitative differences between the highest resolution fixed mesh results in figure 7 and those shown in figure 8, the differences are minimal compared with figure 9.

For quantitative comparison between the fixed and adapted mesh results, figures 10 and 11 show the time evolution of the maximum value of the advected field, for which the analytical solution preserves a value of 10, and the relative *L*^{2} norm (Hanert *et al.* 2004) of the difference between the numerical and the reference solution whose construction is described above. The results from figures 8 and 9 correspond to the lowest and middle resolution adaptive results in figures 10 and 11.

The fixed mesh results can be seen to diffuse more than the adaptive mesh when the tracer is close to the western boundary current, even for the highest resolution mesh where over 150 000 nodes are used. The adaptive mesh results lose some of the initial magnitude in the tracer field at the start of the simulation owing to the resolution of the initial mesh used before the first adapt, i.e. in the first 20 time steps. Improved results would be obtained if a higher resolution initial mesh was used or if the mesh was adapted at the start of the simulation to better represent the initial condition. At the end of the simulated time, the magnitude of the tracer field for the highest resolution mesh (157 045 nodes) can be seen to have dropped to below that of the middle adaptive resolution simulation, which corresponds to figure 9, i.e. 4303 nodes.

In the adaptive mesh results, the number of nodes used increases as the tracer field enters the western boundary current which is due to the gradients (and curvatures) increasing as the tracer field is stretched, and the number of nodes then decreases as the tracer field leaves the western boundary current region and the simulation becomes less computationally demanding. This demonstrates an important capability of adaptive methods: the ability to remove resolution where it is no longer required.

The loss in the magnitude of the tracer can be seen to be most significant while the tracer field is in the western boundary current region, which is due to the resolution being comparably lower as the scale of the advected field reduces, hence numerical diffusion increases. The effect is lower in the adaptive simulations as the mesh resolution is increased more significantly in this area at this time.

These results are consistent with the error plots in figure 11. The errors grow more rapidly when the tracer is being advected through the western boundary current. This is the case for both fixed and adaptive mesh results, where again the highest fixed-resolution mesh results match most closely with those obtained from the middle adaptive mesh simulation (figure 9), i.e. using over an order of magnitude fewer nodes.

With both the *L*^{2} error and the maximum value of the advected field, the fixed mesh results actually tend to perform better before and after the tracer is passing through the western boundary current, for example, as can be seen from the growth rates of the plotted lines in figures 10 and 11 at the start and end of the simulation. The adaptive mesh runs perform better, in terms of growth of errors, when the problem becomes more challenging numerically and the model is able to automatically use increased anisotropic resolution close to the western boundary.

The computational cost of the metric construction and the mesh optimization procedure are approximately equal for this problem. When combined they correspond approximately to the time taken for the model to complete three time steps for this problem. Hence, the computational overhead of adaptivity can be estimated at 15 per cent of the total run time. The middle adaptive mesh simulation (corresponding to figure 9) uses a maximum number of nodes, which corresponds closely to the number used by the coarsest fixed mesh simulation. However, the adaptive simulation runs approximately 60 per cent faster than the fixed mesh simulation.

### (c) Wind-driven barotropic gyre

#### (i) Problem description

In this problem, the unsteady velocity field in the domain *Ω*=*L*×*L* (*L*=1000 km) driven by a constant wind stress is computed. The two-dimensional constant-density Navier–Stokes equations in the following form are considered:
3.9a
and
3.9b
where ** u**=(

*u*,

*v*)

^{T}is the two-dimensional velocity vector,

*t*represents time and

*p*is the pressure. The (eddy) kinematic viscosity is given by

*ν*=1000 m

^{2}s

^{−1}. The rotation vector is defined in terms of the beta-plane approximation 2

**=(0,0,**

*Ω**f*

_{0}+

*βy*)

^{T}, where

*f*

_{0}=0 and

*β*=1.8×10

^{−11}. A wind stress is applied through the source term

**where ,**

*f**τ*=2×10

^{−7}and

*f*

_{y}=0.

#### (ii) Discretization

A description of the discretization based on continuous piecewise-linear basis functions for both velocity and pressure is given in Piggott *et al.* (2008*b*). This problem is simulated for 2 years with a time step of 6 h using a Crank–Nicolson method for the temporal discretization. For the adaptive simulations, the mesh is adapted every 80 time steps (20 days). Linear interpolation is used to transfer the solution field between adapted meshes.

#### (iii) Results

This problem has been simulated using a series of fixed uniform isotropic resolution unstructured meshes. The number of nodes in these meshes can be seen from table 1. In figure 12, the pressure field computed after 2 years of simulated time is presented. In figure 13, the magnitude of the computed velocity field is presented. In figure 14, the pressure and velocity magnitude fields and the corresponding adapted mesh are shown in an adaptive mesh simulation that had 2739 nodes. The adaptive mesh simulation can be seen to be recreating the correct dynamics of the problem, with much higher anisotropic resolution being used close to the western boundary where the flow is most intense.

Table 1 presents some diagnostic values for this problem. The adaptive mesh simulation can be seen to be doing a particularly good job in the representation of the maximum values of the north–south velocity component and the vorticity, both of which have their maxima in the western boundary current. They are able to do this with a reduced number of nodes because of the focusing in of the resolution near the western boundary, and the use of anisotropic refinement which further increases the efficiency of the approach and which was also seen in the previous two examples.

## 4. Concluding remarks

This paper has presented an optimization-based approach to anisotropic mesh adaptivity. An important step is the construction of an error metric that allows the size and shape quality of individual elements to be assessed, so that local operations may be applied to improve overall mesh quality. This approach is able to yield anisotropic meshes, and so a metric that incorporates information on the anisotropy of the solution variables being approximated is advantageous. Here an error metric inspired from interpolation theory was described, which is built up from the curvatures of solution fields. Future work will consider the inclusion of more powerful, and ideally problem specific, information in the metric definition (e.g. Micheletti & Perotto 2005; Power *et al.* 2006; Bernard *et al.* 2007). Even with the use of a relatively simple error metric, the advantages of adaptive resolution for three barotropic problems possessing western boundary currents were presented.

The adaptive approach presented here extends naturally to three dimensions. However, there are still open questions regarding the use of adaptive methods for large-scale stratified ocean dynamics: for example, whether the horizontal and vertical should be treated separately, for both the discretization method and mesh structure, and for which classes of problems an adaptive mesh approach coupled to an unstructured mesh model can be competitive with a highly optimized structured, finite-difference-based model. This is the subject of ongoing research by the authors and will be the subject of future publications.

Further related topics of research include the development of efficient load-balanced parallel algorithms for use with time-varying inhomogeneous resolution, and more accurate interpolation routines for passing solution fields between pre- and post-adapted meshes, in particular methods that are high order, bounded and conservative.

Finally, the approach presented here applied local topological operations to the mesh at discrete intervals throughout a simulation. Further increases in adaptive power, and complexity, result when continual mesh movement is also permitted in the model, for example to track isopycnal layers. The use of a spatially varying approximation order in the discretization of solution variables is also possible within a finite-element framework, and this can also be used to impart an adaptive capability in a model: for example, increasing the order of the method in regions where solution variables are smooth in an attempt to benefit from some of the exponential convergence properties of spectral methods. When multiple options are available for varying the approximation power of a model, via adaptive mesh algorithms, an important issue is how to automate the procedure by which the model is able to decide which combination of options is optimal to use at any point and at any time in the calculation. This leads back to the design of reliable metrics that are able to gauge not only what is going on in the simulated dynamics, but also what the implications are for the resulting numerical errors under different adaptive algorithms.

## Acknowledgements

The authors would like to acknowledge the funding of the UK Natural Environment Research Council under grants NE/C52101X/1, NE/F004184/1 and NE/F012594/1, the support of the Grantham Institute for Climate Change and the High Performance Computing Service at Imperial College London.

## Footnotes

One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.

- © 2009 The Royal Society