Numerics and Physics

MagmaThermoKinematics.jl is built around a finite-difference energy solver with semi-Lagrangian advection and tracer-based tracking of injected sill material.

Governing Equation

The thermal state of the magmatic system is described by the advection–diffusion equation for temperature $T$ with latent heat release:

$$\begin{array}{}\text{(A1)}& \rho C_p(\frac{\partial T}{\partial t}+\mathbf{v}\cdot \mathrm{\nabla }T)=\mathrm{\nabla }\cdot (k\mathrm{\nabla }T)+\rho L\frac{\partial {\varphi }_s}{\partial t}+H\end{array}$$

where

Symbol	Description
$\rho$	Rock density (kg m⁻³)
$C_p$	Isobaric heat capacity (J kg⁻¹ K⁻¹)
$T$	Temperature (K or °C)
$t$	Time (s)
$\mathbf{v}$	Solid-state velocity field (m s⁻¹) induced by dike/sill emplacement
$k$	Thermal conductivity (W m⁻¹ K⁻¹)
$L$	Latent heat of crystallisation (J kg⁻¹)
${\varphi }_s=(1-\varphi )$	Solid fraction (dimensionless, 0–1), with $\varphi$ being the melt fraction
$H$	Radiogenic heat production (W m⁻³)

Operator Splitting

Equation (A1) is solved with a first-order operator split that separates advection from diffusion+latent heat. Each time step $\mathrm{\Delta }t$ proceeds in two sub-steps:

Step 1 — advection (semi-Lagrangian):

$$\begin{array}{}\text{(A2a)}& T^{\ast }=T(\mathbf{x}-\mathbf{v}\mathrm{\Delta }t,t^n)\end{array}$$

The departure point $\mathbf{x}-\mathbf{v}\mathrm{\Delta }t$ is found by back-tracking each grid node by one time step along the velocity field. The temperature at the departure point is obtained by bi-linear (2D) or tri-linear (3D) interpolation from the current grid values $T^n$. The semi-Lagrangian scheme is unconditionally stable with respect to the advective CFL criterion.

Step 2 — diffusion and latent heat:

$$\begin{array}{}\text{(A2b)}& \rho C_p\frac{T^{n+1}-T^{\ast }}{\mathrm{\Delta }t}=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }{T}^{n+1}\right)+\rho L\frac{\partial {\varphi }_s}{\partial t}{|}^{n+1}+H\end{array}$$

Effective Heat Capacity Formulation

Because melt fraction $\varphi$ depends on temperature, the latent heat term in (A2b) can be rewritten using the chain rule:

$$\rho L\frac{\partial {\varphi }_s}{\partial t}=\rho L\frac{\partial {\varphi }_s}{\partial T}\frac{\partial T}{\partial t}$$

Substituting into (A2b) and rearranging yields a modified heat capacity:

$$\begin{array}{}\text{(A3)}& \rho C_p^{\text{eff}}\frac{T^{n+1}-T^{\ast }}{\mathrm{\Delta }t}=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }{T}^{n+1}\right)+H\end{array}$$

where the effective heat capacity absorbs the latent heat contribution:

$$\begin{array}{}\text{(A3a)}& C_p^{\text{eff}}=C_p+L\frac{\partial \varphi }{\partial T}\end{array}$$

Where we made use of the fact that $\partial {\varphi }_s/\partial t=-\partial \varphi /\partial t$. This formulation is advantageous because it keeps the structure of the standard diffusion equation while automatically capturing latent heat release whenever $\partial \varphi /\partial T\ne 0$ (i.e., within the two-phase region). As this latter effect is taken into account in an implicit manner, it is numerically more stable.

An important point to keep in mind, though, is that $\partial \varphi /\partial T$ should be continous throughout the domain, including at the solidus and liquidus to prevent numerical instabilities. For this reason, the GeoParams package implements smoothening functions for any melting curve (as many parameterisations in use are discontinuous).

Explicit Time Discretisation

The diffusion operator in (A3) is discretised explicitly:

$$\begin{array}{}\text{(A4)}& T_{i,j}^{n+1}=T_{i,j}^{\ast }+\frac{\mathrm{\Delta }t}{\rho C_p^{\text{eff}}}[{(\mathrm{\nabla }\cdot k\mathrm{\nabla }T)}_{i,j}^n+H_{i,j}]\end{array}$$

Staggered-Grid Finite Differences

The spatial derivatives in (A4) are evaluated on a staggered grid:

• Temperature $T$ and scalar material properties ($\rho$, $C_p$, $\varphi$) are defined at cell centres.

• Heat fluxes $q=-k\mathrm{\nabla }T$ are defined on cell faces (half-integer indices).

• Conductivity $k$ is harmonically averaged to cell faces to conserve flux continuity across material interfaces.

The discrete Laplacian in 2D is:

$${(\mathrm{\nabla }\cdot k\mathrm{\nabla }T)}_{i,j}\approx \frac{k_{i+\frac{1}{2},j}(T_{i+1,j}-T_{i,j})-k_{i-\frac{1}{2},j}(T_{i,j}-T_{i-1,j})}{\mathrm{\Delta }{x}^2}+\frac{k_{i,j+\frac{1}{2}}(T_{i,j+1}-T_{i,j})-k_{i,j-\frac{1}{2}}(T_{i,j}-T_{i,j-1})}{\mathrm{\Delta }{z}^2}$$

