Mixed finite elements for global tide models

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2 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.

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 [5,6,10,22,[30][31][32][33]. 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 [8,26,29].
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 shallowwater 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 [12,27]. 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 parameterized drag term to model the effects of bottom friction, as described in [21]. 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 [34], whilst continuous finite element approaches have been used in many studies ( [11,18,23], for example). The lowest order Raviart-Thomas element for velocity combined with P 0 for height was proposed for coastal tidal modeling in [35]; 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 [21]. 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 timevarying tidal forcing, the equations have a time-varying attracting solution that all solutions converge to as t → ∞. 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.
Although our present work focuses squarely on the shallow water equations, we believe that many of our results will apply to other hyperbolic systems with damping. For example, the model we consider is just the damped acoustic wave equation plus the Coriolis term. Our techniques should extend to other settings where function spaces have discrete Helmholtz decompositions, most notably damped electromagnetics or elastodynamics.
Also, we point out that our main aim here is the theoretical analysis of the damped system. This is the first such analysis of which we are aware demonstrating the strong damping and hence optimal long-time error bounds. We do not assert whether similar results hold for other discretizations, just that they are unkown. While an extended experimental and/or theoretical study of these properties for a wide range of discretizations could be a fruitful project for the tide modelling community, it is beyond the scope of the present work. We do point out that the lowest-order Raviart-Thomas element has been previously proposed for tidal modeling [35], and our framework covers this case as well as the extension to higher-order methods and other mixed spaces.
The rest of this paper is organised as follows. In Sect. 2 we describe the finite element modelling framework which we will analyse. In Sect. 3 we provide some mathematical preliminaries. In Sect. 4 we derive energy stability estimates for the finite element tidal equations. In Sect. 5 we use these energy estimates to obtain error bounds for our numerical solution. Appendix A includes the discussion of embedded manifolds.

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 where u is the nondimensional two dimensional velocity field tangent to Ω, u ⊥ = (−u 2 , u 1 ) is the velocity rotated by π/2, η 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), f is the spatially-dependent 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), C is the (spatially varying) nondimensional drag coefficient and H 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 β 2 ∇η , we assume some F ∈ L 2 , giving our model as It also becomes useful to work in terms of the linearized momentum u = Hu rather than velocity. After making this substitution and dropping the tildes, we obtain A natural weak formulation of these equations is to seek u ∈ H (div) and η ∈ L 2 so that We now develop mixed discretizations with V h ⊂ H (div) and W h ⊂ L 2 . Conditions on the spaces are the commuting projection and divergence mapping V h onto W h . We define u h ⊂ V h and η h ⊂ W h as solutions of the discrete variational problem We will eventually obtain stronger estimates by working with an equivalent secondorder form. If we take the time derivative of the first equation in (5) and use the fact where F = F t . This is a restriction of which is the variational form of to the mixed finite element spaces. Note that here, we have taken the time derivative at the discrete level, and so this is an equivalent formulation to (5), making use of the compatible spaces, and so the problems associated with discretizing the equations in wave equation form (such as the methods described in [19,25]) do not arise. 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 [20], 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 (3) 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 H (div) proves proper long-term behavior of the model.

