Potential vorticity dynamics in the framework of disk shallow-water theory: II. Mixed Barotropic-Baroclinic Instability

# Potential vorticity dynamics in the framework of disk shallow-water theory: II. Mixed Barotropic-Baroclinic Instability

O. M. Umurhan 1-31-3
###### Key Words.:
Hydrodynamics, Astrophysical Disks – theory, instabilities
offprints: O.M. Umurhan 11institutetext: School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, U.K. 22institutetext: School of Natural Sciences, University of California Merced, Merced, CA 95343, USA 33institutetext: Astronomy Department, City College of San Francisco, San Francisco, CA 94112, USA
###### Abstract

Context:Potential vorticity dynamics for astrophysical disks.

Aims:Extend exploration of these instabilities in cold astrophysical disks which are three dimensional whose mean states are baroclinic. In particular, we seek to demonstrate the potential existence of traditional baroclinic instabilities of meteorological studies. We show this in a simplified two-layer Philips Disk Model of a midplane symmetric disk.

Methods:We consider the dynamical normal-mode response of thin annular disks with two potential vorticity defects, one in each of two disk layers. Each disk layer is of constant but differing densities. The resulting mean azimuthal velocity profile shows a variation in the vertical direction implying that the system is baroclinic in the mean state. The stability of the system is treated in the context of disk shallow water theory wherein azimuthal disturbances are much longer than the corresponding radial or vertical scales. The normal-mode problem is solved numerically using two different methods.

Results:The results of a symmetric single layer barotropic model is considered and it is found that instability persists for models in which the potential vorticity profiles are not symmetric, consistent with previous results. The instaiblity is interpreted in terms of interacting Rossby waves. For a two layer system in which the flow is fundamentally baroclinic we report here that instability takes on the form of mixed barotropic-baroclinic type: instability occurs but it qualitatively follows the pattern of instability found in the barotropic models. Instability arises because of the phase locking and interaction of the Rossby waves between the two layers. The strength of the instability weakens as the density contrast between layers increases.

Conclusions:Baroclinic instability is feasible for astrophysical disks but has the character of mixed barotropic-baroclinic type. The instability as explored in this study could be present in protoplanetary disks with weak vertical density stratification.

## 1 Introduction

The phenomenon of baroclinic instability has been studied with intensity for over 70 years. Sometimes referred to as sloping convection (Vallis 2006), it is an instability that relies on the presence of a vertical shear in a zonal flow in a rotating atmosphere. This archetypical sheared flow exists as a result of a thermal-wind relationship balancing Coriolis effects with meridional (equator-to-pole) temperature gradients. It is ultimately a reflection of the mismatch in the isobars and isopycnals in an atmosphere. It is alternatively interpreted as an instability resulting from the interaction of inter-layer Rossby waves (Hoskins et al. 1985, Baines and Mitsudera 1994, Heifetz 2004). Baroclinic instability is important to a variety of geophysical and astrophysical phenonomena. Among others, it plays a well-known role in the generation of terrestrial weather as well as in the dynamics of the Ocean Mixed Layer (Haine and Marshall 1998). It has been conjectured to play an important role in the mixing processes in stars (Knobloch and Spruit 1982, Fujimoto 1988).

The original studies of the instability (e.g. Charney 1947, Eady 1949, Philips 1954) considered the process in the small Rossby number limit in which the Rossby number, Ro, is given by , where the terms are (respectively) the typical speed of a vortex, planetary rotation rate (at latitude) and the typical horizontal length scale of the flow. Baroclinic instability is a process that figures prominently in the geostrophic limit (Ro ) but theoretical considerations have shown that it persists relatively unabated even when Ro is (Stone 1970, Stone 1971, Molemaker et al. 2005). These studies also argue that the instability also persists for significant deviations from hydrostasy.

The possibility that baroclinic instability is present in accretion disk flows, which are characterized by Ro , was examined in the studies of Cabot (1984) and Knobloch and Spruit (1985,1986). Owing to the fundamental inseparability of the resulting linearized problem (even in the relatively simplified Boussinesq limit) Knobloch and Spruit (1986) were only able to conclude that the instability is feasible for radial and vertical disturbances comparable to the local scale height of the disk. The feature complicating the analysis and rendering it inseparable is the dominance of the Keplerian shear making traditionally employed small Ro analyses of meteorological studies here unusable.

Adding to the challenge: because the Keplerian shear is barotropic the resulting problem for a disk is fundamentally of mixed baroclinic-barotropic type. This type of mixed problem has been examined in the context of quasigeostrophic flows by McIntyre (1970), Ioannou & Lindzen (1986), James (1987). These studies generally agree that the structure of baroclinic waves will get increasingly confined in the shearwise direction with increasing strength of the barotropic shear. The confinement, due to the barotropic governor mechanism (James 1987), appears to have the effect of weakening growth rates.

