z˜ -Coordinate, an Arbitrary Lagrangian–Eulerian coordinate separating high and low frequency motions

Ocean Modelling 37 (2011) 139–152

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

~z-Coordinate, an Arbitrary Lagrangian–Eulerian coordinate separating high and low frequency motions Matthieu Leclair a,⇑, Gurvan Madec a,b a b

Laboratoire d’Océanographie et du Climat: Expérimentation et Approches Numériques (LOCEAN), 4, place Jussieu, Paris, France National Oceanography Centre Southampton (NOCS), United Kingdom

a r t i c l e

i n f o

Article history: Received 11 May 2010 Received in revised form 21 January 2011 Accepted 1 February 2011 Available online 13 February 2011 Keywords: ALE Vertical coordinate NEMO Spurious diapycnal mixing

a b s t r a c t A new geopotential based vertical coordinate, named ~z, is presented, that treats external motions as a z⁄coordinate, internal low frequency motions in a Eulerian way and high frequency ones in a Lagrangian way. To that purpose, a prognostic equation for the low frequency horizontal flux divergence is added in the system in order to discriminate between low and high frequency parts of vertical motions. This coordinate keeps the intrinsic advantages of the z-coordinate while enabling a drastic reduction of the spurious diapycnal mixing induced by vertical and lateral schemes truncation errors associated with high frequency motions in level models. An idealised experiment illustrates this property in an internal wave propagation framework. The reduction of Eulerian vertical velocities enables a proportional reduction of the associated advection schemes truncation errors. Following high frequency motions, the ~z-coordinate lets lateral schemes operate along tracer anomalies created by those motions.  2011 Elsevier Ltd. All rights reserved.

1. Introduction The choice of a vertical coordinate system is one of, if not the most, important aspects of an ocean model’s design (Griffies et al., 2000a). Until recently, ocean models based on primitive equations almost exclusively used one of the three main vertical coordinates: the geopotential, terrain following or isopycnal coordinate. A detailed review of those coordinates is given by Griffies et al. (2000a). The z-coordinate models have some key advantages linked to the predominance of lateral over vertical transport in the ocean and hence to its essentially horizontal structure. The highly nonlinear equation of state is easily and accurately represented, geopotential levels are a natural framework for the computation of the horizontal pressure gradient, avoiding truncation errors and the surface mixed layer is naturally parametrized. z-Models have also the non-negligible advantage of easy and efficient implementation. However, some drawbacks are also inherent to this vertical coordinate. The stepwise representation of topography is not accurate, although some alternative methods like shaved cells (Adcroft et al., 1997) or partial cells (Pacanowski and Gnanadesikan, 1998) successfully address this issue, especially when used in association with an energy–enstrophy conserving momentum advection scheme (Penduff et al., 2007; LeSommer et al., 2009). The representation of down-slope flows or overflows is difficult (Legg ⇑ Corresponding author. Tel.: +33 476825045. E-mail address: [email protected] (M. Leclair). 1463-5003/$ - see front matter  2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2011.02.001

et al., 2009) and horizontal levels are unnatural for the description of diffusion processes occurring along inclined iso-neutral surfaces in the ocean interior, leading to spurious diapycnal mixing (Griffies et al., 2000b). Isopycnic models are a very efficient remedy to these disadvantages enabling notably a perfect control of the diapycnal mixing. Unfortunately, this quality is obtained at the expense of the main advantages of z-models. A q-coordinate is inappropriate for the representation of the surface mixed layer, considering an accurate empirical equation of state is difficult and so-called truncation errors in the computation of the horizontal pressure gradient may occur, though satisfactory addressed (Adcroft et al., 2008). Isopycnal and geopotential coordinates are thus quite complementary. The last main classical vertical coordinate used in ocean models is the terrain-following or r-coordinate. The advantages of this model class are mainly concentrated at the ocean bottom. Topography is smoothly represented and allows a good description of the bottom boundary layer. However, the spurious diapycnal mixing problem is more critical than in z-models (Marchesiello et al., 2009), and although some methods reduce them (Shchepetkin and McWilliams, 2003; Marsaleix et al., 2009), truncation errors in the computation of horizontal pressure gradient are still an issue. There is obviously no perfect vertical coordinate. Hence, the recent general tendency of the ocean modelling community is to develop generalised vertical coordinate transformations also called Arbitrary Lagrangian–Eulerian (ALE) vertical coordinates. Two strategies are used for the practical realisation of such coordinates.

140

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

We call them Quasi Lagrangian (QL) and Quasi Eulerian (QE) algorithm hereafter. Each of them is adapted, and used, when developing an ALE coordinate starting respectively from an initially Lagrangian (q-coordinate) or Eulerian (z, r-coordinates) model. The (QL) algorithm was introduced by Hirt et al. (1974) and consist in three phases:  A Lagrangian phase where vertical motions are described by vertical displacements of computational surfaces.  A Regridding phase where the new vertical grid is computed.  A Remapping phase where the fields previously updated in a Lagrangian way are remapped on the new grid. The (QE) algorithm (Kasahara, 1974) can be decomposed in two stages:  A Regridding stage during which the new vertical grid is calculated. The Eulerian velocity, i.e., the vertical velocity relative to the mesh, is also obtained during this stage.  A Eulerian stage during which vertical motions are described by advection through computational surfaces, the advection velocity being the Eulerian velocity previously computed. The HYCOM model, mostly based on a q-coordinate, was the first model that tried to use a mix of z and q vertical coordinates within a single framework (Chassignet et al., 2003) and is the prototype of ocean models using the QL algorithm: primitive equations are first solved in a fully Lagrangian way and results are then re-mapped on a vertical grid that tends to be isopycnal in the ocean interior and evolves toward geopotential and terrainfollowing coordinates respectively in the vicinity of ocean surface and bottom (Bleck, 2002; Halliwell, 2004). More recently, White et al. (2009) developed a vertical ALE coordinate using this strategy with high order regridding and remapping schemes. Several vertical coordinates have been developed using the QE algorithm. The first of them is the r-coordinate which is not fully Eulerian. r levels are rescaled proportionally to the sea surface height, so that barotropic motions are treated in a Lagrangian way. This principle has been applied to the geopotential coordinate, giving birth to the z⁄-coordinate (Adcroft and Campin, 2004; Stacey et al., 1995). Burchard and Beckers (2004) and recently Hofmeister et al. (2010) developed an adaptive r-based coordinate allowing vertical resolution to follow strong shear and stratification regions. The new vertical coordinate proposed here enters the class of the ALE coordinates using the QE algorithm. This coordinate, which we name ~z, differs from the previously mentioned coordinates as it distinguishes between the Eulerian and Lagrangian parts in the temporal rather than in the spatial domain. It aims at following high frequency vertical motions like internal and barotropic waves in a Lagrangian way, as an isopycnal model would do, and at allowing cross-level velocities, i.e., a z-model like Eulerian description of the vertical flow, for low frequencies. It thus can be viewed as a generalisation of the z⁄-coordinate to high frequency baroclinic motions. The coordinate is designed so that it keeps the zcoordinate advantages at low frequencies and benefits of the isopycnal coordinate qualities at high frequencies. Staying in a zcoordinate like framework at low frequencies allows the truncation error on the hydrostatic pressure gradient to remain negligibly weak and preserves the easy handling of model outputs (no interpolation onto z-level is required) while the Lagrangian behaviour at high frequencies addresses the spurious diapycnal mixing issue present in z-level models. Moreover, recent studies showed a general predominance of high frequencies in the vertical velocity spectrum, especially in the deep ocean (Komori et al., 2008). The article is structured as follows. The primitive equations are first briefly described in an ALE coordinate (Section 2). Then the

~z-coordinate is defined in the continuous time and space domain (Sections 3.1 and 3.2) followed by its discretization (Section 3.3). A numerical experiment illustrates the reduction of spurious diapycnal mixing in Section 4. Section 5 summarises the results and concludes. 2. The primitive equations in a general vertical coordinate The primitive equations are solved in the NEMO model in a horizontally orthogonal curvilinear coordinate system. A time varying vertical coordinate, like the r- or z⁄- coordinate can be used, so that the model equations (1a)–(1e) are treated in the framework of a general vertical coordinate (Madec et al., 1996). Note, that between the two possible options, the so-called ‘‘vector invariant’’ and ‘‘flux’’ forms, we present here the first one, the most frequently used with this model, that enables enstrophy conservation (Madec et al., 1991; LeSommer et al., 2009; Barnier et al., 2009). Contrary to Madec (2008) that derives the primitive equations in an ALE vertical coordinate in the tensorial formalism for the NEMO model, they are presented here in a more ‘‘standard’’ way, as for instance in Adcroft and Hallberg (2006), in order to avoid requiring a supplementary effort from the reader. Finally, we do not consider here the treatment of the external mode, which is not in the scope of the present study and does not interfere with the ~z-coordinate formulation.

@uh 1 w @u 1  ðrh p ¼ ð f k þ rh  uh Þ  uh  rh juh j2  2 e3 @k q0 @t þ qg rh zÞ þ Duh þ Fuh

