Modeling thermal evolution

Thermal evolution of lithosphere is often one of the key components of the long-term tectonics and is modeled by solving the heat equation:

$$\rho c_p\dot{T}=k{\nabla }^{2}T,$$

where $T$ is the temperature field while $c_p$ and $k$ are the heat capacity and the thermal conductivity of the lithosphere material. Multiplying by a weighting function on both sides and integrating by parts over the domain, we get

$$C_a{\dot{T}}_a^{t+\Delta t}=-\sum _e^{a\in e,b\in e}\left(k{D}_{ab}{T}_b^t{\Omega }_e\right)+\sum _e^{a\in s,s\in \partial {\Omega }_e}\left(\frac{1}{M-1}{\mathbf{q}}_{\mathbf{s}}\cdot {\mathbf{n}}_s{L}_s\right),$$

where the diffusion matrix

$$D_{ab}=\sum _{a,b\in e}\sum _i\frac{\partial N_a^e}{\partial x_i}\frac{\partial N_b^e}{\partial x_i}$$

is evaluated at the barycenter of each element since we use constant strain triangles (linear finite elements on simplexes). The lumped thermal capacitance (mass) is given by,

$$C_a=\sum _e^{a\in e}\left(\frac{1}{M}\rho c_p{\Omega }_e\right),$$

and ${\mathbf{q}}_{\mathbf{s}}$ is the prescribed boundary heat flux on a segment $s$. Then, the temperature is updated explicitly as:

$$T_a^{t+\Delta t}=T_a^t+\Delta t{\dot{T}}_a^{t+\Delta t}.$$

The stability condition for the explicit integration of temperature is usually satisfied by the time step size determined by the scaled wave speed, but if a stable time step size for heat diffusion is smaller, it becomes the global time step size.