Keplerian disks are notoriously difficult to finesse into a geostrophic-hydrostatic framework familiar in traditional meteorology. Knobloch and Spruit (1985,1986) note that it is only in the long stremwise wavelength limit (azimuthal in disk geometry) of a Keplerian flow will some of the formalism of quasigeostrophy carry over. Umurhan (2008) developed a scaling analysis in which a local annular disk section can be viewed in a semi-geostrophic, quasi-hydrostatic framework (see also Balmforth and Spiegel 1996). This, in turn, permits a dynamical analysis purely in terms of the potential vorticity, typical of meteorological studies. This construction was used to re-examine the Rossby wave instability (Li et al. 2000, Meheut et al. 2010) by illustrating it in a setting stripped down to its essential components (Umurhan 2010).

We use the limit developed in Umurhan (2008) to examine baroclinic instability in an annular section of a cold accretion disk. We circumvent the difficulties encountered in Knobloch and Spruit’s calculation by constructing a two-layer Philips model (Philips 1954). The Philips model is an effective simplification of the governing equations that allows for analytical tractability without losing the underlying physics at play in the baroclinic instability (see an extended discussion of this in the text of Vallis 2006).

This work is organized as follows. In Section 2 we develop a two-layer Philips disk model from the disk shallow water equations found in Umurhan (2008) utilyzing a number of simplifications and assumptions. In Section 3 we set up the general normal mode problem for this two-layer, mid-plane symmetric, model system including the boundary conditions and numerical methods employed. We describe the model whose stability we test including a description of the features we include which is designed to strip away the complications arising from critical layers. In Section 4 we consider the stability properties of a purely barotropic (single-layer) configuration and interpret the results from the perspective of interacting Rossby waves developed in previous studies. In Section 5 we build upon the picture developed for the barotropic calculation and describe a baroclinic model configuration and examine its stability properties. In Section 6 we offer a summary of our principle results and a short litany of reflections and remarks.

## 2 Derivation of the Two-Layer Philips Disk Model

The analysis of this study assumes dynamics occur on azimuthal length scales which are much longer than the system’s vertical and radial scales. The equations appropriate for such a thin annular setting were developed from a scaling analysis detailed in Umurhan (2008). However they can be found in various other forms in older studies as well, e.g. Knobloch and Spruit (1985,1986) and Balmforth and Spiegel (1996).

These equations describe dynamics in a thin annular section of a disk centered at some radius rotating with the local Keplerian rate . If the local soundspeed of the disk is given by then the ratio forms a natural small parameter of the disk equations. Formally this means we will take . The equations are therefore derived by assuming a number of scalings, namely that the azimuthal velocity scales as and that the vertical and radial velocities scale as . Together with the assumption that the radial and vertical scales are a factor of smaller than the azimuthal scale , the velocity assumptions here imply a long time scale . It is important to note that this long timescale necessarily filters out acoustic and gravity wave type motions from these dynamics. The combination of these scalings into the fundamental equations of motion results in the following set of reduced non-dimensionalized equations

 −2Ω0V = −1ρ∂xP−qΩ20x, (1) dVdt+2Ω0u = −1ρ∂yP, (2) 0 = −1ρ∂xP+Ω20z, (3)

together with

 1ρdρdt+∇⋅u = 0, (4) ρTdSdt = Q, (5)

in which

 ddt=u∂x+V∂y+w∂z.

is the pressure, where and respectively denote density and temperature. The velocity vector is comprised of components in the radial (), azimuthal () and disk vertical () directions respectively. The lengths and represent the radial and vertical scales scaled by while the azimuthal scales being scaled by . The entropy is where is the ratio of specific heats. represents the non-dimensionalization of the usual Coriolis terms that emerge when one moves into a rotating frame. In this formulation it is formally ‘’ but we keep it around for the sake of latter clarity. For Keplerian disk systems .

The lack of the usual inertial terms in the resulting radial and vertical momentum equations expresses the extremeness of the scalings describing the thin annular disk region. The resulting equations describe dynamics which are vertically hydrostatic but, moreover, geostrophic in the radial direction. In meteorological studies this is sometimes referred to as semi-geostrophic because the azimuthal flow dynamics are not in fundamental geostrophic balance.

To derive a Philips model that is tractable and transparent the following assumptions are made:

1. Dynamics are symmetric about the midplane ,

2. Two layers are examined with differing uniform densities (and the configuration is stable to buoyant instabilities),

3. The flow in each layer is incompressible.

The assumption of incompressibility means that the continuity equation is replaced, instead, by

 ∇⋅u=0, (6)

in each layer. A pictoral representation of the system envisioned for this analysis is shown in Figure 1. From here on forth, quantities in the upper layer are designated with a subscript ‘’ while in the lower layer they are represented by the subscript ‘’. Thus, for example, the uniform density of the upper layer is while in the lower layer it is . Because incompressibility has been assumed for each layer there is no longer any need for the energy equation (5). Instead, the evolution of each layer’s thickness will be followed (see below).

