Five groundwater models, built one on top of the other, running live in your browser. In each module you first work a number out by hand, then type it in, then run the model and watch your number appear in the water budget.
How this works. The solver behind every module is a finite-difference groundwater engine — it solves the same equations as MODFLOW, the U.S. Geological Survey's standard groundwater modeling code, on a 21 × 11 grid of 100 m cells. Every reference number here was cross-checked against real MODFLOW 6.8.0 runs of the same models (open the browser console to see the self-test). Checkpoints are self-assessment only; nothing is recorded, and progress resets if you reload.
One row of 21 cells, water held at 10 m on the left and 0 m on the right. Small enough to solve on paper — which is exactly why we start here: if the model can't reproduce your algebra, nothing bigger can be trusted.
Hydraulic head h is the level water would rise to in a pipe stuck into the aquifer — a direct readout of the water's energy, in metres. Groundwater always moves from higher head to lower head, and Darcy's law says how fast:
K is hydraulic conductivity — how easily the sediment lets water through (fine sand ≈ 1–10 m/d). b is the saturated thickness, W the width of the slab, and Δh/L the gradient — the slope of the energy line. In this module the layer is confined: b is fixed at the full 10 m no matter what the head does, so everything stays beautifully linear.
Our aquifer: K = 10 m/d, thickness b = 10 m, one row W = 1 m wide. The two end cells are constant-head boundaries (MODFLOW's CHD package): head is pinned at 10.0 m and 0.0 m at the cell centres, which sit 2000 m apart.
How much water flows through this aquifer, in m³/day? Use Darcy's law with the numbers above.
The 21 cell centres run from x = 50 m to x = 2050 m. If head drops linearly from 10 m to 0 m, what head should the model compute at the middle cell (cell 11, x = 1050 m)?
This is the exact model, solved by finite differences the way MODFLOW does it. If your two checkpoint answers were right, they should appear — to the decimal — in the tiles below.
Two things to try. Double K and watch the flow double while the heads don't move at all — in a purely head-driven problem the profile is set by the boundaries; K only scales the throughput. Then raise the east boundary and watch the whole line pivot. The budget-error tile is the modeler's seatbelt: in ≈ out, always, or something is wrong.
Before you touch the slider: if K doubles from 10 to 20 m/d, what happens to the head at cell 11?
Real aquifers are modelled with the same ingredients you just used: a grid, a K value per cell, and boundary conditions. In MODFLOW's vocabulary the pinned end cells are the CHD (constant-head) package. A MODFLOW 6.8.0 run of this exact model returns the linear profile to within 10⁻¹⁵ m and a water-budget discrepancy of 0.00% — reproducing a hand solution to machine precision is the standard first test of any modeling tool, and now you have performed it yourself.
The same aquifer, widened to 11 rows — a 2.1 km × 1.1 km slab — with a well pumping from dead centre. The head field becomes a landscape, and the budget starts answering the only question that matters: whose water is the well taking?
In 2-D it gets tedious to carry K and b separately, so hydrogeologists multiply them once and give the product a name:
With K = 10 m/d and a 10 m thick confined layer, what is the transmissivity T?
Before any well: the aquifer is now 11 rows × 100 m = 1100 m wide. How much water flows through the whole slab, west to east, in m³/d? (Same Darcy's law as Module 1 — only W changed.)
Drag the rate up and watch two things at once: the shape (heads dimple into a cone around the well) and the budget (the boundaries rearrange themselves to feed it).
Read the budget at 500 m³/d. Inflow rises 550 → 800 (+250) and outflow shrinks 550 → 300 (−250). The well's 500 is exactly half "new water pulled in" and half "water that used to leave" — the split is symmetric because the well sits dead centre. Hydrogeologists call this capture. Steady state has no storage: every drop the well takes is a drop some boundary gives up.
At steady state, where does a pumping well's water come from?
In MODFLOW a pumping well is the WEL package: a single line of input naming the cell and the rate, negative for extraction — everything else you watched happen (the cone, the capture) is the physics responding to that one line. The MODFLOW 6.8.0 run of this model gives a well-cell head of 0.993 m and the 800 / −300 budget you just reproduced. Agreement to the third decimal between two independent solvers is what "verified" means in modeling.
Flip one switch — the layer is now unconfined — and the model turns nonlinear. The saturated thickness is the water table, so pumping erodes the very property that delivers water to the well. This module is about what that feedback does.
Unconfined, the transmissivity is T = K · h — thickness equals the water-table height above the aquifer bottom. Drop the water table and T drops with it. Two consequences:
1 · The natural profile curves. With head 10 m on the west and 2 m on the east, the east side is thin and slow, so it needs a steeper gradient to carry the same flow. The result is the classic Dupuit parabola:
2 · Wells can fail suddenly. Pumping thins the aquifer at the well → less T → less inflow → more drawdown → thinner still. Past a threshold the loop runs away and the cell dries out.
Why the east boundary moved to 2 m. In Modules 1–2 it was 0 m — the aquifer bottom. Unconfined, that would mean zero saturated thickness at the outlet: T = 0, a wall no water can cross. A water-table boundary must stay wet, so this model holds the east side at 2 m.
Use the Dupuit formula to predict the natural (no pumping) water-table height at the well column — x = 1000 m along the 2000 m between boundaries, with h₁ = 10 m, h₂ = 2 m. Give h, not h².
The solver below uses MODFLOW's standard formulation, where drying is abrupt. Raise the rate slowly and watch the delivered water and the well-cell head. Somewhere there is a cliff.
The confined aquifer of Module 2 delivered 500 m³/d without complaint. Why does the same aquifer, unconfined, fail near 280?
In MODFLOW, confined versus unconfined is a single flag per cell (ICELLTYPE: 0 confined, 1 convertible) — the entire nonlinearity you just wrestled with switches on with one input value. A MODFLOW 6.8.0 bisection on this model put the maximum sustainable yield at 282.5 m³/d — with a band just above it where the solver fails to converge at all, and a false-calm zone above ~294 where it "converges" to a dry cell and a well delivering zero. Professional water-table models use the Newton–Raphson method to ride through drying gracefully; the solver behind Modules 4 and 5 does the same.
Steady state is a convenient fiction — real aquifers take years to get there, and the journey is paid for out of storage. This module turns time on: the well starts at day zero against the natural water table, and you watch where each day's water actually comes from.
When the water table drops, the sediment above it doesn't empty completely — capillary forces keep some water behind. The fraction that actually drains is the specific yield Sy (≈ 0.10–0.25 for sands). It converts water-table movement into volume:
One model cell is 100 m × 100 m and Sy = 0.15. How many m³ of water are released if the water table under that one cell drops by 1 m?
Press Run: the solver computes 100 growing time steps across 10 years, then plays them back. The chart splits every day's pumping into its two possible sources.
The lesson in the crossover. Early on, storage supplies almost everything — the boundaries haven't felt the well yet. Over months to years the cone reaches them and capture takes over; at true steady state, storage's share is zero. At 500 m³/d something else happens too: near year 3 the water table reaches the aquifer bottom at the well and physics forces the rate down toward ~400 m³/d — the model version of a well that "can't make its rate" anymore. Overpumping is invisible exactly as long as storage is picking up the bill.
If you cut Sy from 0.15 to 0.05 (a finer-grained aquifer) and rerun, the day the well first fails to make its full rate will…
In MODFLOW, turning on time means declaring a transient stress period and supplying storage properties through the STO package — Sy becomes an input alongside K, and it is the parameter a modeler calibrates against observed drawdown. Water-table transients also need the Newton–Raphson formulation, which turns wells down smoothly instead of crashing as cells approach dry. The MODFLOW 6.8.0 run of this model delivers the full 500 m³/d until day ~1064, then reduces to 403 m³/d by year 10 — the same curve your run traces.
A stream crosses the aquifer 300 m east of the well. It talks to the groundwater through its bed — both directions. This is the module where "sustainable yield" stops being a number and becomes a question: sustained at whose expense?
MODFLOW's river (RIV) boundary exchanges water at a rate proportional to the head difference across the streambed:
The conductance C bundles the streambed's geometry and permeability. For one model cell's worth of stream:
If the water table ever drops below the streambed bottom the stream disconnects: leakage caps at C · (stage − bottom) no matter how much deeper the aquifer is pumped.
Our stream: bed conductivity 0.1 m/d, reach length through one cell 100 m, width 10 m, bed thickness 0.5 m. Compute the conductance per cell.
Before pumping, the model's water table at the stream is 5.064 m and the stream stage is 5.0 m. How much water does ONE stream cell gain from the aquifer, in m³/d? (Give the magnitude.)
Press Run. The chart tracks the stream–aquifer exchange on a log time axis — the way stream-depletion curves are always plotted — while the map shows the cone reaching the stream.
Three acts. (1) For weeks the exchange barely moves — the well is living on storage and the stream hasn't noticed. (2) Near day 140 the curve crosses zero: the whole reach flips from gaining to losing. (3) It levels off with the stream donating ~217 m³/d and the well finally making its full 500 — the stream's support is precisely what makes the "unsustainable" rate of Module 3 sustainable. Total depletion ≈ 142 lost baseflow + 217 induced leakage ≈ 359 m³/d, 72% of the pumping, invoiced to whoever is downstream.
A regulator asks: "The well is 300 m from the stream — is the stream affected?" Based on this model, the honest answer is:
In MODFLOW the stream is the RIV (river) package, and its inputs are exactly the three numbers you computed by hand: stage, conductance, and bed bottom. MODFLOW 6.8.0 numbers for this model: natural gain 141.6 m³/d, flip at day 142, equilibrium leakage 216.6 m³/d — check them against your tiles.
Where modelers go from here: a streamflow-routing package that tracks the river's own discharge shrinking reach by reach; moving the well and watching the flip day scale with distance²; recharge, seasons, more layers. Each is one more package on the five-model backbone you just climbed.