ð1aÞ

@ðe3 hÞ @ðhwÞ ¼ rh  ðe3 huh Þ  þ Dh þ F h @t @k

ð1bÞ

@e3 @w ¼0 þ rh  ðe3 uh Þ þ @k @t

ð1cÞ

@p ¼ qge3 @k

ð1dÞ

q ¼ qðT; S; pÞ

ð1eÞ

Here k is the local upward vertical unit vector and k 2 R is a nondimensional vertical coordinate defined such that its integer values match the discrete tracer point of the C-grid when (1) is discretized. It varies from kb, the bottom value, to ks, the surface value. e3 is the vertical scale factor, also called the thickness in reference to its discrete counterpart. It is given by e3 = @ kz, where z is the local depth, and depends on the three-dimensional space coordinates and on time. rh = (@ x, @ y, 0) is the lateral derivation operator,1 i.e., acting along k = cst surfaces. uh is the horizontal velocity, w is the Eulerian velocity across k = cst surfaces, h could be either the potential temperature T or salinity S, q is the in situ density and p the pressure. q0 is a constant reference density, f is the Coriolis factor and g the acceleration of gravity. Duh and Dh are the diffusive terms, Fuh and Fh the forcing terms. Finally, we also introduce here the lateral local transport divergence or equivalently the lateral thickness weighted divergence, D, defined by

D ¼ rh  ðe3 uh Þ

ð2Þ

which will be widely used hereafter. 1 The h subscript indicates the horizontal direction and is abusively kept here for simplicity even if the operator does not exactly operate in the horizontal direction.

141

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

3. The ~ z-coordinate 3.1. The ~z-coordinate concept Let us first describe how vertical motions are treated in the different models. In a generalised coordinate framework, the vertical velocity is the sum of a ‘‘Lagrangian’’ part, @ tz, the velocity of the k-surfaces and a ‘‘Eulerian’’ part, w, the velocity across the ksurfaces. Note that these two velocities are not strictly speaking the Lagrangian velocity, i.e., the velocity of a fluid parcel, and the Eulerian velocity, i.e., the velocity relative to a fixed coordinate system. That is the reason why we call them ‘‘Lagrangian’’ part and ‘‘Eulerian’’ part here. However, in the following, they will be called abusively ‘‘Lagrangian’’ velocity and ‘‘Eulerian’’ velocity for simplicity reasons. They both appear in the continuity equation (1c) by their derivation with respect to k, i.e., @ te3 = @ t(@ kz) and @ kw. Eulerian and Lagrangian vertical velocities are entirely determined by the association of the continuity equation and a second equation, thereafter called the coordinate equation, that depends on the vertical coordinate used:  In pure vertically Lagrangian models, like isopycnal ones, k-surfaces are material surfaces. There is no Eulerian vertical velocity, so that the coordinate equation is w = 0 and (1c) is used as a prognostic equation for e3: @ te3 = D, where D is defined by (2). D is entirely absorbed by the vertical scale factor evolution.  In pure vertically Eulerian models, like z-coordinate ones, k-surfaces remain fixed in time so that there is no vertical Lagrangian velocity. The coordinate equation is @ te3 = 0 and (1c) is used as a diagnostic equation for w: @ kw = D. D is entirely counterbalanced by the vertical Eulerian velocity.  In ALE vertical coordinate models, an entire process we call here coordinate equation2 determines the vertical coordinate. In models using the QE algorithm, @ te3 is deduced from the coordinate equation and (1c) is then used to diagnose the Eulerian velocity: @ kw = D  @ te3, so that D is absorbed partly by the Lagrangian vertical velocity and partly by the Eulerian one. Note that the determination of w is not required in models using the QL algorithm, as the vertical advection is implicitly done during the remapping phase. The z⁄-coordinate is a particular example of ALE coordinate using the QE algorithm. Indeed, it can be formulated as follows.3 The lateral local transport divergence is split into a barotropic and a baroclinic part, D = D⁄ + D0 , where D⁄ is deduced from the vertical integral of D and D0 is its departure from D⁄:

D ¼

e3 H0 þ g

D0 ¼ D  D

Z

ks

D dk

ð3Þ

kb

ð4Þ

where H0 = z(kb) is the ocean depth and g = z(ks) the sea surface height (ssh hereafter). We precise here that the definition of the word barotropic is personal as the resulting ‘‘barotropic’’ divergence D⁄ does depend on k. The depth averaged horizontal divergence, which is also the divergence of the barotropic transport divided by the water column thickness, is indeed multiplied by the local vertical scale factor. This ‘‘barotropic’’ divergence is thus more precisely defined as the local contribution to the barotropic transport. Then, the vertical scale factor and vertical Eulerian velocity are respectively defined by the coordinate equation, @ te3 = D⁄, and the continuity equation, @ kw = D  @ te3 = D0 . In other words, ver2 The so-called coordinate equation can be a single equation, as in the z⁄-coordinate case, but also a much more complex regridding algorithm. 3 In this subsection, we ignore the fresh water forcing for simplicity reasons.

tical displacements associated with barotropic motions are treated in a Lagrangian way, while baroclinic ones are treated in a Eulerian way. The ~z-coordinate we introduce here is also a particular case of the QE algorithm and can be viewed as a generalisation of the z⁄coordinate to high frequency baroclinic motions. Let us further e and a low frequency part, split D0 into a high frequency part, D, 0  0 e D , so that D ¼ D þ D þ D . The philosophy of the ~z-coordinate LF

LF

is that vertical motions associated with both the barotropic and high frequency baroclinic parts of the flow are treated in a Lagrangian way while the remaining low frequency baroclinic part is treated in a Eulerian way. The resulting equations are:

@e3 e ¼ D  D @t @w @e3 ¼ D  ¼ D0LF @k @t

ð5Þ ð6Þ

e and ensure the robustness The questions are now how to define D of the coordinate in a realistic global ocean context. This is the topic of the next subsection. 3.2. Practical formulation of the ~z-coordinate In order to discriminate between high and low frequencies, D0LF is introduced in the system as a new prognostic variable. Many ways can be considered to derive D0LF . We choose here a simple restoring equation towards the total baroclinic transport lateral divergence D0 , which can also be seen as a first order low-pass filter:

 @D0LF 2p  0 D  D0 ¼ @t s~ LF

ð7Þ

e is then sim~, the cutoff period of the filter, is a parameter. D where s ply deduced from D0LF by:

e ¼ D0  D0 D LF

ð8Þ

In practise, the formulation (5) and (6) is slightly complicated in order to ensure the robustness of the coordinate and a clean separation between high frequency internal variations and external variations of e3. Let us decompose the vertical scale factor e3 as

e3 ¼ e03 þ f e3 þ e3

ð9Þ

where:  e03 is the scale factor associated with the background z-coordinate. e03 is horizontally homogeneous and does not depend on time. Its vertical integral is the ocean depth H0. f e3 is the part of e3  e03 that evolves with internal high frequency motions, it oscillates around zero and its vertical integral is zero.  e3 is the part of e3  e03 that evolves with external motions, it remains small ðe3  e03 Þ, as in a z⁄-coordinate, and its vertical integral is equal to g, the ssh. e03 being fixed in time, two prognostic equations have to be established for f e3 and e3 to entirely define the vertical coordinate. Let us first consider the f e3 prognostic equation. Starting from the schematic overview given previously, the prognostic equation for e Nevertheless, a more robust coordinate f e3 should be @ t f e3 ¼  D. is obtained by the addition of two other terms. The actual prognostic equation of f e3 becomes:

@ðf e3 Þ e3 þrh  ½j ~ rh ð f e 2spz f e3 Þ ¼ |{z} D |fflffl{zfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} @t r1

r2

r3

ð10Þ

142

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

