On the divergence-form finite-difference approximation to the momentum advection in curvilinear orthogonal coordinates

Ocean Modelling 1 (1999) 95±99

www.elsevier.com/locate/oadm

On the divergence-form ®nite-di€erence approximation to the momentum advection in curvilinear orthogonal coordinates Vladimir M. Kamenkovich * Department of Marine Science, Institute of Marine Sciences, The University of Southern Mississippi, Room 249, Bldg. 1103, Stennis Space Center, MS 39529, USA Received 3 June 1999; received in revised form 11 October 1999; accepted 18 October 1999

Abstract The traditional ®nite-di€erence schemes for the dynamical equations in curvilinear orthogonal coordinates have a basic ¯aw: they conserve only energy, but not momentum. A ®nite-di€erence approximation to these equations is suggested that conserves both energy and momentum. Ó 2000 Elsevier Science Ltd. All rights reserved.

The purpose of this note is to suggest a ®nite-di€erence scheme for the ocean-dynamics equations in curvilinear orthogonal coordinates that conserves both energy and momentum. By the conservation of momentum (energy) is meant that the ®nite-di€erence approximation to the advection of momentum (energy) does not generate any computational internal sources of momentum (energy). The curvilinear orthogonal coordinates are widely used in oceanology. We can refer to modeling the global ocean circulation in spherical coordinates or modeling the coastal ¯ows in a specially designed curvilinear coordinate system. The ®nite-di€erence schemes traditionally applied for approximating the dynamical equations in such coordinates conserve only energy, but not momentum. This basic ¯aw cannot be remedied if one considers the momentum equations in the component form. We will show that the required scheme can be constructed based on the approximation of the momentum equation in the vector form. For simplicity's sake we restrict ourselves to the analysis of the shallow water equations. From what follows it will be clear how to generalize the results for 3D case. Introduce curvilinear orthogonal coordinates n1 and n2 with the unit coordinate vectors a1 and a2 and metric scale factors m1 and m2 such that

*

Tel.: +1-228-688-3091; fax: +1-228-688-1121. E-mail address: [email protected] (V.M. Kamenkovich).

1463-5003/00/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 1 4 6 3 - 5 0 0 3 ( 9 9 ) 0 0 0 0 9 - 8

96

V.M. Kamenkovich / Ocean Modelling 1 (1999) 95±99

ds2 ˆ m21 …dn1 †2 ‡ m22 …dn2 †2 ;

…1†

where ds is an element of length. The velocity vector v is represented as v ˆ v1 a1 ‡ v2 a2 :

…2†

The shallow water equations in coordinates n1 and n2 have the following form: ov1 v1 ov1 v2 ov1 v2 om2 v1 v2 om1 g og ‡ ÿ 2 ‡ ÿ fv2 ˆ ÿ ; ‡ m1 on1 ot m1 on1 m2 on2 m1 m2 on1 m1 m2 on2 ov2 v1 ov2 v2 ov2 v2 om1 v1 v2 om2 g og ‡ ÿ 1 ‡ ‡ fv1 ˆ ÿ ; ‡ m2 on2 ot m1 on1 m2 on2 m1 m2 on2 m1 m2 on1   oh 1 o o ‡ …m2 v1 h† ‡ …m1 v2 h† ˆ 0; ot m1 m2 on1 on2

…3† …4† …5†

where f is the Coriolis parameter; g is the gravity acceleration; h ˆ H ‡ g, H is the depth, g is the sea surface elevation. The nonlinear terms in Eqs. (3) and (4) represent the momentum advection. By invoking the mass conservation equation (5) one commonly rewrites Eqs. (3) and (4) in the form:   o 1 o o gh og …hv1 † ‡ …m2 v1 hv1 † ‡ …m1 v2 hv1 † ÿ … f ‡ F †hv2 ˆ ÿ ; …6† ot m1 m2 on1 on2 m1 on1   o 1 o o gh og …hv2 † ‡ …m2 v1 hv2 † ‡ …m1 v2 hv2 † ‡ … f ‡ F †hv1 ˆ ÿ ; …7† ot m1 m2 on1 on2 m2 on2 where F ˆ

  1 om2 om1 v2 ÿ v1 m1 m2 on1 on2