Mathematical preliminaries
For the velocity space V h , we will work with standard H (div) mixed finite element spaces on triangular elements, such as Raviart-Thomas (RT), Brezzi-Douglas-Marini (BDM), and Brezzi-Douglas-Fortin-Marini (BDFM) [3,4,28]. We label the lowestorder Raviart-Thomas space with index k = 1, 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 k = 1 as well. We will always take W h to consist of piecewise polynomials of degree k − 1, not constrained to be continuous between cells. In the case of domains with boundaries, we require the strong boundary condition u · n = 0 on all boundaries.
In the main part of this paper we shall present results assuming that the domain is a subset of R 2 , i.e. flat geometry. In the Appendix, we describe how to extend these results to the case of embedded surfaces in R 3 .
Throughout, we shall let · denote the standard L 2 norm. We will frequently work with weighted L 2 norms as well. For a positive-valued weight function κ, we define the weighted L 2 norm If there exist positive constants κ * and κ * such that 0 < κ * ≤ κ ≤ κ * < ∞ almost everywhere, then the weighted norm is equivalent to the standard L 2 norm by A Cauchy-Schwarz inequality holds for the weighted inner product, and we can also incorporate weights into Cauchy-Schwarz for the standard L 2 inner product by 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 V h we consider, the divergence maps V h onto W h . Also, the spaces of interest all have a projection, : H (div) → V h that commutes with the L 2 projection π into W h : for all w h ∈ W h and any u ∈ H (div). We have the error estimate when u ∈ H k+1 . Here, σ = 1 for the BDM spaces but σ = 0 for the RT or BDFM spaces. The projection also has an error estimate for the divergence for all the spaces of interest, whilst the pressure projection has the error estimate Here, C and C π are positive constants independent of u, η, and h, although not necessarily of the shapes of the elements in the mesh. We will utilize a Helmholtz decomposition of H (div) under a weighted inner product. For a very general treatment of such decompositions, we refer the reader to [2]. For each u ∈ V , there exist unique vectors u D and u S such that u = u D + u S , ∇ · u S = 0, and also 1 H u D , u S = 0. That is, H (div) is decomposed into the direct sum of solenoidal vectors, which we denote by and its orthogonal complement under the 1 H ·, · inner product, which we denote by Functions in N (∇·) ⊥ satisfy a generalized Poincaré-Friedrichs inequality, that there exists some C P such that We may also use norm equivalence to write this as Because our mixed spaces V h are contained in H (div), the same decompositions can be applied, and the Poincaré-Friedrichs inequality holds with a constant no larger than C p .

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 (5), and find that Since u ⊥ h · u h = 0 pointwise, we add these two equations together to find Hence, we have the following.

Proposition 1
In the absence of damping or forcing, the quantity is conserved exactly for all time.
Now suppose that F = 0 still but that 0 < C * ≤ C ≤ C < ∞ pointwise in Ω. The same considerations now lead to so that

Proposition 2 In the absence of forcing, but with
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 3 With nonzero F, we have that for all time t,
Proof We choose w h and v h as without forcing, and find that Cauchy-Schwarz, Young's inequality, and norm equivalence give 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 (6), we begin with vanishing forcing and damping terms, putting v h = u h,t to find which simplifies to so that the quantity is conserved exactly for all time.
If C is nonzero, we have that which implies that E(t) 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 u h = u D h + u S h in the 1 Hweighted decomposition. We let 0 < α be a scalar to be determined later and let the test function v in (6) and we rewrite the left-hand side so that We use the fact that This has the form of an ordinary differential equation where and By showing that for suitably chosen α, both A(t) and B(t) are comparable to E(t) defined in (28), we can obtain exponential damping of the energy.

Lemma 1 Suppose that
Then Proof We bound the term 1 H u h,t , u D h , with Cauchy-Schwarz, Poincare-Friedrichs (19), and weighted Young's inequality with δ = √ β : So, then, we have and the result follows thanks to the assumption (36).
Showing that B(t) is bounded above by a constant times E(t) is straightforward, but not needed for our damping results.

