Convective "bubbles"

Chuck Doswell

Website version created: 18 May 1998 Updated: 24 September 2006 - minor revisions and moved over from its former site.

Started: 03 October 1995

For a faucet, using Bernoulli's Theorem, the quantity B, where

is a constant along a streamline, where F is the gravitational potential [], p is pressure, r is density and u is the flow velocity. For a faucet, p/r is assumed constant, so the velocity u increases as the potential F decreases to maintain constant B. This means that the flow is speeding up as it falls in the gravitational field. The mass continuity equation then requires, if the flux is to be constant, that the cross sectional area A of the flow must decrease, since mass continuity requires ruA=constant along the flow. It is this narrowing of the flow channel, coupled with surface tension and some other physical effects (according to Shu 1992 pp. 71-74),^[1] that causes drops to pinch off from the steady flow.

For all practical purposes, in the atmosphere, we can regard to be constant. For a convective plume, it would be unthinkable to let p/r be considered constant, however, so the problem is somewhat more complex. Assume a uniform steady inflow, which by mass continuity requires an updraft. During its initial ascent, air from the inflow layer is assumed to be negatively buoyant, such that the buoyancy effect is inhibiting the updraft. Some energy must continually be supplied to the acending parcels if they are to continue their ascent. Ignore for the moment where this energy might come from; let's just say it's there somehow. For the sake of simplicity for the moment, assume that the inflowing air has an LCL and an LFC somewhere above the surface.

Again for the sake of an argument, let's say that the negative buoyancy opposes the forced ascent such that the rising parcels ascend at roughly constant w. The appropriate Bernoulli constant is

so if w is approximately constant, this means that F+p/r is approximately constant along the flow. As the gravitational potential increases (remember, the parcel is rising, not falling), the ratio p/r is decreasing.

Abruptly, the parcel crosses its LFC, at which point its buoyancy becomes positive, perhaps increasingly so. Therefore, the parcel begins to accelerate upward relatively rapidly. Suddenly, along the streamline, the Bernoulli equation says that the pressure must be decreasing very rapidly. If the mass flux is to be constant, as with the dripping faucet, the area of the updraft must decrease. The large increase in upward velocity is going to create substantial horizontal vorticity, especially since the distance across the updraft is decreasing while the velocity difference across it is increasing. This enhances a toroidal vortex ring. The vortex ring advects environmental air into the center, creating a spherical bubble of buoyant air that "pinches off" from the narrowing flow. The steady inflow should create a series of bubbles, analogous to drops of water from a faucet.

If this physical picture is basically correct, then depending on the parameters of the situaion, it may well turn out that even with a smooth, invariant forcing flow on scales above that of the convective cloud, the convection itself will have an inherent unsteadiness. It may even be possible to develop a time scale associated with the frequency of bubble generation. The bubble-development process may be chaotic, as is the dripping of water from a faucet, but that remains to be seen.

The trick will be to develop a model that illustrates the behavior. It might work with a 1 1/2-D or a 2-D model (probably circularly symmetric). A flat, non-rotating earth, of course. The inflow layer should have steady inflow from the sides, with a steady outflow along the boundaries above the inflow that balances the mass inflow. No-slip bottom boundary seems superfluous. There will have to be condensation and the flow nonhydrostatic (probably anelastic). The inflow layer might be forced to stay constant in some way. There should be an LCL below the LFC, and I suspect it would be best to have initial conditions involving signficant CAPE, shear seems superfluous, no worries about condensation physics, so it will be a body-force buoyancy effect perhaps or let the condensate be massless and irrelevant to the dynamics; its only impact is on the thermodynamics .. pure moist adiabatic ascent (not reversible moist adiabatic) above the LCL. Perhaps a tropopause or at least an EL with a strong stable layer at the top (do I want to worry about the damned gravity waves though?); the idea is to let the bubbles pass out of the domain without screwing up the works in the convective layer.

I think an appropriate model should show the effect of "bubbles" from steady forcing on scales larger than the cloud.

