`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 of10 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

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

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_{1} and a dewpoint lapse rate given by
g_{d1}, to a height of
h_{1}. 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

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

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_{1} 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.

- : 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.
- : 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).
- : when w' is decreasing with z, this is "vertical convergence" that acts to increase the positive vorticity associated with w increasing with r.
- : 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

negativelycorrelated. In effect, when this correlation is negative, relatively large (negative) values of r_{o}u' are pushed upward, or small values are pushed downward, both leading to positive pseudo-vorticity3. : the density term, by virtue of being squared, is positive-definite. Now this term is positive when u' and w are

positivelycorrelated, 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_{o}u' is becoming less negative as one goes up.

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

This will still need some work, to be put in terms of the streamfunctions, y

** **

**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_{1}, 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.

** **

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.

[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.]