…8†

The traditional ®nite-di€erence schemes are based on the approximation of the momentum Eqs. (6) and (7). Notice that only a part of momentum advection is represented here in the divergence form. The other part …ÿFhv2 and Fhv1 † caused by the curvature of the coordinate system n1 ; n2 is represented analogously to the Coriolis acceleration. Such a representation is useful in constructing the energy conserving scheme (Arakawa & Lamb, 1977), but it cannot ensure the conservation of momentum. It appears that in curvilinear coordinates Eqs. (6) and (7) always produce some internal sources of hv1 and hv2 conditioned by the momentum advection. In the di€erential formulation these sources are balanced because the advection of momentum vector does not generate any internal sources. But using the ®nite-di€erence approximation inevitably violates this balance and consequently produces some computational internal sources of momentum vector. Our task is to suggest a divergence-form ®nite-di€erence approximation to the whole expression of momentum R advection. The important point is that S hvi dS; i ˆ 1; 2; has nothing to do with the total momentum of a volume S  h. That is why, by considering Eqs. (6) and (7) separately, one fails to represent in a divergence form the sum of all nonlinear terms. In other words, one cannot consider separately the advection of hv1 and hv2 in curvilinear coordinates. Therefore the component forms of the

V.M. Kamenkovich / Ocean Modelling 1 (1999) 95±99

97

momentum equations (6) and (7) are not useful in the case considered and we rewrite the momentum equation in the vector form. Invoking the known formulas for oa1 =oni and oa2 =oni (see, for example, Batchelor, 1967; Appendix 2) and Eq. (5) yields the relevant ¯ux form of the momentum equation   o 1 o o …hv† ‡ …m2 v1 hv† ‡ …m1 v2 hv† ‡ f k  hv ˆ ÿghrg; …9† ot m1 m2 on1 on2 where k is the unit vertical vector. From Eq. (9) we see that the advection of momentum vector does not not produce any internal sources. The divergence form of the energy equation follows from Eqs. (9) and (5). It is derived by adding Eq. (9) multiplied scalarly by v and Eq. (5) multiplied by g. Recalling that h ˆ H ‡ g, we ®nd after some manipulations:    v
v i v
v i o v
v g2 1 o h o h m2 v1 h m1 v 2 h ‡ ‡g h ‡ gg ‡ ‡ gg ˆ 0: …10† 2 ot 2 m1 m2 on1 2 on2 2 We will formulate now a ®nite-di€erence scheme based on the approximation of Eqs. (9) and (5). We will require that the corresponding ®nite-di€erence equations should not produce: (1) any internal sources of momentum vector due to the momentum advection and (2) any internal sources of total energy. It should be recognized that because we are dealing with the vector form of momentum equation, we consider only those lattices, for which both velocity components are de®ned in the same grid points. Therefore widespread lattice C appears to be not appropriate for the construction of momentum conserving schemes in curvilinear coordinates (see the classi®cation of lattices in Mesinger & Arakawa, 1976, p. 47). Note that the de®nition of both horizontal velocity components in the same grid points is also widely applied in ocean modeling. The Bryan± Cox±Semtner global model (Bryan, 1969; Semtner, 1986) and all of the subsequent advanced versions of this model are using lattice B. We will start with lattice A. We have:      o 1 n1 n n2 n 1 2 …hv† ‡ m2 v 1 h v ‡ m1 v 2 h v ‡ f k  hv ˆ ÿghrg; …11† n1 n2 ot m1 m2      oh 1 n1 n2 ‡ m2 v1 h ‡ m1 v 2 h ˆ 0; …12† n1 n2 ot m1 m2 where rg ˆ

1 n1 1 gn1 a1 ‡ gnn22 a2 ; m1 m2

and the notation           1 d d 1 d d ni ani ˆ a ni ‡ ÿ a ni ÿ ; a ˆ a ni ‡ ‡ a ni ÿ d 2 2 2 2 2

