Energy Conservation Equation
The conservation of energy is a fundamental physical principle stating that energy can neither be created nor destroyed, but only transformed. In geodynamical modeling, this principle is typically expressed in terms of temperature, which evolves due to diffusive and advective transport processes. Assuming radioactive heat sources only, the general energy equation is defined as:
\[\begin{equation} \left(\frac{\partial E}{\partial t} + v_j\frac{\partial{E}}{\partial{x_j}}\right) + \frac{\partial q_{i}}{\partial x_{i}} = \rho H, \end{equation}\]
where energy is defined as $E = c_p \rho T$. Here $c_p$ is the specific heat capacity [J/kg/K], $\rho$ is the density [kg/m³], $T$ is temperature [K], $\partial/\partial t$ is the time derivative, $t$ is time [s], $v_j$ is the velocity [m/s] in direction $j$, $q_i$ is the heat flux [W/m²] in direction $i$, $\partial/\partial x_i$ is a directional derivative in $i$, $H$ is the internal heat production per unit mass [W/kg]. Repeated indices imply summation.
The heat flux $q_i$ is described by Fourier’s law:
\[\begin{equation} q_i = - k \frac{\partial{T}}{\partial{x_i}}, \end{equation}\]
where $k$ is the thermal conductivity [W/m/K]. The flux is directed opposite to the temperature gradient and represents the amount of heat passing through a unit surface per unit time.
Substituting Fourier’s law into the energy equation yields the temperature conservation equation in Eulerian form:
\[\begin{equation} \rho c_p \left(\frac{\partial T}{\partial t} + v_j\frac{\partial{T}}{\partial{x_j}}\right) = -\frac{\partial q_i}{\partial x_i} + \rho H. \end{equation}\]
This equation captures temperature changes due to diffusion (right-hand side) and advection (left-hand side). For simplicity and assuming a spatially constant internal heat production, these processes can be split using an operator splitting technique, solving the advection and diffusion steps sequentially. If internal heat production varies spatially, a more advanced advection scheme is required to account for source term integration.
Heat Diffusion Equation
GeoModBox.jl provides several finite difference (FD) schemes to solve the diffusive component of the time-dependent or steady-state temperature equation—including optional radioactive heating—in both 1-D and 2-D. Available methods include:
- Forward Euler
- Backward Euler
- Crank–Nicolson
- Alternating Direction Implicit (ADI)
See the documentation for the 1-D and 2-D solvers for detailed descriptions of each method.
Currently, only Dirichlet and Neumann boundary conditions are supported. Most implementations assume constant thermal properties, although certain 1-D and 2-D solvers allow for variable parameters. See the HeatEquation source directory for implementation details.
Examples
- 1-D oceanic geotherm
- 1-D continental geotherm
- Comparison of FD schemes on a Gaussian anomaly
- 2-D resolution test with Gaussian anomaly
- 2-D Poisson equation resolution test
For explanations, see the examples documentation and the full example directory.
Exercises
- 1-D forward Euler
- 1-D backward Euler
- 2-D Poisson equation
- 2-D transient plume heating
- 2-D transient sill heating
Heat Advection Equation
To solve the advective component of the temperature equation, GeoModBox.jl offers several schemes:
- Upwind scheme
- Staggered leapfrog scheme
- Semi-Lagrangian scheme
- Passive tracers/markers
See the advection documentation and associated source code for method-specific details. Tracer-related functionality is located in src/Tracers, whereas other advection schemes are implemented in src/AdvectionEquation.
Examples
See the examples documentation for further details.