gwsource tutorials#

Jupyter notebooks that teach cell-mixing flow tracing with gwsource, from a first capture zone to transient tracking and parameter ensembles. Each notebook builds its own small MODFLOW 6 model and runs in a few seconds to about three minutes (the table below gives the run time of each). The notebooks are committed with their outputs, so they read on GitHub without running them.

Setup and launch#

Create the gw-source environment and download MODFLOW 6 into bin/ (once). You need Miniforge (mamba and conda) or micromamba. The commands are the same on Linux, macOS and Windows. On Windows, run them in the Miniforge Prompt. Run them from the repository root, because environment.yml installs gwsource from this folder:
```
mamba env create -f environment.yml
mamba activate gw-source
python scripts/setup_env.py
```
scripts/setup_env.py installs FloPy from its develop branch and downloads MODFLOW 6, triangle and PEST++ into bin/. With conda or micromamba, write conda or micromamba in place of mamba. To update an existing environment, run mamba env update -n gw-source -f environment.yml instead of mamba env create.

The gw-source environment includes JupyterLab and the notebook tools (environment.yml). To add them to an existing environment: mamba install -n gw-source -c conda-forge jupyterlab ipykernel nbclient nbformat nbconvert.

Tutorial 04 and section 5 of tutorial 11 also need vorflow (Hugman and Sanchez 2026), which builds their Voronoi grid of the tutorial aquifer. It is an optional dependency of the tutorials only, not of gwsource. environment.yml installs it. To add it to an existing environment: pip install vorflow.

Tutorial 09 also needs PEST++ (bin/pestpp-ies) for its history-matching section. scripts/setup_env.py downloads it into bin/ with pyemu’s get-pestpp command. If it is missing, the notebook downloads it the same way.

Launch JupyterLab from the repository root:

mamba activate gw-source
jupyter lab tutorials/

Open a notebook and run all cells. Model workspaces are written to tutorials/_work/<tutorial>/ (git-ignored; safe to delete).

The notebooks share one support module, herebedragons.py. It holds the tutorial model builders, an MF6 GWT tracer runner, wrappers around the FloPy model splitter, and the table and map helpers. It is tutorial code, not part of the gwsource API. The notebooks themselves are the narrative. To check that every notebook still runs: pytest -m tutorials (deselected from the default test run).

Order#

Start with 01 and 02: they introduce the outputs and the accounting that every later tutorial uses. After that, 03 to 11 can be read in any order; the list below is the suggested one.

