We built a numerical solution to our full model using a finite difference approach. We will started by implementing a finite difference scheme for a simple advection equation in one dimension with constant wavespeed and no forcing. In one dimension we use the domain $[0,1]$ which represents a column in the atmosphere. The velocity $u:[0,1]\times[0,T]\rightarrow\mathbb{R}$ must satisfy

$\frac{\partial u}{\partial t} + v \frac{\partial u}{\partial x} = 0, $

where v is a constant wavespeed. We will implement zero Dirichlet boundary conditions, that is, $u(0,t) = u(1,t) = 0$ for all $t \in [0,T]$. We know that this equation has solution $u(x,t)=g(x-vt)$ for any $g \in C^1([0,1])$ that satisfies the boundary conditions. Knowing a solution provided a useful starting point for our numerical solution as we were able to test the correctness/accuracy of the method.

For the finite difference approximation we introduce a grid on our domain with N points with equal spacing $\Delta x$ and a timestep $\Delta t$. The approximation of the solution at a grid point $i\in\{0,1,2,...,N\}$ at timestep $n\in\{0,1,2,...,T/\Delta t\}$ is denoted $u_i^n \approx u(i\Delta x, n\Delta t)$. We will use an upwind scheme for advection, this provides the finite difference approximation of the spatial derivative $\frac{\partial u}{\partial x}$ as

\begin{align} \frac{u_{i}-u_{i-1}}{\Delta x} &\mbox{ for } v \geq 0,\\ \frac{u_{i+1}-u_{i}}{\Delta x} &\mbox{ for } v < 0, \end{align}

and we use an explicit Euler method in time. We note that for stability of the solution we require a condition on the timestep,

\[ \Delta t < \frac{\Delta x}{v}. \]

We test this method with initial velocity $u(x,0) = \sin(2\pi x)$ and periodic boundary conditions. The numerical solution displays numerical diffusion. This is an inherent flaw in the approximation method, however since we will eventually be introducing diffusion into the equation which will be larger than the numerical diffusion, hence it's effect on the solution will be insignificant. We continued building up the model by adding in diffusion and forcing terms and non-constant wavespeed to give Burger's equation. Then we moved into two dimensions and incorporated incompressibility giving us a simulation of Navier-Stokes in 2D. Then it was relatively simple to include the density and temperature equations for the full model. The full stability condition for our model, coming from the velocity and temperature equations, is

\[ \Delta t \leq \min\left\{\frac{\Delta x^2}{2(4\nu + \sup_{i} |u^n_i| \Delta x)}, \frac{\Delta x^2}{2(4\mu + \sup_{i} |u^n_i| \Delta x)} \right\}. \]

Below are animations of the simulations. For the first one, we set the initial conditions as follows, the cloud density q_c is set to zero everywhere and is measured in hPa kg J^-1 or equivalently kg m^-2 / 100. The water vapour density q_v is set to an exponential of temperature, specifically 90% of the saturation density. Physically this is a large amount of water vapour, we use this initial vapour content to ensure some cloud behaviour; we expect a fairly large cloud formation.

In the next animation we start with a pocket of saturated air near the ground. The initial motion and formation is somewhat different from the first. The concentrated pocket of saturated air quickly condenses to form a dense cloud, then this cloud is spread around over time. However in the first simulation, the water hasn't reached saturation density at any point initially, therefore the cloud slowly condenses over time.

Unexpectedly the buoyancy force causes the vertical velocity to become negative and forces the formation into the ground. We found that this forcing was very sensitive to the ambient air density parameter $\rho$ and we could not find a stable value. The below animations use a forcing term on the horizontal velocity and no forcing on the vertical velocity to illustrate the model further.