where  r1 is the dominant term of the equation, enabling the vertical coordinate to follow high frequency baroclinic motions as an isopycnal coordinate would do.  r2 is a restoring term of f e3 towards 0, i.e., a restoring of the total coordinate towards the z⁄-coordinate, with the restoring time scale sz being a parameter. This term prevents any potential long term drift of the coordinate. In the real ocean, a single point will indeed not necessarily see an entire number of periods of e term. As a typthe high frequency motions contained in the D ical example, a front passing at that point will cause the coordinate to move to a new position without coming back. The r2 term ensures that the ~z-coordinate will remain close to the background coordinate at time scales larger than sz. This restoring term amounts for an additional first order filtering of vertical motions so that the high frequency motions captured by the f e3 evolution result in the end from a second order filtering (see Appendix A).  r3 is an optional term, the so-called thickness diffusion term, ~ . This term is also present which eddy diffusivity coefficient is j in layer models and enables the dissipation of energy cascading down to small space scales which could otherwise lead to grid point noise in the k-surfaces. Note that f e3 is oscillating around zero, so that, in opposition to the isopycnal coordinate, the diffusion operator is not required to be strictly positive and a more scale selective operator, like a bi-harmonic one, could be used. As mentioned before, f e3 is of baroclinic nature, i.e., its vertical integral is zero. This property ensures a perfect separation of the treatment of vertical motions associated with either the high frequency baroclinic part of the flow or its barotropic part. In particular, such a separation prevents the restoring and thickness diffusion terms, r2 and r3, from affecting the ssh. The prognostic ~ equation (10) preserves the baroclinic nature of f e3 as soon as j does not depend on k and the thickness diffusion fluxes are set to zero through the domain boundaries. As shown in Appendix B, Rk one can form a prognostic equation for kbs f e3 dk whose right hand side terms are a restoring and a diffusive term. The horizontally integrated variance of that variable is thus necessarily decreasing. Rk As it starts from zero, one easily shows that "(x, y, t), k s f e3 ðx; y; k; b tÞ dk ¼ 0. The prognostic equation of e3 requires some special care. In a pure z⁄-coordinate framework, e3 is defined as the redistribution of the0 ssh proportionally to the constant background scale factor: e e3 ¼ H30 g. In the ~z-coordinate framework, the relative proportions of vertical scale factors are not constant in time (due to f e3 evolution), therefore e3 must be defined by its time tendency equal to the redistribution of the ssh time tendency proportionally to the scale factor:

@ðe3 Þ e3 @ g ¼ @t H0 þ g @t

ð11Þ

Rk Integrating (11) vertically, it is obvious to show that kbs @ t ðe3 Þ dk ¼ R ks @ t g. Since @ t g ¼ kb D dk þ Q , where Q is the fresh water forcing, (11) can be rewritten as @ t e3 ¼ D þ He03þQg . In other word, treating the external mode as in a z⁄-coordinate model means representing barotropic vertical motions in a Lagrangian way, except for fresh water forcing. ~Þ ! 0, it There are two limit cases for the ~z-coordinate. If ðsz ; s can be easily demonstrated that f e3 ¼ 0 and e3 is given by   e3 ¼ e03 þ e3 ¼ e03 1 þ Hg0 , i.e., corresponding to the z⁄-coordinate. ~Þ ! 1, (7) yields D0LF ¼ 0 and the r2 term in Conversely, if ðsz ; s (10) vanishes. Furthermore, if both fresh water forcing and

e so thickness diffusion terms are ignored, @ t e3 ¼ D and @ t f e3 ¼ D, that equations (8), (10) and (11) eventually reduce to @ t e3 ¼ e þ D ¼ D. The ~z-coordinate becomes a pure @ t ðf e3 þ e3 Þ ¼ D ~Þ can thus be seen as a pair of tuning Lagrangian coordinate.4 ðsz ; s parameters defining which part of the baroclinic temporal spectrum is considered as high frequency and thus treated in a Lagrangian way and which part is low frequency and thus subjected to a z-coordinate like treatment, the whole barotropic spectrum being treated in a z⁄ manner (except when fresh water forcing is taken into account). The ~z-coordinate is now entirely defined by the coordinate equations (10) and (11). The vertical Eulerian velocity should thus be simply given by the continuity equation but a slight modification is done here. Indeed, the thickness diffusion term r3 in (10) is designed to smooth the coordinate when k-surfaces tends to be too noisy but we do not want this term to generate a compensating vertical Eulerian velocity. To that purpose, r3 is reformulated as a lateral advective correction (Gent and McWilliams, 1990; ~h, Bleck et al., 1992): introducing a thickness diffusion velocity, u the prognostic equation (10) becomes:

j~ rh ð f e3 Þ e3 @ðf e3 Þ e  2p f ~hÞ e  rh  ðe3 u ¼ D @t sz 3

~h ¼  u

ð12Þ ð13Þ

~ h is taken into account in the Then, the thickness diffusion velocity u modified continuity equation:

@w @e ~ h Þ  3 ¼ rh  ½e3 ðuh þ u @k @t

ð14Þ

This behaviour is inspired by isopycnic models where there is no Eulerian velocity even if a thickness diffusion term is applied. Eventually, satisfying the local consistency between volume and tracer budgets (Griffies et al., 2005; Leclair and Madec, 2009) the tracer equation (1b) is modified to take the same advective correction as in (14) into account:

@ðe3 hÞ @ðhwÞ ~ h Þ  ¼ rh  ½e3 hðuh þ u þ Dh þ F h @t @k

ð15Þ

The repartition of the local lateral transport divergence either absorbed by the vertical scale factor variation (Lagrangian part) or balanced by the across k-surface velocity (Eulerian part) is now entirely defined. In order to demonstrate the filtering properties of the ~z-coordinate, a study of the coordinate frequency response is conducted in a linear and baroclinic context and the concepts of gain and phase originating from the electronics domain are introduced (Bode and July, 1940). We define the complex amplification factors K1, K2 and K3. K1 is the ratio between low frequency lateral divergence and full lateral divergence. K2 represents the part of the baroclinic lateral divergence that is absorbed by the scale factor evolution, i.e., the ratio between scale factor evolution in the ~z-coordinate and in a fully isopycnic one. And symmetrically, K3 is the part of the baroclinic divergence absorbed by the Eulerian cross-level velocity or, in other words, the ratio between Eulerian velocity in the ~z-coordinate and in a standard z one.   The gain and phase Gdb i ; Ui ¼ ð20log10 ðjKi jÞ; ArgðKi ÞÞ are then used to obtain the so-called Bode diagrams (Bode and July, 1940). The method used to calculate them is detailed in Appendix A and they are plotted in Fig. 2. The two upper panels constitute the Bode diagram of K1. The characteristic 20 db per decade slope in the Amplitude diagram and  P2 phase for high frequencies 4 Note that the pure Lagrangian coordinate is just considered as an extreme case. It does not mean that the ~z-coordinate could be used in an isopycnal configuration as non-trivial processes like vanishing layers are not taken here into account.

143

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

w f u f

f

v v

u

T w

f

confirm that D0LF is defined by a first order low-pass filter. The Bode e3 diagram of K2 is given by the two middle panels. It shows that @ t f indeed absorbs the high frequency part of the lateral divergence. The asymptotic behaviour for low frequencies shows a characteristic second order high-pass filter. Finally the two lower panels constituting the Bode diagram of K3 show a first order low-pass filter, confirming that @ kw absorbs the low frequency part of the baroclinic lateral flow divergence. 3.3. Discretization of the ~z-coordinate

Fig. 1. Repartition of the calculation points in the staggered C-grid. T-points are located at (i, j, k), U-points are at (i + 1/2, j, k), V-points at (i, j + 1/2, k), W-points at (i, j, k + 1/2) and F-points (the vorticity points) at (i + 1/2, j + 1/2, k).

Although the ~z-coordinate can be used in any time and space discrete framework, we present here the particular case of the

Period [day]

Period [day]

Period [day]

Period [day]

Period [day]

Period [day]

Fig. 2. Bode diagrams presenting the spectral analysis of the ~z-coordinate. Amplitude and phase diagrams are respectively plotted in the first and second column. The first row shows the amplification factor of D0LF , the low frequency baroclinic flux divergence, compared to D0 , the full baroclinic flux divergence. The second row represents the amplification factor of the vertical scale factor variation @ te3 compared to a full layer model, where the total divergence is absorbed in the vertical scale factor variation. Eventually, the third row, shows the amplification factor of the Eulerian velocity w compared to a z-model where the total horizontal divergence generates Eulerian vertical velocity. The x axes are graduated in period in days, the gain axes are in decibel (db) and the phase axes are in degrees. For all curves, the cut period of the low frequency 2p e ¼ 100 m2 s1 . The different curves ~ ¼ 5 days, the wave number is k ¼ 200 baroclinic horizontal divergence, has been set to s km and the thickness diffusion coefficient is K correspond to different settings of sz: in black when sz = 100 days, in blue when sz = 30 days and in red when sz = 5 days. The Bode diagram corresponding to D0LF contains only one curve as the amplification factor does not depend on sz. The vertical dashed lines correspond to the cut frequencies of each amplification factor. The frequency of the waves generated in the test case of Section 4, 0.5 day, is represented on each plot with a dotted line. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

144

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

NEMO model which is discretized in space on a three-dimensional Arakawa-C-type grid using a modified Leap-Frog time stepping scheme associated with a Robert–Asselin time filter (Leclair and Madec, 2009). 3.3.1. Space discretization Equations (13) and (11) defining the prognostic equation for the vertical scale factor e3 are computed at T-points on the grid (Fig. 1). The interpolation of e3 at U, V and W-points is determined by energetic considerations. It can be demonstrated that in the nondissipative primitive equations without forcing and with a linear equation of state,  Property 1: the integrated kinetic energy variation over the domain is equal to the integrated power of the pressure gradient forces.  Property 2: the integrated potential energy variation over the domain is equal to the opposite of the integrated power of the pressure gradient forces. As shown in Madec (2008), preserving property 1 in the discrete equations defines the interpolated vertical scale factors at U and V-points as follows:

