Reducing spurious diapycnal mixing in ocean circulation models

(1)

Reducing spurious diapycnal mixing in

ocean circulation models

DISSERTATION

of

Margarita Smolentseva

at the University of Bremen and the Alfred Wegener

Institute for Polar and Marine Research

Supervisors: Prof. Dr. Thomas Jung, Prof. Dr. Sergey Danilov

Reviewers: Prof. Dr. Thomas Jung, Prof. Dr. Jo

̈rn Behrens

Date of submission: 23.08.2020

(2)

(3)

1

Abstract

Spurious diapycnal mixing of water masses occurs in ocean circulation models as an artifact of numerical algorithms used to advect tempera-ture and salinity. Most of the ocean models used in climate research are based on geopotential vertical coordinates (z- coordinates), which inter-sect isopycnal surfaces. The non-alignment of coordinate surfaces with isopycnals causes spurious diapycnal mixing during horizontal advection of a water-parcel by high-order upwind transport schemes. The growth in the potential energy of the system appears without any sources. This behavior is physically incorrect and leads to an energetic inconsistency and incorrect water mass transformation. Therefore, spurious diapycnal mixing in ocean models is one of the reasons that lead to the incorrect hydrological state of the ocean basins after some integration time. Im-provements are required which would reduce spurious mixing in ocean models.

There are several ways that can potentially reduce numerical mixing and spurious diapycnal mixing in particular. Three of them are considered in the current work. The first way is the design of more accurate advection schemes with the intention to achieve a reduction in truncation error which leads to a decrease in numerical mixing in a system. The second option is the stabilization of central high-order advection schemes by isoneutral diffusion. And the last approach is a choice of the right mesh. It is a question to investigate whether meshes can cause spurious mixing due to their structure, regularity, and other properties.

The current work deals with the problem of spurious diapycnal mix-ing. It analyses the stability of numerical implementation of isoneutral diffusivity on triangular meshes of FESOM2. It proposes a new compact advection scheme characterized by a reduced truncation error compared to other finite volume schemes in FESOM2. It shows that the application of isoneutral diffusion to stabilize central schemes can reduce spurious di-apycnal mixing in models, however, it requires special tuning for every initial state of a model. It is also found out that mesh irregularity does not necessarily imply an enhanced numerical mixing in a system, however, it might depend on the type of triangles.

(4)

Chapter 1 Introduction

A significant part of our knowledge about the variability of the Earth climate relies on numerical modeling. Numerical climate models help researchers to understand the climate of the past and make projections for the future. Ocean circulation models are one of the main components of Earth climate models. Ocean significantly contributes to the meridional heat transport, and ocean-atmosphere interactions provide boundary conditions for atmospheric models over a large part of the Earth’s surface. Ocean transport processes are very slow. It takes centuries for water masses to complete the whole circulation cycle. Ocean has a huge impact on climate, being a part of the planetary energy cycle (Rahmstorf [2002]), and on biogeochemical cycles, exchanging gases with the atmosphere. It has a higher heat capacity than atmosphere and land, influencing climate fluctuations on daily, seasonal, annual, and interannual time scales (Clark et al. [2002]).

Despite ongoing efforts, the existing climate models demonstrate a wide range in their projections concerning the behavior of climate trajectories under differ-ent carbon dioxide increase scenarios (e.g., see Community [2020], Vivek and et [2019], Karen and Richard [2007]). Moreover, model outputs demonstrate quite large biases from historical observations which also bring particular uncer-tainty into their projections. Many possible reasons for such results have been investigated, and a significant part of these reasons arises from numerical errors. Numerous modeling groups over the world make significant efforts to improve the results of the models. Nevertheless, climate models, and particularly global ocean circulation models, still suffer from insufficient resolution and inaccuracies coming from parameterizations and numerical assumptions. With the growth of computer facilities, finer computational meshes become affordable, allowing

(6)

4 CHAPTER 1. INTRODUCTION

ter representation of dynamical processes. However, currently, even the highest resolution that can be used in models for climate research (about 1/10 degree) is insufficient for the representation of many processes. Therefore, these processes are parameterized. The imperfection of numerical methods leads to a lack of consistency in simulated balances and energy pathways.

One of the problems that researchers face is spurious numerical mixing ac-companying the simulation of ocean motion of the ocean water masses. Unlike the atmosphere, the deep ocean is nearly adiabatic, with small background mix-ing on the level of 10−5 m2/s provided largely by internal-gravity waves through their nonlinear evolution and final breaking. The spurious mixing caused by the implementation of numerical operators in ocean models can easily reach values comparable to background mixing or even exceed them as estimated by Ilicak et al. [2012] and Megann [2018]. As explained below, spurious numerical mixing generally occurs when a fluid parcel is advected in ocean models. For numerical stability or the need for positive-definite solutions for the transported quantities, advective operators are, as a rule, equipped with build-in numerical dissipation that depends on the detail of the implementation. Spurious mixing of tempera-ture and salinity generally implies also spurious mixing of density. The biggest concern is due to a diapycnal (occurring across isopycnals) component of mixing which leads to a direct rearrangement of isopycnal surfaces. Spurious isopyc-nal (along isopycisopyc-nal surfaces) mixing also has consequences and may affect the ocean density structure through nonlinearities of the ocean water equation of state (cabeling and thermobaricity).

Mixing in the ocean is the major source of the deep water masses transfor-mation. For example, it is well known that the lowest part of the meridional overturning circulation is driven by mixing that is balanced by slow upwelling of the deep waters (e.g., Kuhlbrodt et al. [2007]). The fact that numerical mixing in ocean models can be as high or even sometimes higher than the physical mixing implies that spurious mixing can strongly affect water mass transformations on long time scales. That might be one of the reasons for biases developing in these models over the long integration time. Furthermore, mixing raises the center of mass and potential energy of the water, so, it cannot happen without exter-nal sources of energy. Background mixing parameterizations such as IDEMIX (Olbers and Eden [2013]) aims to take into account this kind of sources as an important part of the wider idea of energy-consistent modeling. The presence of spurious mixing also means the loss of energy balance.

(7)

5

sensitive to the intensity of simulated eddy motions (to the grid-scale Reynolds number) as it was found by Ilicak et al. [2012]. Nowadays resolution in modern climate models is becoming finer, which allows scientists to perform simulations that resolve mesoscale eddy motions. Mesoscale eddies are an indispensable component of ocean dynamics responsible for vertical and horizontal exchanges in the ocean. Resolving them also means that the simulated eddy kinetic energy becomes much larger, which also might imply spurious mixing in the deep ocean. The existing analyses for 1/4 degree global ocean models show that it can indeed be the case (see, e.g., Ilicak et al. [2012]) but can be controlled in other cases (Gibson et al. [2018]). It remains to be seen how this behavior will be modified with finer resolution and in longer simulations.

