FINITE VOLUME SIMULATION OF THE ... - Semantic Scholar

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT IN A ROTATING SHALLOW WATER SYSTEM∗ ´ ´ † MANUEL J. CASTRO , JUAN ANTONIO LOPEZ , AND CARLOS PARES Abstract. The goal of this article is to simulate rotating flows of shallow layers of fluid by means of finite volume numerical schemes. More precisely, we focus on the simulation of the geostrophic adjustment phenomenon. As spatial discretization, a first order Roe-type method and some higher order extensions are developed. The time discretization is designed in order to provide suitable approximations of inertial oscillations, taking into account the hamiltonian structure of the system for these solutions. The numerical dispersion laws and the wave amplifications of the schemes are studied and their well-balanced properties are analyzed. Finally, some numerical experiments for 1d and 2d problems are shown. Key words. Shallow water equations, geostrophic adjustment, well-balanced schemes, highorder methods, finite volume methods, Roe’s methods. AMS subject classifications. Primary, 74S10, 65M06; Secondary, 35L60, 35L65

1. Introduction. This paper is devoted to the development of finite volume numerical schemes that correctly simulate the geostrophic adjustment in rotating flows of shallow layers of fluid. This process, in which the Coriolis force balances the pressure forces, involves two different types of motion: high-frequency inertia gravity waves and low-frequency quasi-geostrophic motion. Therefore, in order to accurately simulate this phenomenon, a numerical method has to be able to handle these two kinds of motions. When rotation is present the shallow water system is dispersive: waves of different wave-length travel at different velocities. In order to correctly simulate the highfrequency inertia gravity waves, the discrete dispersion law has to be close enough to the exact one. Some comparisons between the discrete and the exact dispersion laws can be found in [2] for finite difference schemes and in [17], [18], [19] for finite element methods. Since quasi-geostrophic flows are close to equilibrium, the ability of the numerical schemes to approximate stationary solutions is also a fundamental aspect. This issue appears already for non rotating shallow water systems with nontrivial bottom topography. In [3] the condition called C-property was introduced: a scheme is said to satisfy this condition if it solves correctly the steady-state solutions corresponding to water at rest. This idea of constructing numerical schemes that preserve some equilibria, which are called in general well-balanced schemes, has been studied and extended by many authors: see [4] for a review on this topic. The design of wellbalanced high-order methods is nontrivial and some progress has been made recently: see [7], [10], [23], [24]. The well-balanced properties of finite volume schemes in the presence of rotation has been previously analyzed in [5] and [22]. The goal of this article is to obtain high-order numerical schemes for rotating shallow-water system with correct dispersion laws and good well-balanced properties. The article is organized as follows: in Section 2 a 1d simplified model is considered for ∗ This

research has been partially supported by the Spanish Government Research project MTM2006-08075. † Departamento de An´ alisis Matem´ atico. Facultad de Ciencias. Universidad de M´ alaga. 29071 M´ alaga, Spain. ([email protected], [email protected], [email protected]). The second author’s research was supported by a FPU predoctoral grant (MEC). 1

2

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

which the description of the numerical schemes and the analysis of their properties is easier. In Section 3, a first order semi-discrete numerical scheme for this simplified model is presented, which is based on a generalized Roe linearization (see [30]). The dispersion law of the method is calculated and compared with the exact one. Next, a semi-discrete high-order extension of the Roe method, based on the use of a reconstruction operator is presented (see [7], [26]). The dispersion laws corresponding to a second and a third order numerical schemes based on the use of WENO reconstrucions are analyzed. In Section 5, the discretization in time of the previously introduced semi-discrete methods is considered. It will be seen that, for low order methods, in order to obtain good dispersion relations, the time discretization has to preserve, in a sense to be precised, the Hamiltonian structure induced by the Coriolis term. Again, the dispersion laws corresponding to the full-discretized methods are calculated and compared with the exact one. In Section 6, the well-balanced properties of the numerical schemes are studied. Next, some numerical tests are presented in order to verify in practice these properties, as well as the ability of the methods to simulate the geostrophic adjustment for 1d flows. In Section 8, the 2d case is considered: the extensions of the schemes to this case are briefly described and a simulation of a 2d geostrophic adjustment is presented. Finally, in order to examine how the introduced schemes behave when the Coriolis parameter varies with latitude, they are applied to a numerical test proposed in [19]. 2. The 1d model. We consider the rotating shallow water equations with no dependence on the along-front coordinate y: Wt + F (W )x = S1 (W ) + S2 (W )Hx ,

x ∈ R, t > 0,

(2.1)

where



 h W =  q1  , q2





q1

  2  q1 1 2 F (W ) =   h + 2 gh   q1 q2 h

 0 S1 (W ) =  f q2  , −f q1



   ,   



 0 S2 (W ) =  gh  . 0

Here, h is the depth of water, H is the bottom depth from a fixed reference level, f is the Coriolis parameter, and g is the gravity acceleration. Variables q1 (x, t) and q2 (x, t) represent the mass flow in the x and y directions respectively, i.e. q1 = hu and q2 = hv, where u and v are the horizontal components of the averaged velocity. As it happens for the bed slope source term, a centered discretization of the Coriolis term, either implicit or explicit, may produce a numerical scheme which is not well-balanced. In order to correctly upwind this source term, it will be interpreted

3

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

as a nonconservative product. To do this, we introduce the function σ(x) = x and use the identity S1 (W ) = S1 (W )σx .

(2.2)

Adding the equations σt = 0 and Ht = 0, the system can be rewritten as follows: ft + A( eW f )W fx = 0, W

where

x ∈ R, t > 0,

(2.3)



 W f =  H , W σ

(2.4)