ator, preserve a quasi second order accuracy for the vertical advection scheme. This choice however breaks the energetic consistency considered above. The results and implications of these two formulations or discussed in Section 4. 3.3.2. Time discretization The NEMO model is discretized in time using a leap-frog time stepping scheme associated with a Robert–Asselin time filter. The time stepping sequence corresponding to the non-dissipative continuous equations (3)–(8) and (11)–(15) without forcing is described hereafter. The f lower script will denote time filtered quantities and the n upper script is the time step. First the baroclinic lateral transport low frequency divergence e is deduced: ðD0LF Þ is computed and the high frequency one ð DÞ

Dn ¼

en3 H 0 þ gn

Z

ks

Dn dk

ð19Þ

kb

D0n ¼ Dn  Dn

ð20Þ

0n1 D0n  2p LF ¼ DLF

Dt 

sd

D0n1  D0n LF



e n ¼ D0n  D0n ; D LF

ð21Þ ð22Þ

e3iþ1=2 ¼

ð16aÞ

so that the baroclinic part of the vertical scale factor can be advanced in time:

e3jþ1=2

ð16bÞ

~h ¼  u

 1 1 ðe1 e2 e3 Þi þ ðe1 e2 e3 Þiþ1 ðe1 e2 Þiþ1=2 2 i 1 1h ðe1 e2 e3 Þj þ ðe1 e2 e3 Þjþ1 ¼ ðe1 e2 Þjþ1=2 2

where e1 and e2 are the horizontal scale factors, in other words, the grid spacing in the two horizontal directions. i and j are the indexes in these two directions. The interpolation of e3 at W-points is more complex as it is linked to the interpolation of density at those points. We propose here two formulations favouring either energetic consistency or the accuracy of the vertical advection. The first one is obtained by preserving property 2. To our knowledge, the only solution found so far in the literature is to associate the following interpolations for density (tracers) and vertical scale factors at W-points (see for example Marsaleix et al., 2008):

 1 e3 þ e3kþ1 2 k  1 q e3 þ qkþ1 e3kþ1 ¼2 k k e3kþ1=2

e3kþ1=2 ¼

ð17aÞ

qkþ1=2

ð17bÞ

However, this conservative solution is obtained at the expense of the accuracy of the vertical advection scheme: when vertical resolution is not constant, the classical second order centred advection scheme degenerates into a first order one (Yin and Fung, 1991). Marti et al. (1992) and Tréguier et al. (1996) indeed explain that second order accuracy is maintained as soon as T- and W-points depths and scale factors are derived from an analytical continuous function and tracers are interpolated at grid interfaces with a simple half averaging. The second solution given by (18a) and (18b) is designed in order to preserve a quasi second order advection scheme:

e3kþ1=2 ¼ e03 ðk þ 1=2Þ þ 1 2

qkþ1=2 ¼ ½qk þ qkþ1 

 1 e3  e03 ðkÞ þ e3kþ1  e03 ðk þ 1Þ 2 k

ð18aÞ ð18bÞ

where k#e03 ðkÞ is an analytical continuous function defining the background z-grid. Within the ~z-coordinate framework, vertical scale factors and depths oscillate around their background values. Hence (18a), which aims at interpolating e3 at W-points around its initial analytical value and (18b), the direct half averaging oper-

f e3 nþ1

j~ en3



e3 n1 rh f f



ð23Þ



 n  n1 e n  2p f ~ ¼f e3 n1 þ 2 D t  D  r  e e u 3f h f 3 h

sz

ð24Þ

Then the barotropic part e3 , and thus the entire scale factor e3, are updated:

e3nþ1 ¼ e3 n1 þ f

 en3  nþ1 g  gn1 f n H0 þ g

¼ e03 þ f e3 nþ1 þ enþ1 enþ1 3 3

ð25Þ ð26Þ

Even if f e3 and e3 are prognostic variables, only one of them needs to be stored. As the entire vertical scale factor e3 has to be stored anyway, the second one is recovered from (9). Here we chose to store f e3 . Then w can be computed and tracers advanced in time:

   enþ1  en1 3f ~h  3 @ k wn ¼ rh  en3 unh þ u 2Dt     nþ1 n1 ~ h hn  @ k ðwn hn Þ enþ1 ¼ en1  2Dt½rh  en3 unh þ u 3 h 3f hf

ð27Þ ð28Þ

Eventually, the Robert–Asselin time filter is applied to the vertical scale factors and tracers:

  f e3 nf ¼ f e3 nþ1  2f e3 n þ a f e3 n þ f e3 n1 f    2en3 þ en1 en3f ¼ en3 þ a enþ1 3 3f   nþ1 n1  2en3 hn þ en1 en3f hnf ¼ en3 hn þ a enþ1 3 h 3f hf

ð29Þ ð30Þ ð31Þ

where a is a small non-dimensional parameter that controls the magnitude of the time filter. In the presence of forcing, this filter has to be modified in the way described by Leclair and Madec (2009) in order to ensure global and local tracer conservation. A remarkable point highlighted by the time stepping sequence is that only the very simple prognostic equations (21) and (23) imply additional computations compared to the already available z⁄-coordinate option in NEMO, so that the model’s cost does not significantly increase.

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

145

4. Numerical experiments: reduction of the spurious diapycnal mixing The expected benefit of this coordinate that follows high frequency vertical motions in a Lagrangian way and low frequency ones in a Eulerian way is to reduce the spurious diapycnal mixing. This property is illustrated in this section in a test case that aims at showing the behaviour of the ~z-coordinate compared to a classical z-coordinate in an internal wave propagation framework. The geometrical, physical and numerical settings of this configuration are first described. The results in two particular situations are then detailed in order to exhibit vertical advection schemes truncation errors. Finally, a quantification of the spurious diapycnal mixing is proposed and the influence of wave amplitude and lateral diffusion is addressed. 4.1. Test case configuration The test case presented here is an f-plane 600 m deep flat bottom basin, 1500 km wide, closed in the zonal direction and periodic (i.e., infinite) in the meridional direction. The latitude of the f-plane is 30N so that the inertial period is 1 day. The resolution is Dx  Dz = 2.5 km  20 m. The equation of state is linear and density depends only on temperature. The initial temperature field is horizontally uniform. Its profile is given in Fig. 3: the initial temperature is constant for the first 200 m then linearly decreases for the next 200 m and again constant over the last 200 m. An ‘‘internal wave generator’’ is placed at the centre of the basin. It consists in the addition of a momentum forcing in the zonal direction, i.e., an acceleration, A(x, z, t) = A0M(t)F(x, z). The spatial distribution of this forcing and its temporal modulation are given by (32b) and represented in Fig. 4:

8      < sin 2p x sin p zþH0  1 jx  x0 j 6 k2x kx 2 gþH0 Fðx; zÞ ¼ :0 jx  x0 j P k2x

ð32aÞ

  81 1  cos 2p 6Tt > 2 >  > t <1    MðtÞ ¼ sin 2p 1 T > 1 þ cos 2p t6T > 6T > :2 0

ð32bÞ

Fig. 4. Illustration of the wave generator. The top panel shows the spatial distribution of the zonal momentum forcing confined in the [25 m, 25 m] interval. The thick arrows represent the acceleration vector and the thin lines represent the model grid. The temporal modulation of momentum forcing is plotted on the bottom panel. The time axis is graduated in number of internal wave periods. The amplitude is smoothly growing during the first 3 periods, then kept constant during the 3 following ones and smoothly decreasing again during 3 periods.

0 6 t 6 3T 3T 6 t 6 6T 6T 6 t 6 9T t P 9T

where T = 0.5 days is the wave period and kx = 50 km is the horizontal wave length. The two first tests performed in order to validate the implementation of the ~z-coordinate are the two extreme behaviours ~ ¼ Dt and described in Section 3.2: when setting 2s~p ¼ 2spz ¼ 0 or s sz = 2Dt, the model actually matches respectively an isopycnic or a z⁄-coordinate one. Fig. 5 shows the layer displacements when

Fig. 5. Illustration of the wave trains propagation after 7 days. The experiment has been performed with a full Lagrangian setting of the ~z-coordinate, i.e., restoring ~ are infinite. The model level displacements are plotted in the weak periods sz and s amplitude case on the top panel and in the strong amplitude case on the bottom e ¼ 20 m2 s1 is applied in one. For strong amplitude waves a thickness diffusion of K order to remove the grid-point noise arising in the unstratified layers in the forcing zone. Results obtained for the ~z-coordinate with settings described in Section 4.1.3 are very close to the full Lagrangian results shown here.

Fig. 3. Initial temperature profile in the numerical test case of Section 4. The temperature is homogeneous over the top and bottom 200 m and linearly stratified in between the temperature field is homogeneous in the horizontal direction.

the model is used as an isopycnic one and illustrates the two wave trains that are generated at the centre of the basin and travel zonally. Their amplitude is determined by A0: the isopycnal

146

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