For the reasons mentioned above, there is a considerable amount of research aimed at diagnosing spurious mixing and exploring measures to minimize it. Numerous papers devoted to the analysis of water mass transformation (see e.g. Xu et al. [2018] and references therein) concentrate on estimates of diapycnal ve-locities which are produced by both physical and spurious mixing. The difficulty here lies in the need to compare the diagnosed total diapycnal transport with the transport which we would have if all mixing were physical. This can be done in a simplified way (e.g., zonally averaged), as demonstrated by Lee et al. [2002] and Megann [2018], but is difficult otherwise. Nevertheless, Xu et al. [2018] pre-sented an example of situations where the presence of spurious mixing is clearly identifiable. Klingbeil et al. [2014], based on Burchard and Rennau [2008] pro-pose an approach that estimates the discrete variance decay associated with the implemented advective operator. These methods provide three-dimensional maps of spurious mixing, however, they do not distinguish between dia- and isopycnal directions. Approaches for estimating effective spurious diapycnal dif-fusivities were proposed in Griffies et al. [2000] based on adiabatic sorting of water parcels and in Getzlaf et al. [2010] and Hill et al. [2012] based on releasing passive tracers. The method of Griffies et al. [2000], similarly to Megann [2018], provides a basin-averaged view on spurious diapycnal mixing, which is not nec-essarily related to the mixing caused by numerically simulated water parcels (see Ilicak et al. [2012]). Finally, Ilicak et al. [2012] suggested to consider just the increase in the so-called reference potential energy (RPE) proposed by Winters et al. [1995] and also used in Griffies et al. [2000] as a part of their analysis. Since then, the concept of RPE is used by many model development teams to identify their spurious mixing (see Ilicak et al. [2012], Petersen et al. [2015], Mohammadi-Aragh et al. [2015], Gibson et al. [2018], Kärnä et al. [2018]). Its

(8)

advantage is the ease of diagnostics. Once again, it is a global measure, charac-terizing the entire basin. Ilicak [2016] proposes an extension of this concept to quantify the spatial distribution of mixing. The concept of RPE will be used in this work and will be explained in detail in chapter 2.

Ocean modelers are mostly interested in the reduction of diapycnal spurious mixing because it affects the density structure of the ocean and, thus, the struc-ture of the entire ocean circulation. It is also well understood that motions in nearly adiabatic deep ocean happen along isopycnal surfaces: The exchange of two water parcels along isopycnals does not change gravitational potential en-ergy and thus does not require extra sources and sinks. One of the main reasons why the density structure of the deep ocean is not preserved by numerical ocean dynamics is that their vertical discretization often does not follow isopycnals coordinate. Most models used for global ocean simulations use geopotential or

z-coordinates in vertical. Although isopycnals have a very small slope angle

across wide regions of the ocean, they still cross z-model level surfaces, and the angle of crossing can be substantial in frontal areas, e.g. nearby shore in polar regions. Numerical operators in z-coordinate models are implemented in horizontal and vertical directions. They are causing numerical errors. There are two issues, one of which is related to diffusion. One can think of physical diffusion as a combination of mixing along isopycnals with much smaller mixing across isopycnals. Because of this difference, even though isopycnals slope at small angles to the horizontal z-surfaces outside the mixed layer, the approach of splitting mixing into vertical and horizontal components leads to a spurious cross-isopycnal component which is not negligible. This problem was recognized long ago, and Redi [1982] proposed to introduce a diffusivity tensor, which ap-proaches diffusion in such a way that it occurs along isopycnals with isopycnal diffusivity Ki and across isopycnals with much smaller diffusivity Kd defined

by physical parameterizations. Numerical implementation of these (rotated or Redi) diffusivities faced certain problems before Griffies et al. [1998] proposed a numerical solution based on the variational principle. A recent discussion of related numerics and stability analysis for quadrilateral meshes has been pre-sented in Lemarié et al. [2012]. The appearance of models working on triangular meshes requires to reconsider these analyses. For Finite volumE Sea-ice Ocean Model (FESOM) the variational implementation was described in Danilov et al. [2017], however, its stability remained unexplored. Analysis of the stability is provided in chapter 3.

(9)

7

tracers follows a high-order upwind scheme, it, as a rule, causes spurious mix-ing. Upwind schemes are widely used in ocean circulation models, for example, Regional Ocean Modeling System (ROMS, see Moore et al. [2011]), Nucleus for European Modelling of the Ocean (NEMO, see Madec and Team [2016]), MIT General Circulation Model (MITgcm, see Marshall et al. [1997]). Indeed, in this case, it can be shown that the numerical advection operator has the main truncation error of diffusive type (see, e.g., Lemarié et al. [2012]). To stabi-lize central, or so-called mixed (combination of central and upwind), schemes, additional methods such as slope or flux limiters are used. In particular, the flux corrected transport method described by Zalesak [1978] is one of them. Limiters are applied to ensure that solutions stay positive (in case of concentra-tions) or monotone. These methods introduce either local mixing or unmixing (Mohammadi-Aragh et al. [2015]). They are utilized in ocean circulation mod-els such as Finite volumE Sea-ice Ocean Model (FESOM2, see Danilov et al. [2017]) or Model for Prediction Across Scales-Ocean (MPAS-O, see Ringler et al. [2013]). In all cases this mixing will have a cross-isopycnal component which is fully spurious because continuous operators of advection cannot change tracer variance, i.e., create mixing. Numerical advection operators are believed to be the major source of spurious mixing in general circulation ocean models with

z-coordinate models. Therefore, the reduction of spurious mixing is an

impor-tant problem in numerical ocean modeling, which becomes even more significant with the increase of model resolution. There are several ways of handling this problem. One is the design of more accurate advection schemes because a re-duction in truncation error means the rere-duction of spurious effects due to the dissipativity of this error (Hill et al. [2012], Mohammadi-Aragh et al. [2015]). The other approach is the use of high-order centered schemes stabilized with isopycnal (rotated) diffusion, as advised by Lemarié et al. [2012]. A new advec-tion scheme with a reduced truncaadvec-tion error is proposed and analyzed compared to other schemes in chapter 4. The application of isoneutral diffusion is analyzed in chapter 5.

The last topic considered in this thesis relates to meshes. Regular meshes and among them the meshes composed of equilateral triangles are expected to be associated with smaller numerical approximation errors than irregular or distorted meshes. However, variable-resolution meshes are always distorted to some degree. Thus, there is a question if the lack of regularity and structure of meshes can influence numerical mixing in the system and in which way. This is also investigated in the last section of chapter 5.

(10)

All the methods described in this thesis were implemented in FESOM version 2. A concise outlook of the model is provided in chapter 2.

(11)

Chapter 2 Methods

In this chapter, I give an overview of the required basis which the current work is founded on. I describe methods used for reducing spurious mixing and its analysis by researchers. Also, a short description of the model FESOM version 2, used in this work, is provided.

2.1 FESOM2

FESOM2 is a finite-volume sea ice-ocean circulation model. In the model, the equations of motion, continuity, and tracer balance are integrated with respect to time. The main principles of FESOM2 are described in this subsection. More detailed information about FESOM2 can be found in Danilov et al. [2017].

FESOM2 is based on unstructured triangular meshes of variable resolution. The mesh is defined by a grid of triangles in the horizontal directions and a system of horizontal levels in the vertical direction which partition the com-putational domain into triangular prisms. In the horizontal plane, the scalar discrete variables of FESOM2 are placed at vertices of triangles, discrete hor-izontal velocities are placed at centroids of triangles. A control volume for velocities is defined by a triangular prism in the centroid of which this velocity value is defined. Control volumes of scalars are defined by median-dual prisms. A median-dual prism for a scalar T is obtained by linking centroids of all the triangles that have the vertex of T through the mid-edge points. Figure 2.2 shows schematically a control volume for velocities and scalars.

In the vertical direction, the horizontal velocities and scalars are placed in the mid-layers, and the vertical velocity is at full levels. Horizontal gradients of scalar quantities are located at the centers of triangles, and vertical gradients

(12)

10 CHAPTER 2. METHODS

vertex shaded with light blue colors. The control volume for velocity in the horizontal plane is represented by a triangle (indicated by light red colors for

v5). The light purple color is the area where both these control volumes intersect.

