VIEWS: 6 PAGES: 168 POSTED ON: 2/9/2012
Geophysical Fluid Experiments with the Princeton Ocean Model Lie-Yauw Oey Princeton University Page 2 COVER ILLUSTRATION A three-dimensional surface of near-inertial energy = 0.03 m2s-2 on Sep/03/12:00, approximately one week after the passage of the disastrous hurricane Katrina, Aug/25- 30/2005, in the Gulf of Mexico, USA, simulated by the Princeton Ocean Model. This shows penetration of intense energies to deep layers due to the presence of the warm-core ring and the Loop Current, both represented by the dark contours in the three cut-away xy-planes. The location of an observational mooring where extensive model- observational analyses have been conducted is shown as vertical dashed line [from: Oey et al. 2008: “Stalled inertial currents in a cyclone,” in press, Geophysical Res. Lett., with permission to reproduce]. Page 3 PROLOGUE: About This Book “.. I have no special talents. I am only passionately curious..” (Albert Einstein) This book is about exploring various aspects of fluid motions on a rotating earth using a popular, relatively simple yet powerful numerical model - the Princeton Ocean Model (POM). It is a book aimed primarily for advanced undergraduates and graduates in physical oceanography; but the book should also be useful for researchers or anyone fascinated by and intensely curious about oceanic fluid motions. Some knowledge of fluid mechanics is assumed, equivalent to a two-semester course which usually covers up to the derivations of conservation laws of viscous fluid motions, including the boundary layer theory; these prerequisites should not be an impediment to those interested enough to want to open this book. On the other hand, we have strived to make the book more or less self-contained by building each chapter from simple to more advanced concepts. By conducting geophysical fluid experiments on a computer and analyzing the results, we hope that the student will gain a solid understanding of geophysical fluid dynamics (GFD) and physical oceanography; the student will also learn to critically examine the results (be a skeptic) and to attempt connecting them to the real-world phenomena. And of course, the student will learn to use POM, numerical methods and, we hope, to keep exploring long after he or she surpasses this book. Why POM? One reason is our familiarity with the model. More importantly, however, it is because the model is relatively easy to tweak, say to suit a different flow problem, or to change the model physics. This makes it an excellent educational tool because it gives the student a hands- on experience. It is much like the difference between being a driver and being a driver as well as a mechanics. Most of us are the former, but only the privileged few belong to the latter. POM shares many of the same features of other popular ocean models now available in the scientific community, and it uses basically the same equations of motions and conservation of mass etc. Therefore, an accomplished student of this book should find it easy to transit to other models should he or she so desire. The book begins in chapter 1 with the near-surface ocean response to wind built upon Ekman’s fundamental work. This is not the common Page 4 approach of a book in physical oceanography. On the other hand, most of us experience the ocean perhaps through a summer swim in a coastal sea or in a lake, and wind-driven surface motions are the most apparent. For a majority of flows on the rotating earth, we find that by directly dealing with viscous boundary-layer (i.e. Ekman’s) flows, one can quickly grasp why most of the ocean’s interior is nearly inviscid yet why the boundary layers are so important for the interior flows. POM and related files for exercises outlined in this book can be downloaded from: ftp://aden.princeton.edu/pub/lyo/pom_gfdex/ Download latest public-released POM from: http://www.aos.princeton.edu/WWWPUBLIC/PROFS/waddow nload.html Download POM User Guide other releases from: http://www.aos.princeton.edu/WWWPUBLIC/htdocs.pom/ Page 5 CHAPTER 1: Wind-Driven Ocean Currents The book begins with the near-surface ocean response to wind. We will for the moment assume that the ocean density ( 1035 kg m-3) is constant. Throughout this book, unless explicitly stated, we will use the SI units. Most of us experience the ocean perhaps through a summer swim in a coastal sea or in a lake. Wind-driven surface motions in the form of waves and surfs are familiar. However, these motions can involve large vertical accelerations and are of very small scales (order of centimeters to meters, O(cm-m)) which we will not model for now, at least not directly. Instead, we focus on motions which are of sufficiently large scales in the horizontal (scale L ~ O(1-1000 km)) compared to vertical (scale D ~ O(10- 1000 m)), D/L << 1, that the pressure field satisfies the hydrostatic equation: , (1-1) where a Cartesian coordinate system xyz is used such that, conventionally, x and y are the west-east and south-north axes respectively, and z is the vertical axis with z = 0 placed at the mean sea level (MSL; Figure 1-1); x, y and z are measured in m. The pressure p is in N m-2 (also called Pascal, Pa), and g ( 9.806 m s-2) is the acceleration due to gravity. In other words, the motions of interest occur in such thin fluid layers (with D/L << 1) that as far as the vertical force balance is concerned the fluid may be treated as being static.1 For such an atmosphere-ocean system, the pressure at some depth z = z’ in the water is just the weight per unit area of fluid (water and air) above, or p = [mass per area].g = [ (-z’) + <eazea> ].g (1-2) 1 The thinness of the oceanic (or atmospheric) layer relative to its horizontal extent is comparable to that of the thickness of paper of this book relative to its width. Page 6 Figure 1-1. A sketch of water of mean depth H(x,y) and free surface (x,y,t). Here, z = (x,y,t) is the free surface; and/or z’ can be positive or negative, and z’ ( 0) represents the distance from the sea-surface to the point z = z’ in the water. The term <eazea> represents the mass per area of the atmosphere, it is given by <eazea> = (1-3) where a is the air density. Thus ea and zea may be thought of as the effective density and height of the “air-containing” atmosphere. The main bulk (about 80%) of the atmosphere’s mass is within the lowest layer (the troposphere) of about 10 km thick such that the atmospheric pressure at the sea surface, pa = <eazea>.g 105 N m-2 = 1 bar.2 Equation (1-2) can be written as 2 An easy way to remember this is to use ea 1 kg m-3, zea 10 km and g 10 m s-2. Note also the familiar weather reporting of “1000 millibar” etc, which is 1 bar, the approximate sea- level pressure. Page 7 p = g(-z’) + pa(x,y,t) (1-4) which is also obtained by integrating (1-1) from z = z’ to z = ; in general, the pa is a function of the horizontal position and time. If both a and pa do not vary with (x, y) and either are steady or at most only vary so slowly with time that the fluid remains in hydrostatic equilibrium, i.e. (1-1) is satisfied, then the fluid remains motionless provided that initially it is also at rest. A simple system to imagine is water in an annulus channel (Figure 1-2). If the axis of the annulus is far from the center, so that r/Ranu << 1, where r is the width of the channel, and Ranu is the radius of curvature of the annulus, then one can approximate the (motionless) water to be in a straight channel with x directed along the channel axis, y across the channel (Figure 1-2) and a vertical slice along the axis gives Figure 1-1. Modelers refer to this kind of channel a “periodic channel” since, starting from a yz-section at any x-location and going around along the axis, the field variables return to the same values. When conducting geophysical fluid experiments, a periodic channel is a convenient configuration to use because it often allows the extraction of essential flow physics while at the same time alleviates the modeler from having to formulate more complicated boundary conditions. Page 8 Figure 1-2. A circular annulus channel of radius Ranu much larger than the width of the annulus r. This is used to illustrate along-channel periodic fluid motion within the annulus. 1-1: A Simple Shearing Flow by the Wind Consider therefore the straight channel (Figure 1-1) in which we will further assume that all variables are independent of the cross-channel axis “y,” and that the earth’s rotation is nil. The flow is described by the three components of velocity (u, v, w) and the pressure p. Since there can be no flow across the channel wall, i.e. v = 0, at y = 0 and y = r, t. (1-5) the cross-channel velocity v must be zero everywhere; the symbol t means “at all time.” The channel is further idealized by letting its depth H = constant. Initially, at t = 0, the water is at rest. A wind stress (N m-2) is then applied uniformly at the surface, so that is a function of time only: = (t) (1-6) Similarly, we could stipulate that pa is also a function of time only; but for simplicity we will set pa = 0. (1-7) It follows that, the along-channel velocity u cannot vary with x. Therefore, at any point, there can be no accumulation (convergence or divergence) of mass, and since w = 0 at z = H, the vertical velocity is nil everywhere, and the sea-surface remains flat. Thus, u = u(z, t), w = 0 = /x = p/x. (1-8) Under these very specialized conditions, a parcel of fluid experiences acceleration only due to the vertical shear stress. The momentum balance is then: (1-9) Page 9 where zx denotes shear stress (force per unit area) in the x-direction acting on the fluid elemental face that is perpendicular to the z-axis, and D(.)/Dt is the material derivative which for a fluid property S is given by . (1-10) For the present specialized case, the last three terms are nil, and with S = u, we have = . A loose analogy is a stack of poker cards which are glued with a special adhesive that is never dry and thus remains sticky. The top card is then “pulled” parallel to the card’s surface a small distance and it drags upon the card below it. In our fluid system, the wind is doing the “pulling” by transferring air momentum onto the water’s surface (and we assume that no ripples or wind waves are produced!). Momentum is vertically transferred from the surface to fluid layers below – upper fluid drags upon the lower fluid which in turn drags upon the layer further below, and so on (Figure 1-3). Figure 1-3. The ocean is treated as a (vertical) stack of thin layers each of thickness z : shown are layers “1” and “2;” layer “0” is air. The Fmn is (shearing) force per unit (horizontal) area due to layer “m” on layer “n.” The formulae show that the net force per unit area in layer “1” is F1net = F01 + F21 =z.zx/z. Therefore, the net shearing force per unit mass is F1net x Page 10 y/(xyz) = (zx/z)/. This, in the idealized case described in the text, is = u/t, the acceleration in the x-direction. For the so called Newtonian fluid such as water, the shear stress zx is proportional to the vertical rate of change of fluid velocity (i.e. the vertical “rate of strain”): zx = (1-11) where is the viscosity which to a good approximation may be taken as a constant and moreover is rather small; it is 10-3 kg m-1 for water and 1.7×10-5 kg m-1 for air. Equation (1-11) is valid for laminar or slowly- moving viscous flows void of turbulent movements we usually see in say, a mountain stream or in swirling air vortices around a house. However, as will be seen below, one often model oceanic and atmospheric flows using a similar formula, except that an “eddy viscosity” e is used to represent the aggregated effects of small eddies on the large-scale flows. The eddy viscosity is much larger than the molecular viscosity and is moreover not a constant – it depends on the flow. The “ e” will be ascribed below but for now we will proceed with our model formulation using equation (1-11). Equation (1-9) then becomes: (1-12) where e = e/, the kinematic eddy viscosity. This is a simple “heat- diffusion” equation that can be easily solved subject to some initial conditions and boundary conditions at z = 0 and z =H; for example: u = 0, t=0 (1-13) u = 0, z = H (1-14a) Page 11 , z = 0, (1-14b) where o is a constant wind stress. Exercise 1-1A: Use POM to solve (1-12) subject to (1-13) and (1-14) in a periodic straight channel. Compare your solution with exact (analytical) solution. Experiment with different grid sizes as well as with more complicated initial and/or boundary conditions. Analyze and discuss your results (e.g. use Matlab or IDL etc). Exercise 1-1B: Repeat Exercise 1-1A using a e-profile that decays with depth, e.g. e = eo exp(z), where eo and are constants. Exercise 1-1C: Repeat Exercise 1-1A using an annulus channel (Figure 1- 2) instead of the straight (periodic) channel. Use a curvilinear grid (in POM) for this exercise. Experiment with various values of Ranu and r. Compare and discuss your results. Are the results substantially different from the straight-channel case when Ranu is “not too large” (define what this means), why? Is equation (1-12) still valid then? Why (or why not)? 1-2: Effects of Earth’s Rotation In the above straight-channel flow example, we saw that the cross- channel flow is nil. This is no longer true if the channel is placed on a turntable and rotates about a vertical axis. We now expand the study of wind-driven shearing flow to the case when the earth’s rotation cannot be ignored. We can imagine that the earth is the “turntable.” If our channel is placed at the (celestial) pole, then it rotates about the earth’s rotation axis once every day (0.997258 day to be more exact) anticlockwise at the north pole and clockwise at the south pole (Figure 1-4). Thus the rotational frequency = 2/(86163.09 s) = 7.292 10-5 s-1. If the channel is placed on the (celestial) equator, then there is no rotation (though the channel is being carried west to east at a great speed = Rearth 1673 km/hour; Rearth 6371 km is the radius of a sphere having the same volume as the earth). At a given latitude, the rotational frequency is: Rot() = .sin(), /2 +/2 (1-15) Page 12 with the convention that positive (negative) is anticlockwise (clockwise). Figure 1-4. The Earth's axial tilt (or obliquity) and its relation to the rotation axis and plane of orbit (from Wikipedia). We are anchored (by gravity) to the earth and therefore also rotate as the channel does. It is most convenient to describe phenomenon referenced to this rotating coordinate frame rather than to a fixed, inertial frame. For scalar quantities such as the temperature or salinity of a (generally moving) fluid parcel, a change in the frame of reference clearly will not change their value. The parcel’s momentum is a different story. While the phenomena are still independent of the frame of reference, our perception of the phenomena as well as their descriptions are altered. An example is offered in Figure 1-5. Page 13 Figure 1-5. A turntable rotating at a constant radians/s with two observers “A” and “B” on it. At time t o “A” throws the ball (red) to “B.” A time-interval “t” later, the relative position of both observers are unchanged (and neither perceive any change). Let the ball’s speed be such that it arrives at the original B’s position in the same time-interval. To “A,” the ball has been forced to the right of “B.” This apparent force is called the Coriolis force. It is only perceived by “A” (and “B”) in the rotating frame but not felt by an outside observer in an inertial frame. We now derive the momentum equation on a rotating frame. We examine first the rate of change of an arbitrary vector (e.g. the momentum) B = B1i1 + B2i2 + B3i3 in a rotating Cartesian frame of reference with unit vectors (along x, y and z) i1, i2 and i3 [Pedlosky, 1979]. By chain rule: (1-16) where repeated index k means summation over all k = 1, 2 and 3, and the subscript I denotes rates of change as seen from an observer in the inertial frame, respectively. Clearly, the first term on the RHS of (1-16) is (1-17) i.e. it is the rate of change of Bk when ik is fixed, i.e. when the observer is himself rotating, hence is not able to sense the rotation of the coordinate axes. The second term on the RHS of (1-16) involves , the rate of change of a vector of fixed length (i.e. a unit vector), hence its change is completely due to its changing direction as the coordinate axes rotate. From Figure 1-6, we see that (1-18) where is the angular velocity at which the constant-length vector A rotates. Setting A = ik and using (1-17) and (1-18), equation (1-16) gives (1-19) Page 14 Figure 1-6. Rate of change of a vector A of fixed length caused by it being rotated by . This rate of change is seen from the sketch to be a vector of length |||A|sin(), where is the angle between A and , directed perpendicular to the plane formed by A and , i.e. it is simply the vector product A. Thus the rate of change of the vector B in the inertial (non-rotating) system consists of two parts. The first part, , is the rate of change as seen by an observer (that is us) “glued” to the rotating frame of reference (that is the earth). To this must be added the second part, B, in order that the correct rate of change, , can be used in Newton’s law to describe fluid motion. This equates the rate of change of momentum following a fluid element in an inertial frame of reference to the sum of forces acting on the element: Page 15 (1-20) where the subscript I again reminds us that the equation is valid only in an inertial frame of reference. The RHS is the sum of the pressure gradient force p, the body force , and the viscous force F(uI) which for Newtonian fluids with molecular viscosity is (1-21) These equations are derived, for example, in Batchelor [1967]. Applying (1-19) to the position vector r of the fluid element, and also to uI (i.e. setting B = r and B = uI in turn), we obtain, (1-22a) (1-22b) where in (1-22a) we have set and . As seen in the inertial frame (from a “fixed star” outside the planet earth), the fluid velocity uI consists of the velocity uR measured by an earth-bound observer, plus an amount equal to the rate r (in m/s) at which the fluid element is being carried along by the rotating earth. Since uR is the velocity we directly observe, we substitute (1-22a) into (1-22b), (1-21) and (1-20) to express the equation of motion entirely in term of uR. Thus, (1-22b) becomes: , (1-22b) where we have assumed that is small compared to other terms. Figure 1-7 explains that (r) is a centrifugal acceleration that can be expressed as the gradient of a centrifugal potential c; it therefore can be grouped with the body force on the RHS of (1-20), i.e. ( + c) = g, the acceleration due to gravity. The direction of g is conveniently taken to align with the z- Page 16 axis, positive upwards. Since |(r)| < ||2 Rearth 3.410-2 m s-2, the centrifugal acceleration makes but a small correction to the radial direction that is “vertical” were the earth not rotating. Note that (r) is zero at the poles and is maximum at the equator; it accounts for the slight bulge of the earth around the equator. Substituting (1-22a) into (1-21), one can show that F(uI) = F(uR). Thus substituting (1-22b) into (1-20) gives: (1-23) which is entirely in terms of quantities in the rotating frame of reference. As written, equation (1-23) is valid also for non-constant . Henceforth the subscript R will be dropped. Figure 1-7. (a) Sketch showing the centrifugal acceleration vector (r) = ||2|r|sin()nc of a fluid material volume element with Page 17 position r on the earth’s surface (at latitude = /2 - ), directed along the unit vector nc perpendicular to and away from the axis of rotation of angular momentum . For convenience, the origin of r is at the earth center, but this is irrelevant. (b) Sketch showing the vector sum ( + c), of (i) , the gravitational acceleration vector directed towards the center of the earth, where is the gravitational potential, and (ii) c =(r) =||2r [see sketch (a) for meaning of r], where c is the centrifugal potential. Let r = xk ik , then c = (||2x2k )/2 as direct differentiation using = ik (/xk) shows. Note also that c = (||2|r|2)/2 = |r|2/2. In geophysics, the summed vector ( + c) = g is called the acceleration due to gravity and the direction of g (i.e. “upwards”) is conveniently aligned with the z-axis. The Equation of Motion Applied to Large-Scale Oceanic Flows: Written in component form, (1-23) is: (1-24a) (1-24b) (1-24c) In arriving at these equations, we have used the incompressibility condition (1-25) to eliminate the second term on the RHS of (1-21). Equation (1-25) is also often called the continuity equation, a term that we will also use, though this is correct only if effects of density change on mass balance are negligible. Instead of , we have also used the more conventional notation to denote the substantive derivative: Page 18 (1-26) i.e. the rate of change of a field variable G following the fluid’s material volume element. This is permitted since the vector r (Fig.1-7) denotes the position of the fluid element. A localized coordinate has been chosen such that the direction of 23 points in the z-direction and its magnitude is (c.f. Figure 1-7) 2sin() = f, say, where is the latitude of the fluid element. The “f” is called the Coriolis parameter; it is zero at the equator, and is a maximum (=2) at the pole as noted previously when discussing equation (1-15). For large-scale, thin-fluid flows (i.e. D/L << 1) that we will mostly focus on, it can be shown [e.g. Pedlosky, 1979; Gill, 1982] through a scaling analysis that the terms in (1-24) are generally quite small compared to other terms in the same equation; these terms will therefore be dropped. The molecular viscous terms ( u) are also very, very small; though their importance increases as the scale of motion becomes very small and eventually they provide the necessary energy dissipation for the fluid system. In models, effects of turbulence are usually parameterized by invoking an eddy or effective viscosity which in the vertical (z) is denoted by M (or diffusivity KH for heat and salt, chapter 3; two excellent texts are Hinze, 1961, and Schlichting, 1963); these have much enhanced values than their molecular counterparts: KM >> , a result of “eddy mixing” on the small scales.3 Unlike the molecular viscosity, the eddy viscosity depends on the flow itself, hence is a function of time and space, and in general takes different forms in the vertical and horizontal directions (i.e. non-isotropic). Various models of turbulence are available for KM (and KH). For example, the Mellor and Yamada’s [1982] parameterization scheme is very popular and is used in many oceanic and atmospheric circulation models including 3 A proper treatment involves Reynolds averaging – see Hinze, 1961. Equation (1-24) is then called the Reynolds equation. Page 19 the Princeton Ocean Model. While parameterizations for the vertical eddy viscosity and diffusivity are relatively well-known, those for enhanced horizontal eddy viscosity and diffusivity are in comparison still a subject of much debate. Unless otherwise stated, we will assume that horizontal viscosity and diffusivity are small, and include them only for the purpose of numerical computations: stabilization of the numerical schemes and/or elimination of grid-point (“2x” oscillatory) noise [e.g. Roach, 1972]. Excellent reviews of turbulence mixing in geophysical flows can be found in the various chapters of the book edited by Baumert et al. [2005]. With the above simplifications, equation (1-24c) becomes (1-27) which is the hydrostatic equation (1-1) except that it is seen to be valid also for non-uniform . Indeed, (1-24) admits a “static” solution u = 0, = o(z) and p = po(z) (1-28) which however is not very interesting. More interesting flows are produced as a result of deviations of the pressure field from this hydrostatic state: = o(z) + ’, p = po(z) + p’ (1-29) where the ’ and p’ are perturbation density and pressure respectively. In the ocean o oc = 1025 kg m-3 say, where oc is a constant density, and ’ is generally a small fraction of this: ’/o 0.01 and smaller. Subtracting (1-28) from (1-24), letting o in the acceleration terms, and dividing through by o, we obtain: (1-30a) (1-30b) Page 20 (1-30c) with an error in (1-30a,b) that is at most O(’/o). 1-3: Wind-Driven Homogeneous Flows with Rotation Equations (1-25) and (1-30) are insufficient to solve for u, v, w, ’ and p’. Additional equations involving the heat and salt equations as well as the equation of state are required. These will be discussed in chapter 3. Here, we study the case of a homogeneous fluid of constant density (i.e. ’ is known, ’ = 0), for which equations (1-25) and (1-30) are then complete and may be solved for u, v, w, and p’, given appropriate initial and boundary conditions. For this homogeneous fluid system, while ’ = 0, p’ is not and it may be written in terms of the free surface . Assume then a free-surface z = (x, t) in an idealized ocean of depth z = H but with no lateral boundaries. There can be no flow across the surface and bottom, i.e. z = and z = H are material surfaces, and D(z)/Dt and D(z+H)/Dt are both = 0: (1-31a) (1-31b) These equations constitute the upper and lower boundary conditions for our problem. From (1-30c) with ’ = 0, we have p’ = p’(x, y, t); also the first of (1-29) gives = o, a constant. Integrating (1-27) and using the second of (1-29) yield: p = ogz + function(x, y, t) = po(z) + p’(x, y, t), (1-32a) Page 21 from which we may let po = ogz (1-32b) Across the air-sea interface, z = , the total pressure p must be continuous and equal to pa, so that (1-32) gives: p’ = og + pa. (1-33) The incompressibility condition (1-25) can now be used to relate (or p’) to the velocity. Integrating (1-25) from z = H to z = and applying (1-31): (1-34a) (1-34b) Equation (1-34) can also be re-written in term of a vertically-averaged velocity , where , so that . The form used in (1-34) shows explicitly however that the equation is valid for general (u, v) that may also be a function of z (as well as of x, y and t). Substituting (1-33) into (1-30a,b), the resulting equations and (1-34) are grouped together here for convenience: (1-35a) (1-35b) . (1-35c) These constitute three equations to be solved for (u, v, ) given appropriate initial and boundary conditions. These are in fact the equations solved numerically in POM. After (u, v) are solved, the incompressibility condition (1-25) is then used to solve for w. Before we describe the numerical solution however, it is instructive to examine more closely how one can approach the problem analytically. Page 22 1.3.1 Boundary-Layer and Interior Flows: For analytic treatment, we closely follow Pedlosky [1979]. We study a problem in which the ocean is confined in the vertical by a free surface z = (x,y,t) and a bottom z = H(x,y), but is unbounded in the horizontal. It is more direct to use (1-25) and (1-30), which are repeated here: (1-36a) (1-36b) (1-36c) (1-36d) Under each term, we indicate its order of magnitude by using the following scale for each variable: (x,y) L; z D; (u,v) U; KM K; (1-37a) t L/U; w W ~ UD/L; (1-37b) p’ P ~ oLfoU; f fo (1-37c) where “” means “have/has the scale of.” Expression (1-37a) sets the basic scales for the horizontal (L) and vertical (D) distances, and also the horizontal velocity (U) and the eddy viscosity (K). The time scale (L/U) in (1-37b) is chosen to be the “advective time scale,” i.e. the time taken for a parcel of fluid to cover a distance of about L. If f is constant, then: f/fo = ±1 (1-37d) where the plus (minus) sign is for the northern (southern) hemisphere. The vertical velocity scale W is chosen as follows. In order to satisfy the continuity equation (1-36a), horizontal divergences or convergences Page 23 must be balanced by vertical motion such that their net sum is zero; thus we set terms (1) and (2) in (1-36a) to be of the same order, and W ~ UD/L. It follows then that the scale of term (2) of equations (1-36b,c), WU/D ~ U2/L, which is also the scale of term (1) of these equations; although |w| is in general much smaller than |u| or |v|, it multiplies the vertical gradient or which is generally much larger than the horizontal velocity gradients such as . The scale for pressure in (1- 37c), P, is such that the horizontal pressure gradients balance the Coriolis acceleration (see below). It is easy to see now that when the variables in (1-36) are non-dimensionalized by the above scales, the equations become [after dividing (1-36b,c) through by the scale (foU)]: (1-38a) (1-38b) (1-38c) (1-38d) where now all variables are dimensionless, f = ±1 for constant Coriolis and (1-39) are the two parameters of the problem (the factor “2” in Ev is for convenience, below). The non-dimensionalized kinematic boundary conditions at z = and z = H have the same form as equations (1-31a,b). Equations (1-38b,c) show that the Coriolis and pressure gradient terms are of the same order, both of these terms in (1-38b,c) are multiplied by “1.” Thus by choosing the scales as those in (1-37c), we implicitly assume that the Coriolis acceleration terms are important, and are in fact of the same order of magnitude as the pressure-gradient terms. Clearly this choice is only valid if “f” is not zero; the scaling (and balance of terms) will need to be changed close to the equator. Observations of most oceanic (and atmospheric) flows do indicate that away from the equator, pressure gradients approximately balance Coriolis accelerations. From (1-38b,c) we Page 24 see that the conditions for this to be true are that the Rossby and Ekman numbers are small. The former requires that either the flow is not strong and/or the horizontal scale L is sufficiently large – for example, U 0.1 m/s, L > 10 km for fo 10-4 s-1. The latter requires that the Ekman scale (2K/fo)1/2 is much smaller than the ocean depth D; if K 10-4~10-1 m2 s-1, this requires that D > 10~50 m. However, close to a boundary, no matter how small the eddy viscosity KM is, at some point the shear terms become important. This is most apparent at the bottom boundary where the fluid parcel cannot have a relative motion with respect to that boundary, i.e. the “no-slip” condition must be satisfied: nt u = 0, or approximately u = 0, at z = H(x,y) (1-40) where nt is the unit vector tangential to the bottom. Therefore, unless the flow is trivially zero, the interior (i.e. away from the bottom) velocity of a real-fluid (i.e. fluid with viscosity) flow must decrease in value until it becomes zero at the bottom.4 Since Ev is small, the only way the viscous terms can become large enough to balance the Coriolis acceleration term is if the z-derivatives of (u, v) become very large, i.e. if the decrease of velocity from the interior to the bottom occurs within a thin layer in the immediate neighborhood of the bottom. This thin layer is called a boundary layer, or a bottom boundary layer in the present case. A similar thin layer also exists near the surface where wind stresses are applied. Geostrophic Nearly-Inviscid Interior Flows: Far away from boundaries, then, one expects that the flow is nearly inviscid and one may approximate equation (1-38b,c) by dropping the term 4 An inviscid flow needs only to satisfy the no-normal flow condition equation (1-31b), so that a fluid parcel can ‘slip’ past the bottom. Page 25 involving Ev. Just exactly how far is “far” will be determined below. If is also small, equations (1-38b,c) then give (for f = ±1): (1-41a) (1-41b) Oceanographers call the horizontal velocity (ug, vg) that satisfy (1-41) geostrophic velocity, and the equation geostrophic relation. In the northern hemisphere, a higher pressure in the south (east) than in the north (west), p’/y < 0 (p’/x > 0) would produce an eastward (northward) geostrophic velocity ug > 0 (vg > 0). The geostrophic flow around a center of high (low) pressure is therefore clockwise (anticlockwise) or anticyclonic (cyclonic). The direction of the flow is reversed in the southern hemisphere (f < 0). For homogeneous fluid flow, since p’ is independent of z (see 1-38d), (ug, vg) is therefore also independent of z. From (1-41), a geostrophic flow also has zero horizontal divergence: , (1-42) so that the corresponding vertical velocity wg satisfies: (1-43) All three components of the geostrophic velocity (ug, vg, wg) are therefore independent of z, the rotation axis. This is the Taylor-Proudman theorem which is valid for inviscid homogeneous slow (small ) flow. If the upper or lower boundary is flat, so that wg = 0 there, then wg = 0 through the water column. The motion is then two dimensional, (ug, vg) 0, but wg = 0. Figure 1-8 shows an example of this type of flow. (See Exercises). Strictly speaking, the model solution in Fig.1-8 is still evolving in time, so it is not steady and does not exactly satisfy equations (1-41), (1-42) and (1-43). This is an important point – purely geostrophic flows satisfying exactly (1- 41) through (1-43) are indeterminate. Any p’(x,y) satisfies the hydrostatic relation (1-38d), and (ug, vg) can be constructed to satisfy the two (approximate) momentum equations (1-41a,b); the unknown (ug, vg) in turn Page 26 is horizontally non-divergent (equation 1-42) which simply requires that the vertical velocity wg is independent of z. We know all these information (hence the Taylor-Proudman theorem) but cannot pin down what the flow is! To arrive at the unique model solution shown in Fig.1-8, we actually included the small terms on the left hand side of momentum equations (1- 38b,c), the time-dependent terms in particular.5 Another way to pin down a unique solution is to include the shear terms which of course have the added benefit that we can then also satisfy the physically relevant no-slip bottom condition (1-40) and the conditions at the sea surface where windstress may be applied. 5 There were additionally horizontal viscous terms on the RHS of (1-38b,c) included for numerical stability, but these were small and unimportant to the arguments presented here. Page 27 Figure 1-8. Nearly-steady homogeneous (constant-density) flow in an x- periodic channel, 600km by 300km, of depth 200m except at the channel’s center where a cylinder rises 50m above the bottom. Color is sea-surface height in meters and vectors are velocity at (A) z=0m (i.e. surface), (B) z=90m and (C) z=180m. This model calculation was carried out for 100 days when the flow has reached a nearly steady state. The experiment is meant to illustrate Taylor-Proudman theorem: flow below the cylinder’s height goes around the cylinder while above it the flow also tends to go around as if the cylinder extends to the surface. The velocity does not vary with “z” and the vertical velocity (which is not shown) is nearly zero. Viscous Boundary-Layer Flows: Page 28 It was mentioned before that no matter how small the Ekman number Ev is, at some point close to the boundary the shear terms become large enough that they balance the Coriolis and pressure gradient terms which by our scaling (see equation 1-38) are of O(1); purely geostrophic flows then break down. Thus in a thin layer close to the boundary, the small Ev multiplies to make their product O(1). We express this mathematically by changing the coordinate z to , where (1-44) so that the shear terms become , (1-45) which is of O(1). The heuristic interpretation is, since is large, of O(Ev-1), inside the boundary layer because a small change in z results in a large change in (u, v), the change is made more gradual in the stretched coordinate (i.e. ‘one moves more slowly’ in than in z) and one can resolve or “see” the rapid velocity change so that the boundary condition (1-40) (and a corresponding one at the surface, see below) can be satisfied. The “z” measures distance changes outside the boundary layer, i.e. “z” is an outer coordinate (see equation 1-37a, in which the dimensional z* (where the asterisk denotes dimensional) is non-dimensionalized by the water depth D >> boundary layer depth), and the geostrophic relation (1-41 through 1- 43) is an outer solution, though at the moment an incomplete one. Inside the boundary layer, vertical distances are measured by scaled by a much smaller (localized) length scale l* which from (1-44) is . (1-46) The is called an inner coordinate, and the solution to be derived below is called an inner solution. Indeed, equation (1-38) belongs to a class of problems that are amendable to solution by the singular perturbation Page 29 methods, the early development and usage of which is in aerodynamics [van Dyke, 1964]. With the new vertical scale l*, we have , which would appear to be much larger than (and therefore cannot be balanced by) the horizontal divergence term in the continuity equation (1-38a). This apparent dilemma is resolved by recognizing that in actuality the vertical velocity w is very small. Indeed, inside the boundary layer, the continuity equation requires that W*/l* ~ U/L, where W* is the scale of w inside the boundary layer; thus (1-47) where we have used (1-46) and (1-37b; for W). Therefore, because of the thinness of the boundary layer, hence a much smaller aspect ratio l*/L than the D/L previously assumed in (1-37), the appropriate vertical velocity scale W* is correspondingly much smaller, by O(Ev1/2), than the scale W that we have assumed. In the absence of other processes, the reader should satisfy himself or herself to show that w* is the scale of the vertical velocity that is valid both inside and outside the boundary layer. If we now write the newly rescaled vertical velocity as = w*/W* = w.W/W* = w/Ev1/2 (1-48) using (1-47), the terms on the RHS of the momentum equations (1- 38b,c) become: (1-49) and equation (1-38) becomes: (1-50a) (1-50b) (1-50c) (1-50d) Page 30 where (1-45) has been used in the viscous shear terms in (1-50b,c), f = ±1 for constant Coriolis and the hydrostatic relation is unchanged but for a trivial change of the vertical coordinate. Because of (1-50d) (or 1-38d), the thinness of the boundary layer means that the pressure-gradient force (per unit mass), (p’/x, p’/y)/o, may be represented by their values just outside the boundary layer in the geostrophic interior of the water column, i.e., by f(vg, ug) from (1-41) (f = ±1). Now that equation (1-50) is properly scaled for flows inside the boundary layer, we may simplify it (as we did with equation (1-38) when deriving the geostrophic relation) by dropping the terms involving : (1-51a) 0 (1-51b) 0 (1-51c) where (1-41) has been used in (1-51b,c) to replace the pressure-gradient terms by the geostrophic velocity. Also, for constant Coriolis, f = ±1, where the first (second) or top (bottom) sign of or applies to the northern (southern) hemisphere with positive (negative) Coriolis. Note that since (ug, vg) is independent of z (or ), the (u, v) appearing in the - derivatives in equations (1-51b,c) may be replaced by (uE, vE) = (u-ug, v-vg), (1-52) the Ekman velocity; this is in fact what we will do below. A wind stress is applied at the surface, and the boundary condition in dimensional form is: , at z* = 0 (1-54) where the kinematic wind stress (Xs*, Ys*) = (o*/o) (Xs, Ys), and o* is the scale of the wind stress in N m-2. In non-dimensional form, (1-54) is: , at =0 (1-55a) Page 31 where s = (Xs, Ys) is the non-dimensionalized kinematic wind stress, and . (1-55b) The boundary condition at the bottom is the no-slip condition: (u, v) = (0, 0), at z = H (1-56) Surface Boundary Layer with Constant Eddy Viscosity: We assume constant eddy viscosity KM (= 1, non-dimensional), and also that the depth of the channel is sufficiently deep that the surface and bottom boundary layers do not interact. Equations (1-51b,c) are solved by first re-writing them in terms of A = uE + ivE, where i = (1). Thus {(1- 51c) i (1-51b)} gives (for f = ±1): . (1-57) We first consider the surface boundary layer, for which the solution that satisfies (1-55a) and that vanishes as ~ (so that u ~ ug and v ~ vg as one moves into the geoestrophic region outside the boundary layer) is: (1-58) or, (1-59a) (1-59b) In dimensional form, equation (1-59a,b) are: (1-60a) (1-60b) where . Equation (1-59) shows that over a distance of about E* = l* = (1-61) Page 32 the velocity profiles (u, v) ~ (ug, vg), the geostrophic velocity in the ocean’s interior outside the boundary layer. The E* is called the Ekman depth; it measures the thickness of the boundary layer which is therefore thicker for larger K and/or smaller fo. At the surface, = 0, the vector product (uE, vE) (Xs, Ys) gives the angle w between the velocity (uE, vE) and wind-stress vector (Xs, Ys) as: w = sin-1[cos(/4)] = /4; (1-62) the positive (negative) sign is for the northern (southern) hemisphere and it means that the Ekman velocity is directed to the right (left) of the wind stress, at an angle of 45o. To obtain w, we use the continuity equation (1-51a), and also (1-59): (1-63) which gives upon integration and setting : (1-64a) In dimensional form, using equation (1-48), we have: (1-64b) where = z*/(2K/fo)1/2. Outside the surface boundary layer, ~ , the vertical velocity is: (1-65) In terms of the original variables “w” and “z,” we note from equation (1-44) that the ~ means a finite but small z-distance below the surface, i.e. ~ is equivalent to z ~ 0 and Ev1/2 << 1, (1-66) so that: (1-67a) Page 33 In dimensional form, this is: (1-67b) Equation (1-67) is the main result from the analysis of surface boundary layer: the effect of the wind is to induce a small vertical velocity which is directly proportional to the wind stress curl: (1-68) For positive wind stress curl, there is upwelling (downwelling) in the northern (southern) hemisphere. The vertical velocity is small, but as we shall see in a simple example, it has a profound effect on the interior flow field. Before we give this example, it is necessary to conduct a similar analysis on the bottom boundary layer near z =H. Bottom Boundary Layer with Constant Eddy Viscosity: The analysis for the bottom boundary layer follows much the same way as for the surface layer. Instead of the wind-stress condition (1-55), we now apply the no-slip condition (1-56). Let zb = z + H (1-69) where the subscript “b” (here and below) is used to denote the bottom boundary layer, to distinguish the solution from that obtained previously for the surface boundary layer. The bottom boundary layer is near zb > 0. Equation (1-56) then becomes: (u, v) = (0, 0), at zb = 0 (1-70) The boundary-layer coordinate (1-44) becomes: (1-71) Equation (1-57) is the same, but a solution that decays as b ~ + is chosen: uEb + i vEb = (ug + i vg) exp[(1 i)b ] (1-72) which also satisfies the no-slip condition (1-70). Thus, (1-73a) (1-73b) Page 34 and the corresponding dimensional form simply has the u’s and v’s replaced by corresponding variables with asterisks: (1-74a) (1-74b) The vertical velocity is again obtained from the continuity equation (1-51a), and also (1-73): (1-75) which gives upon integration and setting : (1-76a) In dimensional form, using equation (1-48), we have: (1-76b) where b = (z*+H*)/(2K/fo)1/2. In terms of the original variables “w” and “z,” we note from equation (1-71) that the b ~ + means a finite but small z-distance above the bottom, i.e. b ~ + is equivalent to zb ~ 0+, or z = H+ and Ev1/2 << 1, (1-77) so that (1-76a) and (1-48) give: (1-78a) In dimensional form, this is: (1-78b) Equation (1-78) is the main result from the analysis of bottom boundary layer. Unlike the surface boundary layer, it is the interior geostrophic flow that is driving the bottom boundary layer. In the northern (southern) hemisphere, a positive curl in the geostrophic velocity produces upwelling. Page 35 A Simple Problem: Quasi-Geostrophic Interior Solution: Equipped with the above boundary-layer formulae, we are now ready to revisit the interior flow and attempt to close the otherwise indeterminate (geostrophic) solution (equations 1-41 through 1-43). The correctly scaled equation that governs the interior flow is (1-38). We will find an approximate solution which as will become clear is called quasi-geostrophic. We already saw that a first approximation is the geostrophic relation, and indicated that the solution in Figure 1-8 was obtained through the inclusion of small terms on the left-hand-side of the momentum equations. These small terms represent small imbalances to the purely geostrophic relation. Equations (1-38b,c) suggest then that the “degree” of this imbalance, or ageostrophy, is of O( ). It is reasonable to assume then that the solution consists of predominantly the geostrophic part, plus a small correction that is of O( ): u = ug + u1 + O( 2); v = vg + v1 + O( 2); etc. (1-79) and similarly for other variables also. Higher order correction terms , O( 2) and higher are also indicated in equation (1-79), since there is no reason to expect that the geostrophic relation plus the O( ) correction (i.e. u1 etc) will provide the complete solution. It is important to understand that just as the (ug, vg, …) are properly scaled of O(1), the (u1, v1, …), (u2, v2, …), etc are also all scaled to be of O(1). Only then it makes sense to assume the form of (1-79) for (u, v, …). Equation (1-79) is called a perturbation expansion, which by the way may not be necessarily convergent [van Dyke, 1964]. From the discussion prior to Figure 1-8 (Taylor-Proudman Theorem) the geostrophic velocity components are all independent of z: (ug,vg,wg)/z = (0, 0, 0). In particular, for “w,” we already know from the boundary-layer solutions at either the surface or the bottom boundary (see 1-67a and/or 1- 78a) that w ~ O(Ev1/2). Therefore, since Ev1/2 << 1 ~ O( ) say, the corresponding expansion (i.e. 1-79) for w = wg + w1 + … shows that wg = 0, for all z; (1-80) Page 36 otherwise, the solution for w will be inconsistent with its small values near the surface and bottom boundaries. In quasi-geostrophic approximations, the vertical velocity component is very small, O( ; Ev1/2), and represents only the ageostrophic portion of the flow. However, this ageostrophy is very important. Substituting (1-79) into (1-38b,c), and collecting and equating terms that are of O(1) and O( ) separately we obtain for the O(1) terms the geostrophic relation as before, while the O( ) terms give: (1-81a) (1-81b) (1-81c) (1-81d) Note that the viscous terms and the vertical advection terms are missing in (1-81) because of course, but in the latter case also because wg = 0 from (1-80). Thus the O( ) dynamics are also hydrostatic and incompressible, but are not in perfect geostrophic balance. The ageostrophy is due to the unsteadiness and advective terms by the geostrophic flow (i.e. the left-hand-side of 1-81b,c). We eliminate the pressure gradients from (1-81) by subtracting the y- derivative of (1-81b) from the x-derivative of (1-81c), i.e. by [(1- 81c)/x(1-81b)/y], and using (1-42) and also (1-81a): (1-82) Page 37 where g = vg/xug/y is the relative vorticity of the geostrophic flow (not to be confused with the scaled vertical coordinate of the boundary layer, e.g. equations (1-44) or (1-71)). Since g is also independent of z, equation (1-82) can be integrated to yield, upon using (1-67a) and (1-78a): (1-83) Given the wind stress curl, s, this equation now allows the determination of the geostrophic vorticity field g, hence also the geostrophic velocity (ug, vg). It is the consideration of the viscous boundary layers near the surface (providing the wind stress curl) and also near the bottom (providing the ‘damping’, g) that allows the existence of a small vertical velocity and the determination of the O(1) geostrophic field. A simple example further illustrates these ideas. A simple problem: Consider a homogeneous fluid in a channel of constant depth = H that extends to in the x-direction; all fields (and forcing) are assumed independent of the y-direction. A y-directed kinematic wind stress s = (Xs, Ys) is applied at the surface: (Xs, Ys) = (0, x) (1-84a) The solution is seeked only in the region L* x* +L*, where L* = 500 km, or in terms of the non-dimensional “x,” in 1 x +1. Taking the scale of the kinematic wind stress o*/o = 10-4 m2 s-2: ; (1-84b) Page 38 we then have wind stresses of about 0.1 N m-2 at x* = 500 km. Since the flow is independent of y, the geostrophic relation (1-41) requires that ug = 0. Assuming also steady state, the LHS of equation (1-83) is then zero, yielding the simple solution: (1-85a) or in dimensional form: (1-85b) The solution (equation 1-85) states a simple balance between the applied wind stress at the surface and the frictional damping dues to the viscous effects at the bottom. Knowing the geostrophic velocity in the interior, the complete solution is supplemented near the surface by equations (1-59) and (1-64a) (or in dimensional forms, equations 1-60 and 1-64b) and near the bottom by equations (1-73) and (1-76a) (or in dimensional forms, equations 1-74 and 1-76b). Some details of this solution are given in Figures 1-9 and 1-10 (see exercises). Exercise 1-3A: Use POM to solve the simple problem (above) of wind- driven flow in a channel. (See Figures 1-9 and 1-10). Modify the code to solve also the case with no-slip conditions at the bottom, and replot the figures for this latter case. Exercise 1-3B: Use POM to solve the Taylor-Proudman column (Figure 1-8). Note that in this case POM is modified without vertical viscosity. Solve also the case of a seamount instead of the cylinder. Exercise 1-3C: Use POM to solve a case that has wind as well as the seamount: i.e. a combination of exercises 1-3A and 1-3B above. Derive the corresponding analytical solution. Page 39 Figure 1-9. Nearly-steady homogeneous (constant-density) wind-driven rotational flow in a channel of depth 200m and theoretically unbounded in x, but in practice an x-length = 1000 km is used in the solution. The flow is assumed independent of y. The wind stress is specified in y-direction only: (0, 10-4).(x/500km 1) m2 s-2. Panel (A) is for the y-directed velocity and panel (B) the x-directed velocity. In (B), profiles of u are also plotted at the four indicated x-locations with scales 0.1 m/s given along the bottom for the first and fourth locations. Detailed plots comparing the profile at x=205 km (i.e. the first x-location) with analytical solution are given in Figure 1-10. The Coriolis parameter is constant, fo = 610-5 s-1, and the vertical eddy viscosity is also constant, K = 510-3 m2 s-1. This numerical solution is at time = 100 days, and differs slightly from the analytical one discussed in the text; instead of the no-slip condition at the bottom, a Page 40 matching of the velocity to the law-of-the-wall log layer is used [see the POM manual by Mellor, 2004]. Page 41 Figure 1-10. Profiles of (A) u, (B) v and (C) w at the x = 205 km location described in Figure 1-9b. Solid line indicates analytical solution assuming no-slip condition at the bottom, while dashed line the corresponding numerical (i.e. Figure 1-9) solution assuming log-layer at the bottom. Black line indicates the Ekman velocity while red (in panels (A) and (B)) indicates the total (Ekman + geostrophic) velocity. Note that in panel (B), total v/4 is plotted (so the profile fits in the same panel). Page 42 CHAPTER 2: Wind-Driven Homogeneous Ocean: Boundary Currents In this chapter we study the effects of wind on the general circulation of the ocean. We will examine how the large-scale oceanic gyre may be explained by the theories of Sverdrup [1947], Stommel [1948] and Munk [1950]. The circulation consists of a slow equatorward drift in the mid-ocean and over the eastern portion of the ocean basin, and a swift and narrow poleward current along the western boundary. We then numerically simulate these large-scale features with simple model experiments using POM. We will compare the numerical results against theories. We still assume that the fluid is homogeneous ( = constant), and that the motions have horizontal scales much larger than the vertical scales; in other words, we still assume that the hydrostatic equation (1-30c) is satisfied. In fact the scales are so large (~1000 km’s) that the Coriolis parameter “f” is no longer a constant, i.e. equation (1-37d) is not valid, and f = 2.sin(), /2 +/2 (2-1) can be a function of the latitude 2-1: The -plane & the QG-Potential Vorticity Equation For the moment, unprimed symbols are dimensional variables; primed symbols (such as u’ etc) denote dimensionless variables. Suppose we focus on a certain latitude, = o where f = fo, and let “y” be the northward (distance) coordinate measured from this latitude (so that y = yo at = o), then for “small” distances away from this latitude, we have: f fo + [df/dy][yyo] = fo + [(df/d)/ro][yyo] = fo + [2.cos(o)/ro][yyo] (2-2) or, f fo + o(yyo) (2-3a) where o = 2.cos(o)/ro (2-3b) Page 43 Equation (2-3) constitutes the so called “-plane,” and is valid for |yyo| not too large (i.e. a few thousands of kilometers, for o 2×10-11 m-1s-1 near mid-latitude.) We now derive a vorticity equation similar to equation (1-82) but taking into account that f is a function of y, i.e. is given by (2-3). In terms of dimensionless variables, denoted by primes, (2-3) is: f’ f’o + ’(y’y’o), ’ = (oL/|fo|) << 1. (2-4) where f’o = fo/|fo| = ±1 denotes northern or southern hemisphere. Note that ’ is small; we will in fact let ’ , and from (1-39): = ’/ = oL2/U O(1). (2-5) With (2-4), the Coriolis terms in equation (1-38) (which is dimensionless) becomes: f'v’ = [f’o + ’(y’y’o)][v’g + v’1 + O( 2)] = f’ov’g + f’ov’1 + ’(y’y’o) v’g + O( ’, 2 ) (2-6a) f'u’ = f’ou’g f’ou’1 ’(y’y’o) u’g + O( ’, 2 )] (2-6b) With these expressions, substituting (1-79) into (1-38b,c), and collecting and equating terms that are of O(1) and O( , ) separately we obtain for the O(1) terms the geostrophic relation as before, while the O( , ) terms give (c.f. equation 1-81): (2-7a) (2-7b) (2-7c) (2-7d) where fo = ±1, = ’/ , and dimensionless variables are now unprimed (except for the perturbation pressure p’1). Taking the curl, i.e. [(2-7c)/x(2-7b)/y] now gives (instead of 1-82): Page 44 (2-8a) where (2-8b) Since (ug, vg) hence g are independent of z, equation (2-8) can be integrated to yield, upon using (1-67a) and (1-78a) (c.f. 1-83): (2-9) In dimensional form [with asterisks *; using (1-39) and (1-55b)], this is: (2-10) where is the z-component of the curl of the kinematic wind stress (see equation 1-55a) and ; (2-11) is the bottom friction coefficient; l* is the Ekman depth (from 1-46); note that an asterisk has been placed on D, i.e. D* to avoid confusion with non-dimensional variables. 2-2: A Model of Large-Scale Ocean Gyre & Western Boundary Current In this section, unless otherwise stated, all variables will be dimensional. Consider a “box” ocean with sides Lx and Ly on a -plane (Figure 2-1). A zonal wind stress is applied: (2-12) so that Page 45 (2-13) where wo (which will be set = 10-4 m2s-2 below) is the magnitude of the (kinematic) wind stress. Initially, the fluid is at rest: (u, v, w) = (0, 0, 0), and also, free surface = 0 (2-14) At the four side walls, both the normal and tangential components of the horizontal velocities are zero: (u, v) = (0, 0) at x = 0, x = Lx and also at y = 0, y = Ly. (2-15) At the ocean bottom and at the free surface, equations (1-31) are satisfied: (2-16a) (2-16b) The continuity equation and the x and y momentum equations are then (c.f. 1-35): (2-17a) (2-17b) . (2-17c) where (2-18) is the depth-integrated transport vector (per unit length). These equations differ from (1-35) in that the nonlinear advective terms in momentum equations (2-17b,c) have been dropped. Page 46 Figure 2-1. A “box” ocean circulation driven by eastward (northern half of basin) and westward (southern half) wind stress distribution. Depth-integrate (2-17b,c), and again drop nonlinear terms such as u etc: (2-19a) (2-19b) where w and b are wind and bottom stress vectors respectively. These equations together with equation (2-17a) constitute there equation for three unknowns (, U, V). For our purpose, it is more convenient to write them in terms of the depth-averaged velocity: (2-20) so that (2-17a) and (2-19) become: Page 47 (2-21a) (2-21b) (2-21c) where r is the bottom friction coefficient with dimension m/s, we have modeled the bottom friction as: ( , (2-22) and we have also added a source term Q to the right hand side of (2-21a). This we will use later to model the net effect of precipitation ( ) minus evaporation ( ), i.e. Q=( – )/o (unit is m/s). (2-23) Let H = constant, a vorticity equation similar to (2-10) may now be derived from (2-21). Taking the curl of (2-21b,c), we obtain: (2-24a) Apart from the terms involving Q and /t, this equation is equivalent to (2-10). The term /t does not appear in (2-10) because it is neglected when we derived the vertical velocity at the surface (see equation 1-67). We will explain in the next section when this term may be deleted; assuming that it can be, then (2-24a) becomes: (2-24b) The appearance of -Qf together with curlz(w) means that as far as the large-scale ocean circulation is concerned, effects of wind stress curl and net precipitation minus evaporation are equivalent. The reason is that in both cases, the ocean layer beneath the surface boundary layer (where wind and precipitation etc act) ‘sees’ only vertical mass flux coming in (or out) of the surface layer. In the case of wind stress curl, this is “Ekman pumping” (Chapter 1) caused by convergence or divergence in the surface layer, Page 48 while in the case of precipitation etc, it is more direct – due to the flux of water from the atmosphere. We will model these effects and give a physical interpretation below. Also, while H is the ocean depth, it may be (and for our purpose will be) interpreted as the depth of the main thermocline in which the predominant large-scale currents reside. In that case, the friction is then considered as the friction at the base of the main thermocline. Sverdrup Relation: In the open ocean away from side boundaries, the relative vorticity is assumed weak, in fact , and the main vorticity balance in steady state is, from (2-24b) (omitting Q for simplicity): (2-25) Figure 2-2. A schematic to explain the equatorward Sverdrup transport V caused by a negative wind stress curl, curlz(w) < 0. Page 49 Equation (2-25) is the Sverdrup relation; it shows that in steady motion, a negative curlz(w) produces southward transport. Figure 2-2 explains the physics involved. The negative curlz(w) produces convergent flow in the thin surface (Ekman) layer (Chapter 1) and hence a downward pumping that squashes the vortex tube in the deeper interior; this is equivalent to a divergent subsurface flow. Squashing increases the cross-sectional area of the tube, and if it stays at the same latitude a negative relative vorticity (i.e. a reduced spin) would be produced, . To keep the same 0, and assuming that the tube remains approximately vertical (and it does), the tube must move south along the curved earth surface where it can stretch and retain its initial zero spin. The curved surface in this case is the physical equivalence of the -effect, i.e. the Coriolis parameter f varies with latitude. Gill [1982; p.465] gives another, equivalent explanation. The divergent subsurface flow beneath the negative curlz(w) increases the cross-sectional area of the vortex tube, so the absolute vorticity (i.e. f + ) decreases. But since f >> , the only way that this decrease can be accomplished is that the tube moves southward where f is smaller. Figure 2-2 illustrates two other features. One is that the effect of (positive) Q [= ( – )/o] is similar to (negative) curlz(w): this is because Q ‘pumps’ water into the subsurface layer and again squashes vortex tube. The other one is that while the surface is slightly elevated ( > 0) it is not important in the argument given above; in steady state, it is in fact = 0. If as shown in Figure 2-1 the Sverdrup transport occurs over most portion of the ocean basin, then the total southward transport is: SV (m3 s-1). (2-26) Page 50 Figure 2-3. Annual mean distribution of Ekman pumping, positive = upwelling, in unit of 10-6 m s-1 (or 0.1 m/day). From Tomczak and Godfrey. Figure 2-3 shows a global map of Ekman pumping – note large areas of downwelling in the subtropical regions, suggesting equatorward Sverdrup transports there. From (2-26), we have roughly SV curlz(w).Lx/o 2 Lx/(Lyo) (2×10-4/10-11)(Lx/Ly) m3 s-1 20 Sv (1 Sv = 106 m3 s-1) for Lx Ly, 10-4 m2 s-2 and o 10-11 m-1s-1. The corresponding Ekman pumping (downwelling) is from (1-64b) = curlz(w)/fo 2 /(Lyfo) 1~2×10-6 m s-1 for Ly 2000 km and fo 6×10-5 s-1, roughly equal to the values shown in Figure 2-3. This Ekman velocity (or Figure 2-3) may be compared with contours of Q [= ( – )/o] shown in Figure 2-4, which gives local maxima of roughly 10-7 m s-1, about 10 times weaker than the Ekman pumping. Page 51 Figure 2-4. Annual mean distribution of precipitation minus evaporation, Q = ( – )/o, in m yr-1 ( 3×10-8 m s-1); shaded regions are positive. Stommel’s Model of the Western Boundary Current (WBC): In our idealized closed basin (Figure 2-1), the southward Sverdrup transport must be returned northward. Since this return cannot be in the mid-ocean (i.e. open ocean) because of the Sverdrup constraint (as seen above), it must occur near the western and/or eastern boundaries where the Sverdrup balance is broken. In other words, other terms in the potential vorticity equation (2-24b), such as friction and/or nonlinearity (which is neglected in 2-24b), become important near the boundary. Stommel assumes that (bottom) friction is important so that (2-24b) becomes, in steady-state: . (2-27) Including curlz(w) adds unnecessary complexity without adding physical insight so it is omitted. Also, near the western or eastern boundary, . The solution is: (2-28) where C(y) is an arbitrary function of y to be determined. Equation (2-28) shows a northward jet which is maximum at the boundary (for an eastern boundary, the “x” is Page 52 replaced by “xLx”) and either decays (if western boundary) or grow (if eastern boundary) exponentially away from the boundary. Therefore, the northward return flow must be along the western boundary and crucially depends on a non-zero o; the width of this return jet may be estimated as the e-1 x-scale: xwbc = r/(Ho) (2-29) This is 50 km for r 10-4 m s-1, H 200 m and o 10-11 m-1 s-1. The C(y) is determined by requiring that the total northward transport carried by the WBC is equal to the southward Sverdrup transport at the same latitude (i.e. same “y”). Thus, , i.e. (2-30) With the aforementioned values of SV 20 Sv, o 10-11 m-1 s-1 and r 10-4 m s-1, C 2 m s-1. The Stommel WBC solution (2-28) therefore consists of a jet concentrated against the western boundary. Since this jet has to return northward all the Sverdrup transport across a narrow width xwbc, the jet’s velocity is large. The Sverdrup velocity scale is, from (2-25), curlz(w)/(Ho) 2 /(HLyo) (2×10-4/10-11)/(HLy) m s-1 0.1 m s-1 for H = 200 m and Ly = 2000 km; while the Stommel velocity scale is C 2 m s-1. Stommel’s jet is maximum at the coast x = 0 and decays exponentially offshore to become very weak in the ocean’s interior where the southward Sverdrup transport prevails. Note that the Stommel’s jet, being wholly positive, cannot join smoothly to the Sverdrup flow, nor is it expected to since it is based on a frictional equation that ignores the Ekman pumping by the wind; but this is not important. Nor is it important that the Stommel’s jet cannot satisfy the no-slip condition at the coast. The important thing is that: (1) beta (o) is crucial for the existence of the western-intensified boundary current, and (2) the WBC jet provides a mechanism that closes the ocean-basin mass balance. Munk’s Model of the Western Boundary Current (WBC): Page 53 A slightly more complicated model can be constructed to satisfy the no-lip condition at x = 0. Instead of the bottom friction we include horizontal viscosity terms on the RHS of the momentum equations: (2-31a) (2-31b) (2-31c) where AH is the horizontal (eddy) viscosity coefficient (unit is m2 s-1), assumed a constant. If we now take the curl etc as in the derivation of (2-24b) we obtain: (2-32) Define a stream function (unit m2 s-1) such that: (2-33) and therefore (2-34) In the WBC, we now let the beta term balances the horizontal viscosity term, equation (2- 32) in steady state is then (c.f. 2-27): (2-35a) where = (o/AH)1/3x. (2-35b) Page 54 Equation (2-35a) is a fourth-order ordinary differential equation (with “y” as the ‘floating’ parameter not directly involved), and a solution of the form exp(m ) maybe assumed. Four roots are: m1 = 0, m2 = 1, m3 = exp(i/3)=(1+i3)/2, m4 = exp(i2/3) =(1i3)/2 (2-36) The solution then is of the form = C1 + C2exp(m2 ) + C3exp(m3 ) + C4exp(m4 ) (2-37) where the C’s are arbitrary function of “y” yet to be determined. The m2 root must be deleted - by choosing C2 = 0, because it leads to a solution which exponentially grows with (or x). The root m1 leads to C1 – which is chosen to be C1 = I(y), where the subscript “I” denotes “interior” (of the ocean away from the boundary) and from the interior Sverdrup solution (2-25) with , (2-38) is the stream function of the interior’s Sverdrup flow; o(y) is an arbitrary function. This choice is so that ~ I for ~ , since m3 and m4 give exponentially decaying solutions for large . Substituting m3 and m4 from (2-36) to (2-37), we obtain then: = I + [C3 + C4 ] (2-39) the C3 and C4 here differ from those in (2-37) but the detail is not important; they are still arbitrary functions of “y.” The corresponding velocity components from (2-33) are: (2-40a) (2-40b) The conditions that there can be no normal flow ( ) and no slip ( ) at = 0 give: Page 55 C3 = I(0,y) + K, where K = constant, and (2-41a) . (2-41b) Therefore, (2-39) gives: = I(x, y) I(0,y) [ + ] +K [ + ] (2-42) From equation (2-35b), we note that (AH/o)1/3 has the unit of length (m). Its value is 100 km (for AH = 103 m2 s-1 and o 10 -11 m-1 s-1) which is thin compared to the width of the basin 2000 km (say); it is the e-folding decay scale of the western boundary current near x 0. We are therefore making only a small error if I(0,y) in (2-42) is replaced by I(x,y). Another way of saying the same thing is that, since I varies over a large scale >> (AH/o)1/3, I(x, y) I(0,y) near the western boundary. To determine K, we note that, since our basin is closed, the stream function around the boundary is zero. Therefore, for to be also = 0 at x = 0, (2-42) shows that K must also = 0. Therefore, (2-42) becomes: = I(x, y) {1 [ + ]} (2-43) The northward velocity is /x: , = (o/AH)1/3x, (2-44) in which we have neglected the interior velocity contribution I/x, which is much smaller. Note that the jet form is again apparent by the exponential decay term. But in contrast to the Stommel’s solution equation (2-28), because of the sine function, the northward velocity is zero at x = 0. Also because of the sine function the velocity becomes slightly negative (i.e. reversed flow) at the offshore edge of the jet. Class: estimate the magnitude of from (2-44). Page 56 Dear Jane, Eda and Alu, here (section 2-3) is a derivation of the quasi-geostrophic potential vorticity equation used in Sheremet’s paper (his eqn.1). He used a so called reduced-gravity form of the equation. Pls. read and work out the steps as much as you can, and hopefully I can help you on Thursday. Have a good weekend. leo 2-3: The Equivalent Reduced-Gravity Model In this section non-dimensional variables will be denoted by subscripts “n.” The above model has one layer and involves the free surface . We now show that with slight reinterpretations of the variables, the equations of that model are exactly the same as those governing the motion of a two-layer ocean in which one of the layers is infinitely thick and is quiescent. Before we do this, we find out first the condition when the term involving /t in equation (2-24a) may be neglected. This is necessary in order to understand the restrictions on the spatial and temporal scales under which the resulting equation(s) apply. Instead of (2-24a), we will non-dimensionalize the vorticity equation that uses the horizontal viscosity instead of the bottom friction, i.e. we use (2-32) with f(/t)/H retained: (2-45) Define the following non-dimensional variables: w = wno, Q = Qn Qo, AH = AHn AHo. (2-46) Note that the scale for , i.e. fUL/g, is determined from balancing the pressure gradient and Coriolis terms in equation (2-31b,c). Substitute (2-46) into (2-45): (2-47) where R = (gH)1/2/f is the Rossby radius of deformation. We see that dropping the term involving n/tn amounts to saying that R >> L, i.e. the spatial scale of motion is much smaller than the Rossby radius. For H 1000 m and f 6×10-5 s-1, R 1500 km, Page 57 so that motions with scales smaller than about 1000 km will not be much affected by the production of vorticity due to stretching associated with the motion of the free-surface. This kind of motion is called a barotropic motion. On the other hand, as we will now show, R can be quite small for the reduced-gravity model, in which case the stretching term becomes important. The Reduced-Gravity Model from the Two-Layer Model: We now derive the reduced-gravity (or 1.5-layer) equations from a two-layer model as sketched schematically in Figure 2-5: Figure 2-5. An illustration in the xz-plane of an ocean model with two (inmiscible fluid) layers with subscripts “1” and “2” denoting the upper and lower layers respectively. Symbols are: 1 and 2 = constant densities of the upper and lower layers respectively, 2 > 1; H1 = constant undisturbed upper-layer depth; H2(x,y) = undisturbed lower-layer depth which may be a function of (x,y) for variable bottom topography; h(x,y,t) = layer depth; u(x,y,t) = the x-component velocity; v(x,y,t) = the y-component velocity; p(x,y,z,t) = pressure; 1(x,y,t) = disturbed free-surface height measured from z = 0; 2(x,y,t) = disturbed inter-layer height measured from z = -H1. Note that 1 and 2 can be positive (upward) or negative (downward), and that h1 = H1 + 1 2. Page 58 The caption of Figure 2-5 explains all the variables to be used. Since in each layer the density is constant, the continuity (conservation of volume) equation takes exactly the same form as that for a homogeneous fluid (see equation 2-31a): , i = 1, 2. (2-48) where i = 1 is for layer-1 and i = 2 is for layer-2, and for simplicity the source/sink term Q is omitted. (The subscripts ‘i’ in equation 2-48 are not to be confused with x and y components, and repeated subscripts do not mean summation!). Except for the pressure gradient terms, the momentum equations (2-31b,c) will also remain essentially the same. To evaluate the pressure gradients, we will assume as before that the horizontal scales are much larger than the vertical scales, so that the pressure is hydrostatic (see equations 1-1 and 1-2). Thus the pressure at some depth z = z’ in the water is just the weight per unit area of fluid above. For a point in layer 1, the pressure p1(x,y,z,t) is: p1(x,y,z,t) = pa(x,y,t) + g1(1 z), for 1 z H1+2, (2-49a) where pa = atmospheric pressure (i.e. = weight per unit area of air above z = ). Similarly, for a point in layer 2, the pressure p2(x,y,z,t) is: p2(x,y,z,t) = p1(at z =H1+2) + g2(H1 + 2 z), for H1+2 z (H1+H2), (2-49b) Then the horizontal pressure gradients are (for example /x): p1/x = pa/x + g11/x, for 1 z H1+2, (2-50a) p2/x = p1/x + g(1h1)/x, for H1+2 z (H1+H2), (2-50b) where = (2 1), and we have used (2-49a) to evaluate p1(at z =H1+2), and also 2 = H1 + 1 h1 in obtaining the last term in (2-50b) (see Figure 2-5). The momentum equations (2-31b,c) applied to layers 1 and 2 separately lead to (2-51b,c) and (2-52b,c) below: (2-51a) Page 59 (2-51b) (2-51c) (2-52a) (2-52b) (2-52c) Note that we have generalized /t to D/Dt in the momentum equations. These are the equations for our two-layer ocean. More complicated forms with additional terms exist in the literature; for our purpose, the above suffices. Note that the lower layer experiences both the layer-1 pressure gradient (p1/x)/o as well as the additional force (per unit mass) (g/o)(1h1)/x; and similarly for the y-component. In the reduced-gravity model, we let one of the layers to be infinitely thick. Supposing we let layer-2 to be infinitely thick, then from the continuity equation (2-52a) the deep- layer currents (u2, v2) must approach zero, and h2/t ~ 0 also. Subtracting (2-52) from (2-51) then gives: (2-53a) (2-53b) (2-53c) in which we have put back the “Q” on the RHS of the continuity equation, and have made a further approximation that 1x << h1/x and similarly 1y << h1/y, which means that the variations of the free surface is negligible in comparison to the variations of the upper-layer depth. Equations (2-53) are identical to the equations that govern the motion of a single-layer (homogeneous) fluid of constant depth H1 (i.e. the topography) with a free surface except that the acceleration due to gravity g is reduced to: g’ = g/o (2-54) Page 60 which appropriately is called the “reduced gravity.” To see the equivalence, set h1 = H1 + , (2-55) where now the “” is interpreted as the departure of the upper-layer depth from its ‘mean’ of H1. For typical /o 10-4~10-2, the corresponding shallow-water baroclinic phase speed Ci = (g’H1)1/2, hence also the baroclinic Rossby radius of deformation Ri = (g’H)1/2/f, is some 10 to 100 times smaller. Therefore, for reduced- gravity (baroclinic) motions, the time-dependent stretching term (L/Ri)2 h1n/tn that produces vorticity in equation (2-47) can be significant. In fact, this term overwhelms the other time-dependent term 1n/t in the limit of (L/Ri)2 >> 1. (2-56) The reduced-gravity equivalence of (2-45) is: (2-57) in which the nonlinear advection term u1. 1 (from the D/Dt terms in 2-53b,c) is included. But, consistent with the quesi-geostrophic approximations of Chapter 1, we may replace the u1 in this nonlinear term by its geostrophic values: H1v1 = +(Ci2/f)h1/x, H1u1 = (Ci2/f)h1/y (2-58) Equation (2-58) suggests a stream function 1 such that: 1/x = H1v1 and 1/y = H1u1 (2-59) By comparing (2-58) and (2-59), we then have: 1 = (Ci2/f)h1 = fRi2 h1 (2-60) Equation (2-57) then becomes: Page 61 (2-61) where J(A,B) = A/x.B/yB/x.A/y is the Jacobian (of the functions A and B). Eda, Jane and Alu, this equation (2-61) is identical to Sheremet’s equation (1), except that his is actually my 1/H1, and he sets the wind stress curl and Q to zero. 2-4: An Application of the Reduced-Gravity Model: Jet & Eddy Interaction Figure 2.6a shows an eddy, either cyclonic or anticyclonic, drifting westward and approaching a straight meridional jet. The jet’s velocity is a maximum along its center, and exponentially decays symmetrically on either side. Such a velocity profile is obtained by specifying a constant potential vorticity (PV) jump across the jet. The PV is defined as an anomaly from a rest state. For a reduced gravity model an upper layer of depth h = H + (where H = undisturbed layer depth and = perturbed layer height) and velocity (u, v) overlying an infinitely deep resting lower layer the PV is: PV = (f + )/h fo/H = ( fo/H + y)/(H + ), for a -plane, (2-62a) or, PV = ( fo/H)/(H + ), for an f-plane. (2-62b) Then the PV-distribution of the jet is: PVJ = +QJ/2, x > 0, = QJ/2, x < 0. (2-63) Note PV-jump|Jet = QJ. Similarly, for the eddy, assumed to be of a circular shape with radius REddy, we have, PVEddy = QE, r < REddy, = 0, r > REddy, (2-64) Page 62 and the PV-jump|Eddy = QE. Here, r denotes the radial distance from the center of the eddy. Inverting the PVJ and PVEddy to Obtain the Corresponding Velocity Fields: Figure 2-6a. A sketch of the jet with PV-jump|Jet = QJ, and the approaching eddy with PV-jump|Eddy = QE. Figure 2-6b. A sketch of a weak (cyclonic) eddy and/or a strong jet, QJ > QE. The eddy does not cross the jet and is being advected along with the jet. For a weak eddy and/or a strong jet, such that QJ > QE (Figure 2.6b), the eddy introduces only a small perturbation on the jet. The eddy does not cross the jet and is being advected along with the jet [Stern and Flierl, 1987; Bell, 1990]. Page 63 Figure 2-6c. A sketch of the evolution at four indicated times for the combined system of a strong (cyclonic) eddy and/or a weak jet, QE > QJ. If the eddy’s strength is increased and/or the jet’s strength is weakened, QE > QJ, where > 1, the eddy interacts nonlinearly with the jet. In other words, the eddy’s velocity field can significantly influence the jet, and the jet’s velocity field in turn affects the eddy. The evolution of the combined system of a cyclone and jet is illustrated in Figure 2.6c for four time instances, t1 < t2 < t3 < t4. Time t1 is when the jet begins to “feel” the presence of the eddy “E,” when the latter comes within 2Ro of the jet’s axis, where Ro = (g’H)1/2/fo is the Rossby radius. The sketch at time = t2 shows that the eddy’s velocity field downstream (upstream) pushes the jet westward (eastward) to develop the downstream (upstream) meander “DM” (“UM”). Vorticity inside “DM” (“UM”) is anticyclonic (cyclonic). The “dipole” formed by DM and E moves southwestward due to their mutually-induced velocity shown by the green vector, while UM and E also have a mutually-induced velocity shown by the blue vector. The resultant vector is westward indicated by the red vector in a direction that tends to cross the jet. Figure 2-6c shows how this jet-crossing may occur at later time = t3, during which the E and UM merge, but the DM-E dipole continues to “peal” itself from the jet. This “pealing” process is complete at time = t4 when the original eddy E has completely crossed the jet. Page 64 If eddy “E” is an anticyclone, then the same argument as Figure 2-6c will show that the eddy cannot cross the jet because the resultant (red) vector would be eastward. However, the jet may become baroclinically unstable (this is not modeled by the reduced- gravity model because baroclinicity requires an active lower layer), and jet-crossing may still occur. References Bell, G.I., 1990: Interaction between vortices and waves in a simple model of geophysical flows. Phys. Fluids A 2 (4), 575-586. Stern, M.E. and G.R. Flierl, 1987: On the interaction of a vortex with a shear flow. J. Geophys. Res, 92 (C10), 10733-10744. Eda, start here please. Thank you. Page 65 Pressure Compensation In Smith and O’Brien [1983; JPO, 13, 1681-1697], you will hear many times the word “compensation.” All this means is that the lower-layer has very small velocity. In my note “pom_v14.pdf” Fig.2-5, this just means that the 1 and 2 are of opposite signs. To see this, we use lower-layer equation (2-52b,c) and set the pressure gradient terms on the RHS to zero. After integration, and note that there is an arbitrary constant “C” (has to be a constant – why?), we get: 1 + C = ( or 1 = ( choosing C = (. When the free surface convex upwards, 1 > 0, the isopycnal separating layers 1 and 2 concave downward, 2 < 0. 2-5: An Application of the Two-Layer Model Isobe [2004; JPO, 34, 1839-1855; equations will be prefixed with an”I”] provides a nice application of the 2-layer model. The 2-layer system was previously shown in Fig.2-5 (of the “GFD book”), copied here: Figure 2-5. An illustration in the xz-plane of an ocean model with two (inmiscible fluid) layers with subscripts “1” and “2” denoting the upper and lower layers respectively. Symbols are: 1 and 2 = constant densities of the upper and lower layers respectively, 2 > 1; H1 = constant undisturbed upper-layer depth; H2(x,y) = undisturbed lower-layer depth which may be a function of (x,y) for variable bottom Page 66 topography; h(x,y,t) = layer depth; u(x,y,t) = the x-component velocity; v(x,y,t) = the y-component velocity; p(x,y,z,t) = pressure; 1(x,y,t) = disturbed free-surface height measured from z = 0; 2(x,y,t) = disturbed inter-layer height measured from z = -H1. Note that 1 and 2 can be positive (upward) or negative (downward), and that h1 = H1 + 1 2. The equations of a two-layer system were derived in section 2-3, equations (2-51) and (2- 52). The x and y-momentum equations (2-51b,c) for layer-1, ignoring wind stresses, atmospheric pressure and horizontal viscosity, are: (2-5.1a) (2-5.1b) Taking the curl: (2-5.1c)/x - (2-5.1b)/y] then gives: 1/t + ( 1u1)/x + ( 1v1)/y + f .u1 = 0 (2-5.2) For layer-2 (see equation (2-52)), the x and y-momentum equations are: (2-5.3a) (2-5.3b) where g’ = g/o is the reduced gravity. Taking the curl, we obtain exactly the same equation (2-5.2) with subscript “1” replaced by “2.” Thus, we write: k/t + .( kuk) + f .uk = 0 (I2.1) (2-5.4) The continuity equations are from section 2-3, equations (2-51a) and (2-52a): hk/t + .(hk uk) = 0 (I2.2) (2-5.5) Split each field variable into time-mean and fluctuation: Page 67 u = <u> + u’, h = <h> + h’, = < > + ’ (2-5.6) where <.> denotes a mean and primes denote fluctuations such that <A’> = 0 for any variable A. Substitute (2-5.6) into the continuity equation (2-5.5), then take time- averaging: .(<hk> <uk>) = .(<hk‘uk’>) (I2.4) (2-5.7) The LHS represents divergence of the mean mass flux field, and this is balanced by the divergence of the mass flux due to eddies – the RHS. Analogous to heat diffusion by turbulent mixing, say -<w’T’> which we model as: -<w’T’> = (KHT/z)/z (2-5.8) see Fig.2-5.1, Fig.2-5.1 On average, turbulent velocities bring warmer water downward and cooler water upward, as shown in the sketch. Therefore, -<w‘T’> > 0. This is often called down-gradient diffusion (down the gradient of T – i.e. heat goes from warm to cold – and is commonly often modeled as <w‘T’> = KH dT/dz, where KH is the turbulent eddy diffusivity. we also model the RHS of (2-5.7) as: -<hk‘uk’> = k <hk > (I2.5) (2-5.9) where k is the diffusivity for layer-thickness (which you can think of also as “heat” since thicker h means warmer water). Equation (2-5.7) then becomes: Page 68 .(<hk > <uk>) = .(k <hk >) (I2.6) (2-5.10) By using the idea of down-gradient diffusion, i.e. using “diffusivity,” we avoid the very difficult task of computing turbulence. In the case of large-scale geophysical fluid flows (e.g. Kuroshio meanders and frontal eddies), we then avoid directly analyzing ‘eddies’ – which are called “geostrophic turbulence.” Recall in GFD course that eddies in geostrophic turbulence are rigidly constrained by rotation because of the Taylor- Proudman theorem, so that their motions are “planar” in the horizontal plane perpendicular to the rotation axis. In the case of Kuroshio frontal eddies, equation (2- 5.9) models in a crude way how eddies transport heat from the Kuroshio to the (cooler) shelf, and how they also transport cooler shelf water to the Kuroshio – i.e. mixing. The development of these eddies are say, via baroclinic instability. In Isobe (2004), these eddies (or some idealized form of eddies) are actually numerically computed. Similar to time-averaging the continuity, the time-average of the vorticity equation (2- 5.4) gives: .(< k> <uk>) + .(< k‘uk’>) + f .<uk> = 0 (I2.7) (2-5.11) (note that the time-average of the time-dependent term k/t is 0). The second term .(< k‘uk’>) can again be modeled using some kind of eddy viscosity times the gradient of the mean vorticity. But it is usually small, and is neglected. Eddy-mixing of relative vorticity in (2-5.11) is not as important as eddy-mixing of heat in equation (2- 5.10) because | k’| is generally smaller than |f| (and also smaller than |< k>|). Physically, what is important is the mixing of heat across the Kuroshio, and less so for the relative vorticity. Mixing of heat (h’) causes layer to change which intuitively is also important because the change directly affects potential vorticity: e.g. (f + 1)/h1 f/h1. The f .<uk> in (2-5.11) can be obtained from (2-5.10): f .<uk> = f<uk> <hk>/<hk> + f .(k <hk >)/<hk> = +<uk><hk>. (f /<hk>) + f .(k <hk>)/<hk>;(assume f –plane) = +<Uk>. (f /<hk>) + f .(k <hk>)/<hk>;(let <Uk> = <uk><hk>) Page 69 Substituting this into (2-5.11), we obtain: .(< k><uk>) + <Uk>. (f /<hk>) = f .(k <hk>)/<hk> (I2.8) (2-5.12) Now the down-gradient portion on the RHS of the previous equation, “ <hk>”, can be approximated by “thermal wind” for this 2-layer system which can be obtained from the geostrophic form (time-averaged) of (2-5.1) and (2-5.3). We first do (2-5.1) minus (2- 5.3) (i.e. upper-layer minus lower-layer), in vector form: z×<ug1-2>=(g’/f) (<h1>-<1>)(g’/f) <h1> (since |1| << h1 see fig.2-5) so that <h1> = z×<ug1-2>(f /g’) (2.9) (2-5.13) Substituting this into (2-5.12) for the upper layer (k = 1), and noting that: .( z×<ug1-2>)=<ug1-2>. zz. ×<ug1-2> (Use .(a×b)=b.( ×a)-a.( ×b)) = 0 Curlz(<ug1-2>)(Since z = 0) we obtain from (2-5.13) and (2-5.12): .(< 1><u1>) + <U1>. (f /<h1>) = f21g’-1 Curlz(<ug1-2>)/<h1> (I2.10)(2-5.14) To derive the lower-layer vorticity equation, we similarly get the “thermal-wind” equation by doing (2-5.3) minus (2-5.1) (i.e. lower-layer minus upper-layer), and noting from Fig.2-5 that (1-h1) = 2 = (h2-H2), where H2(x,y) is the bottom topography: z×<ug2-1>=(g’/f) ( <1>-<h1>)=(g’/f) (<h2>H2) (2-5.15) from which we obtain: <h2> = z×<ug2-1>(f /g’) + H2 (2-5.16) Substituting this into (2-5.12) for the upper layer (k = 2), we obtain: .(< 2><u2>)+<U2>. (f/<h2>)=f22g’-1Curlz(<ug2-1>)/<h2>f .(2 H2)/<h2> Page 70 (2-5.17) However, if we assume in (2-5.15) that <h2> is >> H2, i.e. that the topographic slope varies slower than the interfacial slope, then: .(< 2><u2>) + <U2>. (f /<h2>) = f22g’-1 Curlz(<ug2-1>)/<h2> (I2.11)(2-5.18) 2-6: Heuristic derivation of the quasigeostrophic (QG) potential vorticity (PV) equation for a stratified fluid The inviscid x, y & z-momentum and continuity equations are: u/t + uu/x + vu/y + wu/z -fv = -r-1(p/x) (2-6.1a) v/t + uv/x + vv/y + wv/z +fu = -r-1(p/y) (2-6.1b) p/z = -g (2-6.1c) u/x + v/y + w/z = 0 (2-6.1d) Also, a beta-plane is assumed such that: f = fo + oy (2-6.2) where fo is the Coriolis parameter at y=0 (around where later we will develop the instability analysis), o = df/dy, and it is also assumed that |oy| << |fo|. For o 10-11 m-1 s-1, and |fo| 5×10-5 s-1 or larger (poleward of 20oN/S), we restrict |y| 1000 km or less. The symbols have the usual meanings, and r = reference density which is a function of z only. The basic QG-approximation idea is that the motions are close to being geostrophic, so we will use the geostrophic velocity assuming constant f fo (which is ‘zeroth-order’) to evaluate the “difficult” terms which in (2-6.1a,b) are the non-linear and “o” terms. To the zeroth-order approximation, the motions are nearly geostrophic, so that: -fovo = -r-1(po/x) +fouo = -r-1(po/y) (2-6.3a,b) Page 71 so that uo/x + vo/y = -wo/z 0 (2-6.4) Therefore, wo 0 since it is 0 at the surface. In any case, the zeroth-order or geostrophic approximation of w = wo is small compared to |uo| or |vo|.6 Substitute the (uo, vo, wo0) into the non-linear and beta parts of (2-6.1), we then obtain the following equations for the next-order term with subscripts “1” (u1, v1, p1): do/dt{uo} – fov1 oyvo = -r-1(p1/x) (2-6.5a) do/dt{vo} + fou1 + oyuo = -r-1(p1/y) (2-6.5b) where do/dt = /t + uo/x + vo/y. (2-6.6) Taking the (z-component of the) curl of (2-6.5) to eliminate the pressure, and using the first of (2-6.4) that uo/x + vo/y 0, we obtain: do/dt{o} + ovo + fo H.u1 = 0, where H=i/x+j/y and u1=(u1,v1); or, do/dt{o + oy} = fow1/z (2-6.7) after using (2-6.6) and (2-6.1d), i.e. u1/x + v1/y + w1/z = 0. For = constant, by vertically integrating (from z=-H to z=0) equation (2-6.7) we recover the QGPV equation for a homogeneous fluid, which we studied previously in section 2-1 (see equation (2.10)). For stratified fluid, we need to evaluate w1/z using the density (or buoyancy b = -g/r) equation [see section 11.1.1 of Marshall and Plumb, 2008] which in the absence of heat/salt sources and sinks on the RHS is: /t + u/x + v/y + w/z = 0 (2-6.8) 6 Typically, w = Ekman pumping near the surface ×o/(fo) 0.1(N/m2)/[104(m).10-4(s-1).103(kg/m3)] 10-4 m/s 10 m/day, for a strong wind stress curl of 0.1 N/m2 over a distance of 10 km. Thus comparing to typical |uo| 5×10-2 m/s, the “w” is indeed small. In the open ocean, the “w” is typically smaller than 10 m/day. Page 72 The density is given by: = r(z) + ’(x,y,z,t) (2-6.9) and the hydrostatic equation is then assumed also for the perturbation density ’ (and pressure po): po/z = -g’ (also of course pr/z = -gr) (2-6.10a) or (po/z)/r = -g’/r = b’ (2-6.10b) Using the QG-approximation on (2-6.8), we get: do/dt{’} + w1(dr/dz) = 0, (2-6.11a) or multiplying by –g/r, we get do/dt{-g’/r} + w1(-gdr/dz)/r = 0, i.e. do/dt{b’} + w1N2(z) = 0 (2-6.11b) where N2(z) = (-gdr/dz)/r is the squared buoyancy (or Brunt-Vaisala) frequency. We see that (in the absence of sources and sinks) the vertical velocity w1 is related to the isopycnal movement due to the geostrophic flow: w1 = do/dt{b’/N2(z)}, or also = do/dt{’}/{-dr/dz}. (2-6.11c) The last form shows, since -dr/dz > 0, that deepening isopycnal (i.e. do{’}/dt < 0) corresponds to downwelling, and shallowing isopycnal corresponds to upwelling. We now use w1 in (2-6.7): do/dt{o + oy} = fow1/z = fodo/dt{[b’/N2(z)]/z} or upon using the hydrostatic equation (2-6.10b): do/dt{o + oy} = fodo/dt{[(po/z)/(rN2)]/z} i.e., do/dt{o + oy + fo[(po/z)/(rN2)]/z} = 0 (2-6.12a) Page 73 or, do/dt{( 2 po)/(for) + oy + fo[(po/z)/(rN2)]/z} = 0 (2-6.12b) since from (2-6.3a,b), we have for the geostrophic vorticity: o = vo/xuo/y = 2 po/(for) (2-6.13) For the ocean, r can be assumed to be constant in (2-6.12b), so that equation becomes: do/dt{ 2 (po/for) + oy + [(po/for)/z)(fo /N2)]/z} = 0 (2-6.14) From the geostrophic relation (2-6.3), the po/(for) is equivalent to geostrophic stream function: = po/(fo r) (2-6.15) so that (2-6.14) is usually expressed in terms of : do/dt{ 2 + oy + [(/z)(fo /N2)]/z} = 0 (2-6.16a) The quantity inside the {..} is the QGPV for a stratified fluid: Q= 2 + oy + [(/z)(fo /N2)]/z (2-6.16b) Non-dimensionalization: We use the following scales: Table 2-6.1 Scales Variables So that… (primes denote nondimensional variables) L (x,y) (x,y) = L (x’,y’) D z z = Dz’ U (u,v) (u,v) = U(u’,v’) L/U t t = (L/U)t’ fo rUL po po = fo rULpo’, using geostrophic eqn.(2-6.3) UL = UL’, using eqn. (2-6.15) Ns N N = NsN’, for some typical Ns value of N Then the non-dimensionalized form of (2-6.16) is (dropping all primes): Page 74 do/dt{ 2 + y + [S-1(/z)]/z} = 0 (2-6.17a) and Q= 2 + y + [S-1(/z)]/z (2-6.17b) where = oL2/U, S = (LD/L)2, and LD = NsD/fo is the baroclinic Rossby radius based on ocean’s depth D. Note that Ns2D ~ g/r ~ g’, the reduced gravity, so that Ns2D ~ g’D ~ baroclinic phase speed squared. The nondimensional = o/[(U/L)/L] = o/(o/L) is a measure of the ratio of gradient of planetary vorticity “o” to the gradient of flow vorticity “o/L”. By assuming that ~ O(1) in (2-6.17), we are examining an important GFD case in which the planetary vorticity gradient contributes equally with the relative vorticity gradient to the overall vorticity balance. If the scale of motion is small compared to the Rossby radius, L << LD so that S >> 1, the effect of vortex-tube stretching to the PV-balance (i.e. the [S-1(/z)]/z term) is small, and the flow’s relative vorticity (i.e. 2) becomes important. If on the other hand the scale of motion is large, L >> LD, then S << 1, and the flow’s vorticity becomes relatively unimportant, i.e. the flow looks horizontally more uniform. 2-7: Baroclinic instability Fig.2-7.1 suggests that baroclinic instability may occur along special paths “(iii)” across slanting isopycnals which of course imply the existence of vertical velocity shear by the thermal-wind equations (from geostrophic eqn.(2-6.3a,b) and hydrostatic eqn.(2-6.10b)): k × fo uo/z = b’ (2-7.1a) or fo uo/z = k × b’ (2-7.1b) We will now find the condition – what type of vertical shears can produce baroclinic instability? Or, in general, what type of vertical and horizontal shears can produce baroclinic and barotropic instabilities? We will use the non-dimensionalized QGPV equation (2-6.17) so that all variables are from here on non-dimensional. Dimensional variables will have subscript “*”. Consider an initial or background state of purely zonal flow Uo(y,z) with stream function (y,z): Uo(y,z) = /y (note that Vo = /x = 0) (2-7.2) Page 75 Introduce perturbation “ ”, so that the total, time-dependent stream function is: (x,y,z,t) = (y,z) + (x,y,z,t) (2-7.3) Fig.2-7.1. Schematized diagrams showing various scenarios when exchange of fluid parcels “A” and “B” along the blue path line either leads to increased potential energy so that the mass center of the system moves up (stable), is unchanged (neutral) or to decreased potential energy so that the mass center of the system moves down (maybe unstable). It is easy to see that only for the special path within the “wedge” formed by the isopycnal (red line) and the horizon (dashed line) in “(iii)” can the exchange possibly lead to instability (such that after the exchange the two parcels move away from each other). Such an instability is called baroclinic instability. The function is perturbation to the initial state ; it represents the structure of the evolving perturbation field. Substitute (2-7.3) into (2-6.17): (/t + Uo/x y/x + x/y){q + 2/y2 + y + [S-1/z]/z} = 0 Page 76 or (/t + Uo/x)q + J( ,q) + x/y =0 (2-7.4) where q(x,y,z,t) is the perturbation PV: q= 2 + [S-1( /z)]/z (2-7.5a) and o = 2/y2 + y + [S-1/z]/z (2-7.5b) is the background or initial state’s PV, and its meridional gradient is (using eqn(2-7.2)): o/y = 2Uo/y2 + [S-1Uo/z]/z (2-7.5c) Notice how the thermal-wind relation leads to the vertical shear Uo/z that then appears in this last equation for o/y. How does the structure of Uo(y,z) determine the evolution of the perturbation field ? That is, given a particular background or initial state Uo(y,z), will the perturbation “injected” on the flow grows or decays? If grow, then the initial state is unstable with respect to the perturbation . To show that Uo is stable, we must check all possible ’s. On the other hand, to show instability, we only need to find one perturbation to which the initial state Uo is unstable. Linear Stability Analysis: We assume that | | << 1 so that the J( ,q) in (2-7.4) is dropped: (/t + Uo/x)q + x/y =0 (2-7.6) The boundary conditions are that vertical velocity w* = wo* + w1* + O( 2) is zero at z=0 (surface) and z=-1 (ocean’s bottom). Here = U/(foL) is the Rossby number which is assumed to be small, and we already noted previously (see eqn.(2-6.4) and discussion) that wo* = 0, so that w* w1* The non-dimensional w1 can then be expressed in terms of using equations (2-6.11c), the hydrostatic relation (2-6.10b), and (2-6.15) which gives: b*’ = fo*/z*, hence w1* = do*/dt{fo*/z*/N2} (2-7.7a,b) and w1 = do/dt{ S-1/z} = 0 at z=0 & -1. (2-7.8) Page 77 A perturbation energy equation can be derived and it can be shown that: s (o/y)(<2>/t) dydz = 0 (2-7.9) where is the meridional displacement of fluid elements defined by: /t + Uo/x = /x (2-7.10) and <.> is a zonal-averaging. From (2-7.9), we see that if there is to be a growth in the displacement of fluid elements in time, i.e. if <2>/t > 0, then o/y must be somewhere positive and somewhere else negative in the yz-plane, or o/y can also be identically zero everywhere. In other words, o/y must vanish on a line in the yz- plane. Clearly, this condition (i.e. that o/y must vanish somewhere) is not sufficient. --------------------------- In Tulloch et al.’s paper (JPO, 2011), they use a that is a function of x,y & z: (x,y,z) = V(z)x U(z)y Substituting this into (2-6.17) then gives their eqn.2 (instead of eqn.2-7.6 above). Their second equation is just equation (2-7.8). Their equation (1) is just our equation (2-7.5c) above with 2Uo/y2 = 0 and their last term (in their eqn.1) is a consequence of applying the boundary condition (2-7.8) at the surface z=0. Page 78 Simple Ideas of Frictional Spindown and Buoyancy Shutdown Ocean is spun down by friction frictional spindown. But with stratification and bottom slope, it is possible that bottom mixing (w/stratification+slope) produces horizontal density gradients, hence thermal-wind shear in the mixed layer so that the interior flow “sees” a slippery bottom and Ekman pumping & friction are shutdown buoyancy shutdown. Rough derivations: To compute /y due to bottom slope h/y and bottom mixing (downwelling case), we consider an isopycnal =constant; then = 0 = /z h + /y y so that /y = -/z.h/y so that (g/o)/y = N .hy, where N = -g(/z)/o. Or, 2 2 with b = -g/o we have: b/y = -N2hy where N2 = b/z. By geostrophy within the mixed layer, we have fu = -p/y/o, so that fu/z = +(g/o) /y = -b/y = N2hy from previous equation. So that integrating across ML, from z=- h to z=-h+, where = ML thickness: ui – u = (-h+-z) N2hy assuming that N2 is constant. Therefore: u=ui-(-h+-z) N2hy. Note: (a) buoyancy shutdown arises only if bottom has slope h/y0; Page 79 (b) if hy is =0, then u = ui extending all the way to bottom, and the interior flow is then subject to frictional spindown (because then b = -rui operates); (c) basically, at buoyancy shutdown, the ui is canceled by thermal-wind part and b=0. References: Garrett, MacCready & Rhines, Ann. Rev. Fluid Mech, 1993, 291-323. Chapman, JPO, 2002, 336-352. Page 80 Chapter 3. Barotropic and Baroclinic Instabilities Barotropic Instability: Rayleigh-Kuo's Theorem If , then the jet is stable (For eddy, if ) Consider potential vorticity equation: (I1) Consider parallel basic flow in x direction, U=U(y) uT U u vT 0 v T Z (I2) dU Z dy (I3) After linearization and neglect quadratic and high order terms, (I1) becomes: Z U v v 0 (I4) t x y (I4)*ζ: 2 U 2 ( ) 2v (I5) t Z y Z y x Integral over a wave length, 2L=wave length Page 81 L 2 U L L L Z y )dx Z y 2 vdx (I6) 2 ( t L L The 2nd term of LHS is zero, because of the periodicity, the RHS will be: v u 1 v 2 uv v 1 v 2 uv 1 u 2 v v( ) ( u ) ( ) x y 2 x y y 2 x y 2 x (I7) L 1 2 L L Lvdx L 2 x (v u ) y (uv) dx L y (uv)dx 2 Thus, (I6) becomes: 2 L t Z y )dxdy 2L y (uv)dydx ( (I8) since uv=zero at y=±b are walls, so RHS of (I8)=0 2 ( ) t Zy dxdy 0 (I9) 2 For there to be an instability, >0, then (I9) only vanish if β+Zy changes sign somewhere in "y" t Suppose that ( Z y ) 1 -1 2 y5 4 1 y 2 (I10) -1 0 y 1 0 otherwise Page 82 Then, clearly α changes sign at some "y". Therefore, if ∂/∂t(ζ2)>0, flow is unstable, then ∫α∂/∂t(ζ2)dy=∫αdy =0, so that (I9) is satisfied. In summary, we say that if the flow is unstable, i.e. if ∂/∂t(ζ2)>0, it is necessary for α to change sign; or we can also say that "a change in sign in α" is necessary for flow to be unstable. But, "a change in sign in α" is not sufficient for the flow to be unstable. For example, the integral ∫α∂/∂t(ζ2)dy=-∫αdy =0, then, ∂/∂t(ζ2)<0, it is stable for the flow! Again, we have (I4): Z U v( ) 0 (I4) t x y u v 0u , v (I11) x y y x where ψ=perturbation stream function (note: ζ= 2ψ) so that, (I4) becomes: o 2 U v( 2 )0 (Rayleigh-Kuo Equation2) (I12) t x y o Z where (Pedlosky) U yy (I13) y y Let real{ y e ik ( x ct) } (real part of it) (I14) Page 83 be a disturbance of wave number k in x (periodic), and σ=-ikc=-ik(cr+ici)=ikcr+kci (I15) which is -i*kc(frequency), such that kcr=real frequency (I16a) kci=growth rate (I16b) Thus, flow is unstable if the imaginary part of c, i.e. ci>0 From (I14), only take the real part: u y y e ik ( xct) v x ik e ik ( xct) (I17a, b, c) 2 k 2 yy e ik ( x ct ) (I12) then becomes: ikc k 2 yy Uik k 2 yy ik y 0 (I18) oy (U c) k 2 yy y 0, or k 2 yy (U c) 0 Multiply (I18) by *, which is the conjugate of , so that * =abs( )2=real, and also 2 * ( * ) , then integrate , we have y yy y y dy (U c) Y Y Y Y y ( * 2 y ) dy k2 dy 0 2 oy 2 y (I19) Y Y The 1st term on LHS is zero, since we will assume either Page 84 (i) Y and -Y are walls so that = *=0 there, or (ii) Y and -Y are infinite and decays → 0 as y →±∞ so (I19) becomes: dy U c Y Y (U cr ici ) 2 k2 dy 0 2 oy 2 y 2 Y Y where U c (U cr ici )(U cr ici ) real , so 2 U c Y Y oy ci dy i 2 k2 (U cr ) dy 0 2 2 2 y oy (I20) Y U c 2 2 Y It means real=0, and imaginary=0. We focus on the imaginary part since we want to find out more about "ci" (which determines instability etc.), so Y y U c dy 0 2 ci 2 (I21) Y Since "c" is complex, it occurs in pair in conjugate, so that it has both positive ci (growing) and negative ci (decaying) parts. So all we need to show for instability is that ci≠0 Y y U c dy 0 2 (I21) can only be satisfied if either ci=0 or 2 (I22) Y If there is an instability, i.e. ci≠0, then Y oy U c dy 0 2 2 (I23) Y But since both U c are nonzero, oy must change sign somewhere in y: -Y≤ y≤ Y. 2 2 and This is the necessary condition that we obtained in (I9). Page 85 Y oy On the other hand, if oy does not change sign anywhere in "y", then U c 2 2 dy cannot be 0, so Y that ci must be 0, and the flow is stable. So the condition that " oy does not change sign in y" is sufficient to ensure stability. (I24) This proves the Rayleigh-Kuo Theorem. Page 86 Chapter 4 (Gill's Chapter 4; Gill’s Figures and Equations are denoted by G*): 1. What is a material element? [p.63]. [material element = infinitesimally small fluid sample that retains its identity.] 2. Why is it that if molecular diffusion is =0, the material element will consist of the same particles? What are "particles" in this case? [p.64]. [If diffusion=0, then molecules (=particles) cannot spread, so material element consists of the same particles.] 3. Why is (G4.2.1) valid (i.e. why can the LHS be broken up into the 3 terms on the RHS)? [Use d(ab)/dk = adb/dk+bda/dk).] 4. Derive (G4.2.3), (G4.2.4) & (G4.2.5). 5. What is the difference between advection and convection? [p.66] 6. Derive (G4.1.8) from (G4.3.2) and (G4.2.4) (== (G4.2.3)). 7. Derive (G4.2.5) [Use Fig.G4.2]. 8. Explain the flux vector F in (G4.3.3) - as advective + diffusive fluxes. 9. Note how (G4.3.4) [or G4.3.7] is similar to the Salinity equation used in POM. Is it identical? Why (yes or not)? 10. The entropy of a material element is fixed if the element is moved (a) without exchanging heat with its surrounding and if (b) it retains its fixed composition. The motion is then called isentropic. 11. So, the salinity or 'salt' equation is fairly straight-forward, since the Qv on p.67 = rho*s = mass of salt per unit volume; and the "mass of salt" is a 'measurable' scalar quantity - i.e. a "mass." 12. Temperature is somewhat different however. It is still a scalar of course, but it is closely tied to the agitation of molecules and therefore also to the entropy or "heat content." Temperature is derived from the "heat equation." To understand this, we need to explore a little thermodynamics. Page 87 13. The "heat" equation is obtained from the 2nd law of thermodynamics: dE = Td pdvs (G3.2.1) E = internal energy per unit mass Vs = specific volume (volume per unit mass) T and p = temperature and pressure respectively = specific entropy (entropy per unit mass) Energy increases because of (i) agitation (molecules) d, and (ii) compression dvs < 0. We assume no phase change and salinity is fixed – i.e. the fluid is of fixed composition; so = (p, T). Td change in heat content per unit mass pdvs = work done by p in compressing fluid’s volume per unit mass. 14. Eqn.(G3.2.6) is used often later – so we will derive it. From (G3.2.1): Td = dE + pdvs (N4.1) So, increase in heat content for a change in T while keeping p constant is then, from (N4.1): Cp T(/T)p = (E/T)p + p(vs/T)p (N4.2) The Cp is called the specific heat at constant pressure. Similarly, we obtain the following for increase in heat content for a change in p while keeping T constant: T(/p)T = (E/p)T + p(vs/p)T (N4.3) Now, recall that = (p, T), the Td (i.e. LHS of (G3.2.6)) is: Td = T(/T)p dT + T(/p)T dp = CpdT + T(/p)T dp (N4.4) where the definition of Cp from (N4.2) has been used for the first term. For the second term, T(/p)T, we may use (N4.3), but we can get rid of the “E” as follows. We can [T(N4.3) p(N4.2)]: Page 88 (/p)T +T( /pT)T( /Tp) = ( E/pT)( E/Tp)+p( vs/pT)p( vs/Tp)(vs/T)p so that, (/p)T = (vs/T)p (G3.2.5) which can be used in (N4.4) to give (G3.2.6): Td = CpdT T(vs/T)pdp. (G3.2.6) 15. A process is reversible if it changes very, very slowly due to infinitesimal changes in some property of the system, so that no dissipation occurs. For example, if a gas in an insulated cylinder is compressed by pushing a piston (see Fig.N4.1) and there is no friction between the piston and the cylinder wall, then if the piston is let go the compressed gas would expand back so the piston is reversed to its original position. If there is friction, then heat is permanently loss and the process is irreversible. In practice, no process is reversible. But if the applied changes are slow and/or if the system’s response is fast, then the process may be considered as being reversible. Gill assumes this [p.42, approx. line 6]. piston gas cylinder Fig.N4.1 16. Since Td = increase in heat content (of the system), then if the process is isentropic and reversible, then it is also called adiabatic. In our case, we assume “reversible processes,” so isentropic = adiabatic. Thermodynamic Definition: An adiabatic process a process in which no heat is transferred to or from the working fluid. 17a. When a parcel of fluid sinks adiabatically, say into the deep ocean, it experiences increased pressure and therefore warms up. Similarly, when air parcel rises (to mountain top) it experiences less pressure and cools down. In this case, “adiabatic” means that there is no turbulent or molecular diffusion and no radiation across the parcel’s boundary, and also no change in composition or phase. Page 89 Thus, in the case of adiabatic sinking or rising fluid parcel, because of pressure compression (expansion), DT/Dt 0; it would be nice to define a quantity such that the D(.)/Dt = 0. Also, in both cases (sinking or rising parcel), the (vertical) column of fluid (air) is apparently unstable (i.e. “warm” is under “cool”). But this is not true – i.e. the fluid is actually stable. To avoid this apparent confusion, and to make D(.)/Dt = 0, the quantity potential temperature is introduced. The potential temperature “” is the temperature of the parcel at pressure p would have if the parcel were brought adiabatically (i.e. isentropically; i.e. without exchange of heat with its surrounding) to a reference (i.e. common) pressure pr. By using , we eliminate the apparent warming or cooling due to changes in pressure. Since now the various temperatures at different depths have been brought to a common pressure, we can compare them without worrying that the ’s represent false temperatures. The variation of with height can therefore be used to judge if the vertical column is statically stable or not. The situation is a bit similar to trying to judge if person (A) on Jupiter is heavier than person (B) on Pluto (Fig.N4.2). It would appear that (A) is heavier than (B). But to more fairly judge them, we can move them (without loosing even a single hair; constant entropy) all to a reference place, say the earth, then re-weigh them. The weight measured on earth may then be called “potential weight.” (A) I’m 2 tons (B) I’m 2kg Jupiter Pluto Fig.N4.2 17b. An example (for air) and some exercises will make the concept of potential temperature clearer. It is convenient to do this with air than with water because air behaves like an “ideal gas” and explicit formulae may then be obtained. Suppose we follow an air parcel from some height with (p,T) which sinks adiabatically (i.e. isentropically) from a high mountain to the ground. For isentropic process, we have from (G3.2.6) with d = 0: Page 90 dT = T(vs/T)p,s dp/Cp = T(/T)p,s dp/(2Cp). (N4.5) Assume that air obeys the ideal-gas law: = p/(RT) (N4.6) [note that high p gives rise to high T and also smaller volume (per unit mass) vs = -1], we have then, (/T)p,s = p/(RT2) (N4.7) and therefore (using (N4.6) again): (/T)p,s /(2Cp) = +(R/Cp)/p = (2/7)/p = /p, (N4.8) where = (R/Cp) = 2/7 for dry air from Gill’s (G3.3.2). Equation (N4.6) then becomes: dT/T = dp/p (N4.9) Integrating (N4.9) from a height where pressure is p and temperature is T to the reference level where the pressure is pr and temperature is, by definition, the potential temperature , we have ln()ln(T) = [ln(pr)ln(p)], or ln(/T) = ln(pr/p)2/7, or /T = (pr/p)2/7. (N4.10) It is clear that potential temperature is only meaningful when we have prescribed a reference level; the most common one is the surface, where the pressure is 1013 mb. It is convenient to simply let pr = 1000 mb. Exercise 17b1: Air parcel “A” at a height where p = pA = 100 mb has an in situ temperature TA = 175 K. Another parcel “B” closer to the ground where p = pB = 400 mb has a TB = 220 K. Is the air column stable? [Hint: calculate the corresponding .] Please do not read beyond this line before you attempt exercise 17b1. Page 91 More Thoughts on Exercise 17B1: In the above exercise, TA < TB and according to (N4.6), the density A seems to be greater than B, i.e. A > B, i.e. air column appears to be unstable. But (N4.6) cannot be used this way to determine the densities because the pressures pA and pB are different. In fact, since pA < pB, it may be possible for A to be less than B. This is why we introduce the potential temperature by “moving” both “A” and “B” adiabatically to the common level where p = pr. By doing this “moving” adiabatically, we ensure that both parcels’ heat-gain (and rise in temperature) as they are moved closer to the ground is “kept inside” the parcels. Therefore, from (N4.10), A 338 K > B 286 K, so that the air column is in fact stable. Exercise 17b2: If we increase TB, the air column would then tend to the state of “less stability.” (we will learn the precise definition later). Calculate the value of TB such that the air column is neutrally stable. Exercise 17b3: In Exercise 17b1, is the air column stable if TB = 300 K? Exercise 17b4: Figure N4.3 plots T as a function of p for three values of ’s. Think about this in view of the exercises you have just done. By choosing a few numbers, convince yourself that you have understood the concept of . Fig.N4.3. The drop in local temperature (T) with height (or pressure) for various potential temperatures = 240, 280 and 320 K, computed from (N4.10), rewritten as T = (p/pr)2/7, pr = 1000 mb. Page 92 Answers to Exercises 17B2&3: TB_neutral = 260 K; unstable. 18. Static Stability (section G3.6). Suppose the fluid is at rest; to judge a parcel’s static stability, we move the parcel up or down isentropically, calculate its density change, d, then compare this change with the change in the surrounding (i.e. ambient) density, da. If the parcel is moved up, then for stability, d must > da so that the parcel can sink back to its original depth when let go. The general formula for the change in density d(p,T,s) of the fluid (ambient or parcel fluid) as a result of a change in height “dz” is, by chain-rule: d(p,T,s) (d/dz).dz = [ (/p)T,s (dp/dz) + (/T)p,s (dT/dz) + (/s)p,T (ds/dz) ].dz (G3.6.8) where s = salinity. For a parcel of fluid that is moved isentropically, with = (p, T), and fixed s; see “13”, (G3.6.8) becomes: d = (/p)T,s dp + (/T)p,s dT (N4.11) We can replace “dp” and “dT” with “dz” as follows. From (G3.2.6), for isentropic process (d = 0): dT = T(vs/T)p,s dp/Cp = T(/T)p,s dp/(2Cp) = Tdp/(Cp) = dz (G3.6.1) where = thermal expansion coefficient and = adiabatic lapse rate: = -1 (/T)p,s > 0 (G3.6.2) = gT/Cp, (G3.6.5) and the last equality in (G3.6.1) (involving ) is obtained by using the hydrostatic equation: dp/dz = g (N4.12) Substituting “dT” from (G3.6.1) and “dp” from (N4.12) into (N4.11), which becomes: d = [g(/p)T,s + ] dz (G3.6.7) Page 93 Reminder: equation (G3.6.7) is the change in density, dpar, of the fluid parcel as it is moved. On the other hand, equation (G3.6.8) is general, and gives (for example) the change, damb, of the ambient (i.e. surrounding) fluid density with z. Suppose that the parcel is moved upward (i.e. dz > 0), we must have, for stability, that dpar > damb, so that the parcel can move back (i.e. downward) to its original position when let go: [ + dT/dz] s ds/dz > 0 (G3.6.9) or Cp-1 2gT + dT/dz s ds/dz > 0 (G3.6.10) where s = -1 (/s)p,T > 0 (G3.6.11) is the expansion coefficient for salinity (I use s to avoid possible confusion with the planetary beta ). Exercise 18.1: Show that the same stability condition is obtained if similar arguments as above are applied but the parcel is moved downward instead of upward. Answers: If parcel is moved down, then for stability we must have dpar < damb. Then from (G3.6.7) and (G3.6.8): [g(/p)T,s + g(/p)T,s + + dT/dz sds/dz] dz < 0. But since dz < 0, the […] must be > 0, and (G3.6.9) is obtained. 19. The buoyancy frequency (or BruntVaisala frequency; section G3.7) squared is: N2 = g [ ( + dT/dz) s ds/dz] = g[Cp-1 2gT + dT/dz s ds/dz]. (G3.7.1) When N2 > 0 (< 0), then the medium is stable (unstable). In a stable medium, N is real and represents the frequency of (vertical) oscillation of a fluid parcel – closely linked to internal waves. In the ocean N 10-4 to about 10-2 s-1, or periods of about 10 minutes to a few hours. 20. Since is the temperature that a parcel has if it is moved isentropically to a reference pressure (see “17a”), then a particular has a particular entropy . Therefore = constant when = constant. So if s = fixed, = (), i.e. is a function of only. The function can be obtained by fixing p = pr (reference pressure) in (G3.2.6) which then becomes, since T is then = : Page 94 d/d = Cp(pr, )/ (G3.7.6) 21. Assume that s = constant (i.e. fixed composition). Instead of p and T, we can use p and , so that = (p, ). Then for a parcel that is moved isentropically (i.e. d = 0), the change in its density is (c.f. N4.11) dpar = (/p),s dp + (/)p,s d (N4.13) The change of the surrounding density is from (G3.6.8) with T replaced by : damb(p, ,s) (d/dz).dz = [ (/p),s (dp/dz) + (/)p,s (d/dz) + (/s)p,T (ds/dz) ].dz (N4.14) For stability, dpar damb > 0 (if dz > 0) as in “18”, then: g-1N2 = ’(d/dz) ’s (ds/dz) > 0 (G3.7.9) where ’ = -1(/)p,s and ’s = -1(/s)p, . In the form written as (G3.7.9), if s = fixed, then the stability condition becomes very simple, i.e. the fluid medium is stable if the potential temperature increases upwards. This fact was used in Exercise 17b.1-4. 22. The Heat Equation (isentropic motion): If fluid elements do not exchange heat with their surroundings, and further more if they retain a fixed composition, then their motions are adiabatic or isentropic, and we have: D/Dt = 0 = Cp(pr, )-1 D/Dt (G4.4.1) the last equality is from (G3.7.6). Thus for adiabatic motion, the heat equation takes on this simple form: D/Dt = 0, (N4.15) Page 95 i.e. the potential temperature is constant following the parcel’s movement. The heat equation looks a little more complicated if written in terms of the in situ temperature. From (G3.2.6), we have TD/Dt CpDT/Dt (T/)Dp/Dt = 0. (N4.16) The heat equation can also be written in terms of the internal energy E (per unit mass; why? See the pressure work term in eqn.G3.2.1). From (G3.2.1): TD/Dt DE/Dt + pDvs/Dt = 0. (N4.17) This last equation can be re-written by making use of the continuity equation (G4.2.4) which is: /t + .(u) = 0, (G4.2.4) so that DE/Dt + E*LHSof(G4.2.4) = t + u. E + E*LHSof(G4.2.4) = (E)/t + .(uE); or: DE/Dt = (E)/t + .(uE). (G4.3.6) Then, times the second equation of (N4.17) gives: DE/Dt (E)/t + .(uE) = pvs-1Dvs/Dt = p .u, (G4.4.3) in which (G4.3.6) and (G4.2.2) have been used. Equation (G4.4.3), (E)/t + .(uE) = p .u, (N4.18) says that the rate of increase of energy per unit volume following a fluid parcel, D(E)/Dt > 0, is equal to the rate at which work is being done due to the parcel’s compression (i.e. convergence) .u < 0. Similarly, there is a decrease in E if the parcel expands (i.e. divergence) .u > 0. By the way, one can use Figure G4.2 and its derivation of the time rate of change, /t, of mass per unit volume (= ), to apply to energy per unit volume = E. 23. The Heat Equation (non-isentropic motion): Page 96 Then additional terms are added to account for various heat transfers and sources: Frad, k T and QH. Similar to the “advective heat flux density” Eu, we can identify Frad and k T as the “radiative and conductive heat flux densities” respectively. In other words, they are “fluxes” that ‘flow’ in and out of the side faces of the rectangular volume element of Figure G4.2 (see the comment of the last sentence of “22” above). By contrast, the “QH” (= rate of heating per unit volume) is a heat source (or sink) inside the volume element and is therefore grouped together with the pressure work p .u. The result is therefore equation (G4.4.4): (E)/t + .(Eu + Frad k T) = QH p .u. (G4.4.4) An alternate form using the in situ temperature T is obtained by equating times (N4.16) to times (N4.17), and then using (G4.3.6) and (G4.4.4), to obtain: CpDT/Dt T Dp/Dt = .( k T Frad) + QH. (G4.4.6) Yet another alternative form using the potential temperature is simpler, using (G4.4.1), (N4.16) and (G4.4.6): TCp(pr, )-1 D/Dt = .(k T Frad) + QH. (G4.4.7) The Princeton Ocean Model (POM) uses (so do most other ocean models), hence it uses equation (G4.4.7). In fact an approximate form is used, in which QH = 0, the k is replaced by the Smagorinsky’s horizontal diffusivity coefficient Asmag (why??), and Frad represents the radiative penetration of light from the ocean’s surface and is significant only near the surface z > 30 m (where T): D/Dt = .(Asmag ) (Frad)/z.(Cp)-1, (N4.19) where .(Asmag ) is supposed to represent mixing effects by small-scale eddies; we replace T by , and also lump the (Cp)-1 into Asmag (so that it has a unit of m2s-1) thus acknowledging that we don’t know much about this small-scale eddy mixing process, and therefore the errors due to these “replacement” and “lumping” of terms are inconsequential. In actual numerical simulation, we set Asmag to be very small, in fact often zero. 24. Momentum Balance (or the Equation of Motion): This is (G4.5.20) [or pom_v13.pdf (or *.doc) equations (1-24a,b,c)]: Page 97 (Du/Dt + 2 × u) = p’ ’g + .( u), (G4.5.20) where u = (u, v, w); = (1, 2, 3); g = (0, 0, g). (N4.20) 25. Mechanical Energy (Kinetic Energy) Equation: This is (G4.6.1) [with (G4.6.2)]; it is derived by taking the scalar product of “u” with eqn.(G4.5.20) i.e. u.(G4.5.20). Term by term, we have: u.Du/Dt = u.(u/t + uu/x + vu/y + wu/z) = D(u.u/2)/Dt [decompose into components] u.(2 × u) = 0 [the 2 vectors are perpendicular, so their dot product = 0] u. p’ = .(u p’) + p’ .u [by chain rule] u. (’g) = wg’ u.[ .( u)] u.[ .( u)] = .(u. u) ( u . u) = .( (u.u/2)) [(u/x)2+(u/y)2+(u/z)2] Putting all these together, we have then for u.(G4.5.20): D(u2/2)/Dt = wg’ .(u p’) + .( (u2/2)) [(u/x)2+(u/y)2+(u/z)2] + p’ .u (N4.21) or, D(u2/2)/Dt = wg’ + .(u p’ + (u2/2)) + p’ .u (G4.6.1) where = [(u/x)2+(u/y)2+(u/z)2] = dissipation rate (G4.6.2) 26. It is also useful to think of the rate of change of K.E. for a fixed volume – Gill’s equations (G4.6.3- G4.6.4). Also understand (G4.6.5) and interpretations based on Fig.G4.6. Page 98 27. For a large volume, we have eqn.(G4.6.7): dK/dt + Fn’dS = dxdydz, (G4.6.7) where Fn’ = n.F’ for F’ = energy flux density given by (G4.6.4), and n is the unit vector normal to the bounding surface dS. 28. Importance of Viscosity! We normally say (or hear the statement) that viscosity is not important on the large scales (say a few kilometers or larger) typical of the energy input (say by wind). But let’s think about this more carefully. Let’s take the simple case of a homogeneous ocean. Input of energy is due to the pressure and viscous stress (see equation G4.6.4) at the free-surface only, since no energy can pass through the solid bottom and side walls. Yet despite this constant energy input, we know that the ocean’s energy is not continually increasing, i.e. dK/dt 0 in (G4.6.7). Therefore, in that equation, the large-scale energy flux input Fn’dS can only be balanced by the small-scale energy dissipation dxdydz by molecular viscosity! How small is small? Let’s call this small scale = l. The l can be assumed to depend only on and : l = ( 3/ )1/4 O(10-3 m) (N4.22) Therefore, because of the very disparate scales, the only way that the large-scale energy can be dissipated by the small-scale (l ) viscous dissipation is that there will be a transfer of energy from large to small. This transfer is accomplished by the nonlinear terms (e.g. uu/x etc; if u ~ eikx, then uu/x ~ ei2kx etc) – “breaking” down large eddies into smaller ones, which are in turn broken down again into even smaller ones etc until the smallest eddies have scales of O(l ) when molecular viscosity can act to dissipate the (kinetic) energy into heat. The fact that we have disparate scales creates big problem for numerical models which typically must use some kind of eddy viscosity >> molecular viscosity in order to dissipate input of large-scale (say wind) energy. However, this “eddy viscosity” parameterization does not guarantee that we are removing the energy in a realistic way – thus we have (for example) Mellor-Yamada scheme or the Smagorinsky viscosity to try to model this “energy-removing” process. If we have a super-super computer that allow a grid resolution of 10-3 m, then we will not need to use Mellor-Yamada or Smagorinsky! Page 99 Chapter 5. Centrifugal and Slantwise Instabilities: Applications to Mixed Layer Recall that when warm fluid “A” is below cold fluid “B,” the system is unstable if the potential temperature of “A” is higher than that of “B.” In this case, the system is unstable because of gravity – i.e. lighter fluid tends to rise and heavier fluid tends to sink. This is called convective instability. This topic was discussed when we studied Gill’s Chapter 4. An interesting analogous instability occurs in rotating flows. Recall that on a rotating earth, centrifugal acceleration may be “lumped” together with gravity such that a gravitational potential may be defined (see Fig.1-7). This suggests that in rotating flow, centrifugal force acts like gravity and there may exist a similar kind of instability. In section 5-1 we study the nature of this instability for a homogeneous fluid (constant density) – centrifugal instability. We then extend the ideas to the case with density stratification in section 5-2 slantwise or symmetric instability (Emanuel, 1994). In section 5-3, we then utilize the same parcel theory used in section 5-1 to re-derive the growth rates of both instabilities, and also to outline same for baroclinic instability. Finally, in section 5-4, we discuss the applications of the theories to instabilities that can happen in a mixed layer where horizontal stratification can exist and the vertical stratification is weak. 5-1 Centrifugal Instability (CI) We examine axis-symmetric flow – a circular eddy for example – in which there exists an azimuthal velocity v(r,t), positive cyclonic and a radial velocity u(r,t), positive outward from the center. Here r is the radius from the axis of rotation and t is the time. All variables are independent of the azimuthal angle (Fig.5.1). Page 100 Figure 5-1. Axis-symmetrically rotating homogeneous fluid with angular velocity = v/r. The radial displacement of a circular ring is also sketched. The motion is characterized by and M: = angular velocity = v/r (5-1.1) M = angular momentum = rv = r2 (5-1.2) The here should not be confused with the same symbol used for earth’s rotation rate in the previous chapters. Assuming inviscid fluid, the azimuthal momentum equation is: Dv/Dt + uv/r = o (p/)/r (5-1.3) and the radial momentum equation is: Du/Dt v2/r = o (p/r) (5-1.4) where o = 1/o (not to be confused with the thermal expansion coefficient used in Chapter 4, = -1(/T)p,s). [If you are familiar with the y and x-momentum equations involving the Coriolis parameter “f,” then it is easy to remember (5-1.3,4), just replace the “f” by “” and use equation (5-1.1); also replace “y” by “r” and “x” by “r.”] We now show that, following a fluid parcel, the angular momentum “M” is conserved, i.e. DM/Dt = 0, for axis-symmetric flow. To show this, we multiply (5-1.3) by r and expand the D/Dt written in (r,) coordinate [i.e. D/Dt /t + u/r + (v/r)/]: rv/t + ruv/r + vv/ + uv = o (p/) (5-1.5) We then note that DM/Dt = D(vr)/Dt is: DM/Dt = M/t + uM/r + (v/r)M/ = rv/t + uv + ruv/r + vv/ (5-1.6) which is the same as the LHS of (5-1.5). Therefore, (5-1.5) gives : DM/Dt = o (p/) = 0, (5-1.7) for axis-symmetric flow, since then “p” (also all other variables) is independent of . Page 101 When we analyzed convective instability (Chapter 4), we examined a parcel’s movement in the vertical (z) within a background of (ambient) density stratification. We then moved the parcel vertically and adiabatically to determine if at the new position the parcel is lighter or heavier than the ambient fluid. A similar analysis is now carried out for a parcel undergoing axially symmetric displacements. Instead of “z” we now have “r,” and instead of potential temperature, we now have that the angular momentum (M) is conserved. The parcel in the present axially symmetric situation is a “ring.” So, what is the ambient, equivalent “hydrostatic” state? This is obtained from (5-1.4) by setting Du/Dt = 0 (similar to setting “Dw/Dt = 0” in hydrostatic case): Mo2/r3 = o (dpo/dr) (5-1.8) where Mo = rvo, and po is the pressure distribution that maintains the balanced azimuthal flow having an angular momentum = Mo. Suppose now that the flow is initially at this balanced state, and consider the ring-parcel radially displaced by r from r = ro. For small displacement, we assume also that “p” remains at its initial pressure distribution po. Equations (5-1.4) and (5-1.8) then give, for the ring-parcel: Du/Dt = [M2 Mo2(ro+r)]/r3 (5-1.9) Here, u and M represent the fluid state after the displacement at r = ro+r. Since M is conserved (i.e. equation 5-1.7), we have: M2 = Mo2(ro) (5-1.10) since the ring-parcel initially (i.e. before being displaced) has M = Mo. By Taylor’s expansion about ro assuming that r is small, we also have, Mo2(ro) Mo2(ro+r) [(dMo2/dr)|ro].r (5-1.11) where the vertical-bar symbol “|ro” means “at ro.” Equations (5-1.10) and (5-1.11) then give M2Mo2(ro+r) [(dMo2/dr)|ro], and equation (5-1.9) then becomes: [Du/Dt]|(ro+r) [(dMo2/dr)|ro].r/ro3 (5-1.12a) or simply (with the “at” symbol “|” being understood): Page 102 Du/Dt (dMo2/dr).r/ro3 (5-1.12b) Equation (5-1.12) is the result we want. It shows that as the ring-parcel is displaced outward, r > 0, it is further accelerated outward (and therefore moves further away from its initial position ro), Du/Dt > 0, if the squared ambient angular momentum Mo2 decreases with r, i.e. if (dMo2/dr) < 0. Thus, if (dMo2/dr) < 0, then the flow is unstable (5-1.13a) This kind of instability is called centrifugal instability (CI). On the other hand, if (dMo2/dr) > 0, then the flow is stable (5-1.13b) Suppose the ambient aximuthal velocity varies with “r” as: vo = rn (5-1.14) then Mo2 = r2(n+1), and dMo2/dr = 2(n+1).r2n+1 < 0 if n < 1, (for vo = rn) (5-1.15) Therefore, an (axially symmetric) eddy which “does not spin fast enough,” i.e. whose velocity distribution decays faster than 1/r, is unstable to axially symmetric perturbations. Exercise: Go to a playground with your friend(s) or someone you may trust. Stand on a “merry- go-round (MGR; see fig.5-2).” Ask the trusted friend to spin the MGR, then jump off while the MRG has spun to an almost steady rotational speed. Describe and explain what you experience the moment you touch the firm (almost stationary) ground. By the way, why is the ground “almost stationary”? (Warning: this exercise may not be suitable for teachers at or over 55). Page 103 Figure 5-2. A merry-go-round. Localized Centrifugal Instability (LCI): The above theory requires axially symmetric perturbations. But as the Merry-Go- Round exercise suggests, localized perturbation (the exerciser jumping off) may also induce localized centrifugal instability (LCI). This is in fact true for ocean (and atmosphere). To develop the relevant formulae, we now define localized conditions for CI based on the inviscid equations of motion on an f-plane. We will develop the theory carefully, step by step. In equation (1-24), we (a) ignore the horizontal components of the earth’s rotation: 1 and 2; (b) align the z-axis so that it is perpendicular to the plane tangent to the local earth’s surface where the latitude is o, i.e. so that z is locally “upward” and 23 = 2sin(o) = f, where = magnitude of the earth’s angular velocity; (c) assume motions of scales small enough that “f” does not vary significantly, i.e. the motions are confined to the neighborhood of latitude “o”; and, (d) ignore the viscosity . Page 104 We then divide (1-24) through by the density , and set = 1/; so that equation (1-24) becomes: (5-1.16a) (5-1.16b) (5-1.16c) Let us assume that in our flow field, there exists a direction yg such that p/yg = fug = constant (could be zero!), we can then align the “y” in (5-1.16) in this direction, i.e.: p/y = fug = constant (5-1.17) where ug is therefore the (constant) geostrophic flow in the (new) x-direction. Equation (5-1.16) then becomes: (5-1.18a) ) (5-1.18b) (5-1.18c) We can now translate the coordinate system in the x-direction at the constant speed ug, so that with respect to this translating coordinate system the new “u” and its material derivative are given by: unew = u ug, and Dunew/Dt = Du/Dt (5-1.19) The other velocity components, v and w, in the directions perpendicular to the x-axis, are clearly not affected by the uniform x-translation (assuming negligible Einstein’s relativistic mechanics!). Also xnew = x ug.t, so that /xnew = /x. Equation (5-1.18) is then (dropping the “new” etc): (5-1.20a) Page 105 (5-1.20b) (5-1.20c) In this new aligned and translating coordinate system, we note that the y-momentum equation (5-1.20b) has no pressure forces, so it becomes: D(v + f x)/Dt = DM/Dt = 0 (5-1.21) where M=v+fx (5-1.22) is the absolute momentum. Therefore the absolute momentum M is conserved for inviscid motions on an f-plane on which one component of the geostrophic velocity (ug in our case) is constant. While ug is constant, the geostrophic velocity in the y-direciton, i.e. vg is quite free, and in general is a function of x and z, and is defined as (see fig.5-3): vg = (p/x)/f (5-1.23) Page 106 Figure 5-3. A sketch showing ug and vg, and also the tube-parcel elongated in the y- direction and given a perturbed displacement x in the x-direction. We consider then the flow situation in Fig.5-3. Suppose a perturbation is introduced. In order that M (= v + f x) remains conserved, the perturbation must be such that the ug remains constant. From (5-1.17), it follows then that p’/y = 0 (5-1.24a) where p’ is the perturbation pressure. It then follows that (perturbation)/y = 0 (5-1.24b) and the perturbation must be highly (in fact infinitely) elongated in the y-direction as the tube-parcel of fluid sketched in Fig.5-3 shows. When the tube-parcel is displaced in the x-direction, the governing acceleration is obtained from (5-1.20a) which upon using (5- 1.23) becomes: (5-1.25) where the approximation is from neglecting the perturbation pressure from p/x (i.e. the pressure remains at its initial pressure distribution that produces vg) and Mg = vg + f x (5-1.26) is called the geostrophic absolute momentum (in analogous to the absolute momentum in equation 5-1.22). In (5-1.25), the u, M and Mg are the values at the tube-parcel’s displaced position x + x. Then comparing (5-1.25) with (5-1.9) and following the analysis that leads to (5-1.12b), we now have: Du/Dt (Mg/x).x.f (5-1.27) which is the result we want. It shows that as the tube-parcel is displaced to the right, x > 0, it is further accelerated to the right, i.e. Du/Dt > 0, if (Mg/x) < 0. Thus, if (Mg/x) < 0, i.e. if vg/x + f < 0, then the flow is unstable (5-1.28) Page 107 The growth rate may be estimated by noting that u x/t, so that (5-1.27) then gives: cf2 = (t)-2 f.( vg/x + f) (5-1.29) Equation (5-1.28) shows that only anticyclones can have centrifugal instability, especially at lower latitudes. Page 108 Gill's Chapter 5: Adjutment under Gravity in a Nonrotating System Section 5.1: Many fluid mechanical phenomena (including GFD) involve studying how the fluid system (e.g. the atmosphere-ocean system) tends to adjust to equilibrium; i.e. how it approaches a steady state once it is initially perturbed. What are the characteristics of the equilibrium (e.g. a bucket of water or a small lake at rest; the Kuroshio or the atmospheric jet-stream in perfect geostrophic balance; etc)? How long does the adjustment take? How may the adjustment process be understood? In this chapter, we study the adjustment without the complications due to the earth's rotation and curvature, and furthermore we only consider small departures from the hydrostatic state - i.e. we will linearize the equations. Fig.5.1 -- left (side X) has density 2, right (side Z) has density 1 < 2. The response is similar to the familiar "estuarine circulation;" Fig.5.2a: pA pB = c gh; (N5.1) Fig.5.2b: Let "F" be the z-position of the interface that is on the left side of the partition. Then: pA = pF + 2 gh, and pB = pF + 1 gh; so that: pA pB = (2 1) gh = 2 g’h (N5.2) where g' = (2 1) g/2 is the reduced gravity. Exercise N5.1. Why is the adjustment slower in the case of Fig.5.2b when compared with that of Fig.5.2a? Section 5.2: Consider an ocean of constant depth H and constant density c with an atmosphere of zero density above, so that the in-situ density of the system is: = c for H < z 0 Page 109 = 0 for z > 0. A Cartesian coordinate (x, y, z) is chosen such that z = 0 is at the free surface at rest, x points eastward and y northward. The z-momentum equation applied to this equilibrium (i.e. at rest) fluid is then: dpo/dz = g, or po(z) = gz. (5.2.1) Section 5.3 -- Surface Waves: The governing equation is just the Laplace Equation (5.2.7). Note that this is NOT a wave equation – it is not even time-dependent! It is in fact just a statement of incompressibility condition (5.2.4). The “wave” part comes from the surface boundary condition equation (5.2.10) which gives the time-dependent motion to the water beneath to generate waves. The bottom boundary condition (5.2.8) restricts wave motion at the ocean’s floor. But actually, this condition is less important physically especially when the ocean is very deep – why? Also, why are there no boundary conditions on the sides? What implicit assumption is Gill making by ignoring side boundaries? To see this, derive the dispersion relation (5.3.8) in section 5.3. 2 = g.tanh(H). (5.3.8) Note that the phase speed is c = / = [(g/)(tanhH)]1/2, so that long waves with small travels faster than short waves with large . For H << 1, we have tanh(H) ~ H, so that: c = / (gH)1/2 (N5.3) For H >> 1, we have tanh(H) ~ 1, so that: c = / (g/)1/2. (N5.4) Page 110 What is the meaning of H << 1 and H >> 1? Derive the expressions for u and v: Equation (5.2.5) and (5.3.6) give: u/t = gkosin(kx+lyt)cosh[(z+H)]/cosh(H), and similarly for v; so that: (u,v) = (k, l).(go/)cos(kx+lyt)cosh[(z+H)]/cosh(H) (N5.5) Let’s consider only waves propagating in the xz-plane (i.e. the y-wavenumber l = 0). Also consider infinitely deep ocean, i.e. H ~ . Then (N5.5) becomes: u (kgo/).cos(kxt).exp(kz) (N5.6a) and (5.3.7) gives: w (kgo/).sin(kxt).exp(kz) (N5.6b) Exercise: Plot on the xz-plane paths experienced by fluid parcels due to the wave motion. Hint: x ~ u.dt and z ~ w.dt. Some Useful Things to Remember for Deep-Water Waves: Condition for “deep water” is: H >> 1, or H >> /(2) (N5.7) i.e. water depth deeper than about one wavelength . For example, at z = /2, equation (N5.6) gives the amplitudes for u and w to be approximately exp(2) 0.04, or 4% of the amplitude at the surface, and at z = , exp(2) 0.002, which is only 0.2% of the amplitude at the surface. From (N5.4), we have for deep-water wave phase speed in m/s: Page 111 c = / (g/)1/2 = [g/(2)]1/2 1/2 1.25 1/2 (N5.8a) and also for period tp in seconds: tp = /c 0.80 1/2 (N5.8b) The Table below summarizes the above results. Table 5.1 Typical values of deep-water waves (m) c (m/s) tp (s) Water deeper than H (m) ~ /2 1 1.25 0.8 0.5 10 4 2.5 5 100 12.5 8 50 1000 40 25 500 Page 112 Gill's Chapter 6: It is perhaps best to approach the concept of normal modes in a continuous system, then we will go back to the two layer model. I will rely on your knowledge about the Fourier theory. Mr. Fourier basically asserted that a continuous (maybe piecewise) function can be represented by a simple sum of infinite number of terms comprising of the sines and cosines: (6.1) where the an and bn are the Fourier coefficients. The sines and cosines are called the basis functions; together they form a complete basis – which loosely means that together they are sufficient to represent F(x). The sines and cosines are solutions of many boundary/initial-value differential equations – simple harmonic oscillator, seiches in channel etc. We will find that normal modes for motions in a stratified sea of constant depth play the same role as the sines and cosines in these other problems. Consider the following linearized equations: u/t fv = P/x + x/z (6.2a) v/t + fu = P/y + y/z (6.2b) (p/z) /o = g/o = +b (6.2c) /t + wdo/z = 0, or b/t + N2w = 0 (6.2d) u/x + v/y + w/z = 0 (6.2e) where P = p/o is the “kinematic” pressure, x and y are also kinematic stresses, i.e. stresses divided by o, N2 = (g/o)(do/dz) is the buoyancy frequency, o(z) = background density, and (u, v, w, p, , b) are all small perturbations representing departures from a static fluid at rest. We consider that the motion is initiated from rest, so that (u, v, w, p, , b) = 0, at t=0 (6.3) Equation (6.2d) then suggests that we define (6.4) Page 113 so that b = N2 , and (6.2c) becomes (p/z) /o = N2 (6.5) Now, the total pressure PT is the sum of a hydrostatic part and the perturbation: PT = gz + P (6.6) At the free surface z = (x,y,0,t) = (x,y,t), PT is the atmospherics pressure, taken to be zero. So that (6.6) gives, P = g, and P/t = g/t, at z 0 (6.7) Integrating (6.5) w.r.t. z from z to 0 then gives: (6.8) In addition to the boundary condition at z=0, we also have at the bottom: w=0 at z = H (H=constant) (6.9) Normal Modes: We seek separable solutions of the homogeneous part of (6.2) (i.e. without forcing ): (u, v) = (Uk(x,y,t), Vk(x,y,t)).dQk/dz, k=1,2,3,… (6.10a,b) where “Q(z)” is a function of z only. From (6.2e), we have: (Uk/x + Vk/y).dQk/dz = w/z (6.11) which suggests that the “z” part of w be “Q”, or since w is related to through (6.4), that = Sk(x,y,t).Qk(z) (6.10c) Page 114 Substitute (6.10) into /z of (6.2a) (for example; with = 0) gives, using also P/z from (6.8): LHS of (6.2a) [Uk/t fVk].d2Qk/dz2 (6.12a) RHS of (6.2a) [Sk/x].N2Qk(z) (6.12b) For both sides to be consistent, the “z” part must match to an arbitrary constant factor, say “k”, so that: d2Qk/dz2 + N2k.Qk(z) = 0 (6.13) This is a simple ODE, and is supplied with 2 BC’s one at z=0 (i.e. 6.7) and the other one at z=H (i.e. 6.9). A particularly simple case is, instead of (6.7), we also take w = 0 at z = 0: Qk = 0 at z = 0 and at z = H (6.14) The solution of (6.13) consists of sines and cosines, which are the normal modes or eigenvectors of the linearized, shallow-water equations (6.2) in an ocean of constant depth. Substituting (6.10) into (6.2) [/z the momentum eqns. first], and making use of (6.13): Uk/t fVk = ck2 Sk/x (6.15a) Vk/t + fUk = ck2 Sk/y (6.15b) Uk/x + Vk/x = Sk/t (6.15c) where we have written ck2 = 1/k. These equations are exactly the same as the equations the depth-integrated equations of a homogeneous fluid in an ocean basin of constant depth. In fact, we can identify Sk to be the “” and ck2 = gHk (6.16) which can be named the square of the equivalent shallow-water wave phase speed for the mode “k” in an ocean of equivalent depth Hk. Thus the solution to this stratified fluid problem consists of just solving a set of shallow-water equations, one for each “k”, then summing them all up to get the complete solution. Page 115 In Fourier series, we calculate the coefficients an and bn by using the orthogonality property of of sine and cosine: ; (6.17a) ; (6.17b) . (6.17c) The first equality on the LHS of (6.17a) is derived by first noting that , then using ei(m+n)x = [cos(mx).cos(nx)sin(mx).sin(nx)] + i[sin(mx).cos(nx)+ cos(mx).sin(nx)] to equate the real part. The second equality, i.e. that the integral is = 0, is shown by integration by parts, say of , which gives , so that in order to be consistent with the first equality we must have, if m n, each of the two integrals = 0. Equating the imaginary part then gives (6.17b). Equation (6.17c) is easily derived. We first have by integration by parts. Also, since cos(m+m)x = cos2mx sin2mx = 1 2sin2mx, we have ; and also . Using equations (6.17), the Fourier coefficients an and bn can be determined for the function F(x) in (6.1). The requirement is that the basis functions (sines and cosines) be complete and orthogonal in the interval [-, ]. An exactly parallel approach can be used for the normal modes. We now show that they are orthogonal. Multiply (6.13) by Qm, take the integral from H to 0, integrate by parts: Qmd(dQk/dz) + N2kQmQkdz = 0, or, (dQm/dz).(dQk/dz)dz + N2kQmQkdz = 0. Page 116 Interchange “m” and “k” in the above equation and subtracting the two resulting equations: (km). N2 QmQkdz = 0, or N2 QmQkdz = 0, if k m. (6.18a) so that from the line before “Interchange “m” and “k”..”, we also have (dQm/dz).(dQk/dz)dz = 0, if k m. (6.18b) Comparing this with (6.17), we see that the normal-mode eigenfunctions Qm (or dQm/dz) are orthogonal in the interval z = [H, 0]. Thus any (piecewise continuous) function in “z” may be expressed as a sum of Qm’s with coefficients determined similar to the way the Fourier coefficients are determined. In particular, the frictional force, x/z and y/z may be expressed in terms of the dQk/dz’s: x/z = g k k x . dQk/dz, y/z = g k k y . dQk/dz (6.19) where the kx,y can be functions of (x,y,t). Given the functions x/z and y/z, the x,y k are obtained just like we obtained the Fourier coefficients an and bn: x,y k = x,y/z. dQk/dz dz/(g( dQk/dz )2dz). (6.20) Thus instead of (6.15), we now have: Uk/t fVk = ck2 Sk/x + g k x (6.21a) Vk/t + fUk = ck2 Sk/y + g k y (6.21b) Uk/x + Vk/x = Sk/t (6.21c) This completes the formulation of the linearized shallow-water equations (6.2) in water of constant depth H in terms the normal modes. Instead of the velocity (u,v) and density, we have transformed the problem into one that requires the solution of the “transports” (Uk, Vk) and “elevation” Sk for each of the mode “k.” We will now show that in a two- layer system that we studied in Gill’s chapter 6, we use similar “transport” and “elevation” idea. Page 117 Chapter 7: Notes on Island Rule Island rule (Godfrey, 1989) is derived by an application of the Kelvin’s circulation theorem. We will first learn this theorem, following Pedlosky’s book (1979), then we will learn island rule by choosing a very special circuit. Kelvin’s Theorem Take a curve C enclosing a surface A in the fluid. Then the line integral of the velocity u around C is: (7.1) Note that dr is on C. When C is a material circuit, meaning that it is composed of the same fluid elements at all time and therefore it moves with the fluid. Then we can d/dt eqn.(7.1) and bring the d/dt inside the integral: = (7.2) Since u.du = d|u2|/2, the integral . Now we substitute the momentum equation: du/dt + 2 u = p/ + g + F (7.3) into (7.2) which then gives: (7.4) Meaning of Each Term: 1. --- consider a u that is perpendicular to . (7.5a) 2. (7.5b) 3. --- this is the frictional term. (7.5c) Page 118 In application to the Island Rule, the fluid is single-layer and homogeneous, so the friction F is given by: F = (o b)/H, (7.6) where H is water depth and o and b are wind and bottom stresses respectively. The velocity u is then a depth-averaged velocity (u, v). Since the motion is barotropic, the p and surfaces coincide, then the p = 0, and the contribution from in (7.5b) vanishes. The x and y momentum equations can be written as: u/t + (+f) k u = p/ + |u2|/2) + (o b)/H. (7.7) where k is the unit upward vector and = v/x u/y is the z-component of the relative vorticity. This equation looks somewhat different from the usual form. But by expanding (k u) and (|u2|/2): (k u) = (iv + ju) = i v(v/xu/y) + j u(v/xu/y) (7.8a) and (|u2|/2) = i [(u2+v2)/2]/x + j [(u2+v2)/2]/y = i [uu/x + vv/x] + j [uu/y + vv/y] (7.8b) and adding them (and including the f also), we have: (+f) (k u) + (|u2|/2) = i [uu/x + vu/y fv] + j [uv/x + vv/y + fu] (7.9) which (in (7.7)) gives the usual form. Integrate (7.7) now around a circuit CI that surrounds an island (as in 7.2 and 7.4): , on circuit CI. (7.10) Page 119 The circuit-integrals of other terms, (+f) k u and p/ + |u2|/2), vanish, first because the normal component of the velocity to the island is zero, and secondly because the circuit integral of (…) = 0. Consider that the island is Taiwan at longitude xI and latitudes ys to yn; the island is located between China at xw and America at xe. Then consider a new circuit CR in the Pacific Ocean between Taiwan and America: (xe, ys) (xe, yn) (xI, yn) (xI, ys) (xe, ys). The Bernoulli’s terms again integrate to zero, but the vortex terms are not zero along the two latitude lines (xe, yn) (xI, yn) and (xI, ys) (xe, ys): , on CR, (7.11) where n = unit normal to CR, and the circuit integral is positive cyclonic (hence the minus sign in front of fn). Now, (7.12) is just the Sverdrup transport + Kuroshio transport and is in fact the stream function on the island, and is equal to the transport between China and Taiwan = Taiwan Strait transport. To see this, consider the definition of /x = v, and integrate to calculate for from China to Taiwan then to America – try it. Adding the circuit integrals on CI and CR (7.10) and (7.11), we then have: , on circuit CI + CR. (7.13) which in steady state and assuming that b 0: , on circuit CI + CR. (7.14) which is the Island Rule. Page 120 3. Homogenization of 2 Fluid Layers by Wind L.-Y. Oey (lyo@princeton.edu) We wish to derive Eqn.2 of Oey et al. 2006 (Loop Current warming by hurricane Wilma, GRL 33, L08613, doi:10.1029/2006GL025873, 2006). The followings have variously been derived by several people: J. Turner (original?), John Simpson (also original?), Chris Garrett, Peter Rhines and perhaps others. Consider two fluid layers that are mixed into one layer by turbulence due to, e.g. wind, as shown schematically by the following figure: Figure 3-1 Left: two layers of different densities and - ( > 0) before the onset of mixing (by wind). Right: wind mixes the fluids and homogenizes the two layers to one layer of some intermediate density = am. The physical situation may be that in summer, warm fluid overlies cold fluid (Figure 3-1 left). A tropical cyclone (typhoon or hurricane) comes and mixes the two fluids into one layer (Figure 3-1 right). In this idealization, we assume that there is little agitation of the fluid below the lower layer of Figure 3-1 (left). Page 121 What is am, the final density of the fluid after mixing, in terms of the original densities of the two fluids before mixing? If h1 = h2 then clearly am will be a simple average of and , i.e. am = /2. In general, am is a weighted-average of and . In other words, by mass conservation (i.e. mass after mixing = mass before mixing), we have per unit horizontal (xy) area: am (h1 + h2) = .h2 + ().h1, (3-1a) Which shows that am is a weighted-average of and : am = [.h2 + ().h1]/[h1 + h2] (3-1b) Wind pumps energy into the fluid by doing work through turbulent mixing, and thereby increases the potential energy of the fluid. The rate (per unit area) at which the wind does work is simply equal to the frictional stress (the drag force per unit area acted by the ocean surface on the air above) = aCd|W|2 times the air parcel’s speed (assume that the wind speed >> ocean current speed), thus: (Wind Work)/Area = aCd|W|3 in J/(m2.s) or Watts/m2 (3-2) where J = Joules (= Newton meters). If the wind blows for a time period = , then the energy (per unit area) supplied by the wind to the water column, or the power dissipation by the wind, is: , in J/m2. (3-3) Here, the factor accounts for the fact that not all the wind work is used to produce mixing in the water column – some of it is used to create waves, sprays etc. In fact, only a relatively small fraction ~ 10% (or less) is used to produce mixing in the fluid below. However, for our purpose (see below), we need not worry about this. To see that the wind work raises the potential energy (PE; again per unit area) of the fluid, we need to compare the PE after mixing, PEAM, with the PE before mixing, PEBM. The PE is simply the work done in raising the fluid through the water column, i.e. ~ the integral with respect to z of the weight of the fluid parcel as it is raised through the water column: (3-4) Page 122 where the primes denote general or dummy variables. Evaluating the last two integrals in (3-4), we obtain for the potential energy of the fluid, per unit area, before mixing: PEBM = (g/2)[h22 + 2h1h2 + ( )h12 (3-5) Similarly, (3-6) where equation (3-1a) has been used to substitute for “am (h1 + h2)” in the final expression. The change in PE, PE = PEAM PEBM, is then given by: or, , in J/m2, (3-7) which is positive, confirming that the PE of the mixed fluid is indeed increased. The ratio, (3-8) which is equation (2) of Oey et al. (2006), is then a ratio between the amount of energy required to thoroughly mix two stratified fluid layers and the amount of energy from wind that is actually available to do the mixing. If is “large,” then either the wind power is not sufficiently strong (or the wind has not acted long enough) to thoroughly mix the fluid layers, or that the stratification between the two layers is too strong, or the layer thicknesses are too large for the given wind to mix the layers. The opposite is true if is “small,” wind is very strong, or stratification is weak, or layers are thin. Page 123 It is difficult to exactly say just what value of is to qualify it as “large” or “small.” (One of the difficulties is of course our very poor knowledge of the value of .) However, one can estimate , , h1 and h2, and calculate (and plot) the for particular situation, say as a time series, and attempt to correlate (visually is the simplest!) it with, say, the time-series of sea- surface temperature or better still the temperature a few meters below the surface. If you suspect a correlation, you can then try to estimate what the ‘critical” is for your particular problem. Exercise: Carry out the integrations in (3-4) to obtain (3-5). Page 124 Gill’s Chapter 7: Effects of Rotation The principle of conservation of angular momentum Why do global winds display bands of easterlies and westerlies (see Gill’s fig.2.2)? To answer, we use the principle of conservation of angular momentum which is valid when there is no friction. Let τx(φ) = eastward wind stress = force/area acting on the earth’s surface. This stress transfers a tourque = rate of transfer of angular momentum per unit area of the earth’s surface: torque = rE.cos(φ).τx(φ) Multiply this by an elemental area of zonal strip = 2πrE.cos(φ).rE.dφ, inegrate from south pole to north pole to cover the entire earth’s surface, then set the integral to zero because there can be no net torque (otherwise the earth’s spin may vary with time and it will be hard for us to plan a daily schedule): or, . (2.3.1) Therefore, since cos2(φ) ≥ 0, τx(φ) must change sign, i.e. must alternate between westerly and easterly (figure). The principle of conservation of angular momentum is used over and over again in GFD. Page 125 7.2 The Rossby Adjustment Problem Previously in chapter 5 we studied the outward propagation of an initial discontinuity in free surface when there was no rotation. By method of characteristics we were able to calculate that the discontinuity spreads at a propagation speed = (gH)1/2. We now examine what happens to the discontinuity if there is rotation. In that case, we have eqn.(7.2.4). Exercise 1a: Using hydrostatic equation ∂p/∂z = −ρg, show that p = patm + po(z) + ρgη, where po(z) = −ρgz is the equilibrium (static) pressure and patm is the atmospheric pressure. (Often we set patm = 0, or at least patm = 0. In numerical modeling of storm-generated currents, the patm may be specified from weather-chart and can be important). Exercise 1b: State all conditions and assumptions under which equations (7.2.1), (7.2.2) and (7.2.3) are valid. Exercise 2: Derive eqn. (7.2.4): ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + fHζ = 0 (7.2.4) For f = 0, this equation can be solved given initial and boundary conditions for η. For f ≠ 0, we need another equation that relates ζ and η. (Alternatively, of course, one can solve for (u,v,η) from the 3 Page 126 equations (7.2.1)−(7.2.3); but it is often easier to solve just 1 equation (7.2.4) for η when f=0, or derive another one when f≠0 to solve for η and ζ). Exercise 3: Derive eqn. (7.2.8), the conservation of (linearized) potential vorticity: ∂(ζ/f – η/H)/∂t = 0 (7.2.8) Thus for a homogeneous rotating fluid, the quantity Q’ is invariant with time: Q’ = (ζ/H – fη/H2) (7.2.9) Note that Q’ may be and in general is a function of the spatial coordinates (x,y). But at each point (x,y), it retains its initial value for ALL time: Q’(x,y,t) = Q’(x,y,0). (7.2.10) Exercise 4: The full PV = (f + ζ)/(H + η). Show that in linearized form (i.e. ignoring products such as (ζη) terms and higher), this becomes f/H + perturbation, where the perturbation is Q’. Eqns.(7.2.9) and (7.2.10) give: fHζ = fH2Q’(x,y,0) + f2η (N7.1a) or after dividing by f2H: ζ/f – η/H = HQ’(x,y,0)/f = ηosgn(x)/H (N7.1b) Substituting this into (7.2.4) then gives: ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = − fH2Q’(x,y,0) (N7.2) Consider now our initial condition of a discontinuous free surface: η(x,y,0) = −ηo.sgn(x), and u = v = 0, (7.2.11) Page 127 so that from (7.2.9) we have Q’(x,y,0) = fηosgn(x)/H2 which after substituting into (N7.2) gives: ∂2η/∂t2 – c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = − fH2Q’(x,y,0) = − f2ηosgn(x) (7.2.13) Exercise 5: Carefully follow the above steps, then derive the appropriate equations on the β-plane, i.e. for the case when f = fo + βy. Then look back notes discussed previously on Sverdrup transport etc, compare and discuss. 7.2.2 The Steady Solution: Geostrophic Flow We wish to know the steady state that the initial discontinuity in elevation (eqn. 7.2.11) may evolve into after a long time. It seems reasonable that we find the steady solution from the steady equations – i.e. from equations (7.2.1)−(7.2.3) with all the ∂/∂t terms set to zero − so let’s try that: fv = g∂η/∂x, fu = −g∂η/∂y (7.2.14) Exercise 6: Describe what you find after you try it, and are you then able to find the steady-state solution for (u,v,η)? Why yes (if you could find the solution) or why not (if you could not)? Exercise 7: Based on the conclusions of Exercise 6, look back pom-book’s chapter 2 and find out what we did when we encountered a similar dilemma with the Ekman problem. In Exercise 6, you should find that the geostrophic relation (7.2.14) gives a velocity field (u,v) that identically satisfies the (steady-state) continuity equation (7.2.3). Thus any η(x,y) would produce (u,v) (from 7.2.14) that satisfies all the steady-state equations. In other words, we cannot find a unique solution. In Exercise 7, you should find that for the steady Ekman problem of pom-book’s chapter 2, we relied on viscosity (i.e. Ekman layer) to overcome the difficulty and found the unique Ekman-layer solution. In the present case, we do not have viscosity. But recall that in the Taylor-Proudman column problem we were able to find the solution numerically, and we mentioned then that that was accomplished by solving the pom-model by a time-dependent integration. So now we can do the same, but we do it analytically. Actually, we have done so when we integrated (7.2.8) in time to obtain (7.2.10), and finally we obtained (N7.1b). Page 128 Exercise 8: In (N7.1b), approximate the ζ using the geostrophic velocity from (7.2.14), then show that: (g/f2)(∂2η/∂x2 + ∂2η/∂y2) – η/H = ηosgn(x)/H, so that: −c2(∂2η/∂x2 + ∂2η/∂y2) + f2η = −f2ηosgn(x) (N7.3) which is the steady-state form of (7.2.13). Equation (N7.3) governs the form of η, given some conditions on η for large |x|. Since the RHS (which is the ‘forcing’) is not a function of ‘y’ we may look for a solution η(x), function of ‘x’ only, so that: −c2d2η/dx2 + f2η = +f2ηo, for x < 0, and (N7.4a) −c2d2η/dx2 + f2η = −f2ηo, for x > 0 (N7.4b) The solution that is finite at large |x| is then: η/ηo = 1.exp[+x/(c/f)] + 1, for x < 0, and (N.7.5a) η/ηo = 2.exp[−x/(c/f)] − 1, for x > 0. (N.7.5b) Matching the two solutions at x =0 where the η(0) needs to be continuous, =0, we choose 1 = −1 = − 2, so that finally: η/ηo = −exp(+x/a) + 1, for x < 0 = +exp(−x/a) − 1, for x > 0, (7.2.22) where a = c/f = (gH)1/2/f (7.2.23) is the Rossby radius of deformation. From the geostrophic relation, we obtain: u = 0, and v = −(gηo/fa).exp(−|x|/a) (7.2.24) Page 129 The flow is not in the direction of the pressure gradient, but is perpendicular to it along line of constant η (i.e. contours of η), as shown in Gill’s fig.7.1. That geostrophic flow is along η-contours can be shown by noting that u. η = (g/f).[−∂η/∂y.∂η/∂x + ∂η/∂x.∂η/∂y] = 0, so that the velocity vector u and the direction normal to the contours η are perpendicular; i.e. u and η- contours are parallel. An alternative derivation: I mentioned in the class that the reason we could obtain a unique solution is that we retain the small ∂/∂t terms, which then leads to (7.2.10). We will now do this more systematically by first scaling the equations. Scaling reveals which terms are small and which are larger. Let U be the velocity scale, H the scale for η, L the scale for (x,y), and εL/U the scale for time, where ε = U/fL is the Rossby number. The reason we scale the time “t” in this way will become clear shortly. We then define nondimensional variables (denoted by subscript n) as follows: u = un.U, v = vn.U, x = xn.L, y = yn.L, t = tn/(εf), η = ηn.H (N7.6) By scaling, we anticipate that all non-dimensionalized variables are of O(1). We see that, since ε is expected to be small, tn ~ O(1) ~ ε.ft, and we therefore require that the time “t” varies very slowly over many inertial periods ~ 1/(εf). Then equation (7.2.1) becomes (dropping the subscript “n”): fUε (∂u/∂t) – fU.v = −(gH/L).∂η/∂x, or ε ∂u/∂t – v = −(R2/εL2).∂η/∂x, ε ∂u/∂t – v = −(1/ε) ∂η/∂x, (N7.7a) where R = c/f is the Rossby radius (Gill uses “a” for this), and we have chosen L = R. Similarly for eqn.(7.2.2), we have: ε ∂v/∂t + u = −(1/ε) ∂η/∂y. (N7.7b) The continuity equation (7.2.3) becomes: Page 130 ∂η/∂t + (∂u/∂x + ∂v/∂y) = 0. (N7.7c) We try to find the solution to (N7.7a-c) by using power series expansion in ε: u = uo + ε.u1 + ε2.u2 + …; v = vo + ε.v1 + ε2.v2 + …; where un & vn ~ O(1). (N7.8a,b) For η, since the flow is expected to be nearly geostrophic, (N7.7a,b) suggests that: η = ε.η1 + ε2.η2 + …; where ηn ~ O(1). (N7.8c) Substituting (N7.8) into (N7.7), we obtain for the zeroth-order equations: vo = ∂η1/∂x, uo = −∂η1/∂y, and ∂uo/∂x + ∂vo/∂y = 0 . (N7.9a,b,c) The leading (i.e. zeroth) order then satisfies the degenerate geostrophic balance and the corresponding velocity (uo, vo) cannot be uniquely determined because it satisfies the (leading order) continuity equation (N7.9c); i.e. we only have two equations for three unknowns (uo, vo, η1). We therefore examine the O(1) equations: ∂uo/∂t − v1 = −∂η2/∂x, (N7.10a) ∂vo/∂t + u1 = −∂η2/∂y, (N7.10b) ∂η1/∂t + ∂u1/∂x + ∂v1/∂y = 0. (N7.10c) Use (N7.10a,b) to eliminate u1 and v1 from (N7.10c): ∂[η1 – ζo]/∂t = 0, (N7.11) where ζo = ∂vo/∂x − ∂uo/∂y = 2η1 from (N7.9a,b). Equation (N7.11) is the PV-conservation equation (7.2.8) in non-dimensionalized form. Integrating (N7.11): ζo – η1 = Q(x,y,t) = Q(x,y,0), or Page 131 − 2η1 + η1 = −Q(x,y,0) (N7.12) which is the steady-state form of (7.2.13) (i.e. eqn.7.2.20 when ∂/∂y = 0) in non-dimensionalized form (check it!). We see clearly that through our systematic scaling and series expansion, we obtain the same steady-state equation (N7.12) which is not degenerate, and can be solved. What have we learnt from this systematic approach to solve large-scale, slowly-varying flows with rotation? 1. The scale for length is the Rossby radius, L ~ R = (gH)1/2/f; 2. The scale for time is 1/(εf) ~ many inertial periods if ε << 1, which is the case; 3. η*/H is small ~ O(ε); but our chosen velocity scale ‘U’ is such that u*/U ~ O(1); here the asterisks *’s denote dimensional quantities. 4. To a good approximation (i.e. to the leading order), the velocity (uo, vo) is geostrophic; 5. The vorticity ζo can be evaluated assuming geostrophy; but the equation that relates it’s temporal variation (i.e. ∂ζo/∂t) to the slow time-dependent variation in η*/H (i.e. ∂η1/∂t) can only be obtained by considering the O(ε) balance in the series expansion. This leads to the important PV- conservation equation (N7.12) that may be used to solve for η1, after which the (uo, vo) may be computed through geostrophy – equations (N7.9a,b). The above approach is general, and the resulting approximation is called quasi-geostrophic approximation (or QG-approx). “Quasi” means “almost like but not exactly so” – in our case, we used geostrophy for (u,v) but we did not use it when we wanted to relate η with ζ. Another important thing to learn is that we can in fact include nonlinear advection terms in the momentum equations and we will obtain a very similar equation as (N7.12) and geostrophy for (uo, vo) (see exercise 9). The thing is, we never had to directly deal with ∂u/∂t or ∂v/∂t, just ∂ζ/∂t. Exercise 9: Show that if nonlinear advection terms are included, we obtain instead of (N7.7a,b): ε Du/Dt – v = −(1/ε) ∂η/∂x, (N7.13a) ε Dv/Dt + u = −(1/ε) ∂η/∂y. (N7.13b) and if H is still kept constant (it is not hard to also include a non-constant but slowly- varying H, i.e. gentle slope ~ O(ε)) the continuity equation (7.2.3) becomes: Dη/Dt + (1 + η)(∂u/∂x + ∂v/∂y) = 0. (N7.13c) Page 132 Now show that by substituting the power series expansion (N7.8a,b,c) as before, we obtain (N7.9a,b,c) as the leading order as before, and then the O(ε) equations are still (N7.10a,b,c) as before. Therefore, the PV-equation (N7.12) remains unchanged. Page 133 Lecture 8. Boundary Waves Kelvin Wave: Consider a homogeneous fluid (or the equivalent normal-mode equation, or reduced gravity), the x and y momentum equations are (Gill’s (7.2.1) and (7.2.2)): u/t fv = g /x, (8.1a) v/t + fu = g /y. (8.1b) Also, the mass conservation equation (H = constant): /t + H u/x = 0. (8.1c) Let y = 0 (or y = constant) be the straight coast. Kelvin discovered a type of rotational wave that can propagate along the coast. He reasoned that, near the coast, v 0, so that (8.1a,b) becomes: u/t = g /x, (8.2a) fu = g /y, (8.2b) Thus the y-momentum equation (8.2b) is geostrophic, i.e. the along-shore velocity “u” is always in geostrophic balance with the cross-shore sea-level gradient. But the alongshore momentum (8.2a) is highly ageostrophic. In fact, the x-momentum “u” is driven by pressure gradient in the same direction (as the momentum) – as if the flow has no rotation. Thus alongshore momentum (8.2a) and mass conservation (8.1c) in every y = constant plane (i.e. parallel to the coast) are exactly the same as in non-rotating flow, and may be combined to give: t2 = c2 x2 (8.3) Noting that “y” is kept constant in (8.3), so its solution is: = F(x+ct,y) + G(x-ct,y), (8.4a) which can be substituted in (8.2a) to give u: u = (g/H)1/2 [F(x+ct,y) G(x-ct,y)]. (8.4b) Page 134 Substitute and u from (8.4) into the y-momentum equation (8.2b) now gives, since “F” and “G” are independent solutions: F/y = F/R, and G/y = G/R, (8.5a,b) where R = (gH)1/2/f is the Rossby radius of deformation. For f > 0 and y points seaward, only (8.5b) with the exponentially decaying solution “G” is allowed, so that: G = A(x-ct).exp(y/R) (8.6a) which represents a wave propagating in the positive x-direction – see fig.8.1. If f < 0 and y again points seaward, then (8.5a) with the “F(x+ct,y)” solution must be chosen: F = A(x+ct).exp(+y/R) (8.6b) Figure 8.1 Kelvin wave. Substitute (8.6a) into (8.4) then gives for and u: = A(x-ct).exp(y/R) (8.7a) u = (g/H)1/2 A(x-ct).exp(y/R) = (g/c). (8.7b) where c = fR. If we think R as being always > 0, then c > 0 in northern hemisphere, and c < 0 in southern hemisphere. Then (8.7a) and the second form of (8.7b) are valid for both f > 0 and f < 0. In summary: 1. Kelvin wave travels at the speed “c”, exactly the same as for non-rotating shallow-water waves; Page 135 2. The wave form “A” is the same on every y=constant plane (parallel to the coast), except that its amplitude decreases exponentially with distance from the coast; 3. Therefore the wave has significant amplitude only within a distance R c/f from the coast; 4. The velocity normal to the coast is zero: v = 0; 5. The velocity parallel to the coast, u, is in geostrophic balance with the sea level (show this from (8.7)). Therefore, from the second form of (8,7b), u is > 0 (< 0) in the northern (southern) hemisphere where coastal sea-level |y=0 is > 0. And the signs of u are reversed where the |y=0 is < 0. Things to discuss on Friday Apr/30/2010: 1. Kelvin wave energy flux.. 2. UE/t .. total atmos+ocean mass transport = 0; 3. Why T/z << u/z .. mixed layer for temperature? 4. Why barotropic geostrophic flows must be parabathic? 5. Read Gill 10.9 (Storm Surge) for next week. Coastal Setup (or Setdown)by Wind: Consider an infinite coast bounding an ocean of uniform depth H (Fig.8.1) on a rotating earth with constant f. At t = 0, a uniform kinematic wind stress ox blows in x-direction. Since coast is infinite and wind is uniform, we consider the case /x = 0. Then instead of (8.1), we have:7 u/t fv = (oxbx)/H (8.8a) v/t + fu = g/y by/H (8.8b) /t + Hv/y = 0 (8.8c) where (bx,by) is bottom frictional stress assumed = 0. Wind blows in x-direction, so there is Ekman transport, hence v 0. We can (8.8b)/t, then use (8.8a) for u/t, and use (8.8c)/y for 2/yt: 2v/t2 + f2v c22v/y2 = fox/H (8.9) After a time longer than inertial period, t >> 1/f, transient inertial oscillations die away, the 2/t2 term is small, i.e. steady state.8 The solution consists of exponentials, and we get rid of the growing, e+ part and also make the solution satisfy v = 0 at y = 0: 7 Derive this by integrating the linearized equations with vertical shear stress terms, and assuming = constant. Page 136 v = ox/(fH) [1 ey/R]. (8.10) Divergence (u/x + v/y) of water near the coast is (since /x = 0), from (8.10): v/y = ox/(RfH)ey/R (8.11) If ox > 0, we then have convergence i.e. v/y < 0 and sea-level rises from (8.8c). If ox < 0, we have divergence i.e. v/y > 0 and sea-level falls. From (8.8c), the rise and fall is: = ox/(Rf)ey/R.t (8.12) Therefore, although v is steady, is not, and linearly increases with time. Since from (8.8b) the “u” is in geostrophic balance with the sea-level, it too linearly increases with time: u = +(ox/H).ey/R.t (8.13) The PV-equation is: ( f/H) = 0 (8.14) so that the (relative) vorticity also increases indefinitely. In practice, as “u” increases, friction eventually becomes large on the RHS of (8.8a).9 Kelvin wave energy flux: For a fixed y, Kelvin wave behaves as in flow without the earth’s rotation – see eqns.(8.3) & (8.4). So we derive first the energetic for f = 0, and then apply it to Kelvin wave. The equations are: /t + (Du)/x + (Dv)/y = Q (8.15a) u/t fv = g/x ru (8.15b) 8 To see this, define non-dimensional variables: tn = ft, yn = y/R, vn = v/U, and onx = ox/(fUH), where << 1, then vn/tn2 + vn 2vn/yn2 = onx, and the first term can be neglected. 2 2 9 Effects of friction are seen by considering very near-coast region where “v” is small. Then (8.8a) gives ox bx, etc. Page 137 v/t + fu = g/y rv (8.15c) Here, D = H+, Q = source (in m s-1) and r = constant bottom friction coefficient. Multiply (8.15a) by g: (g2/2)/t + g[(Du)/x + (Dv)/y] = gQ (8.16) Note that g2/2 = PEA (8.17a) is potential energy per unit (horizontal) area, since it is = gz dz’, where the integral is from z’ = 0 to . For our case, = constant since the fluid is assumed to be homogeneous. Equation (8.16) is therefore the potential energy equation. Multiply (8.15b) by Du and (8.15c) by Dv, and then sum the resulting 2 equations, we obtain, with KE = (u2+v2)/2: D.(KE)/t = g(uD/x+vD/y) rD(u2+v2) = g{ [(uD)/x+(vD)/y] [(Du)/x+(Dv)/y] } 2rDKE which after using (8.16) and (8.17a) and also approximating D.(KE)/t = (D.KE)/t KE./t (D.KE)/t, gives: (KEA + PEA)/t = g[(uD)/x+(vD)/y] + gQ 2rDKE (8.18) where KEA = D.KE = D. (u2+v2)/2 (8.17b) can be seen to be the kinetic energy per unit (horizontal) area of the column of water (of depth “D”). For Kelvin wave, the solution is (8.7a,b). Let’s consider the special case when the amplitude is of a simple cosine form: A(x-ct) = cos(x-ct) = cos( ), = x-ct. (8.19) So (8.7a,b) become: Page 138 = o cos( ) exp(y/R) (8.20a) u = o (g/H)1/2 cos( ).exp(y/R) = (g/c). (8.20b) Let <.> denotes the mean over all “ ”, so that since cos2 .d /(2 where the integration is from 0 to 2, we have from (8.17a,b) that the mean PE or KE per unit x-length (PEx or KEx) is: PEx = g<2>/2.dy = go2R/8 = H<u2>/2.dy = KEx (8.21) where the integral over ‘y’ …dy is taken from y= 0 (coast) to y = . Therefore, the mean kinetic and potential energies per unit x-length are equal. This is the same conclusion reached for f = 0 step- problem in Gill’s Chapter 5 that the energy of motion is partitioned equally between KE and PE. The mean energy flux at any “x” along the coast is given by (see the 1st term on the RHS of 8.18): Efluxx = gH<u>dy = g2Ho2/(4|f|) (8.22) For Kelvin wave propagating along the coast, the 2nd term on the RHS of (8.18), (vD)/y, is small compared to the 1st term, so that the energy flux is Efluxx given by (8.22). Now, if friction can be neglected, we have by integrating the total energy equation (8.18) over an x-section of the coast where Kelvin wave passes, say from x1 to x2, that: (<KEA> +< PEA>)dxdy}/t = 0 = gH[<u>|x1<u>|x2].dy (8.23) The LHS is = 0 since the total energy of Kelvin wave over (x1, x2) is constant. Therefore, <u>|x1 = <u>|x2 and the mean energy flux Efluxx is constant. Therefore, if at any “x2” the water depth H2 is shallower than the section previously x = x1 < x2 and H2 < H1, then the amplitude at x2 must increase: [o]|2/[o]|1 = (|f2|/|f1|).(H1/H2)1/2 > 1, if |f| does not change much. (8.24) Therefore, flooding and drying can occur when a Kelvin wave (generated say by a distant storm) propagates into a shallow coastal region; made worse if propagation is poleward (so that |f2| > |f1|). Summary: 1. Step-function alongshore wind produces inertial oscillations that radiate away; 2. The wind then drives steady Ekman transport towards or away from the coast; 3. The convergence or divergence causes coastal sea-level to rise or fall linearly with time; Page 139 4. The alongshore velocity “u” also increases linearly with time, in geostrohpic balance with ; 5. Both “u” and “” decays exponentially away from the coast over a distance ~ R. 6. Eventually bottom friction is important, and balances the wind. Problem 8.1. Modify the plume-on-shelf problem in http://pom-training.pbworks.com/ (File: runpom08_01_assignment5b_horcon050_nadv1 in ftp://aden.princeton.edu/pub/lyo/pom_gfdex/wmo09training/anIntroCourseNumOceanExpsUsingPOM ). Domain is a square or rectangular basin, coastal walls on all sides. Then ramp (in time) an alongshore wind in a finite coastal strip 50km alongshore and 10km cross-shore say at the center portion of the western boundary. Examine the resulting Kelvin wave. a. No rivers, no tides, and ocean is initially homogeneous, bottom is flat; b. Same as ‘a’ but ocean has shelf/slope all around; c. Same as ‘a” but ocean is vertically stratified; d. Same as ‘b” but ocean is vertically stratified; e. Same as ‘a’ but a circular cyclonic wind system (typhoon) propagating from open ocean (east; i=im) to land (west; i=1) along the center latitude (j=jm/2) of the domain at a speed of 5 m/s. Problem 8.2. Design an experiment to verify (8.24). Page 140 9. Dynamics of Long Rossby Wave (Propagation) Forced by Wind Stress Curl Consider the linearized reduced-gravity equations: h/t + (hu)/x + (hv)/y = Q h (9.1a) u/t fv = g’h/x + ox/H (9.1b) v/t + fu = g’h/y + oy/H (9.1c) where (x, y, t) = (west-east, south-north, time), h = upper-layer depth, H = mean upper-layer depth, Q = a source term, = Newtonian cooling coefficient, (u,v) = horizontal velocity, g’ = g/o is reduced gravity, f = Coriolis parameter, and (ox, oy) = wind stress vector. For slow motions with time scales much larger than 2/f ~ days, the /t terms underlined in (9.1b,c) are small compared to the Coriolis terms, so may be neglected. Then the motion is nearly geostrophic: vg (g’/f) h/x (9.2a) ug (g’/f) h/y (9.2b) where the subscript “g” denotes “geostrophic.” From (9.2), we have ugh/x + vgh/y = 0. (9.3) Equation (9.1a) is: h/t + (uh/x + vh/y) + h(u/x + v/y) = Q h (9.4) The term “(uh/x + vh/y)” is nearly zero by (9.3), while the “h(u/x + v/y)” term is obtained by taking the curl of the momentum equations (9.1b,c), i.e. by “(9.1c)/x (9.1b)/y”, which gives: h(u/x + v/y) H(u/x + v/y) = Hv/f + ×o/f Hvg/f + ×o/f = [(Hg’)/f2]h/x + ×o/f = R2 h/x + ×o/f (9.5) where = df/dy and R = (Hg’)1/2/f is the baroclinic Rossby radius. Using (9.5), equation (9.4) then becomes, with CR = R2, the long Rossby wave speed: h/t CR h/x = ×o/f + Q h. (9.6) Page 141 Eqn.(9.6) describes the westward propagation of (long) Rossby wave which without forcing (RHS=0) has the solution of the form h = G(x+CRt) (9.7) But the wave is also being continuously forced by some distribution of sources and sinks on the RHS given by (i) ×o/f = Ekman pumping (note water is added to upper-layer if this is <0), and (ii) Q representing some waters entering (>0 i.e. entrainment) or leaving (<0 i.e. detrainment) the upper layer, Note that mathematically the source/sink terms ×o/f and Q are indistinguishable from each other. At the same time the wave is also damped by (iii) h = Newtonian cooling or heating, ‘cooling’ or ‘heating’ because “h” can physically be imagined as representing the upper-layer warm water. It is “damping” because of the minus sign, so that when h, this term tends to produce h/t, and vice versa. Note that also the wind stress acts on a thin surface Ekman layer (~E), and its effects on the upper- layer motion is mainly through Ekman pumping as a result of differential Ekman transports inside E (as we will shortly see). Page 142 Lecture 10. El Niño, La Niña and the Southern Oscillation 0. Historical Introduction: El Niño is the name given to the irregular development of warm ocean surface waters along the coast of Ecuador and Peru. It develops when there is change in the circulation of the atmosphere across the tropical pacific. Globally, the development of El Niño is also associated with a number of other climatic changes in other parts of the world. To explain El Niño, it is necessary to explain how the ocean adjusts to the changes in the surface winds. Responses (wind and ocean currents) near the equator are rapid. Recall Gill’s adjustment with step problem – non-rotational gravity waves (gH)1/2, and also Rossby waves R2 (for example) are faster because R = Ci/f is larger. So ocean responds fast to winds, and air-sea coupling is particularly efficient: slight change in one leads to immediate changes in the other and vice versa. Gilbert Walker (1923~1937): knew that interannual pressure fluctuations over the Indian Ocean and eastern tropical Pacific are out of phase: “When pressure is high in the (eastern) Pacific Ocean it tends to be low in the Indian Ocean from Africa to Australia.” He called this irregular fluctuation the Southern Oscillation (SO). The SO is correlated with major changes in the rainfall patterns and wind fields over the tropical Pacific and Indian Oceans, and with temperature fluctuations in southeastern Africa, southwestern Canada, and southeastern USA. El Niño is that phase of the SO when the trade winds are weak and when pressure is low over the eastern and high over the western tropical Pacific. Page 143 La Niña is apposite when sea surface temperatures in the central and eastern tropical Pacific are unusually low and when the trade winds are very intense. http://www.physicalgeography.net/fundamentals/7z.html Page 144 CHAPTER 2???: Internal Waves and Rays in the Ocean In chapter 1 we studied how wind generates ocean currents and sea- level variations assuming homogeneous fluids ( = constant). We also focused on motions having horizontal scales that were much larger than the vertical scales; in other words, we assumed also that the hydrostatic equation was satisfied. These simplifications allowed us to examine rudimentary boundary-layer dynamics and to understand their profound effects on the large-scale (interior) oceanic flows. In this chapter, we remove these restrictions. We examine internal fluid motions in a stratified fluid ( is not constant, and generally a function of the vertical coordinate z). We do not assume hydrostasy, but we restrict ourselves to motions which are of small amplitudes that the equations are linear, and which are also of small horizontal scales and short time scales that the Coriolis force can be neglected. We will study rays – which are paths (generally in three dimensions) along which internal wave energy propagate. We then attempt to numerically model these waves (and rays) in two-dimensional vertical- section plane (xz-plane, say). We conclude the chapter by showing how the study of internal-wave rays may be applied to an entirely different physical system in the ocean: that of the study of ray propagation of hydrostatic and rotational topographic Rossby waves. We show that under certain conditions, the two systems are mathematically equivalent, and demonstrate a surprising physical consequence of rotation. 2-1: The Spring-Mass System - Generalized Inertia & Stiffness An oscillatory system such as the wave motion can be viewed as a spring-mass system: (2-1) Page 145 where m is the mass, is the spring’s stiffness, and z is distance by which the spring has been pulled or compressed (z = 0 is where the spring is in its natural position with zero tension). The solution to (2-1) consists of cosines and sines and the frequency squared is: 2 = /m (2-2) Equation (2-1) can be rewritten as (2-3) Now, the potential energy (PE) is = z2/2 while the kinetic energy (KE) is = , so that (2-3) is simply a statement that the total energy is constant in a system without dissipation. Equations (2-2) and (2-3) suggest that: (2-4) We will now use this idea to derive the frequency of oscillations for surface- interfacial waves, then internal waves also. 2-2: Surface & Interfacial Waves Consider then a thin (horizontal) interface that separates a heavier fluid of uniform density 2 below from a lighter fluid of uniform density 1 above. We assume that that the interface is far from the free surface above and also from the ocean bottom below, so that the free surface z = zT and ocean bottom z = h do not interfere too much with the interfacial motion; Page 146 here z = 0 at the interface. Let the interface be disturbed by a (small) displacement . The PE of the lower layer per unit area is: (2-5a) Thus the lower-layer PE is increased by 2g2/2 from its undisturbed state (i.e. when = 0). Similarly, the PE of the upper layer per unit area is: (2-5b) and in this case the upper layer has its PE decreased, 1g2/2 from its undisturbed state. Therefore the total PE per unit area due to is: (2-6) and the generalized stiffness, the coefficient of 2/2, is seen to be: generalized stiffness = g(2 1) (2-7) The disturbed interface displacement produces currents and therefore kinetic energy. To derive an expression for it, we need to examine surface (or interfacial) waves. Free-Surface Waves on Water with Constant Density: Page 147 Consider then a homogeneous fluid with a free surface (say we think of the lower layer only with uniform density 2, and in the followings think of “o” as “2”). Equations (1-24a,b,c) written in vector form and dropping the nonlinear, Coriolis and viscous terms become: o u/t = p’ (2-8) where p’ is the pressure in which a hydrostatic (undisturbed) pressure distribution po(z) (integrating equation 1-27 from z to 0; see equation (1- 32b)): po(z) = patmos ogz (2-9) has been removed (from the total pressure p): p’ = p patmos + ogz (2-10) and the that multiplies the acceleration on the LHS of (2-8) has been approximated by o. Taking the curl of (2-8) shows then that the vorticity field u is independent of time: ( u)/t = 0 (2-11) Therefore, the rotational part of the velocity field induced by this stationary vorticity field is steady, hence its corresponding pressure is zero by (2-8), and it cannot disturb . The remaining part of u is irrotational (i.e. it has Page 148 zero vorticity) and therefore can be written as (for simplicity we still call it u, and ignore the rotational part): u= (2-12) the gradient of a scalar function, so called the velocity potential. Only this part exhibits the fluctuations associated with wave propagation. From (2- 12), we see that, since water is incompressible (equation 1-25), we have: 2 =0 (2-13) This Laplace equation cannot by itself support wave motions if the boundaries are steady. But it can if the free surface is undulating in time. From (2-10), we have, setting p = patmos at z = : p’ = og, at z = 0 (2-14) the last approximation is because has been assumed to be small compared to the total water depth. Using (2-12) and (2-8) (giving o /t = p’), equation (2-14) gives: /t = g, at z = 0 (2-15) At the free surface, equation (1-31a) also holds (i.e. the free surface is a material surface), thus: /t = w = /z, at z = 0 (2-16) Equations (2-15) and (2-16) then give: Page 149 2 /t2 = g /z, at z = 0 (2-17) Equation (2-13) describes the motion in the fluid interior away from boundaries, and with equation (2-17) the system is almost complete. We assume that side boundaries are sufficiently far (in fact they are infinitely far) that waves propagate horizontally unimpeded. Besides the surface condition (2-17) we then also have the condition at the ocean bottom (1-31b) which for a flat bottom is simply that the vertical velocity is zero: /z = 0, at z = h (2-18) We will, however, relax the matter even further: we will assume that the bottom is also sufficiently deep that a solution (satisfying 2-13 and 2-17) that decays sufficiently rapid with depth will be adequate. Consider horizontal propagation in the x-direction. The Laplace equation admits a separable solution of the form: = (z) exp[i(tkx)] (2-19) where is frequency and k wave number. Substituting into (2-13) and solving for (z): (z) = o ekz (2-20) This solution decays with depth z ~ . Substituting (2-20) into (2-17) now gives the dispersion relation connecting and k: Page 150 2 = gk (2-21) The solution (2-20) allows us now to compute the total kinetic energy: (2-22) where equation (2-13) and the divergence theorem are used to obtain the second integral over the bounding surface S consisting of the surface and the bottom only. The bottom contribution vanishes, either because of (2-18) or because the solution decays exponentially with depth according to (2-20). At the surface, /n = /z = k-1(/t)2 after using (2-19), (2-20) and then (2-16). Equation (2-22) then gives for the kinetic energy per unit area: KE = (o/2) k-1(/t)2 (2-23) Comparing with equation (2-4), then, for this free-surface system, the generalized inertia is: generalized inertia = o k-1 (2-24) We also have, from (2-6), PE = go2/2, so that generalized stiffness = og (2-25) Equation (2-4) then gives 2 = gk in agreement with (2-21). Dispersion Relation for Interfacial Waves: Page 151 We saw from (2-6) that the disturbance gives rise to the PE available for wave motion. The same disturbance gives rise to kinetic energy given by (2-23) with the o replaced by the lower-layer density “2,” and similarly in the upper layer with o replaced by “1.” The total KE = [(1 + 2)/2] k-1(/t)2 (2-26) Equation (2-4) then gives, 2 = gk(2 1)/ (2 + 1) (2-27) Comparing with the frequency of the wave on a free surface, equation (2- 21), we see that the interfacial wave frequency (and hence also the corresponding phase speed = /k) is smaller, by a factor proportional to the square root of the density difference between the two layers. Note that (2- 27) leads to (2-21) if 2 >> 1. 2-3: Internal Waves Extensions of the above ideas to the case when the density o(z) is continuous require an exposition to the idea of a vertically stable stratified fluid when the density decreases with height. If a fluid parcel at z is moved up by a small amount to z + , the density of the surrounding fluid is now decreased: o + (o/z), with o/z < 0, > 0. (2-28) Page 152 Through the hydrostatic equation (1-27), po(z) = o(z)g, pressure of the surrounding fluid is also reduced to: po(z) p = po(z) o(z)g. (2-29) This pressure drop pressed upon the fluid parcel results in its density being reduced. Assuming isentropic process, the parcel’s density is reduced to: o(z) o(z)g/co(z)2, (2-30) where co(z)2 = (p/)s is the square of the speed of sound, and the subscript “s” denotes constant entropy. The parcel’s excess weight (per unit volume) over the surrounding buoyancy is then g[(2-30) (2-28)], or, g[o(z)g/co(z)2 + (o/z)] = oN2(z) (2-31) where N2 = [g2/co(z)2 + g(o/z)/o] (2-32) In order that there is a restoring gravitational force on the fluid parcel, this excess weight must be positive, which means that N is real and also that the vertical density gradient must be (from 2-31): o/z < o(z)g/co(z)2 (2-33) This shows that for stable stratification, it is not sufficient to merely require that o/z < 0. Page 153 For N2 > 0, the restoring force (2-31) is like the force exerted by the spring in the spring-mass system described above. The is also similar to the perturbed interfacial displacement described previously, but since the fluid is continuously stratified, the may be viewed as a generalized vertical coordinate for each fluid parcel. Then the potential energy per unit volume that results from is: (2-34) and generalized stiffness = o(z)N2(z) (2-35) What is the corresponding kinetic energy per unit volume? Equation (2-23) gives KE per unit area for a surface wave having a discontinuity in density, and a penetration depth scale ~ k-1. With continuous stratification, we do not know a priori what the penetration depth scale is (in fact there is none(!) as we shall see, there exists wave-propagation in the vertical also). It is more useful, therefore, to consider KE per unit volume which from (2-23) is: KE/volume = (o/2)(/t)2 (2-36) With surface or interfacial waves, the “” can only move vertically, and “the vertical” therefore represents the only generalized coordinate. The generalized inertia of these purely vertical oscillations is from (2-36) = o(z). But this must be the minimum inertia in the case of internal waves in continuously stratified fluid. This is because in a continuously stratified fluid, it is possible for a fluid parcel to oscillate along a line that makes an Page 154 angle T to the horizontal (or at = /2 T to the vertical); the corresponding displacement is , say (Figure 2-1). These types of oscillations clearly increase the kinetic energy which becomes ( o/2)( /t)2, hence also the generalized inertia (the coefficient of (/t)2/2), which is now increased to a value = o(z)/sin2(T). The purely vertical oscillations, for which T = /2, therefore represent a state of minimum inertia, hence maximum possible frequency. Figure 2-1. Parcel fluid oscillations (displacement = , say) inclined at an angle T to the horizontal, or at = /2 T to the vertical. Thus = /sin(T). In general then, (2-36) is replaced by: KE/volume = (o/2)(/t)2/sin2(T) (2-37) Page 155 and equation (2-4) gives, with (2-34) and (2-37): 2 = (generalized stiffness)/(generalized inertia) = N2 sin2 (2-38) or, = N.sin(T) = N.cos(), = /2 T (2-39) Figure 2-2 illustrates schematically how an IW excited in a thermocline may be trapped within it. It makes use of the idea expressed in (2-39) that the maximum possible frequency of oscillations cannot exceed the local N. Figure 2-2. Left: a typical N-profile in the ocean’s thermocline. Right: an internal wave (IW) with = N1 excited within the thermocline cannot escape to depths (beyond upper and lower boundaries) where the local N < N1. Exercises: 2-3A: Find out what the typical values of N are in the ocean. Then estimate typical (minimum) IW periods and compare them with the Coriolis periods. Page 156 2-3B: Re-derive (2-33) by letting < 0 instead of positive. 2-3C: Study the following alternative derivation of = N sin(T) from the equations of motion. Then close the book, and re-derive yourself. Derivation of (2-39) Using the Equations of Motion: In equation (1-24), we drop rotation, viscosity and the nonlinear advective terms, and replace the on the LHS (i.e. the acceleration terms) by o (this is the first part of the so called Boussinesq approximation), to obtain in vector form the following momentum equation: o u/t = p + g (2-40) where g = (0, 0, g). Subtract from this the hydrostatic relation (see equations 1-27, 1-28 and 1-29): po = og (2-41) we obtain: o u/t = p’ + ’g (2-42) The conservation of mass (i.e. continuity equation) is (see e.g. “Fluid Dynamics” by G.K. Batchelor, 1966): /t + (u) = 0 (2-43) Page 157 Now write /t + (u) = ’/t + [(o+’)u] and use the second part of the Boussinesq approximation that density fluctuations ’ (which of course affects motions, from equation 2-42) in the conservation of mass are neglected. In other words, the fluid motions of concern evolve sufficiently slowly that the rate of change of ’ as well as its spatial gradients do not affect the mass balance very much. Therefore, (ou) = 0 (2-44) Also, since ’ is the excess density of a fluid parcel over its surrounding o(z) (see equation 1-29), we have from the derivation of (2-31) that ’ = oN2/g (2-45) for the (excess) density when the fluid parcel is displaced a small distance from its initial position; note that (2-45) is valid for positive or negative (exercise 2-3B). We now manipulate (2-42), (2-44) and (2-45) to derive an equation solely in terms of one variable which we will choose to be ow = upward component of mass flux. Take the divergence of (2-42) and use (2-44), 2 p’ = g ’/z (2-46) then eliminate p’ from the z-component of (2-42), which is: Page 158 (ow)/t + p’/z = ’g (2-47) we obtain, 2 (ow)/t = g2’/z2 g 2 ’ = g(2’/x2 + 2’/y2) (2-48) The ’ on the RHS is now eliminated by first taking the time derivative, (2- 48)/t, then noting from (2-45) that g’/t = o(/t)N2 = owN2, (2-49) to finally obtain, 2 (2q/t2) = N2 h 2 q (2-50) where q = ow, and h 2 = (2/x2+2/y2) is the horizontal Laplacian operator. Incidentally, equation (2-49) written as b/t = wN2, b = g’/o = buoyancy, (2-51) may be recognized as being derivable from the density equation D/Dt = 0 (2-52) after linearization about the base (rest) state o(z) etc. (e.g. Gill, 1982). The Dispersion Relation: Page 159 For constant N, we may look for a (separable) wave-like solution of equation (2-50): q = ow = qo exp[i(tkxlymz)], (2-53) where k = (k, l, m) is the wavenumber vector. Substituting this into (2-50) gives, 2 = N2(k2+l2)/(k2+l2+m2) (2-54a) or, = N(k2+l2)1/2/(k2+l2+m2)1/2 (2-54b) which shows, in agreement with (2-38) or (2-39), that the maximum possible frequency of oscillations is N. Page 160 Figure 2-3. Incompressibility means that the fluid parcel oscillation with displacement vector and the internal wave vector k are perpendicular; assuming a plane-wave solution. (See also Figure 2-1). In fact, equations (2-54b) and (2-39) are equivalent (apart from the fact that we did not explicitly assume constant N in 2-39; see however, later). This may be proven as follows. Since the fluid is incompressible, we have from (2-44) that for a plane-wave solution such as (2-53), ku=0 (2-55) meaning that the wavenumber vector k and the fluid velocity u are perpendicular to each other. Since in Figure 2-1, the u is in the same direction (and plane) as the oscillation displacement vector , we have then k and are also perpendicular as shown in Figure 2-3. Equation (2-39) then follows by noting that the angle T made to the horizontal by the fluid parcel’s displacement vector during the oscillations is also the angle spanned by the wavenumber vector k with the vertical axis (Figure 2-3). Expressions for p’ and u in terms of qo: For completeness, we give here the corresponding formulae for p’ and u. For plane wave, we have (c.f. equation 2-53): (p’, u) = (po, uo).exp[i(tkxlymz)]. (2-56) Substitute the p’ into /t[(2-46)] and using (2-49) and (2-53), we obtain, p’ = m(k2+l2)-1qoexp[i(tkxlymz)]. (2-57) Page 161 Substitute (u, v) into the (x, y)-component of (2-42) and using (2-57) gives ou = [km(k2+l2)-1, lm(k2+l2)-1, 1]qoexp[i(tkxlymz)] (2-58) where the ow part just repeats (2-53). A Summary of the Governing Equations for IW’s: Before deriving the corresponding dispersion relation for TRW’s in the next section, it is useful to summarize here the governing equations for IW’s so that one may better see the differences with the TRW equations. In IW’s the momentum equations are (2-42): o u/t = p’/x (2-59a) o v/t = p’/y (2-59b) o w/t = p’/z ’g (2-59c) The conservation of mass (i.e. the “continuity equation”) is /t + (u) = 0 (2-60a) or D/Dt + u = 0 (2-60b) which we used for IW’s in the form of equation (2-44) keeping the o inside the divergence operator ( ) mainly for conciseness so we obtained an Page 162 expression for the vertical mass flux ow. However, for incompressible fluid, it may be shown from (2-60b) that the first term is of O(1/co2) so it is very small and therefore u=0 (2-61) which is the impressibility condition stated previously in equation (1-25). The important point is that in deriving the (oceanic) IW equations we have assumed that the fluid is incompressible. Finally, the density equation is (2-51), repeated here: b/t = wN2, b = g’/o = buoyancy. (2-62) Equations (2-59a,b,c), (2-61) and (2-62) give five equations for the five unknowns (u, v, w, ’, p’). 2-4: Topographic Rossby Waves In this section, we show the curious equivalence for the two- dimensional propagation of internal waves studied in the previous sections with the propagation of hydrostatic bottom-trapped topographic Rossby waves in the horizontal xy-plane. In other words, we show that one can study the propagation of non-hydrostatic internal waves using a hydrostatic model, or vice versa! We derive the corresponding dispersion relation for TRW’s. We consider fluid motions over a sloping bottom. We will assume that the scales are large that the earth’s rotation and hence also the Coriolis parameter f (Chapter 1) are important, but not too large so that we may ignore the variations of f with latitudes, i.e. f/y = 0. The scaling for these large-scale and relatively slow (compared to the earth’s rotation) flows have been previously detailed in chapter 1 (see, e.g. equation 1-38). It is more instructive for the present purpose of pointing out the similarities and Page 163 differences between TRW’s and IW’s to formulate by modifying equations (2-59), (2-61) and (2-62). In TRW’s we also assume that the fluid is incompressible so that equation (2-61) holds. The density equation (2-62) also remains the same (of course!). The x and y momentum equations, however, are modified by adding the Coriolis terms, and the z-momentum equation is reduced to the hydrostatic equation (even though the fluid is in motion!); thus equations (2-59a,b,c) become: u/t fv = P’/x (2-63a) v/t + fu = P’/y (2-63b) P’/z = b (2-63c) u/x + v/y + w/z = 0 (2-63d) b/t = wN2 (2-63e) where P’ = p’/o, the buoyancy b is used in place of (g’/o) in the hydrostatic approximation (2-63c), and the incompressibility and density equations are repeated from (2-61) and (2-62) respectively. We assume fluid bounded only at top and bottom but unbounded in the horizontal x and y directions. The boundary conditions are: w = 0 at z = 0 (2-64a) w = uh/x vh/y (2-64b) Take the curl of (2-63a,b) gives: /t + f (u/x + v/y) = 0 (2-65) Taking the divergence of (2-63a,b) gives: Page 164 /t (u/x + v/y) f = h 2 P’ (2-66) and this becomes after eliminating the horizontal divergence term from (2- 65): 2/t2 + f2 f 2 h P’ =0 (2-67a) or, assuming that the time scales of motions are much larger than the earth’s rotation period, f h 2 P’ = 0 (2-67b) We now (/z) then (/t) equation (2-63c), and then use (2-63e) to substitute for (2b/zt) to obtain: (/t)(2P’/z2) = N2w/z (2-68) Use now the incompressibility equation (2-63d) to replace w/z by the horizontal divergence, and then use (2-65) to obtain: (/t)(2P’/z2) = (N2/f) /t (2-69) Eliminate by using (2-67b) to obtain, finally, (/t)[ h 2 P’ + (f2/N2) 2P’/z2] = 0 (2-70) The boundary conditions are written in terms of P’ by using (2-63c) and (2- 63e): Page 165 2P’/zt = 0 at z = 0 (2-71a) 2P’/zt = (N2/f) (h/y P’/xh/x P’/y) at z = h (2-71b) where also fv = P’/x and fu = P’/y from (2-63a,b) have been used, again assuming that the time scales of motions are much longer than the earth’s rotation period. Substituting as before a separable plane-wave solution: P’ = (z).exp[i(tkxly)] (2-72) (2-70) gives: zz 2 = 0, 2 = (N2/f2)(k2 + l2) (2-73) The solution is: = Ae z + Be z (2-74) where A and B are constants. Equation (2-71) requires zt, which is: zt = i z = i [Ae z Be z] (2-75) Applying (2-71a) at z = 0 gives A = B = o, say. Equation (2-71b) becomes: z = [N2/(f)][hyk hxl] = [N2/(f)][Kh h]|z, at z = h Page 166 (2-76) Equaitons (2-75) and (2-76) then give: (/N) tanh( h) = sign(f) (nk h)|z (2-77) Now, tanh( h) = tanh(h/htrap), where -1 htrap = = |f|/(N|Kh|), where |Kh| = (k2 + l2)1/2 (2-78) so that if h/htrap > 1 (2-79) then tanh( h) 1 and (2-77) becomes: = sign(f) N | h| sin(T) (2-80) where T is the clockwise angle that the wave vector makes with h (Figure 2-4). Equation (2-80) is in the same form as the corresponding IW dispersion relation, equation (2-39). Page 167 Figure 2-4. A sketch that describes the relation between topography h, its gradient h, TRW wave number vector K and group velocity Cg for the northern hemisphere. The T is the clockwise angle the K-vector makes with h. Page 168 References Batchelor, G.K., 1967: Fluid Dynamics. Baumert et al., 2005: Gill, A., 1982: Atmospheric & Oceanic Dynamics Mellor, G.L., 2004: POM User Guide. Mellor G.L. and Yamada, T., 1982: Review of Turbulence Models. Pedlosky, J., 1979: Geophysical Fluid Dynamics. Roach, P., 1972: Computational Fluid Dynamics. Van Dyke, M., 1964: Perturbation Methods in Fluid Mechanics. Stewart, Robert H., 2000: Introduction to Physical oceanography. Columbia University Press, http://www.earthscape.org/r3/stewart/index.html.