A qualitative inspection of the state of affairs depicted in Figure 1 shows that, indeed, the isopycnals and isobars of the atmosphere are misaligned which means that the basic state will be baroclinic with the in-layer azimuthal velocities being different from one another. The misalignment has been placed here by hand and is controlled by the functional form of the heights and . The origins of their respective distributions ultimately rests upon radially dependent thermal process involving radiation and disk chemistry requiring detailed calculations (e.g. Turner et al. 2012) without making the simplifying assumptions we have made. For the purpose of this model, how the atmosphere got arranged into this configuration is not our concern: it is enough that the profile is in steady state.

The azimuthal velocity in both layers are written as a sum of the local Keplerian part plus a deviation where is either ‘’ or ‘’ as in, . Assuming that the pressure of the disk is zero at its top, the vertical momentum equation may be immediately integrated to yield,

 P=Ω2012⎧⎪ ⎪⎨⎪ ⎪⎩ρ0(h20−z2);h1≤z≤h0,ρ0(h20−h21)+ρ1(h21−z2);0≤z≤h1. (7)

Putting this solution for the pressure back into the horizontal momentum equations shows that for each layer

 −2Ω0vi = −∂xΠi, (8) (DDt)ivi+Ω0(2−q)ui = −∂yΠi, (9)

where

 Π0=12Ω20h20,Π1=12Ω20[δh20+(1−δ)h21]. (10)

The quantity may be viewed as the analog of the streamfunction in meteorological studies. Each layer height, which is in effect an integration constant, is designated by . The parameter measures the ratio of the two densities. Henceforth, since interest is in configurations which are buoyantly stable, attention will be given only to those values of .

Inspection reveals that the resulting horizontal momentum equations are independent of the vertical coordinate. As such the time derivative operator in each layer has been written as

 (DDt)i=∂t+ui∂x+(vi−qΩ0x)∂y.

The incompressiblity equations may be integrated in each layer revealing

 w=⎧⎪⎨⎪⎩−z(∂xu0+∂yv0)+¯W0;h1≤z≤h0,−z(∂xu1+∂yv1);0≤z≤h1. (11)

The symmetry condition has been imposed upon the vertical velocity profile by enforcing that its value is zero at . The velocity term is determined by enforcing certain matching conditions when following the motion of the individual interfaces. Since plays no essential role in the subsequent linear theory we do not detail its form here. The motion of each layer edge is viewed from both layers and is accordingly set equal to the vertical velocities found in (11). Combining this with the criterion that there is no layer separation reveals that the dynamics of each layer’s thickness is governed by

 (DDt)1h1=−(∂xu1+∂yv1)h1, (12)

and

 (DDt)0(h0−h1)=−(∂xu0+∂yv0)(h0−h1). (13)

The above equations may be combined and simplified following standard procedures resulting in

 (DQ0Dt)0=0,Q0≡Ω0(2−q)+∂xv0h0−h1, (14)

in the upper layer while, for the lower layer,

 (DQ1Dt)1=0,Q1≡Ω0(2−q)+∂xv1h1. (15)

The above relate the fact that the potential vorticities in each layer, represented by and , are materially conserved by the flow. It should be noted here that even though the are conserved in each of their respective layers, because the same involve quantities that develop in layers outside of their own means that layers are coupled to one another.

Equations (14-15) together with (8-10) completely describe the model system. It is referred to hereafter as the two-layer Philips disk model.

## 3 Steady states, their perturbations and numerical method

As the purpose of this work is to demonstrate the possible normal-mode responses for various configurations of this two-layer Philips model, in this section we develop the formalism for the perturbative study of generalized steady state configurations. This steady states are assumed to be radially dependent, however, perturbations of these states are non-axisymmetric. In the subsequent two sections we will detail the normal mode response for a variety of specific steady state profiles.

For the sake of transparency in the following presentation shall be set to its nominal value of . We study perturbations of already laid down potential vorticity profiles. These steady profiles describe mean azimuthal velocities which are deviations around the Keplerian base flow. There are no mean radial velocities. This means in effect, that the mean PV profiles will depend only upon the radial coordinate . Therefore, assuming and it follows from the relationships between ’s and found in (14-15) that

 2(2−q)+∂2x¯Π0−2¯Q0(¯h0−¯h1) = 0, (16) 2(2−q)+∂2x¯Π1−2¯Q1¯h1 = 0, (17)

in which the relationships between and , found in (10), have been inverted:

 ¯h0=(2¯Π0)1/2,¯h1=21/2(¯Π1−δ¯Π01−δ)1/2; (18)

expressing the mean steady height’s in terms of the steady in-layer pressures . In deriving the above expressions the following relationships have also been explicitly used,

 ¯v0=12∂2x¯Π0,¯v1=12∂2x¯Π1. (19)

Linear perturbations of the Philips model system developed by writing for all dependent quantities ‘

 Fi→¯Fi(x)+F′i(x,y,t) (20)

where the primes denote perturbations. Equations (14-15) and (8-10) are explicitly written out

 (∂t+¯v0∂y−qx∂y)Q′0+u′0∂x¯Q0 = 0, (21) (∂t+¯v1∂y−qx∂y)Q′1+u′1∂x¯Q1 = 0, (22)