are at the vertices and full levels (see Fig.2.3). A triangular prism of the mesh can therefore be split into six triangular sub-prisms with unique combinations of discrete values of scalars and their derivatives (see Fig. 2.4). We will need these sub-prisms to explain the implementation of the isoneutral diffusion operator.

(13)

3j+1 ∂T ∂z

2j+1 ∂T ∂z

1j ∂T ∂z

3j ∂T ∂z

2j

Figure 2.3: A triangular prism element of the mesh. The scalar values are located in the middle of the prism edges, the horizontal gradient of the scalars is placed in the middle of the prism, and their vertical derivatives are in the prism vertices. Six sub-prisms obtained by splitting the triangular prism by the mid-plane (drawn) and faces of median-dual scalar control volumes (not shown) are characterized by triples of scalar value, scalar horizontal gradient, and scalar vertical gradient.

Figure 2.4: A triangular prism divided into six sub-prisms.

Governing equations

Let us introduce vertical layer thicknesses hk = hk(x, y, z) where k changes from

1 to K, and K is the number of layers. Then the continuity equation integrated over the vertical extent of the layer is written as

(14)

Here u is horizontal velocity, W is the water flux leaving the ocean at the surface,

wt

k and wbk denote transport velocities through the top (the superscript t) and

bottom (the superscript b) boundaries of layer k, and ∇h = (∂x, ∂y). Let us

denote an arbitrary tracer as T . The layer-integrated equation for the tracer will be written as

∂t(hT )k+ ∇h· (uhT )k+ (wtTt− wbTb)k+ W TWδk = ∇ · hkK∇Tk. (2.2)

In this equation ∇ = (∂x, ∂y, ∂z) denotes the 3-dimensional operator, i.e., the

right hand side in (2.2) is the divergence of the flux due to 3-by-3 diffusivity tensor K which will be described in details in this chapter. The equation for the surface elevation η (also called the sea surface height (SSH)) is obtained by summing (2.1) over layers

∂tη + ∇h· ∑︂

k

hkuk+ W = 0. (2.3)

The pressure field in hydrostatic approximation is expressed as

p = pa+ gρ0η + ph, ph = ∫︂ 0

z

ρdz (2.4)

where pa is the atmospheric pressure and ρ the deviation of density from its

reference value ρ0, and ph is the hydrostatic pressure due to density variation ρ.

The momentum equation is taken directly, without layer integration, and the layer index k is omitted. It is written below in a vector-invariant form:

∂tu + ω + f h k × uh + (︂ (w∂zu)t+ (w∂zu)b )︂ /2 + ∇h (︂ p/ρ0+ u2/2 )︂ + gρ∇hZ/ρ0 = Duu + ∂z(ν∂zu), (2.5) but also a flux form is available. In this equation gρ∇hZ reflects the fact that

the model layers may deviate from geopotential surfaces. The variable Z stands for the vertical coordinate z of the midplane of the layer with the thickness h. The second term of the equation represents the potential vorticity q = (ω +f )/h.

Finally, the set of equations is completed by the equation of state,

ρ = ρ(T, S, z),

with T the potential temperature, S the salinity, and x the depth which replaces the pressure because of the Boussinesq approximation used in the equations above.

(15)

2.1. FESOM2 13

ALE

FESOM2 provides three options for the treatment of free-surface. One of them is the linear free surface where the change in the thickness of the upper layer due to surface elevation is ignored in the tracer and momentum equations so that the vertical layers are kept fixed in time, and volumes of the prisms stay constant throughout the simulations. Two other options represent the full free surface and use arbitrary Lagrangian-Eulerian (ALE) vertical coordinate (see Donea and Huerta [2003]) to account for the moving surface, which warrants full conservation. The full free surface is represented in FESOM2 with two options: zlevel and zstar. With the zlevel option, only the thickness of the upper layer is changing with time following the variations in SSH. With the zstar option, the thicknesses of all the layers are dynamically updated and the change in the thickness of the fluid column due to the elevation is equally distributed among all the vertical layers. The zlevel option is slightly more computationally efficient since only the thickness of the upper layers is updated. However, this option may impose limitations when the unperturbed thickness of the upper layer is too small: the thickness of the upper layer may become too close to zero, or even negative, which is not allowed. For more details see Scholz et al. [2019], Danilov et al. [2017]. The ALE and full free surface options allow one to get rid of so called virtual salinity fluxes. When the freshwater is added to the surface, the thickness of the upper layer is modified according to the thickness equation 2.1 and the salinity in the upper layer is controlled by that change. There is zero flux of salt, and total salt is conserved. The heat flux is taken into account through the boundary conditions at the surface in the temperature equation 2.2.

Spatial discretization

First, let us consider the notation. The indices c, e and v enumerate cells, edges, and vertices respectively. The notation like v(e) will be used to denote the set of vertices of edge e. The notation e(v) means the set of all the edges coming out of the vertex v. This rule is applied to all other possible sets c, e and v.

In order to discretize the governing equations in terms of finite volumes, these equations are integrated over the control volumes. For the velocities it means

Acuc= ∫︂

c

udS,

and similarly for the scalar tracers

AkvTkv = ∫︂

kv

(16)

Here Acis a cell area, and S is a contour. Akv is a horizontal area of a prism for

the vertex v at the layer k. The layer index k is suppressed everywhere where this causes no ambiguity.

There are three options for the momentum equation in FESOM2. Two of them have a flux form and another one has a vector-invariant form. Cell-vertex discretization on triangular meshes has an excessive number of velocity degrees of freedom (compared to scalar degrees of freedom) which was taken into account in these options. Implementation of the momentum equation requires a certain amount of averaging in order to avoid noise which can occur on a grid scale.

• Vertex velocity option. Velocity at a vertex is calculated by averaging among the velocities defined at the cells which contain the vertex v

Avuvhv = ∑︂

c(v)

uchcAc/3.

Here c(v) is a set of all the prisms containing the vertex v, hc are prism

thickness. c(v) coincides with c(v) in the upper layers, but bottom to-pography excludes some prisms in deeper layers. Using this averaging the diversion of horizontal momentum flux is calculated

Ac (︂ ∇ · (huu))︂ c= ∑︂ e(c) le (︄ ∑︂ v(e) ne· uvhv )︄(︄ ∑︂ v(e) uv/4 )︄ .

Here ne is an external normal at edge e of cell c, le is the length of the

edge e.

• Scalar control volumes. The horizontal part of the momentum flux diver-gence on scalar control volumes is computed as

Av (︂ ∇ · (huu))︂ v = ∑︂ e(v) ∑︂ c(v) uchc· necucdec.

For the vertical part, the flux going through the top boundary of the prism in layer k is

Av(wvut) = wv ∑︂

c(v)

ut_cAc/3.

Here the superscript t indicate the top surface.

• Vector-invariant form. The Coriolis parameter (see (2.5)) has to be defined on the vertices of a prism because the relative vorticity is defined there

(︂

(ω + f )k × (u))︂=∑︂

v(c)

(ω + f )vk × uc/3.

The kinetic energy in (2.5) is first computed at scalar locations, and then its horizontal gradient is computed at triangles.

(17)

2.1. FESOM2 15

Transport equation is discretized with high order schemes are used. FESOM2 uses upwind (3rd order), central (4th order), and mixed schemes of 3rd-4th order. The schemes are based on the estimates of the tracer gradients on an upwind or central stencil extending beyond the nearest neighbors. This and other advection schemes are described in more detail and analyzed in chapter 4.

