The Geodynamic Modelling ToolBox is a Julia package primarily intended for teaching purposes. It provides various finite difference, staggered discretization schemes to numerically solve the governing equations of two-dimensional geodynamic problems. These include the conservation equations of:
GeoModBox.jl includes a series of exercises and examples of geodynamically well-defined problems. The exercises are provided as Jupyter notebooks for students to complete. The theoretical background is documented here.
The solvers for each governing equation can be used separately or in combination for dimensional or non-dimensional problems, with only minimal modifications when calling the functions. For more informations on how to use the individual functions please see the list of functions or individual exmples. Some typical initial conditions, such as a linearly increasing temperature, are predefined and can be called using specific functions. In the following a brief explenation is given regarding the governing equations and the numerical method to solve them within the GeoModBox.jl. For more detailed information see the individual documentations.
To properly solve the governing equations, a staggered finite difference scheme is chosen for the energy and momentum equations. A staggered grid enables a correct and straightforward implementation of boundary conditions and ensures conservation of stress between nodes in cases of variable viscosity. This requires certain parameters to be defined on different grids. For more information regarding the physical and numerical background, please refer to this.
Within the GeoModBox.jl, temperature, density, pressure, normal deviatoric stresses, and heat production rate are defined on the centroids. The deviatoric shear stresses are defined on the vertices, and velocities are defined between the vertices. Viscosity is required on both.
For further details on the implementation in GeoModBox.jl, see the individual documentations for each governing equation.
In geodynamics, the energy is described by the temperature and needs to be conserved within a closed system. Within the GeoModBox.jl, the temperature conservation equation, or temperature equation, is solved using an operator splitting method, that is, first the advective part of the temperature equation is solved, followed by the diffusive part.
GeoModBox.jl provides several finite difference schemes for solving the diffusive part of the time-dependent or steady-state temperature equation, including radioactive heating, in both 1-D and 2-D. The solvers are located in src/HeatEquation. Currently, only Dirichlet and Neumann thermal boundary conditions are supported. Most functions assume constant thermal parameters (with the exception of the 1-D solvers and the 2-D, iterative implicit solver, called iterative defection correction method).
GeoModBox.jl provides various methods to advect properties within the model domain. The routines are structured so that any property defined on centroids (including ghost nodes at all boundaries) can be advected using the described solvers. Using passive tracers, one may choose to advect either the absolute temperature or the phase ID.
On geological timescales, Earth's mantle and lithosphere deform slowly due to their high viscosity, allowing us to neglect inertial forces. This simplifies the Navier-Stokes equation into the Stokes equation. GeoModBox.jl provides two main methods to solve the Stokes equation in 1-D and 2-D: the direct method and the defection correction method, applicable for both constant and variable viscosity fields. Velocity and pressure are defined on a staggered grid, and ghost nodes are included to ensure proper implementation of free-slip and no-slip boundary conditions.
The following are visualizations of selected examples provided by GeoModBox.jl. For further details, refer to the documentation linked in each title.
Figure 1. Gaussian Diffusion. Time-dependent, diffusive solution of a 2-D Gaussian temperature anomaly at a resolution of 100 × 100, using the Crank-Nicholson approach, compared to the analytical solution.
Top Left: 2-D temperature field with numerical isotherms (solid black) and analytical isotherms (dashed yellow).
Top Right: Total deviation from the analytical solution.
Bottom Left: 1-D y-profile along
Bottom Right: Root Mean Square (RMS) total deviation over time.
Figure 2. Resolution test. Maximum RMS error
Figure 3. Rigid-Body-Rotation. Time-dependent advection of a rotating circular temperature anomaly using the upwind (top), semi-Lagrangian (middle), and tracer (bottom) methods on a 100 × 100 grid. Within a circular region, the velocity field follows rigid rotation; outside, it is zero. Temperature for tracers is interpolated to the grid for visualization but not updated on the tracers.
Figure 4. Isoviscous Falling Block. Time-dependent simulation of an isoviscous falling block at 50 × 50 resolution with 9 tracers per cell. The solver handles variable viscosities. Tracers advect the phase ID, which is used to interpolate density and viscosity on centroids and vertices, respectively.
Figure 5. Falling Block Sinking Velocity. Block sinking velocity vs. initial viscosity ratio
Figure 6. Falling Block Benchmark. Tracer distribution at the final stage for selected viscosity ratios
Figure 7. Rayleigh-Taylor Instability. Evolution of two-layered Rayleigh-Taylor instability.
Figure 8. Bottom-Heated, Isoviscous Convection for Ra =
TOP: Transient temperature field with velocity vectors.
BOTTOM: Horizontally averaged temperature–depth profiles at each time step.
Solvers: defect correction (momentum), semi-Lagrangian (advection), Crank-Nicolson (heat diffusion).
Boundary conditions: Dirichlet (top/bottom), Neumann (sides), free-slip (velocity, all sides).
Figure 9. Internally Heated Convection for
Same setup as above, but with Neumann boundary at the bottom (zero heat flux) and constant internal volumetric heat production
Figure 10. Mixed-Heated Convection for Ra =
Combination of the above two setups (bottom heating + internal heating).