#	Tutorial	What it answers	What you learn	Runtime
01	Quickstart	Where does the water pumped by a well come from?	Build and run a model, read a step, select a well, backward trace; fraction through and the F > 0.1 and F > 0.5 capture zones; volume through and volume from; the split by source (recharge, stream, boundary); the closure check; F against flow-weighted MF6 PRT particles (zone areas, mean \|P − F\|); npz and netCDF export	~5 s
02	Accounting	What exactly is being computed?	Full mixing and the linear system; the Foley & Black (2013) 3×3 worked example rebuilt by hand (net accounting, exact match); gross vs net accounting on a weak-sink well (net attributes 6.7× the pumping); residual, closure and duality diagnostics; boundary-control and dominant-sink maps	~5 s
03	Forward mode	Where does the water in each cell, and in the well, come from?	Source fractions of stream water, regional inflow and recharge; the well’s water by source; source fraction as dilution ratio C/C0 and mixing of source concentrations; volume to by sink; a recharge zone (hypothetical fertilised field) reaching the well; forward/backward duality; agreement with MF6 GWT at equilibrium (≤ 4.2e-13)	~5 s
04	Grids	Does it work on my grid?	The same trace code on DIS, a Voronoi DISV grid refined around the well and the stream (built with vorflow, an optional dependency) and DISU; DIS = DISU exactly, DISV close; boundary placement when comparing grids; why the zone edge is wider on the Voronoi grid (squares rotated 45° to the flow and a finer Voronoi grid); which outputs work on which grid type; XT3D flow loops, the loop report and the exact block solver; boundname and zone-array selectors; GeoPackage and netCDF export	~10 s
05	Transient snapshots	How does the capture zone change through the seasons?	Time steps as snapshots and their limits (no travel lag, storage as source/sink); the well’s water by source per step; min/max/mean F, always/sometimes/never, exceedance count, probability of volume through > 1 m³/d; time vs equal weighting; one way to combine the maps into nested zones	~10 s
06	Travel time	How old is the water reaching the well?	Porosity as an explicit input; Eulerian mean travel time and time-of-travel zones; the exact volume identity ΣVF/Q; weight-losing, flow-weighted MF6 PRT particles: overall agreement (largely the volume identity), per-cell differences and why; when to use particles vs Eulerian times	~10 s
07	Networks and multi-model	What if water moves through streams, drains or wells between layers, or across models?	Cells-only vs package-node traces; SFR reaches and MVR transfers as nodes and edges (drain → mover → stream → aquifer → well); MAW borehole across layers; a split simulation matches the single model; required save options and limitations	~3 s
08	Transient tracking	How much of the well water arrived since the simulation started?	Why snapshots are not transport; forward tracking of tagged sources and initial water; backward (adjoint) fate maps and closure; duality; volume modes “gwt” vs “flow”; agreement with MF6 GWT; four transport sub-steps per flow step, and why (early arrival)	~15 s
09	Ensembles	How sure are we about the capture zone?	50 random-K realisations from 100 kriged pilot points; capture probability P(F > 0.1) and P(F > 0.5); mean and quantile maps; spread of the well’s source shares; exact vs histogram quantiles and the memory budget; serial vs parallel runs; parameter spread vs full mixing and face routing; a PEST++ IES history match of the same ensemble (40 heads, run with local workers) and the prior vs posterior capture probability and river share	~1.5 min
10	Unit concentration	Does tracking reproduce a published unit-concentration study?	The Bedekar et al. (2024) example rebuilt in MF6 and the values the paper leaves open; well shares of lateral inflow, recharge, stream and stored water over 25 years; agreement with four MF6 GWT unit-concentration runs (≤ 1.5e-10) and with the paper’s Figure 2 (≤ 1.3 percentage points); the steady long-time limit; transport sub-steps and arrival fronts; backward time-limited contributing areas and where the pumped recharge and stream water entered	~2.5 min
11	Numerical dispersion	How wide is the edge of a capture zone, and how does it change with the grid?	The band where 0.1 < F < 0.9 against cell size (200 to 12.5 m) and grid orientation on two test problems: uniform flow and a well against the analytical capture zone of Javandel and Tsang (1986), and a recharge-fed zone against MF6 PRT on a 12.5 m grid (93,856 particles); DIS and Voronoi DISV grids; the band narrows as about h^0.5 (fitted rates 0.41 to 0.54); F > 0.5 area within 5 % of the reference area for cells of 50 m or less; F > 0.1 as a conservative envelope (30 to 90 % and 8 to 16 % too large); squares turned 45° widen the edge from 121 m to 215 m; face routing; shares by source and mean travel time barely change	~1 min

Thresholds on the fraction through, and the numerical dispersion of full mixing, are measured in tutorial 11.

The tutorial models#

All models share one plan: a 2 km × 2 km unconfined aquifer in 50 m cells with recharge, a constant-head west edge that brings in regional inflow, a stream (SFR) flowing north to south, its reaches traced as source/sink terms of their cells unless a trace reads them as nodes, and a supply well next to it. herebedragons.py builds variants of it:

Builder	Model
`build_quickstart`	steady, one layer (tutorial 01)
`build_layered`	three layers, partially penetrating supply well and a weak-sink well
`build_voronoi`, `build_disu`	DISV Voronoi grid from vorflow, refined around the well and the stream; DISU copy of the DIS model
`build_xt3d`	heterogeneous anisotropic K with XT3D (flow cycles)
`build_transient`	storage, seasonal recharge and pumping
`build_stream_mover`	the quickstart model plus two drains that feed the stream through the water mover; the stream leaks part of that water back to the aquifer near the well
`build_maw`	three layers with a multi-aquifer well that moves water between layers
`build_realisation`	a seeded random log-K realisation (the synthetic truth of tutorial 09)

Tutorial 10 does not use this plan. It rebuilds the published model of Bedekar et al. (2024) with its own builder, build_bedekar in 10_unit_concentration/helpers.py: 101 × 101 cells of 10 m in three layers, lateral inflow on the west edge, a stream on the east edge and a well pumping 700 m³/d from year 0.

Tutorial 11 uses the plan only in its section 5. Its two test problems are built by its own helpers, in 11_numerical_dispersion/helpers.py: problem A, a confined 2 km × 2 km aquifer with uniform regional flow and a well pumping 500 m³/d, whose perimeter cells hold the analytical head, and problem B, one unconfined layer of 1.5 km × 1 km with recharge, a river, a constant-head west edge and a well pumping 800 m³/d. Both are built on DIS and on Voronoi DISV grids (Triangle, bin/triangle) of several cell sizes.