Figure 2. Schematic figure showing the grid, as of 10 October 1995

Thoughts - 14 October 1995

The folllowing development is for the "mesoscale" part of the flow (i.e., u=u_m, w=w_m, y=y_m). Unless I have my algebra wrong, the anelastic mass continuity equation in cylindrical coordinates for a fluid in which the density, r_o, varies only as a function of height, z, is

[Strictly speaking, if r_o=r_o(z) alone, then ]. Therefore, we have that

or simply

. (1)

When the flow is nondivergent, then the l.h.s. vanishes, so that

. (2)

Now, define a streamfunction, y, such that

We note that in this system, the component of that is along the tangential direction (normal to the plane of this z-r cross section) is given by

.

This h is positive for a counterclockwise circulation in the r-z plane, which is normal. According to a proper Laplacian in this space, then, we have

which seems a little odd, but the equation actually describing h has the simplicity of a Poisson equation in Cartesian coordinates.

I note that if I specify the inflow velocity, u, at the surface for the radius r=r_o, then the velocity decreases outward in both radial directions from r=r_o. It decreases outward to account for the non-divergence of the flow outside r_o, and it decreases inward to account for the constant convergence assumed inside r_o. The rate of outward decrease has already been noted above in (2). In the regions of constant divergence (convergence below, divergence above), we must have that

.

In the region of convergence at low levels, u is negative while for convergence, D is negative (negative u is increasingly negative as r increases). The opposite is true in the layer of divergence aloft. Note that since u=0 at r=0, when we have specified u(r_o)=u_o, then D = u_o/r_o. The behavior of the u-component is illustrated in Fig. 3 for some convenient numbers

---

Figure 3. Schematic illustration of variation of u-component magnitude as a function of radial distance, r, when the speed at the radius r_o is 3 m s^-1 and when r_o is 3 km; the constant convergence within r_o is 1.0 x 10^-3 s^-1, with non divergence outside of r_o.

It appears that there may still be some messiness in setting up the "mesoscale" flow regime in view of (1) above, involving the variation of density with height. I think ultimately, once the streamfunction is found, this will prove to be the simplest approach.

More thoughts - 16 October 1995

Most of this material is contained in the discussion about "getting pressure." The idea is that the intial sounding (horizontally homogeneous) is in height coordinates, and specifies q and q (mixing ratio) in the lowest layer, to some height h_o. Above that, there is a stable layer with a lapse rate given by g₁ and a dewpoint lapse rate given by g_d1, to a height of h₁. Then there are M significant levels above that, each associated with a change in lapse rate at the m^th level to g_m, and/or a change in the dewpoint lapse rate to g_dm. At the M-1 level, we reach the tropopause, above which the lapse rates stay the same at stratospheric values of g_s and g_ds. This allows quite a bit of flexibility in designing the intial sounding.

However, given these thermodynamic quantities, we need to know the pressure, in order to specify, say, q at every level. The write-up for that is in "finding pressure," and won't be repeated here. The process as coded will involve a first guess at each level and a single correction to an average pressure halfway between the lower level pressure and the first guess. What I have shown in "finding pressure" is that the first correction gets the vast majority of the answer, so it seems superfluous to go beyond the first iteration. This will be relatively easy to code up in FORTRAN and will allow considerable flexibility in the specification of the initial sounding ("getsdg"). Using the spreadsheet, I have developed the following sample of what it might look like (Fig. 4).

Figure 4. Sample sounding showing various significant levels, starting with a constant (q,q) layer in the lowest 1 km, with a lapse rate change at 8 km and a tropopause at 12 km.

---

More thoughts - 20 October 1995

We have the relation that

,

so for regions of nondivergence, it must be that

. (3)

One solution for this is clearly that w is identically zero everywhere in the column, as in the outer region of the model domain. However, if at the top of the boundary layer convergence zone within r=r_o we have w ‚ 0, then the variation of density with height means that there will be a vertical change of w associated with that variation according to the above. Then, above that, the flow will be diverging and w should decrease to zero at the top of the layer. In order for that to happen, the equation to solve will be