The Gent-McWilliams (GM) parametrization of eddy stirring (Gent and McWilliams [1990]) is implemented in FESOM2 using the algorithm proposed by Ferrari et al. [2010]. The GM parametrization is always applied together with isoneutral diffusion proposed by Redi [1982]. The GM parametrization together with isoneutral diffusion is important for parametrizing eddy-scale motions and diapycnal mixing in the ocean. Isoneutral diffusion and its analysis in FESOM2 is one of the key points considered in the current work.

Isoneutral diffusivities

Isoneutral diffusivity directs mixing of temperature, salinity, and other tracers, caused by unresolved eddies, along isoneutral surfaces in the ocean. It can also be used for the stabilization of central advection schemes.

An isopycnal is commonly understood as the surface of constant potential density defined with respect (referenced) to a particular pressure level. The nonlinearity of the equation of state leads to the fact that this surface deviates from a locally referenced density surface, which is also called a neutral density surface and commonly denoted γ, even if the potential density and isoneutral density coincide at some point. Small and almost energetically neutral exchanges of fluid particles are expected to occur along neutral surfaces and be isoneutral.

The neutral slope vector is computed as

∇γ = −α∇T + β∇S.

Here T and S are the potential temperature and salinity respectively, α is the coefficient of thermal expansion and β is the coefficient of haline contraction. Note that the pressure appears in this equation only through the expansion and contraction coefficients, and therefore these coefficients depend on local pressure (depth in the Boussinesq approximation). γ surfaces can be computed only approximately and this is the reason the language of isopycnal surfaces is often used to express the physics approximately. In most cases, the notions of isopycnal and diapycnal directions below have the precise meaning of isoneutral and dianeutral directions respectively.

(18)

Diffusive flux FT of tracer T is commonly parameterized as being

propor-tional to the gradient of the tracer. Since the flux is a vector, a general connection between it and another vector ∇T is a diffusivity tensor given by 3 by 3 matrix

K in the standard Cartesian representation, FT = −K∇T.

Here the minus sign takes care that the flux is down gradient in the simplest case when K = κI, where I is the identity matrix.

To have the flux vector oriented along the γ surfaces one should write K =

KiP, where Ki is the isoneutral diffusivity and P is a projecting operator that

eliminates components aligned with ∇γ. For the matrix, K this operator is expressed as

P = I − ∇γ∇γ/|∇γ|2,

or in the component form

Pij = δij − ninj, ni = ∂iγ/|∇γ|,

where i, j ∈ {x, y, z}, and δij is the Kronecker delta. Normally a slope vector in

the direction of ∇γ is used instead of the unit vector ni. This vector is defined

as

s = −∇hγ/∂zγ, ∇h = (∂x, ∂y).

The matrix K in terms of the slope vector s

K = Ki 1 + s2 ⎛ ⎜ ⎜ ⎜ ⎝ 1 + s2 y −sxsy sx −sxsy 1 + s2x sy sx sy s2 ⎞ ⎟ ⎟ ⎟ ⎠ , (2.6)

where sx and sy are components of the isoneutral slope vector s, and s is its

magnitude.

We assume that the slope is small which means that the values of the vector

s are also small, |s| ≪ 1. In this case, we can approximate 1+s2 _{by 1 up to small}

errors that are quadratic in small |s|. For the same reason, all quadratic terms in the matrix can be ignored except for the diagonal s2 _{term which corresponds}

to the vertical diffusion. Thus, the expression above is simplified to

K = ⎛ ⎜ ⎜ ⎜ ⎝ Ki 0 sxKi 0 Ki syKi sxKi syKi s2Ki ⎞ ⎟ ⎟ ⎟ ⎠ ,

(19)

2.1. FESOM2 17

which is referred to as small-angle approximation. Finally, we should add the physical dianeutral diffusion Kd that is provided by vertical mixing

parameter-ization schemes. The dianeutral direction is very close to the vertical one, so with sufficient accuracy, this addition affects the lower diagonal term. The over-all form of the diffusivity tensor (in smover-all-angle approximation) will take the form: K = ⎛ ⎜ ⎜ ⎜ ⎝ Ki 0 sxKi 0 Ki syKi sxKi syKi s2Ki+ Kd ⎞ ⎟ ⎟ ⎟ ⎠ .

In the mixed layer, the slope vector is no longer small. However, the fluid motion is diabatic in the mixed layer, and fluid is mixed horizontally across the isoneutral surfaces. For this reason, we taper the off-diagonal components of the diffusivity tensor to zero at the locations in and close to the mixed layer. Thus,

Ki becomes the coefficient of horizontal diffusion in regions where tapering is

done. A practical criterion is the magnitude of the slope vector.

There are two important points to consider while implementing the rotation diffusivity numerically. The first one is that in finite-volume methods the com-ponents of the slope vector are naturally computed at different locations: α and

β are defined at the center of a scalar cell, horizontal and vertical derivatives

are computed, respectively, at lateral and horizontal faces. The slope vector is not defined at a single point unless some averaging is made. However, if this averaging is made, it is not guaranteed that the diffusive operator will lead to variance decay. The rigorous approach requires using variational calculus, i.e., writing a discrete dissipation functional and obtaining discretization by taking variations (differentiating) of this functional, as first proposed by Griffies et al. [1998] (see also discussion in Lemarié et al. [2012] and the explanation below).

The second point is the numerical stability of the time integration. Because of the presence of mixed (horizontal-vertical) derivatives, the explicit time stepping of the isoneutral diffusion operator faces severe restrictions if the aspect ratio of a cell is smaller than the slope. In order to avoid this difficulty, the integration has to be split into explicit and implicit parts. We explored the stability of this split at triangular meshes where scalars are placed at vertices in chapter 3.

The time integration split, tapering, and other numerical details make the isoneutral diffusion not truly isoneutral, although the dianeutral component, introduced by them, is small compared to the horizontal diffusion.

(20)

Analysis of isoneutral diffusion

We consider a diffusion equation

∂T

∂t = ∇K∇T. (2.7)

Here T is a tracer field (temperature, salinity, etc.). This equation has two properties.

• Property 1.

If tracer T is the isoneutral density γ, from definitions above (see equation 2.6)

K∇γ = 0. (2.8)

This happens because K is the projection operator which removes compo-nents aligned with ∇γ.

• Property 2.

We multiply equation (2.7) with T and integrate over volume. The integral on the right-hand side is integrated by parts assuming no flux conditions to obtain ∂t ∫︂ T2/2 dV = ∫︂ T ∂iKij∂jT dV = − ∫︂ ∂iT Kij∂jT dV ≤ 0. (2.9)

The last equality follows from matrix Kij of K being symmetric and

posi-tive if T differs from γ. Here i, j ∈ {x, y, z}.

Let us prove that the matrix K is positive-definite as it is mentioned in the prop-erty 2. For this, let us consider Sylvester’s criteria which says that a symmetric matrix is positive-definite when all its leading principal minors are positive. In-deed, ∆1 = |Kiso| > 0, ∆2 = Kiso2 > 0, and ∆3 = Kiso2 (s2Kiso+ Kd) − s2xKiso3 −

s2

yKiso3 = s2Kiso2 Kd > 0. Thus, the matrix K is positive-definite.

Properties 1 and 2 must be maintained at a discrete level. The difficulty is that in the discrete case T , its horizontal derivatives, and its vertical derivatives lie at different locations and the same concerns the quantities needed to compute the isoneutral gradient vector defining K.

We introduce the dissipation functional F (T ) = −1

2

∫︂

∇T K∇T dΩ. (2.10)

As can be seen, the right-hand side of (2.9) is expressible as 2F (T ). One can also readily see that the right-hand side of (2.7) can be obtained by the calculus

(21)