in which

 (¯h0−¯h1)Q′0 = ∂xv′0−¯Q0(h′0−h′1), (23) ¯h1Q′1 = ∂xv′1−¯Q1h′1, (24)

together with

 2v′0 = ∂xΠ′0, (25) 2v′1 = ∂xΠ′1, (26) h′0 = (2¯Π0)−12Π′0, (27) h′1 = 1212(11−δ)12(¯Π1−δΠ0)−12(Π′1−δΠ′0). (28)

The linear equations resulting from perturbations of (9) are rewritten in order to express the perturbation radial velocity, , in terms of the other perturbation quantities,

 u′0 = −1(¯h0−¯h1)¯Q0{∂yΠ′0+(∂t+¯v0∂y−qx∂y)v′0}, (29) u′1 = −1¯h1¯Q1{∂yΠ′1+(∂t+¯v1∂y−qx∂y)v′1}. (30)

Inserting (23-26) and (29-30) into (21-22) reduces the system to two coupled PDE’s for the variables . This is to say,

 (31)

and

 −2∂yΠ′1⋅Dln¯Q1=0, (32)

where the equations relating the perturbation heights to the perturbation pressures are given in (27-28). The symbol ‘’ denotes differentiation with respect to .

All perturbation quantities are assumed to be periodic in the azimuthal direction and, furthermore, we assume normal mode temporal structure. That is to say, the solution ansatz we adopt is

 f′=^f(x)eik(y−ct)+c.c., (33)

where is the azimuthal direction wavenumber and is the wavespeed. The goal will be to find the admitted values of for the variety of configurations examined in the following sections. Values in which Im correspond to instability.

For the radial boundary conditions it will be assumed that the fluctuation perturbation pressures go to zero at symmetrically placed radial boundaries, i.e. that as . These boundaries will be placed significantly far from the region where there are nontrivial variations in the mean flow. Typical values taken are (unless otherwise specified). We are concerned with uncovering effects which are insensitive to the effects introduced by using such an artificial boundary condition like we have here. As such, in all instances we have verified that the growth rates we determine are in fact insensitive to the position of the boundaries themselves by recalculating the growth rates for several values of significantly greater than . Because all eigenfunctions we generate show exponential decay as one approaches these artificial boundaries, we are confident that there does exist, in fact, a minimum value of beyond which the effects of the boundaries are unimportant in determining the dynamics we study. (Except, however, see the special case discussed at the end of Section 5).

The coupled linearized ordinary differential equations are solved using a two-step procedure. The x-domain is discretized on a Chebyshev grid . The number of grid points is . The solutions are assumed to be represented by the corresponding Chebyshev decomposition on this grid:

1. Step 1. This step involves expressing the governing equations (31-32) and their operators in physical space. For a solution to represented in physical space by the two coupled equations form a single matrix equation

 −icMΠ=LΠ, (34)

where

 Π=[Π(1)0,Π(2)0,⋯,Π(N−1)0,Π(N)0,Π(1)1,⋯,Π(N)1]T,

in which the superscript denotes the vector transpose operation. is a vector of length while and are matrices. and are matrices containing the Chebyshev representation of the linear operators of the problem. The boundary condition that is zero at the endpoints are implemented the following way. The rows corresponding to the boundaries are replace by the equation

 ∂tΠ(1,N)(0,1)=−ϖΠ(1,N)(0,1),

where we choose to be a very large number (typically ). We then determine the full eigenvalue set, using standard packages for determining numerical eigenvalues as found in Matlab. We note that the four eigenvalues of are discarded since they are unphysical.111The normal modes associated with these eigenvalues physically represents a very fast relaxation of the system to the boundary condition. This scheme is used widely in the fluid dynamics literature. The remaining eigenvalues and corresponding eigenvectors represent the responses of the system. The boundary conditions are satisfied with an error of . We assume a value of which typically yields converged solutions.

2. Step 2 We verify and further refine the solutions we determine in Step 1 by solving (34) using a standard Newton-Raphson-Kantorovich (NRK) method. We input into the procedure a solution vector together with its corresponding eigenvalue . The error tolerance required of NRK solution was that it is correct to one part in .

All quoted converged solutions in this study passed both of the outlined steps above. Furthermore, we checked resolution convergence using the NRK method by taking a given solution determined on a grid of size , interpolate it onto a grid of size , and then using this resulting interpolated solution as an initial guess for the NRK method. In rare instances we found solutions determined via Step 1 to have false unstable eigenvalues, i.e. Im. These were identified using the grid refinement procedure: for those spurious eigenvalues, Im asymptotically as .

## 4 Barotropic configurations

In order to anchor our intuition we begin this analysis by reconsidering the dynamics in a single-layer configuration. Under those conditions the shear is barotropic. The linearized equations we will need to solve are those given in the previous section with set to zero. With the upper layer has no dynamical effect upon the lower layer and, as such, the resulting dynamics describe flow of the lower layer entirely isolated from the upper layer. Therefore, the equations for the steady state are (17) with given and , as given in (19). For the linearized dynamics we solve only (32) with (28) to generate solutions for the lower layer perturbation enthalpy . We also assume the normal mode ansatz (33) and enforce that .