. (4)

I believe I can do all this via the spreadsheet since each vertical column will look the same in terms of w. The u-values will look like Fig. 3, above. Once I have the (u,w) values, then the next step is to solve for the density-weighted "mesoscale" streamfunction.

---

Figure 5. Solution of the mass continuity equation where the density follows the indicated curve, derived from the "artificial sounding." This was developed in a spreadsheet using manual iteration to make the mass balance come out in the upper layer to drive w back to zero at the top of that layer.

The mass continuity result shows that the divergence at the top is 0.25485062 times the convergence at the bottom, whereas the depth of the bottom layer is 0.25 times the depth of the top layer. This is the contribution of the density variation, as is the slight increase of w with height in the layer of non-divergence. The r-variation of u with radius in the upper divergent layer looks like the preceding figure (with the appropriate sign, of course), multiplied by the somewhat odd number: 0.25485062. Note that this number results from hand-iterating the divergence value in the spreadsheet "mass continuity" to make the w-value at the bottom of the upper layer come out very small (~10^-10 or so) and slightly positive, with a definite zero assigned to the next higher level.

Next stop ... the streamfunction.

More thoughts - 21 October

The program code is beginning to piece itself together. The processes of (1) adjusting the divergence to complete the w-profile and (2) adjusting the u-component to match the divergence obviously belong in a subroutine ("getmeso"), to develop some flexibility in specifying what the "mesoscale" flow looks like. I might want to change the levels or the divergence value, or whatever. This routine would start with the low-level convergence value, along with the heights of the layers shown in the schematic figure of the model domain, and do all of the mass continuity adjustments to derive the (u,w) field for the"mesoscale" flow.

Looking down the road ... I need to derive the momentum equations and transform to the vorticity equation, including mass continuity along the way. Thermodynamics will be tough. I will be performing a "convective adjustment" and releasing latent heat. I think what I will do is create a "cloud" array that is either 1 or 0 simply to show the outlines of the cloud. The issue of evaporation arises, but I really don't think I want to deal with that ... I may have to in order to have some semblance of realism, but I believe I need to think about that. Also, I probably want to add small perturbations (randomly distributed between 0 and some maximum value) at some locations (involving more randomness) to simulate spatial variation in the heating, so I don't get uniform development of cloud all over the inner radius at the same time.

I also need to find appropriate numerical schemes to use. At the moment, like the idea of centered and space differencing with an Asselin time filter, which looks like

,

for some scalar F, where a = 0.25. I probably will use the 9-point Arakawa Jacobian (if I can find it again) with the stream function for advection. Also, the "sponge" layer will have to be worked out, as well as the lateral (no gradient) boundary conditions to prevent reflected gravity waves.

Another coding thought: start the program with a very short routine that reads a file with all the "tunable" parameters, putting them in COMMON blocks, and then calls the real main program. This would prevent the need to recompile and relink every time I change a parameter. This "toy" model is beginning to look a bit real, now! I have some confidence it will be running in a month. Can I do it?

Thoughts - begun: 27 October

Now that some code has been written to do the preliminaries, let's think about the main deal. Generally speaking, variables all should look like

The details of the mesoscale, ()_m, variables have been provided above. They are hydrostatic, satisfy mass continuity, and the mesoscale flow is time-independent. However, the thermodynamic variables will change with time, as they are acted upon by the combined influence of the mesoscale and perturbation flows.

The total derivatives of the flow in this simulation (axisymmetry and no tangential flow):

,

while the total derivative of q looks like

,

and the anelastic continuity equation implies that there is no local time change of density (except as it expresses itself in buoyancy) - the usual trick, we will be expanding upon below.

---

Derivation of the vorticity equation begins with the momentum equations:

where Coriolis is neglected and the viscous terms (in brackets) only apply in the "sponge" layer. Consider the following: on the right-hand side of each momentum equation, there is a term involving 1/r. This looks like

,