2.2. VARIANCE DECAY AND TRUNCATION ERRORS 19

of variations as the functional derivative

δF = −(1/2) ∫︂ ∇δT K∇T dΩ − (1/2) ∫︂ ∇T K∇δT dΩ = − ∫︂ ∇δT K∇T dΩ = − ∫︂ ∇(δT K∇T )dΩ + ∫︂ δT (∇K∇T )dΩ = ∫︂ δT (∇K∇T )dΩ (2.11) Here δT is a small variation of tracer field satisfying boundary conditions. The first equality in the chain of transformations follows from the symmetry of K, and the second one is due to our assumption that no flux is coming through boundaries (fluxes coming to the ocean are associated with vertical (dianeutral) diffusivity). The chain of transformations implies that the right-hand side of (2.7) is also a multiplier of δT in the last integrand, i.e., it is the functional derivative of the dissipation functional,

δF

δT = ∇(K∇T ). (2.12)

Thus, starting from negative-definite dissipation functional one arrives at the right-hand side of (2.7). Vice versa, it was shown above that the right-hand side of (2.7) leads to the variance decay.

The idea of Griffies et al. [1998] exploits this fact. To get a numerical im-plementation of isoneutral diffusivity that satisfies Property 2, one first needs to write the dissipation functional (2.10) in the discrete form and obtain the expression for the right-hand side of discretized (2.7) by performing a discrete analog of variation computation. Analysis of isoneutral diffusivity, and dissi-pational functional in particular, on triangular meshes is provided in chapter 3.

2.2 Variance decay and truncation errors

Consider a 1D advection–diffusion equation

∂tT + ∂x(uT ) = ∂xκ∂xT (2.13)

in a domain x ∈ [a, b] assuming periodicity of u and T for simplicity to eliminate the consideration of boundary effect. In the equation above T is a scalar field (tracer), u is the advecting velocity and κ is the diffusivity. Integrating this equation over the domain we easily find that ∂t

∫︁b

aT dx = 0, i.e. the total

(22)

order to see the effect of diffusion we should consider behavior of variance T2_.

To learn about behavior of the variance, we multiply the equation above with T and integrate over the interval. After simple manipulations we find:

∂t ∫︂ b a T2dx = −2 ∫︂ b a κ(∂xT )2dx, (2.14)

which means that the variance decay is related only to the diffusion, and not to the advection. Indeed,∫︁b

a T ∂x(uT ) dx = ∫︁b

a∂x(uT2/2) dx = 0 (in a more common

case this result is a consequence of impermeability of boundaries). Thus in the continuous case advection does not create variance decay, i.e. it does not relate to mixing.

In the discrete case, our approximation of the advection operator always contains truncation errors. Let us discretize the interval [a, b] with cells of the size h. The cells will be indexed with n, and it will be assumed that Tn is

the discrete value of T at the center of cell n. Assume for simplicity that u is constant and positive over the interval. With the first order upwind method, the following approximation is used

u∂xT |n ≈ u(Tn− Tn−1)/h. (2.15)

In this formula, we assume that the velocity u is positive. We use the Taylor expansion for Tn−1 around the center of cell n to obtain

Tn−1 = Tn− h∂xT |n+ (h2/2)∂xxT |n+ O(h3),

i.e.,

u∂xT |n≈ u(Tn− Tn−1)/h = u∂xT |n− (uh/2)∂xxT |n+ O(h2). (2.16)

The leading term in the truncation error is diffusive and will lead to mixing that is rather high (with κA= uh/2).

Acting in the same way we will find that for the second-order central differ-ences

u∂xT |n≈ u(Tn+1− Tn−1)/(2h) = u∂xT |n+ (uh2/6)∂xxxT |n+ O(h4). (2.17)

Here the main truncation error is coming with odd derivative, and the leading error is dispersive. The order of the method here is associated with the power of h in the leading truncation term.

The difference between the two methods lies in both their order and the type of their truncation term. This behavior is preserved for higher-order methods.

(23)

2.2. VARIANCE DECAY AND TRUNCATION ERRORS 21

For example, the standard third-order upwind method has the truncation error (uh3_/12)∂

xxxx (see, e.g., Lemarié et al. [2012]), which is a biharmonic diffusion.

The numerical illustrations for the errors of the fourth-order method are given in chapter 4.

The type of truncation error is important because a dispersive error does not lead to variance decay. Take, for example, the case of second-order centered differences. In this case, the contribution of the error to the variance decay is proportional to ∫︂ b a T ∂xxxT dx = ∫︂ b a ∂x(T ∂xxT )dx − ∫︂ b a ∂x(T ∂xT )dx = 0.

The integrals above take zero values due to the periodic boundary conditions of the flux. One will get the same zero result for higher order even methods. There might be some boundary effects in a general case, but their role is expected to be small if the domain is large.

This is the reason why one expects that even-order methods do not create mixing on their own. If velocity is variable, this statement is not necessarily correct. Nevertheless, it is natural to expect that even-order methods will be nearly non-dissipative. The drawback of even-order methods is that, because of dispersion, the numerical propagation speed of perturbations becomes a func-tion of their wavelength even if u is constant. As a result, oscillafunc-tions will be developing around frontal features, leading finally to code instabilities. An even-order scheme has to be supplemented with an appropriate dissipation which will suppress oscillations. Thus, in the end, both odd-order schemes with build-in dissipation and even-order schemes with added dissipation imply some numeri-cal dissipation, i.e. spurious mixing. The conceptual advantage of an even-order scheme is that dissipation can be treated independently from advection. In the 3D case, dissipation can be turned along isoneutral surfaces to minimize or even avoid dianeutral mixing (Lemarié et al. [2012]).

If the order of scheme is increased, the truncation term becomes asymptoti-cally smaller (assuming that the tracer field is smooth enough), and will either introduce less dissipation or will need less dissipation to counter dispersion. For this reason, methods reducing the magnitude of truncation errors are expected to need less dissipation for stability. Chapter 4 presents a fourth-order compact scheme developed in the framework of this thesis that has a reduced truncation error.

The simple analysis above explains the motivation behind particular parts of this work. There are, however, complications in the practical realization, which

(24)

may partly modify the simple arguments described above. One of the compli-cations is the discrete time stepping. In FESOM2 advection is implemented with second-order Adams–Bashforth method. Using the 1D example above and assuming κ = 0, the procedure of FESOM2 is modeled by

T_nm+1− Tm n = −(uτ /h)(T AB2 n+1 − T AB2 n−1)/2, (2.18)

where m is the time step index, τ is the time step, u = const and

T_nAB2 = (3/2 + ϵ)T_nm− (1/2 + ϵ)Tm−1

n . (2.19)

Here, ϵ is a small parameter (0.1) needed for stability. To proceed, the fields are expanded around (m, n) point in time and space in a Taylor series, and then time derivatives are expressed in terms of space derivatives with account for the order of approximation. As a result, we get a modified equation that shows what is actually solved. The procedure is called a modified equation method and is described in many textbooks (see, e.g., Donea and Huerta [2003]). The result for the case considered is

∂tT + u∂xT = ϵu2τ ∂xx− uh2∂xxxT (1/6 + (uτ /h)2(−5/12 − ϵ + ϵ2)) + O(τ, h)3

(2.20) The time stepping stabilization through ϵ introduces in the leading order diffu-sion with the equivalent diffusivity coefficient ϵu2_{τ , i.e. spurious mixing.}

