Documentation of the simple dynamical core part 2: Continuous equations

September 12, 2022

The system of continuous partial differential equations of motion for the atmosphere (relating the variables $\mathbf u = (u,v,w), \rho, T, p$) in its general form consists of:

Three equations for the conservation of momentum, usually written as a vector eqution. In the most general case, this is the Navier-Stokes equations on a rotating planet (neglecgting viscosity terms as these are unimportant for large-scale flow) $$ \frac{\partial \mathbf u}{\partial t} + \mathbf u\cdot \nabla\mathbf u+ 2\mathbf\Omega \times\mathbf u=-\nabla\left (\Phi -\frac{1}{2}(\mathbf \Omega \times \mathbf u)^2\right )-\frac{1}{\rho}\nabla p.$$
The continuity eqution for the conservation of mass $$ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf u) = 0.$$
The thermodynamic equation for conservation of energy $$ \frac{\partial\theta}{\partial t} + \mathbf u\cdot\nabla(T) = \frac{\theta}{T}\frac{Q}{c_p}$$ where $\theta = T(\frac{p_0}{p})^{(R/c_p)}$. In addition, a relation between pressure, density and temperature is needed, like the ideal gas law.

This system of equations is complicated to solve numerically due to the presence of fast sound waves propagating in three dimensions. Typically, a series of approximations are made to get the set of equations called the hydrostatic primitive equations in the traditional approximation. I chose this set of equations to implement in my model since this is the typical set of equations used in most global climate models (including CESM which I used for my PhD). The approximations are:

The hydrostatic approximation:
- No acceleration in the vertical, gravitational term balances vertical pressure gradient. This turns the vertical equation of motion into the hydrostatic equation.
The shallow atmosphere approximation:
- $r = a+z$ where $a$ is the radius of the Earth, and $a»z$. This means that the vertical scale of motion is much smaller than the horizontal scale. This means that we can replace $r$ with $a$ except when differentiating with respect to $r$ (when $z$ is commonly used instead).
The traditional approximation, which needs to be made together with the shallow atmosphere approximation:
- Metric and coriolis terms involving the vertical velocity are neglected (metric terms come from the expression of the total derivative of a vector in spherical geometry. I will point out the remaining ones later).

After these approximations, the equations of motion are:

Conservation of momentum: $$\frac{\partial \mathbf u}{\partial t} + \mathbf u\cdot \nabla\mathbf u+ f\mathbf k \times\mathbf u=-\frac{1}{\rho}\nabla p$$ where $\mathbf u=(u,v)$ now is the horizontal velocity vector. The vertical equation of motion has turned into the hydrostatic equation $$\frac{\partial p}{\partial z}=-\rho g$$
Conservation of mass is unchanged when written in terms of the $\nabla$-operator and the velocity vector, though some terms have changed under the hood.
Conservation of energy is also unchanged in the same sense.

In all of the above the implicit coordinate system is a form of tue “geoscientific sherical coordinate system” $(\lambda,\phi,z)$ with $\lambda\in[-\pi,\pi)$ representing longitude, $\phi\in[-\pi/2,\pi/2]$ representing latitude, and $z\in[0,z_T)$ reprenting height above sea level in meters.

Since the system of equations lacks an equation for $w$, numrical solution of it is complicated. It can be made more tractable by changing from height to pressure as the vertical coordinate: Instead if $z$ we take $p\in[p_T,p_s]$ where $p_T, p_s$ are the pressure at the model top and the Earth’s surface, respectively. The variables in this system are $\mathbf u = (u,v),\omega, \rho, T, z$. This change implies replacing the pressure gradient force with a term representing the gradient of the geopotential of a pressure layer. The continuity equation takes the simple form

$$\nabla\cdot \mathbf u + \frac{\partial \omega}{\partial p}=0.$$

At this point, we need to talk about upper and lower boundary conditions. At the top of the model ($p=p_T=const.$), we have $\omega=\frac{dp}{dt}=0$, assuming zero mass transfer across the model top. At the surface ($p=p_s$) we have

$$\omega=\omega_s=\frac{\partial p_s}{\partial t} + \mathbf u_s \cdot \nabla p_s.$$

The continuity equation and upper boundary condition can be used to calculate the vertical pressure velocity $\omega$ by integrating from the top of the model, down to a specific level: $$\omega=-\int_{p_T}^{p}\nabla\cdot\mathbf u dp.$$

Extending this integration to the surface we have another expression for $\omega_s$: $$\omega_s=-\int_{p_T}^{p_s}\nabla\cdot\mathbf udp.$$

Substituting this expression for $\omega_s$ into the lower boundary condition and rearanging gives us the surface pressure tendency equation: $$\frac{\partial p_s}{\partial t}=-\mathbf u_s \cdot\nabla p_s-\int_{p_T}^{p_s}\nabla\cdot\mathbf udp.$$

Using the ideal gas law to emilimate $\rho$ from the hydrostatic equation and integrating from the surface ($z=z_s$) yields the following expression for the geopotential height of a given layer:

$$gz = gz_s+R\int_{p}^{p_s} T d(\ln p).$$

So to summarize, to implement the atmospheric model I used the primitive equations in the form described above:

Equation of motion for zonal wind $u$ (prognostic, i.e. involving time derivatives): $$\frac{\partial u}{\partial t} + \mathbf u\cdot \nabla u + \omega\frac{\partial u}{\partial p}+ fv + \frac{uv}{a}\tan\phi=-g\nabla z + F_\lambda$$
Equation of motion for meridional wind $v$ (prognostic): $$\frac{\partial v}{\partial t} + \mathbf u\cdot \nabla v + \omega\frac{\partial v}{\partial p} - fu + \frac{u^2}{a}\tan\phi=-g\nabla z + F_\lambda$$
Thermodynamic equation (prognostic): $$\frac{\partial T}{\partial t} + \mathbf u\cdot\nabla T + \frac{\omega T}{\theta}\frac{\partial \theta}{\partial p} = \frac{Q}{c_p}$$
Surface pressure tendency (prognostic): $$\frac{\partial p_s}{\partial t}=-\mathbf u_s \cdot\nabla p_s-\int_{p_T}^{p_s}\nabla\cdot\mathbf udp$$
Vertical pressure velocity equation (diagnostic, i.e. involves no time derivatives): $$\omega=-\int_{p_T}^{p}\nabla\cdot\mathbf u dp$$
Vertical pressure velocity equation at the surface (diagnostic): $$\omega_s=-\int_{p_T}^{p_s}\nabla\cdot\mathbf udp$$
Equation for the geopotential height (diagnostic): $$gz = gz_s+R\int_{p}^{p_s} T d(\ln p)$$

In the equation of motion and the thermodynamic equation the terms $F_{\lambda,\phi}$ and $Q$ represent all other sources of force (typically friction) or heating (absorption of sunlight, moist processes), respectively. The terms involving $\tan\phi$ are the remaining metric terms.

The gradient of a scalar in the spherical coordinate system assuming a shallow atmosphere is: $$\nabla\psi=\frac{\hat i}{a\cos\phi}\frac{\partial\psi}{\partial\lambda}+\hat j \frac{1}{a}\frac{\partial\psi}{\partial\phi}.$$ The divergence of a vector is: $$\nabla\cdot\mathbf A=\frac{1}{a\cos\phi}\left[\frac{\partial A_\lambda}{\partial\lambda}+\frac{\partial (A_\phi\cos\phi)}{\partial\phi}\right].$$