when r'/r_o<<1, a typical assumption. Thus, we have (ignoring the viscous terms, now)

Note that, since because r_o is a function of z alone, we have

Using this, we find that

Expanding the first terms

Thus, we have the difference,

This means that the lhs of the vorticity equation, L, is given by

Substituting for L₁ in the above and doing some minor rearranging, we have

If we make use of the following expansions,

then we find there is a bit of cancellation among terms, and r-derivatives of r_o (and its z-derivative) drop out, so we end up with

We can make use of the definition of our pseudo-vorticity to show that

or solving for the actual vorticity in terms of our pseudo-vorticity

Therefore, the lhs of the pseudo-vorticity equation looks like

This rather ugly looking mess is the price we pay for a simple elliptic equation to solve for the stream function based on this particular form of the pseudo-vorticity.

We have deferred the consideration of the rhs of the equation but now is the time to consider it. We observe that

where it is clear that , so that the rhs of the equation, R, is given by

At this point, we have to consider the component parts of the pressure, p. As with the other variables, it is made up of a mesoscale component, p_m, and a perturbation component, p', where

Therefore, the rhs becomes

Now it is customary in these situations to ignore the effects of the perturbation pressure, and I think I am going to do so here in the interest of avoiding too much complication to the physics. Hence, terms involving the derivatives of p' are going to be neglected, leaving us with

Basically, this means that the primary non-conservative term in the psuedo-vorticity equatoion is a buoyancy term associated with perturbation density gradients. If buoyancy increases with radius (which will be the case for a warm convective bubble), this means that h' is becoming more negative, which is consistent with the normal sign convention on vorticity. Since the prediction equations will be for q and q, a calculation of the perturbation density from these predicted variables will be simple. Putting L=R, and solving for the time tendency of the psuedo-vorticity, we have

The other terms not associated with advection or buoyancy are associated with (a) messiness arising from the fact that the basic state density is factored into the definition of the pseudo-vorticity, and (b) interaction among the velocity derivatives. Let's consider the interaction terms, one by one.

1. : when w' is increasing with r, this means positive vorticity; when u is increasing with r, both terms are positive, but the overall effect is to reduce the vorticity because the increasing u will act to "spread apart" the upward motions, leading to a reduction in the positive vorticity.
2. : this is analogous with the first term, but when u' is increasing with z, the vorticity is negative; spreading apart the shear with increasing upward motion with height reduces the negative vorticity (which is equivalent to an increase in positive vorticity).
3. : when w' is decreasing with z, this is "vertical convergence" that acts to increase the positive vorticity associated with w increasing with r.
4. : when u' increases with r, this is "horizontal divergence" that acts to decrease the negative vorticity associated with w increasing with z.

What about the density terms? Let's consider them, one by one.

1. : since density decreases with height, this term is negative for rising motion and positive pseudo-vorticity. Decreasing density means that the individual terms in the definition of the pseudo-vorticity are decreasing for ascending parcels.

2. : the rate of density decrease decreases with height, so the second derivative is positive. This term is positive for positive vertical motions when the radial velocity perturbation is negative and the converse; i.e., when u' and w are negatively correlated. In effect, when this correlation is negative, relatively large (negative) values of r_ou' are pushed upward, or small values are pushed downward, both leading to positive pseudo-vorticity

3. : the density term, by virtue of being squared, is positive-definite. Now this term is positive when u' and w are positively correlated, which suggests an antisymmetry with the other term. If u' is constant (and negative, as in the inflow layer) with height, then in general rou' decreases in magnitude with height, implying negative psuedo-vorticity, because r_ou' is becoming less negative as one goes up.

---

Figure 6. Schematic showing how the product r_ou' varies with height under conditions where u' is a constant, and illustrating how the contributions to the pseudo-vorticity arise.

Current thoughts on thermodynamics:

There are three prognostic variables: y', q, and q. A "cloud water" pseudo-variable will be carried, but it is really simple because it is zero outside of cloud and unity inside of cloud. Mass continuity will be assured through the use of the streamfunction, I believe. There may be some question about using r_o instead of r' in y' -- this will need to be evaluated. The advection terms will be of the form:

,

where

For the record:

whereas

Thoughts - begun on 05 November

As noted in Pellerin et al. (1995) ^[2] - an appropriate set of thermodynamics equations can be written as

where C is the condensation rate, found at each time step as the difference between the parcel's q-value and its q_s-value, assuming that q > q_s at some point during its ascent; the q_c is the cloud water content. L is the latent heat of condensation, and the ()_o variables represent the initial parcel values. In principle, I could evaporate cloud water via this system, but at the moment, I am undecided whether or not to allow this to happen. At the moment, the w-momentum equation does not take the accumulating cloud water content into account, and I'm not really convinced I want to mess with this water loading term. If I don't worry about it, I don't think it is a devastating problem re the physics I want to elucidate.

If I ignore the effects of changing the cloud water content, then once cloud is formed, it is simply advected by the flow, probably to leave the domain out the sides. I would guess that I am going to simplify the problem by having L and c_p as constants, ignoring any changes with temperature and ignoring ice physics, clearly.

As I now see the topic, it seems to me that at sufficiently strong "mesoscale" forcing, a steady flow will be produced ... the convection will be plume-like. As the mesoscale forcing is decreased, there will be some critical value below which the convection will begin to "drip" like a faucet; i.e., to become bubble-like. There may be hysteresis, if the mesoscale flow changes during the simulation ... a problem I am not going to consider at this time. The critical value of the mesoscale flow will certainly depend on the initial sounding parameters: the overall CAPE, the strength of the "lid" at the top of the inflow layer, etc. As I have envisioned the physics, then, this is a fairly complex issue. It may be possible, once the solutions have been analyzed some, to develop a non-dimensional parameter that will describe the relative contributions to "dripping" and have some critical value, comparable to the onset of instability with the Richardson number. It remains to be seen precisely how this is going to work, but I am pretty sure we can develop a way to predict the character of the convection in terms of its being "plume-like" or "bubble-like." If possible, this would have some fairly profound implications, I think.

I also am considering a semi-Lagrangian formulation of the advection. It appears that such a form, perhaps with semi-implicit time differencing, might be stable enough to permit relatively coarse time resolution. I need to look into this some, but I don't want to get too hung up over the numerics at this point.

More thoughts ... if the forcing is in some way related to the steadiness of the convection, this is reminiscent of the Binger tornadic storm on 22 May 1981 ... there, the storm clearly was pulling negatively buoyant air up through a restraining inversion. The mesocyclone aloft, in combination with the huge CAPE, was enough to maintain the storm against the restraining inversion and produce a quasisteady updraft. But above the LFC, there still appeared to be convective bubbles superimposed on the flow. Thus, although the flow might have legitimately considered "plume-like" it still had "bubble-like" aspects. I hope that the model will show this sort of behavior for large mesoscale forcing. Perhaps this unsteadiness superimposed on the plume is responsible for the cyclic behavior of supercells ... in rare instances, a very nearly steady state might be reached, but it would require delicate balancing and also might be "convolved" with the sorts of issues pursued in the papers with Harold about the relative strengths of the mesoscyclone and the storm-relative flow in the layer bearing the precip. [Thought added on 18 May 1998: The effect of Kessler-like "condensation oscillations"^[3] resulting from precip being shed from the updraft might also make steadiness an idealization that is nearly impossible to achieve.] In effect, there may be a number of factors that would need delicate balancing to create a very nearly steady storm ... hence, their rarity.

The vorticity equation in the sponge layer:

Within this layer, on the right-hand sides of the momentum equations are the additional terms

For the sake of simplicity, I have not allowed the viscosity coefficients in the horizontal and vertical to be different. The viscosity coefficient, n, does not vary with r, but obviously does so with z (a functionality yet to be specified, except to note that it increases rapidly with height in some way beginning at the bottom of the "sponge" layer).