Spuri-ous mixing can be essential in frontal regions where the ocean velocity can reach 1 m/s. Note also that the amplitude of the dispersion term is modified by time stepping. The combination uτ /h is called the Courant number, which is com-monly small in ocean applications. For example, FESOM2 in global simulations on a 1/4 degree mesh is running with τ = 1200 s, giving uτ /h ≤ 0.1. For this reason, the compensation of dispersive error by time stepping is mostly weak, and errors of spatial discretization are as a rule dominant.

It should be reminded once again that the analysis above is only valid for uniform velocity. In realistic application, velocity is variable and depends on the fields of active tracers such as temperature and salinity. However, it can be expected that the behavior of advection schemes remains qualitatively similar to that described above.

2.3 RPE analysis

To estimate spurious mixing in the model, the concept of reference potential energy (RPE) will be applied in a way how it was described by Winters et al.

(25)

2.3. RPE ANALYSIS 23

[1995]. RPE is the minimum potential energy of the system obtained by sorting density without mixing. RPE changes only through dianeutral transport. When mixing occurs in a model, RPE generally increases with time (Ilicak et al. [2012]). For calculation RPE all the water parcels in a model have to be sorted in a way that the heaviest parcels are located at the bottom and the lightest ones -at the top of a reservoir. The RPE is calcul-ated as an integral by volume over density-weighted geopotential of the new state of a model with sorted parcels:

RP E = g

∫︂ ∫︂ ∫︂

ρ∗zdV. (2.21)

Here ρ∗ stands for density of the state with sorted parcels. It is important to remember that after this operation when a parcel is placed at another depth under different pressure, its density will change. This behavior causes additional problems in the RPE computation for the ocean with the full equation of state, requiring new sorting after each time a parcel for the next depth is sought. However, Ilicak et al. [2012] shows that using an isopycnal density instead of in situ density in this procedure allows one to do sorting only once, for incurring errors are small. Basically speaking, RPE can be considered as unavailable potential energy of the system. It would not change with time in models with switched off explicit dissipation and forcing unless the effects of spurious mixing. When mixing exists in the system, e.g. from dissipation brought by advection schemes, RPE will grow.

Hence in an ocean domain with no boundary fluxes of buoyancy, the evolution of the RPE directly reflects a nonzero dianeutral transport (Ilicak et al. [2012]). To avoid unnecessary information about RPE due to its small changes, I consider relative changes of RPE with respect to initial state:

RP E = RP E(t) − RP E(t = 0)

RP E(t = 0) . (2.22)

Following Ilicak et al. [2012] even after adiabatic sorting we avoid diagnosing local, or isopycnally averaged, diffusivity from the sorted profile due to its noisy character and because it can be potentially misleading. Exploring the potential energy of mixing through the RPE analysis is both consistent and relevant to the global energy of the ocean.

(26)

(27)

Chapter 3 Isoneutral diffusion on triangular meshes

In deep ocean mixing of water mass properties by eddy motion, up to small di-aneutral mixing, occurs along isoneutral surfaces. This is taken into account by the diffusivity tensor K defined above. Numerical implementation of isoneutral diffusivity is a challenge because it has to satisfy several properties mentioned above. The existing discrete implementations were formulated and analyzed for quadrilateral meshes (see, e.g., Griffies et al. [1998], Lemarié et al. [2012]). Addi-tional analyses are required for triangular meshes. The description of discretiza-tion appropriate for triangular meshes of FESOM2 was provided in Danilov et al. [2017]. I worked on the implementation of isoneutral diffusivities in FE-SOM2 which was then tested in Scholz et al. [2019]. This chapter deals with the analysis of limitations on the time step required for numerical stability of isoneutral diffusion. First, it considers a 2D case to explain the procedure, and then continues with the analysis of 3D discretization of FESOM2.

3.1 2D discretization

The intention of this section is to explain computations using a 2D example. A fragment of 2D mesh is shown in Fig. 3.1. Indices k and i enumerate mesh cells in the vertical and horizontal directions respectively. The black dots show the placement of discrete tracer variables. Using the finite volume method, we can discretize equation 2.7 as follows

Vk,i∂tTk,i= ∫︂ k,i ∇K∇T dΩ = ∫︂ k,i (K∇T )ndS. (3.1)

In this equation Vk,i stands for the volume of the cell {k, i} (area in

two-dimensional case); S stands for the surface of the volume Vk,i, and n is the

(28)

26 CHAPTER 3. ISONEUTRAL DIFFUSION ON TRIANGULAR MESHES V_kul_−1,i−1 Vdl k−1,i−1 Vur k−1,i−1 Vdr k−1,i−1 k + 1, i k, i− 1 k, i k, i + 1 k_{− 1, i} h_k₋1 2 i− 1 2 i + 1 2 d_i₋1 2

Figure 3.1: Discretization on a quadratic mesh with the step h in the vertical direction and step d in the horizontal direction. Subvolumes V_k,ij are shown in purple color in the left upper volume.

outer normal vector to the surface. In two-dimensional case the K matrix will have the following form

K = Kiso ⎛ ⎝ 1 s s s2 ⎞ ⎠. (3.2)

Discrete vertical derivatives of the tracer will be at horizontal faces (edges) of the cells, i.e at k + 1/2 and so on. Horizontal derivatives will be at i + 1/2 and so on. Since K combines vertical and horizontal derivatives, fluxes through the cell faces in (3.1) on each face will depend on the derivatives defined at other places. However, there is no immediately straightforward way to average vertical derivatives to vertical faces and horizontal derivatives to the horizontal faces so that it will ensure that the discrete variance decays.

First, let us introduce the discrete form of the dissipation functional F : F = −1 2 ∑︂ k,i ∫︂ k,i ∇T K∇T dΩ. (3.3)

Since horizontal and vertical derivatives are defined at different places, Vk,i will

(29)

3.1. 2D DISCRETIZATION 27

right parts of the volume. Thus, the volume Vk,i will be split into four smaller

volumes each of which takes the fourth part of the volume of Vk,i, denoted as

Vul

k,i, Vk,iur, Vk,idl, Vk,idr as it is shown on the Fig. 3.1.

Using these sub-volumes, we split the integral in (3.3) in four contributions F = −1 2 ∑︂ k,i 4 ∑︂ j=1

(∇T )j_k,iK_k,ij (∇T )j_k,iV_k,ij . (3.4) Here j is array of indexes from 1 to 4 where 1 stands for the combinations ul, 2 for ur, 3 for dl and 4 for dr. Each of four sub-volumes of Vk,i is characterized

by unique combination of tracer Tk,i, and horizontal and vertical derivatives,

producing a positive contribution to the sum in (3.4) because of the positivity of the matrix K_k,ij . The small variation of the tracer now becomes the vector of perturbations of the discrete values δTi,k. Since (3.4) is a quadratic form in Ti,k,

the variation of (3.4) can be computed as

δF =∑︂ k,i δTk,i ∂F ∂Tk,i . (3.5)

On the other hand, in analogy with (2.11), we write

δF =∑︂

k,i

δTk,iRk,iVk,i. (3.6)

Here Rk,i is used to represent the right hand side of discrete diffusion equation

(3.1) as ∫︁

k,i∇K∇T dΩ = Vk,iRk,i. Hence, we will get

Vk,iRk,i =

∂F ∂Tk,i

. (3.7)

Thus, the right hand side of (3.1) is computed by differencing of the discrete dissipation functional. It can be seen that it ensures discrete variance decay (Property 2). Indeed,

∑︂

i,k

Tk,iVk,iRk,i= ∑︂ k,i Tk,i ∂F ∂Tk,i = 2F .

The last equality is once again the consequence of F being a quadratic form in

Ti,k. Property 1 will be ensured if the slope is computed separately in every four