with the analogous 3D extension for the $y$-direction.

CFL Stability Criterion

For the explicit thermal diffusion step the time step must satisfy the standard parabolic CFL condition:

$$\begin{array}{}\text{(CFL)}& \mathrm{\Delta }t\le \frac{1}{2}\frac{min(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z{)}^2}{max(\kappa )}\end{array}$$

where the thermal diffusivity is $\kappa =k/(\rho C_p^{\text{eff}})$. In practice the code evaluates $\kappa$ at every grid point and uses the global maximum to set a conservative $\mathrm{\Delta }t$.

Nonlinear Iterations (Picard)

Because $C_p^{\text{eff}}$ itself depends on $T^{n+1}$ through $\partial \varphi /\partial T$, equation (A4) is nonlinear. The solver uses damped Picard (fixed-point) iterations to resolve this nonlinearity within each time step:

Predictor — evaluate material properties at the current iterate $T^{(k)}$:

$$\begin{array}{}\text{(A5)}& C_p^{\text{eff},(k)}=C_p\left(T^{(k)}\right)-L\frac{\partial \varphi }{\partial T}{|}_{T^{(k)}}\end{array}$$

Corrector — advance temperature with the updated effective heat capacity:

$$\begin{array}{}\text{(A6)}& T_{i,j}^{(k+1)}=T_{i,j}^{\ast }+\frac{\mathrm{\Delta }t}{\rho C_p^{\text{eff},(k)}}[{(\mathrm{\nabla }\cdot k^{(k)}\mathrm{\nabla }{T}^{(k)})}_{i,j}+H_{i,j}]\end{array}$$

Damping — to aid convergence the update is damped:

$$T^{(k+1)}\leftarrow \alpha T^{(k+1)}+(1-\alpha )T^{(k)},0<\alpha \le 1$$

Convergence — iterations continue until the maximum pointwise temperature change falls below a prescribed tolerance $\epsilon$:

$$\underset{i,j}{max}|T_{i,j}^{(k+1)}-T_{i,j}^{(k)}|<\epsilon$$

Typically three to five Picard iterations are sufficient per time step.

Melt Intrusion Algorithms

Three distinct emplacement modes are available for adding new melt to the system:

Mode	Description
Geneva sill (under-accretion)	New magma is injected at the base of the existing melt body and solid host-rock is displaced downward.
UCLA-HD (central injection)	Melt is added at the centre of the intrusion; host rock is displaced radially outward while conserving volume.
Elastic dike	An elliptical dike geometry is used and host rock is displaced by an analytically prescribed elastic displacement field.

All three modes inject a Dike object, update the velocity field $\mathbf{v}$ for one time step (used in the semi-Lagrangian advection step A2a), and add tracer particles at the new melt location.

Tracers and Temperature–Time Paths

Passive tracer particles are advected with the host-rock velocity field using the same semi-Lagrangian scheme as the temperature field. At each time step each tracer records its current temperature, producing a continuous $T$–$t$ path. These paths are the primary input to the ZirconGrowth integration.

Tracers are kept on the CPU, as they generally use a lot of memory.

Dimensions and Geometry

Supported configurations include:

• 2D Cartesian.

• 2D axisymmetric (via specific workflows in examples).

• 3D Cartesian.

Material Parameter Handling

Material properties are integrated through GeoParams-based parameterizations and package-level helper routines for conductivity, density, heat capacity, melt fraction, and related quantities.

One of the highlights is that melting parameterizations and material properties can be changed very easily by using GeoParams (see ).

Changing Melting Parameterizations

Material setup is defined with SetMaterialParams, so changing physics is extremely simple and only requires you to change your MatParam struct. For example, you can switch between the Caricchi melting parameterization and a simple 4th-order polynomial melting curve with:

Caricchi

MatParam = (
	SetMaterialParams(Name="Rock", Phase=1,
		Density=ConstantDensity(ρ=2700kg/m^3),
		HeatCapacity=ConstantHeatCapacity(Cp=1000J/kg/K),
		Conductivity=T_Conductivity_Whittington_parameterised(),
		Melting=MeltingParam_Caricchi()
        ),
)

Smooth Polynomial

MatParam = (
	SetMaterialParams(Name="Rock", Phase=1,
		Density=ConstantDensity(ρ=2700kg/m^3),
		HeatCapacity=ConstantHeatCapacity(Cp=1000J/kg/K),
		Conductivity=T_Conductivity_Whittington_parameterised(),
		Melting=SmoothMelting(MeltingParam_4thOrder())
        ),
)
# Example alternatives:
# Conductivity=ConstantConductivity(k=3.3W/K/m)

This design makes it straightforward to test sensitivity to melt models and thermophysical assumptions without changing solver internals.

Performance and Parallelism

Backend selection through environment! configures CPU threads, or CUDA execution, and the package composes with ParallelStencil finite-difference modules.

• CPU and CUDA workflows are typically run with Float64 precision.

When your own script uses ParallelStencil macros directly (for example @zeros or @parallel), call @init_parallel_stencil(...) in script scope after environment!(...).

If you run this on a CPU, you can run it in parallel by starting julia in multi-threading mode:

$julia -t auto