eW f ) is the 5 × 5 matrix whose block structure is and A(   A(W ) −S1 (W ) −S2 (W ) eW f) =  0 0 0 , A( 0 0 0 0 000

(2.5)

where A(W ) is the jacobian matrix of F (W ):   0 1 0 A(W ) =  −u2 + gh 2u 0  . −uv v u f take values in the set: The states W

f = [h, q1 , q2 , H, σ]T ∈ R5 , Ω = {W

(2.6)

h > 0},

eW f ) are: The eigenvalues of A(

f ) = u − c, λ2 (W f ) = u, λ3 (W f ) = u + c, λ4,5 (W f ) = 0, λ1 (W

p where c = gh, with associated eigenvectors:    1   λi     e2 (W f) =  ei (W f ) =  v  , i=1,3, R R      0  0 

  e f R4 (W ) =   

where Fr =

1 0 v 1 − Fr2 0



  ,  



0 0 1 0 0

fv g



  ,  

    0 e5 (W f) =  R  f v2  − f h(1 − Fr2 )   g  0 1 − Fr2

|u| is the Froude number. c

(2.7)



     ,    

(2.8)

4

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

The system is hyperbolic provided that c > 0 and Fr 6= 1. Furthermore, the e1 (W f ) and R e3 (W f ) are genuinely nonlinear: characteristic fields R ej (W f ) · ∇λj (W f ) 6= 0, R

f ∈ Ω, j = 1, 3, ∀W

e2 (W f ), R e4 (W f ) and R e5 (W f ) are linearly degenerate: while R ej (W f ) · ∇λj (W f ) = 0, R

f ∈ Ω, j = 2, 4, 5. ∀W

3. First order space discretization. Let us consider (2.1) with initial condition 

 h0 (x) W (x, 0) = W0 (x) =  q1,0 (x)  , q2,0 (x)

x ∈ R.

(3.1)

To discretize in space the system, it is first rewritten in the form (2.3) with initial conditions:   h0 (x)  q1,0 (x)    f (x, 0) =  q2,0 (x)  , x ∈ R. W (3.2)    H(x)  x Next, the computational domain is decomposed by using computing cells Ii , Ii = [xi−1/2 , xi+1/2 ]. For simplicity, let us assume that the cells have constant size ∆x and that xi+ 12 = i∆x. We define xi = (i − 1/2)∆x as the center of the cell Ii . We denote fi (t) the approximation of the cell averages of the exact solution provided by the by W numerical scheme, that is, Z xi+1/2 1 f (x, t) dx. fi (t) ∼ W W = ∆x xi−1/2 The initial conditions for the numerical scheme are given by: Z xi+1/2 fi (0) = 1 f0 (x) dx, ∀i. W W ∆x xi−1/2 fi (t): The following notation will be used for the components of W   hi (t)    q1,i (t)  Wi (t)   fi (t) =  Hi (t)  =  q2,i (t)  . W    Hi (t)  σi (t) σi (t) eW f ) in the sense The spatial discretization is based on a Roe linearization of A( defined in [30]: first a family of paths Ψ in the set Ω has to be chosen, i.e. a locally Lipschitz continuous map Ψ : [0, 1] × Ω × Ω → Ω satisfying fL , W fR ) = W fL , Ψ(0; W

fL , W fR ) = W fR , Ψ(1; W

(3.3)

5

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

and f, W f) = W f, Ψ(s; W

(3.4)

as well as certain regularity conditions (we refer to [9] for further details). Next, an itermediate matrix Aei+1/2 (t) is chosen at every inter-cell xi+1/2 at every time t that satisfies: 1. Aei+1/2 (t) is diagonalizable. fi (t) = W fi+1 (t) = W f , then Aei+1/2 (t) = A( eW f ). 2. If W 3. The following equality is satisfied:   fi+1 (t) − W fi (t) = Aei+1/2 (t)· W Z

1

0

e fi (t), W fi+1 (t))) ∂Ψ (s; W fi (t), W fi+1 (t)) ds. A(Ψ(s; W ∂s

(3.5)

Then, the semi-discretized numerical scheme can be written as follows (see [7]):
fi − W fi−1 ) + Ae− fi+1 − W fi ) , f ′ = − 1 Ae+ · (W (W W i i−1/2 i+1/2 ∆x

(3.6)

where the dependency on t has been dropped for simplicity. As usual, e e± e −1 Ae± i+1/2 = Ki+1/2 Li+1/2 Ki+1/2 ,

e where Le± i+1/2 is the diagonal matrix whose coefficients are the eigenvalues of Ai+1/2 : i+1/2

λ1

i+1/2

, . . . , λ5

,

e i+1/2 is a 5 × 5 matrix whose columns are associated eigenvectors. and K In this work we shall consider the family of paths composed by the segments: fL , W fR ) = W fL + s(W fR − W fL ), Ψ(s; W

and the Roe linearization given by  Ai+1/2 Aei+1/2 =  0 0 0 000

Ai+1/2 

S1,i+1/2

−S1,i+1/2 0 0



0 =  −u2i+1/2 + ghi+1/2 −ui+1/2 vi+1/2 0

   hi+1 + hi  g = 2    0

    ,   

S2,i+1/2

(3.7)

 −S2,i+1/2 , 0 0 1

2ui+1/2 vi+1/2 



s ∈ [0, 1],

0 0 ui+1/2

(3.8) 

,

0

  q2,i + q2,i+1   f 2 =   q1,i + q1,i+1   −f 2

(3.9)



    .    

(3.10)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

6 where

p p √ √ vi hi + vi+1 hi+1 ui hi + ui+1 hi+1 hi + hi+1 p p hi+1/2 = , vi+1/2 = . , ui+1/2 = √ √ 2 hi + hi+1 hi + hi+1 (3.11) Notice that Ai+1/2 is the usual Roe matrix for the conservation law system obtained by neglecting the source terms in (2.1). It can be easily verified that (3.8) if a Roe linearization of (2.3) for the family of paths (3.7). For this particular choice of intermediate matrices, the equations for Hi (t) and σi (t) in (3.6) are trivial and thus: Z xi+1/2 1 Hi (t) = Hi = H(x) dx, ∀i, ∆x xi−1/2 σi (t) = xi ,

∀i,

while those for Wi (t) reduce to: 1 P + (Ai−1/2 (Wi − Wi−1 ) − Si−1/2 ) ∆x i−1/2
− +Pi+1/2 (Ai+1/2 (Wi+1 − Wi ) − Si+1/2 ) ,

Wi′ = −

where

Si+1/2 = S1,i+1/2 (Hi+1 − Hi ) + S2,i+1/2 ∆x, ± Pi+1/2 = Ki+1/2



I ± sign(Li+1/2 ) 2



−1 Ki+1/2

(3.12)

(3.13)

(3.14)

and Ki+1/2 is a matrix whose columns are the eigenvectors of Ai+1/2 . In the absence of rotation, (3.12) coincides with the Q-scheme of Roe upwinding the source term introduced in [3]. Let us study the dispersion law of this numerical scheme, i.e. the relation between the frequency and the wave number for wave solutions of the linearized scheme. This discrete dispersion law will be then compared to the exact one. ˆ q1 = 0, q2 = 0, To do this, we linearize (2.1) around the constant state h = h, under the assumption of a flat bottom:  ∂h ∂q1   + = 0,   ∂t ∂x      ∂q1 ˆ ∂h − f q2 = 0, (3.15) + gh  ∂t ∂x         ∂q2 + f q1 = 0. ∂t The wave solutions of (3.15) verify the following dispersion law:  2 ω = 1 + (Rd k)2 , f

(3.16)

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

7

q ˆ is the Rossby deformation radius, ω is the wave-frequency and k where Rd = g h/f is the wave number. We consider now the spatial semi-discretization of the linearized problem (3.15): Wi′ (t) = −


1 P + (A(Wi − Wi−1 ) − Si−1/2 ) + P − (A(Wi+1 − Wi ) − Si+1/2 ) . (3.17) ∆x

where: 

0 ˆ A =  gh 0

 1 0 0 0 , 0 0

P±

 q ˆ 0 g h 1 ±1/  1 q , =  ˆ  1 0  2 ± gh 0 0 1

and the numerical source term is given by:  Si+1/2

After some  dhi   =   dt        dq1,i     dt =                dq2,i = dt



0

  q2,i + q2,i+1  ∆x  f 2  =  q1,i + q1,i+1   −f ∆x 2

(3.18)



    .    

(3.19)

algebraic calculations, (3.17) can be rewritten as follows: −

Rd f 1 1 (q2,i−1 − q2,i+1 ), (−hi−1 + 2hi − hi+1 ) − (q1,i+1 − q1,i−1 ) + 2∆x 2∆x 4Rd

−

Rd2 f 2 Rd f (hi+1 − hi−1 ) − (−q1,i−1 + 2q1,i − q1,i+1 ) 2∆x 2∆x

f + (q2,i−1 + 2q2,i + q2,i+1 ), 4 f − (q1,i−1 + 2q1,i + q1,i+1 ). 4

(3.20)

We look for solutions of (3.20) of the form hi (t) = ϕh (t)eikxi ,

q1,i (t) = ϕ1 (t)eikxi ,

q2,i (t) = ϕ2 (t)eikxi .

(3.21)

Introducing these expressions in (3.20), the following linear system of differential equations can be easily obtained:  i i Rd f  sin(k∆x)ϕ2 , (1 − cos(k∆x))ϕh − sin(k∆x)ϕ1 − ϕ′h = −    ∆x ∆x 2R d      iR2 f 2 Rd f f ϕ′1 = − d sin(k∆x)ϕh − (1 − cos(k∆x))ϕ1 + (1 + cos(k∆x))ϕ2 ,  ∆x ∆x 2         ϕ′ = − f (1 + cos(k∆x))ϕ . 1 2 2 (3.22)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

8

The eigenvalues of the matrix of the linear system (3.22) are: α1 = 0, α2,3

Rd f =− (1 − cos(k∆x)) ± i ∆x

r

Rd2 f 2 f2 sin2 (k∆x) + (1 + cos(k∆x))2 . 2 ∆x 4 (3.23)

If we define Rd f (1 − cos(k∆x)), ω = µ= ∆x

r

Rd2 f 2 f2 sin2 (k∆x) + (1 + cos(k∆x))2 , 2 ∆x 4

(3.24)

all of the solutions (3.21) of (3.20) have terms of the form e−µt ei(kxi −ωt) and, therefore, the numerical dispersion law is given by  2 ω R2 1 = d2 sin2 (k∆x) + (1 + cos(k∆x))2 , f ∆x 4

(3.25)

Rd . ∆x In Figure 3.1 we compare the numerical and the exact dispersion laws in the interval 0 < k < π/∆x, with fixed parameter Rd /∆x = 2. Some important differences between both graphs can be observed. The wave-frequencies of the numerical and exact solutions are similar only for small values of the wave number k. Nevertheless, although waves with a greater wave number cannot be correctly represented, they are filtered by the scheme, as they appear multiplied by a factor e−µt , close to zero except for small wave numbers (see Figure 3.1). which only depends on the parameters k∆x and

Dispersion law

Wave amplification

7 Exact Roe

Exact Roe

6

1

5 0.8

ω/f

4 0.6 3 0.4 2 0.2

1

0

0

0.5 k ∆ x / π

1

0

0

0.5 k ∆ x / π

1

Fig. 3.1. Left: comparison of the exact (solid line) and the numerical (dashed line) dispersion law corresponding to the semi-discrete Roe scheme. Right: comparison of the ∆x wave amplification as a function of the wave number at time T0 = √ . Rd /∆x = 2. ˆ gh

4. High-order semi-discrete methods based on reconstruction of states. In this section we develop higher order extensions of the Roe method introduced in Section 2 in order to improve the dispersion properties. To do this, we consider a reconstruction operator of order s, that is, an operator that, associates to a sequence

9

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

fi } of cell values, two new sequences of reconstruction at the inter-cells {W ± }, {W i+1/2 f , i.e. so that if the cell values are the averages of a smooth function W fi = 1 W ∆x

then:

Z

Ii

f (x) dx, W

(4.1)

s f± f W i+1/2 = W (xi+1/2 ) + O(∆x ).

Reconstruction operators are usually built by means of some approximations functions Pei whose calculation is usually made by means of interpolation techniques based on the values at some of the neighbor cells (the stencil ): fi−l , . . . , W fi+r ). Pei (x) = Pei (x; W

Once these functions have been calculated, the reconstructed states are given by: f+ W i−1/2 =

lim

x→x+ i−1/2

Pei (x),

f− W i+1/2 =

lim x→x− i+1/2

Pei (x).

Let us assume that the approximation functions satisfy the following properties: if fi } satisfy (4.1) for a smooth function W f then the cell values {W f (x) + O(∆xs1 ), Pei (x) = W

d e f ′ (x) + O(∆xs2 ), Pi (x) = W dx

∀ x ∈ Ii ,

(4.2)

∀ x ∈ Ii .

(4.3)

The semi-discrete method based on a high-order extension of a generalized Roe method for the system (2.3) can be written as follows (see [7]):  fi′ (t) = − 1 Ae+ f + (t)−W f − (t))+Ae− f+ f− W (t)·(W i−1/2 i−1/2 i+1/2 (t)·(Wi+1/2 (t)−Wi+1/2 (t)) ∆x i−1/2 Z xi+1/2  e Peit (x)) d Peit (x) dx , (4.4) A( + dx xi−1/2 f ± (t)} are the reconstructed states associated to the sequence {W fi (t)}, where {W i+1/2 Pe t the approximation functions: i

fi−l (t), . . . , W fi+r (t)), Peit (x) = Pei (x; W

and Aei+1/2 (t), the generalized Roe matrix of problem (2.3) associated to the states f ± (t). In [7], it has been proved that this scheme is of order at least γ = W i+1/2 min(s, s1 + 1, s2 + 1). We will use the following notation for the components of the intermediate states and the approximation functions:     ± Wi+1/2 Pit   ± f± Hi+1/2  , Peit =  ptH,i  , W i+1/2 =  ± ptσ,i σi+1/2

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

10 where:

± Wi+1/2

 h± i+1/2   ± =  q1,i+1/2  , ± q2,i+1/2 



 pth,i Pit =  ptq1 ,i  . ptq2 ,i

Again, the equations for the components Hi (t) and σi (t) are trivial. If the functions Pit are assumed to be continuous, the equation for Wi (t) can be written as follows:
1  + + − Wi′ = − Pi−1/2 · Ai−1/2 (Wi−1/2 − Wi−1/2 ) − Si−1/2 ∆x
− + − + Pi+1/2 · Ai+1/2 (Wi+1/2 − Wi+1/2 ) − Si+1/2 Z xi+1/2 d (4.5) − + S1 (Pit (x)) ptH,i (x) dx + F (Wi+1/2 ) − F (Wi−1/2 )− dx xi−1/2 Z xi+1/2  d S2 (Pit (x)) ptσ,i (x) dx , − dx xi−1/2 where + − + − Si+1/2 = S1,i+1/2 (Hi+1/2 − Hi+1/2 ) + S2,i+1/2 (σi+1/2 − σi+1/2 ).

(4.6)

Let us assume that the reconstruction operator is exact for first degree polynomials, as it will be the case for all the operators considered in this article. As σi (t) = xi ,

∀i,

we have: ± ptσ,i (x) = x, σi+1/2 = xi+1/2 ,

∀t, i.

Therefore, (4.5) reduces to:
1  + + − Pi−1/2 · Ai−1/2 (Wi−1/2 − Wi−1/2 ) − Si−1/2 ∆x
− + − + Pi+1/2 · Ai+1/2 (Wi+1/2 − Wi+1/2 ) − Si+1/2 Z xi+1/2 d − + S1 (Pit (x)) ptH,i (x) dx + F (Wi+1/2 ) − F (Wi−1/2 ) − dx xi−1/2 Z xi+1/2  S2 (Pit (x)) dx , −

Wi′ = −

(4.7)

xi−1/2

where now: + − Si+1/2 = S1,i+1/2 (Hi+1/2 − Hi+1/2 ).

(4.8)

In the particular case of WENO reconstructions (see [29]), two approximation functions are calculated at every inter-cell: − Pei+1/2 (x) =

r−1 X

k=0

e k (x), ωk− Q

+ Pei−1/2 (x) =

r−1 X

k=0

e k (x), ωk+ Q

11

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

e k is the derivative of an interpolation polynomial based on the values of the where Q fi at the stencil sequence W Skr = {xi−k , ..., xi−k+r−1 }

and the weights ωk± satisfy ωk± ≥ 0,

r−1 X

ωk± = 1.

k=0

These weights are calculated so that, on the one hand, the reconstruction operator is of order 2r − 1 and, on the other hand, the weight ωk is near to zero when the data on the stencil Skr indicate the presence of a discontinuity. We denote this operator as r-WENO, and the resulting scheme as r-WENO-Roe. For the details about WENO interpolation, see [14], [20], [28], [29]. We consider the following approximation functions at the cells:  Pe + (x) if x ∈ [xi−1/2 , xi ), i−1/2 Pei (x) = (4.9) Pe − (x) if x ∈ (xi , xi+1/2 ]. i+1/2 These functions are in general discontinuous, as

f − = Pe + (xi ) 6= Pe− (xi ) = W f+. W i i i−1/2 i+1/2

Due to this fact, an extra term has to be added at the right hand side of (4.5): F (Wi+ ) − F (Wi− ). Nevertheless the jumps at xi are of order O(∆xr ) and, due to this, there is no need to consider Riemann problems at the center of the cells. In fact, in [7] an alternative approach was introduced in which the approximations functions Pei were continuous, but in practice the numerical results given by the two versions of the scheme were similar being the continuous approach more costly. WENO reconstructions satisfy (4.2) and (4.3) with s1 = r and s2 = r − 1. Therefore, the order of accuracy of the WENO-ROE method is r while for conservative problems its order is 2r − 1. 4.1. 2-WENO-Roe and 3-WENO-Roe dispersion law. The expresion of the numerical scheme (4.5) corresponding to the linearized equations (3.15) is the following: 1  + + − + − A (Wi−1/2 − Wi−1/2 ) + A− (Wi+1/2 − Wi+1/2 ) ∆x Z xi+1/2  − + S2 (Pit (x)) dx . +A(Wi+1/2 − Wi+ + Wi− − Wi−1/2 )−

Wi′ = −

xi−1/2

Here, A± = P ± A,

(4.10)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

12

where A, P ± are given by (3.18). To calculate the dispersion law for the semi-discrete 2-WENO-Roe method and the semi-discrete 3-WENO-Roe method the procedure of Section 3 is carried out. As the solutions of the linearized semi-discrete scheme considered in this procedure are regular, the values of the weights in the WENO reconstructions are those corresponding to smooth solutions, that is: ω0+ = 2/3, ω1+ = 1/3, ω0− = 1/3, ω1− = 2/3, and ω0+ = 3/10, ω1+ = 3/5, ω2+ = 1/10, ω0− = 1/10, ω1− = 3/5, ω2− = 3/10, for the 2 and the 3-WENO reconstructions, respectively. It can be shown that the solutions of the form (3.21) satisfy the following system of ordinary differential equations: ϕ′ = M ϕ, where ϕ = [ϕh , ϕ1 , ϕ2 ]T and M is the matrix whose coefficients mij are q  ˆ gh 1 4 1 + cos(2k∆x) − cos(k∆x) , m11 = − 2∆x 3 3 m12 = − m21

ˆ g hi =− ∆x



i ∆x



 4 1 sin(k∆x) − sin(2k∆x) , m13 = 0, 3 6

 4 1 sin(k∆x) − sin(2k∆x) , m22 = m11 , m23 = f, 3 6 m31 = 0, m32 = −m23 , m33 = 0,

for the 2-WENO-Roe method and q  ˆ 1 1 gh 1 1 m11 = − − cos(k∆x) + cos(2k∆x) − cos(3k∆x) , 2∆x 3 2 5 30 m12

i =− ∆x m21



 17 22 1 − sin(2k∆x) + sin(k∆x) + sin(3k∆x) , m13 = 0, 30 15 30

ˆ g hi =− ∆x

  17 22 1 − sin(2k∆x) + sin(k∆x) + sin(3k∆x) , 30 15 30

(4.11)

(4.12)

(4.13)

(4.14) (4.15)

(4.16)

(4.17)

(4.18)

m22 = m11 , m23 = f,

(4.19)

m31 = 0, m32 = −m23 m33 = 0,

(4.20)

for the 3-WENO-Roe method. The corresponding dispersion laws and wave amplifications are depicted in Figure 4.1 and compared to the exact ones and those corresponding to Roe: it can be clearly seen how the improvement increases with the accuracy of the scheme.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

Dispersion law

13

Wave amplification

4 Exact 3−WENO−Roe 2−WENO−Roe Roe

ω/f

3

1 0.8 0.6

2 Rd / ∆ x = 1

Rd / ∆ x = 1

0.4 1 0

0.2 0

0.5 k∆x/π

1

0

0

0.5 k∆x/π

1

Wave amplification

Dispersion law 40 1

ω/f

30

0.8

20

0.6

Rd / ∆ x = 10

Rd / ∆ x = 10 0.4

10 0

0.2 0

0.5 k∆x/π

1

0

0

0.5 k∆x/π

1

Fig. 4.1. Left: comparison of the exact and the numerical dispersion laws corresponding to the semi-discrete Roe, 2-WENO-Roe, and 3-WENO-Roe methods. Right: comparison of ∆x the wave amplifications as a function of the wave number at time T0 = √ . Rd /∆x = 1, ˆ gh

Rd /∆x = 10.

5. Time discretization. 5.1. First order schemes. The dispersion laws and wave amplifications obtained in the previous Section change dramatically when the system is also discretized in time. Therefore, the choice of a suitable time discretization is a crucial aspect for

14

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

the problem under study. For example, the use of the classical forward Euler scheme produces nonphysical amplifications of the inertial oscillations, which are the solutions of the following simplified system:    dq1 = f q2 ,   dt (5.1)  dq  2   = −f q1 , dt corresponding to H = const, h = const, q1 = q1 (t), q2 = q2 (t). To illustrate the inconvenience of using the explicit Euler method we consider, following [32], the system (5.1) in compact form as: ∂Q = −if Q, ∂t

(5.2)

where Q = q1 + q2 i, and write the method as follows: Qn+1 = (1 − if ∆t)Qn .

(5.3)

For an initial state Q0 , the numerical solution is Qn = λn Q0 ,

(5.4)

where λ = 1 − if ∆t. Observe that |λ| > 1, which implies the amplification of inertial oscillations. Notice that (5.1) is a Hamiltonian system, whose Hamiltonian, H(q1 , q2 ), is given by H(q1 , q2 ) = − 21 f q12 − 21 f q22 . As symplectic schemes deal correctly with these oscillations (see [1], [13] for details), the idea here is to consider, if possible, explicit time discretizations that reduce to a symplectic scheme whenever H = const, h = const, q1 = q1 (t), q2 = q2 (t). The fact that the Hamiltonian is separable (i.e., H(q1 , q2 ) = T (q1 )+ U(q2 )) allows us to design explicit first order symplectic schemes, as the partitioned Runge-Kutta method: q1n+1 = q1n + ∆tf q2n , (5.5) q2n+1 = q2n − ∆tf q1n+1 . In order to design a time discretization that reduces to this scheme for the simplified system (5.1), we write formally the space-discretization as follows: Wi′ (t) = X (i; W (t)), where W (t) represents the vector of approximations of the cell-averages at time t. We propose the following first order time discretization: (1)

Wi

= Win + ∆tBX (i; W n ), (1)

Win+1 = Wi

∀i,

+ ∆tCX (i, W (1) ),

(5.6) ∀i,

15

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT (1)

where W n = {Win }, W (1) = {Wi }, and 

1 0 B= 0 1 0 0

 0 0 , 0



0 0 C= 0 0 0 0

 0 0 . 1

(5.7)

Let us study the dispersion law of the first order scheme obtained by applying this time discretization to the semi-discrete Roe scheme (3.12)-(3.14). To do this, (3.20) is discretized in time by using (5.6). Then, we look for solutions of the form: hni = ϕnh eikxi ,

n q1,i = ϕn1 eikxi ,

n q2,i = ϕn2 eikxi .

(5.8)

Some straightforward calculation show that ϕn = [ϕnh , ϕn1 , ϕn2 ]T has to satisfy: ϕn+1 = M · ϕn ,

(5.9)

M = (I + ∆t · B · M ) · (I + ∆t · C · M ),

(5.10)

where:

with 

    M =    

−

Rd f (1 − cos(k∆x)) ∆x −

iRd2 f 2 sin(k∆x) ∆x 0

− −

i sin(k∆x) ∆x

Rd f (1 − cos(k∆x)) ∆x

f − (1 + cos(k∆x)) 2

 i sin(k∆x)  2Rd    f (1 + cos(k∆x))  . 2    0 −

(5.11) Associated to every eigenvalue α0 of M one can find discrete wave solutions whose phase frequency satisfy: arg(α0 ) arg(α0 ) Rd ω = = · , f f ∆t CF L ∆x

(5.12)

where arg(α0 ) represents the argument of the complex number α0 , and whose ampliq ˆ fication factor after a time ∆x/ g h is |α0 |1/CF L .

The eigenvalues of M have been calculated numerically. Figure 5.1 shows the dispersion law and the amplification factors for some values of the CFL parameter for Rd /∆x = 2. Values of Rd /∆x greater than 2 may produce a numerical dispersion law which is even closer to the exact one. Nevertheless, small values of this parameter lead to real eigenvalues or eigenvalues with modulus greater than one, which prevents the scheme from representing waves correctly. Observe that, for small CFL numbers, the dispersion law of the fully discretized scheme is close to the semi-discretized one; CFL values near 1 or near 0.5 produce a dispersion law which is very close to the exact one; finally CFL values between 0.5 and 1 can lead to numerical waves which are faster than the exact ones. Concerning the amplification factor, notice that waves are not damped for great CFL values, while short waves are filtered for the rest of the values.

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

16

Dispersion law

Wave amplification

10

1

5 CFL = 0.20 0

0

0.5 CFL = 0.20 0.5

1

10

0

0.5

1

ω/f

0

0

0.5

1

0

0.5

1

0.5

1

0.5

1

0.5 CFL = 0.50 0.5

1

10

0

0

1

5 CFL = 0.70

0

0.5 CFL = 0.70 0.5

1

10

0

0

1

5 CFL = 0.98 0

1

1

5 CFL = 0.50

0

0.5

0.5 CFL = 0.49

10

0

0

1

5 CFL = 0.49 0

0

0

0.5 CFL = 0.98 0.5

k∆ x / π

1

0

0

k∆ x / π

Fig. 5.1. Left: comparison of the exact (solid line) and the numerical (dashed line) dispersion law corresponding to the Roe method with time discretization (5.6). Right: q comˆ parison of the wave amplifications as a function of the wavenumber at time T0 = ∆x/ g h.

Rd /∆x = 2.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

17

5.2. Second order schemes. As T (q1 ) = − 12 f q12 and U(q2 ) = − 21 f q22 are quadratic functions, it is possible to construct a second order symplectic scheme for (5.1) based on a Nystr¨om method (see [13] for details): (1)

q2 = q2n − q2n+1

=

q2n

∆t n f q1 , 2

−

(1)

q1 = q1n +

(1) ∆tf q1 ,

q1n+1

=

q1n

∆t (1) f q2 , 2

+

(5.13)

(1) ∆tf q2 .

In order to design a numerical discretization of the 2-WENO-Roe semi-discrete method that reduces to this scheme for (5.1), we denote again by X (i; W ) the operator corresponding to the spatial discretization and consider the three steps algorithm: (1)

= Win +

(2)

= Wi

Wi

Wi

(1)

∆t BX (i; W n ), 2

+

∆t CX (i; W (1) ), 2

Win+1 = Win + ∆tX (i; W (2) ), where



1 0 B= 0 0 0 0

 0 0 , 1



0 0 C= 0 1 0 0

∀i, (5.14)

∀i, ∀i,  0 0 . 0

(5.15)

The solutions of the form (5.8) satisfy now (5.9), where M is now the matrix given by: M = I + ∆tM (I +

∆t ∆t CM )(I + BM ), 2 2

and M is the matrix whose coefficients are given by (4.12)-(4.15). Figure 5.2 shows the corresponding dispersion laws and the wave amplification factors. If, instead, the usual second order TVD Runge-Kutta method (see [12]) is used: (1)

Wi

= Win + ∆tX (i; W n ),

Win+1 =

∀i,

1 n 1 (1) ∆t W + Wi + X (i; W (1) ), 2 i 2 2

(5.16) ∀i,

the matrix M for system (5.9) would be M=

1 1 ∆t I + (I + M )2 . 2 2 2

The corresponding dispersion law for Rd /∆x = 2 is shown in Figure 5.3. Notice that the scheme produces some numerical waves which are faster than the exact ones, as the numerical dispersion law is above the exact one. Futhermore, certain waves can be amplified. In fact, if (5.16) is applied to (5.1) the numerical solutions have the form (5.4) with: r (∆tf )4 . (5.17) |λ| = 1 + 4 As |λ| > 1, inertial oscillations are amplified for any time step ∆t.

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

18

Dispersion law

Wave amplification

8 6 ω / f

1

Exact Semidiscrete 2−WENO−Roe Symplectic 2−WENO−Roe

4

0.5

CFL = 0.2

CFL = 0.2

2 0

0

0.2

0.4 0.6 k∆x/π

0.8

1

8

0

0

0.2

0.4 0.6 k∆x/π

0.8

1

0.4 0.6 k∆x/π

0.8

1

0.4 0.6 k∆x/π

0.8

1

1

ω / f

6 4

0.5

CFL = 0.7

CFL = 0.7

2 0

0

0.2

0.4 0.6 k∆x/π

0.8

1

8

0

0

0.2

1

ω / f

6 4

0.5

CFL = 0.9

CFL = 0.9

2 0

0

0.2

0.4 0.6 k∆x/π

0.8

1

0

0

0.2

Fig. 5.2. Left: comparison of the exact and the numerical dispersion laws corresponding to the semi-discrete 2-WENO-Roe scheme and the 2-WENO-Roe scheme with time discretization (5.14). Right: comparison of the wave amplification factors as a function of q ˆ the wavenumber at time T0 = ∆x/ g h. Rd /∆x = 2.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

Dispersion law

19

Wave amplification

7

1.4 Exact Semidiscrete 2−WENO−ROE 2−WENO−Roe / RK2

6

1.2

1

4

0.8

3

0.6

2

0.4

1

0.2

ω/f

5

0

0

0.2

0.4 0.6 k∆x/π

0.8

1

0

0

0.2

0.4 0.6 k∆x/π

0.8

1

Fig. 5.3. Left: comparison of the exact and the numerical dispersion laws corresponding to the semi-discrete 2-WENO-Roe scheme and the 2-WENO-Roe scheme with time discretization (5.16). Right: comparison of the wave amplification factors as a function of q ˆ the wavenumber at time T0 = ∆x/ g h. Rd /∆x = 2.

5.3. Third order schemes. It has not been possible for us to find explicit time discretizations for the 3-WENO-Roe scheme which are suitable for this problem and that reduce to third order symplectic schemes for (5.1): the presence of negative intermediate time steps lead to non-physical oscillations. We have thus considered the usual third order TVD Runge-Kutta method (see [12]): (1)

= Win + ∆tX (i; W n ),

(2)

=

Wi Wi

∀i,

3 n 1 (1) ∆t W + Wi + X (i; W (1) ), 4 i 4 4

Win+1 =

1 n 2 (2) 2∆t W + Wi + X (i; W (2) ), 3 i 3 3

∀i,

(5.18)

∀i.

By applying this scheme to (5.1), it can be easily shown that inertial oscillations are not amplified under the constraint: f 2 ∆t2 < 3.

(5.19)

Now the solutions of the form (5.8) satisfy (5.9), where M is the matrix given by:   1 2 1 3 M = I + (I + ∆tM ) I + (I + ∆tM )2 , 3 3 4 4 and M is the matrix whose coefficients are given by (4.16)-(4.20). Figure 5.4 shows the dispersion laws and the wave amplification factors corresponding to this scheme.

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

20

Dispersion law

Amplification

8

1

Exact Semidiscrete 3−WENO−Roe 3−WENO−Roe / RK3 4 CFL = 0.2 2

ω / f

6

0

0

0.2

0.4 0.6 k ∆x/π

0.8

0.5

1

8

0

CFL = 0.2

0

0.2

0.4 0.6 k∆x/π

0.8

1

0.4 0.6 k∆x/π

0.8

1

0.4 0.6 k∆x/π

0.8

1

1

ω / f

6 4

0.5

CFL = 0.7

CFL = 0.7

2 0

0

0.2

0.4 0.6 k ∆x/π

0.8

1

8

0

0

0.2

1

ω / f

6 4

0.5

CFL = 0.9

CFL = 0.9

2 0

0

0.2

0.4 0.6 k ∆x/π

0.8

1

0

0

0.2

Fig. 5.4. Left: comparison of the exact and the numerical dispersion laws corresponding to the semi-discrete 3-WENO-Roe scheme and the 3-WENO-Roe scheme with time discretization (5.18). Right: comparison of the wave amplification factors as a function of q ˆ the wavenumber at time T0 = ∆x/ g h. Rd /∆x = 2.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

21

6. Well-balanced properties. As it was mentioned in the Introduction, wellbalancing is a fundamental requirement for the numerical scheme in order to correctly simulate the geostrophic adjustment. It can be easily verified that the smooth stationary solutions of (2.1) are either of the form:  q = 0,  Z  1 f h(x) = H(x) + v(x)dx, (6.1)  g  q2 (x) = v(x)h(x), where v is an arbitrary smooth function, or   q1 = Q 6= 0, q2 (x) = h(−f x + V ),  h(x) is, for each x a root of the polynomial:

  p(h, x) = gh3 + (−f x + V )2 /2 − gH − C h2 + Q2 /2,

(6.2)

(6.3)

where C is a fixed real number. According to a general result stated in [25], a Roe scheme (3.6) based on a linearization that satisfies (3.5) for the family of segments (3.7) solves up to the second order every smooth stationary solution of (2.3). Moreover, it solves exactly all the f (x) satisfying the following property: stationary solutions W eW f (x), W f (y)) · (W f (y) − W f (x)) = 0, A(

∀x, y,

(6.4)

eW f (x), W f (y)) is the Roe matrix corresponding to the states W f (x), W f (y). In where A( effect, if this property is satisfied, it can be trivially checked that the vector given by f: the point values of W fi = W f (xi ), W

∀i,

is a critical point for the ODE system (3.6). f (x) of the form (6.1) with v ≡ const Let us verify that a stationary solutions W satisfies (6.4): some straightforward calculations allow us to obtain the equality:     f v/g 1     0  
 0  2  f f   W (y) − W (x) = H(y) − H(x)  v  + (y − x)  f v /g   , ∀x, y,  1    0 0 1 and the two vectors appearing in the last expression are eigenvectors associated to 0 eW f (x), W f (y)): see (2.8) and take into account that Fr = 0 in of the Roe matrix A( this case. Therefore the semi-discrete and fully discretized Roe schemes introduced in Sections 3 and 5 solve exactly this family of stationary solutions. In particular, they solve exactly water at rest solutions, i.e. solutions of the form (6.1) with v = 0. In fact, when q2 = 0, the numerical scheme coincides with the Q-scheme of Roe upwinding the source term introduced in [3] which is known to satisfy the C-property. On the other hand, the semi-discrete or fully discretized high-order extensions introduced in Sections 4 and 5 solve any arbitrary smooth stationary solution up to

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

22

the order of the method. Moreover, they solve exactly the stationary solutions for which the reconstruction operator is well-balanced in the following sense: Definition 6.1. The reconstruction operator is said to be well-balanced for a f (x) if: stationary solution W • the approximation functions {Pei (x)} associated to the sequence of cell averf are also stationary solutions of the system, i.e. ages of W e Pei (x)) d Pei (x) = 0, A( dx

∀x ∈ Ii , ∀i;

• the reconstructed states verify:

f+ f− Aei+1/2 (W i+1/2 − Wi+1/2 ) = 0.

Again, it can be easily verified that the vector of cell averages of a stationary f is a critical point for the ODE system (4.4) if the reconstruction operator solution W is well-balanced for that solution. A reconstruction operator which is well-balanced for solutions of the form (6.1) with v ≡ const can be obtained as follows: first, consider a reconstruction operator that, given a family of scalars cell values {wi }, produces a family of approximation functions: pi (x; wi−l , ·, wi+r ). Let us assume that this operator is exact for linear functions and that, given any sequence of cell values {wi } and a real number λ, the following equality holds: pi (x; λ wi−l , ·, λ wi+r ) = λ pi (x; wi−l , ·, wi+r ). Notice that reconstruction operators based on high-order polynomial interpolation satisfy both properties. Next, given a sequence [hi , q1,i , q2,i , Hi ]T of cell averages, consider the new sequence [hi , q1,i , q2,i , ηi ]T with ηi = hi − Hi . Then, apply the reconstruction operator to obtain the approximation functions: ph,i ,

pq1 ,i ,

pq2 ,i ,

pη,i .

Finally, define pH,i = ph,i − pη,i . Let us verify that this reconstruction operator is well-balanced for stationary solutions of the form (6.1) with v ≡ const. Notice that the cell averages of these solutions satisfy: q1,i = 0,

q2,i = v hi ,

ηi =

vf xi + C, g

∀i,

for some constant C. As a consequence: pq1 ,i ≡ 0,

pq2 ,i = vph,i ,

pi,η (x) =

vf x + C, g

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

23

and then, Pei = [ph,i , pq1 ,i , pq2 ,i , pH,i , x]T

is also a stationary solution of the form (6.1) with v ≡ const. Finally, it can be easily verified that the reconstructed states satisfy:   1  
 0  + − f+ f−  v , W − W = H − H  i+1/2 i+1/2 i+1/2 i+1/2   1  0 and the vector appearing in the right-hand side of this equality is again an eigenvector ei+1/2 associated to the eigenvalue 0. of A 7. Numerical experiments.

7.1. Inertia gravity waves. In order to illustrate the behavior of the schemes for inertia gravity waves, we consider the linearized system (3.15) with initial conditions:  Q if |x| ≤ a, h0 (x) = const, q1 (x, 0) = 0, q2 (x, 0) = (7.1) 0 if |x| > a, where Q is an arbitrary constant, and f = 10−4 rad/s,

Rd = 2, a

Rd = 5, 10. ∆x

(7.2)

In Figures 7.1, 7.2, 7.3, and 7.4 the numerical results are compared with the exact solution which is computed following [6]. Together with the numerical schemes discussed above, we have also considered a second order scheme based on a a MUSCL reconstruction operator (see [31]) and a third order scheme based on a hyperbolic reconstruction operator (see [21]). The analysis of the dispersion law for this latter scheme has been omitted due to the technical difficulties involved. Observe in Figures 7.2 and 7.3 that first and second order schemes based on symplectic time integrators capture better the wave velocities than those based on Runge-Kutta methods. Nevertheless, for second order schemes the differences diminish as the parameter Rd /∆x increases. Notice that, for small values of this parameter, the second order schemes do not capture correctly the geostrophic regime, which is not the case for third order schemes (see Figure 7.4). 7.2. Well-balanced properties. Two numerical tests are shown to verify in practice the well-balanced properties of the numerical schemes. In both cases, we have solved the system with the Roe method with time discretization (5.6), the 2WENO-Roe method with time discretization (5.14), and the 3-WENO-Roe method with time discretization (5.18). 7.2.1. Test 1. We consider (2.1) in the interval [−1/2, 1/2], with a flat bottom (H = 1), f = 2, g = 9.81. The initial condition is the stationary solution of the form (6.1) with v ≡ 1. Free boundary conditions are considered. We have considered four different meshes with 20, 40, 80, and 160 cells. The error in L1 norm between the initial condition and the numerical steady state solution reached by using the three schemes is zero up to machine accuracy. As en example, in Table 7.1 the absolute errors for the 3-WENO-Roe method are shown.

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

24

0.1

0.15 Roe / Symplectic, Rd / ∆ x = 10

Roe / Euler Rd / ∆ x = 10

Roe / Symplectic Rd, / ∆ x = 5

0.08

Roe / Euler, Rd / ∆ x = 5 0.1

Exact

Exact

0.06 0.05 ^ h−h

^ h−h

0.04

0.02

0

0

−0.05

−0.02 −0.1 −0.04

0

5

10

15 x / Rd

20

25

30

0

5

10

15 x / Rd

20

25

30

Fig. 7.1. Inertia gravity waves: free surface at time t=80 hours for the linear system (3.15) with initial conditions (7.1). Comparison of the exact and the numerical solutions obtained with the Roe method with time discretization (5.6) (left) or forward Euler (right).

0.1

0.1 MUSCL / Symplectic, Rd / ∆ x = 10 MUSCL / Symplectic, Rd / ∆ x = 5 Exact

0.08

MUSCL / RK2, Rd / ∆ x = 10 MUSCL / RK2, R / ∆ x = 5

0.08

d

Exact

0.04

0.04 ^ h−h

0.06

^ h−h

0.06

0.02

0.02

0

0

−0.02

−0.02

−0.04

0

5

10

15 x/R

20

25

−0.04

30

0

5

10

d

15 x/R

20

25

30

d

Fig. 7.2. Inertia gravity waves: free surface at time t=80 hours for the linear system (3.15) with initial conditions (7.1). Comparison of the exact and the numerical solutions obtained with a second order scheme based on a MUSCL reconstruction operator with time discretization (5.14) (left) or (5.16) (right).

0.1

0.1 2−WENO−Roe / Symplectic, Rd / ∆ x = 10

2−WENO−Roe / RK2, Rd / ∆ x = 10

2−WENO−Roe / Symplectic, Rd / ∆ x = 5

0.08

2−WENO−Roe / RK2, Rd / ∆ x = 5

0.08

Exact

Exact

0.04

0.04 ^ h−h

0.06

^ h−h

0.06

0.02

0.02

0

0

−0.02

−0.02

−0.04

0

5

10

15 x/R

d

20

25

30

−0.04

0

5

10

15 x/R

20

25

30

d

Fig. 7.3. Inertia gravity waves: free surface at time t=80 hours for the linear system (3.15) with initial conditions (7.1). Comparison of the exact and the numerical solutions obtained with the 2-WENO-Roe method with time discretization (5.14) (left) or (5.16) (right).

25

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

0.1

0.1 3−WENO−Roe / RK3, R / ∆ x = 10

Hyperbolic rec. / RK3, R / ∆ x = 10

d

d

3−WENO−Roe / RK3, R / ∆ x = 5

0.08

Hyperbolic rec. / RK3, R / ∆ x = 5

0.08

d

d

Exact

Exact

0.04

0.04 ^ h−h

0.06

^ h−h

0.06

0.02

0.02

0

0

−0.02

−0.02

−0.04

0

5

10

15 x/R

20

25

−0.04

30

0

5

10

d

15 x/R

20

25

30

d

Fig. 7.4. Inertia gravity waves: free surface at time t=80 hours for the linear system (3.15) with initial conditions (7.1). Comparison of the exact and the numerical solutions obtained with 3-WENO-Roe method with time discretization (5.18) (left) and a third order scheme based on a hyperbolic reconstruction operator with time discretization (5.18) (right).

n part. 20 40 80 160

error h 8.43e-16 1.19e-15 2.73e-15 5.48e-15

error q1 1.21e-15 1.84e-15 3.63e-15 8.10e-15

error q2 2.47e-15 3.24e-15 7.45e-15 1.42e-14

Table 7.1 Numerical results for the 3-WENO-Roe method with time discretization (5.18) in Test 1 at time t = 1.

7.2.2. Test 2. We consider now (2.1) in the interval [−5, 5], with a flat bottom (H = 1), f = 10, g = 1. The initial condition is the stationary solution of the form (6.1) with: v(x) =

gx −x2 e . 2

In Tables 7.2-7.4 the absolute errors in L1 norm between the initial condition and the numerical steady state solution are depicted for the three schemes using four different meshes with increasing number of cells. In every case, the order of accuracy agrees with the theoretical analysis. n part. 50 100 200 400

error h 1.14e-02 1.44e-03 4.39e0-4 1e-04

error q1 3.59e-04 8.23e-04 1.49e-04 4.33e-05

error q2 1.91e-04 4.31e-03 1.02e-03 2.54e-04

order h 2.55 2.14 2.13

order q1 2.12 2.09 2.15

order q2 2.14 2.07 2.00

Table 7.2 Numerical results for the Roe method with time discretization (5.6) at time t = 30 in Test 2.

7.3. Geostrophic adjustment simulations. Following [5] and [22], we consider system (2.1) with initial conditions: h0 (x) = H,

q1,0 (x) = 0,

q2,0 (x) = HV NL (x),

(7.3)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

26 n part. 50 100 200 400

error h 7.82e-01 1.83e-01 2.52e-02 3.22e-03

error q1 2.40e-03 5.91e-04 1.29e-04 2.16e-05

error q2 4.67e-01 1.19e-01 1.69e-02 2.15e-03

order h 2.09 2.86 2.96

order q1 2.02 2.19 2.57

order q2 1.97 2.81 2.97

Table 7.3 Numerical results for the 2-WENO-Roe method with time discretization (5.14) at time t = 30 in Test 2.

n part. 50 100 200 400

error h 8.62e-02 3.03e-03 7.20e-05 1.01e-05

error q1 2.01e-03 8.23e-03 4.91e-06 6.71e-07

error q2 6.12e-01 2.28e-02 8.68e-04 8.30e-05

order h 4.83 5.39 2.83

order q1 4.61 4.06 2.87

order q2 4.76 4.71 3.46

Table 7.4 Numerical results for the 3-WENO-Roe method with time discretization (5.18) at time t = 30 in Test 2.

where H, V are two positive constants and NL is the function given by: NL (x) = We write the system           

where:

(1 + tanh(4x/L + 2))(1 − tanh(4x/L − 2)) . (1 + tanh(2))2 in adimensional form:

∂h ∂q1 + = 0, ∂t ∂x   1 ∂ q12 1 Bu 2 ∂q1 = h q2 , + + 2 ∂t ∂x h 2 Ro Ro         ∂q2 1 ∂  q1 q2    = − q1 , + ∂t ∂x h Ro V , Ro = fL

R2 Bu = 2d , L

√ gH Rd = , f

(7.4)

(7.5)

are, respectively, the Rossby number, the Burger number Bu , and the Rossby deformation radius. Figure 7.5 shows the approximation of the height h along an inertial period T = 2π f for Bu = 0.25, Ro = 1, CFL=0.9 and Rd /∆x = 50, obtained with the 2-WENO-Roe method and time discretization given by (5.14). Observe that fast inertia gravity waves are emitted as a consequence of the initial momentum disturbance and two shocks arise in the wave front. Although all the methods considered in previous Sections provide good results in this test, higher order methods are more suitable for long time simulations. A comparison between several methods of different order at time t = 10 · T is depicted in Figure 7.6. Figure 7.7 shows the geostrophic balanced between f v and ghx . Although the adjustment is rapidly achieved in the jet core, small oscillations with frequencies close

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

27

t / T = 0.0

t / T = 0.2

h

t / T = 0.4

t / T = 0.6

t / T = 0.8

t / T = 1.0

−8

−6

−4

−2

0

2

4

6

8

10

x / Rd

Fig. 7.5. Geostrophic adjustment: evolution of the free surface provided by the 2-WENORoe method with time discretization (5.14). Initial condition: (7.3). Parameters: Ro = 1, Bu = 0.25, CF L = 0.9, Rd /∆x = 50.

to f remain for a long time in the vicinity of the jet as it was noticed in [15], [5], or [22]. 8. The 2d model. Finally, we show some results for the bidimensional case. The system is now as follows:  ∂q2 ∂h ∂q1   + + = 0,   ∂t ∂x ∂y          ∂q ∂ q12 1 2 dH ∂  q1 q2  1 + + gh + , = f q2 + gh ∂t ∂x h 2 ∂y h dx           ∂ q22 dH ∂  q1 q2  1  ∂q2  + + + gh2 = −f q1 + gh .  ∂t ∂x h ∂y h 2 dy

(8.1)

As in the 1d case, (8.1) can be written as a non-conservative hyperbolic system by considering the functions σ1 (x, y) = x and σ2 (x, y) = y. Adding the equations ∂H = 0, ∂t

∂σ1 = 0, ∂t

∂σ2 = 0, ∂t

(8.1) can be written as ft + A f1 (W f )W fx + A f2 (W f )W fy = 0, W

(8.2)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

28

Roe 2−WENO−Roe 3−WENO−Roe Ref. solution

1.6 1.4

h

1.2 1 0.8 0.6 0.4 −8

−6

−4

−2

0 x / Rd

2

4

6

8

−2

0 x / Rd

2

4

6

8

0.8 Roe

0.6

2−WENO−Roe 3−WENO−Roe

0.4 q

2

Ref. solution

0.2 0 −0.2 −8

−6

−4

Fig. 7.6. Geostrophic adjustment: numerical solution at time t = 10T provided by the Roe method with time discretization (5.6), the 2-WENO-Roe method with time discretization (5.14), and the 3-WENO-Roe with time discretization (5.18). Initial condition: (7.3). Parameters: Ro = 1, Bu = 0.25, CF L = 0.9, Rd /∆x = 14.28.

where 

   f W =   

h q1 q2 H σ1 σ2

fi (W f ), i = 1, 2 are the 6 × 6 matrices : and A



   ,   

(8.3)

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

29

t / T = 22 gh

0.4

x

fv

0.2 0 −0.2 −20

−15

−10

−5

0 x/R

5

10

15

20

d

t / T = 22.25 0.4

gh

0.2

fv

x

0 −0.2 −20

−15

−10

−5

0 x/R

5

10

15

20

d

t / T = 22.5 ghx

0.4

fv

0.2 0 −0.2 −20

−15

−10

−5

0 x/R

5

10

15

20

d

Fig. 7.7. Geostrophic adjustment: geostrophic balance provided by the 2-WENO-Roe method with time discretization (5.14). Initial condition: (7.3). Parameters: Ro = 1, Bu = 0.25, CF L = 0.9, Rd /∆x = 50.



   e f A1 (W ) =    

0 c2 − u21 −u1 u2 0 0 0

1 2u1 u2 0 0 0

0 0 u1 0 0 0

0 −c2 0 0 0 0

0 − 21 f q2 1 2 f q1 0 0 0

0 0 0 0 0 0



   ,   

(8.4)

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

30 and



   e f A2 (W ) =    

0 −u1 u2 c2 − u22 0 0 0

0 u2 0 0 0 0

1 u1 2u2 0 0 0

0 0 −c2 0 0 0

0 0 0 − 21 f q2 0 12 f q1 0 0 0 0 0 0



   ,   

(8.5)

√ where ui = qi /h, i = 1, 2 and c = gh. In order to obtain a first order Roe scheme for system (8.2), following [8], a finite volume mesh of the computational domain is constructed and piecewise constant approximations of the solution are considered. These approximations are updated by considering the method of lines: at any time level, a family of projected Riemann Problems in the normal direction to each edge of the mesh is considered. These projected Riemann problems are then linearized by using a Roe linearization as in Section 3. The approximated solutions of these 1d linear Riemann problems are finally averaged in the cells to obtain the new piecewise constant approximation of the solution. High-order schemes are next obtained by extending to the 2d case the procedure developed in Section 4 for 1d problems: a reconstruction operator is considered, i.e. an operator that, given a family of constant values at the cells of the mesh, provides two functions at the edges, in such a manner that if the values at the cell are the averages of a regular function, then these edge functions are high-order approximations of of the traces of the regular function (see [8] for more details). In this section a MUSCL-type second order reconstruction is used (see [31]). Finally, the time discretizations considered are either the usual first and second order TVD Runge-Kutta methods, or the first and second order symplectic numerical schemes described in Sections 5.1 and 5.2. 8.1. Geostrophic adjustment simulation. We consider a test proposed in [16] corresponding to the following initial condition:   q √ √ ( λx)2 + (y/ λ)2 − Ri A0   , 1 − tanh  h(x, y, 0) = 1 + 2 RE

(8.6)

q1 (x, y, 0) = 0, q2 (x, y, 0) = 0,

using the parameters A0 = 0.5, λ = 2.5, RE = 0.1 and Ri = 1. The gravity and the Coriolis parameter have been fixed to g = 1 and f = 1 respectively. We consider a quadrilateral 400 × 400 mesh for the domain [−10, 10] × [−10, 10] and the CFL number is fixed to 0.9. Figure 8.1 shows the numerical results obtained with a MUSCL-type second order reconstruction and the second order symplectic time discretization described in Section 5.2. The initial elliptical mass imbalance evolves in an nonaxisymmetric way. Two shock waves propagate leaving behind an elevation which is slowly spinning clockwise. These results are similar to those appearing in [16]. Finally, similar results are obtained if the usual second order TVD Runge-Kutta method is used.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

Fig. 8.1. 2d model: evolution of the free surface for the initial condition (8.6).

31

32

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

8.2. Anticyclonic eddy propagation. We consider a test proposed in [19] whose objective is to study the influence of the numerical scheme on the propagation of planetary (Rossby) waves. Rossby waves are due to the variation of the Coriolis parameter with latitude and they propagate information and energy westward across the basin. In the large-scale oceanic circulation, they are responsible for the westward intensification associated with western boundary currents. Following [19], a lineari-

(a) Free surface

(b) Velocity field

Fig. 8.2. Anticyclonic eddy propagation: initial condition. The minimum and maximum value are specified at the bottom right corner of each panel (only the maximum value for the velocities).

ˆ q1 = 0, q2 = 0 under the assumption zation of (8.1) around the constant state h = h, of a flat bottom is considered:  ∂h ∂q1   + +   ∂t ∂x       ∂q1 ˆ ∂h + gh  ∂t ∂x       ∂q2   ˆ ∂h + gh  ∂t ∂y

∂q2 = 0, ∂y (8.7)

= f q2 , = −f q1 ,

where f = f0 + βy, that is, the β-plane approximation (see [11]) is considered. A rectangular basin of 2000 km × 1200 km is considered. As initial condition, a gaussian distribution of the free surface centered at the origin of the domain is prescribed together with a velocity field which is in geostrophic balance (f k × u = −g∇h): ˆ + η(x, y, 0), h(x, y, 0) = h

η(x, y, 0) = Ae−(x

2

+y 2 )/B 2

, 2

+y 2 )/B 2

q1 (x, y, 0) = h(x, y, 0)u1 (x, y, 0),

g y −(x u1 (x, y, 0) = 2A f0 +βy B2 e

q2 (x, y, 0) = h(x, y, 0)u2 (x, y, 0),

g x −(x u2 (x, y, 0) = −2A f0 +βy B2 e

2

,

+y 2 )/B 2

(8.8) ,

ˆ = 1.631 m, A = 0.95 m, B = 130 km. f0 and β are evaluated at 25o N, where h f0 = 6.1635 × 10−5 s−1 and βq = 2.0746 × 10−11 m−1 s−1 . The radius of deformation ˆ 0 ≈ 65 km. at the midbasin is thus Rd ≡ g h/f

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

(a) Free surface after four weeks of simulation

33

(b) Velocity field after four weeks of simulation

(c) Free surface after eight weeks of simulation (d) Velocity field after eight weeks of simulation Fig. 8.3. Anticyclonic eddy propagation (linear model): evolution of the free surface and the velocity field when a Roe first order scheme is used. The minimum and maximum value are specified at the bottom right corner of each panel (only the maximum value for the velocities).

Free boundary conditions are considered and ∆x = ∆y = 20 km. The CF L parameter is set to 0.9. The model is integrated for a time of eight weeks. The initial condition is shown in Figure 8.2. First, the initial condition adjusts to the β-plane balance of the model and, after this initial adjustment, the anticyclonic vortex is supposed to evolve westward. In Figure 8.3 we show the evolution of the surface elevation and the velocity field predicted by using the Roe method with the symplectic time discretization described in Section 5.1. In the figures, black bullets indicate the starting position of the centers of both the elevation and the velocity field at t = 0. The eddy has almost disappeared at the end of the simulation. Notice that the well-balanced properties of the numerical scheme for the case of a constant Coriolis parameter is not enough to properly capture the slow regime. In order to obtain good results the mesh size must be dramatically reduced or the accuracy of the scheme must be increased. For instance, Figure 8.4 shows the evolution of the surface elevation and the velocity field predicted by a second order scheme based on a MUSCL-type reconstruction and the symplectic time discretization described in Section 5.2. These results agree with those obtained in [19]. Similar results are obtained if the usual second order TVD Runge-Kutta method is used. Finally, the same experiment is reproduced for the fully non-linear system (8.1). Figure 8.5 shows the evolution of the free surface and the velocity field predicted by using a second order scheme based on a MUSCL second order reconstruction and the symplectic time discretization. REFERENCES

34

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

(a) Free surface after four weeks of simulation

(b) Velocity field after four weeks of simulation

(c) Free surface after eight weeks of simulation (d) Velocity field after eight weeks of simulation Fig. 8.4. Anticyclonic eddy propagation (linear model): evolution of the free surface and the velocity field when a MUSCL-type second order reconstruction is used. The minimum and maximum value are specified at the bottom right corner of each panel (only the maximum value for the velocities).

[1] L. Abia, M.J. Sanz-Serna, Partitioned Runge-Kutta methods for separable Hamiltonian problems. Math. Comput. 60 (1993), pp. 617-634. [2] A. Arakawa, V. R. Lamb, Computational design of the basic dynamical processes of the UCLA General Circulation Model. Methods in Computational Physics, 17. J. Chang(ed), Academic Press (1977), pp. 173-265. ´ dez, M. E. Va ´ zquez, Upwind methods for hyperbolic conservation laws with source [3] A. Bermu terms. Comput. & Fluids, 23 (1994), pp. 1049–1071. [4] F. Bouchut, Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources. Birkh¨ auser, 2004. [5] F. Bouchut, J. Le Sommer, V. Zeitlin, Frontal geostrophic adjustment and nonlinear wave phenomena in one dimensional rotating shallow water. Part 2: High-resolution numerical simulations. J. Fluid Mech. 514 (2004), pp. 35-63. [6] A. Cahn, An investigation of the free oscillations of a simple current system. J. Meteor. 2 (1945), pp. 113-119. es, High order finite volume schemes based on recons[7] M.J. Castro, J.M. Gallardo, C. Par´ truction of states for solving hyperbolic systems with non conservative products. Applications to shallow water systems. Math. Comp. 75 (2006), pp. 1103-1134. ´ ndez-Nieto, A.M. Ferreiro-Ferreiro, J.A. Garc´ıa-Rodr´ıguez, [8] M.J. Castro, E.D. Ferna es, High order extensions of Roe schemes for two dimensional nonconservative C. Par´ hyperbolic systems. Submitted. [9] G. Dal Maso, P.G. LeFloch, F. Murat, Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74 (1995), pp. 483-548. [10] J.M. Gallardo, C. Par´ es, M.J. Castro, On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas. J. Comput. Phys. 227 (2007), pp. 574–601. [11] A. E. Gill, Atmosphere-Ocean Dynamics. International Geophysics Series, Academic Press, 31. (1982). [12] S. Gottlieb, C.W. Shu, Total variation diminishing Runge-Kutta schemes. Math. Comp. 67 (1998), pp. 73–85.

FINITE VOLUME SIMULATION OF THE GEOSTROPHIC ADJUSTMENT

(a) Free surface after four weeks of simulation

35

(b) Velocity field after four weeks of simulation

(c) Free surface after eight weeks of simulation (d) Flow speed field after eight weeks of simulation Fig. 8.5. Anticyclonic eddy propagation (non-linear model): evolution of the free surface and the velocity field when the Roe scheme is used. The minimum and maximum value are specified at the bottom right corner of each panel (only the maximum value for the velocities).

[13] E. Hairer, PP. Norsett, G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems. Springer-Verlag. (1987). [14] G. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126 (1996), pp. 202-228. [15] A. C. Kuo, L. M. Polvani, Time-dependent fully nonlinear geostrophic adjustment. J. Phys. Ocean. 27 (1997), pp. 1614-1634. [16] A. C. Kuo, L. M. Polvani, Nonlinear geostrophic adjustment, cyclone/anticyclone asymmetry, and potential vorticiy rearrangement. Phys. of Fluids. 12 (2000), pp. 1087-1100. [17] D. Y. Le Roux, Dispersion relation analysis of the P1NC - P1 finite-element pair in shallowwater models. SIAM J. Sci. Comput. 27 (2005), pp. 394-414. [18] D. Y. Le Roux, V. Rostand, B. Pouliot, Analysis of numerically-induced oscillations in 2D finite-element shallow-water models, Part I: Inertia-gravity waves. SIAM J. Sci. Comput. 29 (2007), pp. 331-360. [19] D. Y. Le Roux, B. Pouliot, Analysis of numerically-induced oscillations in 2D finite-element shallow-water models, Part II: Free planetary waves. To appear in SIAM J. Sci. Comput. [20] X.D. Liu, S. Osher, T. Chan, Wighted essentially nonoscillatory schemes. J. Comput. Phys. 115 (1994), pp. 200-212. [21] A. Marquina, Local piecewise hyperbolic reconstuction of numerical fluxes for nonlinear scalar conservation laws. J. Sci. Comput. 15, (1994), pp. 892-915. ´c ˇova ´ -Medvido ´ va ´ , S. Noelle, M. Kraft, Well-balanced finite volume evolution [22] M. Luka Galerkin methods for the shallow water equations. J. Comput. Phys. 221 (2007), pp. 122147. [23] S. Noelle, N. Pankratz, G. Puppo, J.R. Natvig, Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows. J. Comput. Phys. 213 (2006), pp. 474-499. [24] S. Noelle, Y. Xing, C.-W. Shu, High order well-balanced finite volume WENO schemes for shallow water equation with moving water. J. Comput. Phys. 226 (2007), pp. 29-59. es, M. Castro, On the well-balance property of Roe’s method for nonconservative [25] C. Par´ hyperbolic systems. Aplications to Shallow-Water systems. ESAIM: M2AN. 38 (2004), pp.

36

´ ´ M.J. CASTRO, J.A. LOPEZ AND C. PARES

821-852. [26] C. Par´ es, Numerical methods for nonconservative hyperbolic systems. A theoretical framewok. SIAM J. Num. Anal. 44 (2006), pp. 300-321. [27] C. Rossby, On the mutual adjustment of pressure and velocity distributions in certain simple current systems II. J. Mar. Res. 1 (1938), pp. 239-263. [28] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes. J. Comput. Phys. 77 (1998), pp. 439-471. [29] C.-W. Shu, Essentially non-oscillatory and weighted essentiallly non-oscillatory schemes for hyperbolic conservation laws. ICASE Report 97-65, 1997. [30] I. Toumi, A weak formulation of Roe’s approximate Riemann solver. J. Comp. Phys. 102 (1992), pp. 360-373. [31] B. Van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method. J. Comput. Phys. 32 (1979), pp. 101-136. [32] J. Wang, M. Ikeda, Inertial stability and phase error of time Integration schemes in ocean general circulation models. Monthly Weather Review 125 (1996), pp. 2316-2327.