Atmospheric and Oceanic Modeling - Quasi-Geostrophic Models 2

1 Introduction We describe here the quasi-geostrophic ocean model that was developed by John Marshall, George Nurser and Roger Brugge - see Marshall, J., A.J.G. Nurser and R. Brugge: (1988) On the time-averaged °ow of quasi-geostrophic wind-driven gyres, J. Geophys. Res. (Oceans), 93, 15427-15436. It integrates a prognostic (potential vorticity) equation for each layer in a closed basin with (linearized) bottom topography and can be driven by both wind and buoyancy forces. The model can be integrated in complex geometries; basin and channel versions are available. Both wind and buoyancy driven °ows can be studied, but within the con¯nes of quasi-geostrophic dynamics. 2 The model The multi-layer quasi-geostrophic (Q-G) basin model integrates the potential vorticity equation !"! + $(%! & "! ) = '! (1) !# in each layer (, where the Jacobian $()& *) is given by $()& *) = !) !* !) !* ¡ !+ !, !, !+ (2) It functions as follows: (I) Given the Q-G potential vorticity " and streamfunction % in each layer and known forcing functions '! , the tendency ¯eld !".!# is calculated and used to update the " ¯eld. (II) This revised " ¯eld is inverted with appropriate boundary conditions to give the new streamfunction ¯eld (III) These new " and % ¯elds are then used to calculate the tendency in (I) and so on. The equations are ¯nite-di®erenced on a regular grid using standard second order ¯nite di®erences and are stepped forward in time using a leap-frog scheme. Boundary conditions upon the baroclinic modes required for the elliptic problem are deduced from the constraint that the total quantity of °uid in each layer is known from the continuity equation. A high order r6 % (enstrophy-destroying) friction is employed demanding two extra boundary conditions. We have chosen to apply the boundary conditions r2 % = r4 % = 0. Other boundary conditions (for example no-slip) can easily be introduced. Before going on to describe the numerical implementation we discuss the dynamical setup, the potential vorticity equation and its non-dimensionalisation. 1 Figure 2.1: The ocean basin and its topography. (See text for meaning of symbols.) 2.1 The physics of the model The equations are solved in an /" £ /# rectangular basin (typically /" & /# = 3000 km) with bottom topography 0(+& ,), see Fig. 2.1. The ocean is taken to have total depth 1 (typically 5000 m) with 2 layers of thicknesses 11 & - - - & 1! & - - - & 1$ , (see Fig. 2.2) and densities 31 & - - - & 3! & - - - & 3$ . We denote the density jumps 3!+1%! = 3!+1 ¡ 3! between the layers. The 4!+1%! , the displacements from equilibrium of the interfaces between the (( + 1)th and (th layer, and also the bottom topography 0(+& ,) are required [for consistency with the quasi-geostrophic scaling] to be small relative to the layer thicknesses. Note that the double su±x e.g. 4!+1%! denotes a quantity evaluated at the interfaces between layers (( + 1) and ( | the `thermodynamic levels' in the equivalent level formulation. The Q-G potential vorticity equation may be derived from: 1) the vorticity equation for each layer !6 0 (r2 % + *,) = 50 0# !7 (3) which in the (th layer takes the ¯nite di®erence form 50 0 (r2 %! + *,) = (6!%!¡1 ¡ 6!+1%! ) 0# 1! (4) (here 6!+1%! is the vertical velocity at the ( + 1/(th interface, r2 %! is the relative vorticity, *, is the planetary vorticity and 50 is the value of the Coriolis parameter at the southern boundary), and 2) the equations for the movement of the interfaces between the layers 6!+1%! = 8!+1%! + 2 0 4!+1%! 0# (5) Figure 2.2: Schematic representation of the vertical model structure. (See text for meaning of symbols.) where 4!+1%! is the deviation of the interface between layers (( + 1) and (. Thus 04.0# represents the vertical velocity of the interface and the cross-interface velocities 8!+1%! represent the transformation of °uid from one density to another, and may be considered as arising from heating terms in the temperature perturbation equation to which, applied at the thermodynamic levels, Eqn. (5) is equivalent. The thermal wind relation gives the interface displacements in terms of the streamfunctions in the layers above and below thus 4!+1%! = 50 30 (%!+1 ¡ %! ) 9¢3!+1%! (6) where 9 is the gravitational acceleration and 30 is the mean density of water. Substitution of Eqn. (6) into Eqn.(5) gives (in fact proportional to) the `buoyancy equation' 0 9¢3!+1%! (6!+1%! ¡ 8!+1%! ) = 0 (7) (%! ¡ %!+1 ) + 0# 30 50 or, expanding the substantial derivative ! 9¢3!+1%! (6!+1%! ¡ 8!+1%! ) = 0(%! ¡ %!+1 ) + $(%! & %! ¡ %!+1 ) + !# 30 50 (8) Elimination of vertical velocities in Eqn. (4) with Eqn.(8) gives the Q-G potential vorticity equation !"! (9) + $(%! & "! ) = '! !# 3 where "1 "! &02 '0 (% ¡ %2 ) () n 1 ¢'212 1 & ' &02 '0 0 0 r2 %! + *, ¡ ()! ¢' (%! ¡ %!+1 ) + ()! ¢' (%! !+1"! !"!¡1 2 &0 '0 r2 %$ + *, ¡ ()# ¢' (%$ ¡ %$¡1 ) + &)0#* #"# ¡1 = r2 %1 + *, ¡ = "$ = 9 > = o > ¡ %!¡1 ) are the quasi-geostrophic potential vorticities in each layer and &0 (60 ¡ 821 ) ¡ :r6 %1 )1 &0 (8!%!¡1 ¡ 8!+1%! ) ¡ :r6 %! )! &0 8$%$¡1 ¡ :r6 %$ ¡ ;r2 %$ )# '1 = '! = '$ = 9 > = > ; > > ; (10) (11) are the potential vorticity forcing functions. Note 1) 60 is the Ekman suction velocity imposed at the top of the surface layer and can be related to the imposed wind stress (< ) curl thus: 60 = =->?@A(< .30 5)- (12) 2) The vertical velocity at the bottom of the ocean is given by 6+,--,. = 6/ + 0 0 0# (13) where 0 = 0(+& ,) is the (assumed small) bottom topography and 6/ is the rate of ejection of °uid from the bottom Ekman layer. This has been parameterised in Eqn. (11) by assuming 6/ = ;1$ 2 r %$ 50 (14) 3) It may be shown that the (equivalent) perturbation temperature B!+1%! at the interface is given by 50 (%! ¡ %!+1 ) B!+1%! = (15) 9) 0-5(1! + 1!+1 ) where ) = (¡!3.!B ).30 is the coe±cient of thermal expansion. The cross-interface velocity between the (th and (( + 1)th layers may then be related to the total heat input/unit area H0 experienced between the middle of the (( + 1)th and (th layers by 8!+1%! = ) H0 C1 ¢3!+1%! (16) Values of ) = 2 £ 10¡4 K¡1 , C1 = 4000 Jkg¡1 K¡1 appropriate to near surface waters gives (in ms¡1 ) (17) 8!+1%! = 5 £ 10¡8 0 H0 .¢3!+1%! if ¢3!+1%! is speci¯ed in kgm¡3 and H0 in Wm¡2 . 4 2.2 The non-dimensional equations We choose to scale + and , by /# , and 1! and 0 by 1, the total depth of the ocean. We write the depth of the layer as a fraction of the total depth D! = 1! .1, E = 0.1, and the density jumps ¢3!+1%! .30 = ¢F!+1%! £ 10¡3 , where ¢F!+1%! is the di®erence in density expressed in F units. If time is scaled by (*/# )¡1 then for /# = 3000 km and * = 2 £ 10¡11 , this timescale (*/# )¡1 » 10 hours. We scale % by G2 /# where G2 is a characteristic velocity, and scale q by */# . Then Eqn. (9) becomes 1 !"! + $(%! & "! ) = '0! (18) H !# where H = G2 .(*/2# ) is a Rossby number for the vorticity equation, 9 n o 342 1 ¡62 > "1 = Hr2 %1 + , ¡ 51 4$2 6¢7 > > 21 % n = o 2 34$ 6! ¡6!¡1 6! ¡6!+1 2 "! = Hr %! + , ¡ 5! 42 ¢7!"!¡1 + ¢7!+1"! (19) % n > o > 2 34 > ¡6 6 # # ¡1 "$ = Hr2 %$ + , ¡ 5# 4$2 ¢7 + 94&0$85# ; #"# ¡1 % the `pseudo Rossby radius' /: = q () 103 &02 (of magnitude 70 km for 1 = 5000 m and 50 = 10¡4 s¡1 ). The forcing terms (scaled by *G2 ) are '01 = '0! = '0$ = 1 ;& (60 ¡ 821 ) ¡ 94< 5 r6 %1 51 ;' $ ! ;& < 6 (8 ¡ 8 !%!¡1 !+1%! ) ¡ 945 r %! 5! ;' $ $ ;& 8$%$¡1 ¡ 94< 5 r6 %$ ¡ 94= $ r2 %$ 5# ;' $ where the Sverdrup velocity is given by 9 > = > ; G> = 50 60 .(*1)- (20) (21) After scaling, the thermodynamic equation becomes G> /2: ! ¢F!+1%! (6!+1%! ¡ 8!+1%! ) = 0(%! ¡ %!+1 ) + H$(%! & %! ¡ %!+1 ) + !# G2 /2 (22) We have scaled the vertical velocities 8 by 60 , a typical value of the Ekman pumping chosen to be 30 m yr¡1 (10¡6 ms¡1 ). G> is the depth-averaged value of the meridional velocity implied by the Sverdrup balance resulting from this Ekman pumping. In oceanographic studies we have taken the characteristic velocity G2 = G> ; for the values of *, 50 , 1 and 60 chosen above, G> » 1mms¡1 . 2.3 Numerical integration of the quasi-geostrophic potential vorticity equation The major steps in the numerical integration are now outlined. 5 A centred di®erence leapfrog scheme is used for the forward time-stepping of the ¯nite di®erence version of the potential vorticity equation (18). Thus ª © -+1 -¡1 (23) ¡ "?%@ = 2¢#H $?%@ + '0-¡1 "?%@ ?%@ is the Jacobian and the where the ?%@ su±x denotes the value at the i,jth gridpoint, the $?%@ -¡1 "?%@ are the "'s from the preceding timestep. is calculated using values of " and % from the current timestep. The The Jacobian $?%@ Arakawa Jacobian (Arakawa, 1966) formulation is used. This ensures that there is no net advection out of the region of interior points over which time-stepping is performed of 1) potential vorticity 2) enstrophy 12 " 2 or 3) energy, so long as % is constant along the boundary X X X i.e. $?%@ = "?%@ $?%@ = %?%@ $?%@ = 0(24) (!)*+(,+ 1,?!-A (!)*+(,+ (!)*+(,+ 1,?!-A 1,?!-A This stepping forward is only performed at the interior gridpoints. Values of % are then found everywhere by solving an elliptic equation with appropriate boundary conditions. The inversion of the Helmholtz equations does not require knowledge of " along the boundary. The new " along the boundary is, however, required in the calculation of the Jacobian at the points immediately adjacent to the boundary. It is determined by imposing the condition that the relative vorticity r2 % = 0 on the boundaries. Values of " and % from the previous timestep are used to calculate the frictional contributions (Ekman | r2 %$ , and biharmonic | r6 %! ) to the forcing terms '0! . Calculation of r2 % is straightforward: the biharmonic friction is evaluated by repeated application of the ¯nite di®erence r2 operator to the above determined r2 % using the condition that r4 % = 0 on the boundary. The mechanical forcing (appearing as an Ekman suction velocity generated by the wind stress curl) and thermal forcing (imposed as interfacial velocities) are then calculated and added to the frictional terms to give the '0! 's. 2.4 Inverting the potential vorticity Eqn. (19) may be written in matrix form, eg. " = *, + r2 % ¡ I% (25) where I is a positive de¯nite tridiagonal stretching matrix, not necessarily orthogonal, and ", , and % are column vectors comprised of the values of ", , and % on the layers. The matrix I is de¯ned by 9 I1%1 = J: .D1 (F2 ¡ F1 ) > > > > I1%2 = ¡I1%1 > > 9 > > I?%?¡1 = ¡J: .D? (F? ¡ F?¡1 ) = > > = I?%? = ¡I?%?¡1 ¡ I?%?+1 K = 2& - - - 2 ¡ 1 (26) ; I?%?+1 = ¡J: .D? (F?+1 ¡ F? ) > > > > I$%$¡1 = ¡J: .D$ (F$ ¡ F$¡1 ) > > > > I$%$ = ¡I$%$ ¡1 > > ; I?%@ = 0 jK ¡ Lj ¸ 2 6 where 103 G2 502 (27) *91 and 1 is the total model depth, G2 is the characteristic velocity and 50 is the value of the Coriolis parameter at the southern boundary. To solve for %, we require that (i) I is non-singular [i.e. EM#I 6= 0(9 I st I I¡1 = 1)] and (ii) all eigenvalues of I are real. There are then 2 (not necessarily di®erent) non-zero eigenvalues N. , with corresponding eigenvectors M. such that J: = I M. = N. M. & (28) . . : where M. = (M. 1 & M2 & - - - M$ ) (a column vector). So if a function 9 = (91 & 92 & - - - & 9$ ), can be written 9 = 9M. (29) (i.e. its values in di®erent layers are in the same proportion as those of the eigenvectors M. ), then (30) r2 9 ¡ I 9 = r2 9 ¡ N. 9We will use this to ¯nd a solution of Eqn. (25). Because I is non-singular, these M. are `complete', i.e. we can express any 9 P$ . . > % = % M = .=1 P$ . . , = ,1 = , .=1 1 M P$ > ; . . " = .=1 " M (31) So, expressing Eqn. (25) in terms of the eigenvectors M. ". = *,1. + r2 % . ¡ N. % . & O = 1& - - - & 2- (32) This results in 2 independent linear elliptic equations, which are readily solved (using a Helmholtz-equation solver) for % . . It is then easy to recover % = (%1 & %2 & - - - & %$ ) from the % . by using Eqn. (31). The method of solution is thus (i) project " on to the eigenvectors M. , and ¯nd ". and *,1. , (ii) use a Helmholtz equation solver to ¯nd % . , (iii) recover % from % . . Consider Eqn. (25); suppose we pre-multiply by P ¡1 , where 0 M11 M21 - - - M$ 1 B M1 M2 - - - M$ 2 2 B 2 P = (M1 & M2 & - - - & M$ ) = B .. .. .. @ . . . 1 2 M$ M$ - - - M$ $ 7 1 C C C& A (33) to give P ¡1 " = *,P ¡1 1 + r2 P ¡1 % ¡ P ¡1 I P P ¡1 % (34) which can be written as Q = *,¨ + r2 ª ¡ 0 ª& (35) 0 = P ¡1 I P- (36) where Q = P ¡1 ", ¨ = P ¡1 1, ª = P ¡1 % are the column vectors composed of the projections of ", 1 and % respectively on to the 2 eigenvectors M. , and 0 B B Now P ¡1 P = B @ 1 0 .. . 0 --1 --.. . 0 0 .. . 0 0 --- 1 From Eqn. (33) we can express 1 C C C and P = (M1 & M2 & - - - & M$ ) : the columns are eigenvectors. A P ¡1 where the rows are M¤. , with 0 1 ¢ ¢ ¢ M¤1 ¢ ¢ ¢ B ¢ ¢ ¢ M¤2 ¢ ¢ ¢ C B C =B C .. @ A . ¤$ ¢¢¢ ¢¢¢ M (37) M¤. -M! = D.%! - (38) These M¤. are known as the complementary eigenvectors. Now consider 0 = P ¡1 I P : I P = I(M1 & M2 & - - - & M$ ) = (N1 M1 & N2 M2 & - - - & N$ M$ )- (39) Hence 0 1 0 ¢ ¢ ¢ M¤1 ¢ ¢ ¢ B ¢ ¢ ¢ M¤2 ¢ ¢ ¢ C B B C B 1 2 $ M & N M & & N M ) = 0=B (N C 1 B .. 2 $ @ A @ . ¢ ¢ ¢ M¤$ ¢ ¢ ¢ N1 0 - - - 0 0 N2 - - - 0 .. .. .. . . . 0 0 - - - N$ 1 C C CA (40) 0 is a diagonal matrix composed of eigenvalues, and % (the values of % on each of the layers) may be found from the modal projections ª using 0 1 1 % B %2 C C % = P ª = (M1 & M2 & - - - & M$ ) B (41) @ ¢¢¢ A%$ 8 2.4.1 Boundary conditions We now consider the boundary conditions to be imposed on the side boundary values of % . when evaluating them from the " . [by inversion of the Helmholtz equation (32)]. For the barotropic mode, it is clearly acceptable to require % 1 = 0 along the side boundaries. However, it is less clear what are the appropriate boundary conditions on the baroclinic modes. By considering the horizontally-averaged mean interface displacements, conservation of mass requires "# o 9¢3 ! n "# "# ¢ !+1%! ¡ %!+1 ¡ %! 6!+1%! ¡ 8!+1%! = !# 30 50 (42) where R denotes a horizontal average of R over the basin, 6!+1%! is the vertical velocity at the interface and 8!+1%! is the cross-interface velocity resulting from buoyancy forcing. Note that in the horizontal averaging, half weight must be given to points on the boundary and quarter weight to corner points, since the condition % = 0 is applied at the gridpoints. The horizontal advection terms must sum to zero since there can be no °ow through the side boundaries. Initially, we will require 6!+1%! = 8!+1%! = 0, which represents the condition of no net vertical motion or buoyancy forcing. Since the mean interfacial displacements are propor"# "# "# tional to ¥!+1%! = %!+1 ¡ %! , it is clear that these (2 ¡ 1) ¥!+1%! (there are 2 layers) are determined if we know the mean vertical velocities and buoyancy forcings. Let 0 P$ ¤1 "# 1 !=1 M %! "# C B ³ "# ´ "# "# : C B %2 ¡ %1 1 =B S = % & ¥21 & - - - & ¥$%$¡1 C .. A @ . "# %$ ¡ %$¡1 10 0 "# 1 ¤1 %1 M "# C B ¡1 B 1 0 --- 0 C C B %2 "# C B "# C B B 1 --- 0 C = B 0 ¡1 C B %3 C = T % C C B .. B .. .. .. . @ . . . . A @ .. A "# 0 0 0 --- 1 %$ As long as T is non-singular (detT 6= 0), we can write % "# = T ¡1 S- (43) But we can express the modal projection in terms of levels thus: ª = P ¡1 %& "# so the mean modal projections ª "# ª (44) must satisfy = P ¡1 T ¡1 S = P 2 S = U9 (45) Thus the continuity equation (42) has given us the U. , the horizontally-averaged values "# of the baroclinic modes. These will generally be di®erent from the %0. , the spatial averages of the solutions derived by inverting r2 %0. ¡ N. %0. = " . (46) with %0. = 0 on the boundaries. The solution is thus found by adding a multiple > of the homogeneous solution %+. which obeys (47) r2 %+. ¡ N. %+. = 0 with %+. = 1 on the boundaries, such that %0. + >%+. "# = U. . ³ ´ "# "# Thus (i) Eqn. (46) is used to calculate %0. , (ii) > is set equal to U. ¡ %0. .%+. (the %+. ¯elds do not change and therefore need only be calculated once), and (iii) >%+. is added to the %0. ¯elds to give the new % . ¯eld. 2.4.2 Summary In the model the following procedure is adopted: (1) We de¯ne the matrix I . (2) The eigenvalues N. of I and the matrix P are found; the column vectors of P are the eigenvectors of I, de¯ned by I M. = N. M. - (48) (3) The inverse matrix of P, P ¡1 , is computed. (4) The matrix product of 0 M¤1 1 B ¡1 B T = B .. @ . 0 T and P is calculated by multiplying 0 1 ¤1 M¤1 M M11 M21 - - 2 $ B M1 M2 - - 1 --- 0 C 2 B 2 C = & P B .. C .. .. .. @ . . . A . 0 --- 1 M1$ M2$ - - - the two matrices 1 M$ 1 C M$ 2 C (49) .. C . A M$ $ (5) The inverse of this product is evaluated to give P 2 P 2 = (T P)¡1 = P ¡1 T ¡1 "# (6) All the elements of ¥ (50) are set to zero initially (7) %+. is set to 1 on the boundaries for modes 2 to 2. The interior values are then "# calculated (these remain ¯xed throughout the simulation) and %+. is found by horizontal averaging of %+. . The steps which are taken in inverting the " ¯eld to give the % ¯eld each timestep now follow: 10 "# (8) ³ ¥ is found ´for levels 2 to 2 by determining the mean interfacial displacements "# %? ¡ %?¡1 . (9) %0. is set to zero on all the boundaries. (10) After subtracting *, and the topography from ", " . is found by projecting " on to the modes, eg. $ X . (PB¡1 ). "B (51) " = B=1 where = is the level number, and O is the mode number. (11) The Helmholtz equation solver is used to invert the " . to give the %0. for all the "# internal gridpoints %0. is then calculated by horizontal averaging. "# "# (12) For the ¯rst (barotropic) mode,¥ is equated to %01 (i.e. the ¯rst element in vector S) while for the other (baroclinic) modes U. is found by evaluating the product of P 2 and S. The coe±cient > is then found using "# "# > = (U. ¡ %0. ).%+. - (52) (13) % . is then evaluated using %0. + >%+. . (14) The streamfunction in each layer of the model is then determined from the modes using Eqn. (31). 11