1 Introduction
The primitive equations are the constitutive system of equations in ocean modeling. This system is composed of a momentum equation, an equation for the ocean thickness, equations for the transport of tracers, as well as an equation of state. The momentum and thickness equations describe the dynamics, i.e. the change in time of the velocity and thickness of the water. The tracer equation describes the transport of tracers, such as temperature, salinity or chemicals. Depending on the number of tracers considered, several equations may be added to and coupled to the dynamics system. The coupling between the tracers and the dynamics depends on the nature of tracers. Temperature and salinity impact the density of the layers of ocean water, influencing the dynamics, and for this reason they are called active tracers. Hence, between the dynamics and active tracers there is a twoway coupling. However, for most tracers the coupling is only oneway, with such tracers being called passive [1].
To this day, the numerical solution of the primitive equations remains a challenging task. One of the main issues is the presence of multiple timescales, where different processes (e.g., external and internal gravity waves, eddies, biochemical reactions) take different times to be completed. Since the primitive equations arise from hyperbolic conservation laws ([2, 3]), explicit time integrators and RungeKutta schemes would hypothetically be good choices for solving it. However, these schemes are not capable of efficiently handling multiple timescales, because the timesteps restrictions are too severe, resulting in a significant degradation of performance. To better handle multiple timescales, methods have been developed with the help of mathematical and algorithmic techniques such as splitting strategies and semiimplicit approaches. Due to the different properties of the dynamics and the tracer equations, schemes have been created to separately deal with the two subsystems ([4, 5]). Methods for the dynamics need to take into account the different wave speeds of the model, with the goal of having a model for which time step sizes are governed by the slow wave speeds or the speed of advection and not by the fast wave speeds as is the case for standard explicit schemes. Examples of timestepping schemes for the dynamics include implicit ([6]) and splitexplicit ([7, 4]) methods. More recently, exponential time differencing (ETD) methods, also known as exponential integrators, have gained attention in the ocean modeling community due to their stability properties that allow time steps considerably larger than those dictated by the CFL condition. In [8], an ETD scheme has been developed for the rotating shallow water equations with multiple horizontal layers, which correspond to a vertical discretization of the primitive equations dynamics in an isopycnal vertical coordinate system. The main idea behind exponential integrators is a splitting of the righthand side term of an equation into a linear part and a remainder, i.e.
with an appropriate choice of the linear operator . For a review of exponential integrators we refer to [9].
In this work, we devise an ETD method for the tracer equation. A tracer is supposed to satisfy a conventional advectiondiffusion equation of the form,
(1) 
where is the tracer studied, is the velocity of water, which is usually split in to the horizontal velocity and the vertical velocity , is a diffusion term, and represents interior sources or sinks. The equation (1) has to be solved on a threedimensional domain, split into a twodimensional horizontal and a vertical coordinate , and . We note that, the tracer equation can usually also be written as
(2) 
with the horizontal diffusion , and the vertical diffusion coefficient . To avoid problems with shocks, the velocity field is usually chosen to be divergence free and tangential to the boundary.
When dealing with the tracer equation, a key point is appropriately including vertical mixing. Tracer vertical mixing usually occurs on small timescales and can be induced by density differences and/or by turbulent motions. For explicit time stepping schemes, the time step requirement is usually set by the horizontal advective CFL condition, hence very small time steps may need to be used with explicit time stepping methods to include realistic vertical mixing of tracers. To avoid this issue, in the timestepping schemes used by popular ocean models, the vertical diffusion term is treated implicitly. In POP [5], the vertical tracer diffusion term is treated with an implicit Euler algorithm, whereas the remaining terms of the equation are treated with a leapfrog algorithm. In MPASOcean [4], tracer equations are stepped forward with the midtime velocity values and this process is repeated in a predictorcorrector way. Implicit vertical mixing of tracers completes each timestep, where, as in POP, the vertical tracer diffusion term is treated with an implicit Euler algorithm.
Another important phenomenon to deal with, besides vertical mixing, is when an inflow of cold water near the coast leads to cold water on top, scenario that creates density and pressure variations and downward motion. In an isopycnal configuration, there is no vertical transport, however, if a mostly Eulerian coordinate system is employed (zstar, zlevel), rapid variations in the pressure/density may induce a fast vertical flow. This phenomenon is amplified by the usage of fine meshes in the vertical, with much smaller vertical than horizontal spacing, for example [4] employs –km horizontal, –m vertical. As we said, many ocean models treat the vertical advection explicitly and so, when a fast transport of water among different layers occurs, the model may be unable to appropriately capture this behavior. Consequently, instabilities in the simulation may occur, causing the need to decrease the timestep. To resolve this issue, we propose an ETD solver where all the vertical terms, i.e. vertical advection and diffusion, are treated with a matrix exponential, whereas the horizontal are dealt with in an explicit way. This means that we are splitting the linear operator into two parts: that accounts for all vertical terms, and that contains all horizontal terms. is then incorporated in the remainder , so the actual linear operator we are working with is . This operator splitting has two advantages. First, by treating exponentially terms related to fast timescales, bigger time steps can be taken, and so computational speedups that can be obtained over other explicit methods. Compared to semiimplicit methods, we expect higher accuracy due to an exact treatment of the fast scales. Second, by including only the operator in the exponential, not the whole operator , the computational cost is reduced, effectively lowering the total computational cost and time of the whole method.
Another important challenge to face when dealing with tracer equations, and in general with primitive equations, is in the numbers of tracers considered. The amount of tracers in an ocean simulation is usually around 40, but may increase up to 70, causing a significant computational load. Hence, efficiently solving multiple tracers equations is an important task in an ocean model. When an ETD scheme is used, efficiency corresponds to a lowcost evaluation of functions. For a given matrix , assembling is generally prohibitive in the large scale context; iterative methods are used instead, i.e. Krylov subspace algorithms ([10, 11]). Scaling and squaring is used usually for dense matrices ([12, 13, 14, 15, 16]) However, it has also been used in the context of multiwavelet Galerkin methods, where a certain level of sparsity can be maintained throughout the iterations of the method and the approximation of ([17]). In this work, together with the usage of a Krylov method, we pursue a second approach based on scaling and squaring relations, focusing on the fact that there is no communication in ocean due to the vertical exponential, since the domaindecomposition is horizontal, and the exponential is vertical. The proposed approximation scheme is based on polynomials of moderate degree that result in a consistent and stable approximation of with low bandwidth. To efficiently solve multiple tracer equations, we intend to preassemble instead of computing for every right hand side . Additionally, since the current ocean models use only 40100 layers, even a full storage of exponential is feasible.
Finally, our solver is compared with existing ocean models, to make sure that it is able to reproduce similar results under the same physical conditions. To do so, the whole primitive equations are solved and two tests from [18] are performed. The tracer equations is coupled with the dynamics, and to solve this latter system a second ETD solver (Exponential Rosenbrock Euler) is used. The results obtained with our ETD solvers are compared with those obtained with three other codes: MPASOcean, MITgcm and MOM.
The paper is organized as follows. Section 2 describes the discretization of the tracer equation both in the vertical and in the horizontal. Section 3 focuses on exponential integrators and their implementation. For the computation of functions, a restarted Krylov subspace method is described, and a scheme based on scaling and squaring relations is proposed. For the latter scheme, an error analysis is provided. In section 4, the proposed ETD solver for the tracer equation is presented, and the properties of the operator splitting are described. The numerical tests are shown and discussed in section 5, while the conclusions follow in section 7.
2 Discretization of the Tracer Equation
In this section, the discretization of the tracer equation (1) in the vertical and in the horizontal is described. In most ocean models, an hydrostatic condition is assumed, leading to describing the primitive equation by the incompressible Boussinesq equations in hydrostatic balance. Following this assumption, the tracer equation in continuous form can be rewritten as
(3) 
where is the tracer, is the horizontal velocity, is the vertical velocity and is the pseudodensity. The variable represents the vertical coordinate and it defined positive upward. and indicate the horizontal and vertical diffusion term, respectively. These terms are defined as
(4)  
(5) 
where and are the horizontal and vertical diffusion, respectively. Without loss of generality, we assume that no forcing term is present, i.e. . This assumption is equivalent to consider that no external factors have an influence on the tracer behavior.
To discretize equation (3), we employ zlevel coordinates in the vertical [19] and a finitevolume method using a Cgrid staggering in the horizontal [20]. For more details about the discretization of all primitive equations, please refer to [4].
2.1 Tracer Equation with vertical discretization
As in MPASOcean [4, 19], the vertical coordinate we use is Arbitrary LagrangianEulerian (ALE). With ALE, several coordinate systems can be specified depending on the application. Common choices for the vertical coordinates include zlevel, where all layers have a fixed thickness except for the top layer, zstar, where all layer thicknesses vary in proportion to the sea surface height, and isopycnal, where there is no vertical transport between layers.
The tracer equation with vertical discretization is
(6) 
where indicates the vertical layer, is the top layer and increases downward up to ; is the mean elevation of the free surface, and the coordinate is positive upward. The variable indicates the transport of fluid from layer to , i.e. across the top interface of layer . The pseudodensity has been replaced by the , which is the layerthickness. The operator , on a generic variable , is the vertical average between the layer and the above layer , i.e.,
(7) 
Finally, and indicate the discretized horizontal and vertical tracer diffusion terms, respectively, and are defined as
(8) 
The discrete operators and , on a generic variable , are defined as
(9)  
(10) 
The vertical tracer diffusion term can actually be rewritten without introducing the operator as .
The choice of the vertical coordinate system is enforced in the computation of the vertical transport. can be found by solving the thickness equation for . The thickness equation discretized in the vertical has the form
(11) 
For a Boussinesq fluid, this equation represents the continuity equation for the pseudodensity . To obtain from (11), all variables at the previous time step must be known, in particular the time derivative of at layer , , must be known. For this purpose, a new quantity, named , is introduced. represents the desired thickness for the new time, and it is used to compute using a firstorder finite difference approximation. In this way, can be found as
(12) 
For isopycnal simulations, the vertical transport is set to zero in (12). For other coordinate systems, like zlevel and zstar, the way is computed determines the type of coordinates chosen. For example, for zlevel vertical coordinates
where is the layer thickness when the ocean is at rest, and is the sea surface height defined as . For zlevel coordinates (and for ztype coordinates in general) the resting thickness is considered constant in each horizontal layer, but for other coordinate systems, like sigma coordinates, varies horizontally in proportion to the column’s total depth. The simulations presented in section 5 use zlevel vertical coordinates. Please refer to [19] for more details about the computation of for other coordinate systems.
2.2 Tracer Equation with horizontal discretization
The horizontal discretization is a Cgrid, finitevolume method applied to a spherical centroidal Voronoi tessellation (SCVT) mesh. Height, tracers, pressure and kinetic energy are defined at centers of the convex polygons, and the velocity is located at cell edges. Vorticity (curl of velocity) is defined at cell vertices. In the following, the subscripts and indicate the discretized variables through cell centers and edges, respectively. Since we are focusing on the discretization of the tracer equation only, we will not work with variables and operators defined at cell vertices.
The tracer equation with horizontal discretization is
(13)  
(14) 
Each variable now has two subscripted indices, the first indicating the vertical layer, and the second indicating its position on the horizontal grid, namely either or . Colons in subscripts may be places as second index to indicate that multiple edges or cell centers are used in computing the horizontal operator. For a generic variable , the symbol represents the averaging of the variable from two adjacent centers to the corresponding edge. We would like to point out that the vertical transport through the sea surface and at the bottom surface is zero, i.e. and . Moreover, we consider on all boundary edges.
For a generic vector field
and variable , the discrete horizontal operators and are defined as(15)  
(16) 
indicates the Voronoi cell area, is the distance between cell centers, is edge length and represents the sign of the vector at edge with respect to cell . The sets are the edges about cell , and the sets are the cells neighboring edge . Thus, the divergence moves from edges to cellcentered quantity, while the gradient moves from cell centers to edges.
3 Exponential Time Integration
This section describes exponential time differencing methods, that later will be used to solve the tracer equation. ETD methods have already been employed to solve the single layer ([21, 22, 23]) and multilayer ([8]) shallow water equations.
3.1 Exponential Integrators
Let
be a system of partial differential equations (PDEs), where
denotes the vector of the solution variables for , and is the righthandside term. The interval refers to one time step. The main idea behind exponential integrators is a splitting of the righthandside term into a linear part and a remainder, i.e.(17) 
where represents a linear operator, and denotes the remainder, which in general is nonlinear. Applying the variation of constants formula to equation (17), the solution at time , i.e. , is obtained as
(18) 
At this point, to build a concrete exponential integrator, an approximation of must be considered. By substituting with its Taylor expansion truncated at (), namely
the solution can be approximated by
The above expression can be rewritten as
(19) 
by introducing the, so called, functions defined by
(20)  
(21) 
By performing the change of variable , (21) can be obtained from (20). For , we have that .
The parameter in (19) indicates the number of stages of the method. By taking , the exponential Euler method is obtained as
(22) 
where , with
indicating the identity matrix. For a generic linear operator
, this method is firstorder accurate, but if is the Jacobian matrix of the system evaluated at , namely , then the method becomes secondorder accurate [9, 23]. In this case, the scheme (22) is called Exponential Rosenbrock Euler. By taking , a secondorder singlestep method with two stages is obtained as(23)  
(24) 
where the first derivative of is approximated with a firstorder finitedifference approximation
ETD2RK fulfills the stiff order conditions of order two described in [24], section 5.1. However, by relaxing some of these conditions, different schemes can be obtained that are more convenient for computations. For instance, the following twostage predictorcorrector scheme makes use only of the function.
(25)  
(26) 
This scheme fulfills the nonstiff order conditions up to order two, and the stiff order conditions up to order one. From an implementational point of view, it is simpler than ETD2RK, since does not need to be assembled, and a further advantage in terms of computational time can be obtained if the matrix can be precomputed and stored efficiently.
3.2 Computation of the functions
The computation of the functions is of major relevance for ETD methods. Let us assume that is in , where and is the index of the function. Fist of all, let us consider the special case , for which we have . The most popular method for computing the matrix exponential is the scaling and squaring algorithm [25, 15]. In [14], the scaling and squaring algorithm used by the MATLAB function expm is described, and a variation of this method that alleviates overscaling problems is presented in [12]. The Expokit package [26] uses a scaling and squaring method as well for computing matrix exponentials. An alternative to scaling and squaring methods is given by Krylov subspace projections [27, 11]. Krylov schemes can be used for evaluating any function with index . With Krylov schemes, matrixvector products of the form are computed, without the actual construction of the matrix . To overcome some memory issues associated with standard Krylov methods, restarted schemes have been developed [10]. The C/C++ library SLEPc [28] uses the restarted algorithm presented in [10].
In the following, a brief description of the algorithm presented in [10] is given, since this method is used in the numerical results section, and scaling and squaring algorithm for indexes
is presented, for which an error estimate is derived.
3.2.1 Krylov subspace approximation
The Krylov subspace approximation has been found very efficient to compute matrixvector products due to the optimality of the matrix polynomials produced by Krylov methods. The idea behind a Krylov subspace approach is to approximate the vector , which lives in , in a smaller space of dimension . The Krylov approximation of is based on an Arnoldi decomposition of
(27) 
where is a matrix whose columns form an orthogonal basis for the Krylov subspace of dimension
is an upper Hessenberg matrix, and denotes the th vector of the standard basis of . The matrices and are such that
(28) 
therefore can be seen as the projection of the action of to the Krylov subspace .
Standard Krylov methods requires the storage of , i.e. the storage of vectors of size , which may be costly for moderate to large values of . Therefore, a standard Krylov subspace method may be impractical because of the storage requirements associated with . To overcome this issue, the Arnoldi approximation could be modified in a way that allows the construction of successively better approximations of based on a sequence of Krylov spaces of small dimension. Methods based on this approach are called restarted Krylov subspace methods, and they use an Arnoldilike decomposition where a sequence of ascending (not necessarily orthonormal) basis vectors is introduced. We now focus on the restarted Krylov subspace algorithm presented in [10], which can be summarized as follows. An Arnoldilike decomposition of is constructed
(29) 
where ,
, and
, .
Since (29) is an Arnoldilike decomposition, the columns of are only blockwise orthonormal.
The matrices ,
and the scalars
are obtained from proper Arnoldi decompositions.
Setting
the approximation of after restart cycles is given by
(30) 
Therefore, the approximation is obtained from the previous approximation plus a correction term. Only has to be stored from the previous cycle of the algorithm, and the matrix (together with ) can be discarded after computing . An efficient implementation of this algorithm can be found in [29], where it is also shown how to stably compute the coefficient vector .
3.2.2 Scaling and Squaring
We now present a scaling and squaring method for the computation of . We are going to base the scaling and squaring method on the recursive relations for the functions; cf. [30, 31, 32]. They are given as (see [32, Lemma 3]):
(31)  
(32) 
Note that in both expressions the matrix function of the matrix is expressed as a product of two matrix functions of plus some correction terms. The first identity is simpler, but the second one reduces the number of correction terms, which is slightly more efficient for computations; however [32] mentions that the first form is more stable in numerical experiments.
The first formula follows from the definition (20) using , and a splitting of the integral into the intervals and
where the second integral was shifted to . Now, using the binomial identity in the second term and the definition of and shows (31). Now, to derive the second identity, we use in a first step
Proceeding iteratively times with the first term, we obtain
using this identity for and respectively, we obtain (32).
In order to evaluate the matrix functions, we are going exploit these recursive relations to reduce the computation to matrix functions of a scaled matrix . Then, we use a polynomial approximation for this matrix. In order to ensure accuracy and stability, we do this in the following way: First, let
be a polynomial of degree that approximates the exponential function up to order . Here, is the Taylor approximation to , and is a remainder. Then, for all we define the consistent polynomial approximations to the functions as
(33) 
Now, we define the higher order recursive approximations for to using (31) as
(34) 
This definition ensures that the resulting functions have similar properties as the original functions.
Proposition 1.
For any and there holds
(35)  
(36) 
Proof.
For this holds according to definition. For higher we follow an induction argument.
First, we show (35). We use the recursive definition (34) to obtain
and use the induction hypothesis to obtain
using again the recursive definition (34) for and the wellknown summation formula of binomial coefficients. Dividing by yields (35). Concerning (36), we note that is suffices to repeatedly apply (35) for , , …, . ∎
Remark 1.
Remark 2.
In (45)(46), which is the ETD method we actually use to later solve the tracer equation, the only function needed is . For this specific case, the recursion formula (31) looks like
and
Using the above relations, the following algorithm computes from .

Define

For , Given and , compute as

Given , compute as
Finally, we provide an error estimate for the approximation . First, we consider the polynomial approximation on a subset of the complex plane. To prepare for the general case, we let be some compact subset of the negative complex plane shifted by , and assume that the underlying polynomial fulfills the stability assumption
(37) 
Moreover, since is compact, there exists such that
(38) 
This will be the basis of the error estimates.
Remark 3.
There are many situations where this assumption is fulfilled. In the case that is the Taylor polynomial and , can be chosen as the intersection of the negative halfplane and the wellknown stability region of a RungeKutta scheme of order . In this case, due to the relation and we have . For an overview over known results on polynomials with optimal stability properties for various forms of and a computational approach to determine them, we refer to Ketcheson [33].
We first investigate the case of the matrix exponential .
Proposition 2.
Let the polynomial fulfill (37). Then for all we have the stability and approximation properties:
(39)  
(40) 
which hold for all , for the enlarged stability radii .
Comments
There are no comments yet.