Mixed finite elements for global tide models

# Mixed finite elements for global tide models

Colin Cotter Department of Mathematics, Imperial College London, South Kensington Campus, London SW7 2AZ. This work was supported by NERC grant NE/I016007/1    Robert C. Kirby Department of Mathematics, Baylor University, One Bear Place #97328, Waco, TX 76798-73278. (robert_kirby@baylor.edu). This work was supported by NSF grant CCF-1117794.
###### 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 second-order formulation for the linearized momentum, we prove long-time 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 forced-dissipative rotating shallow-water equations, which are used for predicting global barotropic tides. The main point of interest here is whether the discrete solutions approach the correct long-time solution in response to quasi-periodic 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 divergence-free 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 divergence-free velocity fields in the finite element space; this approach was extended to develop finite element methods for the nonlinear rotating shallow-water 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 forced-dissipative linear rotating shallow-water equations have the correct long-time 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., depth-averaged) 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 shallow-water 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 Raviart-Thomas 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 long-time 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 time-varying tidal forcing, the equations have a time-varying 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):

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

2. This attracting solution converges to time-varying 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

 ut+fϵu⊥+βϵ2∇(η−η′)+Cu=0,ηt+∇⋅(Hu)=0, \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 spatially-dependent non-dimensional 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

 ut+fϵu⊥+βϵ2∇η+Cu=F,ηt+∇⋅(Hu)=0. \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

 1Hut+fHϵu⊥+βϵ2∇η+CHu=F,ηt+∇⋅u=0. \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 second-order form. If we take the time derivative of the first equation in (LABEL:eq:discrete_mixed) and use the fact that , we have

 (1Huh,tt,vh)+1ϵ(fHu⊥h,t,vh)+βϵ2(∇⋅uh,∇⋅vh)+(CHuh,t,vh)=(˜F,vh), \hb@xt@.01(2.6)

where . This is a restriction of

 (1Hutt,v)+1ϵ(fHu⊥t,v)+βϵ2(∇⋅u,∇⋅v)+(CHut,v)=(˜F,vh), \hb@xt@.01(2.7)

which is the variational form of

 1Hutt+fHu⊥t−βϵ2∇(∇⋅u)+CHut=˜F, \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 first-order 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 first-order 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 second-order equation in proves proper long-term 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 Raviart-Thomas (RT), Brezzi-Douglas-Marini (BDM), and Brezzi-Douglas-Fortin-Marini (BDFM) [28, 4, 3]. We label the lowest-order Raviart-Thomas space with index , following the ordering used in the finite element exterior calculus [1]. Similarly, the lowest-order Brezzi-Douglas-Fortin-Marini and Brezzi-Douglas-Marini 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 positive-valued weight function , we define the weighted norm

 ∥g∥2κ=∫Ωκ|g|2dx. \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

 √κ∗∥g∥≤∥g∥κ≤√κ∗∥g∥. \hb@xt@.01(3.2)

A Cauchy-Schwarz inequality

 (κg1,g2)≤∥g1∥κ∥g2∥κ \hb@xt@.01(3.3)

holds for the weighted inner product, and we can also incorporate weights into Cauchy-Schwarz for the standard inner product by

 (g1,g2)=(√κg1,1√κg2)≤∥g1∥κ∥g2∥1κ. \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 :

 (∇⋅Πu,wh)=(π∇⋅u,wh) \hb@xt@.01(3.5)

for all and any . We have the error estimate

 ∥u−Πu∥≤CΠhk+σ|u|k \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

 ∥∇⋅(u−Πu)∥≤CΠhk|∇⋅u|k \hb@xt@.01(3.7)

for all the spaces of interest, whilst the pressure projection has the error estimate

 ∥η−πη∥≤Cπhk|η|k. \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

 N(∇⋅)={u∈V:∇⋅u=0}, \hb@xt@.01(3.9)

and its orthogonal complement under the inner product, which we denote by

 N(∇⋅)⊥={u∈V:(1Hu,v)=0, ∀v∈N(∇⋅)}. \hb@xt@.01(3.10)

