Atmospheric and Oceanic Modeling - Barotropic Models

Lecture 12 Barotropic equations 8.1 Barotropic equations The barotropic equations describe “balanced” motion in a two-dimensional incrompessible homogeneous fluid. A quick (and dirty) derivation involves imposing a rigid-lid on the shallow water equations; the vorticity and conti­ nuity equation are then: ρt (Hq) + ρx (Hqu) + ρy (Hqv) = 0 ρx (Hu) + ρy (Hv) = 0. (8.1) (8.2) Hq = f + ψ is the absolute vorticity. If we further assume a flat bottom then we can re-write the equations as ρt ψ + J(τ, f + ψ) ψ v u = = = = 0 ∀2 τ ρx τ −ρy τ (8.3) (8.4) (8.5) (8.6) where ψ is the relative vorticty and τ is a stream-function. J(p, q) is the Jacobian operator J(p, q) = ρx pρy q − ρy pρx q which may also be written in flux form as J(p, q) = ρx(pρy q) − ρy (pρxq) or as J(p, q) = ρy (qρx p) − ρx (qρy p) Because of these flux forms, if either p or q is either �zero or constant on � the boundary of a domain then the domain integral of J(p, q) dx dy = 0. Further, 1 < pJ(p, q) > = < J( p2 , q) >= 0 2 1 < qJ(p, q) > = < J(p, q 2 ) >= 0 2 These properties imply that domain integrated kinetic energy, |∀τ|2 , and enstrophy, |f + ψ|2 , are conserved. We’ll consider three approaches to solving the barotropic equations nu­ merically. For pragmatic purposes we must include friction and viscosity; hte most typical form of the barotropic equation is ρt ψ + J(τ, ψ) + ωρx τ = F − ζψ + ∂∀2 ψ (8.7) where F is a general forcing and other terms on the R.H.S. are arbitrary forms of dissipation; ζ is a (bottom) drag coefficient and ∂ is a lateral viscosity. The model is forced by a wind-stress curl, F = ∀ � � which is prescribed. The model domain may be either a doubly periodic domain, a channel or a bounded box with no normal flow on the sides which means that the stream function is constant on the boundary: τ = 0 on boundaries (8.8) and no-slip at the side walls: ∀τ · n̂ = 0 on boundaries 8.2 (8.9) Finite difference barotropic model We’ll describe a finite difference approximation to the above equations for the purposes of comparison. We’ll use the leap-frog method for advection (the Jacobian term) and a forward method for the dissipation terms: � ψ n+1 = ψ n−1 + 2�t F − J n − B n + D n−1 � (8.10) where F is the forcing term, J n and B n are the Jacobian and beta terms evaluated at the mid-point, n, and D n−1 is the dissipation terms evaluated at the past point, n − 1. The dissipation terms are discretized in space using second order finite differences: D n = −ζψ n + 3 3 βM βM n β ψ + βjj ψ n ii �x2 �y 2 (8.11) The beta term is a centered second order difference: Bn = 1 i βi τ n �x (8.12) The vorticity is discretized: ψ n = ωyij + 1 1 βii ψ n + βjj ψ n 2 2 �x �y (8.13) The treatment of the Jacobian term will be desribed in the next section. The complete algorithm then is as follows: • from ψ n , compute τ n by inverting ∀2 τ n = ψ n , • compute the Jacobian J n = J(τ n , ψ n ) • compute the dissipation D n−1 • compute ψ n+1 = ψ n−1 + 2�t(F − J n − B n + D n−1 ) • increment n and repeat cycle and at some interval, say every 50 steps, after computing ψ n+1 we will filter out the leap-frog computational mode by setting ψ n = (ψ n+1 + ψ n−1 )/2. 8.3 The Arakawa Jacobian Fig. 8.1 shows a grid with labels for the stencil of points involved in the simplest finite difference Jacobians proposed by Arakawa, 1960. Using these labels, the Jacobians evaluated at point 5 are defined �x�yJ++ = (τ8 − τ1 )(ψ7 − ψ4 ) − (τ6 − τ3 )(ψ7 − ψ0 ) �x�yJ+× = τ8 (ψ7 − ψ4 ) − τ1 (ψ0 − ψ2 ) − τ6 (ψ7 − ψ0 ) + τ3 (ψ4 − ψ2 ) �x�yJ×+ = ψ6 (τ7 − τ0 ) − ψ3 (τ4 − τ2 ) − ψ8 (τ7 − τ4 ) + ψ1 (τ0 − τ2 ) 13 12 11 0 6 7 10 1 5 8 9 2 3 4 Figure 8.1: The stencil of points for evaluating the Arakawa Jacobians; evalu­ ation at point 5 involves contributions from points 0 through 7 and evaluation at point 7 involves contributions from points 5 through 13. The blue cross indicates the stencil of points of the “+” terms and the red indicates the stencil of points of the “×” terms, both at point 5. The continuum energy equation is found by multiplying the vorticity equation by τ and noting that < τJ(τ, ψ) >= 0. Consider the contributions to the energy equation arising from τ5 J×+5 and τ7 J×+7 : τ5 J×+5 = τ5 ψ6 τ7 − τ5 ψ6 τ0 − τ5 ψ3 τ4 + τ5 ψ3 τ2 −τ5 ψ8 τ7 + τ5 ψ8 τ4 + τ5 ψ1 τ0 − τ5 ψ1 τ2 τ7 J×+7 = τ7 ψ12 τ11 − τ7 ψ12 τ13 − τ7 ψ8 τ9 + τ7 ψ8 τ5 −τ7 ψ10 τ11 + τ7 ψ10 τ9 + τ7 ψ6 τ13 − τ7 ψ6 τ5 We see that the underlined terms appear in both contributions with the opposite sign so that they cancel; summing over all points, all terms cancel in a similar fashion. It follows that J×+ globally conserves energy and it follows by symmetry that J+× globally conserves enstrophy. However, similar consideration of the contributions from ψJ×+ shows that there is no such cancellation: ψ5 J×+5 = ψ5 ψ6 τ7 − ψ5 ψ6 τ0 − ψ5 ψ3 τ4 + ψ5 ψ3 τ2 ψ7 J×+7 −ψ5 ψ8 τ7 + ψ5 ψ8 τ4 + ψ5 ψ1 τ0 − ψ5 ψ1 τ2 = ψ7 ψ12 τ11 − ψ7 ψ12 τ13 − ψ7 ψ8 τ9 + τ7 ψ8 τ5 −ψ7 ψ10 τ11 + ψ7 ψ10 τ9 + ψ7 ψ6 τ13 − ψ7 ψ6 τ5 There is only one term with the quantity ψ5 ψ6 τ7 and so nothing can cancel it. However, the J++ conserves neither but it happens that The most obvious discretization of the Jacobian: J++ (τ, ψ) = � 1 � i j j i βi τ βj ψ − β j τ βi ψ �x�y (8.14) corresponds to the continuum form, J(τ, ψ) = ρx τρy ψ − ρy τρx ψ. The flux form Jacobian � j i 1 j i βi (τβj ψ ) − βj (τβi ψ ) J+× (τ, ψ) = �x�y (8.15) corresponds to the continuum form, J(τ, ψ) = ρx (τρy ψ) − ρy (τρx ψ) and the other flux form Jacobian � j i 1 i j βj (ψβi τ ) − βi (ψβj τ ) J×+ (τ, ψ) = �x�y (8.16) corresponds to the continuum form, J(τ, ψ) = ρy (ψρx τ) − ρx (ψρy τ). In the fashion shown above, Arakawa, 1960, showed that • J+× conserves the integral of enstrophy, ψ 2 , • 1 (J++ 2 + J×+ ) conserves the integral of enstrophy, ψ 2 , • J×+ conserves the integral of energy, ∀τ · ∀τ, • 1 (J++ 2 + J+× ) conserves the integral of energy, ∀τ · ∀τ, so that the famous combination of 1 J n = (J++ + J×+ + J+× ) 3 conserves both the enstrophy and energy and therefore also conserves mean wave number; this prevents nonlinear instability occurring. 8.4 Solving elliptic boundary value problems In the preceding finite difference model, we had to solve a problem of the form ∀2 τ = ψ given some boundary conditions on τ (e.g. τ = 0). One approach is to note that the solution to the above elliptic problem is the steady state solution to the following time-dependent ρ t τ = ∀2 τ − ψ This is simply a diffusion problem with a source term and can be solved use simple explicit method with an appropriately chosen time-step. Consider the FTCS discretization τ n+1 = τ n + � �t � n n n n τ + τ + τ + τ − 4τ i,j − �tψ i+1,j i,j−1 i,j+1 �x2 i−1,j This system can be integrated forward in a stable fashion so long as �t � 1 �x2 . If we chose the largest time-step then the system simplifies to 4 τ n+1 = � 1� n �x2 n n τi−1,j + τin+1,j + τi,j−1 + τi,j+1 − ψ 4 4 which is known as Jacobi’s method. The method converges but converges slowly. A slightly more efficient iterative method is the Gauss-Seidel method which uses the latest estimates on the RHS as the computations procede. A commonly used GS method is known as the red-black Guass-Seidel method where the points on the grid are labelled in a checker-board fashion (red­ black). The Jacobi method is applied to the red and black points seperately, using the latest values in each consecutive copmutation. A more efficient iterative method is successive over-relaxation (SOR) which is described in some detail in Numerical Recipes, chapter 19. For the recast problem n n ai,j τin−1,j + bi,j τin+1,j + ci,j τi,j−1 + di,j τi,j+1 + ei,j τi,j = ψ the SOR algorithm is � 1 � n n ψ − ai,j τin−1,j − bi,j τin+1,j − ci,j τi,j−1 − di,j τi,j+1 ei,j � n = φτi,j + (1 − φ)τi,j τ� = n+1 τi,j where φ0 = 1 1 φ 1 = (1 − �2J acobi )−1 2 1 φ n = (1 − φ n−1 �2J acobi )−1 n = 2, 3, . . . 4 and �J acobi is a function of a, b, c, d and e. The original elliptic problem can be solved directly by Fourier analysis. Suppose we have Dirichlet boundary conditions (τ = 0 on the boundaries) then a natural basis function is sine waves. The solution is represented by the Fourier series τi,j = N −1 M −1 � � ak,l sin νik νjl sin N M N −1 M −1 � � bk,l sin νik νjl sin N M k=1 l=1 and ψi,j = k=1 l=1 This allows us to solve the elliptic problem by first computing bk,l = −1 M −1 � 2 2 N� νik νjl ψi,j sin sin M N M i=1 j=1 N and then for each Fourier mode � νjl ν(i + 1)k νjl ν(i − 1)k sin + sin sin M N M N νik ν(j − 1)l νik ν(j + 1)l + sin sin + sin sin M N M N νjl νjl νjl νjl = bk,l sin sin −4 sin sin M M M M ak,l sin or � νjl νik νk sin − 2) sin M N N νjl νl νik +(2 cos − 2) sin sin M M N νjl νjl νjl νjl = bk,l sin sin sin −4 sin M M M M ak,l (2 cos or simply ak,l = 8.5 �k N 2 cos bk,l + 2 cos �k −4 N Finite element barotropic model We’ll now solve the barotropic equations using the Galerkin approximation. The basis function must now be two-dimensional. Since we are in a rect­ angular domain, we can use a square based function known as the pagoda function: |x − xi | |y − yi | )(1 − )) �ij (x, y) = max(0, (1 − �x �y which gives bi-linear interpolation between values on a quadrilateral grid. The result of the Galerkin approximation is: Ax Ay ρt ψ + J + Ay ω β3 β3 i βi τ = −ζAx Ay ψ − M2 Ay βii ψ − M2 Ax βjj ψ (8.17) �x �x �y The operator Ax Ay is a two-dimensional operator and thus potentially ex­ pensive to invert but since it is separable into two tri-diagonal operators we can invert it relatively efficiently. The Jacobian term is the complicated term to calculate. The Galerkin approximation is: J = � � � � �ij ρx τρy ψ dx dy − −� −� � � −� � � −� Substituting in the series expansions τ = � τk,l �k,l (x, y) k,l ψ = � m,n ψm,n �m,n (x, y) �ij ρy τρx ψ dx dy into the first integral in J gives � � dx −� � � dy �ij −� � τk,l ρx �k,l (x, y) � ψm,n ρy �m,n (x, y) m,n k,l Since the pagoda functions only interact with the eight neighboring elements the summations reduce to sums over k = i − 1, i, i + 1 ; l = j − 1, j, j + 1 ; m = i − 1, i, i + 1 ; n = j − 1, j, j + 1 we can anticipate up to 9 × 9 = 81 coefficients arising from the double sum. Fortunately, many of these are zero since some of the pairs of basis functions, (k, l) and (m, n) don’t overlap (e.g. (k, l) = (i − 1, j − 1) and (m, n) = (i + 1, j + 1)). Others can be deduced to be zero by symmetry, for � � example, dx dy�i,j ρx �i,j+1 ρy �i,j+1 = 0. However, rather than compute all these coefficients we can make an inspired guess based on the structures of the other terms: i j j �x�yJ = Ay βi τ Ax βj ψ − Ax βj τ Ay βi ψ i If we use the leap-frog scheme for the Jacobian and beta terms and forward method for the dissipation terms then our model algorithm is: • compute D n−1 = −Ax Ay ζψ n−1 − • compute B n = −Ax βi τ n • compute J n = 1 �x�y � 3 �M Ay βii ψ n−1 �x2 − 3 �M Ax βjj ψ n−1 �y 2 i i j j Ay βi τ Ax βj ψ − Ax βj τ Ay βi ψ • invert Ax Ay ρt ψ = F − J n − B n + D n−1 i � • step forward for ψ n+1 = ψ n−1 + 2�tρt ψ • increment n and repeat cycle An alternative approach is to invert the equation so that each term re­ quires a matrix solve. Multiply through by the inverse of Ax Ay : ρt q = −J � − Ax−1 1 β3 β3 i βi τ − ζq − M2 Ax−1 βii q − M2 Ay −1 βjj q �x �x �y Now, solve Ax ρ x τ � = Ay ρy τ � = Ax ρ x q � = Ay ρ y q � = Ay ρxx q � = Ay ρyy q � = 1 i βi τ �x 1 j βj τ �y 1 βi q i �x 1 βj q j �y 1 βii q �x2 1 βjj q �y 2 (8.18) (8.19) (8.20) (8.21) (8.22) (8.23) And then the Jacobian can be appriximated as J � = ρ x τ � ρy q � − ρ y τ � ρx q � and similarly the beta term is B � = ρx τ � and dissipation term 3 D � = βM (ρxx q � + ρyy q � ) This method is the method of compact differencing and is gives fourth order accuracy. The prognostic equation then becomes � q n+1 = q n−1 + 2�t F − J n� − B n� + D n−1 8.5.1 � � Spectral and pseudo-spectral barotropic model For simplicity, we will limit ourselves to a doubly periodic domain or di­ mensions 2ν by 2ν. This reduces the number of coefficients in the following model. Using fourier basis functions and the Galerkin method we can obtain the differential equation for the spectral coefficients (k 2 + l2 )ρt akl + � � ωik + ζ(k 2 + l2 ) + ∂(k 2 + l2 )2 akl � � 1 = �k (x)�l (y)(F − J ) dA Ikkll where Ikkll = Ikk Ill and Ink = � �n (x)�k (x) dx noting that Ink = 0 � n ≤= k because �k (x) = eikx form an orthogonal set of basis functions. In the spectral method, the fourier transform of J leads to a convolution and is relatively expensive to evaluate. The transform method is the more efficient approach in which J is evaluated in physical space and the transformed to wave-number space: τx� (xi , yj ) = �� k τy� (xi , yj ) = = = J(xi , yi ) = ilakl e−ikxi e−ilyj l �� k ψy� (xi , yj ) l �� k ψx� (xi , yj ) ikakl e−ikxi e−ilyj −ik(k 2 + l2 )akl e−ikxi e−ilyj l �� −il(k 2 + l2 )akl e−ikxi e−ilyj k l � � τ x ψy − τy� ψx� 1 �� J(xi , yj )eikxi eilyj NM i j 1 �� = F (xi , yj )eikxi eilyj NM i j Jkl = Fkl so that we can now integrate � � (k 2 + l2 )ρt akl + ωik + ζ(k 2 + l2 ) + ∂(k 2 + l2 )2 akl = Fkl − Jkl Note that the left hand side can be integrated in time analytically but that some explicit and stable time-stepping method must be used for the right hand side. 8.5.2 Aliasing and nonlinear aliasing Discrete representations of functions always allow aliasing of unresolved wave numbers onto resolved numbers. Fig. 8.2 shows an example using the regular finite difference grid xj = 2νj/N with N = 8. Evaluating the function cos(5x) at xj (black) only on the finite difference grid gives the eight values indicated by the circles. The highest mode that can be resolved on this grid is cos(3x), cos(7x) & cos(11x) 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 0.1 0.2 0.3 0.4 0.5 x/(2*π) 0.6 0.7 0.8 0.9 1 Figure 8.2: Aliasing of high wave number to a low wave number on a finite resolution grid: the grid has 8 nodes and on that grid, cos(5x) (solid black)is indistinguishable (circles) from cos(3x) (dashed blue). The highest mode that could be represented is cos 4x. The mode cos(11x) (red) is also aliased to cos(3x). cos( N2 x) = cos(4x) and all higher mode number waves are mis-represented on this grid. In this case, cos(5x) is aliased to cos(3x); evaluating the function cos(3x) on the finite difference grid gives the exact same values as from cos(5x). More generally, the relation 2�n 2�n eikxj = eikj�x = ei2�n eikj�x = ei(k+ �x )j�x = ei(k+ �x )xj � n = 0, ±1, ±2, . . . means that a resolved mode k can result in aliasing an infinite set of modes k + 2νn/�x. In the above example, with �x = 2ν/8, n = −1 corresponds to cos(5x) and n = 1 corresponds to cos(11x). Non-linear instability can occur when the result of nonlinear interactions between two waves result in higher mode number waves that can not be represented on the grid. The aliasing of these unresolved waves implies a false source of energy at lower wave numbers. Thus, if nonlinear interactions cos(3x), cos(4x) & cos(3x)cos(4x) 1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 0.1 0.2 0.3 0.4 0.5 x/(2*π) 0.6 0.7 0.8 0.9 1 Figure 8.3: Nonlinear instability occurs when the product of two resolved waves leads to a high mode number component that can not be resolved by the grid and is aliased to lower resolved waves. Here, we show the same 8­ point grid with the functions cos(3x) (blue circles) and cos(5x) (red circles) and their product cos 3x cos 5x = 21 cos(8x) + 21 cos(2x) (black squares) which has a cos(8x) component that aliases to cos(2x). produces a wave number k1 + k2 > N/2, then it will falsely appear as k3 = N − (k1 + k2 ). We illustrate this in Fig. 8.3 in which the product of cos(3x) and cos(4x) on an 8-point grid results in a higher mode signal with