Figure (2) depicts the PV profile we shall examine here. We refer to it as the ‘double-tanh’ profile and is given by the form

 ¯Q1 = A1[1+tanh(Δ−xξ)]+ (35) B1[1+tanh(x−Δξ)]+C1,

where the constants are specified below. This functional form roughly represents a two-step PV-profile with the steps located at where is some parameter. We will refer to twice this parameter as quantifying the width of this type of PV profile. The parameter represents the width of the transition in moving from one region of PV to the other. In the limit the transition from one region of PV to the other (e.g. near in Figure 2) is proportional to a Heaviside step function. One one hand we numerically we try to avoid the limit, yet on the other hand we want to have be small in order to relate the results of the model to our simplified understanding of the system. As such we typically choose . Thus, as depicted in Figure 2, the limit makes the PV-profile look like two successive step functions or ‘defects’ in the PV. We interpret the length as quantifying the separation of these nominal defects of the PV.

We must now relate the constants to physically intuitive features of the model. Because we envisage a PV profile that has an induced velocity profile (as a deviation from the background Keplerian flow, i.e. ) that is zero as , inspection of (15) or (17) requires of us that

 limx→±∞¯Q1=2−qh1(±∞). (36)

In other words, we design this profile to yield constant mean heights, , in the limit of . Of course, these heights may be different from one another but the constant asymptotic form of it ensures that the mean induced flow correspondingly asymptotes to zero as . As such this implies that

 2A1+C1 = 2−qh1(−∞), 2B1+C1 = 2−qh1(∞).

To complete this description we suppose that depth of the depression in , located in the region , to be some factor of the asymptotic value of the PV as , i.e. . This amounts to

 2A1−C1=ϕ¯Q1(−∞)=ϕ2−qh1(−∞).

In summary, then, the barotropic model (35) is characterized by the five parameters: the asymptotic heights , the symmetric positioning of the steps from , , the relative depth of the PV-dip in the middle-zone of this model, , and the width of the transition zones . Note that the middle-zone has a width equal to . The normal-mode system is additionally described by the shear parameter .

A model similar to this one was considered in Umurhan (2010). Unlike here, however, the assumed PV-profile adopted in that previous study were (symmetric) step functions chosen to facilitate analytical tractability. In comparison to the model considered here, the step function used in Umurhan (2010) would correspond to the limit of this model in which . Note also that the symmetry of the profiles studied in Umurhan (2010) is recovered here by choosing . We have therefore introduced the parameter to control the thickness of the transition from one zone to the next.

The rationale of introducing this feature into the model is to cleanly represent the propagation and interaction of Rossby waves without enduring the complexities of critical layers. Because Rossby waves represent disturbances that propagate in regions where the PV changes, setting to be small ensures that the waves are localized to those places where the PV changes the greatest. In this model, for example, there will be two Rossby waves, each of which are localized in at with a radial width of scale . We typically use values of and this usually represented the Rossby wave dynamics without the appearance of critical layers. The robustness some solutions was evaluated by also generating solutions with values of (which necessitated using a higher resolution grid in order to resolve the transition zone with at least 5-7 grid points). Critical layers tended to appear for values of . However, we do not examine these situations in this study.

The barotropic incarnation of the Rossby wave instability was interpreted in Umurhan (2010) to arise as a result of the interaction of the two Rossby waves propagating along the azimuthal directions at the two locations, consistent with the more generic picture of interacting edgewave disturbances (Baines & Mitsudera 1994). An analogous form of this dynamic is present in this model and is especially pronounced for . For example, in reference to the PV profile shown in Figure 2, we observe that the jump in the PV is positive at the step located at . Therefore, in the local frame of reference of the flow at and when taken in isolation, the supported Rossby wave propagates in the negative direction. Similarly, since the average change in the PV profile at is negative, it can support a Rossby wave propagating in the positive direction when viewed from the local reference frame of the flow at . When the laboratory measured wavespeeds of both waves are nearly equal then instability may manifest itself, i.e. when the Hayashi-Young criterion is satisfied (Hayashi & Young 1987). Instability exists when the Rossby waves are essentially counterpropagating along the flow (Heifetz 1999) because it is only when the local propagating tendency of the Rossby wave is against the mean flow that the Hayashi-Young criteria of instability can be met.

In the next section we will also refer to an effective surface density. We define this here as

 Σ1(x)=∫h1(x)0ρ1dz=ρ1h1(x). (37)

As can be seen is equivalent to in the single layer model, but will not be so in the two layer model. Henceforth we take .

A representative summary of the results of this section is displayed in Figure 3. Even though we solve for the wavespeed we depict the growth rates Im, together with having here (and henceforth) set .222Disk-shallow water theory is such that the streamwise wavenumber drops out of all calculations owing to it being, essentially, a long-wavelenth theory. See also Knobloch & Spruit (1985,1986). The growth rates are quoted in time units used to derive the disk shallow water equations (see ), namely, time is scaled in the longtime units .

