.left margin 13 .right margin 90 .tab stops 35,75 .nojustify .spacing 2 .figure 5 .center;WIND-DRIVEN SPINUP IN EDDY-RESOLVING OCEAN MODELS .center;FORMULATED IN ISOPYCNIC AND ISOBARIC COORDINATES .s 4 .center;by .s 4 .center;Rainer Bleck and Douglas Boudra .s 2 .center;Rosenstiel School of Marine and Atmospheric Science .center;University of Miami .center;Miami, Florida 33149 .s 4 .center;December 1985 .page .center;A B S T R A C T .p Wind-driven spinup of the four-layer, quasi-isopycnic, eddy-resolving primitive equation model of Bleck and Boudra (1981) is compared to that obtained with a (numerically dissimilar) 'pure' isopycnic coordinate model and an isobaric (i.e., quasi-Cartesian) coordinate model. In particular, the onset of hydrodynamic instabilities in the flow forced by a double-gyre wind stress pattern is studied. The spin-up processes associated with the isopycnic and quasi-isopycnic model are found to be similar, whereas the flow pattern produced by the quasi-Cartesian model deviates in the direction of Holland's (1978) and Holland and Lin's (1975a,b) two-layer solutions. .page .p $1. Introduction$ .p One of the principal guidelines in approximating partial differential equations by algebraic finite difference expressions is numerical $consistency$. Testing for consistency is equivalent to testing whether the truncation error goes to zero with decreasing mesh size, and the rate at which this happens is often used as the primary indicator of the $accuracy$ of a finite difference scheme. However, not in all situations is the validity of a consistent solution determined by the asymptotic behavior of the truncation error. Exact conservation of first or second moments of the numerical solution, for example, may be equally important in long-term integrations of the hydrodynamic equations (Arakawa, 1966). Another factor limiting the validity of a numerical solution is the spurious transport of scalar properties across isopycnic surfaces in three-dimensional models of stratified, incompressible fluids. In models of finite resolution this diapycnic transport apparently can only be suppressed by treating the isopycnals themselves as coordinate surfaces. Attempts to overcome this problem in Cartesian coordinate models by using higher order advection-diffusion schemes (Redi,1982) are bound to fail in the long run because they cannot remove from the solution the dispersive effects of nonzero phase error in the horizontal advection terms. Isopycnic coordinate models are not exempt from this phase error (which is known to be the main cause for the gradual deterioration of finite difference solutions of hyperbolic systems in general) but do not let it result in material transport across isopycnic surfaces. .p Another desirable property of models with coordinate surfaces tied to the density stratification is that frontogenetic processes are resolved by the distortion of the coordinate system itself. This feature has been found useful in the modeling of atmospheric cyclone waves (Shapiro, 1975), but fundamental problems concerning intersections of coordinate surfaces with the ground have prevented a widespread adoption of atmospheric isentropic coordinate models. .p Until recently, the outcropping of layer interfaces at the sea surface has also been a major impediment to multi-layer isopycnic coordinate ocean modeling. One of the techniques proposed for overcoming this problem is the 'hybrid' coordinate approach of Bleck and Boudra (1981), hereafter referred to as BB. In the BB model, coordinate surfaces coincide with isopycnals when- and wherever this is possible without violating a minimum layer thickness requirement; otherwise, the coordinate surfaces assume a quasi-Cartesian character. Some modelers may react with skepticism to the notion of allowing an individual grid point to float up and down with an isopycnal at one time while holding it in a fixed position at another time. However, the prognostic equations permitting this kind of flexibility are shown in BB to be hardly more complex than the equations in conventional primitive-equation models. .p In meteorology, the relative merits of different model representations are usually tested against observed data, but this avenue is not readily available to ocean modelers due to the relative lack of observations. It is therefore natural to begin by documenting, and perhaps understanding, the differences generated when several models are run under the same initial and boundary conditions. .p BB tested a four-layer version of their fully nonlinear eddy-resolving model in the situation of a double-gyre, wind-driven spinup patterned after the two-layer experiments of Holland and Lin (1975a,b) and Holland (1978). We shall refer to those papers jointly with the acronym HLH. Since the primary purpose of BB's work was to demonstrate the viability of a novel vertical coordinate system under extreme forcing conditions, they chose a 50% higher wind forcing than HLH, thereby sacrificing the model intercomparison aspect of their work. Substantial differences in the resulting motion fields were indeed observed, leading naturally to the question of whether one of the two models was deficient in representing the underlying dynamic processes. One must be careful, here, not to invoke similarity with observed oceanic conditions, since neither the geometry (particularly the absence of bottom topography) nor the forcing was chosen to be particularly realistic. To address the intercomparison aspect, an 18-year experiment with half the forcing amplitude has since then been performed, which exhibits circulation patterns very similar to the double-gyre experiments of HLH. Results from this experiment will be reported elsewhere. .p In the present paper we focus attention on one particular aspect of model behavior. We wish to examine the initial wind-driven spinup of models with three different coordinate representations: the quasi-isopycnic coordinate model of BB, a newly developed 'pure' isopycnic coordinate model, and an isobaric coordinate model, similar to the Bryan (1969) z-coordinate representation, but with eddy-resolving horizontal grid spacing. In particular, we wish to look at how the three models represent the initial onset of hydrodynamic instabilities, i.e., the breakdown of the free jet in the middle of a double-gyre wind forced circulation pattern. It will become apparent that the differences between the two material coordinate-oriented models in the character of the spinup process are minor compared to the differences between quasi-Lagrangian and Cartesian coordinate representation. .p In the layout of the paper, Sections 2, 3 and 4 deal with the technical aspects of the models used. There follows a discussion of energy diagnostics in the different coordinate systems (Section 5), while Section 6 describes the experimental setting and the characteristics of the three sets of model runs. .p The model differences in representing the breakdown of the free jet are explained largely in terms of the associated energy conversions. The question of what causes the meridional displacement of the separation point in the two material coordinate models remains unanswered, as does the perhaps related tendency for strong anticyclonic eddy formation near the western boundary. These topics are beyond the scope of the current paper, but will be addressed in the future, when more extensive long-term experiments with eddy-resolving primitive equation models are economically feasible. .s .test page 9 .p $2. The non-hybrid isopycnic coordinate model$ .p This model takes advantage of the recently developed 'Flux Corrected Transport' (FCT) method (Boris and Book, 1973; Zalesak, 1979) for calculating mass flux divergences in finite difference form. The impact of this method on isopycnic modeling is perhaps best understood by viewing the three-dimensional isopycnic model as a stack of shallow fluid layers, each one characterized by a constant value of density. The motion within each layer is governed by equations resembling the shallow water equations. The layers interact primarily through pressure forces. Frictional stresses along the interfaces may be present but should not be regarded as the primary means of communication among the layers. .p The main numerical difficulty with this modeling approach is that individual fluid layers must be allowed to 'dry up' anytime and anywhere in the grid domain. In other words, the numerical algorithm for handling the horizontal mass exchange among adjacent grid boxes must be capable of allowing areas of zero layer thickness to co-exist with areas of non-zero layer thickness without creating negative thickness values in the transition zone. The FCT algorithm prevents this otherwise common overshooting by breaking up the mass flux divergence calculation into two steps. In the first step, mass fluxes across grid box boundaries are calculated using a diffusive, low-order scheme incapable of causing overshooting. In the second step, these diffusive fluxes are compared pointwise to fluxes obtained from a more accurate, less diffuse, but ripple-prone method (accuracy measured here in terms of the order of the truncation error). The difference between the low- and high-order flux calculations is aptly named 'antidiffusive' flux. The essential step in the FCT algorithm is to reduce the antidiffusive flux values pointwise until they can be added to the previously computed diffusive solution without creating new maxima or minima, i.e., new ripples. It is easy to see that this procedure regenerates steep gradients lost in the diffusive treatment without permitting overshooting. (As an aside, the procedure is also known to exaggerate gradients in smoothly varying fields.) An important aspect of the FCT algorithm is that it can be implemented efficiently on modern vector-processing computers. .p While the above algorithm for the first time permits unrestricted multi-layer isopycnic modeling in the true sense of the word, it should not be regarded as superior in every respect to the hybrid, quasi-isopycnic approach. Since the layer thickness in the 'pure' isopycnic model is not bounded away from zero, quadratic-conservative finite difference expressions derived from the flux form of the prognostic equations cannot be used to evaluate the nonlinear terms in that model in regions where the layer thickness may decrease to zero in the course of the integration. (For the same reason, use of the flux form itself is inappropriate.) As discussed in BB, the hybrid model conserves layer $potential$ vorticity and $potential$ enstrophy at least in those grid regions where the coordinate surfaces coincide with isopycnals (which is generally the case in the ocean interior). On the other hand, the finite difference operators in the purely isopycnic model, where division by layer thickness is not allowed, are quadratic-conservative only in the limited sense that they prevent false generation of enstrophy by the $nondivergent$ part of the flow field. .p Specifically, the finite difference forms of the momentum equations in the isopycnic model (omitting friction for the time being) are .spacing 1 .i 15 &#&#|x##&#&#|y .i 5 ^{&\q&u}|{ \qt} + \d|x[^{&1}|{ 2}(u^{\2} + v^{\2} )] - v\%^{^{xy}} \h\%^{^y} = - \d|xM (1a) .s 2 .i 15 &#&#|x##&#&#|y .i 5 ^{&\q&v}|{ \qt} + \d|y[^{&1}|{ 2}(u^{\2} + v^{\2} )] + u\%^{^{xy}} \h\%^{^x} = - \d|yM (1b) .spacing 2 .s where \h = \d|xv-\d|yu+f is the absolute vorticity and M \Z o\%\/^{^s} + \ap\%^{^s}. These equations follow directly from (BB-2a,b) by canceling all references to layer thickness in the term involving Q, dropping the vertical advection term and combining the horizontal pressure force terms into one involving the Montgomery potential M. Note that M, contrary to the geopotential o\/, is used in the model as a $layer$ variable. Starting from the hydrostatic equation as used in the BB model, \d|so\/=-\a\d|sp, it can be easily shown that M satisfies the hydrostatic equation in the form .p \d|sM = p \d|s\a (2) .p Thus, changes in M coincide with changes in \a. As these only happen at layer interfaces, (2) implies a depth-independent pressure force in each coordinate layer which in turn is consistent with the assumption of a depth-independent velocity field. .p There is no thermodynamic equation to solve as long as diabatic processes, compressibility effects, and temperature and salinity conservation equations are excluded. The continuity equation is identical to (BB-3), except for the omission of the vertical flux term: .p ^{&\q}|{ \qt}\d|sp + \d|xU + \d|yV = 0 (3) .s The reader should keep in mind, however, that due to the use of the FCT algorithm the mass flux components (U,V) in (3) are no longer given simply by (u\d|sp\%^{^x},v\d|sp\%^{^y}) as in the quasi-isopycnic model, but are a mixture of low-order and high-order fluxes. We presently use the upstream or 'donor cell' scheme to calculate the low-order mass fluxes, and the fourth-order expressions .p U = u[^{&9}|{ 8}\d|sp\%^{^x} - ^{&1}|{ 8}\d|sp\%^{^{3x}}] .break .indent 5 V = v[^{&9}|{ 8}\d|sp\%^{^y} - ^{&1}|{ 8}\d|sp\%^{^{3y}}] .s for the high-order fluxes. .p An important question is how to assign grid point values of u and v in regions where the layer depth is zero, because some of these data points are obviously needed to evaluate inertial and viscosity terms at grid points adjacent to 'massless' regions. The strategy followed at present is to solve the momentum equations at $all$ grid points regardless of the layer depth. Note that the hydrostatic equation (2), if integrated across massless isopycnic layers at the ocean surface (p=0), yields an M field that does not vary from one massless layer to the next. Consequently, the balance of forces in the massless layers should always adhere closely to the balance in the topmost layer exhibiting a nonzero layer depth. The same degree of vertical coherence should then be found in the u and v field. .p Massless regions at the $bottom$ of the model domain do not occur in the present experiment, but must be taken into account if the model is run with variable bottom topography. Since p>0 at the bottom, M will change now from one massless layer to the next according to (2). Therefore, an important question to ask is whether the hydrostatically determined M field in layers $above$ the bottom is affected in any way by the existence of an arbitrary number of massless isopycnic layers sandwiched between the ocean bottom and the first layer exhibiting a nonzero thickness. Fortunately, the M field turns out to be invariant. However, since p>0, the balance of forces and the ensuing (u,v) field $within$ the massless layers at the bottom will differ now from what is found in the first 'real' layer above. A possible remedy to this problem and the related one of how to calculate the horizontal pressure force in an isopycnic layer intersecting the ocean bottom is discussed in the atmospheric context in Bleck (1984). .s .test page 9 .p $3. The isobaric model$ .p This model is a special case of the hybrid BB model. To obtain it from the latter requires no code changes other than keeping the layer thickness \d|sp constant in time. As a consequence, the momentum equations (BB-2a,b) reduce to .spacing 1 .i 37 &#&#&#&#&#|s .i 15 &#&#|x##&#&#|y##############^{&#}x .i5 ^{&\q&u}|{ \qt} + \d|x[^{&1}|{ 2}(u^{\2} + v^{\2} )] - v\%^{^{xy}} \h\%^{^y} + s\A \d|su = - \d|xo\/\%^{^s} (4a) .s 2 .i 37 &#&#&#&#&#|s .i 15 &#&#|x##&#&#|y##############^{&#}y .i 5 ^{&\q&v}|{ \qt} + \d|x[^{&1}|{ 2}(u^{\2} + v^{\2} )] + u\%^{^{xy}} \h\%^{^x} + s\A \d|sv = - \d|xo\/\%^{^s} (4b) .spacing 2 .s with s\Zp. The hydrostatic equation \d|so\/=-\a\d|sp remains unaffected by the change to isobaric coordinates while the continuity equation (BB-3) assumes the simple form .p \d|xu + \d|yv + \d|ss\A = 0 (5) .s The thermodynamic equation (BB-5) reduces to .spacing 1 .i 10 &#&#&#&#|x##&#&#&#&#|y##^{&#&#&#&#}s .i5 ^{&\q&\a}|{ \qt} + u\d|x\a + v\d|y\a + s\A\d|s\a = 0 (6) .spacing 2 .s Note that we are using the specific-volume conserving form (BB-5), as opposed to the energetically consistent form (BB-4). It was stated in BB that retention of the pressure weighting of s\A in the vertical advection term had little effect on the outcome of the hybrid model calculation. However, since vertical motion is more prevalent in a Cartesian than in a quasi-material coordinate model, the isobaric model solutions did turn out to be sensitive to the choice of the vertical advection term. We found that the gradual shift in model mean density caused by the pressure weighting of s\A was too serious to permit the use of the energetically consistent formulation. Once this decision was made, it seemed reasonable to convert the thermodynamic equation in the hybrid model to the specific-volume conserving form (BB-5) as well. .p Note that the horizontal advection terms in (4a,b) are identical to those in (1a,b). Thus, the isobaric and the isopycnic model have the same quadratic conservation properties which are inferior to those of the hybrid model. .s .test page 9 .p $4. Numerical Details Common to All Model Versions$ .p Lateral mixing of momentum in all three model versions, and of specific volume in the isobaric model (Eq.(6)), is handled by the biharmonic lateral diffusion terms (BB-6) in conjunction with the deformation-dependent eddy viscosity (BB-7). Note that the grid resolution in HLH and BB is fine enough to resolve a substantial scale-range of baroclinic eddies. An 'eddy' viscosity is still needed in these calculations to simulate the cascading of enstrophy past the grid scale. .p The weighting of the viscous fluxes in (BB-6) by the layer thickness \qp/\qs, which is done for the purpose of conserving the first moment of the quantity being diffused, evidently cannot be carried out rigorously in the isopycnic model version where \d|sp is not bounded away from zero. In order to avoid division by zero we do not permit \d|sp in (BB-6) to drop below one-tenth of the initially specified layer thickness value. Fluid elements located near the point where their coordinate layer intersects the ocean surface are thus subjected to a slight drag from neighboring, 'massless' grid points. We estimate this drag to be small compared to the balance of forces at those points. .p The proportionality factor linking the eddy viscosity to the total deformation was chosen to be 0.04 as in the low-viscosity experiments reported in BB. The same formulation is used for the lateral mixing of \a in (6) and (BB-5). .p A rigid-lid approximation is used in all three model versions to remove barotropic gravity waves from the solutions, allowing a model time step of slightly less than one-half hour. As described in Appendix D of BB, those waves are filtered by removing the divergent component from the vertically integrated mass flux field after each time step. Caution must be exercised in the isopycnic model to prevent the FCT algorithm from perturbing the nondivergent character of the barotropic flow. As outlined in Section 2, horizontal mass transfer in the isopycnic model is accomplished by a low-order, diffusive flux component supplemented by an antidiffusive component whose magnitude often is such that it barely avoids the generation of negative layer thickness values. If we now define the barotropic mass flux as the vertical integral over the sum of the diffusive and antidiffusive flux components, negative layer thickness values may result if the divergence is removed from the resulting flow field. .p We presently deal with this problem in the isopycnic model by (a) imposing the nondivergence condition on the vertically integrated $high-order$ fluxes, rather than the fluxes modified by the FCT algorithm; (b) keeping track of the amount by which the FCT algorithm reduces the antidiffusive fluxes in a given grid column; and (c) distributing a mass flux equal and opposite to the amount recorded in step (b) evenly among all layers in the grid column. In order to avoid the generation of negative layer thickness values in the process, the corrective flux is prorated among the coordinate layers according to the distribution of layer thickness values in the $upstream$ grid column (upstream relative to the direction of the corrective flux). .p Domain size, vertical density structure, and wind stress pattern in the current experiment are identical to those used in BB. The model ocean is spun up from rest by an east-west surface wind stress pattern of the form .p \t = \t|{\0}cos[2\p(^{&y}|{ Y} - ^{&1}|{ 2})] .s where Y = 2400km is the north-south basin size and \t|{\0} = 1.5 X 10^-^4 m^2s^-^2. The model consists of four layers whose initial thickness values, counting from the top on down, are 233, 333, 433, and 4000 dbar respectively. Initial values of \a in the four layers are 1.000, 0.999, 0.998, and 0.997 X 10^-^3 m^3 kg^-^1 respectively. Using Lighthill's (1969) algorithm we have determined that this stratification implies equivalent depth values of 1.20, 0.217, and 0.102m for the three baroclinic modes. Using the basin-wide average of f, 0.83x10^{-\4}s^{-\1}, these numbers translate into 38, 16, and 11km for the corresponding radii of deformation. The wave length 2\pR|d marking the threshold between dispersive and nondispersive baroclinic Rossby waves is therefore 238, 101, and 69km respectively for the three modes. We may conclude that the nondispersive part of the Rossby wave spectrum is resolved, at least marginally, for all three baroclinic modes, whereas only those dispersive waves associated with the $first$ baroclinic mode can be regarded as being partly resolved in our 25km mesh. .p Other parameters taken directly from BB are the east-west basin size, X = 1200km, the Coriolis parameter at y=Y/2, f|{\0} = 0.83x10^{-\4}s^{-\1}, and \b = 2x10^{-\1\1}s^{-\1}m^{-\1}. As in BB, the wind stress is assumed to decrease linearly to zero over a depth of 100 dbar. This is tantamount to a uniform forcing of the top layer until its thickness becomes less than 100 dbar. If this happens, the linear stress law yields a prescription of how to distribute the wind forcing among coordinate layers found within 100 dbar of the surface. There is no sidewall friction, but the lowest of the four coordinate layers is subjected to a linear bottom drag with a drag coefficient of 5x10^{-\8}s^{-\1}. Interfacial stresses, other than those related to the wind stress in the uppermost 100 dbar, are zero. .p In BB, the first five years of model spinup were performed on a 50km mesh, after which time the fields were interpolated to a 25km mesh. Integration then continued on the finer mesh for several more years. This switch in grid size has been eliminated: the present experiments are carried out on a 25km mesh from the very beginning. The spinup phase is therefore not directly comparable between the present experiments and BB. .p We conclude this section by listing in Table 1 the most relevant properties of the three model versions. .s .spacing 1 .test page 20 .literal T A B L E 1 ----------------------------------------------------------------------------- I Isobaric I Hybrid I Isopycnic I I Model I Model I Model I ----------------------------------------------------------------------------- mass fluxes I second-order I second-order I Flux Corrected I in continuity I space-centered I space-centered I Transport (up- I equation I I I (stream/4th order)I ----------------------------------------------------------------------------- quadratic I vorticity and I pot.vorticity and I vorticity and I conservation I enstrophy (by I pot.enstrophy (in I enstrophy (by I properties I nondiv.velocity) I ocean interior) I nondiv.velocity) I ----------------------------------------------------------------------------- layer I fixed (233,333, I variable, but no I variable, I thickness I 433,4000 dbar) I less than 10 dbar I unrestricted I ----------------------------------------------------------------------------- .end literal .spacing 2 .s .test page 9 .p $5. Model Energetics$ .p The three model versions outlined in the previous sections were integrated over five years (in one case seven years) of model time. While the solutions are similar in many respects, there are differences which may be explained in terms of the truncation properties of Cartesian vs. material vertical coordinates. Energy diagnostics play a major role in this analysis. .p In order to study energy transformations in the numerical solutions, all model output fields must be decomposed into time mean and eddy fields. In the case of variable layer thickness (including the special case of completely deflated layers), this decomposition is most appropriately carried out by subjecting all variables carried on layer interfaces (i.e., p,o\/,s\A) to an unweighted time average, denoted in the following by an overbar, while subjecting all 'layer' variables (u,v,\a) to a mass-weighted time average defined as .spacing 1 .i 9 &#&#&#&# .i 5 q\~ = q\d|sp / \d|sp\%, .spacing 2 .s The respective eddy components are then defined as q*=q-q\% and q'=q-q\~. As an example, the mean and eddy components of kinetic energy are represented in this framework by (\d|sp\%)v\Y\~^{\2}/2 and (\d|sp)v\Y'^{\2}/2 respectively. Note that there is no difference between q\% and q\~ in isobaric coordinates, and thus no difference between q* and q'. .p In an attempt to compare the energetics of isobaric and isopycnic coordinate models we encounter a problem having to do with the definition of eddy potential energy (P|E). The rigorous definition of this quantity in terms of the variance of pressure on isopycnals has no numerically tractable, yet exact, counterpart in the isobaric reference frame. In particular, it does not agree with the customary definition of P|E in terms of the density variance on p surfaces originally introduced by Lorenz (1955). This by itself would already complicate the interpretation of energetic differences between the two model types. However, there is the additional problem that the conversion terms between the eddy and mean components of potential and kinetic energy, if viewed in the isopycnic framework, do not follow the classical pattern (Lorenz, 1967) which disallows direct exchange between P|M and K|E and between P|E and K|M. Instead, a conversion term appears in the isopycnic framework suggestive of an exchange process between P|E and K|M while there no longer is an obvious term linking P|M and P|E. .p A detailed derivation of the energetics in isopycnic coordinates is given by Bleck (1985) who also suggests that an energy cycle formally similar to Lorenz' cycle can be constructed by 're-routing' certain exchange processes through intermediate (virtual) energy states. However, these manipulations do not do justice to the physical makeup of the terms in question. Given our present level of understanding, it seems wisest to avoid the breakdown of P into P|M and P|E in the present study altogether and to consider energy exchange processes in the 'triangle' formed by P, K|M, and K|E only. Since isopycnic and isobaric energetics agree in the sense that neither permits a direct exchange between P|M and K|E, the exchange between P and K|E, which forms one side of the triangle, may be interpreted as being equivalent in all model versions to the exchange between P|E and K|E associated with the release of baroclinic instability. Likewise, the barotropic instability process associated with K|M-to-K|E conversion is accounted for in another side of the triangle. .p Following the generalized coordinate notation used by Bleck (1985), the P-to-K|E conversion term in both isobaric (s=p) and isopycnic (s=\a) coordinates can be written in the form .spacing 1 .s .i 8 &#&#&#&#&#&#&#&#&#&#&####&#&#&#&#&#&#&#&#&#&# .i 5 C(P,K|E): o\/V|s\:(^{&\q&p}|{ \qs} v\Y') - v\Y'\:^{&\q&p}|{ \qs} \aV|sp (7) .s 2 whereas the P-to-K|M conversion term appears as .s .i 8 &#&#&#&#&#&#&##########&#&#&#&#&#&#&# .i 5 C(P,K|M): o\/V|s\:(^{&\q&p}|{ \qs} v\Y\~) - v\Y\~\: ^{&\q&p}|{ \qs} \aV|sp (8) .s 2 The form of the K|E-to-K|M conversion term is .s .i 8 &#&#&#&#&#&#&#&#&#&#&#&#&#&#&#&#&#&#&# .i 5 C(K|E,K|M): v\Y'^{&\q&p}|{ \qs}\:(v\Y'\:V|s + s\A'^{&\q&#}|{ \qs}) v\Y\~ (9) .spacing 2 .s Needless to say, the second term in (7) and (8) vanishes in isobaric coordinates (s=p). The finite difference analogs of the above expressions are straightforward. .s .test page 9 .p $6. Results$ .p We begin our discussion of the three model solutions by showing in Fig.1 the basin-wide average of mean flow kinetic energy (K|M) plotted against model time. (Throughout this section, 'mean' fields represent 360-day running averages.) The graph suggests that, as far as the K|M field is concerned, statistical equilibrium is reached in all three model versions after three to four years of spinup time. No true steady state is achieved. (We wish to mention in passing that the previously referenced 18-year isopycnic model run conducted with half the wind forcing of the present experiment, 0.75 X 10^-^4m^2s^-^2, shows unabated multi-year vascillations between asymmetric, eddy-dominated flow regimes and relatively symmetric double-gyre circulations with a weakly meandering free jet. The energetic signature of these vascillations, a discussion of which is beyond the scope of this paper, can be easily distinguished from the signature of the initial spinup depicted in Fig.1.) Before reaching statistical equilibrium, the isobaric model shows signs of excessive K|M buildup. This overshooting is related, as we will demonstrate shortly, to a somewhat delayed appearance of hydrodynamic instability processes in the isobaric model version. .p The K|M curve for the isopycnic model version has the opposite appearance: it reaches a plateau after 1.5 years, but goes up another step near the end of the third year. Thereafter, the K|M level in the isopycnic model solution is noticeably higher than in the other two solutions; the excess amount turns out to be evenly distributed among the four coordinate layers. Since energy input through wind stress is the same in all experiments, a higher steady-state kinetic energy level can only be caused by a less effective energy dissipation mechanism. Closer inspection of the isopycnic and hybrid model solutions indeed reveals that the isopycnic solution, despite its more energetic character, is slightly less affected by lateral dissipation than the hybrid solution. (The bottom drag, on the other hand, is slightly higher in the isopycnic case due to the higher kinetic energy level in the bottom layer.) Since lateral dissipation (particularly the bi-harmonic formulation used here) is highly scale-selective, we may assume that the isopycnic model produces somewhat larger-scale eddies than the other two model versions. .p As shown by the values plotted at 0.5 years in Fig.2, the dominant energy conversion process in all three model versions during the first year of spinup time is the conversion from K|M to P representing a gyre-scale adjustment of the mass field to the gradual strengthening of the wind-forced circulation. Thus, it should be regarded as a K|M-to-P|M conversion. Soon after the first year, this process becomes secondary to the K|M-to-K|E conversion process (Fig.2b) signaling the onset of barotropic instability processes in the free jet. The most interesting aspect of this changeover is the timing. Whereas the jet in the isopycnic and the hybrid (quasi-isopycnic) model becomes barotropically unstable sometime in the middle of the second year, the K|M-to-K|E conversion in the isobaric model version does not jump up until the beginning of year 3. .p In view of the rather coarse time resolution in Fig.2 (which is the result of an inappropriate, but irreversible archiving decision) the significance of this timing difference should be checked against the instantaneous surface velocity fields which are available at 10 day intervals. Since the thickness of the top coordinate layer differs markedly among the three model versions, we will show velocity vectors averaged over a fixed depth. The chosen averaging interval is 300 m. .p During the first two years the surface flow pattern in the isobaric simulation is dominated by a virtually straight jet extending all across the basin (1200 km). Weak eddies begin to appear in the return flow on either side of the jet after about six months of integration. At the beginning of year 3 some of these eddies are intense enough to interact with the jet during their westward migration, causing strong meander growth in the process. The first noticeable event of this kind, which begins around day 750, is documented in Fig.3. (In this and the following figures only a part of the grid domain consisting of one-third of the cyclonic and five-sixths of the anticyclonic wind stress domain is shown.) The jet is seen to extend meanders simultaneously toward a cyclonic and an anticyclonic eddy. The end result of the episode is a total breakdown of the jet in the eastern half of the domain, a regime change which turns out to be irreversible. .p The jet in the isopycnic model simulation, which initially also extends to the eastern boundary, breaks down in a similar fashion, but it does so shortly after day 610, almost five months before the isobaric jet. The difference in timing agrees rather well with what we were able to deduce from Fig.1. .p One of the earliest indications of differences in model behavior is the inability of the free jet in the hybrid model to maintain a straight course and to reach the eastern boundary during the 'quiet' first 15 months of integration. We are presently unable to offer an explanation for this difference. The hybrid model is the only one of the three that conserves potential vorticity and potential enstrophy (at least in the case of non-surfacing isopycnals), and it is tempting to declare this to be the reason for the jet's behavior. However, recent three-layer model intercomparison studies focussing directly on the effect of replacing (-v\%^{^{xy}} \h\%^{^y},u\%^{^{xy}} \h\%^{^x}) by (-v\%^{^{xy}} \h\%^{^y},u\%^{^{xy}} \h\%^{^x}) in the momentum equations give no indication that the potential vorticity/potential enstrophy constraint reduces the stability of a free jet. In fact, this constraint is generally thought of as one that inhibits, rather than promotes, transfer of energy toward smaller scales (Sadourny, 1975). .p As a consequence of the above difference in the jet's behavior, the first meandering episode comparable in magnitude to the event shown in Fig.3 takes place in the hybrid model simulation as early as day 450. This is five months before the breakdown of the isopycnic and ten months before the breakdown of the isobaric model jet. .p The retarded appearance of instabilities in the isobaric model jet can be easily explained in terms of the truncation properties of the coordinate systems involved. Since the vertical resolution in the isobaric model is fixed, any jet-like structure that is trying to form near the surface is being diluted over a considerable vertical range (233m in the present case). In the isopycnic coordinate models, on the other hand, grid points automatically migrate vertically into the region where densely packed, sloping isopycnals are needed to provide thermal support for a developing jet stream. Another process delaying the formation of sharp horizontal density gradients in the isobaric coordinate model is the horizontal mixing of \a in the thermodynamic equation (6). It is thus hardly surprising that barotropically unstable jets form more quickly in a coarse-resolution isopycnic than a coarse-resolution Cartesian coordinate model. Nevertheless, the magnitude of this effect was not anticipated. .p The K|M-to-P conversion process (Fig.2a) does not vanish entirely with the onset of barotropic instability. In fact, in the isopycnic model version it recovers significantly after the third-year minimum associated with the peak in K|M-to-K|E conversion, while continuing on a fairly steady level in the other two model versions. The role of the K|M-to-P conversion process during the later years of the experiment becomes clear after comparing it to the P-to-K|E conversion process shown in Fig.2c: it functions as the energy source for baroclinic instability which, as pointed out in the previous section, is represented in our analysis by the P-to-K|E term. .p Consistent with the high mean kinetic energy level in the isopycnic model and the associated strong baroclinity of the mean density field (Fig.1), baroclinic instability is seen in Fig.2c to play a larger role in this model than in the other two. In fact, during the fifth model year K|E is generated by baroclinic conversion at such a high rate that a net conversion of K|E to K|M results. This anomaly, incidentally, was the reason for continuing the integration of the isopycnic model for another two years. While the isopycnic model is seen to 'quiet down' in the sense that its energy conversion processes become comparable in magnitude to those in the other models, the gradual buildup of mean kinetic energy during years 6 and 7 (Fig.1) makes it seem likely that instability processes in the isopycnic model will eventually intensify again. .p Figs. 4 and 5 are intended to illustrate the spatial distribution of mean and eddy kinetic energy in the three model versions at various times during the spinup process. As before, these figures only show 33% of the northern (cyclonic) and 83% of the southern (anticyclonic) wind stress domain. Thus, the latitude of maximum eastward wind stress, marked by a dotted line, is located approximately one-third down from the top of each frame. The fields shown are those for the second (Fig.4) and third year (Fig.5) of model time, corresponding to data points plotted respectively at 1.5 and 2.5 years in Fig.1 and 2. .p The circulation in the isobaric model during model year 2 (Fig.4, top) is characterized by two western boundary currents merging into an extremely straight and steady free jet which is positioned approximately 100km south of the wind stress maximum. Judging from the top left panel of Fig.4, there is hardly any eddy motion in the isobaric model at this time. .p The free jet in the hybrid and isopycnic model (Fig.4, center and bottom), on the other hand, appears to spawn some eddy activity during model year 2, as indicated by the presence of eddy kinetic energy contours in the second and third panel on the left. The jet in the hybrid model seems to be slightly less steady in the eastern half of the basin than the jet in the isopycnic model, and the eddy kinetic energy level seems slightly higher, particularly in the recirculation region of the southern gyre. However, these differences are minute compared to the amount by which both solutions deviate from the isobaric one. .p The situation during model year 3 is shown in Fig.5. At this time, the free jet in the isobaric model (top) no longer follows a steady path across the basin, but starts to meander a short distance out from the western coastline. It is tempting to interpret the patchiness in the eddy kinetic energy field (top left) as an indication of a preferred wavelength for instabilities along the free jet in the isobaric model. This wavelength is approximately 250km. Similar, but weaker, patterns can be detected in the eddy kinetic energy fields obtained from the hybrid and isopycnic model during year 2 (second and third panel on the left of Fig.4). The preferred wavelength in the isopycnic model (bottom panel) appears to be near 400km. .p The mean kinetic energy fields for the hybrid and isopycnic model during year 3 (second and third panel on the right of Fig.5) continue to show a steady western boundary current extending south from the cyclonic gyre domain, but its anticyclonic counterpart and the free eastward jet have become too unsteady to leave a trace in those fields. The only persistent feature in the southern domain is the strong anticyclonic eddy near the western boundary which also dominated the solutions obtained by BB. .p We conclude the discussion by showing in Fig.6 snapshots of the near-surface flow in the three models at the end of the five-year integration. The flow patterns at that time are seen to be in near-perfect agreement with the patterns suggested by the third-year energy plots in Fig.5, strengthening our earlier contention that the models have reached a quasi-steady state well before the end of the fifth year. The isobaric model (Fig.6a) apparently manages to maintain a rudimentary jet, about 500 km in length, which separates from a point just south of the latitude of zero wind stress curl. The two isopycnic model simulations (Fig.6b,c), on the other hand, show no sign of a free jet at all, but reveal a chaotic eddy field whose only quasi-permanent features appear to be the southward flowing western boundary current extending far into the anticyclonic wind stress domain and the notorious anticyclonic eddy in the southern gyre region close to the western boundary. .s .test page 9 .p $7. Summary$ .p The main results of the model intercomparison described above can be summarized as follows. .p (a) There is close agreement between the solutions obtained with the quasi-isopycnic BB model and the (numerically dissimilar) 'pure' isopycnic model. This suggests that circulation features found in the quasi-isopycnic model solutions but not in the solutions obtained by HLH are caused by the choice of geophysical parameters and vertical model resolution, rather than by numerical peculiarities inherent in the hybrid coordinate treatment. .p (b) Solutions obtained with the isobaric (i.e., essentially Cartesian) model version bear closer resemblance to the HLH two-layer solutions than do the isopycnic model solutions. This is consistent with our finding that the currents in the isobaric model require more time to become hydrodynamically unstable and that instability processes in general are more subdued in the isobaric than the isopycnic model version. This, in turn, supports the long-held notion that on a grid point-by-grid point basis material coordinate models have better vertical truncation properties than Cartesian models and thus are more capable of resolving baroclinic structures. .p (c) While the hybrid and the isobaric model show similar kinetic energy levels after reaching statistical equilibrium, the purely isopycnic model is somewhat more energetic than the other two. There are two possible reasons for this: the eddy size in the isopycnic model may be slightly larger than in the other models, or the solutions may be less noisy on the grid scale. In either case, lateral dissipation will be a less effective energy sink. .s 2 .spacing 1 .test page 10 $Acknowledgements:$ We wish to acknowledge very helpful discussions with Prof. Claes Rooth. This work was supported by the National Science Foundation under Grant No. OCE80-26010. Computations were carried out at the National Center for Atmospheric Research in Boulder, Colorado. NCAR is sponsored by the National Science Foundation. .page .spacing 2 .center;R E F E R E N C E S .s 3 .left margin 19 .indent -6 Arakawa, A., 1966: Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. $J.Comput.Phys., 1$, 119-143. .s .indent -6 Bleck, R., 1984: An isentropic coordinate model suitable for lee cyclogenesis simulation. $Riv.Meteorol.Aeronaut., 44$, 189-194. .s .indent -6 ----, 1985: On the conversion between mean and eddy components of potential and kinetic energy in isentropic and isopycnic coordinates. $Dyn.Atmosph.Oceans., 9$, 17-37. .s .indent -6 ----, and D.B. Boudra, 1981: Initial testing of a numerical ocean circulation model using a hybrid (quasi-isopycnic) vertical coordinate. $J.Phys.Oceanogr., 11$, 755-770. .s .indent -6 Boris, J.P., and D.L. Book, 1973: Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works. $J.Comput.Phys., 11$, 38-69. .s .indent -6 Holland, W.R., 1978: The role of mesoscale eddies in the general circulation of the ocean - numerical experiments using a wind-driven quasi-geostrophic model. $J.Phys.Oceanogr., 8$, 363-392. .s .indent -6 ----, and L.B. Lin, 1975a: On the generation of mesoscale eddies and their contribution to the oceanic general circulation. I: A preliminary numerical experiment. $J.Phys.Oceanogr., 5$, 642-657. .s .indent -6 ----, and ----, 1975b: On the generation of mesoscale eddies and their contribution to the oceanic general circulation. II: A parameter study. $J.Phys.Oceanogr., 5$, 658-669. .s .indent -6 Lighthill, M.J., 1969: Linear theory of long waves in a horizontally stratified ocean of uniform depth (Appendix). $Phil.Trans.Roy.Soc, 265$, 85-92. .s .indent -6 Lorenz, E.N., 1955: Available potential energy and the maintenance of the general circulation. $Tellus, 7$, 157-167. .s .indent -6 ----, 1967: $The nature and theory of the general circulation of the atmosphere$. World Meteorol.Org., Geneva, 161 pp. .s .indent -6 Redi, M.H., 1982: Oceanic isopycnal mixing by coordinate rotation. $J.Phys.Oceanogr., 12$, 1154-1158. .s .indent -6 Sadourny, R., 1975: The dynamics of finite-difference models of the shallow water equations. $J.Atmos.Sci., 32,$ 680-689. .s .indent -6 Shapiro, M.A., 1975: Simulation of upper-level frontogenesis with a 20-level isentropic coordinate primitive equation model. $Mon.Wea.Rev., 103$, 591-604. .s .indent -6 Zalesak, S.T., 1979: Fully multidimensional flux-corrected transport algorithms for fluids. $J.Comput.Phys., 31$, 335-362. .page .spacing 2 .center;FIGURE LEGENDS .s 3 Figure 1: Mean flow kinetic energy (based on vertically integrated 360-day averages) plotted against time for the three model versions. .s 2 Figure 2: Basin averages of K|M-to-P conversion (top), K|M-to-K|E conversion (center) and P-to-K|E conversion (bottom) plotted against time for the three model versions. .s 2 Figure 3: Instantaneous near-surface velocity fields (0-300 dbar layer averages) from the isobaric model for days 750, 760, 770, and 780. Triangular pointer marks latitude of maximum wind stress. Only part of the grid domain is shown. .s 2 Figure 4: Areal distribution of eddy (left) and mean field (right) kinetic energy for the isobaric (top), hybrid (center), and isopycnic (bottom) model version during model year 2. Thin dotted lines: mean barotropic flow stream function. Heavy dotted line: location of maximum eastward wind stress. Only part of the grid domain is shown. .s 2 Figure 5: As in Fig.4, but for model year 3. .s 2 Figure 6: Near-surface velocity fields (0-300 dbar averages) at day 1800 from the isobaric (left), hybrid (center), and isopycnic model (right). Triangular pointer marks latitude of maximum wind stress. Only part of the grid domain is shown.