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) Mass	conservation or continuity
ρ-ρ0=(P-P0)c-2	Constitutive 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+...

or

c-2d2 P/dt2 = ∇2P+R1

or

c-2d2v/dt2 = c2∇2v +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

∇²P  = - 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.

Summary

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.