Milestone 6 - Fully Discrete EBM

Recap Semi-discretization

We derived the final form of the $2D$ EBM with heat diffusion in milestone 5 in the section about planetary heat transfer

where we here use the generic version in Cartesian coordinates for simplicity.

The so-called semi-discretization of the PDE (discretization in spatial coordinates), derived in milestone 5, is given for a grid node $(j,i)$ as

where $L_{j,i}$ denotes the spatial discretization of the diffusion operator $\vec{\nabla} \cdot (D\vec{\nabla} T)$. Recall, that for the spatial discretization, due to the pole problem of the spherical coordinates, a special treatment of the diffusion operator in the pole region is necessary. This form is called semi-discretization, as only the spatial directions are discretized, i.e., the formulation does not depend on the spatial coordinates $x$ anymore, but still on the temporal coordinate $t$.

The resulting problem is a (system) of ODEs that need to be integrated in time. In particular, we start with the equivalent form that we derived taking into account the linearity of the diffusion operator and the linearity of the radiative heat transfer,

Due to the linearity of the semi-discretization in the unknown variable (temperature $\mathbf{T}$ in case of the EBM), we can represent the spatial EBM operator $\mathbf{R}(\mathbf{T})$ as a constant matrix times the vector of unknown temperatures, $\mathbf{R}(\mathbf{T}) = \underline{\mathbf{ A}} \mathbf{T}$. The matrix $\underline{\mathbf{ A}}\in \mathbb R^{\texttt{NDOF}\times\texttt{NDOF}}$ is the so-called Jacobian matrix of the spatial operator and its computation is described in milestone 5.

Time Integration

We can now apply the two time-integration methods that we learned in milestone 3 to all equations of the linear system (3). For the evolution from time step $n$ ($t = t^n$) to time step $n+1$ ($t = t^{n+1}$), we discussed two options:

• If we use the forward Euler method, the fully discrete scheme reads as

and the solution at the time step $n+1$ can be obtained explicitly as

or eqvivalently

i.e. in the case of the explicit method, no Jacobian computation is necessary, instead the operator can be used directly. Thus, per time step, the update with the explicit Euler method is very quick - but the extreme time step restriction for (von Neumann) stability of the fully-discrete scheme of the EBM is to restrictive in practise: although one explicit time step is cheap, the amount of time steps to reach the equilibrium is so high (due to the restrictive time step stability) that the overall simulation takes a very long time.

• As a more efficient alternative, we discussed the backward Euler method, where the fully discrete scheme reads as

and the solution at the time step $n+1$, $\mathbf{T}^{n+1}\in \mathbb R^{\texttt{NDOF}}$, must be obtained by solving the system of linear equations

where $\underline{\mathbf{ I}}\in \mathbb R^{\texttt{NDOF}\times\texttt{NDOF}}$ is the identity matrix, $\underline{\mathbf{ M}}\in \mathbb R^{\texttt{NDOF}\times\texttt{NDOF}}$ is the system matrix of the linear equation system, and $\mathbf{Y}^{n+1}\in \mathbb R^{\texttt{NDOF}}$ is the right-hand side vector of the linear system.

Equilibrium simulation

Similar to the task in the pure ODE case, see milestone 3, we are seeking to simulate the equilibrium of the temperature field throughout the year. We emphasize again that we are not looking for a single steady-state temperature field - due to the annual dependence of the solar forcing, no global steady state solution can be reached. Instead, it is possible to reach an annual equilibrium state, where the solution changes in time, but repeats itself every year, such that we call it an equilibrium.

Again, the task is to define a proper stopping criterium when performing the temporal EBM simulations. One could again compute an annual average temperature and stop when the change of said temperature is smaller than a given threshold. One could also compare directly the full 2D temperature field throughout the year and define a proper norm to define a stopping criterium.

In summary, analogous to the pure ODE case, we have to compare the yearly temperature solutions with each other, until the difference from one year to the next is smaller than a given tolerance. Among the different choices of norms to compute the yearly solutions, the simplest option is to compute the Euclidean norm of the data vectors in time.