subvolumes by using the same expressions for vertical and horizontal derivatives as for the tracer.

Let us expand the term (∇T )j_k,iK_k,ij (∇T )j_k,i in (3.4) using directional deriva-tives and suppressing indices for briefness

(∇T )j_k,iK_k,ij (∇T )j_k,i= (∂xT, ∂zT )(Kiso(∂xT + s∂zT ), Kiso(s∂xT + s2∂zT )) =

Kiso((∂xT )2+ 2∂xT ∂zT s + s2(∂zT )2).

(30)

28 CHAPTER 3. ISONEUTRAL DIFFUSION ON TRIANGULAR MESHES

Here Kiso stands for isoneutral diffusivity in the matrix K.

Let us have a closer look at the terms of the equation above. We will be writing down the contributions coming from the volume Vk,i into the dissipation

functional (3.4). For simplicity of explanations, we will assume that the slope

s is the same for the entire volume Vk,i. This assumption is compatible with

Property 2 but violates Property 1. Term 1: (∂xT )2Kiso.

The contribution from volume Vki and term 1 differs for the left and right

sub-volumes; all sub-volumes lead to the expression

− [︄ Vki 2 ⎛ ⎝ Tk,i+1− Tk,i di+1/2 )︄2 +Vki 2 (︄ Tk,i− Tk,i−1 di−1/2 )︄2⎤ ⎦ Kiso 2

Here di is the step in x (horizontal) direction. Therefore, the contribution to

the right hand side according to equation 3.7 becomes

Vk,iRk,i: Kiso [︄ 1 di+1/2 Tx,k,i+1/2+ 1 di−1/2 Tx,k,i−1/2 ]︄ 1 2Vk,i,

Vk,i+1Rk,i+1: −Kiso

1

di+1/2

Tx,k,i+1/2Vk,i,

Vk,i−1Rk,i−1: Kiso

1

di+1/2

Txk,i−1/2Vk,i.

In these formulas, Tx is the shortcut for the discrete x-derivative at locations

given by further subscripts. It is important to note that if Kiso is varying as a

function of x, it should be taken in the same way as Vk,i (with indices k, i).

We can proceed by putting the result as written to the right-hand sides of three cells {k, i−1}, {k, i} and {k, i+1}. We can also combine the contributions from the cells {k, i − 1}, {k, i} and {k, i + 1} into the right-hand side of the equation for Tk,i as follows

Vk,iRk,i : KisoVk,i 2 [︄ 1 di+1/2 Tx,k,i+1/2− 1 di−1/2 Tx,k,i−1/2 ]︄ +KisoVk,i+1 2 1 di+1/2 Tx,k,i+1/2− KisoVk,i−1 2 1 di−1/2 Tx,k,i−1/2 = Kiso [︄ Tx,k,i+1/2 2di+1/2 (︂ Vk,i+ Vk,i+1 )︂ −Tx,k,i−1/2 2di−1/2 (︂ Vk,i+ Vk,i−1 )︂ ]︄ . (3.9)

This equation expresses the result in a flux form. We see the flux entering the cell {k, i} through its right and left faces. The fluxes are, however, weighted with volumes of cells on both sides of faces.

(31)

Term2: 2∂xT ∂zT sKiso.

The contribution of the cell {k, i} to the F is: −sVk,iKiso 4 [︄ Tk−1,i− Tk,i hk−1/2 Tk,i− Tk,i−1 di−1/2 +Tk−1,i− Tk,i hk−1/2 Tk,i+1− Tk,i di+1/2 +Tk,i− Tk+1,i hk+1/2 Tk,i− Tk,i−1 di−1/2 + Tk,i− Tk+1,i hk+1/2 Tk,i+1− Tk,i di+1/2 ]︄ .

In this equation hk+1/2, hk−1/2 stand for the distance in the z-direction between

the respective cell centers. Note that here all four sub-volumes contribute dif-ferently. The slope s is a common multiplier because of our simplification. It is computed separately for each sub-volume otherwise. As can be seen, five discrete tracer values are involved, so, that this term will generate five contributions to the right-hand sides. They are computed by taking derivatives

Vk,iRk,i : sKisoVk,i 4 [︄ 1 hk−1/2 Tx,k,i−1/2− 1 di−1/2 Tz,k−1/2,i 1 hk−1/2 Tx,k,i+1/2+ 1 di+1/2 Tz,k−1/2,i − 1 hk+1/2 T, xk,i−1/2− 1 di−1/2 T, zk+1/2,i − 1 hk+1/2 T, xk,i+1/2+ 1 di+1/2 Tz,k+1/2,i ]︄ , Vk−1,iRk−1,i : sKisoVk,i 4 [︄ − 1 hk−1/2 Tx,k,i−1/2− 1 hk−1/2 Tx,k,i+1/2 ]︄ , Vk+1,iRk+1,i : sKisoVk,i 4 [︄ 1 hk+1/2 Tx,k,i−1/2+ 1 hk+1/2 Tx,k,i+1/2 ]︄ , Vk,i−1Rk,i−1 : sKisoVk,i 4 [︄ 1 di−1/2 Tz,k+1/2,i+ 1 di−1/2 Tz,k−1/2,i ]︄ , Vk,i+1Rk,i+1 : sKisoVk,i 4 [︄ − 1 di+1/2 Tz,k+1/2,i− 1 di+1/2 Tz,k−1/2,i ]︄ .

Here Tz is the shortcut for vertical derivative, and indices point to its location.

Similarly to computations above, the contributions from five cells into Rk,i can

be grouped to reveal that one deals with weighted fluxes through the faces of the cell {k, i}. I do not write them down here.

Term3: (∂zT )2Kiso.

The contribution of the cell {k, i} to the functional F for the term 3 is −s 2_K isoVk,i 4 ⎡ ⎣ (︄ Tk−1,i− Tk,i hk−1/2 )︄2 + (︄ Tk,i− Tk+1,i hk+1/2 )︄2⎤ ⎦.

(32)

In terms of variational calculation it will contribute to the term 3 as

Vk,iRk,i: s2_K isoVk,i 2 [︄ 1 h2 k−1/2 (︂ Tk−1,i− Tk,i )︂ − 1 h2 k+1/2 (︂ Tk,i− Tk+1,i )︂ ]︄ , Vk−1,iRk−1,i: − s2KisoVk,i 2 Tk−1,i− Tk,i h2 k−1/2 , Vk+1,iRk+1,i: s2_K isoVk,i 2 Tk,i− Tk+1,i h2 k+1/2 .

Here the derivative abbreviations are skipped because, as it turns out, this term has to be treated implicitly. For the same reason we have to assemble all the parts of the term 3 which contribute to Vk,iRk,i together:

Vk,iRk,i: s2Kiso h2 k−1/2 (︂ Tk−1,i− Tk,i )︂ [︄ Vk,i 2 + Vk−1,i 2 ]︄ − s 2_K iso h2 k+1/2 (︂ Tk,i− Tk+1,i )︂ [︄ Vk,i 2 + Vk+1,i 2 ]︄ .

Now all the terms are considered. The consideration of the 3D case on trian-gular meshes follows the same procedure. The difference is that each triantrian-gular prism will be split into six sub-prisms and that many triangular prisms con-tribute to a given scalar point.

Stability analysis

The time stepping of isoneutral diffusion will be analyzed in this section. As it is common in such studies, our final intention will be to apply the Fourier analysis. For this reason, I will do some further simplifications. Let us assume the slope scalar s and isoneutral diffusivity Kiso from the matrix K defined in

3.2 to be constant. The mesh will be considered to be uniform with Vi,k = hd,

