Mixed finite elements for global tide models
Abstract
We study mixed finite element methods for the linearized rotating shallow water equations with linear drag and forcing terms. By means of a strong energy estimate for an equivalent secondorder formulation for the linearized momentum, we prove longtime stability of the system without energy accumulation – the geotryptic state. A priori error estimates for the linearized momentum and free surface elevation are given in as well as for the time derivative and divergence of the linearized momentum. Numerical results confirm the theoretical results regarding both energy damping and convergence rates.
Key words. Finite element method, global tidal models, H(div) elements, energy estimates.
AMS subject classifications. 65M12, 65M60, 35Q86
1 Introduction
Finite element methods are attractive for modelling the world’s oceans since implemention with triangular cells provides a means to accurately represent coastlines and topography [34]. In the last decade or so, there has been much discussion about the best choice of mixed finite element pairs to use as the horizontal discretization for atmosphere and ocean models. In particular, much attention has been paid to the properties of numerical dispersion relations obtained when discretizing the rotating shallow water equations [10, 6, 5, 21, 31, 23, 30, 22]. In this paper we take a different angle, and study the behavior of discretizations of forceddissipative rotating shallowwater equations, which are used for predicting global barotropic tides. The main point of interest here is whether the discrete solutions approach the correct longtime solution in response to quasiperiodic forcing. In particular, we study the behavior of the linearized energy. Since this energy only controls the divergent part of the solution, as we shall see later, it is important to choose finite element spaces where there is a natural discrete Helmholtz decomposition, and where the Coriolis term projects the divergent and divergencefree components of vector fields correctly onto each other. Hence, we choose to concentrate on the mimetic, or compatible, finite element spaces (i.e. those which arise naturally from the finite element exterior calculus [1]) which were proposed for numerical weather prediction in [7]. In that paper, it was shown that the discrete equations have an exactly steady geostrophic state (a solution in which the Coriolis term balances the pressure gradient) corresponding to each of the divergencefree velocity fields in the finite element space; this approach was extended to develop finite element methods for the nonlinear rotating shallowwater equations on the sphere that can conserve energy, enstrophy and potential vorticity [29, 8, 26]. Here, we shall make use of the discrete Helmholtz decomposition in order to show that mixed finite element discretizations of the forceddissipative linear rotating shallowwater equations have the correct longtime energy behavior. Since we are studying linear equations, these energy estimates then provide finite time error bounds.
Predicting past and present ocean tides is important because they have a strong impact on sediment transport and coastal flooding, and hence are of interest to geologists. Recently, tides have also received a lot of attention from global oceanographers since breaking internal tides provide a mechanism for vertical mixing of temperature and salinity that might sustain the global ocean circulation [27, 12]. A useful tool for predicting tides are the rotating shallow water equations, which provide a model of the barotropic (i.e., depthaveraged) dynamics of the ocean. When modelling global barotropic tides away from coastlines, the nonlinear advection terms are very weak compared to the drag force, and a standard approach is to solve the linear rotating shallowwater equations with a parameterised drag term to model the effects of bottom friction, as described in [20]. This approach can be used on a global scale to set the boundary conditions for a more complex regional scale model, as was done in [14], for example. Various additional dissipative terms have been proposed to account for other dissipative mechanisms in the barotropic tide, due to baroclinic tides, for example [16].
As mentioned above, finite element methods provide useful discretizations for tidal models since they can be used on unstructured grids which can seamless couple global tide structure with local coastal dynamics. A discontinuous Galerkin approach was developed in [32], whilst continuous finite element approaches have been used in many studies [18, 24, 11, for example]. The lowest order RaviartThomas element for velocity combined with for height was proposed for coastal tidal modeling in [33]; this pair fits into the framework that we discuss in this paper.
In this paper we will restrict attention to the linear bottom drag model as originally proposed in [20]. We are aware that the quadratic law is more realistic, but the linear law is more amenable to analysis and we believe that the correct energy behavior of numerical methods in this linear setting already rules out many methods which are unable to correctly represent the longtime solution which is in geotryptic balance (the extension to geostrophic balance of the three way balance between Coriolis, the pressure gradient and the dissipative term). In the presence of quasiperiodic timevarying tidal forcing, the equations have a timevarying attracting solution that all solutions converge to as . In view of this, we prove the following results which are useful to tidal modellers (at least, for the linear law):

