Wednesday, October 20, 2010

A ramble on GCMs Navier-Stokes equations.

This post is an outgrowth of a discussion at the Air Vent on GCM's and condensation. It's a kind of response to comment claims that GCN's don't allow for inversion, or volume change during condensation etc. I found myself saying often that it all comes out in the solution of the Navier-Stokes equations - you don't have to do something special.

It's often said that the N-S equations are at the heart of GCM's, but if you look at a set of equations actually solved, as in this set from CAM3, it isn't obvious. And it's made less obvious by the funny coordinates they use.

But it's true, and here I'll try to set out the reasons.

The Navier Stokes equations

These go back to Navier in 1822, although his physics was eccentric, based on a stretched version of recent atomic theory. GG Stokes in about 1845 gave a definitive derivation and solved most of the special cases that could be done theoretically. In the mean time, Poisson and Saint-Venant also made big contributions. In fact, what I'll be talking about mainly goes back to Euler.

I'll write out a version for slightly compressible fluids. An abstract form is:
ρDv/Dt = -∇.τ Momentum equation
Dρ/Dt = -∇.(ρv) Massconservation or continuity
ρ-ρ0=(P-P0)c-2Constitutive Relation
Update - I've added a P0 to that last equation (I'm used to thinking of ambient pressure as 0, but of course in the atmosphere that won't do).  At Lucia's suggestion I've used the more orthodox term constitutive relation (not eqn of state). I'll also mention that ρ is density, v velocity, P pressure and t time.
D/Dt represents a total or advective derivative - in a frame moving with the fluid. ∇ is a vector derivative operator (∂/∂x,∂/∂y,∂/∂z), and τ is a stress tensor. Don't worry about what that means if it's unfamiliar - it's about to disappear.

I've chosen this "abstract" style to show the structure. The momentum equation is just Newton's second Law -   ma = F . Continuity just says that if mass is moved by advection (being carried with the flow), the density will change. And this constitutive relation just says that to change the density from a standard (local hydrostatic) value, you have to apply a pressure P, and c is the speed of sound. You could have other constitutive relations for different materials.

The stress tensor could be quite complex, but for now, let's just say there's a pressure P, gravity g, and a viscous stress (coef μ). Atmospheric flows are always turbulent, and turbulence models can work by modifying the viscosity. But I mainly want to talk about the pressure-velocity interaction, so I'll just keep a viscous term as an example of the other things that could come in.

OK, so now we have:

ρDv/Dt = -∇P + ρg + ∇.(μ∇v) Momentum equation
DP/Dt = -ρc2∇.(v) Mass conservation or continuity

Time Scales and numerical solution

Suppose you wanted to solve these numerically. You'll have the space divided up into regions, with the field variables (P, v etc) represented on nodes - maybe the corner points. You'll have linear expressions for the space derivatives in terms of nodal values.

So a naive approach is to say - suppose you calculate the right hand sides at a point in time. just calculate the rates of change of v and P, and assume that holds for a timestep. Update and proceed. That's called an explicit method (forward Euler).

But it won't work, except for very small time steps. The reason is that, while you've got starting information within your little regions, what's going on in one region can affect others during that step. You're limited by the fastest process that can propagate physical effects. And that is the speed of sound, which is the way pressure is propagated.

The reason is that buried in the N-S equations is the acoustic wave equation:
ρdv/dt = ∇P dP/dt = -ρc2∇.v+...
c-2d2 P/dt2 = ∇2P+R1
c-2d2v/dt2 = c22v +R2
I've replaced D by d, because only fast processes are retained, and the difference between D and d is an advective term, which is also slow relative to sound speed. R1 and R2 denote the slow terms.

By small time steps, I mean small enough so that a sound wave will only move part way across one of the small regions. In GCM's, the vertical divisions can be just a few hundred metres. That's a second or so of sound time. That would be ridiculously short for GCM modelling, where steps of about half an hour are common. Incidentally, these are constrained by the speed of other "weather-like" processes.

