gwsource.prt#
Thin helper to check Eulerian travel times against forward particle tracking with MF6 PRT.
What it does
------------
1. Contributing cells: cells whose backward fraction F to the destination group exceeds a
threshold (under full mixing F = 0 means no flow path, so no particle could arrive either).
2. Release points: one set per source term (positive external flow, optionally restricted to a
source Selection) in a contributing cell, placed where PRT's velocity field brings that
term's water in. Every point carries an equal-area (or equal-length) share of q_term, so
sums over particles are flux-weighted.
Plan points: the cell's plan polygon (DIS or DISV) is split into a fan of triangles from its
centroid, each triangle into k x k equal sub-triangles, and a point is placed at each
sub-triangle's centroid (weight = sub-triangle area / cell area).
Distributed terms (IFLOWFACE = 0, or no IFLOWFACE auxiliary): each plan point is released
at z_levels local z = (i + 0.5) / z_levels (fractions of the saturated thickness), so the
points fill the cell volume. That is how MF6 represents a boundary source by default:
recharge, river or drain leakage and edge CHD/GHB inflow are distributed over the cell, and
PRT's velocity field has the matching uniform divergence. Releasing everything at one
height is not consistent with that field. In particular release_z=1 (the water table)
traps the particles of the top layer at the water table: the top face carries no flow, so
the vertical velocity there is zero and they never descend (in a 3-layer test model all of
them ended in layer 1 and none reached a well in layer 3).
Face terms: a package can assign its flow to a cell face with the auxiliary variable
IFLOWFACE (AUXILIARY IFLOWFACE in the GWF package; the value is saved as a column of its
budget record, which is where PRT reads it, see gwsource.mf6io). Values (MF6 PRT source:
prt-fmi.f90, MethodDis, MethodDisv): 0 = distributed; -1 = top; -2 = bottom; i = 1..nvert =
lateral face i, the edge from vertex i to vertex i + 1 of the cell polygon (DISV: the
CELL2D vertex order; DIS: 1 west, 2 north, 3 east, 4 south). PRT then treats the term as a
flux through that face (a boundary face flow, no divergence inside the cell), and the
consistent release is on that face: top or bottom face terms at the plan points with local
z = 1 or 0 (one height); a lateral face term at points along its edge (as many as the
cell has plan points, at the midpoints of equal pieces) x the z_levels heights, a
midpoint rule for a flux spread uniformly over the saturated face. Recharge on the top
face is the uppermost water of the column and descends only as recharge is added on top
of it downstream (in the 3-layer test model with IFLOWFACE = -1 no recharge particle
reaches the bottom-layer wells; volume seeding there credits them with 6.9 m3/d of it).
3. PRT: a separate simulation (own folder) with PRT DIS/DISV copied from the GWF input, MIP
porosity, PRP with EXTEND_TRACKING (the steady period is extended until particles
terminate), FMI reading the GWF head and budget files, and EMS. FloPy writes the input; the
MF6 executable is passed explicitly. PRT needs DATA-SPDIS (NPF SAVE_SPECIFIC_DISCHARGE).
The GWF model must have a single steady stress period: PRT moves the particles in the flows
of every period, but their weights and the destination shares come from one FlowStep.
Particles run until they terminate (strong sink, boundary face); there is no stop zone, and
track_particles() asserts that every particle terminated.
4. Tracks: release, cell-exit and termination records (CSV track output). icell is the 1-based
user node; the records give each particle's sequence of cells and the time it entered them.
Weak sinks and the destination (stop_at_weak_sink is deliberately left off). PRT does not
report which boundary package removes water, and with STOP_AT_WEAK_SINK a particle would stop
in every weak sink cell (sinks plus outflow faces), also those of other boundaries (e.g. the
gaining river cells between recharge and a well), which truncates paths. Instead each particle
carries the full-mixing accounting of the gross method along its advective path: in every cell
j it visits (the release cell included, each cell once), its remaining weight W gives
d_j / O_j x W to the destination (d_j = the group's sink flow in j, O_j = the gross outflow of
j) and keeps (1 - Snk_j / O_j) x W (Snk_j = all sink terms of j, the destination's included).
The captured weight is the sum over the path; its arrival time is the entry time into each j.
Nothing is dropped: water that passes a destination cell, or a well that takes only part of
its cell's outflow, is followed on to the destination cells further down.
With thin_weak_sinks=False the particles lose only the destination's share (1 - d_j / O_j), as
in a plain particle run; that over-counts the water reaching the destination, because the
water that other sinks remove on the way is counted too.
Consistency. If the release points resolve the sources and the cell survival splits the
weight among the outflow faces as the flow does, the captured flow equals the destination
inflow (the duality of the gross method, now along advective paths). Observed:
- 3-layer test model (wells in the bottom layer, a weak well upstream of a strong one):
402.6 of 405 m3/d with k = 3 and 2 levels, 405.3 with k = 4 and 16 levels. Both parts of
the method matter: releasing at the water table and stopping particles in the first
destination cell, which drops the weight a weak well passes on, captures 266 m3/d; volume
seeding with that stop rule about 306 m3/d.
With recharge on the top face (IFLOWFACE = -1, same GWF flows): 405.8 released on the top
face, 412.7 seeded over the volume. With the CHD edges on their lateral faces instead (1
west, 3 east; not in the test suite): 398 to 412 m3/d for k = 3..6 and 2..8 levels, the same
scatter as volume seeding of the CHD (395 to 403): each edge cell brings ~126 m3/d, so a few
particles on the capture-zone edge move the sum.
- Neversink (4 layers, 1.68M cells, eight wells, k = 1, 4 levels): 3,736 of 3,671 m3/d,
every well within 11 % of its rate. Stopping particles in the first destination cell gives
1,940 with release at the water table and 3,005 at mid-cell. Following the weight on, a
single mid-cell height gives 3,722, so the z levels matter most where the alternative
would be the water table. Over five seeding layouts the total scattered
between 0.0 and +1.8 % without a trend (quadrature, not the thinning factor); the per-well
errors shrank with the number of z levels (RMS 18 % with one, 6 % with four).
Which sources supply the water differs between advection and mixing (Neversink: 51 % stream
leakage along the paths, 43 % under full mixing), and so do per-cell capture fractions; that
is the method difference, not an error.
What is measured. Water is counted as arriving when it ENTERS a destination cell (at once for
a particle released in it). The arrival time is therefore comparable with
TravelTimeResult.mean_time_to_exit_cell (the time until the water enters the cell through
whose sink it leaves); mean_time also includes the residence in that cell.
PRT supports DIS and DISV only; for DISU the Eulerian travel time is the only option.
Network steps (gwsource.network: several models, or package nodes) are not supported: the
destination shares and the particle weights assume that the destination sinks and every node
of the graph are cells of the one GWF model PRT tracks in, and particles do not follow package
networks (an SFR reach, a MAW well bore) or cross into other models. track_particles() raises
NotImplementedError for them.
- class gwsource.prt.PrtConfig(porosity, fraction_threshold=1e-06, subdivisions=2, z_levels=4, release_z=None, thin_weak_sinks=True, exit_solve_tolerance=1e-05, stoptraveltime=None)[source]#
Bases:
objectOptions of a PRT run. porosity: scalar or per-cell (required; use the same value as for the Eulerian times). fraction_threshold: seed only cells with backward fraction above this. subdivisions: k; each fan triangle of a cell gets k * k points (4 k^2 points per DIS cell). z_levels: n; each plan point of a distributed term is released at the n local z = (i + 0.5) / n, i = 0..n-1, so the points fill the saturated cell volume (see module docstring); lateral face terms use the same heights. A one-layer grid has no vertical flow and gets one level. Top and bottom face terms (IFLOWFACE -1, -2) are always released on their face (local z 1 or 0). release_z: optional single local z instead of the levels (0 = cell bottom, 1 = water table / cell top); release_z=1 traps particles of distributed terms at the water table of the top layer. thin_weak_sinks: multiply particle weights by the full-mixing survival of every cell they pass, so other sinks remove their share (see module docstring). exit_solve_tolerance: PRP EXIT_SOLVE_TOLERANCE. stoptraveltime: optional PRP STOPTRAVELTIME.
- Parameters:
- class gwsource.prt.ReleasePoints(cell, x, y, z, weight, package)[source]#
Bases:
objectParticle release points: 0-based cell (user node), x, y, local z, flow weight, package.
- class gwsource.prt.PrtResult(release, end_cell, end_status, arrival_time, captured_weight, mean_arrival, captured_flow, released_flow, cell_mean_arrival, cell_capture_fraction, package_mean_arrival, ws)[source]#
Bases:
objectPer-particle end points and flux-weighted arrival statistics at the destination. arrival_time: capture-weighted mean time from release to entering the destination cells, per particle (NaN if not captured). captured_weight: particle weight x the share of it captured along its path (see module docstring; 0 if it passes no destination cell). mean_arrival: flux-weighted mean arrival time of the captured particles. captured_flow: sum of captured weights (compare with the group's inflow from these sources). cell_mean_arrival, cell_capture_fraction: per release cell (NaN where nothing was released); package_mean_arrival: flux-weighted mean arrival time by source package.
- class gwsource.prt.PrtGrid(grid_type, dis_kwargs, polygons, ncpl)[source]#
Bases:
objectGWF discretisation copied for PRT: grid type, FloPy keyword arguments, cell polygons.
- gwsource.prt.load_prt_grid(ws, model, dis_type)[source]#
Load the GWF DIS or DISV input with FloPy and return what PRT needs from it.
- gwsource.prt.fan_points(polygon, k)[source]#
Return (points (n, 2), area weights summing to 1) spread uniformly over a convex polygon. The polygon is split into triangles from its centroid; each triangle into k * k equal sub-triangles, with one point at each sub-triangle's centroid.
- gwsource.prt.edge_points(polygon, face, n)[source]#
Return (points (n, 2), equal weights summing to 1) along lateral face `face` (1-based). Face i is the polygon edge from vertex i to vertex i + 1 (MF6 PRT numbering); the points are the midpoints of n equal pieces of it.
- gwsource.prt.term_points(polygon, iflowface, k, local_z)[source]#
Return (plan points (n, 2), local z, weight shares summing to 1) of one source term. IFLOWFACE 0: the cell volume (fan points x local_z); top (-1) or bottom (-2) face: the fan points at local z 1 or 0; lateral face i: as many points along edge i as the fan has, x local_z.
- gwsource.prt.release_points(step, grid, contributing, source, subdivisions, local_z)[source]#
Return flux-weighted release points for the source terms in contributing cells. Each term is released over its cell volume or over its IFLOWFACE face (term_points).
- gwsource.prt.assert_single_steady_period(files)[source]#
Assert that the GWF model has one stress period and that it is steady. PRT tracks through the whole GWF time history, while the particle weights and destination shares come from one FlowStep; they agree only if that step is the one steady flow field.
- Parameters:
files (ModelFiles)
- Return type:
None
- gwsource.prt.write_prt(prt_ws, gwf_ws, grid, release, config, exe, model=None, budget_file=None, stop_at_weak_sink=False, track_release=True, track_exit=True, clear_ws=False)[source]#
Write a PRT simulation tracking the release points until they leave the model. budget_file: GWF budget file for FMI instead of the model's own (e.g. a copy with edited IFLOWFACE columns); it must hold DATA-SPDIS. stop_at_weak_sink: PRP STOP_AT_WEAK_SINK (off by default, see the module docstring). track_release, track_exit: also write release and cell-exit records (cell_visits() needs both); termination records are always written. clear_ws: delete prt_ws first, so no output of an earlier run is left in it.
- gwsource.prt.run_prt(sim)[source]#
Run a written PRT simulation and return the CSV track file.
- Parameters:
sim (MFSimulation)
- Return type:
- gwsource.prt.read_tracks(track_csv)[source]#
Return particle, cell, status, reason and times of a PRT CSV track file (record order).
- gwsource.prt.endpoints(tracks, n_particles)[source]#
Return, per particle (release point order), its end cell (0-based), status and time.
- gwsource.prt.assert_terminated(ends, release)[source]#
Assert that every particle has a termination record (end status >= 0).
- Parameters:
release (ReleasePoints)
- Return type:
None
- gwsource.prt.cell_visits(tracks)[source]#
Return the cells each particle passes, in order, with the time it entered them. A particle is released in its first cell (entry time 0); every cell-exit record ends one visit, and the next record's cell is the next visit, entered at the exit time. Only the first visit of a cell counts (a particle that re-enters a cell, e.g. by grazing a face, is not counted twice).
- gwsource.prt.capture_along_paths(visits, n_particles, share, survival)[source]#
Return, per particle, the captured share of its weight and the mean capture time. In each visited cell j the particle's remaining weight W gives share_j x W to the destination and keeps survival_j x W; the capture time is the entry time into j.
Return d_k / O_k per cell: the share of a cell's outflow leaving through the group.
- gwsource.prt.weak_sink_survival(step)[source]#
Return 1 - Snk / O per cell (1 where O = 0): the share of the outflow through faces.
- gwsource.prt.arrival_statistics(step, release, ends, captured_share, arrival, ws)[source]#
Return per-particle captured weights and flux-weighted arrival statistics. captured_share and arrival are the per-particle results of capture_along_paths().
- gwsource.prt.check_cells_only(step, destination)[source]#
Raise NotImplementedError for network steps (package-node destinations, several models).
- gwsource.prt.track_particles(gwf_ws, step, result, group, config, exe, prt_ws, source=None, model=None)[source]#
Seed, run and evaluate forward PRT particles for one backward destination group. step and result are the GWF step and its backward trace (gross). Particles are released from the source terms (all, or those in source) of cells with F > fraction_threshold. The GWF model must have a single steady stress period, and every particle must terminate.