For the mixed finite element methods that we consider, the spatial semidiscretization also has an attracting solution in the presence of timevarying forcing.

This attracting solution converges to timevarying attracting solution of the unapproximated equations.
Global problems require tidal simulation on manifolds rather than planar domains. For simplicity, our description and analysis will follow the latter case. However, our numerical results include the former case. Recently, Holst and Stern [15] have demonstrated that finite element analysis on discretized manifolds can be handled as a variational crime. We summarize these findings and include an appendix at the end demonstrating how to apply their techniques to our own case. This suggests that the extension to manifolds presents technicalities rather than difficulties to the analysis we provide here.
The rest of this paper is organised as follows. In Section LABEL:se:model we describe the finite element modelling framework which we will analyse. In Section LABEL:se:prelim we provide some mathematical preliminaries. In Section LABEL:se:energy we derive energy stability estimates for the finite element tidal equations. In Section LABEL:se:error we use these energy estimates to obtain error bounds for our numerical solution. Appendix LABEL:ap:bendy includes the discussion of embedded manifolds.
2 Description of finite element tidal model
We start with the nondimensional linearized rotating shallow water model with linear drag and forcing on a (possibly curved) two dimensional surface , given by
\hb@xt@.01(2.1) 
where is the nondimensional two dimensional velocity field tangent to , is the nondimensional free surface elevation above the height at state of rest, is the (spatially varying) tidal forcing, is the Rossby number (which is small for global tides), is the spatiallydependent nondimensional Coriolis parameter which is equal to the sine of the latitude (or which can be approximated by a linear or constant profile for local area models), is the Burger number (which is also small), is the (spatially varying) nondimensional drag coefficient and is the (spatially varying) nondimensional fluid depth at rest, and and are the intrinsic gradient and divergence operators on the surface , respectively.
We will work with a slightly generalized version of the forcing term, which will be necessary for our later error analysis. Instead of assuming forcing of the form , we assume some , giving our model as
\hb@xt@.01(2.2) 
It also becomes useful to work in terms of the linearized momentum rather than velocity. After making this substitution and dropping the tildes, we obtain
\hb@xt@.01(2.3) 
A natural weak formulation of this equations is to seek and so that
\hb@xt@.01(2.4) 
We now develop mixed discretizations with and . Conditions on the spaces are the commuting projection and divergence mapping onto . We define and as solutions of the discrete variational problem
\hb@xt@.01(2.5) 
We will eventually obtain stronger estimates by working with an equivalent secondorder form. If we take the time derivative of the first equation in (LABEL:eq:discrete_mixed) and use the fact that , we have
\hb@xt@.01(2.6) 
where . This is a restriction of
\hb@xt@.01(2.7) 
which is the variational form of
\hb@xt@.01(2.8) 
to the mixed finite element spaces.
We have already discussed mixed finite elements’ application to tidal models in the geophysical literature, but this work also builds on existing literature for mixed discretization of the acoustic equations. The first such investigation is due to Geveci [13], where exact energy conservation and optimal error estimates are given for the semidiscrete firstorder form of the model wave equation. Later analysis [9, 17] considers a second order in time wave equation with an auxillary flux at each time step. In [19], Kirby and Kieu return to the firstorder formulation, giving additional estimates beyond [13] and also analyzing the symplectic Euler method for time discretization. From the standpoint of this literature, our model (LABEL:eq:thepde) appends additional terms for the Coriolis force and damping to the simple acoustic model. We restrict ourselves to semidiscrete analysis in this work, but pay careful attention the extra terms in our estimates, showing how study of an equivalent secondorder equation in proves proper longterm behavior of the model.
3 Mathematical preliminaries
For the velocity space , we will work with standard mixed finite element spaces on triangular elements, such as RaviartThomas (RT), BrezziDouglasMarini (BDM), and BrezziDouglasFortinMarini (BDFM) [28, 4, 3]. We label the lowestorder RaviartThomas space with index , following the ordering used in the finite element exterior calculus [1]. Similarly, the lowestorder BrezziDouglasFortinMarini and BrezziDouglasMarini spaces correspond to as well. We will always take to consist of piecewise polynomials of degree , not constrained to be continuous between cells. In the case of domains with boundaries, we require the strong boundary condition on all boundaries.
In the main part of this paper we shall present results assuming that the domain is a subset of , i.e. flat geometry. In the Appendix, we describe how to extend these results to the case of embedded surfaces in .
Throughout, we shall let denote the standard norm. We will frequently work with weighted norms as well. For a positivevalued weight function , we define the weighted norm
\hb@xt@.01(3.1) 
If there exist positive constants and such that almost everywhere, then the weighted norm is equivalent to the standard norm by
\hb@xt@.01(3.2) 
A CauchySchwarz inequality
\hb@xt@.01(3.3) 
holds for the weighted inner product, and we can also incorporate weights into CauchySchwarz for the standard inner product by
\hb@xt@.01(3.4) 
We refer the reader to references such as [3] for full details about the particular definitions and properties of these spaces, but here recall several facts essential for our analysis. For all velocity spaces we consider, the divergence maps onto . Also, the spaces of interest all have a projection, that commutes with the projection into :
\hb@xt@.01(3.5) 
for all and any . We have the error estimate
\hb@xt@.01(3.6) 
when . Here, for the BDM spaces but for the RT or BDFM spaces. The projection also has an error estimate for the divergence
\hb@xt@.01(3.7) 
for all the spaces of interest, whilst the pressure projection has the error estimate
\hb@xt@.01(3.8) 
Here, and are positive constants independent of , , and , although not necessarily of the shapes of the elements in the mesh.
We will utilize a Helmholtz decomposition of under a weighted inner product. For a very general treatment of such decompositions, we refer the reader to [2]. For each , there exist unique vectors and such that , , and also . That is, is decomposed into the direct sum of solenoidal vectors, which we denote by
\hb@xt@.01(3.9) 
and its orthogonal complement under the inner product, which we denote by
\hb@xt@.01(3.10) 
Functions in satisfy a generalized PoincaréFriedrichs inequality, that there exists some such that
\hb@xt@.01(3.11) 
We may also use norm equivalence to write this as
\hb@xt@.01(3.12) 
Because our mixed spaces are contained in , the same decompositions can be applied, and the PoincaréFriedrichs inequality holds with a constant no larger than .
4 Energy estimates
In this section, we develop in stability estimates for our system, obtained by energy techniques. Supposing that there is no forcing or damping (), we pick and in (LABEL:eq:discrete_mixed), and find that
\hb@xt@.01(4.1) 
Since pointwise, we add these two equations together to find
\hb@xt@.01(4.2) 
Hence, we have the following.
Proposition 4.1
In the absence of damping or forcing, the quantity
\hb@xt@.01(4.3) 
is conserved exactly for all time.
Now suppose that still but that pointwise in . The same considerations now lead to
\hb@xt@.01(4.4) 
so that
Proposition 4.2
In the absence of forcing, but with , the quantity defined in (LABEL:eq:firstorderenergy) satisfies
In the presence of forcing and dissipation, it is also possible to make estimates showing worstcase linear accumulation of the energy over time.
Proposition 4.3
With nonzero , we have that for all time ,
\hb@xt@.01(4.5) 
Proof. We choose and as without forcing, and find that
CauchySchwarz, Young’s inequality, and norm equivalence give
The result follows by dropping the positive term from the lefthand side and integrating.
However, linear energy accumulation is not observed for actual tidal motion, so we expect a stronger result to hold. Turning to the second order equation (LABEL:eq:secondorderdiscrete), we begin with vanishing forcing and damping terms, putting to find
\hb@xt@.01(4.6) 
which simplifies to
\hb@xt@.01(4.7) 
so that the quantity
\hb@xt@.01(4.8) 
is conserved exactly for all time.
If is nonzero, we have that
\hb@xt@.01(4.9) 
which implies that is nonincreasing, although with no particular decay rate.
Now, we develop more refined technique based on the Helmholtz decomposition that gives a much stronger damping result. We can write in the weighted decomposition. We let be a scalar to be determined later and let the test function in (LABEL:eq:secondorderdiscrete) be . This gives
\hb@xt@.01(4.10) 
and we rewrite the lefthand side so that
\hb@xt@.01(4.11) 
We use the fact that
and also that is orthogonal to to rewrite the lefthand side as
\hb@xt@.01(4.12) 
This has the form of an ordinary differential equation
\hb@xt@.01(4.13) 
where
\hb@xt@.01(4.14) 
and
\hb@xt@.01(4.15) 
By showing that for suitably chosen , both and are comparable to defined in (LABEL:eq:Edef), we can obtain exponential damping of the energy.
Lemma 4.4
Suppose that
\hb@xt@.01(4.16) 
Then
\hb@xt@.01(4.17) 
Proof. We bound the term , with CauchySchwarz, PoincareFriedrichs (LABEL:eq:pf), and weighted Young’s inequality with :
\hb@xt@.01(4.18) 
So, then, we have
\hb@xt@.01(4.19) 
and the result follows thanks to the assumption (LABEL:eq:alpha1).
Showing that is bounded above by a constant times is straightforward, but not needed for our damping results.
Lemma 4.5
Suppose that
\hb@xt@.01(4.20) 
where
\hb@xt@.01(4.21) 
Then
\hb@xt@.01(4.22) 
Proof. We use Cauchy Schwarz, the bounds and , and Young’s inequality with weight to write
\hb@xt@.01(4.23) 
Next, it remains to select and to make the coefficients of each norm positive and also balance the terms. First, we pick
and calculating that
we have that
\hb@xt@.01(4.24) 
We let be the solution to
so that
\hb@xt@.01(4.25) 
If we pick , then we have the lower bound for is exactly . However, we are also constrained to pick in order to guarantee that the lower bounds for is positive as well. If we have , then
and so we also have
\hb@xt@.01(4.26) 
We combine these two propositions to give our exponential damping result.
Theorem 4.6
Let and be defined by (LABEL:eq:alpha1) and (LABEL:eq:alpha2), respectively. Then, for any , and any , we have
\hb@xt@.01(4.27) 
Proof. In light of (LABEL:eq:ode), (LABEL:eq:Bbelow), and the lower bound in (LABEL:eq:AequivE), we have that
\hb@xt@.01(4.28) 
so that
\hb@xt@.01(4.29) 
Using the upper and lower bounds of in (LABEL:eq:AequivE) gives the desired estimate.
This result shows that the damping term drives an unforced system to one with a steady, solenoidal velocity field, in which the Coriolis force balances the pressure gradient term, i.e. in a state of geostrophic balance. Using the second equation in (LABEL:eq:discrete_mixed), we also know that the linearized height disturbance is steady in time in this case. These facts together lead to an elliptic equation for the steady state
\hb@xt@.01(4.30) 
It is easy to see that this problem is coercive on the divergencefree subspaces and thus is wellposed. Hence, with zero forcing, both and equal zero is the only solution. The zeroenergy steady state then cannot have a nonzero solenoidal part. Moreover, the exponentially decay of toward zero forces to reach its steady state quickly, driving both and toward zero at an exponential rate. Finally, since almost everywhere, the exponential damping of also forces toward its zero steady state at the same rate.
Now, we turn to the case where the forcing term is nonzero, adapting this damping result to give longtime stability. The same techniques as before now lead to
\hb@xt@.01(4.31) 
Theorem 4.7
For any and
\hb@xt@.01(4.32) 
we have the bound
\hb@xt@.01(4.33) 
Proof. We bound the righthand side of (LABEL:eq:ode2) by
\hb@xt@.01(4.34) 
We put to find
\hb@xt@.01(4.35) 
This turns (LABEL:eq:ode2) into the differential inequality
\hb@xt@.01(4.36) 
Using (LABEL:eq:Bbelow), we obtain
\hb@xt@.01(4.37) 
At this point, we specify so that, with (LABEL:eq:AequivE) we have
\hb@xt@.01(4.38) 
This leads to the bound on
\hb@xt@.01(4.39) 
Using (LABEL:eq:AequivE) again gives the desired result.
These stability results have important implications for tidal computations. Theorem LABEL:thm:dampedstable shows longtime stability of the system. Our stability result also shows that the semidiscrete method captures the threeway geotryptic balance between Coriolis, pressure gradients, and forcing. Moreover, we also can demonstrate that “spinup”, the process by which in practice tide models are started from an arbitrary initial condition and run until they approach their longterm behavior, is justified for this method. To see this, the difference between any two solutions with equal forcing but differing initial conditions will satisfy the same (LABEL:eq:secondorderdiscrete) with nonzero initial conditions and zero forcing. Consequently, the difference must approach zero exponentially fast. This means that we can define a global attracting solution in the standard way (that is, take , for and as the solution starting from zero initial conditions at and define the global attracting solution as the limit as ), to which the solution for any condition becomes exponentially close in finite time. The error estimates we demonstrate in the next section then can be used to show that the semidiscrete finite element solution for given initial conditions approximates this global attracting solution arbitrarily well by picking large enough that the difference between the exact solution with those initial conditions and the global attracting solution is small and then letting be small enough that the finite element solution approximates that exact solution well.
5 Error estimates
Optimal a priori error estimates follow by applying our stability estimates to a discrete equation for the difference between the numerical solution and a projection of the true solution. We define
\hb@xt@.01(5.1) 
The projections and satisfy the firstorder system
\hb@xt@.01(5.2) 
Subtracting the discrete equation (LABEL:eq:discrete_mixed) from this gives
\hb@xt@.01(5.3) 
By choosing the initial conditions for the discrete problem as and , the initial conditions for these error equations are
\hb@xt@.01(5.4) 
We start with estimates for the height and momentum variables, based on the stability result for the first order system.
Proposition 5.1
For any , provided that ,
\hb@xt@.01(5.5) 
Proof. We apply Proposition LABEL:prop:firstorderstability to (LABEL:eq:firstordererror) to find
\hb@xt@.01(5.6) 
Note that for any ,
Using this, that , and norm equivalence bounds the righthand side above by
and the approximation estimate (LABEL:eq:PiL2) finishes the proof.
Since
we combine this result with the approximation estimates to obtain
Theorem 5.2
If the above hypotheses hold, and also and , we have the error estimate
\hb@xt@.01(5.7) 
Note that our bound on the error equations in Proposition LABEL:prop:L2discerr depend only on the approximation properties of the velocity space, while the full error in the finite element solution depends on the approximation properties of both spaces. Consequently, the velocity approximation using BDM elements is suboptimal. Using RT or BDFM elements, both fields are approximated to optimal order.
Now, we use our estimates based on the secondorder system to obtain error estimates for the time derivative and divergence of the momentum. The projection satisfies the perturbed equation
\hb@xt@.01(5.8) 
As in the firstorder case, we have , and subtracting (LABEL:eq:secondorderdiscrete) from (LABEL:eq:secondorderprojeq) gives