…13†

…14†

is used. The similar approximations were mentioned in Grammeltvedt (1969) but for the momentum equations written in the component form and in the Cartesian rectilinear coordinates. In Eq. (14) a is an arbitrary grid function (scalar or vector) de®ned in points of lattice A with the grid

98

V.M. Kamenkovich / Ocean Modelling 1 (1999) 95±99

spacing d, and i ˆ 1; 2. Operators (14) are de®ned on the ®ner grid with the grid spacing d=2. The sought functions v and h in Eqs. (11) and (12) are de®ned in the points of lattice A. As is customary when the conservative properties of ®nite di€erence equations are considered, we maintain the time derivatives in Eqs. (11) and (12). To prove the conservative properties of Eqs. (11) and (12) we will use the following identities ÿ
v
v 1 Uni ‡ ‰ A…ni † ÿ A…ni ÿ d†Š; …15† v
v ni U n ˆ i 2 2d n

g W nii ‡ gnnii W ˆ

1 ‰ B…ni † ÿ B…ni ÿ d†Š; 2d

…16†

where U and W are arbitrary functions, A…ni † ˆ ‰v…ni †
v…ni ‡ d†ŠU …ni ‡ d=2† and B…ni † ˆ g…ni †W …ni ‡ d† ‡ g…ni ‡ d†W …ni †; and i ˆ 1; 2. In the derivation of Eqs. (15) and (16) we followed Grammeltvedt (1969), where similar identities were developed. Let us suppose that Eqs. (11) and (12) are considered on the entire plane n1 ; n2 . Then taking the sum over all of the grid points yields  X X o X hvm1 m2 d 2 ‡ …hrg†m1 m2 d 2 : …17† … f k  hv†m1 m2 d 2 ˆ ÿg ot We see that ®nite-di€erence momentum equation (11) does not produce any internal sources of momentum vector due to the momentum advection, similar to di€erential momentum equation (9). To prove the conservation of energy we, ®rst, multiply Eq. (11) scalarly by v and use the identity o o  v
v  v
v oh h ‡ ; …18† v
…hv† ˆ ot ot 2 2 ot Eqs. (15) and (12). Second, we multiply Eq. (12) by gg and use (16). Adding the results and summing over all of the grid points gives     o X v
v g2 2 h …19† ‡g m1 m2 d ˆ 0: ot 2 2 Thus, ®nite di€erence Eqs. (11) and (12) do not produce any internal sources of the total energy, similar to di€erential Eqs. (9) and (5) (see Eq. (10)). Notice that if the domain of integration is bounded, near-boundary grid points will provide some contribution to Eqs. (17) and (19). n n2 n1 2 For lattice B we should replace h and g in Eqs. (11) and (12) by …h † and …gn1 † , respectively. The derived conservative properties will remain valid and their development will be absolutely similar.

V.M. Kamenkovich / Ocean Modelling 1 (1999) 95±99

99

Acknowledgements This work has been supported by the National Science Foundation under grant OCE-96-33470. I am grateful to anonymous reviewers for a discussion.

References Arakawa, A., & Lamb, V. R. (1977). Computational design of the basic dynamical processes of the UCLA general circulation model. In Methods of computational physics, vol. 17, Academic Press, New York, pp. 173±265. Batchelor, G. K. (1967). An introduction to ¯uid dynamics. Cambridge University Press, Cambridge, p. 615. Bryan, K., 1969. A numerical method for the study of the circulation of the world ocean. J. Comput. Phys. 4, 347±376. Grammeltvedt, A., 1969. A survey of ®nite-di€erence schemes for the primitive equations for a barotropic ¯uid. Monthly Weather Rev. 97 (5), 384±404. Mesinger, F., & Arakawa, A. (1976) Numerical methods used in atmospheric models. GARP Publications Series No. 17, p. 64. Semtner, A. J. Jr., 1986. Finite-di€erence formulation of a world ocean model. In O'Brien, J. J. (Ed.), Advanced physical oceanographic numerical modelling, D. Reidel, Norwell, Mass., pp. 187±202.