displacement reaches 5 m in the weak amplitude case and 50 m in the strong one. Weak and strong amplitude are differencing two regimes in the geopotential coordinate framework: in the first case, the isopycnal displacement is limited to the adjacent cells, whereas the strong amplitude waves advect water through several grid points. The baroclinic nature of the waves is ensured by the fact Rg that H0 Fðx; zÞ dz ¼ 0. As a consequence, results with a z or z⁄ coordinate are almost identical. Therefore, only the z⁄ results are given here. A set of experiments has been performed to compare the behaviour (especially the spurious diapycnal mixing) of the ~z and z⁄ coordinate models. They are summarised in Table 1 and differ through four settings: the wave amplitude already discussed, vertical advection scheme, lateral diffusion and naturally the vertical coordinate. The two possible discretizations of the ~z-coordinate that are described in Section 3.3 favouring either the energetic consistency (equation (17)) or the vertical advection scheme’s accuracy

Table 1 Summary of the different experiments performed in the internal wave test case of Section 4. The advection scheme column designates both vertical and horizontal advection schemes. CEN2 stands for the centred second order scheme and TVD for the total variation diminishing scheme. The amplitude of the generated waves, either 5 m or 50 m is illustrated on Fig. 5. Experiment

Advection scheme

Wave amplitude (m)

Lateral bilaplacian diffusion (m4 s1)

CEN2_small_NoDiff CEN2_small_Diff CEN2_big_NoDiff CEN2_big_Diff TVD_small_NoDiff TVD_small_Diff TVD_big_NoDiff TVD_big_Diff

CEN2 CEN2 CEN2 CEN2 TVD TVD TVD TVD

5 5 50 50 5 5 50 50

0 1  109 0 1  109 0 1  109 0 1  109

(equation (18)) are respectively referred to as conservative and non-conservative ~z from now on. 4.1.1. Advection schemes The two most used vertical advection schemes available in NEMO that have been used for the comparison are the dispersive centred second order scheme (CEN2 hereafter) and a diffusive TVD scheme. The TVD scheme is actually a centred scheme that is mixed to an upstream scheme to prevent the occurrence of local extrema, ensuring total variance decrease and positivity of the scheme. These two schemes are also described in Lévy et al. (2001). In each experiment the same scheme is used for horizontal and vertical advection. The sensitivity to the horizontal advection scheme has also been explored by comparing the simulations performed with a TVD scheme on the vertical but with different advection schemes on the horizontal. It appears that it does not significantly affect the solution. This is due to the fact that, in the wave propagation framework, zonal velocities and temperature anomalies are in phase. There are thus no zonal velocities across the maximum zonal temperature gradients and, conversely, zonal velocities are maximal where the zonal temperature gradients vanish. Therefore the corresponding experiments are not discussed in this paper. However, horizontal advection schemes would certainly have an impact if the waves propagated in a mean horizontal current. In order to still address the impact of lateral schemes, a lateral diffusion is applied with a varying intensity. 4.1.2. Friction and diffusion An iso-level bilaplacian lateral friction operator is applied to the momentum equation, with a coefficient of 1  109 m4 s1, which is sufficient to ensure model stability in this test case even without any vertical friction and lateral or vertical diffusion on tracers. However, lateral diffusion will be applied through an iso-level bilaplacian operator with a varying intensity. Hence, although the

Fig. 6. Temperature anomalies obtained with the weak amplitude wave trains. The first row corresponds to the CEN2 advection scheme (experiment CEN2_small_NoDiff), the second one to TVD (experiment TVD_small_NoDiff). The first column are results obtained with the z⁄ coordinate and the second one with the conservative ~z-coordinate. The non-conservative version is not shown as it differs very slightly from the conservative one in the weak amplitude case. The colourbar is graduated in degrees Celsius. Model outputs are not rescaled on a pure z-grid before plotting so that the depth axis corresponds to the initial depth of the point. The fact that one can almost see nothing in ~zcoordinate plots when using the same colour scale shows first that the vertical motions are mainly treated in a Lagrangian way and second that the truncation errors appearing in the wake of the wave trains are drastically reduced by the ~z-coordinate. Only the eastern half of the basin is shown because the experiment is exactly symmetrical. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

impact of horizontal advection schemes cannot be studied in this wave propagation framework, lateral diffusion still enables to address the influence of lateral schemes on spurious diapycnal mixing. 4.1.3. Vertical coordinate In the ~z-coordinate case, some parameters have to be chosen. As ~ and sz constitute a explained in Section 3.2, the two time scales s kind of cursor determining the behaviour of the model between a ~ that separates low z⁄ and an isopycnic one. The cutoff period s and high frequencies for the baroclinic lateral flux divergence has been set to 5 days so that the generated waves are considered as high frequency waves by the model (see Fig. 2: the wave frequency is plotted on the Bode diagrams). The restoring period sz for the baroclinic vertical scale factors f e3 has been set to 30 days. Both these periods are typical values that could be used in a realistic simulation. They could have been chosen much closer to 0 so that the model behaves much closer to an isopycnic one, removing more of the spurious diapycnal mixing. However, the purpose of this section is to illustrate the behaviour of the ~z-coordinate in realistic conditions, especially the fact that the coordinate should remain geopotential at low frequency, so that the chosen parameters are realistic ones. Finally, no thickness diffusion is applied in the case of the weak 5 m amplitude waves and it has been set to j~ ¼ 20 m2 s1 for the strong 50 m amplitude waves, a value sufficient to remove the grid point noise arising in the forcing region in the unstratified layers. 4.2. Vertical advection schemes truncation errors: CEN2 versus TVD Of the experiments summarised in Table 1, two are detailed in this section, CEN2_small_NoDiff and TVD_small_NoDiff. Fig. 6

147

shows the temperature anomalies after 7 days that are created by the 5 m amplitude wave trains travelling through the basin. No lateral diffusion is applied and the results with the CEN2 and TVD advection schemes are plotted. Only the conservative ~z-coordinate is considered here because the difference between conservative and non-conservative is not distinguishable when level thickness variations are small. (Note that it would not be the case if the initial scale factors were not homogeneous on the vertical.) Corresponding profiles in the centre part of the basin are plotted on the top left panels of Figs. 7 and 8. They are obtained by zonally averaging the temperature anomaly field after 7 days in the interval [100 km, 100 km]. At this stage of the simulation, wave trains are out of the averaging domain. The most prominent impression when looking at Fig. 6 is that almost no temperature anomaly can be seen in the ~z case. It illustrates two different things. Firstly, the fact that the wave train propagation is invisible with the ~z-coordinate illustrates that model layers are actually following high frequency vertical motions so that almost no Eulerian cross-level velocities are generated and consequently almost no temperature anomalies. Secondly, the signature appearing in the wake of the wave train is also absent from the ~z plots. In this part of the basin, model computational surfaces are flat in both cases. Temperature anomalies in this region are thus due to truncation errors that occurred at the top and bottom of the stratification layer during the wave propagation, generating spurious diapycnal mixing. The latter is hence drastically reduced in the ~z case. In the case of the CEN2 scheme, the unstable temperature maximum created by upward advection is spread up to the surface by the hydrostatic adjustment, creating a spurious diapycnal mixing. The downward advection is thus no more compensated by the upward advection and a negative bias appears at the first point above

Fig. 7. Temperature anomaly profiles obtained with the CEN2 advection scheme. Results are plotted in black for the z⁄-coordinate, in blue for the conservative ~z-coordinate and in red for the non-conservative one. In order to address the influence of wave amplitude and lateral diffusion, the four experiments combining weak or strong wave amplitude and the presence or not of lateral diffusion are treated on the four panels. Each point of the profile is obtained by averaging the temperature anomaly field zonally in the interval [100 km, 100 km] after 7 days of simulation. At this stage, the wave trains are out of the averaging interval so that temperature anomalies are only due to truncation errors. Note the different temperature axis scale for weak and strong amplitudes. Note also that the extreme case of a pure Lagrangian ~z-coordinate is not shown here as, by definition, temperature anomalies are zero in that case. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

148

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

Fig. 8. Same as Fig. 7 with the TVD scheme.

the stratified layer. Hence the gradient has been eroded at the top of the stratified layer. The scheme acts symmetrically under the stratified layer. Hence a spurious diapycnal mixing occurs locally at the stratified layer boundaries but the large scale gradient between top and bottom waters has been reinforced. In the case of the TVD scheme, a spurious vertical diffusion arises at the frontier between stratified and unstratified layers, due to the diffusive upstream scheme that is used to avoid the creation of extrema by the centred scheme. These phenomena, spurious vertical mixing due directly to diffusive vertical advection schemes or indirectly to dispersive schemes leading to convecting adjustment are already known issues (see e.g. Griffies et al., 2000b). When using the ~z-coordinate these two sources of spurious diapycnal mixing are strongly reduced because high frequency Eulerian vertical motions are much smaller than in geopotential models (10 times smaller in this particular case, which is exactly what the frequency analyses of Appendix A shows) inducing weaker truncation errors for vertical advection schemes. 4.3. Evaluation of spurious diapycnal mixing In order to quantify the benefit of the ~z-coordinate in terms of diapycnal mixing and its sensitivity to wave amplitude and lateral diffusion, zonally averaged temperature anomaly and equivalent vertical diffusion coefficient profiles have been computed for each experiment (see Figs. 7–9). These experiments do not include the extreme case of the ~z-coordinate in the pure Lagrangian configuration. As the initial stratification is horizontally homogeneous, model layers are exactly following isopycnal surfaces and obviously no spurious diapycnal mixing could be diagnosed in that case. These equivalent diffusion coefficients are obtained as follows. A vertical temperature profile is computed by laterally averaging the temperature field in the interval [100 km, 100 km] after 7 days of simulation. At this stage wave trains have already propagated outside this interval and the remaining temperature anomalies are only due to truncation errors. The vertical heat fluxes F that explain the laterally averaged temperature variation from