where h, d are the vertical and horizontal mesh steps respectively. Traditionally, the horizontal diffusivity in ocean codes is considered explicitly. Although the isoneutral diffusion is a replacement for the horizontal diffusion, it needs explicit-implicit treatment, as demonstrated in Lemarié et al. [2012]. Here I explain why. Explicit discretization of equation 2.7 in time can be written either as

T_jn+1− T_jn= Kiso∆t[∂x(∂x+ s∂z)Tn+ ∂z(s∂x+ s2∂zTn)], (3.10)

where n is the time step index and the time step length is denoted as ∆t. In the explicit-implicit form everything concerning the horizontal operator will be explicit, and vertical term ∂zz will be computed implicitly

T_jn+1− Tn j = Kiso∆t [︂ (∂x∂x+ ∂xs∂z+ ∂zs∂x)Tn+ s2∂z∂zTn+1 ]︂ . (3.11)

(33)

For the analysis let us represent the spatial dependence of tracer T as a Fourier harmonic

T = ̃T eikx+imz. (3.12) In the continuous case for such a choice of T , ∂xT = ikT and ∂zT = imT . In

the discrete case considered above, as it can be readily seen, they should be replaced with ∂xT = (2i/d) sin(kd/2)T and ∂zT = (2i/h) sin(mh/2)T .

Further-more, there will be additional factors related to averaging that is present in the expressions for fluxes obtained by variational method. For briefness we will omit them here, but full expressions will be used further in the 3D case. Since in the discrete case |k| ≤ π/d and |m| ≤ π/h, the difference here will not be large to change the answer qualitatively. If we insert this representation in equations (3.10) and (3.11), we will get the relations

T_jn+1− Tn j = Kiso∆t [︂ (−k2− 2mks)Tn j − m 2_s2_Tn j ]︂ (3.13) and T_jn+1− Tn j = Kiso∆t [︂ (−k2− 2mks)Tn j − m 2_s2_Tn+1 j ]︂ . (3.14)

For stability analysis I introduce λ = Tn+1_/Tn_{. In order solutions do not diverge}

with time, i.e., for the scheme stability, it should be |λ| ≤ 1. Expressing λ we will get λe= 1 − Kiso∆t(k2+ 2mks + m2s2) (3.15) and λi = 1 − Kiso∆t(k2+ 2mks) 1 + Kiso∆tm2s2 . (3.16)

Let us start from the explicit case. As it was mentioned before, |λ| has to be less than 1. This brings us to the following relation

Kiso∆t|k2+ 2mks + m2s2| ≤ 2 ⇒ Kiso∆t π2 d2(1 + S) 2 _{≤ 2,} (3.17) where S = sd/h.

For typical Kisoand d it does not cause severe limitations on time step if S is

small. The limitation on the time step ∆t becomes much stronger if S is large. For example, if we take d = 10 km, h = 10 m and the slope s = 10−3, then

S = 1. With Kiso = 300m2/s we can evaluate the limiting time step as ∆t ≤ 8

hours. However, if d is larger, or h is smaller, or s is larger, we can get a strong limitation implying that ∆t gets smaller in S2 _{times. For instance, with S = 10,}

the time step approximately has to be ∆t ≤ 9 min which is already quite small. With S = 20 it gets really limiting. In reality, the problem will happen earlier,

(34)

as diffusion will be combined with advection, and the latter will contribute to the time step limitation as well.

Now let us go back to the explicit-implicit case. We rewrite λi as

λi =

1 − C(1 + 2S)

1 + CS2 = 1 − C

(1 + S)2

1 + CS2, (3.18)

where C = Kiso∆tk2 and S = sm_k. Note that S is different here from the

previous case, yet it has similar sense. The condition |λi| ≤ 1 implies that

C(1 + S)2/(1 + CS2) ≤ 2. Considering the left hand side of this inequality as a function of S, we find that its extremum is achieved at S = −1 or CS = 1.

Putting S = −1 in (3.18) we see that λ = 1 which means that there are no limitations in this case. For the other extremum we will get λ = −C. Therefore we have to require that C = Kiso∆tk2 < 1. We see that in the explicit-implicit

case the limitation does not depend on S (and thus on using ik and im as spectral symbols of derivatives). If we want to avoid oscillations as well, we have to require that λ stays positive:

1 − C(1 + 2S) > 0 ⇒ C < 1

1 + 2S if 1 + 2S > 0.

Negative S will not cause any problems, but positive values of S may limit the admissible time step for a large S. Even this limitation is less strong than the limitation of the explicit case.

In summary, this implies that isoneutral diffusion has to be treated in an explicit-implicit way.

3.2 3D discretization on triangular meshes

A 3D mesh used by FESOM2 is based on a triangular surface mesh and a set of horizontal level surfaces as it was described in chapter 2. Each surface triangle defines a triangular prism going vertically down to the bottom. It is cut by the horizontal level surfaces into smaller prisms. A dual set of prisms is created by first constructing the so-called median-dual cells around each surface vertex (see Fig. 3.2). Such a cell is formed by connecting centers of triangles with mid-points at their edges. For a regular surface mesh made of equilateral triangles, these cells are hexagonal cells of dual mesh. Scalar degrees of freedom are located in the middle of such dual prisms in the vertical, and under the vertices in horizontal directions.

(35)

3.2. 3D DISCRETIZATION ON TRIANGULAR MESHES 33

Vertical discretization is therefore similar to the 2D case above. The horizon-tal discretization needs special treatment. When discussing it, we will be using a horizontal plane and saying that scalar degrees of freedom are at vertices.

The main aspects of the discretization of isoneutral diffusion are based on the variational approach on meshes used by FESOM2 (see details in Danilov et al. [2017]). I implemented this discretization in FESOM2 and complemented it by the analysis of its numerical stability presented further. Some aspects of the performance of the isoneutral diffusion implementation are tested in Scholz et al. [2019].

The vertex placement of scalar degrees of freedom in the horizontal plane implies that the natural positions of horizontal gradients are at the centers of triangular prisms. Vertical gradients are therefore located at vertices of triangu-lar prisms (see Fig.2.3). To have a unique combination of vertical and horizontal derivative, each triangular prism is split into six sub-prisms cut by the mid-plane and the faces of a dual prism (see Fig. 2.4).

As the result, the dissipation functional is obtained as a double sum over triangular prisms and six sub-prisms. The resulting right-hand side of the tracer equation is, however, obtained in the same way as in the 2D case above by differentiation of dissipation functional. The expressions are even bulkier than the expressions for the 2D case above and are not presented here.

The goal of the rest of this section is the time stepping stability analysis for the 3D discretization of FESOM2.

Let us consider a regular triangular surface mesh consisting of equilateral triangles. In FESOM2 scalar degrees of freedom are placed at the vertices of a triangular mesh. As it was mentioned above, the natural location for horizontal gradients is at centroids of triangles. However, a triangular mesh includes tri-angles of two different orientations, as shown in Fig.3.2. We will distinguish be-tween triangles pointing upward (u-triangles) and downward (d-triangles). This is needed because one type of triangles cannot be obtained from the other type by translations, and this has to be taken into account in the Fourier analysis. To carry out such an analysis, we assume that the slope vector s is constant, and that the mesh is uniform in horizontal and vertical direction. We also assume for the horizontal part

Tj = ̃T eikxj+ilyj, (3.19)

where j enumerates all vertices of surface mesh, k, l are the zonal and meridional wavenumbers respectively, and xj, yj are the coordinates of vertex j in the

Reducing spurious diapycnal mixing in ocean circulation models