The rigid lid streamfunction method was developed by Bryan (1969). This chapter presents the formulation of the rigid lid in a fashion to closely parallel the derivation given for the free surface in Chapter 7.
6.1 The barotropic streamfunction
With a rigid lid at the ocean surface, the surface heightη=0. Settingηto zero eliminates the very fast (order 200 m s−1in water with depth 4000−5000m) external mode gravity waves. As a result of making the rigid lid assumption, the vertically integrated horizontal velocity satisfies
∇h·U = uh(−H)· ∇hH+ Z 0
−H
dz∇h·uh
= uh(−H)· ∇hH− Z 0
−H
dz wz
= uh(−H)· ∇hH+w(−H)−w(0)
= −w(0). (6.1)
To reach this result, Leibnitz’s Rule (4.65), the continuity equation (4.3), and the bottom kine-matic boundary condition (4.24) were used. Taking the additional assumption
w(z=0)=0 (6.2)
renders
∇h·U=∇h·H(u,v) = 0. (6.3)
As a result, the external mode velocity can be expressed in terms of the external mode streamfunction
u = − 1
Ha
ψφ (6.4)
v = 1
Ha·cosφ
!
ψλ. (6.5)
As a vector equation, this relation takes the form
U=H(u,v)=zˆ∧ ∇hψ. (6.6)
59
6.2 Streamfunction and volume transport
The barotropic streamfunction is specified only to within a constant. As such, only differences are physically relevant. In particular, consider the vertically integrated advective transport between two points
where dl is the line element along any path connecting the points a and b, and ˆn is a unit vector pointing perpendicular to the path in a rightward direction when facing the direction of integration. As written,Tabhas units ofvolume/time, and so it represents a volume transport.
The definition of the barotropic streamfunction and Stokes’ Theorem renders Tab =
is a unit vector tangent to the integration path, pointing in the direction of integration from point a to point b. Therefore, the difference between the barotropic streamfunction at two points represents the vertically integrated volume transport between the two points. It is for this reason that the barotropic streamfunction is sometimes called thevolume transport stream-function. Note that Bryan (1969) defined the barotropic streamfunction with an extra factor of the Boussinesq densityρo, such thanψbryan =ρo ψmom. Hence, the barotropic streamfunction of Bryan has the dimensions ofmass/timerather than volume per time, and so it represents a mass transport streamfunction. Since MOM assumes a Boussinesq fluid, the difference is trivial.
6.3 Hydrostatic pressure with the rigid lid
The hydrostatic equationpz =−ρgcan be integrated from the surfacez=0 to some position z<0 to yield wherepa is the atmospheric pressure,pl is the surfacelid pressure, andpbis the hydrostatic pressure arising from the ocean’s density field. The surface lid pressure is the pressure which
6.4. THE BAROTROPIC VORTICITY EQUATION 61 would be exerted by an imaginary rigid lid placed on top of the ocean. An alternative interpre-tation is that it is the pressure exerted by undulations of a free surface. The latter interpreinterpre-tation is further disscussed in Section 7.1. In summary, the horizontal pressure gradient is given by
∇hp=∇h(pa+pl)+g Z 0
z
dz∇hρ. (6.11)
The horizontal pressure gradient at some depth z has a contribution from gradients in the atmospheric and lid pressure, gradients which act the same for all depths, and the vertical integral of horizontal gradients in the interior density. These latter gradients are due to baroclinicity in the density field, which prompts the often used namebaroclinic pressure gradient, and which motivates the “b” subscript. The vertically integrated horizontal pressure gradient is needed for the development of the barotropic vorticity equation. This gradient is given by
Z 0
Since the horizontal gradient of the atmospheric and lid pressures are independent of depth, Z 0
For cases where the ocean model is driven by a realistic atmospheric model, the atmospheric pressurepa may be determined. But since the lid pressurepl is actually an artifact of making the rigid lid approximation, it is generally not available. To address this fact, Bryan formed the barotropic vorticity equation in which case the lid and atmospheric pressures are eliminated.
Through inverting the vorticity equation, the streamfunction is obtained. As a result, the barotropic velocities can then be determined through the relation (6.6). As will be seen in Section 6.6, the lid and atmospheric pressures will likewise not be needed for determining the baroclinic velocity. This algorithm has been used for solving the rigid lid in MOM up until the formulation of the surface pressure approach of Smith, Dukowicz, and Malone (1992) and Dukowicz, Smith, and Malone (1993) (Section 6.8). In the present section, a derivation of the barotropic vorticity equation is provided.
6.4.1 Tendencies for the vertically averaged velocities
The first step is to integrate the horizontal velocity equations (4.1) and (4.2) vertically over the full depth of the rigid lid ocean
∂tU− f V = − H
The vertical integral of depth dependent quantities (mod the Coriolis force) has been lumped
These equations are implemented in MOM through a straightforward vertical integration of the forcing terms.
To facilitate physical interpretation of the surface terms forcing the vertically integrated momentum, it is useful to perform some manipulations onX0andY0. First, use the bottom kinematic boundary condition (4.24) and the surface kinematic conditionw(0) =0 in order to bring the vertical integral of the convergence of the advective flux to the form
−
whereα=u,v. Second, recall from equations (9.187) and (9.193) that the horizontal momentum friction can be written as a Laplacian acting on the horizontal velocity, plus an extra “metric”
term
F=∇h·(Am∇hu)+Fmetric. (6.20)
As such, the depth integral of momentum friction can be written Z 0
are the surface and bottom stresses (dyn cm−2) due to the winds at the ocean surface and friction and topography at the bottom. The above results render
X0= ∆(τλ)
6.4. THE BAROTROPIC VORTICITY EQUATION 63 Using an overline to denote a vertical column average, and using the expression (6.14) for the horizontal pressure gradient, yields the verticallyaveragedvelocity equations
∂tu = f v− 1
represent the vertical average of the horizontal friction metric terms, advection metric terms, divergence of the horizontal viscous fluxes, and convergence of the horizontal advective fluxes.
6.4.2 The barotropic vorticity equation
In order to eliminate the lid and atmospheric pressures, it is sufficient to form the time tendency of the barotropic vorticity
is the planetary vorticity gradient. The forcing for tendencies in ζ consists of the following terms:
1. Meridional advection of planetary vorticity: −βv.
2. Curl of the Coriolis force, which takes the form of the convergence of the barotropic velocity weighted by the Coriolis parameter: −f∇h·uh. For a flat bottom, this term vanishes with the rigid lid. Combined with the βv term, these provide the forcing
−∇h·(fuh) to the barotropic vorticity.
3. The antisymmetric term proportional to∂λ(∂φpb)−∂φ(∂λpb). This term vanishes in a barotropic model, or in a baroclinic model with a flat bottom.
4. Curl of the depth weighted surface minus bottom stresses: ˆz· ∇ ∧(∆τ/ρoH).
5. Curl of the nonlinear lateral friction terms and advection terms embodied byΓ.
To touch bases with familiar textbook dynamics, note that a steady state ocean withΓ =0 and a flat bottom will result in the familiar barotropic Sverdrup balance
βv=zˆ· ∇ ∧ ∆τ ρoH
!
. (6.33)
6.4.3 Caveat: inversions with steep topography
In order to solve for the barotropic velocityu, it is necessary to invert the elliptical operator appearing in the relation (6.30) relating the barotropic vorticity to the barotropic streamfunction.
The presence of the inverse depth in this operator implies that near regions whereHchanges rapidly, the elliptic operator also changes rapidly. As emphasized by Smith, Dukowicz, and Malone (1992) and Dukowicz, Smith, and Malone (1993), numerical elliptic inversions will potentially have problems converging in the presence of such regions. This form for the elliptic operator is also the central reason for the Killworth topographic instability (Killworth 1987). This instability can result in significant time step constraints in rigid lid models with non-smoothed topography, and these time step constraints can be more stringent than the typical CFL constraints. More problematic from the perspective of very long-term climate modeling, the Killworth instability can sometimes be slowly growing, and may become problematical only after some few hundreds of years. A slow growing instability is arguably more pernicious than rapidly growing instabilities, since a significant amount of computation time can be used before encountering the problem. Experience at GFDL indicates that the Killworth instability is a nontrivial problem with realistic models, and so it provides motivation to avoid the rigid lid streamfunction method in such models.
6.5 Boundary conditions and island integrals
6.5.1 Dirichlet boundary condition on the streamfunction
The no-normal flow boundary condition implies that next to lateral boundaries, ˆ
n·H(u,v) = nˆ ·zˆ∧ ∇hψ
= −tˆ· ∇hψ
= 0, (6.34)
where ˆnis a unit vector pointing outwards from the boundary, with the interior of the closed domain to the left, and ˆt =zˆ∧nˆ is a unit vector parallel to the path traversing the boundary.
6.5. BOUNDARY CONDITIONS AND ISLAND INTEGRALS 65 This constraint says that the streamfunction is a constant along the boundaries. In the parlance of applied mathematics, such boundary conditions are known as Dirichlet conditions. In general, the streamfunction can be a different time dependent constant along the different closed boundaries
ψ=µr(t), (6.35)
where µr(t) is a time dependent number, and r = 1,2, ...R labels the particular boundary (an island label), with Rthe total number of islands. The interpretation afforded by Stokes’
Theorem, discussed in Section 6.1, indicates thatµr(t) represents the time dependent volume transport circulating around the island with labelr.
The presence of Dirichlet boundary conditions on the streamfunction indicates that the streamfunction at one point along the boundary is identical to the streamfunction at another point, which can generally be thousands of kilometres away. Such non-local forms of infor-mation can be quite problematical in the solution of the streamfunction on parallel machines.
The discussions in Smith, Dukowicz, and Malone (1992), Dukowicz, Smith, and Malone (1993), and Dukowicz and Smith (1994) highlight this point.
6.5.2 Separating the streamfunction’s boundary value problem
The purpose of this section is to develop an algorithm for solving the streamfunction’s bound-ary value problem (BVP). To do so, reconsider the horizontal momentum equations written in the form
ut=−∇h(p/ρo)+G, (6.36)
where the vector G represents all the remaining terms given in equations (4.1) and (4.2).
Taking the vertical average of this equation, substituting the definition (6.6) for the barotropic streamfunction, and using the expression (6.14) for the hydrostatic pressure gradient, yields
1
Taking the curl of this equation eliminates the atmospheric and lid pressures
∇h∧
can be separated into two simpler BVPs. The first one is a forced elliptical problem with homogeneous boundary conditions
This system can be solved forψousing some time-stepping scheme, such as leap-frog. The sec-ond BVP is a time-independent unforced elliptical problem with constant boundary csec-onditions on the islands
∇h∧ 1
Hzˆ∧ ∇hψr
= 0 everywhere, except island boundaries (6.43)
ψr = 1 on islands. (6.44)
This BVP can be solved for ψr using some type of an elliptical solver, such as congugate gradient. The full streamfunction is built from the sum
ψ(λ, φ,t)=ψo(λ, φ,t)+ XR
r=1
µr(t)ψr(λ, φ), (6.45)
which can be shown to satisfy the original boundary value problem.
6.5.3 Island integrals for the volume transport
As a final step in the streamfunction solution, it is necessary to formulate a prognostic equation for the volume transportsµr(t). To do so, reconsider the elliptical problem for the streamfunc-tion
Now integrate both sides over the area bounding a particular island with labelr, and employ Stokes’ Theorem. The direction normal to the island’s surface is ˆz, and the area element on the island isdΩr. The left hand side becomes Equating yields the prognostic equations for the volume transports around an island
µ˙r = 1 ρo
H dltˆ·(∇hpb−G) H dlnˆ ·H1∇hψr
. (6.49)
6.6. THE BAROCLINIC MODE 67
6.6 The baroclinic mode
The time tendency for the zonal baroclinic velocity is given by
∂tbu = ∂t(u−u)
= 1− 1 H
Z 0
−H
dz
!
Gu− 1 aρ◦ cosφ
!
(pb+pl+pa)λ
!
(6.50) where the vectorGwas introduced in equation (6.36). Since the atmospheric and lid pressures are depth independent, these pressures have zero deviation from their depth average. Hence, the tendency for the baroclinic velocity is independent of the atmospheric and lid pressures.
In determining the baroclinic velocity, it is therefore sufficient to consider
∂tbu = ∂t(u′−u′) (6.51)
whereu′andu′satisfy the full zonal velocity equation and vertically averaged zonal velocity equation, respectively, but without any atmospheric or lid pressure contributing to the forcing.
Without the atmospheric and lid pressure contributions, all forcing terms in the two “primed”
equations are known. Consequently, the primed velocities can be time stepped in a straight-forward manner without using any tricks such as those needed to eliminate the atmospheric and lid pressures from the streamfunction equation. Time stepping the primed velocities then yields the updated baroclinic velocity through equation (6.51). After obtaining the updated baroclinic velocity, the full horizontal velocity field
(u,v) =(u,v)+(bu,bv) (6.52)
is known at the new time step. Formulation of the rigid lid streamfunction method is now complete.
6.7 Summary of the rigid lid streamfunction method
In summary, the rigid lid streamfunction method allows for the computation of the total velocity field using the following algorithm:
• Baroclinic mode
1. Compute the forcing terms for the “primed” momentum equation, which include all terms except the lid pressure and atmospheric pressure (Section 6.6).
2. Solve the primed momentum equation at each model level to get the primed velocity at the next time step.
3. Subtract the vertically averaged primed momentum from the primed momentum in order to get the baroclinic velocity field at the next time step.
• Barotropic mode
1. Compute the forcing terms in the barotropic vorticity equation (Section 6.4).
2. Solve the elliptical boundary value problem for the time tendency of the barotropic streamfunction.
3. Time step the streamfunction forward and compute the barotropic velocity field at the new time step.
• Add the baroclinic velocity to the barotropic velocity to then have the full velocity field.
6.8 Rigid lid surface pressure method
The rigid lid surface pressure method in MOM was developed by Smith, Dukowicz, and Malone (1992) and Dukowicz, Smith, and Malone (1993). A useful discussion of the analyt-ical issues related to the surface pressure formulation can be found in the paper by Pinardi, Rosati, and Pacanowski (1995). The basic advantages over the streamfunction approach are the following:
1. The elliptical operator has anHfactor, rather than theH−1 found in the streamfunction approach. This property makes the elliptical operator more readily inverted than with the streamfunction approach (see Section 6.4.3). In particular, steep topography is much less of an issue with the surface pressure approach.
2. The boundary conditions for the elliptical problem are Neumann rather than Dirichlet.
As discussed by Smith, Dukowicz, and Malone (1992) and Dukowicz, Smith, and Malone (1993), Neumann boundary conditions are much easier to handle on parallel machines than Dirichlet conditions.
As the formulation of the rigid lid surface pressure method follows quite closely the logic of the streamfunction approach, the details are not reproduced here. Rather, reference should be made to the Smith, Dukowicz, and Malone (1992), Dukowicz, Smith, and Malone (1993) and Pinardi, Rosati, and Pacanowski (1995) papers.