the initial state satisfying @ tT = @ zF are then computed from the bottom to the surface, with the bottom boundary condition F(z = H) = 0. Finally, we assume that the temperature variations are due to a simple vertical diffusion process: @ tT = @ z(Fd) where Fd = Kz@ zT is the vertical diffusive flux and Kz the equivalent diffusion coefficient. Assuming F = Fd yields the vertical diffusion profile K z ¼  @zFT 0 where the vertical temperature gradient is taken at the initial state. Kz can thus only be computed where the initial temperature gradient is non-zero, i.e., in the stratified layer (200 m 6 z 6 400 m). A drawback of this method is that the equivalent diffusion diagnostic is not relevant for the CEN2 scheme. In the CEN2_ small_NoDiff experiment discussed previously, the heat content corresponding to the positive bias under the stratified layer is indeed almost exactly compensated by the negative bias spread down in the water column. The heat flux obtained at the frontier between stratified and unstratified layers when summing up the temperature variations from the bottom is thus automatically zero and so does the co-located diffusion coefficient. Hence, no diffusion is diagnosed precisely where the initial gradient erosion occurs. This is symmetrically observed in the top half of the basin. Moreover, the spurious diapycnal mixing due to convective adjustment arising in the two initially homogeneous layers cannot be diagnosed this way neither as the initial temperature gradients are zero in those layers. Therefore, only the equivalent diffusion coefficient profiles obtained with the TVD experiments (Fig. 9) are discussed here. 4.3.1. Influence of wave amplitude Experiment TVD_small_NoDiff (top left panel of Fig. 9) with a small 5 m wave amplitude shows that, in quite linear conditions, spurious diapycnal mixing is essentially concentrated in the stratification changing zones. As explained in Section 4.2 it is due to the truncation errors of vertical advection scheme. The diagnosed equivalent diffusion coefficient reaches almost 1 cm2 s1 with the z⁄-coordinate and is reduced to about 0.15 cm2 s1 by the ~z-coordinate. No influence of the conservative or non-conservative version

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

149

Fig. 9. Equivalent diffusion profiles with the TVD scheme. They correspond to temperature anomaly profiles of Fig. 8. They are only computed at points where the initial temperature gradient was non 0, i.e., within the stratified layer [200 m, 400 m]. Results are plotted in black for the z⁄-coordinate, in blue for the conservative ~z-coordinate and in red for the non-conservative one. As explained in Section 4, this diagnostic is not relevant for the CEN2 scheme so that only results with the TVD scheme are shown. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

of the ~z-coordinate is evident under these conditions. These two formulations are indeed nearly equivalent when the vertical scale factors are initially homogeneous on the vertical and do not depart far from this initial value. The first information given by the 10 times stronger 50 m amplitude wave experiments (experiment TVD_big_NoDiff, bottom left panel of Figs. 8 and 9) is the diffusive nature of the conservative ~z-coordinate. When vertical scale factors are no longer homogeneous in the vertical, which happens by using the ~z-coordinate when high frequency motions have a strong amplitude or the initial vertical grid is not uniform, the vertical advection scheme corresponding to the discretization given by (17a) and (17b) is indeed a first order diffusive scheme. The equivalent diffusion coefficient is not negligible (it reaches 0.7 cm2 s1 in the stratified layer where the non-conservative discretization stays very close 0) so that only the quasi second order non-conservative formulation is retained. The second information given by the strong amplitude experiments is that the TVD and CEN2 schemes are dispersive within the stratified layer: one can observe oscillations of the temperature anomaly and the equivalent diffusion profile (only for the TVD scheme) around 0 in this region. This is a natural behaviour for the CEN2 scheme but less intuitive for TVD. The latter becomes diffusive only where local extrema would be created by the CEN2 scheme. The strong stratification in the central layer does not permit this local extrema generation so that CEN2 and TVD act in the same way in the stratified layer. One can then notice that the reduction of the spurious diapycnal mixing arising at the stratification changing zones seems less important when using the ~z-coordinate than in the weak amplitude case. The equivalent diffusion, which reaches 3.5 cm2 s1 in the z⁄ case is reduced to 1 cm2 s1 by the use of the ~z-coordinate. The temperature gradient is actually diffused over more than one grid cell in the vertical direction when using the z⁄-coordinate

(see bottom left panel of Figs. 7 and 8) due to the stronger Eulerian vertical velocities. Using the ~z-coordinate reduces the vertical Eulerian velocities enough to hold the temperature anomalies within one grid point around the stratification changing zones, whereas the z⁄-coordinate permits the advection over three grid points around the top and bottom of the stratified layer. 4.3.2. Influence of lateral diffusion In order to address the impact of lateral schemes on spurious diapycnal mixing, a bilaplacian operator is applied with a realistic (1  109 m4 s1) coefficient for this 2.5 km horizontal resolution (see Lévy et al., 2010). Experiments TVD_small_Diff and TVD_big_Diff (right panels of Fig. 9) clearly show that lateral diffusion only influences diapycnal mixing when the z⁄-coordinate is used. The equivalent vertical diffusivity appearing in the stratified layer remains quite low compared to the stratification changing zones in the weak amplitude wave case but becomes totally dominant in the strong amplitude one: the equivalent diffusivity profiles is almost constant around 4.7 cm2 s1. Experiments not shown here with a 10 times stronger lateral diffusion have also been performed. Spurious diapycnal mixing with the ~z-coordinate became sensitive to lateral diffusion but far less than the z⁄ model. The equivalent diffusivity in the centre of the stratified layer reached 0.7 cm2 s1 for the non-conservative ~z-coordinate, 1.25 cm2 s1 for the conservative one and 33 cm2 s1 for the z⁄ coordinate. In the test case presented here, the nearly insensitive character of the spurious diapycnal mixing diagnosed in ~z-coordinate simulations to the presence of a lateral diffusion operator is not surprising. Indeed, starting from a horizontally homogeneous stratification, model layers are almost following isopycnal layers. This would not be a robust result in the real ocean as the slopes between mean state isopycnal and geopotential surfaces are not negligible regarding the impact of lateral operators on diapycnal mixing. However, the variance of any field along a model layer

150

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

due to high frequency vertical motions and, hence, the associated diapycnal mixing (caused either directly by diffusion operators or indirectly by diffusive advection ones), will still be cancelled by the use of the ~z-coordinate. Note also that the creation of local extrema by the lateral bilaplacian diffusion operator can be seen in the experiment TVD_big_Diff (bottom right panel of Fig. 8) in the z⁄ case: there are no longer two opposite biases at the stratified layer boundaries but a single stronger bias in the unstratified layer explained by a downward diffusive flux (the effect of the diffusion operator across moving isotherms) and the appearance of a small opposite bias (as for the CEN2 scheme, it is explained by the local extrema created by the bilaplacian operator spread into the stratified layer by the convective adjustment). 5. Conclusions and perspectives A new geopotential based ALE vertical coordinate called ~z has been presented that aims at reducing the spurious diapycnal mixing associated with high frequency vertical motions in level models. The coordinate distinguishes the Eulerian and Lagrangian parts of the vertical flow in terms of frequencies: low frequency vertical motions are described by the Eulerian method in a z-coordinate framework and high frequency vertical motions are followed in a Lagrangian way. This coordinate only requires about 10% additional cpu time as only two easily handled prognostic equations are introduced in the system, one for the low frequency baroclinic lateral divergence and one for the vertical scale factor. The tilde concept applied here in a geopotential framework could also be used in any r-coordinate model without any complication. The ~z-coordinate has been tested in an internal wave propagation framework. An important spurious diapycnal mixing is diagnosed when the geopotential coordinate is used and several causes are identified. The truncation errors of vertical advection schemes contribute to the diapycnal mixing, either directly when diffusive schemes are used or indirectly via spurious convective adjustments when dispersive schemes are employed. This phenomenon is strongly reduced by the use of the ~z-coordinate that only treats a very small part of high frequency vertical motions in a Eulerian way and hence limits the truncation errors of the vertical advection schemes. The ~z-coordinate also maintains temperature anomalies in the adjacent grid points, even in the case of strong amplitude waves, whereas the geopotential coordinate allows vertical advection through several grid points and thus