The derivation of the vorticity equation by cross-differentiation results in terms that look like

They combine to give an extra term on the ride-hand side, denoted by R₁, where

Note that

for any scalar f. Therefore,

Now, using the definition of the pseudo-vorticity

Therefore, with a little rearrangement of terms

Now, it can be observed that, as we have shown earlier,

Therefore, it can be seen that

This term is added to the vorticity equation on the right-hand side, within the "sponge layer" - to damp out the gravity waves. This will also have to be formulated in terms of the streamfunction y'.

Thoughts - begun 09 November

If it turns out best to use a semi-implicit, semi-Lagrangian scheme for the time integration, it seems that a primitive equation formulation might work out best, avoiding all the streamfunction stuff, including getting out of having to solve an elliptic equation every time step. It also simplifies the development of the base state "mesoscale" situation, since there is no longer a need to solve for the streamfunction there, either. The decision remains unmade, though.

In such a formulation, begin with the primitive equations.

Momentum:

The brackets denote that those terms vanish beneath the sponge layer. If I ignore the contributions to pressure change from all the nonhydrostatic effects (as in the earlier derivation) so that the only physical acceleration following the flow is due to buoyancy, this becomes

Mass continuity:

.

First Law of Thermodynamics:

Equation of State:

Definition of Potential Temperature:

Definition of Virtual Temperature:

Note that C* denotes the special character of the cloud water ... unless I change my mind, once it is set to unity, it becomes a conserved variable; i.e., C* is set to zero. The total derivative operator is defined by:

where u_m and w_m are not time dependent.

My question now is enforcement of mass continuity; this has to take place at some point during each time step. How is this done, and when is it done? If we use the momentum equations to advect momentum, then we must enforce mass continuity right away. Right now, I am unclear about how this is done. Since there is no pressure variable in my grossly simplified equations, there is no elliptic equation to solve (derived from the divergence equation).

I am returning to the belief that a vorticity-stream function approach is preferable. I think it might still be possible to use a semi-Lagrangian scheme for integration. If the vorticity equation is written as

and everything involving velocities on the rhs [called F(r,z)] is converted to streamfunctions, and the streamfunction is used to find the velocities for doing the estimated trajectory on the lhs, then it should be possible to do a semi-Lagrangian form.

I'm also wondering if it might not be easiest to use the ARPS model, suitably modified. I need to look at the ARPS manual.

Note added on 13 November:

I don't think the ARPS can be modified to work without a great deal of effort, since it assumes a horizontally homogeneous basic state.

And what about some relatively simple scheme? The choices one must make are rather extensive and it is not entirely clear how to make them. It seems I have reached a critical point in the effort ... a suitable method for simulating the problem is needed and I feel incompetent to choose among the variety of possible schemes.

I just noted that on the rhs of the above we find, via mass continuity, that:

,

where d is the 2-d divergence. It should be noted that

Clearly, following this formulation, mass continuity is taken care of in the step involving the recalculation of the streamfunction. I continue to inch back toward this sort of formulation.

It now seems like it might be easiest to do something simple, just to get started. Perhaps a simple explicit scheme, with semi-Lagrangian advection and a time filter. I think I will try to develop the code as "object-oriented" as possible, so if simple differencing modules don't work well, they can be unplugged and new ones plugged in.

Current status (as of September 2006):

This project was developed as a grant proposal by Dan Weber and myself, submitted to the National Science Foundation, and we're in the third year of our 3-year grant. See here for some hints about what we're doing.

References

[1] Shu, F.H., Gas Dynamics II, University Science Books, Mill Valley, CA., 476 pp.

[2] Pellerin, P., R. Laprise, and I. Zawadzki, 1995: The performance of a semi-Lagrangian transport scheme for the advection-condensation problem. Mon. Wea. Rev., 123, 3318-3330.

[3] Kessler, E., 1969: On the distribution and continuity of water substance in atmospheric circulations. Meteor. Monogr., 10, No. 32, Amer. Meteor. Soc., 84 pp. [see pp. 72 ff.]