We transform near-local Hamiltonian balanced models (HBMs) describing nearly geostrophic fluid motion (with constant Coriolis parameter) into multi-symplectic (MS) systems. This allows us to determine conservation of Lagrangian momentum, energy and potential vorticity for Salmon's L1 dynamics; a similar approach works for other near-local balanced models (such as the -model). The MS approach also enables us to determine a class of systems that have a contact structure similar to that of the semigeostrophic model. The contact structure yields a contact transformation that makes the problem of front formation tractable. The new class includes the first local model with a variable Coriolis parameter that preserves all of the most useful geometric features of the semigeostrophic model.
The rotating shallow-water model, which describes the behaviour of a single layer of incompressible fluid with a free surface under gravity over a rotating bed, is frequently used in geophysical fluid dynamics as an approximation for the dynamics of the atmosphere. When the Rossby number is small (Ro≪1), a further approximation is to filter out the fast motions, driven by inertia gravity waves, to create so-called ‘balanced models’. Typically, these incorporate a balance between the pressure gradient and the Coriolis force.
Salmon (1983, 1985) showed how Hamilton's principle can be used to derive an important class of balanced models systematically. The idea is to define a constraint corresponding to a balance condition (the geostrophic approximation), and to incorporate this into the Lagrangian for the parent model, which is the two-dimensional shallow-water model. A variational principle is then used to obtain the dynamical equations. The accuracy of the resulting balanced model is judged against the parent dynamics, with the latter being considered as exact. The balanced model inherits a Hamiltonian structure and consequent conservation laws. By considering a class of balance conditions that included Salmon's L1 dynamics and semigeostrophic theory, McIntyre & Roulstone (2002) derived balanced models which go beyond geostrophy: when the Froude number is small (Fr≪1), gravity waves are fast compared with the flow velocity and the geostrophic approximation is not sufficient to yield a good approximation. One, called the -model, improves the accuracy of the approximation to the geostrophic flow (but not the advection velocity) by one order in Ro over semigeostrophic theory. The balance conditions considered in McIntyre & Roulstone (2002) and in this paper are referred to as near-local, because they depend on the geopotential and a finite number of its derivatives only.
Recently, Bridges et al. (2005) showed that the shallow-water model and the semigeostrophic model each have a multi-symplectic (MS) formulation; moreover, these two formulations are very similar. MS systems have a common geometric structure: a vector of closed 2-forms is conserved by the flow. (This generalizes the conservation of a single 2-form in canonical Hamiltonian dynamics.) The MS structural conservation law can be exploited to analyse stability (see Bridges (1997) for details), but in the context of balanced models, we focus on the scalar conservation laws that arise from it, namely conservation of Lagrangian momentum, energy and potential vorticity. Another advantage of the MS structure is that when it is preserved by a discretization, the resulting finite-difference methods normally have extremely good stability properties (see Bridges & Reich (2006) and Ascher & McLachlan (2005) for details). Until now, it has not been clear how to extend the shallow-water MS structure to other near-local balanced models. One purpose of this paper is to address this problem. We derive a MS description of L1 dynamics. McIntyre & Roulstone's other balanced models (including the -model) can be tackled by a similar approach, but as this adds considerable complexity without giving further illumination, we merely summarize the steps that lead to their MS formulation.
Although the semigeostrophic model is only first order, from the geometric viewpoint, it is the richest of the balanced models. It possesses both Hamiltonian and contact structures (see McIntyre & Roulstone (2002) for details), and recently Delahaies & Roulstone (2010) showed that the semigeostrophic model possesses a hyper-Kähler structure. The contact structure is particularly useful because it yields a contact transformation that expands the singularities which occur at fronts, and it is therefore able to describe the formation and evolution of fronts. However, like the other balanced models cited above, the semigeostrophic model has a constant Coriolis parameter. So it is natural to try to generalize this model in a way that retains the contact structure; this is the second aim of our paper. There have been several attempts to find models that preserve the structure of semigeostrophic theory and allow the Coriolis parameter to vary. Using Hamilton's principle, Salmon (1985) made the Coriolis parameter a function of geostrophic coordinates, while Shutts (1989) modelled planetary flow, replacing the geostrophic momentum with its projection onto the equatorial plane. Here, we build on the common features of the MS formulations for the shallow-water and semigeostrophic models, and we derive a class of MS systems with the same contact structure. One such system, which allows the Coriolis parameter to vary with position, has the same accuracy as the semigeostrophic model. By restricting this model to the β-plane, we show that this differs from the planetary model of Roulstone & Sewell (1997) by the inclusion of a single non-obvious term.
The plan of the paper is as follows. Section 2 describes the main features of near-local Hamiltonian balanced models (HBMs) with constant Coriolis parameters. A brief presentation of MS systems is given in §3, together with the MS formulations of the shallow-water and semigeostrophic models. In §4, Salmon's L1 dynamics is written as a MS system; we also outline how to recast the other balanced models in McIntyre & Roulstone's class as MS systems. Section 5 generalizes semigeostrophic theory by describing a class of MS systems that preserve energy, potential vorticity and Lagrangian momentum, and which also have a contact structure that yields analogues of Hoskins' geostrophic coordinates. This class includes the first local semigeostrophic-type model that allows the Coriolis parameter to vary with latitude, while retaining such a contact structure. This paper concludes with a brief discussion in §6.
2. Review of shallow-water balanced models
Following Bridges et al. (2005), we present the shallow-water theory in a Lagrangian formulation. The position and velocity of fluid particles are denoted by x=(x1,x2) and u=(u1,u2), respectively, and we label each fluid particle by its position at t=0, which is denoted by m=(m1,m2); then, all variables are treated as functions of m and t. The total derivatives with respect to t and mα are denoted by subscript ‘t’ or ‘α’ after a comma, for example, 2.1 and partial derivatives with respect to any other variable are written in full.
In shallow-water theory, the motion of a shallow layer of incompressible inviscid fluid over a rotating flat-bottomed bed is approximated by the horizontal momentum equations, 2.2 The Coriolis parameter is f (which is assumed constant, except in §5), the gravitational acceleration is directed downwards with magnitude g, and η(x,t) is the fluid height. The common shorthand k×(a,b)≡(−b,a) is used, and as a consequence of the Lagrangian formulation, we have 2.3 The incompressibility hypothesis requires η to satisfy the relation η dx=η0 dm, where the initial height η0=η(m,0) is assumed to be uniform. We choose coordinates in which η0=1, so that this condition yields the following form of the continuity equation: 2.4 Finally, the gradient with respect to x is given in terms of the Lagrangian formulation by 2.5 From equations (2.2)–(2.5), it is easy to show that the potential vorticity 2.6 is conserved following particles (see Bîlâ et al. (2006) for a thorough analysis). The total energy of a blob of fluid that occupies at t=0 is 2.7 given suitable boundary conditions, is also conserved.
Equations (2.2) and (2.3) are Euler–Lagrange equations for the variational principle where square brackets enclose the quantities that are varied. (In other words, [x,u] denotes the functions x, u and finitely many of their derivatives with respect to m and t.) The shallow-water Lagrangian is 2.8 here, H is the total energy density, 2.9 Conservation of potential vorticity and energy are consequences (via Noether's theorem) of the particle relabelling symmetry and the time invariance of the variational problem (see Salmon (1983)).
The semigeostrophic approximation to the shallow-water equations replaces the acceleration u,t in the momentum equations (2.2) by the Lagrangian time derivative of the geostrophic wind , which is defined by 2.10 where ϕ=gη is the geopotential. So the shallow-water semigeostrophic equations consist of equation (2.10), the momentum equation, 2.11 and the continuity equation (2.4). Equations (2.10) and (2.11) arise from a variational problem In this case, the Lagrangian is 2.12 where H is given by equation (2.9). Once again, the particle relabelling symmetry leads to conservation of potential vorticity, which is now defined by 2.13 Given suitable boundary conditions, the total energy of a fluid blob, 2.14 is also conserved because the variational problem is invariant under time translations.
The semigeostrophic approximation is a balanced model with many useful analytical and geometrical properties, including Hamiltonian and contact structures (see McIntyre & Roulstone (2002) for details). However, as it is formally correct only to first order in the Rossby number, its usefulness is limited. The question of how to build more accurate HBMs while retaining the essential mathematical features of semigeostrophic theory has been much studied. McIntyre & Roulstone (2002) provided one answer to this question using the framework of constrained Hamiltonian dynamics pioneered in Salmon (1983, 1985) and Allen & Holm (1996). The method is to construct a constrained Lagrangian Lc[x] by replacing u in the shallow-water Lagrangian Lsw[x,u] by a constraint velocity field of the form 2.15 In this formulation, uc cannot be varied independently of x; the Euler–Lagrange equation is the momentum equation, which must be supplemented by the constraint (2.15) and the continuity equation (2.4). For each , the resulting model conserves the potential vorticity, 2.16 and (given suitable boundary conditions) the total energy in a fluid blob, which is 2.17 Most attention has been focused on just three values of a. The semigeostrophic model corresponds to a=−1/2; it is a simple calculation to show that amounts to equation (2.13). For a=0, the constraint velocity reduces to the geostrophic wind, uc=ug; this gives Salmon's L1 dynamics. Finally, a=1 yields the so-called -model, which is the most accurate of the balanced models derived by McIntyre & Roulstone. In this case, the potential vorticity (2.16) agrees to second order with an asymptotic expansion of the shallow-water potential vorticity (2.6) in terms of the Rossby number. By contrast, L1 dynamics and the semigeostrophic model are each accurate only to first order (Snyder et al. 1991; Delahaies 2009).
The constrained Lagrangian formulations are incomplete because they do not contain the velocity constraint. In the next section, we show that by using a MS approach (which arises from a complete first-order Lagrangian), balanced models are given a common geometric structure. This structure provides local conservation laws, from which one can determine the necessary boundary conditions for a fluid blob to have a conserved quantity. It also reveals that conservation of potential vorticity is a differential consequence of two more fundamental conservation laws.
3. Multi-symplectic systems and the shallow-water equations
Bridges (1997) called a system of first-order quasi-linear PDEs MS if there is a (symplectic) closed 2-form associated with each independent variable, such that these 2-forms satisfy a structural conservation law (which generalizes conservation of symplecticity for Hamiltonian ODEs). For models based on the shallow-water approximation, we seek MS systems in a (local) Cartesian coordinate system with n dependent variables (z1,…,zn) and three independent variables (t,m1,m2). The structural conservation law is of the form 3.1 where 3.2 are closed. Here, summation is from 1 to n for Latin indices and from 1 to 2 for Greek indices; the functions Wij and are locally smooth. Hydon (2005) showed that the structural conservation law is a differential consequence of a 1-form quasi-conservation law of the form 3.3 for some locally smooth functions and S. In the above, the exterior derivative d acts only on the dependent variables zj; it commutes with the total derivative operators (see Bridges et al. (2010) for details). By comparing the exterior derivative of equation (3.3) with (3.2), one finds that Moreover, by expanding the quasi-conservation law and collecting the coefficients of dzi, one obtains the general form of a MS system (subject to the above restrictions), 3.4 A routine calculation shows that equation (3.4) is the set of Euler–Lagrange equations for the variational principle 3.5 where 3.6 So the set of MS systems is equivalent to the set of variational problems with a first-order Lagrangian that is affine linear in the derivatives of z. (Note: the above construction may be extended to deal with systems for which and S also depend upon t and m, but this is not needed for our purposes—see Bridges et al. (2010) for details.)
Given the above equivalence, it is straightforward to state Noether's theorem for the MS system (3.4), as follows. A vector field of the form generates variational symmetries if, for each sufficiently close to zero, the map leaves the variational problem (3.5) unchanged. Thus, the criterion for X to generate variational symmetries is 3.7 for some functions a and bα, where ⌋ denotes the inner derivative of the differential form by the vector field X. By taking the inner derivative with respect to X of the quasi-conservation law (3.3), we obtain the conservation law associated with X, namely 3.8
In particular, for every MS system (3.4), translations in t,m1 and m2 are variational symmetries (with Qi being , and , respectively); they yield the following three conservation laws: 3.9 3.10 and 3.11 Conservation of energy is represented by equation (3.9), whereas equations (3.10) and (3.11) describe conservation of Lagrangian momentum, that is, the quantity that is canonically conjugate to translations in label space. All three conservation laws may also be regarded as components of the pullback of the quasi-conservation law (3.3) to the space of independent variables. From this viewpoint, the dependent variables zi and their derivatives are treated as functions of t,m1 and m2; then, the coefficients of dt, dm1 and dm2 are equations (3.9), (3.10) and (3.11), respectively. The same approach can be applied to the structural conservation law (3.1), yielding conservation laws that are differential consequences of equations (3.9)–(3.11). In particular, the coefficient of dm1∧dm2 is the difference between the m1-derivative of equation (3.11) and the m2-derivative of equation (3.10), namely 3.12 This is the conservation law that corresponds (by Noether's theorem) to the particle relabelling symmetry, which is the infinite-dimensional pseudogroup of volume-preserving diffeomorphisms of label space. For shallow-water theory and for the balanced models that approximate it, equation (3.12) describes conservation of potential vorticity in terms of the Lagrangian coordinates.1 It is interesting that from the Lagrangian viewpoint, conservation of potential vorticity is merely a differential consequence of the (possibly more fundamental) conservation laws (3.10) and (3.11). However, those two conservation laws, unlike conservation of energy and potential vorticity, do not appear in the Eulerian viewpoint (because they cannot be written without reference to the particle labels).
We now review some relevant details of the MS version of the shallow-water model, which was derived in Bridges et al. (2005). The shallow-water Lagrangian (2.8) is not affine linear as it stands because it contains a multiple of However, an equivalent affine linear Lagrangian can be created by introducing new variables . It is also convenient to write the internal (potential) energy term gη/2 as 3.13
Inserting these new elements into equation (2.8), using Lagrange multipliers to enforce the constraints , we obtain the affine linear Lagrangian 3.14 Then, the Euler–Lagrange equations obtained by varying δx1, δx2, δuα, and in turn (with fixed endpoint conditions) are 3.15 3.16 3.17 3.18 and 3.19 The MS structure emerges when we reorganize this system of first-order partial differential equations as 3.20 Here, 3.21 and W is the 12×12 skew-symmetric matrix where W is the 4×4 block 3.22 and the skew-symmetric matrices Kα, α=1,2, are defined by 3.23 with and 3.24 Using equation (3.18) and (3.19) together with equation (2.5), we see that so equation (3.13) yields 3.25 Inserting these expressions in equations (3.15) and (3.16) leads to the horizontal momentum equations (2.2), which shows that the system (3.15)–(3.19) is equivalent to the shallow-water model.
The MS formulation of the semigeostrophic model is almost the same as the above (with ug replacing u), except that the affine linear Lagrangian is 3.26 and so equation (3.22) is replaced by 3.27
The structural conservation law and the 1-form quasi-conservation law for each of these models is obtained by substituting the components into the general forms (3.1)–(3.3) (see Bridges et al. (2005), Hydon (2005) and Cotter et al. (2007) for details).
4. Multi-symplectic formulation of constrained dynamics
The purpose of this section is to derive a MS version of McIntyre & Roulstone balanced models. In §3, we derived a Lagrangian that enabled us to recast the shallow-water model into MS form. We now seek to insert the constraint u=uc into this Lagrangian. The process is described in full for L1 dynamics and is summarized for McIntyre & Roulstone's class of more general constraints (2.15).
(a) L1 dynamics in multi-symplectic form
To recast L1 dynamics into a MS form, we need to insert the constraint u=ug, written in terms of the MS variables, into the Lagrangian (3.14). For the shallow-water model, equation (3.25) gives where Recall that in this setting, the variables are introduced as the Lagrange multipliers corresponding to the constraints . We want to play the same role in the derivation of L1 dynamics, so that they remain as Lagrange multipliers without being used in the constraint u=ug. Therefore, we introduce auxiliary functions , which depend on the variables only, as follows: The constraint u=ug can now be written as 4.1 it depends only on uα, and first-order derivatives of . Inserting this constraint with a vector of Lagrange multipliers v=(v1,v2) into the Lagrangian defining the parent dynamics leads to the affine linear Lagrangian 4.2 where now Here, The Euler–Lagrange equations corresponding to the variations δxβ, δuβ, , δvβ and are 4.3 4.4 4.5 4.6 4.7 4.8 4.9 and 4.10 As shown in Delahaies (2009), this MS system is equivalent to Salmon's L1 dynamics.
(b) Conservation laws for the multi-symplectic L1 model
Just as for the shallow-water model, the closed 2-forms (3.2) are and The dm1∧dm2 component of the pullback of the structural conservation law is , where is the geostrophic potential vorticity (that is, the potential vorticity defined by equation (2.16) with uc=ug). Similarly, the local energy conservation law, which is the dt component of the quasi-conservation law, is It implies that, provided suitable boundary conditions are given, the total energy , defined by equation (2.17) with uc=ug, is conserved. The components of Lagrangian momentum are derived similarly.
(c) Multi-symplectification for McIntyre & Roulstone's balanced models
To recast McIntyre & Roulstone's other balanced models into an MS form, we apply the same technique as for L1 dynamics, that is, we apply the constraint u=uc to the Lagrangian density Lsw. However, by contrast with L1 dynamics, additional functions are needed for the constraint to fit into the MS structure. Recall that the constraint velocity is 4.11 The -model is obtained when . Expanding equation (4.11) using equations (2.5) and (4.1), the constraint u=uc can be written in components as 4.12 and 4.13 The above expressions involve second-order derivatives (via the second derivatives of ), so as to create an affine linear first-order Lagrangian, we introduce new variables, 4.14 and new functions, 4.15 Then, as 4.16 equations (4.12) and (4.13) amount to 4.17 and 4.18 This ensures that the model obtained by adding the constraints (4.14), (4.17) and (4.18) to the Lagrangian (3.14) can be written in MS form (see Delahaies (2009) for further details).
5. A contact-preserving generalization of semigeostrophic theory
Despite its limited accuracy, the semigeostrophic theory has remarkable geometrical properties, such as a contact structure, Legendre duality and a Monge Ampère structure. These properties have been used to show that solutions exist (see Benamou & Brenier (1998)). In this section, we derive a class of MS systems which generalizes semigeostrophic theory in a way that retains a contact structure and a modified version of the Hamiltonian structure that can be obtained by using Hoskins' geostrophic coordinates. By restricting attention to members of the class that have the same qualitative features as the semigeostrophic model, we find a new model that allows the Coriolis parameter f to vary with latitude. Throughout this section, u denotes the constraint velocity (which is ug when f is constant); it is defined by the pair of diagnostic equations (which lack time derivatives).
(a) Geometric features when f is constant
The semigeostrophic equations with constant f have several geometric features that are revealed by using geostrophic coordinates, 5.1 Hoskins (1975) showed that, from the Eulerian viewpoint, 5.2 where ϕ=gη is the geopotential and 5.3 Note that equation (5.2) follows immediately from the definition of the geostrophic wind by the diagnostic equations 5.4 The consequence of equations (5.2) and (5.4) is that, for each t, the transformation is a strict contact transformation because
From the Lagrangian viewpoint, semigeostrophic dynamics is governed by 5.5 At first sight, this appears to be a canonical Hamiltonian system whose Hamiltonian is f−1Φ. However, each ξα is a length, which contravenes the usual idea in mechanics that the Hamiltonian is energy and that the dependent variables are position and momentum. For semigeostrophic flow, Φ is an energy density (energy per unit mass), so it seems sensible to use this as the Hamiltonian and to have as dependent variables a position and its canonically conjugate momentum density. Bearing in mind that we later intend to allow f to depend on the latitude variable x2, we use ξ1 as the position. Then, with 5.6 equation (5.5) amounts to the canonical system 5.7 Note that, by construction, the Hamiltonian Φ is independent of the Coriolis parameter.
The following observations are also helpful. First, using (A,B) as variables simplifies the time-derivative part of the MS Lagrangian (3.26), so that 5.8 Second, using equations (3.2) and (3.6) with equation (3.21), the time component of the MS 2-form reduces to so that the potential vorticity is Third, as f is a constant, we have not altered the contact transformation property by using (A,B) instead of (ξ1,ξ2).
(b) An extension that preserves the Hamiltonian and contact structures
A major distinction between the MS shallow-water and semigeostrophic equations is the rank of W, which is 4 and 2, respectively. Consequently, for the shallow-water equations, there are four prognostic (dynamic) equations, and two pairs of dependent variables are needed. For semigeostrophic flow, however, one can find a single pair of variables (A,B) that are dependent variables in the prognostic equations; moreover, the diagnostic part of the semigeostrophic equations, as expressed by the contact transformation, requires only the same single pair of variables. These facts lie behind the useful geometric structures that Hoskins' coordinates first revealed.
It makes sense, therefore, to see which other systems of equations have the same features in common. We shall leave the Hamiltonian Φ unchanged, but we allow A and B to be arbitrary independent functions of x and u. Just as for the semigeostrophic model with constant f, the MS Lagrangian is 5.9 because we want Φ to remain unchanged. Then, the Euler–Lagrange equations corresponding to variations of δxβ, δuβ, and are 5.10 5.11 5.12 5.13 5.14 5.15 This system can be written as equation (3.20), except that the non-zero block W is now Equation (5.15) amounts to 5.16 and therefore, system (5.10)–(5.13) is equivalent to 5.17 and 5.18 which can be expressed as 5.19 As before, W is of rank 2, so we can treat Φ as a function of (A,B) for each t. Consequently, the system is again Hamiltonian, 5.20 Substituting equation (5.20) into equations (5.17) and (5.18), we find that once again and so the transformation is a strict contact transformation.
So whenever the MS Lagrangian is of the same form as for constant-f semigeostrophic flow, the contact and Hamiltonian structures are preserved. Furthermore, the time component of the MS 2-form is again and thus there is an analogue of potential vorticity that is preserved by the flow, namely
(c) Variable-f equations
The discussion above is entirely mathematical, without reference to any physical constraints. We now consider how to perturb the constant-f semigeostrophic variables (A,B) in a way that allows the Coriolis parameter to vary with x2. As we have seen, the main problem is that there is a great deal of freedom. Here is one rationale that leads to a solution; we cannot claim that this is the only or best solution, but it has the advantage of introducing errors that are no larger than the errors in the constant-f case.
For constant f, As neither of these depend explicitly on x2, we seek to preserve these equations. A straightforward calculation shows that this requires for some functions F and G. Suppose that, as for constant f, A−u1 is independent of u2 and B−x1 is a multiple of u2. Then, On dimensional grounds, the functions and must have the same dimensions as fx2 and 1/f, respectively. Indeed, for constant f, Even with all of these constraints, there remains a gauge freedom: we can add a constant to without changing the explicit form of the equations. We may exploit this freedom by choosing The integral is indefinite, so A and B are local functions, as required.
This is a new variable-f approximation, which preserves that contact and Hamiltonian structure of the constant-f system. Furthermore, a routine calculation shows that the error introduced by this approximation is of precisely the same order of magnitude as the error owing to the neglect of planetary curvature. Consequently, this model is no worse than the constant-f semigeostrophic model.
If f(x2)=f0+βx2, the systems (5.22)–(5.25) provide a β-plane semigeostrophic model that differs slightly from the β-plane version of the planetary geostrophic equations derived in Roulstone & Sewell (1997), principally because it contains a term which ensures that the contact structure is preserved.
6. Summary and concluding remarks
We have adapted Salmon's approach of incorporating a balance condition into a variational principle to the MS framework and thus obtained a local formulation of near-local balanced models. Starting with Salmon's L1 dynamics, we then considered the more general class of near-local constraints used in McIntyre & Roulstone (2002), including the -model. These descriptions are not unique: instead of inserting the balance condition into the Lagrangian through the use of Lagrange multipliers, we could have replaced the velocity field by the suitably expressed balance condition directly into the Lagrangian. As shown for L1 dynamics and the generalization of semigeostrophic theory, the local formulation paves the way for constructing conservation laws from the structural conservation law.
In this paper, we have only considered near-local constraints, in which the constraint is expressed in terms of the local value of the geopotential and a finite number of its derivatives. Higher accuracy typically requires non-local constraints; examples include the non-local second-order balanced models derived in Allen & Holm (1996) and Vanneste & Bokhove (2002).
We have also proposed a MS system that generalizes the system presented in Bridges et al. (2005) for semigeostrophic theory. This generalization enabled us to present a local extension of the f-plane semigeostrophic theory to variable Coriolis parameters, and we have proved that this generalization carries a contact structure. The generalization of semigeostrophic theory presented in this paper was influenced by our prior knowledge of the formulation of semigeostrophic theory in terms of canonical coordinates, namely Hoskins' geostrophic coordinates. McIntyre & Roulstone (2002) found that complex canonical coordinates exist for the class of near-local constraints considered here; they lead to interesting complex geometries. It would be useful to extend our approach to this more general class of complex canonical coordinates.
We thank Prof. V. Roubtsov and Prof. I. Roulstone for useful discussions and for their very helpful comments on the manuscript. We also thank anonymous referees for their useful comments and suggestions. The work of S.D. was supported (in part) by a NERC studentship.
- Received March 6, 2011.
- Accepted June 20, 2011.
- This journal is © 2011 The Royal Society