diffuses the initial temperature gradient over a larger vertical area. The influence of lateral schemes are then addressed. Although the contribution of lateral advection could only be studied if the waves propagated in a mean current, which we have not considered here, the impact of lateral diffusion has been quantified. This gives a qualitative idea of the influence of lateral processes in general. On the one hand, the study shows that, when using a geopotential coordinate, realistic lateral diffusion does not influence much the solution when the wave amplitude is small and becomes the dominant source of spurious diapycnal mixing in the strong amplitude case (the equivalent diffusivity reaching almost 5 cm2 s1). On the other hand, the ~z-coordinate shows a total insensitivity to lateral diffusion due to the diffusion scheme working in the iso-coordinate direction. The ~z-coordinate follows high frequency motions, so that lateral schemes have almost no effect on internal waves propagating in a horizontally homogeneous fluid. The ~z-coordinate presented in this paper is a first implementation and can be further improved. The filtering strategy chosen to obtain the low frequency lateral divergence could for instance be refined by using a second order filter taking the two previous time steps into account. A more complex issue is the discretization of density and scale factors at W-points in the C-grid. As far as we know, the only energetically consistent solution enabling the transfer between kinetic and potential energy through vertical advection is diffusive in the vertical. A conclusion of the present study is that it becomes too diffusive when strong amplitude waves are propagating and it would be especially problematic if the background grid were not uniform. A diffusive scheme is thus the price to pay for the energetic consistency and the elaboration of more precise energetically consistent schemes is still an open issue. The next question naturally arising about the ~z-coordinate is its ability to operate in a realistic ocean context. Some preliminary tests have been performed in order to address this issue. The ~z-coordinate has been used without any noticeable problem in a global ocean-sea ice model at a 2 horizontal resolution, namely the ORCA2-LIM configuration (see Timmermann et al. (2005) for details), that incorporates more complex physics (non-linear equation of state, isoneutral diffusion for instance) and geometry (partial cell bottom topography and coasts) than in the test case of Section 4. The expected result of a strong reduction of the Eulerian vertical velocity is obtained. This is illustrated on Fig. 10: the use of the ~z-coordinate reduces the squared vertical Eulerian velocity by a factor reaching 16 during the first days of the simulation where a

Fig. 10. Zonal mean of the ratio between the squared Eulerian vertical velocity in the z⁄ and the ~z experiments with the ORCA2-LIM configuration (Timmermann et al., 2005) ~ ¼ 5 days for the averaged over the first 5 days, i.e., during the adjustment phase from rest to the Levitus climatology. The parameters used for the ~z-coordinate tuning are s ~ ¼ 0, i.e., no thickness diffusion is separation between high and low frequency baroclinic motions, sz = 30 days for the restoring towards the background z-coordinate and j used. The contour interval is 2.

151

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

lot of internal waves are generated by the adjustment from rest to the Levitus climatological thermohaline structure. The most obvious application (and planed for a close future) is the influence of the ~z-coordinate on regional configurations at high resolution including tidal forcing. Interactions of tides and bathymetry indeed create internal waves, so that high frequencies are a significant part of the spectrum of total vertical motions in such experiments. Although a weaker impact is expected when considering quasi geostrophic turbulence, the ~z-coordinate will certainly be more efficient than the z⁄-coordinate as associated vertical structures are not purely barotropic. As soon as eddies will be considered as high frequency, which depends on their translation velocity, their characteristic length (the Rossby radius) and the chosen cut off frequency of the ~z-coordinate, associated vertical motions will be treated in a Lagrangian way, thus reducing the spurious diapycnal mixing. Coupled ocean–biogeochemical models could also take advantage of the ~z-coordinate properties as biogeochemical tracer structures present very strong vertical gradients. On the longer term, the ~z-coordinate could also have interesting applications in global coupled climate models, diapycnal mixing being a key issue as it sustains the thermohaline circulation. Although coarse climate model resolutions presumably do not give sufficient temporal variability to benefit much from the ~z-coordinate, non negligible effects should appear on very long term simulations. Acknowledgments With thank the NEMO team at the Labratoire d’Océanographie et du climat: Expérimentations et Approches Numériques – Insitut Pierre Simon Laplace (LOCEAN-IPSL) and the Ocean Modelling and Forecasting (OMF) staff at the National Oceanography Centre Southampton (NOCS) for fruitful discussions and comments. We also thank the reviewers and the editor for their very useful comments and the energy spent in long and pertinent reviews. Funding for this study comes from MyOcean and InfraStructure for the European Network for the Earth System Modelling (IS-ENES) projects. Computational resources were provided by the Institut du Développement et des Ressourses en Informatique Scientifique (IDRIS).

8 ~ ðK1  1Þ > < ixK1 ¼ x ~ k2 i1x K2 K2 ¼ ð1  K1 Þ þ xz i1x K2 þ j > 2 : ~ k i1x K2 K3 ¼ 1 þ K2 þ j 8 K1 ¼ 1þi1 x > > ~ x > < 1  K 2 ¼ jk2 () ð1ixx~ Þ 1ixz þ~ x > > > ~ k2 Þixðx ~ þxz Þ ~ ðxz þj x : K3 ¼ ~ þx þj ~ ðx þ j ~ k2 Þþx ~ k2 Þ x2 ixðx z

z

~ ¼ 2s~p and xz ¼ 2sp. Using the same formalism as in the elecwhere x z tronics domain, we define the gain in decibel as Gdb i ¼ 20log10 ðjKi jÞ and the phase as Ui = Arg(Ki). The Bode diagrams plotted in Fig. 2 are the combination of the Gdb i and Ui plots as a function of x with a log-frequency axis. The asymptotic behaviour of Gdb and Ui as i x ? 0 and x ? 1 shows the filtering properties of the ~z-coordinate. The Bode diagram of K1 shows that D0LF is the result of a classical first order low-pass filter applied to D0 : the asymptotic gain slope and phase when x ? 1 are 20 db per decade and  P4 . Its cut pulsation is obtained when the asymptotic straight crosses the ~. pulsation axis in the amplitude diagram which yields xc1 ¼ x Compared to a pure isopycnal model, which would correspond to K2 = 1, f e3 is a obtained by a second order high-pass filter (40 db per decade slope and P2 phase when x ? 0), whose cutoff pulsation qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ ðxz þ j ~ k2 Þ. Eventually the Bode diagram of K3, which is xc2 ¼ x would be 1 in a pure z-coordinate, shows a first order low-pass filter ~ þ xz . for @ kw. Its cut pulsation is xc3 ¼ x Appendix B. Baroclinicity of ee3 In this appendix, we show how Eq. (10) ensures the baroclinic Rk ~ does E3 ðx; yÞ ¼ kbs f e3 ðx; y; kÞ dk ¼ 0, as soon as j nature of f e3 , i.e., f not depend on k and the thickness diffusion fluxes are set to zero through the domain boundaries. To that purpose, we first form the evolution equation of f E3 by integrating (10) over the vertical:

Z

ks

e dk  2p f ~ rh f D E þ rh  ½j E3  sz 3 |fflfflfflfflfflffl{zfflfflfflfflfflffl}

E3 ¼  @t f

Appendix A. Spectral analysis of the ~ z-coordinate

ðB:1Þ

kb

In order to investigate the ability of the proposed ALE coordinate to discriminate properly between low and high frequency motions (corresponding, respectively, to the Eulerian and Lagrangian part of the vertical coordinate) we will perform a linear spectral analysis. Let us imagine that the model is forced by a horizontal baroclinic flow divergence corresponding to a plane wave propagating in the x direction. Barotropic motions are not considered here (D⁄ = 0) because the ~z-coordinate treats them as the z⁄-coordinate, i.e., in a full Lagrangian way. Let us introduce the complex counterparts of D0 , D0LF , @ t f e3 and @ kw. For simplicity they will still be named by their original names:

8 0 D ¼ D0 eiðxtkxÞ > > > 0 > < D ¼ K1 D0 eiðxtkxÞ LF 0 iðxtkxÞ >@ f > t e3 ¼ K2 D e > > : @ k w ¼ K3 D0 eiðxtkxÞ

by their complex counterparts leads to the following linear system and its solution:

¼0

~ is outside of the vertical integral since it does not depend where j on k. Then we form the associated ‘‘quadratic’’ equation by multiplying (B.1) by f E3 :

1 f2 2p f 2 ~ rh f E þ rh  ½j E3  f @ t E3 ¼  E3 2 sz 3

ðB:2Þ

Finally, (B.2) is integrated over the whole horizontal domain we call here S:

1 @t 2

Z Z

S

Z Z 2p f f E3 2 dx dy ¼  E3 2 dx dy sz S Z Z ~ rh f E3  f rh  ½j E3 dx dy þ

ðB:3Þ

S

ðA:1Þ

where x and k are the pulsation and wave number of the plane wave and Ki are the non-dimensional complex amplification factors. The goal of this analysis is to express each Ki as a function of x and k. For this purpose, we will use the system formed by equations (7), (10) and (14). Replacement of the initial variables