Functions in satisfy a generalized Poincaré-Friedrichs inequality, that there exists some such that

 ∥∥uD∥∥1H≤CP∥∥∇⋅uD∥∥1H. \hb@xt@.01(3.11)

We may also use norm equivalence to write this as

 ∥∥uD∥∥1H≤CP√H∗∥∥∇⋅uD∥∥. \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

 (1Huh,t,uh)+1ϵ(fHu⊥h,uh)−βϵ2(ηh,∇⋅uh)=0,βϵ2(ηh,t,ηh)+βϵ2(∇⋅uh,ηh)=0. \hb@xt@.01(4.1)

Since pointwise, we add these two equations together to find

 12ddt∥uh∥21H+β2ϵ2ddt∥ηh∥2=0. \hb@xt@.01(4.2)

Hence, we have the following.

###### Proposition 4.1

In the absence of damping or forcing, the quantity

 E1(t)=12∥uh∥21H+β2ϵ2∥ηh∥2 \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

 12ddt∥uh∥21H+β2ϵ2ddt∥ηh∥2+∥uh∥2CH=0, \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

 ddtE1(t)≤0.

In the presence of forcing and dissipation, it is also possible to make estimates showing worst-case linear accumulation of the energy over time.

###### Proposition 4.3

With nonzero , we have that for all time ,

 E1(t)≤E1(0)+12C∗∫t0∥F(⋅,s)∥2Hds \hb@xt@.01(4.5)

Proof. We choose and as without forcing, and find that

 ddtE1(t)+∥u(⋅,t)∥2CH=(F,uh).

Cauchy-Schwarz, Young’s inequality, and norm equivalence give

 ddtE1(t)+C∗2∥uh(⋅,t)∥21H≤12C∗∥F(⋅,t)∥2H

The result follows by dropping the positive term from the left-hand 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

 (1Huh,tt,uh,t)+1ϵ(fHu⊥h,t,uh,t)+βϵ2(∇⋅uh,∇⋅uh,t)=0, \hb@xt@.01(4.6)

which simplifies to

 12ddt∥uh,t∥21H+β2ϵ2ddt∥∇⋅uh∥2=0, \hb@xt@.01(4.7)

so that the quantity

 E(t)=12∥uh,t∥21H+β2ϵ2∥∇⋅uh∥2 \hb@xt@.01(4.8)

is conserved exactly for all time.

If is nonzero, we have that

 12ddt∥uh,t∥21H+β2ϵ2ddt∥∇⋅uh∥2+∥uh,t∥2CH=0, \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

 (1Huh,tt,uh,t+αuDh)+1ϵ(fHu⊥h,t,uh,t+αuDh)+βϵ2(∇⋅uh,∇⋅(uh,t+αuDh))+(CHuh,t,uh,t+αuDh)=0, \hb@xt@.01(4.10)

and we rewrite the left-hand side so that

 12ddt∥uh,t∥21H+α(1Huh,tt,uDh)+αϵ(fHu⊥h,t,uDh)+β2ϵ2ddt∥∥∇⋅uDh∥∥2+αβϵ2∥∥∇⋅uDh∥∥2+∥uh,t∥2CH+α(CHuh,t,uDh)=0. \hb@xt@.01(4.11)

We use the fact that

 ddt(1Huh,t,uDh)=(1Huh,tt,uDh)+(1Huh,t,udh,t)

and also that is -orthogonal to to rewrite the left-hand side as

 ddt[12∥uh,t∥21H+α(1Huh,t,uDh)+β2ϵ2∥∥∇⋅uDh∥∥2]+αϵ(fHu⊥h,t,uDh)+αβϵ2∥∥∇⋅uDh∥∥2+∥uh,t∥2CH−α∥∥uDh,t∥∥21H+α(CHuh,t,uDh)=0. \hb@xt@.01(4.12)

This has the form of an ordinary differential equation

 A′(t)+B(t)=0, \hb@xt@.01(4.13)

where

 A(t)=12∥uh,t∥21H+α(1Huh,t,uDh)+β2ϵ2∥∥∇⋅uDh∥∥2 \hb@xt@.01(4.14)

