BumiWasesa

Energi bumi yang lestari

📘 Theory & Documentation

Comprehensive guide to heat diffusion simulation theory, numerical methods, and best practices

← Back to Simulator

Table of Contents

  1. Governing Equation
  2. Numerical Method
  3. Stability & CFL Condition
  4. Material Properties Reference
  5. Horizontal Layer Cake Model
  6. Tips & Tricks
  7. References

1. Governing Equation

The Heat Diffusion Equation

The 2D heat diffusion in geological media is governed by:

∂T/∂t = κ(∂²T/∂x² + ∂²T/∂z²) + Q

Where:

  • T(x, z, t) - Temperature field [°C or K]
  • t - Time [s]
  • x, z - Spatial coordinates (horizontal, vertical) [m]
  • κ - Thermal diffusivity [m²/s]
  • Q - Heat source term [K/s]

Thermal Diffusivity

The thermal diffusivity κ relates to material properties:

κ = k / (ρ · cp)
  • k - Thermal conductivity [W/(m·K)]
  • ρ - Density [kg/m³]
  • cp - Specific heat capacity [J/(kg·K)]

Heat Source Term

Volumetric heat generation S [W/m³] is converted to temperature rate:

Q = S / (ρ · cp) [K/s]

This ensures physically realistic temperature increases. For example, with granite properties (ρ=2700 kg/m³, cp=790 J/(kg·K)), a heat source of 1 W/m³ produces a temperature increase rate of 4.69×10⁻⁷ K/s.

2. Numerical Method: Explicit Finite Difference

Spatial Discretization

The continuous domain is discretized into a uniform grid:

  • Grid points: nx × nz
  • Spatial steps: Δx and Δz [m]
  • Temperature at grid point (i,j): Ti,j

Finite Difference Approximations

Second-order spatial derivatives are approximated using central differences:

∂²T/∂x² ≈ (Ti+1,j - 2Ti,j + Ti-1,j) / Δx²
∂²T/∂z² ≈ (Ti,j+1 - 2Ti,j + Ti,j-1) / Δz²

Explicit Time Integration

Forward Euler method is used for time advancement:

Ti,jn+1 = Ti,jn + Δt · [κ·∇²Ti,jn + Qi,j]

Why Explicit Method for Geological Time Scales?

The explicit method is ideal for geological simulations because:

  • ✅ Low thermal diffusivity (κ ~ 10⁻⁶ m²/s) allows large time steps (~40 years)
  • ✅ Few time steps needed (1000 years ≈ 25 steps with default parameters)
  • ✅ Simple implementation - no matrix inversion required
  • ✅ Low memory footprint
  • ✅ Fast per-step computation

Algorithm Steps

  1. Initialize temperature field: T0 = Tbackground
  2. Calculate stable time step: Δt based on CFL condition
  3. For each time step n = 0, 1, 2, ..., nt:
    • Apply boundary conditions (zero flux at all boundaries)
    • Compute spatial derivatives (∇²T) using finite differences
    • Add heat sources at specified locations and times
    • Update temperature: Tn+1 = Tn + Δt·(κ∇²T + Q)
    • Save snapshots at specified intervals

3. Stability & CFL Condition

CFL Stability Criterion

For the explicit scheme to remain stable, the time step must satisfy:

CFL Condition:
Δt ≤ (Δx² · Δz²) / (2κ(Δx² + Δz²))

For square grids (Δx = Δz):

Δtmax = Δx² / (4κ)

Safety Factor

BumiWasesa uses a 0.5 safety factor to ensure robust stability:

Used in simulation:
Δt = 0.5 × Δtmax

This accounts for:

  • Floating-point rounding errors
  • Spatial discretization errors
  • Source term interactions
  • Boundary condition approximations

Example Calculation

With default parameters (Δx = Δz = 100 m, κ = 1×10⁻⁶ m²/s):

Δtmax = (100)² / (4 × 1×10⁻⁶) = 2.5×10⁹ s ≈ 79 years
Δtsafe = 0.5 × 2.5×10⁹ s = 1.25×10⁹ s ≈ 40 years

For a 1000-year simulation: ~25 time steps - very efficient!

Truncation Error

The numerical scheme has:

  • Second-order accuracy in space: O(Δx²)
  • First-order accuracy in time: O(Δt)

To improve accuracy, reduce grid spacing (Δx, Δz) rather than just the time step.

4. Material Properties Reference

Geological Materials

Material Density ρ
(kg/m³)		 Specific Heat cp
(J/(kg·K))		 Thermal Conductivity k
(W/(m·K))		 Thermal Diffusivity κ
(m²/s)
Granite 2700 790 2.5-3.0 1.0-1.5 × 10⁻⁶
Sandstone 2400 920 2.0-3.0 1.0-1.5 × 10⁻⁶
Basalt 2900 840 1.5-2.5 0.6-1.0 × 10⁻⁶
Limestone 2600 900 2.0-3.0 0.8-1.3 × 10⁻⁶
Clay/Shale 2200 900 1.0-2.0 0.5-1.0 × 10⁻⁶
Quartzite 2650 1000 3.0-6.0 1.1-2.3 × 10⁻⁶
Gneiss 2700 800 2.0-3.5 0.9-1.6 × 10⁻⁶
Slate 2700 760 2.0-2.5 1.0-1.2 × 10⁻⁶
Marble 2700 880 2.5-3.0 1.0-1.3 × 10⁻⁶
Water 1000 4186 0.6 1.4 × 10⁻⁷
Air (20°C) 1.2 1005 0.026 2.2 × 10⁻⁵

Volumetric Heat Capacity

The volumetric heat capacity (ρ·cp) determines how much energy is needed to raise temperature:

Material ρ·cp
(MJ/(m³·K))		 Notes
Granite 2.13 Default in simulator
Sandstone 2.21 Common reservoir rock
Basalt 2.44 Volcanic rock
Limestone 2.34 Sedimentary rock
Water 4.19 High heat capacity

Typical Geothermal Gradients

Setting Gradient (°C/km) Heat Flow (mW/m²)
Continental crust (average) 25-30 60-80
Volcanic zones (high) 50-100 100-200
Oceanic crust (young) 30-50 80-120
Old continental shield 15-20 40-50

5. Horizontal Layer Cake Model

Overview

The simulator supports optional horizontal layering, allowing you to model stratified geological structures with different thermal properties at different depths.

How It Works

  • Layers are stacked from top to bottom: The first layer starts at the surface (z_min), and subsequent layers are stacked below it.
  • Each layer has four properties: thickness (m), thermal diffusivity κ (m²/s), density ρ (kg/m³), and specific heat c_p (J/(kg·K)).
  • The last layer is a halfspace: It extends to infinite depth, representing the deep crust or mantle.
  • Properties vary with depth: The simulator automatically assigns properties based on depth, creating horizontal property maps.

Example Configuration

Three-layer model (sediments → crust → mantle):

  • Layer 1 (Sediments): thickness = 3000 m, κ = 1.0×10⁻⁶ m²/s, ρ = 2400 kg/m³, c_p = 1000 J/(kg·K)
  • Layer 2 (Upper Crust): thickness = 12000 m, κ = 1.2×10⁻⁶ m²/s, ρ = 2700 kg/m³, c_p = 790 J/(kg·K)
  • Layer 3 (Lower Crust/Halfspace): thickness = ∞, κ = 1.0×10⁻⁶ m²/s, ρ = 2900 kg/m³, c_p = 820 J/(kg·K)

When to Use Layers

✅ Use layers when:

  • Modeling sedimentary basins with distinct lithological units
  • Studying heat flow through crust with different thermal properties
  • Investigating geothermal systems in layered geological structures
  • You have depth-dependent property data from geological surveys

ℹ️ Use uniform properties (no layers) when:

  • Modeling homogeneous rock masses (e.g., granite pluton)
  • Doing initial exploratory simulations
  • Properties don't vary significantly with depth in your study area
  • You want simpler model setup and faster computation

Important Notes

  • Layers create horizontal discontinuities in thermal properties
  • The simulator uses the exact depth to determine which layer's properties apply
  • Heat still diffuses vertically across layer boundaries
  • For best results, ensure layer thicknesses are much larger than your grid spacing (dz)
  • The geothermal gradient applies continuously across all layers

6. Tips & Tricks for Parameter Selection

🎯 Spatial Resolution

  • Rule of thumb: 5-10 grid points across smallest feature
  • Heat sources: If your heat source is 500m wide, use Δx ≤ 100m
  • Computational cost: Halving grid spacing quadruples computation time
  • Typical values: 50-200 m for geothermal simulations

⏱️ Time Step Selection

  • Automatic: Let the simulator calculate Δt from CFL condition
  • Time scale: Thermal diffusion time ~ L²/κ
    • L = 1 km, κ = 10⁻⁶ m²/s → ~32 years
    • L = 10 km, κ = 10⁻⁶ m²/s → ~3200 years
  • Simulation duration: Run 2-3× the thermal time scale for near-steady state

🌍 Domain Size Selection

  • Include all sources: Add 2-5 km margin around heat sources
  • Boundary effects: Ensure boundaries don't affect interior solution
  • Check: Look at edge temperatures - should be close to background
  • Typical geothermal: 10-20 km horizontal, 10-30 km vertical

🔥 Heat Source Configuration

  • Intensity: Typical geothermal: 0.1-10 W/m³ (not MW!)
  • Volcanic intrusions: 1-100 W/m³ for active magma chambers
  • Timing: Use t_start and t_end for transient events (intrusions, extraction)
  • Spatial extent: Use elliptical sources for realistic geometries

✅ Quality Checks

  1. Mass/energy balance: Total energy should be conserved (within numerical error)
  2. Stability: No oscillations or unbounded temperature growth
  3. Symmetry: Symmetric inputs should produce symmetric outputs
  4. Boundary independence: Interior solution shouldn't change with larger domain
  5. Resolution convergence: Results shouldn't change significantly with finer grid
  6. Physical bounds: Temperatures should be realistic (not billions of degrees!)

⚠️ Common Pitfalls

  • ❌ Using high thermal diffusivity: κ > 10⁻⁵ m²/s requires too many time steps
  • ❌ Domain too small: Boundary conditions affect interior solution
  • ❌ Grid too coarse: Missing important features, poor accuracy
  • ❌ Heat sources too strong: Unrealistic temperatures (check units: W/m³, not MW!)
  • ❌ Simulation too short: System hasn't reached equilibrium yet

7. References

Heat Transfer Theory

  1. Carslaw, H.S. & Jaeger, J.C. (1959). Conduction of Heat in Solids. Oxford University Press.
    Classic text on analytical solutions to heat conduction problems.
  2. Incropera, F.P., DeWitt, D.P., Bergman, T.L., & Lavine, A.S. (2011). Fundamentals of Heat and Mass Transfer (7th ed.). Wiley.
    Comprehensive engineering textbook covering all modes of heat transfer.
  3. Ozisik, M.N. (1993). Heat Conduction (2nd ed.). Wiley-Interscience.
    Detailed treatment of heat conduction in various geometries and materials.

Numerical Methods

  1. LeVeque, R.J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. SIAM.
    Essential reference for finite difference methods and stability analysis.
  2. Morton, K.W. & Mayers, D.F. (2005). Numerical Solution of Partial Differential Equations (2nd ed.). Cambridge University Press.
    Comprehensive coverage of PDE numerical methods including error analysis.
  3. Strikwerda, J.C. (2004). Finite Difference Schemes and Partial Differential Equations (2nd ed.). SIAM.
    Advanced treatment of finite difference methods for PDEs.
  4. Press, W.H., Teukolsky, S.A., Vetterling, W.T., & Flannery, B.P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press.
    Practical algorithms for scientific computing including parabolic PDEs.

Geothermal Applications

  1. Beardsmore, G.R. & Cull, J.P. (2001). Crustal Heat Flow: A Guide to Measurement and Modelling. Cambridge University Press.
    Essential guide to geothermal heat flow measurements and modeling.
  2. DiPippo, R. (2015). Geothermal Power Plants: Principles, Applications, Case Studies and Environmental Impact (4th ed.). Butterworth-Heinemann.
    Comprehensive treatment of geothermal energy systems and applications.
  3. Grant, M.A., Donaldson, I.G., & Bixley, P.F. (1982). Geothermal Reservoir Engineering. Academic Press.
    Classic text on geothermal reservoir modeling and management.
  4. Turcotte, D.L. & Schubert, G. (2014). Geodynamics (3rd ed.). Cambridge University Press.
    Comprehensive treatment of Earth's thermal structure and dynamics.

Computational Geosciences

  1. Wang, H.F. & Anderson, M.P. (1995). Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods. Academic Press.
    Practical guide to subsurface flow and transport modeling.
  2. Cheng, A.H.-D. & Cheng, D.T. (2005). Heritage and early history of the boundary element method. Engineering Analysis with Boundary Elements, 29(3), 268-302.
    Historical perspective on numerical methods in engineering.
  3. Gerya, T. (2019). Introduction to Numerical Geodynamic Modelling (2nd ed.). Cambridge University Press.
    Modern treatment of numerical methods for geodynamic problems.

Online Resources

← Back to Simulator