Lemma 2 Suppose that
where Then Proof We use Cauchy-Schwarz, the bounds 0 < C * ≤ C ≤ C * and | f | ≤ 1, and Young's inequality with weight δ > 0 to write 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 We let α 2 be the solution to If we pick α = α 2 , then we have the lower bound for B(t) is exactly α E(t). However, we are also constrained to pick α ≤ min{α 1 , α 2 } in order to guarantee that the lower bounds for A(t) is positive as well. If we have α ≤ α 2 , then and so we also have We combine these two lemmas to give our exponential damping result.
Theorem 1 Let α 1 and α 2 be defined by (36) and (40), respectively. Then, for any 0 < α ≤ min{α 1 , α 2 }, and any t > 0, we have Proof In light of (33), (42), and the lower bound in (37), we have that so that Using the upper and lower bounds of A in (37) 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 (5), 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 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 u h and η h equal zero is the only solution.
The zero-energy steady state then cannot have a nonzero solenoidal part. Moreover, the exponentially decay of u t toward zero forces u to reach its steady state quickly, driving both u D and u S toward zero at an exponential rate. Finally, since η t = −∇ · u almost everywhere, the exponential damping of ∇ · u 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 Theorem 2 For any 0 < α ≤ min{α 1 , α 2 } and we have the bound Proof We bound the right-hand side of (51) by We put δ 2 = βδ 1 H * αC P 2 to find This turns (51) into the differential inequality Using (42), we obtain At this point, we specify δ 1 = α 2 so that, with (37) we have This leads to the bound on A(t) Using (37) again gives the desired result.
These stability results have important implications for tidal computations. Theorem 2 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 (6) 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 η(x, t; t * ), u(x, t; t * ) for 0 > t * and t > t * as the solution starting from zero initial conditions at t * and define the global attracting solution as the limit as t * → −∞], 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 t large enough that the difference between the exact solution with those initial conditions and the global attracting solution is small and then letting h be small enough that the finite element solution approximates that exact solution well.

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 The projections u and πη satisfy the first-order system Subtracting the discrete equation (5) from this gives By choosing the initial conditions for the discrete problem as u h (·, 0) = u 0 and η h (·, 0) = πη 0 , the initial conditions for these error equations are We start with L 2 estimates for the height and momentum variables, based on the stability result for the first order system.
It is straightforward to get from here to a bound on the error Theorem 4 If the above assumptions hold, and also u t , u tt ∈ L ∞ ([0, t]; H k+1 (Ω)), then

