Lecture 23 Notes: Vertical Coordinates

Lecture 23 Vertical coordinates We have exclusively used height as the vertical coordinate but there are alter­ native vertical coordinates in use in ocean models, most notably the terrainfollowing coordinate models and isopycnal models. We’ll first transform the equations of motion to a general vertical coordinate, r, before considering the most common specific vertical coordinates. A good reference for the trans­ formation is the book by Haltiner and Williams, 1980, section 1-9. The book by Haidvogel and Beckmann, 1999, describes models using all three vertical coordinates. 10.1 General vertical coordinate Consider a general vertical coordinate, r, which is assumed to be a monotonic function of height, z. Transforming the equations of motion into the coordi­ nates (x, y, r, t) makes z = z(x, y, r, t) a dependent variable. Any dependent variable, A = A(x, y, z, t), that can be described in the original coordinate can be described as a function of the new coordinates: A = A(x, y, z(x, y, r, t), t) The vertical partial derivative in the new coordinate is ρz ρA ρA = ρr ρz ρr or equivalently ρA ρr ρA = ρz ρr ρz A partial derivative with respect to any other ordinate, s, where s is one of x, y or t, is obtained using the chain rule: � � � ρA �� ρA �� ρA ρz �� � = � + � ρs �r ρs �z ρz ρs �r or, substituting in for ρz A from above, � � � ρA �� ρA �� ρA ρr ρA �� � = � + � ρs �r ρs �z ρr ρz ρr �r Here, the vertical bar with suffix indicates the partial derivative holding constant the suffix. We can now write down the horizontal gradient of a scalar A in height coordinates expressed in terms of derivatives in r coordinates: �z A = � r A − ρA ρr �r z ρr ρz Similarly, the horizontal divergence of a horizontal vector, κv can be expressed: �z · κv = �r · κv − ρκv ρr · �r z ρr ρz The vertical velocity in height coordinates, w, can be expressed: w = Dt z = = � � � ρz �� ρz �� ρz �� ρz � + � Dt x + � Dt y + Dt r ρt �r ρx �r ρy �r ρr � ρz �� ρz � +κ v · �r z + r˙ � ρt r ρr Finally, the total derivative of A in z coordinates is the total derivative in r coordinates: Dt A = = = � ρA �� ρA � +κ v · �z A + w � ρt z ρz � � � ρr ρA ρA �� ρz �� v · �r A + w − v · �r z � +κ � −κ ρt �r ρt �r ρz ρr � ρA ρA �� v · �r A + ṙ � +κ � ρt r ρr Now we can use the above relations to transform the terms in the equa­ tions of motion. Hydrostatic balance becomes ρr ρp = −g� ρz ρr or more naturally ρp ρgz = −� ρr ρr The horizontal pressure gradient in the momentum equations becomes ρp �r z ρz = �r p + ��r gz �z p = � r p − where gz is the geopotential. The incompressible continuity (or non-divergence) equation is: �z · κv + ρz w = �r · κv + ρr ρr ρκv · �r z = 0 ρr w − ρz ρz ρr The term ρr w can be found from the definition of w = Dt z above: � ρw ρ = ρr ρr � ρz �� ρz v · �r z + ṙ � +κ � ρt r ρr � ρṙ ρz ρκv ρ �� ρ · �r z + v · �r (ρr z) + ṙ (ρr z) + = � (ρr z) + κ ρt �r ρr ρr ρr ρr ρκv ρṙ ρz · �r z + = Dt (ρr z) + ρr ρr ρr so that continuity in r coordinates is ρz ρz ρṙ �r · κv + Dt + ρr ρr ρr � =0 This can be written equivalently as Dt ln or ρt ρṙ ρz + �r · κv + =0 ρr ρr ρ ρz ρz ρz ṙ = 0 + �r · κv + ρr ρr ρr ρr Combining all the above relations, the inviscid and adiabatic HPEs can then be written in r coordinates as: 1 1 Dtκv + f k̂ � κv + �r p + ��r gz = 0 �¯ �¯ ρ ρp + � gz = 0 ρr � ρr ρz ρz ρ ṙ Dt �r · κv + = 0 + ρr ρr ρr � = �(ω, S, z) Dt ω = 0 Dt s = 0 We’ll now chose specific vertical coordinates and see how the above equa­ tions simplify or otherwise turn out in each case. 10.2 Terrain following coordinates The representation of bottom topography in z-coordinates has historically been crude1 . An elegant approach to incorporate the topography smoothly into models is to normalize the vertical coordinate by the fluid depth. There are many choices of terrain following coordinate but we will demonstrate the simplest, known as �-coordinates (sigma) where z �= H(x, y) This maps the bottom at z = −H(x, y) to � = −1 and so the domain is made square. Using r = � = z/H simplifies the scaling factor ρr ρz 1 = or =H ρz H ρr In the hydrostatic equation, the vertical derivative of gz is trivial since z = rH = �H. The resulting �-coordinate equations are: 1 ρp 1 �� z = 0 Dtκv + f kˆ � κv + �� p − �H �¯ ¯ ρ� 1 The finite volume method (or partial cell representation) has somewhat alleviated this problem. However, there are still advantages to terrain-following coordinates such as fixed domain, the ability to use a modal representation in the vertical, etc... a) 0 0 −10 −10 b) −20 −20 −30 −30 −40 −40 −50 −50 −60 −60 −70 −70 −80 −80 −90 −90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 10.1: Schematic of sigma-coordinates (a) and s-coordinates (b) in physical x − z plane. ρp + �gH = 0 ρ� ρ �� · Hκv + H �˙ = 0 ρ� There are two gradient terms in the momentum equations. To a large degree, these cancel but cannot necessarily be canceled analytically. To see this, imagine a resting fluid of constant density. The pressure is simply a function of height, z. However, wherever the topography, H, slopes, the coordinate surfaces above must slope too. This means that �� p is non-zero and indeed is canceled by g��� z. At a numerical level, it is impossible to make these two terms balance exactly for an arbitrary static initial condition and so leads to “spontaneous” motion. The imbalance of terms is known as the “pressure gradient error” and a lot of effort has been put into reducing this error. With modern schemes, a typical error will be of the order of millimeters per second for common circumstances. Over shallow topography, the effective vertical resolution is increased. This can be an advantage but care must be taken to avoid numerical insta­ bility due to small grid lengths (CFL). The implicit treatment of vertical terms is therefore used in models such as SPEM/SCRUM or SEOM. There are many variants on the terrain following coordinates where dif­ ferent functions of bottom depth and height are used. Some also include the free-surface height in the normalization; the s-coordinate is defined as s= z−� H +� so that all though the free-surface can move the computational domain is fixed (−1 � s � 0). The surface variations introduce relatively minor slopes to the coordinate surfaces compared to those induced by the bottom topog­ raphy. Terrain following coordinates have several advantages over z-coordinate models: • Smooth representation of bottom topography. • Allow concentration of coordinate lines in “boundary layers”. but have the following dis-advantages: • Pressure gradient error can be significant and/or lead to spontaneous motion. • Representation of horizontal or along isopycnal processes is awkward. 10.3 Isopycnal coordinates Where the fluid is adiabatic, (potential) density is conserved and since, under statically stable conditions, density is a monotonic function of height it makes a useful vertical coordinate. The real advantage of isopycnal coordinates is their Lagrangian treatment of vertical motion; explicit advection acts only in the horizontal. This makes the models models very adiabatic allowing them to avoid numerical diffusion in the vertical that can be troublesome in z-coordinate and terrain-following coordinate models. Strictly speaking, the vertical coordinate is potential density, �ρ but we will use � as the vertical coordinate to avoid confusion with the terrain fol­ lowing coordinate. Setting r = �, the continuity equation becomes: ρt ρ� z + �r · ρ� zκv + ρ ρ� z�˙ = 0 ρ� where �˙ is non-zero only where diabatic terms force flow across isopycnals. The continuity equation has become a layer thickness equation for each den­ sity class. The hydrostatic equation is ρ ρp ρp ρ + � gz = + �gz − gz ρ� ρ� ρ� ρ� ρ �M − gz = 0 = ρ� where M = p/�+gz is the Montgomery potential. The momentum equations become 1 Dtκv + f kˆ � κv + �� �M = 0 �¯ Isopycnal models have several advantages over the height and terrain-following coordinates: • Ideal for modeling lateral transfer processes. Adiabatic motions mod­ eled without any spurious diabatic terms. • Smooth representation of topography. The bottom topography is rep­ resented as piecewise-linear and is included in the model through a vanishing of the layer thickness. • Conserves volume of density classes Some dis-advantages are: • Full or non-linear equation of state is difficult. • Non-hydrostatic effects/dynamics are not possible. • Density is not a natural coordinate for representing mixing processes such as the surface BBL (shallow and deep mixed layers). • Vertical and horizontal resolution are tightly connected in regions where isopycnals outcrop. This can lead to inadequate horizontal resolution in regions such as the ACC. a) 0 0 −10 −10 a) −20 −20 −30 −30 −40 −40 −50 −50 −60 −60 −70 −70 −80 −80 −90 −90 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 10.2: A schematic of (a) “step” topography in a z-coordinate model and (b) a partial cell representation. 10.4 z-coordinate models Note that the z-coordinate equations are readily recovered. Setting r = z means that ρr z = 1, Dt ρr z = 0, �r gz = 0 and ρr gz = g. Since we are so familiar with z-coordinates we will refrain from writing them out once more. One of the problems with z-coordinates used to be with the crude rep­ resentation of topography as a series of giant steps (see Fig. 10.2). Now, however, we represent the bottom either via a variable thickness bottom layer (partial cells) or with a piecewise-linear representation. This approach has been tested and favourably compared to results from a sigma-coordinate model. The three main z-coordinate models have some form of partial cells implemented. Some advantages of using z-coordinates are: • Simple to implement and use! • The full equation of state can be used. • Diabatic processes, including the mixed-layer, can be represented easily. • Non-hydrostatic and non-Boussinesq terms can be included. Some dis-advantages are: • Representation of along isopycnal processes are awkward. • Representing the bottom boundary layer is awkward. 10.5 Primitive equations in pressure coordi­ nates The other natural choice of vertical coordinate is pressure. Here, however, the non-Boussinesq equations are the best starting point. This leads to equations very similar to those used in meteorology. They have the advantage of not needing to “simplify” the equation of state dependence on pressure but have a major difficulty in representing the bottom. Use of the Philip’s normalized pressure coordinate (a sigma-coordinate) renders equations very similar to the terrain-following coordinate outlined above. The compressible, non-Boussinesq hydrostatic equations in height coor­ dinates are: 1 Dtκv + f kˆ � κv + �z p � ρp + g� ρz ρw 1 Dt � + �z · κv + � ρz � Dt ω Dt s = 0 = 0 = 0 = �(ω, S, p) = 0 = 0 Using the same transformation rules as before, the non-Boussinesq equations written in a general coordinate “r” become: 1 Dtκv + f k̂ � κv + �r p + �r gz � ρ ρp + � gz ρr � ρr ρz ρz 1 ρ ṙ + Dt � + �r · κv + Dt ρr ρr � ρr � Dt ω Dt s = 0 = 0 = 0 = �(ω, S, p) = 0 = 0 The continuity equation looks unusual. having two total derivatives. A judicious choice of coordinates, r = p, will cause these terms to cancel. With r = p, the hydrostaic equation becomes 1 + �g which is normally written ∂+ ρz =0 ρp ρ� =0 ρp where ∂ = 1/� is the specific volume and � = gz is the geopotential. We can now examine the two total derivatives in the continuity equation: Dt 1 1 ρz ρz 1 Dt � = −Dt ( ) − 2 Dt � + g� g� ρp ρp � 1 1 = Dt � − 2 Dt � = 0 2 g� g� The other major simplification is that the gradient of pressure on a pressure surface is zero (�p p = 0) so that the momentum equation retains only one gradient term. Making the above substitutions we arrive at the primitive equations in pressure coordinates: Dtκv + f kˆ � κv + �p � ρ� ∂+ ρp ρσ �p · κv + ρp ∂ Dt ω Dt s = 0 = 0 = 0 = ∂(ω, S, p) = 0 = 0 where σ = ṗ is the cross coordinate flow in pressure corodinates. The equa­ tion of state is now written in terms of specific volume. Examination of the above equations and comparison with the Boussinesq, hydrostaic equations in height coordinates reveals a one-to-one correspon­ dance of terms and variables; the equations are isomorphic. They take the same structural form and this carries through to the boundary conditions also (de Szoeke and Samelson, 2002; Marshall et al., MWR 2004). Because the equations are isomorphic, all of the algorithmic considerations and coordi­ nate transformations we discussed can be applied to the primitive equations in pressure coordinates. This is an approach used in the MIT GCM to model the atmosphere or non-Boussinesq ocean using the same code that is used to model the Boussinesq ocean. 10.5.1 Atmosphereic equations in pressure coordinates For the atmosphere, we simply need to replace the equation of state with that of the atmosphere (we’ll use the ideal gas equation, p = �RT ) and drop salinity from the equations. There is one simplifying step that can be made due to the ideal gas equation which uses a function called the Exner function, �: �� p � = cp po Here, cp is the specific heat at constant pressure, R is the gas constant and η = R/cp and po is a reference pressure used to define the potential temper­ aute: �−� cp T p ω= =T � po Note that �� η cp η p ρp � = = � p po p and RT ηcp T η�ω = ωρp � = = ∂= p p p so that the primitive equations for an ideal dry atmosphere can be written: Dtκv + f kˆ � κv + �p � ρ� ω+ ρ� ρσ �p · κv + ρp Dt ω = 0 = 0 = 0 = 0 Integrating the continuity equation vertically over the air column and apply­ ing appropriate boundary conditions yields the surface pressure equation ρt ps + � · psκv = 0 10.6 Hybrid coordinates The choice of vertical coordinate is perhaps the single most important fea­ ture that differentiates between models and is still an active area of research. Many groups are examining how to use “hybrid” coordinates where the coor­ dinate may be a function of height in the mixed layer, a function of isentropes in the interior and some function of the terrain in the bottom boundary layer.