The results we plot here are an extension of the symmetric models considered in Umurhan (2010). Here we see that that the models do not have to display symmetric PV profiles in order for instability to surface and that, consistent with the results of Li et al. (2000), asymmetric PV profiles can just as easily lead to instability. We find that for a given set of layer parameters and (the depth of the PV depression) - as the background shear is lessened, the range of separation values that admits instability is shifted toward larger values. In other words, weakening the shear means that profiles with steps further separated from one another will result in instability.

In terms of the picture of interacting Rossby waves, this is easily rationalized in the way outlined earlier in this section: (i) for instability to occur the laboratory measured phase speed of the Rossby waves along each step must be nearly the same, (ii) since the phase speeds of the Rossby waves are governed both by the background flow speed at the location of the step plus the intrinsic Rossby wave speed (as measured in the moving frame at the step) induced by the local change in PV at the transition region it follows that, (iii) with all else held equal, a reduction in the background shear requires the location of the defects at to be set further apart so that the laboratory measured phase speeds may be, once again, nearly equal.

When the step profiles are symmetric, i.e. when instability may occur for separation consistent with previous results in Umurhan 2010. However, as begins to deviate away from the value then there exists a minimum value for the defect separation above which instability may occur. It is also interesting to note that there is a maximum disparity of heights beyond which there is no instability possible. E.g. for values of with there is no instability possible in the model in Figure 3. The interpretation of this stability quality is once again similar to the above: varying beyond some critical value results in individual waves whose speeds cannot possibly match each other. In particular, an increase of the speed of the wave at the step occurs because increasing means that the average PV-gradient increases at step which, in turn, boosts the (in frame) counterpropagating Rossby wavespeed at the step.

The insensitivity of these results on has been verified. It is important to keep in mind that the growth rate values convergence to the values quoted once is greater than about 5. For values of the growth rates weaken. The presence of channel walls in problems like these creates the effect of ’image’ waves (like image charges in electrostatics problems, Drazin & Reid 1981) which act to inhibit phase locking. This behavior has been detailed in a study of the Rayleigh problem with variable wall positions (Heifetz et al. 2009) and the trend for instability suppression with more narrowly set walls is consistent with their findings.

The same character follows in the baroclinic cases studied except for one instance discussed later. Finally, consistent with the Howard Semicircle Theorem, we observe that as the shear is weakened the maximum value of the growth rate is reduced.

## 5 Baroclinic profiles

In this section we solve the two layer problem. We assume that each layer has a single jump in PV but that their positions in are, in general, different from one another. We employ the basic type of model encountered in the previous section. In particular we suppose that the mean PV profiles

 ¯Q0 = C0+A0[1−tanh(x−x0ξ)], ¯Q1 = C1+A1[1−tanh(x−x1ξ)]. (38)

The steps for each layer are located at respectively. Without loss of generality we place these positions symmetrically about the point and assign the position values and . is a parameter as before, however, can take on all values (and is not restricted to being strictly positive). As before, the thickness of the step in each layer is governed by the parameter .

Following the arguments laid out in the previous section, we require of the steady PV states to correspond to zero induced azimuthal flow as . Inspection of (14-15) indicates that the constants must simultaneously satisfy

 C0+2A0=2−qh0(−∞)−h1(−∞), C0=2−qh0(∞)−h1(∞) C1+2A1=2−qh1(−∞),C1=2−qh1(∞),

where the parameters are the asymptotic values of the layer heights in the appropriate limit of . The only requirement we place upon the these are that and .

We solve for the steady configuration (16-19) using the above model for the PV profiles, which are governed by the six parameters . The parameters of the steady configuration are also determined by the three remaining other parameters: the shear , the density ratio and the domain size . A qualitative example of the resulting states is shown in Figure 4. It is evident from the resulting flow that the mean deviation velocity fields are different from one another between layers indicating that the resulting flow state is baroclinic.

Figure 5 depicts the character typifying instability in these models. The parameters chosen here were meant to motivate maximum comparison to the results shown in the top graph of Figure 3. The top layer has a density ratio of the bottom layer. The results are qualitatively similar to the behavior of the barotropic model with certain modification to the region and degree of instability. Since, by design, there are only two waves which can interact - one emanating from the step in the upper layer and the other from the step in the lower layer - the interpretation of the instability is the same as for the barotropic case. Thus, baroclinic instability in this kind of baroclinic problem closely shadows the character of the instability in the barotropic case. The fundamental difference is, of course, the fact that the source of the Rossby waves are from different vertical layers.

