STEADY-STATE HEAT CONDUCTION

1D Layered Geothermal Modeler

Edit layer boundaries, radioactive heat generation, and thermal conductivity. The analytical and finite-volume geotherms update immediately.

Load model

Surface temperature, T_s (°C) Basal heat flow, q_b (mW m⁻²) Total depth (km) Numerical cells Layer count

INPUT MODEL

Layer properties

A k Boundary

Selected layer Layer 1

A (µW m⁻³) k (W m⁻¹ K⁻¹)

Mouse: drag orange or blue vertical handles to change A or k; drag horizontal handles to change interface depth. Click a layer to select it.

MODEL RESPONSE

Temperature–depth profile

Analytical Finite volume

Surface heat flow—

Bottom temperature—

Maximum |ΔT|—

Thermal resistance—

Boundary conditions T(0) = T_s; k dT/dz at the model base = q_b. Depth z is positive downward, and heat-flow values are reported as upward-positive magnitudes.

NUMERICAL VALUES

Layer model

Edits in the table and plots are synchronized.

Layer	Top (km)	Bottom (km)	Thickness (km)	A (µW m⁻³)	k (W m⁻¹ K⁻¹)	Top q (mW m⁻²)	Top T (°C)

FORMULATION

Analytical and finite-volume solutions

The formulation is valid for an arbitrary number n of horizontal layers. Layer properties are piecewise constant, depth is positive downward, and heat flow is expressed as an upward-positive magnitude.

1. Layer notation and governing equation

Interfaces0 = z₀ < z₁ < ··· < z_n = L

Layer thicknessh_i = z_i − z_i−1

Layer propertiesA_i, k_i; i = 1, …, n

Numerical cellsM cells; Δz = L/M

\frac{d}{d z} [k (z) \frac{d T}{d z}] + A (z) = 0

(1)

q (z) = k (z) \frac{d T}{d z}

(2)

The boundary conditions are T(0) = T_s and q(L) = q_b. Temperature and heat flow remain continuous at every interface.

2. Exact piecewise solution for n layers

The heat flow at the top of layer i includes the basal contribution and the heat produced in layer i and all deeper layers:

q_{i}^{t} = q_{b} + \sum_{j = i}^{n} A_{j} h_{j}

(3)

Within layer i, define local depth ξ = z − z_i−1, where 0 ≤ ξ ≤ h_i. The exact solution is

q_{i} (ξ) = q_{i}^{t} - A_{i} ξ

T_{i} (ξ) = T_{i}^{t} + \frac{q_{i}^{t} ξ - \frac{1}{2} A_{i} ξ^{2}}{k_{i}}

(4)

The next layer starts from the values obtained at ξ = h_i:

q_{i + 1}^{t} = q_{i}^{t} - A_{i} h_{i}

T_{i + 1}^{t} = T_{i}^{t} + \frac{q_{i}^{t} h_{i} - \frac{1}{2} A_{i} h_{i}^{2}}{k_{i}}

(5)

3. Cell-centered finite-volume discretization

The interval 0 ≤ z ≤ L is divided into M cells. Each cell center receives the A and k values of its containing geological layer. Conductivity at an internal cell face is the harmonic mean:

k_{p + \frac{1}{2}} = \frac{2 k_{p} k_{p + 1}}{k_{p} + k_{p + 1}}, G_{p + \frac{1}{2}} = \frac{k_{p + \frac{1}{2}}}{Δ z}

(6)

For an interior cell p = 1, …, M−2, flux balance gives

- G_{p - \frac{1}{2}} T_{p - 1} + (G_{p - \frac{1}{2}} + G_{p + \frac{1}{2}}) T_{p} - G_{p + \frac{1}{2}} T_{p + 1} = A_{p} Δ z

(7)

The surface half-cell conductance is G_s = 2k₀/Δz. The first and last rows of the matrix are

(G_{s} + G_{\frac{1}{2}}) T_{0} - G_{\frac{1}{2}} T_{1} = A_{0} Δ z + G_{s} T_{s}

- G_{M - \frac{3}{2}} T_{M - 2} + G_{M - \frac{3}{2}} T_{M - 1} = A_{M - 1} Δ z + q_{b}

(8)

4. Calculation procedure

Validate and order the n layer interfaces, then convert km, µW m⁻³, and mW m⁻² to SI units.
Evaluate the exact piecewise solution by starting with q_b and propagating heat flow and temperature through all layers.
Assign A and k to M finite-volume cells, construct the tridiagonal system from Equations (6)–(8), and solve it with the Thomas algorithm.
Interpolate the analytical solution at numerical cell depths and report max |T_FV − T_analytical|. Increase M to assess grid convergence.

Computational cost: the tridiagonal solve requires O(M) operations and O(M) storage. The analytical calculation requires O(n + P) operations, where P is the number of plotted sample points.

Model equation: d/dz(k dT/dz) + A = 0. This tool represents one-dimensional, steady-state conduction; it does not include advection, transient cooling, phase changes, or temperature-dependent material properties.