Coupling and semi-implicit solution

The problem is that with the explicit method, we have to resolve the acoustic oscillations, which are not of interest.

Looking at the pressure wave equation, we can say that, since c is large, only acoustic solutions will have a significant left side. For the solutions we actually want, that can be dropped. Then we get

2P  = - R1

wher again R1 is a collection of slow terms. This is what is called the Pressure Poisson equation. Once solved, the pressure can be used to get the velocity by going back to the momentum equation.

The essential thing here is taking the pressure and velocity equations together - coupling. Everything else can be held constant during the time step (though you can also do better with staggered time steps etc - there's lots of extra art not relevant here).

So now the time step can be extended until the next fastest process will mean that the computational cell cannot be regarded as (sufficiently) isolated during the timestep. Since that process is generally, unlike sound, one of the processes you do want to resolve, that's basically it.


So that's why the Navier-Stokes equations are the core. They can be seen as an engine that matches pressure and velocity, using sound implicitly, but not requiring sound-limited timesteps. Once you have that sorted out in a time-step, all the other physical processes can be updated in good time.

Heretical afterthought - for practitioners

Actually, implicit solution is not as rosy as it is often painted. You can take long time-steps, but you have to solve a big matrix of equations to do it, instead of just calculating new values at each node. You take advantage of the fact that each node only affects a few neighbors, which gives a sparse matrix. So for big systems you use an iterative solver.

The gloomy fact is that the iterative solver, in effect, is an accelerated version of the original explicit solution. And using something like conjugate gradients, the acceleration is significant, but not super. Preconditioning helps, as does a process called multigrid. But the fact remains that, as you extend the timestep, the matrix becomes more ill-conditioned (an analyst's curse word), and you end up with a number of iterations rather related to the number of explicit steps you would have had to perform.

I might say more about preconditioning in a future post.


  1. Thanks for participating at The Air Vent, WUWT, and other skeptical sites.

    Too many people approach climate blogs with an attitude of "let's win the debate" rather than "let's figure out the science".

    Your posts (and the responses to them) move the blogs more towards science discussions and away from echo chambers endlessly repeating conspiracy theories.

    Charlie A.

  2. Nick, NSE equations have nearly nothing to do with GSMs. This is not only the "funny" terrain-following coordinates, there is much more.

    The "dynamical core" of GSMs is built on _ideal_ (Eulerian) fluid approximations developed in Meteorology 40 years ago when no reasonable commuters existed. As far I understand, the atmosphere is eventually presented as several 2D layers ("levels") of ideal fluid. The core equations do not have the native molecular viscosity of NSE (there is not single reference to "viscosity"), and therefore the initial motion would be sustained forever, which, combined with inevitable errors of discretization (in physical or spectral domains) would lead to divergent solutions. They avoid this blowing up ("divergence") by re-introducing an artificial dissipation between atmospheric layers, so the models would run "stable" (sometimes), see "initial divergence damping".

    Another drawback of this approach is that the vertical component of dynamics does not naturally arise in this scheme, and must be modeled parametrically (see "physics of deep convection"). All this require various tricks/adjustments to maintain various laws of conservation across layers.

    The other caveat of all this is that the initial inviscous equations have different, reduced set of solutions due to truncation of the term with highest spatial derivative. Re-introduction of effective viscosity at different level does not restore the broken original topology of solutions, so that the long-term behavior of GSM has nothing to do with NSE, and therefore cannot be trusted, no matter how good your solver is or could ever be.

  3. Charlie,
    Thanks. Actually, I must say it is great to have an audience for talking about science :)

  4. Al,
    I'm unfortunately off to an early dentist appointment - hope to reply in suitable detail in a few hours.

  5. Don't you need thermodynamics equations to properly describe condensation? Also, even in the context of the hydrodynamic equations, you need additional terms. see for example this.

    I also don't think these harken back to Stoke's day, do they?

    Thanks for an interesting post, and I echo Charles A's comment.

  6. Al,
    Yes, there are obvious adaptions of N-S. Horizontal momentum diffusion (viscosity) is ignored, which makes them rather Euler-like. And as you say, special things are done for vertical plumes.

    But the pressure-velocity coupling is still central. In CAM3, this is described in 3.1.4. I pointed to one set of equations there, but a better complete set is 3.33-3.42. Then as they say, the N-S (or Euler) equations are 3.33,3.34 and 3.37 (they include the energy equation in T as a core equation).

    The semi-implicit solution in CAM3 is interesting - I learnt a lot by re-reading. They do a lot of manilulations on the variables, set up time-splitting, and transform to a spectral space using spherical harmonics (3.12 etc) as a basis. Still, it's a Poisson solver, and is the core step. I hadn't expected them to couple in temperature - that means that they have identified temp transport (probably vertical) as the next fastest process, and are prepared to wear the cost of a larger matrix inversion to get a longet timestep.

    The special treatment of plumes is presumably also time-step related. When the atmosphere is stable processes where layers affect one another are relatively slow. So even though the layers are rather thin, a fairly long timestep can be managed. But instability (storms, hurricanes etc) would destroy this. So if it can be taken out of the process by special treatment, the time step can be extended, though I suspect not too far. I'm trying to figure out more about what CAM3 is doing in sec 4, which is relevant to Jeff's inquiry. I;d be glad of any observations that you have.

    You're right about the removal of viscosity and the introduction of dissipation. But that really means they have an anisotropic viscosity tensor, mainly diagonal, but with two zeroes on the diagonal.

    Thanks for your comment - I think you are more familiar with practice than I am. I'd be grateful for any more interpretations etc that you have,

  7. Carrick,
    Yes, indeed you do. N-S equations are often considered now to have a coupled energy equation, which covers transport of heat, and picks up viscous work etc. Yes, that's wasn't there in 1854.

    In my discussion above I had thought of this as one of the slow processes that could be decoupled. But I see CAM3 thinks it is worth including in the semi-implicit step.

    Euler's equations, which I guess I'm really talking about (as Al says, I think) were characterised by Feynman as what you need to describe the flow of dry water. It's the essence of fluidity (fast interchange of pressure velocity). N-S gives you wet water. Most other things you can just add on.

    Thanks for the kind comments. Despite my occasional cantankerousness, I do appreciate the core of scientific interest around.

  8. Actually, it just dawned on me why they have to couple the energy equation. Radiation is even faster than sound, and has to be treated implicitly. It can't really then be separated from other slower heat transfer processes. D'oh.

  9. Re Al #2, if the GCM's are based on 2D layers and they are "artificially" held stable by adjustments to maintain mass balance etc, then the issue of mass balance re the missing H2O in the rain-out problem being discussed at TAV is likely to be even more significant I would have thought.

  10. Nick--
    Is this discussion for cases with condensation? Or without?

    I'm trying to figure out which questions to ask. Are these equations for the mixture of (gas + liquid water + ice) with gas being (air + water vapor)? Or are they for the gas (i.e. air+water vapor) only?

    If it's for the gas only, why doesn't your first continuity equation include a source term to account for water vapor evaporating?

    Regardless of whether these are for gas only or mixture, wouldn't the condensation issue relate to the questions like, "what is the mathematical formulation for "c" in your equation of state:
    (I switched to primes...)

    That equation is only the equation of state is based on of sound being "c2 = ∂p/∂&ρo|s

    For an ideal gas , that's going to be root(gamma R T), but for a mixture of gas and a condensing fluid, additional terms likely (as in almost certainly) appear.

    Anyway, if you can tell me whether these equations are for the gas (i.e. water vapor + air) or mixture (gas + liquid water + ice) that will help me understand what you are trying to say.

    Which is fine. But doesn't that just push Jeff's question back to: And do you know the equation for "c"? That needs to be coded. And if it's a gas/liquid mixture undergoing condensation, do GCM's contain the correct formulation for the speed of sound?

  11. "...account for water vapor evaporating?"
    Water evaporating or condensing.... (I'm sure I'll see other issues later. :) )

  12. Lucia,
    No, it's de-emphasising things like condensation. The CAM3 documentation that I cited has a long string of equations, and I'm trying to focus on the ones that require close coupling, so they go into the semi-implicit solver. Others are then solved as a catch-up process.

    I've used a linear equation of state, mentioning that there are others. Condensation could be included there. I also restricted at this stage to slightly compressible fluids, which excludes condensation.

    Still, as Al points out, it's a real issue in the models, because they use special models for upward plumes etc, which rather bypasses the N-S solver. And of course it's the Issue Jeff is trying to get at.

    The equation of state isn't just gases - but it is linearised. That's enough to bring out the acoustic wave equation, which is the reason for tight coupling (semi-implicit).

    Other issues - yes, there'll be lots. I'm not trying to do a whole GCM in one post. The motivation is really to have something to say when people ask for all sorts of special cases to be included, when the N-S equation solver handles that.

    I hadn't seen an allowance for condensation coming into the speed of sound particularly. I'm envisaging a more complex eqn of state, of which the ss is a derivative, but with that complexity, the derivative (ss) isn't constant.

    Here my main reason for introducing ss is to show which terms divided by it become small, and drop out to provide the Pressire Poisson Equation.

  13. Nick--
    Ok. Because I the first two sentences of the post gave me the impression this had something to do with the question of whether or not GCM's address condensation properly. Obviously, this discussion falls pretty short of the mark if we want to talk about how effects of condensation would be included or incorporated.

    >> The equation of state isn't just gases - but it is linearised.

    Even your 'equation of state' isn't a linearized form of the real equation of state, is it? The ideal gas law is a real, honest to goodness equation of state and that differs from the one you give. If you wanted to derive the expression for "c" in your 'equation of state', you would need to identify the *real* *honest to goodness* equation of state, add some thermo and also reveal the assumptions required to use your constitute relation connecting pressure to density. (Sorry. I can't bear to call that your third equation "equation of state" without including scare quotes.)

    >>when the N-S equation solver handles that.
    Ok... but I don't get where you are going with this.

    The N-S solver can only handle things that can be handled by the NS. And the set of equations you provide can work in cases where you certain assumptions hold and you can replace the *real* equation of state with the one you suggest-- which only holds under special circumstances. (Maybe they hold under all the cases you care about? Clearly not shock waves though-- but you probably aren't concerned about those.)

    >>I hadn't seen an allowance for condensation coming into the speed of sound particularly. I'm envisaging a more complex eqn of state, of which the ss is a derivative, but with that complexity, the derivative (ss) isn't constant.

    Condensation and evaporation *caused by pressure changes due to flow* should affect the speed of sound-- no doubt. So, it ought to matter when the relative humidity is high enough to anticipate condensation when temperatures drop.

    I don't know how much condensation would change the speed of sound, but it should change it some.

    Another pesky issue: I don't know the upper range of mass loading of water drops in air, but their presence may also affect the speed of sound. (Maybe not much though-- I'd guess an upper bound of 10% if we can assume drops are at the same temperature as surrounding air. More if we can't. But I'm guessing mass loading of water. If you know, I could swag that better. )

  14. Nick--
    Ok. I thought you were dealing with condensation. I'm not sure where this is going.

    Condensation induced by temperature variations arising from flow should affect the speed of sound. If you recall the derivation for idea gases, and think of a multi-phase media, you'll see it has to do so. The speed of sound will vary spatially for this reason.

    I think if you want to use this going forward, you need to elaborate how you got what you call an "equation of state". I'd call it a constitutive relation because the *real* equation of state for gas is the idea gas law, which differs from the one you show. You'll need to explain the assumptions restricting your substitution so people can understand when and why it holds in the cases you think "matter".

    Oh.. btw. Droplets in air can change the speed of sound. I doubt it's much. Do you happen to know the upper bound on mass loading of something like raindrops in air? Or water drops in clouds?

  15. Lucia,
    "I'm not sure where this is going."
    Me either. I did call it a ramble :)

    But your query has helped me to clarify that in my mind. I really want here to lay out some computational structure of the N-S in GCM's. That's why the emphasis on the wave eqn and the transition to the PPE - semi-implicitness.

    There is a large section of code where they do fnncy numerical stuff on this subset of variables (pressure, velocity, temp). Spectral transformation etc. And I suspect it's where they chew up the computer time.

    Then there's the loosely attached code where they catch up with all of what I've called "slow" processes. I'm not so interested in the detail here, at this stage, because you can study each in isolation. And I think condensation is in that category.

    So my next post will probably try to describe the actual CAM3 semi-implicit part. Then maybe another on the parts treated specially - plumes etc. Although I'm reasonably familiar with CFD, I'm just learning about all this.

    So on your specific q's:

    Yes, constitutive relation might have been a better term. Normally for the dynamics of atmospherics, incompressible flow is a useful approximation. It's not that you don't get density variations - obviously with altitude, but also with humidity etc - but rather that the dynamic pressure variation associated with flow velocities is small, and so it's dependence on density can be linearized.

    If heat transfer is assumed slow, then the right invariant expression for adiabatic fluctuations is
    ρ^-γ P=const. ρ=density, γ=cp/cv, P=pressure
    So with this slightly compressible gas, dynamic and acoustic fluctuations satisfy
    P-P0 = (γP0/ρ0)(ρ-ρ0) = c^2 (ρ-ρ0)
    which if you then set equilibrium pressure P0 to be zero (which I forgot to spell out) becomes
    P = c^2 (ρ-ρ0)
    So that's the sense in which it's the equation of state. But I'm happy to call it the constitutive relation - it sounds better.

    I'm not too concerned about what the speed of sound actually is, because the point of writing this out is to eliminate acoustics. We need to know that they are there, and to identify the wave equation sitting in the N-S equations, just so they can be eliminated.

    In saying "the N-S eqn solver handles this" I'm touching on two points
    1. There is a subset of the GCM system which is N-S. It may have source terms in the continuity equation, say, and lots of extra things in the momentum equation. But it does contain the relation between pressure and velocity, and you don't have to put in extra relations to cover that.
    2. The N-S solver is a distinct part of the code (semi-implicit (with energy eqn)). Some effects go into it, some not, and it's necessary to track which is which.

    The condensation effect on ss would then go into the RHS of the Pressure PE.

    The issue of condensed water mass/weight was discussed in Jeff's threads eg here. Much discussed was the question of whether condensation increased or decreased pressure.

  16. Nick
    >> Normally for the dynamics of atmospherics, incompressible flow is a useful approximation
    Of course to leading order, flow is incompressible. The assumption will be useful in many cases. (Since you aren't setting divV=0, I assume by incompressible, you probably mean barotropic or baroclinic, not literally incompressible? )

    >>So with this slightly compressible gas, dynamic and acoustic fluctuations satisfy

    This is where you are inserting an assumption-- which can be important in special cases. So, if your intention is to build on this to explain how things are generalized, you need to reveal t this or risk giving people the impression what you are doing applies in the later special case.

    If there is heat transport or phase change (evaporation or condensation) then, variations in temperature, pressure and density will not vary as we see for sound waves in an idea gas. That's where you need a different relation for ∂P/∂ρ at constant entropy. Or, if you prefer, this equation you substitute into your system doesn't necessarily apply:
    P-P0 = (γP0/ρ0)(ρ-ρ0).

    (Temperature will also not vary according to the corresponding equation for temperature variations. You need to do the thermo! )

    Now, I might not be yammering about this if I meant that equation might only apply *approximately*. But the problem is *it might not even be close*.

    Look at a P-v diagram (i.e P-&1/ρ) for steam/water and find the lines at constant "s". The acoustic relation you are substituting right from the start is fine at temperatures high or pressures low enough to get you away from the "hill". Condensation isn't one of those cases. It's not even close. (I can't find a good P-v diagram on line.)

    >> But it does contain the relation between pressure and velocity, and you don't have to put in extra relations to cover that.
    But the Navier Stokes equations *don't* contain that relation betwen pressure and velocity in many flows. To get that, you need to provide a *real* equation of state, and sometimes use the first and 2nd law of thermo. You are using constitutive equation that has made quite a few assumptions.

    So, if you start where you did and don't reveal your assumptions, you can't "generalize" to special cases that might violate those assumptions. So, I think to meet your goal of going forward, you need to reveal your assumption and explain why they apply to some case that interests you -- or better all cases that could possibly affect the accuracy of a GMC.

    The assumptions made to come up with the equations coded in a GCM must be routinely document in development of climate models, right? We do that in engineering-- routinely. (Otherwise, people would screw up and think they can the stuff you wrote down when solving heat transfer problems. BIG MISTAKE!)

    So, even if you don't normally do that sort of thing yourself, the development must be *somewhere*. (Mind you, I haven't seen it... but then, I don't develop or run GCM codes.)

    (I won't double post if this seems not to show up because of length this time. It seems they do show up.)

  17. Lucia,
    I'm allowing just enough compressibility to get sound - basically to show how the PPE comes in to then eliminate acoustics.

    A general theme is the handling of disparate time scales. A GCM is faced with picosecind IR upwards, and trying to cover centuries. The E/M timescale goes out with ray theory (plus scattering models etc). But the sound is real, and requires PPE semi-implicitness. That's heavy computing. That gets us to about the 20-min scale; then the passage to centuries is just big computers.

    Actually handling pressure needs further investigation, because a lot of it is done with a hydrostatic assumption. But I do think it's implicitly there in the Poisson equation.

    I think condensation may have a role in the inclusion of temperature in the semi-implicit phase. It must be because of something fast, and my firsat thought was E/M. But with rays, that's not really PDE at all. So then I thought that pressure waves might interact with vapor, in effect propagating heat at high speed. That's probably equivalent to your idea of an effect on sound speed.

    I have the Handbook of Heat Transfer Fundamentals by Rosenow, , with plenty of P-v diagrams.

    I'm hoping no-one would want to cite a blog ramble for heat transfer problems :). But I'm also hoping that all will become clearer (to me, anyway) getting more specifically into the CAM3 document.

  18. Nick---
    Of course the sound is real. Yes. there are numerical issues. Few will cite a blog ramble-- especially if it's a ramble. (Though citing a blog doesn't technically violate any offcial rules of academia. It may violate IPCC rules for what to include in their reports-- but that's a bit different.)

    I guess I still don't know what you are trying to show others. Yes. You can sometimes separate stuff.....

  19. Nick,
    NSE are only the momentum equations. You need separate continuity, energy, and eq of state, as well as some thermo to get there from here. In fact, a time varying 3-D high Reynolds number flow is not solvable for useful time scales in the real world even with the best computer. When you include limited initial and boundary data and inadequate grid you have what is known as a joke.

  20. Leonard,
    I included momentum and continuity, and an equation of state, which Lucia said would be more conventionally called a constitutive equation - anyway, it relates pressure and density.

    Energy wasn't in the original N-S equations, and I wanted to separate them here to concentrate on the pressure-velocity time scale (acoustic wave equation). However, as Al said, the GCM's do couple the energy equation in the semi-implicit solution, so it probably should be considered part of the system.