Since the thickness diffusion fluxes are set to zero through the domain boundaries, integrating (B.3) by parts leads to:

1 @t 2

Z Z S

Z Z I 2p f f f ~r f E3 2 dx dy ¼  E3 2 dx dy þ E3 ðj E3 Þ  n dl |fflfflfflfflfflfflfflfflhffl{zfflfflfflfflfflfflfflffl ffl} sz ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl S @S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} ¼0 Z Z

60

~ rh f E3 Þ2 dx dy ðj  S |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 60

ðB:4Þ

152

M. Leclair, G. Madec / Ocean Modelling 37 (2011) 139–152

where @S is the boundary of the domain S, dl is the integration element on @S and n is the normal outgoing unit vector. Equation (B.4) RR f shows that E 2 dx dy is decreasing with time. As it is initially S 3 zero and necessarily positive, it stays equal to zero. Then f E3 2 being positive, a zero integral implies that 8ðx; yÞ 2 S; f E3 2 ðx; yÞ ¼ 0, so that Rk e3 ðx; y; kÞ dk ¼ 0. 8ðx; yÞ 2 S; f E3 ðx; yÞ ¼ kbs f References Adcroft, A., Campin, J.-M., 2004. Rescaled height coordinates for accurate representation of free-surface flows in ocean circulation models. Ocean Modelling 7, 269–284. Adcroft, A., Hallberg, R., 2006. On methods for solving the oceanic equations of motion in generalized vertical coordinates. Ocean Modelling 11 (1–2), 224–233. . Adcroft, A., Hill, C., Marshall, J., 1997. Representation of topography by shaved cells in a height coordinate ocean model. Monthly Weather Review 125 (9), 2293–2315. . Adcroft, A., Hallberg, R., Harrison, M., 2008. A finite volume discretization of the pressure gradient force using analytic integration. Ocean Modelling 22, 106– 113. Barnier, B., Madec, G., Penduff, T., Molines, J.-M., Tréguier, A.-M., Le Sommer, J., Beckmann, A., Biatoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theeten, S., Maltrud, M., McClean, J., De Cuevas, B., 2009. Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution. Ocean Dynamics 56, 543–567. Bleck, R., 2002. An oceanic general circulation model framed in hybrid isopycniccartesian coordinates. Ocean Modelling 4 (1), 55–88. . Bleck, R., Ruth, C., Hu, D., Smith, L., 1992. Salinity-driven thermocline transient in a wind-and thermohaline-forced isopycnic coordinate model of the north atlantic. Journal of Physical Oceanography 22, 1486–1505. Bode, H.W., 1940. Relations between attenuation and phase in feedback amplifer design. Bell System Technical Journal 19 (3), 421–454. Burchard, H., Beckers, J.-M., 2004. Non-uniform adaptive vertical grids in onedimensional numerical ocean models. Ocean Modelling 6 (1), 51–81. . Chassignet, E.P., Smith, L.T., Halliwell, G.R., 2003. North atlantic simulations with the hybrid coordinate ocean model (HYCOM): impact of the vertical coordinate choice, reference pressure, and thermobaricity. Journal of Physical Oceanography 33, 2504–2526. Gent, P., McWilliams, J., 1990. Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography 20, 150–155. Griffies, S.M., Böning, C., Bryan, F.O., Chassignet, E.P., Gerdes, R., Hasumi, H., Hirst, A., Treguier, A.-M., Webb, D., 2000a. Developments in ocean climate modelling. Ocean Modelling 2 (3–4), 123–192. . Griffies, S.M., Pacanowski, R.C., Hallberg, R.W., 2000b. Spurious diapycnal mixing associated with advection in a z-coordinate ocean models. Monthly Weather Review 128, 538–564. Griffies, S.M., Gnanadesikan, A., Dixon, K.W., Dunne, J.P., Harrison, R.G.M.J., Rosati, A., Russell, J.L., Samuels, B.L., Spelman, M.J., Wintonand, M., Zhang, R., 2005. Formulation of an ocean model for global climate simulations. Ocean Science 1, 45–79. Halliwell, G.R., 2004. Evaluation of vertical coordinate and vertical mixing algorithms in the hybrid-coordinate ocean model (HYCOM). Ocean Modelling 7 (3–4), 285–322. . Hirt, C.W., Amsden, A.A., Cook, J.L., 1974. An arbitrary Lagrangian–Eulerian computing method for all flow speeds. Journal of Computational Physics 14 (3), 227–253. . Hofmeister, R., Burchard, H., Beckers, J.-M., 2010. Non-uniform adaptive vertical grids for 3D numerical ocean models. Ocean Modelling 33, 70–86. .

Kasahara, A., 1974. Various vertical coordinate systems used for numerical weather prediction. Monthly Weather Review 102 (7), 509–522. . Komori, N., Ohfuchi, W., Taguchi, B., Sasaki, H., Klein, P., 2008. Deep ocean inertiagravity waves simulated in a high-resolution global coupled atmosphere – ocean gcm. Geophysical Research Letters 35, L04610. Leclair, M., Madec, G., 2009. A conservative leapfrog time stepping method. Ocean Modelling 30 (2–3), 88–94. Legg, S., Ezer, T., Jackson, L., Briegleb, B., Danabasoglu, G., Large, W., Wu, W., Chang, Y., Özgökmen, T.M., Peters, H., Xu, X., Chassignet, E.P., Gordon, A.L., Griffies, S., Hallberg, R., Price, J., Riemenschneider, U., Yang, J., 2009. Improving oceanic overflow representation in climate models: the gravity current entrainment climate process team. Bulletin of the American Meteorological Society 90(5), 657–670. . Le Sommer, J., Penduff, T., Theeten, S., Madec, G., Barnier, B., 2009. How momentum advection schemes influence current–topography interactions at eddy permitting resolution. Ocean Modelling 29, 1–14. Lévy, M., Estublier, A., Madec, G., 2001. Choice of an advection scheme for biogeochemical models. Geophysical Research Letters 28, 3725–3728. Lévy, M., Klein, P., Tréguier, A.-M., Iovino, D., Madec, G., Masson, S., Takahashi, K., 2010. Modifications of gyre circulation by sub-mesoscale physics. Ocean Modelling 34 (1–2), 1–15. . Madec, G., 2008. NEMO ocean engine. Note du Pôle de modélisation, Institut PierreSimon Laplace (IPSL), France, No. 27, ISSN: 1288-1619. . Madec, G., Delecluse, P., Crepon, M., Chartier, M., 1991. A three-dimensional numerical study of deep-water formation in the northwestern Mediterranean Sea. Journal of Physical Oceanography 21, 1349–1371. Madec, G., Lott, F., Delecluse, P., Crepon, M., 1996. Large-scale preconditioning of deep-water formation in the northwestern Mediterranean Sea. Journal of Physical Oceanography 26, 1393–1408. Marchesiello, P., Debreu, L., Couvelard, X., 2009. Spurious diapycnal mixing in terrain-following coordinate models: the problem and a solution. Ocean Modelling 26, 156–169. Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., Pairaud, I., Ulses, C., 2008. Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modelling 20, 61–89. Marsaleix, P., Auclair, F., Estournel, C., 2009. Low-order pressure gradient schemes in sigma coordinate models: the seamount test revisited. Ocean Modelling 30 (2–3), 169–177. . Marti, O., Madec, G., Delecluse, P., 1992. Comment on net diffusivity in ocean general circulation models with nonuniform grids by F.L. Yin and I.Y. Fung. Journal of Geophysical Research 97, 763–766. Pacanowski, R.C., Gnanadesikan, A., 1998. Transient response in a z-level ocean model that resolves topography with partial cells. Monthly Weather Review 126, 3248–3270. Penduff, T., Le Sommer, J., Barnier, B., Tréguier, A.-M., Molines, J.-M., Madec, G., 2007. Influence of numerical schemes on current–topography interactions in 1/ 4 global ocean simulations. Ocean Science 3, 509–524. . Shchepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. Journal of Geophysical Research 108, 1–34. Stacey, M.W., Pond, S., Nowak, Z.P., 1995. A numerical model of the circulation in Knight Inlet, British Columbia, Canada. Journal of Physical Oceanography 25 (6), 1037–1062. Timmermann, R., Goosse, H., Madec, G., Fichefet, T., Ethe, C., Dulière, V., 2005. On the representation of high latitude processes in the ORCA-LIM global coupled sea ice ocean model. Ocean Modelling 8, 175–201. Tréguier, A.-M., Dukowicz, J.K., Bryan, K., 1996. Properties of nonuniform grids used in ocean general circulation models. Journal of Geophysical Research 101, 20877–20881. White, L., Adcroft, A., Hallberg, R., 2009. High-order regridding–remapping schemes for continuous isopycnal and generalized coordinates in ocean models. Journal of Computational Physics 228 (23), 8665–8692. . Yin, F., Fung, I., 1991. Net diffusivity in ocean general circulation models with nonuniform grids. Journal of Geophysical Research 96, 10773–10776.