The coupling between layers is governed by the density contrast . In Figure 6 we display how the growth rates are modified for a model in which we vary the separation and the density contrast . Indeed, as one might expect, both the occurrence of instability and its corresponding growth rate diminish as . When viewed from the perspective of a normal mode problem, the lower layer essentially experiences no dynamical influence of the upper layer when the mass of the upper layer is vanishingly small (the lower layer “knows” of the upper layer due to the pressure fluctuations the upper layer induces upon the lower layer). From the mathematical viewpoint of a normal mode problem, both layers decouple from one another when (see the previous section). Nonetheless even for the Rossby waves in each layer are finite (and non-zero) and propagate in isolation. Thus when is increased the bottom layer feels the upper layer in proportion to while the upper layer feels the bottom layer by the latter’s influence upon the upper layer’s thickness. 333We note in passing that when treated as an initial value problem the situation is different for nearly zero. The lower layer still dynamically evolves free of the influence of the upper layer but the upper layer’s dynamics are influenced by the motion of the lower layer as it shows up in the upper layer’s dynamical evolution as an external forcing term. We do not take up this issue in this study. We also note here that the figures show that growth rates persist as but we are less inclined to infer very much about this limit because as as the fundamental scaling relations leading to the Philips disk model begins to cause a break in the scaling assumptions behind the disk shallow water equations. Nevertheless the propensity for instability in these baroclinic models for not very close to zero suggests that baroclinic instability may be manifest in disks which are weakly stratified. We elaborate more upon this in Section 6.

For a given set of model parameter values, varying the shear acts to reduce the growth rate of the instability and to push the unstable range of step separations toward larger values. In Figure 7 we depict the growth rate trends for for variable values of and . Although the pattern of instability in this parameter diagram is more complex than the previous graphs, we see clearly that as the shear weakens (holding the separation fixed) the instability generally goes away. On the other end of the spectrum, as the Rayleigh limit is approached, i.e. , the instability also appears to become suppressed. However we note here that the scaling ordering that is used to develop the disk shallow water equations begins to breakdown as approaches (see Umurhan 2008). We therefore consider these observed trends (in this particular limit) with more caution. A proper investigation of what happens in the Rayleigh limit should involve re-examining the equations of motion in that particular limit.

For given values of the shear there appears to be no parameter configuration that we could find in which there is both instability with the jets lying exactly on top of each other () and where the instability does not vanish as . Such a configuration might correspond to a flow profile in which the jets model a baroclinically unstable setup relevant to more “traditional” baroclinic flow problems (like for the Earth and other Solar System planets). As we noted earlier, all of the results we report in this study are those in which the growth rates remain unaltered as the artificial domain walls in are made very large . In Figure 8 we show an example set of parameters in which there is instability but where the separation in of the two PV jumps is zero. The values of the top of the atmosphere are represented in a way such that the ratio of remains fixed but the overall height of the atmosphere can vary with the parameter . For the example depicted, and . Thus acts like a parameters varying the lid of the atmosphere. The values of the lower layer parameters are fixed. As such Figure 8 displays the growth rates as a function of and with . Instability appears to exist but for values of so large that the atmosphere is effectively thin and tall. When we then increase the value of (holding fixed the remaining atmosphere parameters) the instability eventually vanishes.

We therefore conclude that the possibility of instability in a model in which the jets lie on top of each other, e.g. as implied by the results shown in Figure 8, is an artifact of the narrowness of the radial boundaries and is not likely to be possible for a real protoplanetary disk.

## 6 Summary and Discussion

Baroclinic instability in the guise familiar in meteorology (in which isopycnals and isobars are misaligned) is feasible for cold protoplanetary disks, albeit in a mixed baroclinic-barotropic form owing to the strength of the underlying Keplerian shear. By using a two-layer Philips model we have further substantiated the conjecture made by Cabot (1984) and Knobloch & Spruit (1984,1985) that baroclinic instability is relevant for accretion disks. Aside from determining growth rates and interpreting the mechanics of the instability in terms of Rossby waves interacting across layers, it is also shown that the instability weakens in strength as the vertical density stratification of the layer is increased. The results of the two-layer model suggests that meteorological-type of baroclinic processes, which wholly depend on the mismatch of isobars and isopycnals in the mean state, is likely for disks that have weak vertical density stratification. Of course, this does not limit the full scope of possibilities and, as such, this must be analyzed more fully in the setting of more realistic disk models.

The use of the simple Philips two-layer model has been shown to be very good at capturing the essence of baroclinic flow dynamics present in continuously layered models of atmospheres (Vallis 2006). Comparison of growth rates determined for both the Eady Model (continuous vertical shear) and the Philips model shows that, although they are not precisely the same, the qualitative behaviors are consistent with one another and it is for this reason why the Philips model is so widely regarded as a useful theoretical platform to examine the physics of these and related processes. In fact, studies have been done in which multiple layer Philips models have been examined for the same problem yielding no significantly different results in classical meteorological investigations. This has partly to do with the geometry of the system being analyzed: in a planetary atmosphere there is usually a bottom boundary and the direction of gravity is the same.

By contrast we have investigated a disk model in which symmetry with respect to the disk midplane is imposed a priori. Thus the dynamics possible here, which allows only for a subset of possible phase relationships and Rossby wave interactions between the midplane layer and the layer above it, is a subset of a wider class of dynamical responses possible. For instance, if the restriction of midplane symmetry were relaxed then there comes into existence the possibility of new types of Rossby wave interactions independently involving three layers (the midplane and the two layers sandwiching the midplane) which will contain the dynamics of the two-layer system studied here. This may yield qualitatively new dynamical possibilities.