Numerical results
In this section we present some numerical experiments that illustrate the estimates derived in the previous sections. In all cases the equations are discretized in time using the implicit midpoint rule. The domain is the unit sphere, centred on the origin, which is approximated using triangular elements arranged in an icosahedral mesh structure (see Appendix A for extensions of the results of this paper to embedded surfaces such as the sphere). All numerical results are obtained using the open source finite element library, Firedrake (http://www.firedrake.org). First, we verify the energy behavior in the absence of dissipation, i.e. C = 0. The variables were initialized with u = 0 and η = x yz, and the equations were solved with parameters = β = 0.1, f = 1, H = 1 + 0.1 exp(−x 2 ), and Δt = 0.01. The energy is conserved by the continuous-time spatial semi-discretization, and is quadratic. Since the implicit midpoint rule time-discretization preserves all quadratic invariants (see [24], for example), we expect exact energy conservation in this case; this was indeed observed as shown in Fig. 1. Upon introducing a positive dissipation constant C = 0.1, we observe both that the energy is monotonically decreasing (as implied by Proposition 2), and is scaling exponentially in time (as implied by Theorem 1). These results are also illustrated in Fig. 1.
Second, we verify the convergence results proved in Sect. 5. This was done by constructing a reference solution using the method of manufactured solutions, i.e. by choosing the solution where we have expressed the velocity in three dimensional coordinates even though it is constrained to remain tangential to the sphere. Here η and u are chosen to solve the continuity equation for η exactly, and F is then chosen so that the u equation is satisfied. We used the parameters = β = 0.1, f = H = 1, C = 1000, Ω = 2, and chose Δt = 10 −5 in order to isolate the error due to spatial discretization only. We ran the solutions until t = 0.3 and computed the time-averaged L 2 error for η. Plots are shown in Fig. 2; they confirm the expected first order convergence rate for V = RT1, Q = DG0, and the expected second order convergence rate for V = RT2, Q = DG1. Finally, we illustrate that this type of discretization excludes the possibility of spurious attracting solutions. In the case of the linear forced-dissipative tidal equations with time-dependent forcing, the continuous equations have the property that the solutions lose memory of the initial conditions exponentially quickly with timescale determined from C and the other parameters (and bounded by α in Theorem 1). As discussed among our stability results, any two solutions with different initial conditions should converge to the same solution as t → ∞. We illustrate this by randomly generating initial conditions for two solutions (u 1 , η 1 ) and (u 2 , η 2 ) with the same time-periodic forcing, Plot of the L 2 difference between two pairs of solutions (u 1 , η 1 ) and (u 2 , η 2 ) with different randomly generated initial conditions but the same forcing, as a function of time. As expected, the difference converges to zero (eventually with exponential rate) as t → ∞, demonstrating the absence of spurious solutions and measuring the difference between them as t → ∞. In performing this test, care must be taken to ensure that η 1 and η 2 both have zero mean as implied by the perturbative derivation of the linear equations (since the dissipation cannot influence the mean component). In this experiment, we used the parameters = β = 0.1, C = 10.0, Δt = 0.01 and we used an icosahedral mesh of the sphere at the fourth level of refinement. We indeed observed that the two solutions converge to each other exponentially quickly in the L 2 norm, as illustrated in Fig. 3.

Conclusions and future work
We have presented and analyzed mixed finite element methods for the linearized rotating shallow equations with forcing and linear drag terms. Our more delicate energy estimates rely on an equivalence between the first order form and a second order form, and this equivalence itself relies on fundamental properties of classical H (div) finite elements. In particular, our estimates show that the mixed spatial discretization accurately captures the long-term energy of the system, in which damping balances out forcing to prevent energy accumulation. Because of the linearity of the problem, our energy estimates also give rise to a priori error estimates that are optimal for Raviart-Thomas and Brezzi-Douglas-Fortin-Marini elements. Numerical results confirm both the stability and convergence theory given.
In the future, we hope to extend this work in several directions. First, we hope to study the more realistic quadratic damping model, which will require new techniques to handle the nonlinearity. Second, our estimates have only handled the semidiscrete case, and it is well-known that time-stepping schemes do not always preserve the right energy balances. Without damping or forcing, the implicit midpoint method preserves exact energy balance, and a symplectic Euler method will exactly conserve an approximate functional for linear problems. It remains to be seen how to give a rigorous fully discrete analysis, either including damping by a fractional step or fully implicit method. Finally, even explicit or symplectic time-stepping will require us to consider linear algebraic problems, as it is typically not possible to perform mass lumping for H (div) spaces on triangular meshes. Implicit methods will require additional care.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
In this case V h ⊂ V , W h ⊂ W . [15] dealt with this problem by constructing operators ι V h : V h → V and ι W h : W h → W such that where and π are projections from V to V h and W to W h respectively; these two operators commute with ∇· defined on M h . In particular, The estimates (14-16) then hold with ι V h • and ι W h •π replacing and π respectively, provided that the polynomial expansion of the element geometries in M h have at the same approximation order as V h and W h . There is also still a discrete Poincaré-Friedrichs inequality for V h . This means that all of our stability results from Sect. 4 hold in the manifold case, and it remains to deal with the error estimates. This is done by introducing further variables u h ∈ V h , η h ∈ W h satisfying This equation is of the form (5) but with a modified inner product on V h . Therefore, all of our stability estimates also hold for this modified equation. We split the error in u and η by writing where We can bound θ h and ζ h by applying Proposition 3 adapted to Eq. (73), i.e. by substituting v = ι V h v h into (4) and rearranging so that it takes the form of (73) with a forcing defined in terms of u, then subtracting (73). Similarly, θ h and ζ h may be bounded by rearranging Eq. (73) into the form of (4), then subtracting (4). Terms appear that are proportional to Id −J where and Id −J is the maximum of the operator norms of Id V h −J V h and Id W h −J W h . [15] showed that Id −J converges to zero as h → 0 with rate determined by the order of polynomial approximation in the isoparametric mapping. Hence we obtain a manifold version of Theorem 3, with u h and η h substituted by ι V h u h and ι W h η h respectively. Similar techniques lead to a manifold version of Theorem 4.