The Semi-Lagrangian method is a numerical method for solving advection equations. For the model we have proposed, this will be important to solve the equations for physical quantities, such as for \( T, q_{c}, q_{v} \) - the temperature, cloud density and vapour density variables respectively. We are interested in numerically solving advection equations of the following form:

Let \(d \in \{ 1,2 \}, T \in \mathbb{R}_{>0} \), \(\Omega \subset \mathbb{R}^{d} \). Given \(v: \Omega \times (0,T) \to \mathbb{R}^{d} \), \( f: \Omega \times (0,T) \to \mathbb{R} \), we look for a solution \( \phi : \Omega \times (0,T) \to \mathbb{R} \) that satisfies:
$$\frac{\partial \phi}{\partial t}(x,t) + v(x,t) \cdot \nabla \phi (x,t) = f(x,t) \quad \text{on} \,\,\,\, \Omega \times (0,T).$$

Here \(v\) is a velocity field that advects the variable \( \phi \) and \( f \) is a forcing term. First we consider solving an advection equation without the forcing term:
$$\frac{\partial \phi}{\partial t}(x,t) + v(x,t) \cdot \nabla \phi (x,t) = 0 \quad \text{on} \,\,\,\, \Omega \times (0,T).$$

We define a characteristic starting at \( x_{0} \) to be the solution of the following ordinary differential equation:
$$\frac{\text{d} \mathbf{X}}{\text{d} t} (x_{0},t) = v( \mathbf{X}(x_{0},t), t) \quad \text{and} \quad \mathbf{X}(x_{0},0) = x_{0}. $$
Along such a characteristic, the quantity $\phi$ is conserved:
$$ \frac{\text{d}}{\text{d} t} \left( \phi ( \mathbf{X}(x_{0},t),t ) \right) = 0. $$

This gives a method for solving numerically solving the advection equation, by tracing the trajectory of a point back to time \( t_{i} \) to estimate the value of \( \phi \) at time \( i+1 \). We now consider the numerical implementation of this method.

We first make the assumption that \( \Omega \) is a rectangular domain. Let \( J \in \mathbb{N} \) and \( dz > 0 \) satisfy \( J dz = L \). We consider the case \( d = 1 \) as it is more concise and the case \( d = 2 \) is also similar. We discretise our domain \( \Omega \) by:
$$ z^{j} = j \, dz, \quad \text{where} \quad j \in \{ 0,1,\ldots,J \} ,$$
we also discretise time using \( N \in \mathbb{N}, \, dt > 0 \) satisfying \( N \, dt = T \):
$$t_{i} = i \, dt, \quad \text{where} \quad i \in \{ 0,1,\ldots, N \}.$$

In one-dimension, we denote the velocity by \( w \). Given a grid point and a time \( (z^{j},t_{i+1}) \), we make the assumption that there exists a characteristic \( \mathbf{X}(z_{0},\cdot) \) satisfying \( \mathbf{X}(z_{0},t_{i+1}) = z^{j} \). This characteristic is denoted by the curve on the graph below, which we then approximate by a straight line.

We approximate the position of the characteristic at time \(i \) given by \( \mathbf{X}(z_{0},t_{i}) \), belonging to the charateristics using:

$$ \tilde{z}^{j}_{i} = z^{j} - w(z^{j},t_{i+1}) \, dt. $$

Here, we approximate the characteristic using the velocity at time \( i+1 \). As the point \( \tilde{z}^{j}_{i} \) is not necessarily a grid point, in general we can not evaluate \( \phi(\tilde{z}^{j}_{i}) \) directly in order to calculate \( \phi(z^{j},t_{i+1}) \). Instead, we linearly interpolate the values of \( \phi \) between the nearest grid points at time \( i \), which we denote by \( z^{j_{i}}, z^{j_{i}-1} \). We also incorporate the forcing terms in the same way:

\[ \phi(z^{j},t_{i+1}) = \left(1 - \alpha^{j}_{i} \right) \, \left( \phi(z^{j_{i}},t_{i}) + f(z^{j_{i}},t_{i}) \, \Delta t \right) + \, \alpha^{j}_{i} \, \left( \phi(z^{j_{i} - 1},t_{i}) + f(z^{j_{i} - 1},t_{i}) \, \Delta t \right). \]

In order to implement the Semi-Lagrangian method, we require that for each gridpoint \( z^{j} \), the point \( \tilde{z}^{j}_{i} \) lies in our domain, \( \Omega = [0,L]^{d} \). This can be enforced using the following timestep restriction:
\[ \Delta t \leq \inf_{j \in \mathcal{J}_{i}} \left \{ \frac{z^{j}}{w^{j}_{i+1}} \chi_{ \{ w^{j}_{i+1} > 0 \} } + \frac{z^{j} - L}{w^{j}_{i+1}} \chi_{ \{w^{j}_{i+1} < 0 \} } \right\}, \]

where \[ \mathcal{J}_{i} = \left \{ \,\, 0 \, \leq \, j \, \leq J \,\, \bigg| \,\, w(z^{j},t_{i+1}) \, \neq \, 0 \,\, \right \}. \]

One of our aims for the Semi-Lagrangian method was to simulate our model using the two following temperature equations and comparing the behaviour of the two results, using a constant velocity profile in two dimensions:

$$ \frac{\partial T}{\partial t} \, = \, - (\mathbf{u} \cdot \nabla) T - \Gamma_{d} w + QC_{c}. \label{1} $$
$$ \frac{\partial T}{\partial t} \, = \, - (\mathbf{u} \cdot \nabla) T +\mu\Delta T + QC_{c}. \label{2} $$

Below are the results for the cloud density \( q_{c} \) that we obtained. The result was obtained using the first equation and the result on the right was obtained using the second form of the temperature equation.

In both of the models, the peak of the cloud density term corresponds to a height of at least 3km, which matches observations. However, in order to generate the results for the second equation (corresponding to the image on the right) using the Semi-Lagrangian method, the thermal diffusivity term \( \mu \) was artificially inflated in order to allow saturation to occur. This is the case as initially, for the second temperature equation, there is initially no forcing acting on the temperature to drive the system. Inflating the value of \( \mu \) helps avoid this problem, at the expense of introducing a non-physical parameter.

Initial values of DeltaT, using the first temperature equation on the left and the second temperature equation on the right.

End values of DeltaT, using the first temperature equation on the left and the second temperature equation on the right.

The second graph shows the forcing at the end time, for both temperature equations. There is some change in the forcing for the first model, whereas for the second model the graph is quite different. Notice that the magnitude of the forcing term grows for the second model, due to the large value of \( \mu \).

One should also note that for the second temperature equation, a timestep restriction of the form \( \frac{\mu dt}{2 {dz}}^{2} < 1 \), which is required to maintain the stability of the \( \Delta T \) term. Combining this timestep restriction with the inflated value of \( \mu \) results in a severe restriction on \( dt \) and \( dz \) for the second model, whereas there is no timestep restriction for the first model.

From the analysis of the system of equations, we obtained the following a-priori estimate for \( q_{c}: \)

\[ || q_{c}(t) ||^{2}_{L^{2}} \, \leq e^{t} \int^{t}_{0} || C_{c}(s)||^{2}_{L^{2}} \text{ d} s . \]

However, this estimate holds when following the quantity \( q_{c} \) as it evolves along characteristics, whereas the Semi-Lagrange approach calculates the values of \(q_{c} \) at gridpoints which are fixed in space. As we have been using a velocity profile that is constant in time, it is possible to calculate exactly the trajectories of the characteristics. This allows us to use a change of variables to compute the norm \( || q_{c}(t) ||^{2}_{L^{2}} \), using the results that we have obtained using the Semi-Langrangian method. This allows us to test our results, and to confirm whether our simulated results correspond to the results that we have shown in the analysis of the model.