We should note here the similarities between traditional ‘meteorological’ baroclinic instabilities to other types of baroclinic processes investigated for both disks (e.g. Klahr & Bodenheimer 2003, Petersen et al. 2007a Petersen et al. 2007b, Lesur & Papaloizou 2010, Lyra & Klahr 2011) and the Sun (Gilman & Fox 1997, Zaqarashvili et al. 2010, to name only a few studies). In these analyses potential vorticity disturbances are effectively baroclinically torqued by some process that is either thermal or magnetic or, in some investigations, a combination of both. In those studies of magnetic Rossby waves in the Sun (e.g. Gilman & Fox 1997, Zaqarashvili et al. 2010) Rossby waves will propagate along those places where there are strong gradients in the solar differential rotation. Treated as a hydrodynamic problem, these Rossby waves are not necessarily unstable on their own, especially for differential rotation profiles inferred to be characteristic of the Sun’s tachochline. (Spiegel & Zahn 1992). However, when a strong toroidal field is added to the mix there appears an effective baroclinic source term (owing to the Lorentz force) which can render the wave unstable. In all of these studies the instabilities are global in character as is the instability investigated in this study since we have examined the fate of disturbances with azimuthal wavelengths on the scale of the disk.

The upcoming Paper III in this series examines an alternative derivation of the Shallow Water Equations for disks that lifts the restrictions placed upon the azimuthal scales as compared to the vertical and radial scales.

## References

• (1994) Baines, P. G., & Mitsudera, H. 1994, J. Fluid Mech., 276, 327
• (1996) Balmforth, N. J., & Spiegel, E. A. 1996, Phys. D, 97, 1
• (1984) Cabot, W. 1984, ApJ, 277, 806
• (1947) Charney, J. G. 1947, J. Meteor., 4, 135
• (1981) Drazin, P.G., & Reid, W.H. 1981 , Hydrodynamic Stability, Cambridge Univ. Press, Cambridge
• (1949) Eady, E. T. 1949 Tellus, 1, 33
• (1988) Fujimoto, M. Y. 1988 A&A, 198, 163
• (1997) Gilman, P. A., & Fox, P.A. 1997 ApJ, 484, 439
• (1999) Godon, P., & Livio, M. 1999, ApJ, 523, 350
• (1998) Haine, T. W. N., & Marshall, J. 1998, J. Phys. Ocean., 28, 634
• (1985) Hayashi, Y.-Y., & Young, W. R. 1987, J. Fluid Mech., 184, 477
• (1999) Heifetz, E., Bishop, C. H., & Alpert, P. 1999 Q. J. R. Meteorol. Soc., 125, 2835
• (2009) Heifetz, E., Harnik, N. & Tamarin, T. 2009 Q. J. R. Meteorol. Soc., 135, 2161
• (1985) Hoskins, B. J., McIntyre, M. E., & Robertson, A. W. 1985 Q. J. R. Meteorol. Soc., 111, 877
• (2003) Klahr, H. H., & Bodenheimer, P. 2003 ApJ, 582, 869
• (1982) Knobloch, E., & Spruit, H. C. 1982 A&A, 113, 261
• (1985) Knobloch, E., & Spruit, H. C. 1986 Geophys. Astrophys. Fluid Dyn., 32, 197
• (1986) Knobloch, E., & Spruit, H. C. 1986 A&A, 166, 359
• (2010) Lesur, G. & Papaloizou, J. C. B. 2010 A& A, 513, 60
• (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000 (LFLC-2000) ApJ, 533, 1023
• (21) Lyra, W. & Klahr, H. 2011 A&A, 527, 138
• (22) Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010 A&A, 516, A31
• (23) Molemaker, M. J., McWilliams, J. C., & Yavneh, I. 2005 J. Phys. Ocean., 35, 1505
• (1987) Pedlosky, J. 1987, Geophysical Fluid Dynamics (2nd ed.), Springer-Verlag, New York
• (2007) Petersen, M. R., Keith J., & Stewart, G. R. 2007a, ApJ 658 1236
• (2007) Petersen, M. R., Keith J., & Stewart, G. R. 2007b, ApJ 658 1252
• (1954) Phillips, N. A. 1954, Tellus,6, 273
• (1992) Spiegel, E. A., & Zahn, J.-P. 1992, A&A 265, 106
• (1970) Stone, P. H. 1970, J. Atm. Sci., 27, 721
• (1971) Stone, P. H. 1971, J. Fluid Mech., 48, 659
• (2012) Turner, N. J., Choukroun, M., Castillo-Rogez J., & Bryden G. 2012, Accepted ApJ arXiv:1110.4166v3[astro:ph]
• (2005) Vallis, G. K. 2006, Atmospheric and Oceanic Fluid Dynamics, Cambridge University Press, New York
• (2008) Umurhan, O. M. 2008, A&A, 489, 953
• (2010) Umurhan, O. M. 2010, A&A, ,521, A25
• (2010) Zaqarashvili, T. V., Carbonell, M., Oliver, R., & Ballester, J. L. 2010 ApJ, 709, 749
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters