All rock falls into three families, defined by how they formed — this matters for modelling because each family produces a different kind of 3-D geometry.

The Law of Superposition is the single most important rule for modelling: in an undisturbed sequence, depth = age. Three companion principles refine it:

A stratigraphic column is simply a vertical list of layers at a location, oldest at the bottom. It's the direct ancestor of a borehole log and is exactly the data structure a model's surfaces list represents.

Geological time itself is organized hierarchically — useful vocabulary when reading literature or maps:

A mineral is a naturally occurring, chemically-defined crystalline solid (e.g. quartz, feldspar, calcite). A rock is an aggregate of one or more minerals. For 3-D geological modelling you almost always work at the rock/lithology level, not individual minerals — mineral composition matters more for geochemistry and ore-grade estimation, which is a separate (though related) discipline from the geometric modelling this course focuses on.

Every tilted surface has two measurements that fully describe its 3-D orientation:

import numpy as np

def dip_to_normal(dip_deg, dip_dir_deg):
    """Convert dip angle + dip direction to a unit normal vector (Gx, Gy, Gz)."""
    dip = np.radians(dip_deg)
    dip_dir = np.radians(dip_dir_deg)
    gx = np.sin(dip) * np.sin(dip_dir)
    gy = np.sin(dip) * np.cos(dip_dir)
    gz = np.cos(dip)
    return (gx, gy, gz)

# Flat layer: dip=0 → (0, 0, 1) — points straight up
# Vertical layer dipping east: dip=90, dip_dir=90 → (1, 0, 0)

Folds form when layers bend under compressional stress rather than fracturing. The two basic shapes:

Fold geometry vocabulary

A fault is a fracture along which blocks of rock have moved relative to each other. Faults are classified by the direction of relative motion:

Key fault vocabulary

• Hanging wall — the block of rock above an inclined fault plane
• Footwall — the block below an inclined fault plane
• Throw — vertical component of displacement
• Heave — horizontal component of displacement
• Fault zone / damage zone — fractured rock surrounding the main fault surface

An unconformity is a buried surface representing a gap in time — older rock was deformed, uplifted, and eroded, then younger sediment was deposited on top of the erosion surface.

Igneous intrusions form when magma pushes into existing rock and solidifies underground — these cut across stratigraphy rather than following it.

Joints are fractures with no appreciable displacement — unlike faults, the rock on either side hasn't moved relative to the other side. They're geologically important (control groundwater flow, rock strength, blasting patterns in mining) but are almost never modelled as discrete surfaces in a regional 3-D model.

A geological map shows which rock unit is exposed at the surface across an area, using colour-coded polygons, plus strike/dip symbols and contact lines (boundaries between units).

A borehole is a vertical (or deviated) hole drilled into the ground; the resulting well log records which rock type was encountered at which depth. This is the single most common and most reliable source of surface_points data.

# Typical well log: a list of (top_depth, bottom_depth, lithology)
well_log = [
    (0,   15,  'Topsoil'),
    (15,  120, 'Limestone'),
    (120, 280, 'Shale'),
    (280, 500, 'Sandstone'),
]

collar_elevation = 1450  # surface elevation at the well (metres)
well_x, well_y = 2400, 1180   # well location

# Convert each TOP of a unit into a (X, Y, Z, surface) input point
# Z must be ELEVATION → collar_elevation minus depth
surface_points = []
for top_depth, bottom_depth, lith in well_log[1:]:  # skip topsoil top = ground level
    z = collar_elevation - top_depth
    surface_points.append((well_x, well_y, z, lith))

print(surface_points)
# [(2400, 1180, 1435, 'Limestone'), (2400, 1180, 1330, 'Shale'), (2400, 1180, 1170, 'Sandstone')]

Seismic reflection surveys send sound waves into the ground and record echoes off rock boundaries, producing a 2-D (or 3-D volume of) image where horizontal bands roughly correspond to layer boundaries. Common in oil & gas exploration.

A geological cross-section is a vertical "slice" through the Earth showing rock layers as they would appear if you cut the ground open along a line. It's the standard way geologists communicate a 3-D interpretation in 2-D, and it's exactly what tools like plot_2d() in GemPy generate.

An outcrop is a place where rock is exposed at the surface. Geologists visit outcrops with a compass-clinometer to directly measure strike and dip by placing the tool against the rock face — this is the most direct, ground-truth source of orientation data.

There are two fundamentally different approaches to building a 3-D geological model:

Most implicit modelling software (including GemPy) uses a form of kriging — specifically universal co-kriging with both surface points and orientation (gradient) constraints. You don't need the full mathematical derivation, but you do need the behavioural intuition:

• The interpolator fits a smooth scalar field where surface points all sit on the same isovalue (e.g. value = 0)
• Orientation data constrain the gradient (steepness/direction) of that field at given points
• Away from data, the field interpolates smoothly — it does not know about geological features you haven't told it about
• Sparse data + uniform orientations → smooth, "boring" surfaces; tight folds and abrupt features need denser sampling near them

A series groups surfaces that share one continuous depositional/deformation event. The order of series, and the relation between adjacent series, is how you encode the geological history (Tier 1 & 2 concepts) into the software:

Geological layers usually extend much further laterally than vertically — a limestone bed might span 10 km horizontally but be only 50 m thick. Anisotropy settings tell the interpolator about this directional bias so it doesn't treat horizontal and vertical distances equally.

Unlike a CAD model of a manufactured object, a geological model is always an interpretation of incomplete data — there is no ground truth to check against except more drilling. Understanding where uncertainty comes from is essential for using a model responsibly.

Structural restoration "undoes" deformation on a model — unfolding folds and closing faults — to check whether the pre-deformation state makes geological sense (e.g. layers should restore to roughly horizontal, without gaps or overlaps, which would indicate an impossible interpretation).

Real basins often have dozens of interacting faults rather than one isolated fault. Key concepts for handling this in software:

# Rows/columns = series, in computation order.
# A 1 at [i, j] means series i offsets/affects series j.
import numpy as np

fault_relations = np.array([
#                FaultA  FaultB  FaultC  Strat
    [0, 1, 1, 1],   # FaultA offsets B, C, and stratigraphy
    [0, 0, 1, 1],   # FaultB offsets C and stratigraphy
    [0, 0, 0, 1],   # FaultC offsets only stratigraphy
    [0, 0, 0, 0],   # Stratigraphy offsets nothing
])
geo_model.set_fault_relation(fault_relations)

Expert practice treats a geological model as a hypothesis to be tested, not a final answer. Standard checks:

1. Cross-validation against held-out data — remove some known points, run the model, check how close the prediction is to the held-out truth
2. Thickness sanity checks — verify modelled layer thickness stays within geologically reasonable bounds across the model
3. Stochastic/Monte Carlo runs — perturb uncertain inputs (e.g. borehole pick depths) within their error bounds and re-run many times to get a probability map of where each lithology likely occurs
4. Geological plausibility review — a domain expert visually checks the model against regional geological knowledge that isn't captured in the raw data

import numpy as np

n_realizations = 100
results = []

for i in range(n_realizations):
    # Perturb each borehole pick within its measurement uncertainty
    perturbed_z = original_z + np.random.normal(0, measurement_error_std)
    geo_model.modify_surface_points(point_id, Z=perturbed_z)
    sol = gp.compute_model(geo_model)
    results.append(sol.lith_block)

# Probability of a given voxel being a particular lithology:
results = np.array(results)  # shape: (n_realizations, n_voxels)
prob_sandstone = (results == SANDSTONE_ID).mean(axis=0)