and

 B(t)=αϵ(fHu⊥h,t,uDh)+αβϵ2∥∥∇⋅uDh∥∥2+∥uh,t∥2CH−α∥∥uDh,t∥∥21H+α(CHuh,t,uDh). \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

 α≤α1≡√βH∗2Cpϵ. \hb@xt@.01(4.16)

Then

 12E(t)≤A(t)≤32E(t). \hb@xt@.01(4.17)

Proof. We bound the term , with Cauchy-Schwarz, Poincare-Friedrichs (LABEL:eq:pf), and weighted Young’s inequality with :

 (1HuDh,t,uDh)≤CP2√H∗[ϵ√β∥uh,t∥21H+√βϵ∥∥∇⋅uDh∥∥2]=CPϵ√H∗β[12∥uh,t∥21H+β2ϵ2∥∥∇⋅uDh∥∥2]=CPϵ√H∗βE(t). \hb@xt@.01(4.18)

So, then, we have

 A(t)≤(1+αCPϵ√βH∗)E(t),A(t)≥(1−αCPϵ√βH∗)E(t), \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

 0<α≤α2≡2C∗1+χ, \hb@xt@.01(4.20)

where

 χ=(2+C2P(1+ϵC∗)2βH∗). \hb@xt@.01(4.21)

Then

 B(t)≥αE(t). \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

 δ=β√H∗ϵCP(1+ϵC∗),

and calculating that

 CP(1+ϵC∗)ϵ√H∗δ=C2P(1+ϵC∗)2βH∗,

we have that

 B(t)≥[2C∗−α(2+C2P(1+ϵC∗)2βH∗)]12∥uh,t∥21H+αβ2ϵ2∥∥∇⋅uDh∥∥2=(2C∗−αχ)12∥uh,t∥21H+αβ2ϵ2∥∥∇⋅uDh∥∥2. \hb@xt@.01(4.24)

We let be the solution to

 2c∗−α2χ=α2,

so that

 α2≡2C∗1+χ. \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

 2C∗−αχ≥2C∗−α2χ=α2≥α,

and so we also have

 B(t)≥αE(t). \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

 E(t)≤3E(0)e−2α3t. \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

 A′(t)+2α3A(t)≤0, \hb@xt@.01(4.28)

so that

 A(t)≤A(0)e−2α3t. \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

 (CHuh,vh)+1ϵ(fHu⊥h,vh)−βϵ2(ηh,∇⋅vh)=0(∇⋅uh,wh)=0 \hb@xt@.01(4.30)

It is easy to see that this problem is coercive on the divergence-free subspaces and thus is well-posed. Hence, with zero forcing, both and equal zero is the only solution. The zero-energy 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 long-time stability. The same techniques as before now lead to

 A′(t)+B(t)=(˜F,uh,t+αuDh). \hb@xt@.01(4.31)
###### Theorem 4.7

For any and

 Kα≡12[1+α2C2Pϵ2βH2∗], \hb@xt@.01(4.32)

we have the bound

 E(t)≤3e−α3tE(0)+Kαα∫t0eα3(s−t)∥∥˜F∥∥2Hds. \hb@xt@.01(4.33)

Proof. We bound the right-hand side of (LABEL:eq:ode2) by

 (˜F,uh,t+αuDh)≤∥∥˜F∥∥H∥uh,t∥1H+αCP∥∥˜F∥∥H∥∥∇⋅uDh∥∥1H≤[H∗2δ1+αCP2δ2]∥∥˜F∥∥2H+δ12∥uh,t∥1H+αCPδ22H∗∥∥∇⋅uDh∥∥2 \hb@xt@.01(4.34)

We put to find

 (˜F,uh,t+αuDh)≤1δ1Kα∥∥˜F∥∥2H+δ1E(t). \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

 A′(t)+αE(t)≤Kαδ1∥∥˜F∥∥2H+δ1E(t). \hb@xt@.01(4.37)

At this point, we specify so that, with (LABEL:eq:AequivE) we have

 A′(t)+α3A(t)≤Kαα∥∥˜F∥∥2H. \hb@xt@.01(4.38)

This leads to the bound on

 A(t)≤e−α3tA(0)+Kαα∫t0eα3(s−t)∥∥˜F∥∥2Hds. \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 long-time stability of the system. Our stability result also shows that the semidiscrete method captures the three-way geotryptic balance between Coriolis, pressure gradients, and forcing. Moreover, we also can demonstrate that “spin-up”, the process by which in practice tide models are started from an arbitrary initial condition and run until they approach their long-term 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

 χ≡Πu−u,ρ≡πη−η,θh≡Πu−uh,ζh≡πη−ηh. \hb@xt@.01(5.1)

The projections and satisfy the first-order system

 \hb@xt@.01(5.2)

Subtracting the discrete equation (LABEL:eq:discrete_mixed) from this gives

 (1Hθh,t,vh)+1ϵ(fHθ⊥h,vh)−βϵ2(ζh,∇⋅vh)+(CHθh,vh)=(fϵHχ+1Hχt+CHχ,wh),(ζh,t,wh)+(∇⋅θh,wh)=0. \hb@xt@.01(5.3)

By choosing the initial conditions for the discrete problem as and , the initial conditions for these error equations are

 θh(⋅,0)=0,ηh(⋅,0)=0. \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 ,

 12∥θh(⋅,t)∥21H+β2ϵ2∥ζh(⋅,t)∥2≤C2πh2(k+σ)C∗H∗∫t01ϵ|u(⋅,s)|2k+σ+|ut(⋅,s)|2k+σ+C∗|u(⋅,s)|2k+σds. \hb@xt@.01(5.5)

Proof. We apply Proposition LABEL:prop:firstorderstability to (LABEL:eq:firstordererror) to find

 12∥θh(⋅,t)∥21H+β2ϵ2∥ζh(⋅,t)∥2≤12C∗∫t0∥∥∥fϵHχ+1Hχt+CHχ(⋅,s)∥∥∥2Hds. \hb@xt@.01(5.6)

Note that for any ,

 ∥∥∥1Hg∥∥∥2H=∫ΩH(1H|g|)2dx=∫Ω1H|g|2dx=∥g∥21H.

Using this, that , and norm equivalence bounds the right-hand side above by

 1C∗H∗∫t0∥∥∥fϵχ(⋅,s)∥∥∥2+∥χt(⋅,s)∥2+∥Cχ(⋅,s)∥2ds≤1C∗H∗∫t01ϵ∥χ(⋅,s)∥2+∥χt(⋅,s)∥2+C∗∥χ(⋅,s)∥2ds

and the approximation estimate (LABEL:eq:PiL2) finishes the proof.

Since

 12∥(u−uh)∥21H+β2ϵ2∥η−ηh∥2≤∥ρ∥21H+∥ζh∥21H+β2ϵ2∥χ∥2+β2ϵ2∥θ∥2,

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

 12∥(u−uh)(⋅,t)∥21H+β2ϵ2∥(η−ηh)(⋅,t)∥2≤C2Πh2(k+σ)H∗|u(⋅,t)|2k+σ+C2πβh2kϵ2|η(⋅,t)|2k+2C2πh2(k+σ)C∗H∗∫t0|ut(⋅,s)|2k+σ+C∗|u(⋅,s)|2k+σds. \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 second-order system to obtain error estimates for the time derivative and divergence of the momentum. The projection satisfies the perturbed equation

 (1HΠutt,vh)+1ϵ(fHΠu⊥t,vh)+βϵ2(∇⋅Πu,∇⋅vh)+(CHΠut,vh)=(1Hχtt,vh)+1ϵ(1Hχ⊥t,vh)+(CHχt,vh)+(˜F,vh). \hb@xt@.01(5.8)

As in the first-order case, we have , and subtracting (LABEL:eq:secondorderdiscrete) from (LABEL:eq:secondorderprojeq) gives

 (1Hθh,tt,vh)+1ϵ(fHθ⊥h,t,vh)+βϵ2(∇⋅θh,∇⋅vh)+(CHθh,t,vh)=(1Hχtt,vh)+1ϵ(fHχ⊥t,vh)+(