Shallow Water Dynamics and Numerical Dispersion

Lecture 11 Shallow water dynamics and dispersion 7.1 Shallow water equations (Quick reference) The shallow water equations describe the dynamics of a hydrostatic, homoge­ nous fluid layer: �t u + u�x u + v�y u − f v + g�x � = 0 �t v + u�x v + v�y v + f u + g�y � = 0 �t � + �x (hu) + �y (hv) = 0 (7.1) (7.2) (7.3) where u and v are the two components of (horizontal) velocity in the di­ rections of x and y, � is the free-surface displacement from mean sea-level and h = H(x, y) + �(x, y, t) is the total fluid depth. H(x, y) is the resting fluid depth (or depth of to­ pography). f = 2� sin(�) is the Coriolis parameter and g is gravitational acceleration. Note that the continuity equation could equally be written �t h + �x (hu) + �y (hv) = 0 The shallow water equations conserve higher moments such as energy and potential vorticity. These properties are easier to obtain from the more con­ venient “vector invariant” form of the equations. Here, we make use o f the two relations u2 + v 2 ) − (�x v − �y u)v u�x u + v�y u = �x ( 2 u2 + v 2 ) + (�x v − �y u)u u�x v + v�y v = �y ( 2 These allow the momentum equations (7.1 and 7.2) to be written �t u − hqv + �x (g� + K) = 0 �t v + hqu + �y (g� + K) = 0 (7.4) (7.5) where 1 1 hq = f + ∂ = f + (�x v − �y u) and K = u2 + v 2 . 2 2 Here, q is the potential vorticity and K is the kinetic energy density. Taking the curl of equations (7.4) and (7.5) gives the vorticity equation �t (hq) + �x (hqu) + �y (hqv) = 0 which combined with the continuity equation gives the potential vorticity equation �t q + u�x q + v�y q = 0. Note that to get here, we made use of the property that curl of grad vanishes; we will need to make sure this happens in discrete analogs. The kinetic energy equation is obtained by multiplying the momentum equations by u and v respectively and adding the resulting equations: �t K + u�x (g� + K) + v�y (g� + K) = 0. Note that the (non-linear) Coriolis terms did not contribute to the kinetic energy equation and we should make sure that this happens in the numerical analog. Multiplying by h and adding the continuity equation multiplied by K and g� gives: 1 �t ( g� 2 + hK) + �x [(g� + K)hu] + �y [(g� + K)hv] = 0. 2 7.2 Inertia-gravity waves in 2-D We previously analyzed inertia-gravity wave propagation on one dimensional staggered grids. We’ll no repeat the procedure for the linear shallow water equations in two dimensions, which are: �t u − f v + g�x � = 0 �t v + f u + g�y � = 0 �t � + H�x u + H�y v = 0 The corresponding dispersion relation of the inertia-gravity waves is: φ 2 = f 2 + gH(k 2 + l2 ) where k and l are the mode numbers in the x and y directions respectively. A) B) uv u v η uv η uv v C) u η v u D) u uv v η v u Figure 7.1: Four natural arrangements of shallows water variables on quadri­ lateral grids. There are four arrangements of the variables on a regular quadrilateral grid (shown in Fig. 7.1) that treat the x and y directions symmetrically (i.e. a coordinate rotation of ±ω/2 does not change the stencil). Treating the time derivative as perfect, the corresponding models, using centered second order spatial differencing, are: • A grid: g βi � i = 0 �x g βj � j = 0 �t v + f u + �y H H �t � + βj v j = 0 βi ui + �x �y �t u − f v + • B grid: g βi � j = 0 �x g βj � i = 0 �t v + f u + �y H H �t � + βj v i = 0 β i uj + �x �y �t u − f v + • C grid: g βi � = 0 �x g βj � = 0 �t v + f uij + �y H H βi u + βj v = 0 �t � + �x �y �t u − f v ij + • D grid: g βi � ij = 0 �x g ij βj � ij = 0 �t v + f u + �y H H �t � + βj v ij = 0 βi uij + �x �y �t u − f v ij + i We saw in section 2.2.2 a pattern where the operators βi � and � always led to the expressions 2i sin k�x and cos k�x in the dispersion relation. Here, we 2 2 will use the following relations for the response functions: R(βi �) = 2i sin k�x = 2isk 2 l�y = 2isl 2 k�x i = ck R(� ) = cos 2 l�y j = cl R(� ) = cos 2 R(βj �) = 2i sin where sk , sl , ck and cl are short-hand for the trigonometric functions. We can now write down the dispersion relations for inertia-gravity waves (by inspection): • A grid: φ2 = f 2 + � or φ f 4gH 2 2 4gH 2 2 s c + s c �x2 k k �y 2 l l 2 = 1 + rx2 s2k c2k + ry2 s2l c2l • B grid: φ2 = f 2 + � or φ f 4gH 2 2 4gH 2 2 s c + sc �y 2 l k �x2 k l 2 = 1 + rx2 s2k c2l + ry2 s2l c2k • C grid: φ 2 = f 2 c2k c2l + � or φ f 4gH 2 4gH 2 s + s �x2 k �y 2 l 2 = c2k c2l + rx2 s2k + ry2 s2l • D grid: φ 2 = f 2 c2k c2l + or � φ f 4gH 2 2 2 4gH 2 2 2 s c c + sc c �x2 k k l �y 2 l k l 2 = (1 + rx2 s2k + ry2 s2l )c2k c2l where rx and ry are resolution parameters: rx = 2Ld 2Ld and ry = �x �y � and Ld = gH/f is the deformation radius. These dispersion relations are plotted in Fig. 7.2 and 7.3. We see that much of the analysis done in one dimension still holds, namely that the D grid has stationary modes regardless of the resolutions and that most grids have some group velocity sign errors. However, previously we concluded that the 1-D B grid did well but we see here that the frequency is underestimated for short two-dimensional waves at both high and low resolution. Indeed, the better grid at high resolution is arguably the A grid but even that has group velocity errors. At high resolution, the only grid that produces a monotonically increasing frequency (i.e. one with group velocity pointed in the right direction) is the C grid. B grid 1.2 1.4 1.1 1.2 ω/f ω/f A grid 1 0 1 0 1 0.5 0.5 s k 1 0 1 0.5 0.5 sl s k C grid 1 0 sl D grid 1 ω/f ω/f 1 0.5 0 0.5 0 0 1 0.5 0.5 s k 1 0 1 0.5 0.5 sl s k 1 0 sl Figure 7.2: Inertia-gravity wave dispersion relations on the Arakawa grids A,B,C and D for coarse resolution. The resolution parameters are rx = ry = 0.5. A grid B grid 2 ω/f ω/f 1.5 1 0 1.5 1 0 1 0.5 0.5 sk 1 0 1 0.5 0.5 sl sk C grid sl D grid 2 ω/f 2 ω/f 1 0 1.5 1 0 1 0.5 0.5 s k 1 0 s l 1 0 0 1 0.5 0.5 s k 1 0 s l Figure 7.3: Inertia-gravity wave dispersion relations on the Arakawa grids A,B,C and D resolving the deformation radius. The resolution parameters are rx = ry = 2. 7.3 Rossby waves on staggered grids The linear shallow water models of the previous section all permit the geostrophic modes or Rossby waves. In the continuum, deriving the Rossby wave disper­ sion relation is usually done by expanding the flow in a small parameter (the Rossby number): u = u o + R o u1 + . . . v = v o + R o v1 + . . . � = � o + R o �1 + . . . Substituting into the beta-plane (f = fo + ηy) SWE’s and keeping only leading order terms gives: f vo = g�x �o f uo = −g�y �o H(�x uo + �y vo ) = 0 Note the time-derivative is assumed of order Ro relative to the Coriolis term. At the next order we have: �t uo − f v1 − ηyvo + g�x �1 = 0 �t vo + f u1 + ηyuo + g�y �1 = 0 �t �o + H(�x u1 + �y v1 ) = 0 Forming the vorticity equation from the order Ro momentum equations: �t ∂o + f (�x u1 + �y v1 ) + ηvo = 0 Combining with the order Ro continuity equation we then get the potential vorticity equation: f �t ∂o − �t �o + ηvo = 0 H or g g f �t ( �2 �o − �o ) + η �x �o = 0 f H f The dispersion relation is then φ= −ηL2d k 1 + Ld2 (k 2 + l2 ) The equivalent analysis can be made on each of the staggered grids used in the previous section. Here, we will simply give the form of the geostrophic balance, PV equation and resulting dispersion relation: • A grid: −g g βi � i ; ug = βj � j f �x f �y � L2d L2d ηL2d ijj ii jj PV : �t � + β � − � + βi � = 0 β ii jj �y 2 �x �x2 −ηLd rx sk ck c2l Dispersion : φ = 1 + rx2 s2k c2k + ry2 s2l c2l Geostrophy : vg = • B grid: −g g βi � j ; ug = βj � i f �x f �y � L2d L2d ηL2d i jj ii β βi � = 0 PV : �t � + β � − � + ii jj �y 2 �x �x2 −ηLd rx sk ck Dispersion : φ = 1 + rx2 s2k c2l + ry2 s2l c2k Geostrophy : vg = • C grid: g βi � f �x −g βj � f �y � L2d L2d ηL2d ijj iijj PV : �t + βi � = 0 β � + � − � β ii jj �y 2 �x �x2 −ηLd rx sk ck c2l Dispersion : φ = 2 2 ck cl + rx2 s2k + ry2 s2l Geostrophy : vg ij = • D grid: ; ug ij = −g g βj � ij βi � ij ; ug ij = f �y f �x � L2d L2d ηL2d ijj PV : �t β � + β � − � + βi � = 0 ii jj �y 2 �x �x2 −ηLd rx sk ck c2l Dispersion : φ = 1 + rx2 s2k + ry2 s2l Geostrophy : vg ij = The denominator of each dispersion relation reflects the inertia-gravity wave dispersion relation. Indeed, it is the deficiencies of the denominator that dominates the behavior of the discrete Rossby wave dispersion. The disper­ sion relations are plotted in Figs. 7.4 and 7.5. 7.4 Higher moments in discrete SW models We derived the energy and potential vorticity equations for the SWEs in section 7.1. These conservation equations and properties do not necessary carry through to the discrete models though with some effort we can generally find a discretization that does exhibit some form of higher order conservation. The natural discretization of the continuity equation on any grid is �x�y�t � + βi (�yU ) + βj (�xV ) = 0 where U = hu and V = hv are the thickness wieghted flow. On a C grid, the natural choice for U and V is i U = hu j V = hv Suppose that we write the momentum equations as j 1 i βi (g� + K) + qV �x i 1 j βj (g� + K) − qU �t v = − �y �t u = − we can then see that the contributions to the total energy equation are �x�y g��t � = −g�βi (�yU ) − g�βj (�xV ) �x�y K�t h = −Kβi (�yU ) − Kβj (�xV ) �x�y U �t u = −�yUβi (g� + K) + �x�yUqV i �y�x V �t v = −�xV βj (g� + K) − �x�yV qU j j i Summing over all points in the domain, we can see due to the symmetry or skew symmetry properties of each term mean that all terms cancel giving � ij �x�y (g��t � + K�t h) + � ij �x�y U�t u + � ij �y�x V �t v = 0 B grid 2 2 0 0 ω/f ω/f A grid −2 −1 −2 −1 1 0.5 0 s k 1 0 1 0.5 0 sl s k C grid sl D grid 2 2 0 ω/f ω/f 1 0 −2 −1 0 −2 −1 1 0.5 0 sk 1 0 1 0.5 0 sl sk 1 0 sl Figure 7.4: Rossby wave dispersion relations on the Arakawa grids A,B,C and D resolving the deformation radius. The resolution parameters are rx = ry = 5. A grid B grid 0.5 ω/f ω/f 0.5 0 −0.5 −1 0 −0.5 −1 1 0.5 0 sk 1 0 1 0.5 0 s sk l C grid s l D grid 0.5 ω/f 0.5 ω/f 1 0 0 −0.5 −1 1 0.5 0 sk 1 0 s l 0 −0.5 −1 1 0.5 0 sk 1 0 s l Figure 7.5: Inertia-gravity wave dispersion relations on the Arakawa grids A,B,C and D for coarse resolution. The resolution parameters are rx = ry = 0.2. i j which holds if K = 12 (u2 + v 2 ). Another way to show this is to consider the total energy equation at the � points which involves averaging the contributions of the kinetic energy equation. For instance i g�βi (�yU ) + �yUβi (g�) = βi (g� i �yU ) which is difference of a flux and so conservative. The above discretization of the shallow water equations is the energy conserving scheme proposed by Sadourny (1975). An alternative form that also conserves energy is i 1 βi (g� + K) = 0 �x j i 1 i ij �t v + q h u + βj (g� + K) = 0 �y �t u − j q ij h v j + but does not seem as natural. Sadourny also proposed a potential enstrophy conserving scheme: 1 βi (g� + K) = 0 �x ji 1 i �t v + q i h u + βj (g� + K) = 0 �y j ij �t u − q j h v + The corresponding vorticity equation is i ji j ij �t (βi (�yv ) − βj (�xu)) + βi (�yq i h u ) + βj (�xq j h v ) = 0 which can be written ij ji i j ij �x�y�t (qh ) + βi (�yq i h u ) + βj (�xq j h v ) = 0 Since ij i ji j ij �x�y�t h + βi (�yh u ) + βj (�xh v ) = 0 then ji ij j j i i ij �x�y�t (q 2 h ) + βi (�yq˜2 h u ) + βj (�xq˜2 h v ) = 0 Image from Sandourny, 1972 (JAS), removed due to copyright. Figure 7.6: The evolution of potential enstrophy (Z) and energy (e) in (left) the energy conserving model and (right) the potential enstrophy conserving model. where the tilde indicates a geometric average rather than an arithmetic av­ erage. The potential enstrophy scheme does not conserve energy in general but happens to conserve energy in the special limit of pure rotational flow. Sadourny (1975) discusses the relative merits of conserving energy and poten­ tial enstrophy. Arakawa and Hsu (1990) describe a much more complicated scheme that conserves energy and dissipates enstrophy. Energy conservation is argued to ensure numerical stability since the flow must remain bounded. Potential enstrophy conservation does not guarantee a bounded solution since it does not prevent growth of the divergent part of the flow. How­ ever, Sadourny found the potential enstrophy conserving scheme to be more stable (see Fig. 7.6); while the energy conserving model does not conserve potential enstrophy in practice, the potential enstrophy conserving model does seem to conserve energy (in practice though not formally).