gwsource.tracking#

True transient mixing-cell tracking: tagged water carried through storage.

Snapshots (gwsource.transient) treat each time step as a steady state in which storage release
is a source and storage gain a sink. Here tagged water is carried from step to step instead:
storage is the memory of the system, not a source. Every node (cell or package node) is a fully
mixed reservoir of water volume V; transport is advective only (no dispersion), upstream
weighted and implicit in time, over the flow model's own time steps. Every time step of the
simulation must be in the budget file (OC BUDGET ALL), so the steps join up in time.

Forward scheme (tagged fraction c = C/C0 per node, step n from t_{n-1} to t_n, length dt)
-----------------------------------------------------------------------------------------
Storage terms (STO-SS, STO-SY, the STORAGE record of package nodes and the CSUB storage terms;
TrackingConfig.storage_texts) are removed from the term table. With q_ij >= 0 the intercell
(and package-network) flow i -> j, s_k > 0 the remaining source terms, d_k = |q_k| the remaining
sink terms, O_i = sum_j q_ij + sum_{k in i} d_k, S_i the storage terms of node i summed (positive
= release) and r_i the node's flow residual (net inflow over all terms, storage included):

    V_i^new c_i^n - V_i^old c_i^(n-1)
        = dt [ sum_j q_ji c_j^n + sum_{k in i, k in S} s_k - (O_i + a r_i) c_i^n ]       (1)

S is the tagged source group and a = 1 if TrackingConfig.flow_error_correction (default): the
residual leaves (r > 0) or enters (r < 0) at the node's own fraction, as MF6 GWT's FMI
FLOW_IMBALANCE_CORRECTION does (Langevin et al. 2022, eq. 2-14/2-15). In matrix form
A^n c^n = B^n c^(n-1) + dt s^n with B^n = diag(V^old) and A^n = diag(V^new + dt (O + a r))
- dt Q^T (Q_ij = q_ij). A^n is an M-matrix: its row sums are V^old + dt Src >= 0 (a = 1), so
the solvers of gwsource.solve apply unchanged (topological substitution when the flow graph is
acyclic, block-triangular LU otherwise; no closed-loop isolation is needed, since V > 0 makes
every recirculation loop non-singular).

Volumes. Both modes satisfy V^old - V^new = dt S within a step, which makes (1) consistent: with
a = 1, c = 1 everywhere solves it whenever every source and the initial water are tagged
(row sum of A^n = V^old + dt Src exactly). They differ in how one step's volumes join the next:
- "gwt": V^new = porosity x area x saturated thickness at the end of the step (package
  nodes: the VOLUME of their STORAGE record) and V^old = V^new + dt S, exactly as MF6 GWT's MST
  (mst_fc_sto) and the advanced transport packages (apt_get_volumes) compute them. The volume
  at the start of step n+1 is then not the volume at the end of step n:
  V^old,n+1 - V^new,n = (Sy - porosity) x area x (fall of the water table in step n+1), plus
  the specific-storage part, whenever porosity differs from the specific yield. At each step
  boundary tagged water c^n x (V^old,n+1 - V^new,n) appears or disappears; it is reported as
  the volume adjustment of the balance. One exception copies MF6 SFT (sft_get_volumes): an
  SFR reach whose storage flow is exactly 0 (SFR without kinematic-wave storage, or a steady
  period) takes the previous step's V^new as V^old (0 at the first step), so its volume is
  carried as in flow mode but its within-step water balance is not closed (reported in
  water_balance_error). This mode reproduces GWT unit-concentration results.
- "flow" (default): the volume is carried by the storage flows themselves. The start volume is
  porosity x saturated volume (package volume) at the end of step 1 plus dt S, then
  V^old,n+1 = V^new,n and V^new = V^old - dt S. The water balance of every node is exact by
  construction and tagged water is conserved across step boundaries. Porosity then only sets
  the initial volume (the mobile water that mixes); the change of volume is the flow model's
  storage change (specific yield and specific storage), so the volume can drift away from
  porosity x saturated volume. A node whose volume would become negative (porosity much smaller
  than the specific yield) stops the run with an error.

Tagged initial water: TrackingConfig.initial_water = {name: "all" or a node mask/indices}
adds one forward column per entry with c^0 = 1 on those nodes and no source: the fraction of
pre-simulation water.

Forward balance per group and step (volumes over the step, all exact to round-off):
    storage_end - storage_start = inflow - outflow - flow_error
with storage_start = sum V^old c^(n-1), storage_end = sum V^new c^n, inflow = dt sum s,
outflow = dt sum_i (sink total of i) c_i (also split by package and by sink group) and
flow_error = dt a sum r_i c_i; volume_adjustment = storage_start^n - storage_end^(n-1) (0 in flow
mode; at step 1 it is 0 and storage_start is the tagged initial water).

Backward (adjoint) scheme
-------------------------
The tagged outflow through a destination group D over the whole simulation is
J = sum_n dt^n d^n . c^n (d^n the destination sink flow per node). Its discrete adjoint runs
backward in time from the end of the simulation, with w^(N+1) = 0:

    (A^n)^T w^n = dt^n d^n + B^(n+1) w^(n+1)                                           (2)

w^n_i is the fraction of the water in node i at the start of step n (time t_(n-1)), or entering
it during step n, that leaves through D before the simulation ends (an exact reading of the
implicit scheme, in which the water present at the start of a step and the water arriving
during it are mixed alike). Multiplying (1) by w^n and (2) by c^n and summing over the steps
gives the duality identity, exact for any volumes and any a:

    J = sum_n dt^n s^n . w^n + sum_i V^old,1_i c^0_i w^1_i                               (3)

(tagged source inflow times the adjoint, plus the tagged initial water times the adjoint at the
start). With every source tagged and c^0 = 1 (and a = 1) the right-hand side splits J by
source term and step: volume_from (per node and step), volume_from_by_package, and the initial
water's share initial_volume_from; their total is the destination outflow (closure). With a = 0
the closure misses sum_n dt r^n . w^n (reported as residual_term).

The backward run also solves (2) for "any sink" (d = every remaining sink) and for the water
still in storage at the end (no d; terminal right-hand side V^new,N). For every node and step
fraction_through[D] + fraction_to_other_sinks + fraction_in_storage_at_end + fraction_unaccounted
= 1, where the unaccounted part is the water removed by the flow-error correction and, in gwt
mode, by the volume adjustments (zero up to the flow residual in flow mode).

Sub-steps (TrackingConfig.substeps = m, default 1): every flow step is split into m equal
transport sub-steps of length dt/m. The flows are constant within a flow step, so this changes
only the time discretisation of the transport, not the flows; the volumes are interpolated
linearly between V^old and V^new (exactly what the constant storage flow gives within the step).
Implicit upstream stepping adds numerical dispersion in time of order dt; sub-steps reduce it
without touching the flow model. Per-step outputs are reported at flow-step boundaries: c at the
end of the step, w at its start, and the step's flow-weighted totals (outflow, volume from,
residual terms) summed over its sub-steps. The duality and closure identities hold unchanged
(the adjoint is that of the sub-stepped scheme). m = 1 reproduces the unsplit scheme exactly.

Limit: with flows held constant for long enough, c^n tends to the steady source_fraction of
gwsource.trace.forward and w^1 to the steady fraction_through of gwsource.trace.backward
(storage terms are zero then, so the snapshot and the tracked systems have the same flows).
class gwsource.tracking.TrackingConfig(porosity, volume_mode='flow', flow_error_correction=True, initial_water=None, solver='auto', storage_texts=('STO-SS', 'STO-SY', 'STORAGE', 'CSUB-CGELASTIC', 'CSUB-ELASTIC', 'CSUB-INELASTIC', 'CSUB-WATERCOMP'), keep_series=True, substeps=1)[source]#

Bases: object

Options of transient tracking (track_forward, track_backward).

porosity: scalar, per-node array or {model: per-cell array} (required).
volume_mode: "flow" (default: volumes carried by the storage flows; exact water balance) or
    "gwt" (volumes from porosity x saturated volume, as MF6 GWT; use it to compare with GWT);
    see the module docstring.
flow_error_correction: the node residual leaves or enters at the node's fraction (as GWT's
    FLOW_IMBALANCE_CORRECTION); makes the closure exact.
initial_water: {name: "all" or a node mask / node indices}: forward columns of tagged
    initial water; backward, the initial water's share of each destination's outflow.
solver: "auto", "topological", "block" or "lu" (gwsource.solve).
storage_texts: budget texts of the storage terms (taken out of the term table).
keep_series: keep every step's per-node outputs (else only the last step's state).
substeps: number of equal transport sub-steps per flow step (1 = the flow model's steps).
Parameters:
class gwsource.tracking.StepFlows(edges, tracked, source, sink, outflow, inflow, storage, residual)[source]#

Bases: object

The flows of one step as transient tracking uses them (storage terms taken out).
Parameters:
class gwsource.tracking.StepVolumes(old, new)[source]#

Bases: object

Water volume of every node at the start (old) and end (new) of a step.
Parameters:
class gwsource.tracking.ForwardTracking(config, groups, steps, totim, delt, source_fraction, initial, balance, outflow_by_package, outflow_to, water_balance_error, volume_jump, template, solver_used=<factory>, notes=<factory>)[source]#

Bases: object

Forward transient tracking: tagged fractions over time and tagged-water balances.

source_fraction[g]: (n_steps, n) fraction c at the end of each step (only the last step if
keep_series is False); initial[g]: c at the start. balance[g][key]: per-step volumes (keys in
BALANCE_KEYS). outflow_by_package[g][package] and outflow_to[sink group][g]: tagged volume
leaving per step. water_balance_error: per step, max |V^new - V^old - dt (net inflow)| over
the nodes; volume_jump: max |V^old,n - V^new,n-1| (0 in flow mode).
Parameters:
cumulative_outflow(group)[source]#
Return the tagged outflow of a group by package, cumulated over the steps.
Parameters:

group (str)

Return type:

dict[str, ndarray]

class gwsource.tracking.BackwardTracking(config, groups, steps, totim, delt, fraction_through, fraction_any_sink, fraction_in_storage_at_end, volume_from, volume_from_by_package, volume_from_group, initial_volume_from, initial_volume_from_group, destination_outflow, residual_term, active, template, solver_used=<factory>, notes=<factory>)[source]#

Bases: object

Backward (adjoint) transient tracking from destination groups.

Row n of every (n_steps, n) array is step n: the water in the node at the start of step n
(time_start[n]) or entering it during step n. fraction_through[g]: fraction leaving
through g before the end; fraction_any_sink: through any sink; fraction_in_storage_at_end:
still in the model at the end. volume_from[g]: water entering through the node's source
terms during the step that leaves through g (by package: volume_from_by_package[g][p], per
step totals); initial_volume_from[g]: initial water of each node that leaves through g;
volume_from_group[g][source group] and initial_volume_from_group[g][name]: the same for the
source groups and initial_water entries given. destination_outflow[g]: water leaving
through g in each step. residual_term[g]: per step dt r . w (closure term when the flow-error
correction is off). active: nodes holding or passing water in the step (the fractions of the
other nodes are 0 and meaningless).
Parameters:
property time_start: ndarray#
Simulation time at the start of each step.
fraction_to_other_sinks(group)[source]#
Return the fraction leaving through a sink outside the group before the end.
Parameters:

group (str)

Return type:

ndarray

property fraction_unaccounted: ndarray#
Return 1 - any sink - storage at end (flow-error correction, volume adjustments).

NaN on nodes that are not active in the step.
gwsource.tracking.storage_mask(step, texts=('STO-SS', 'STO-SY', 'STORAGE', 'CSUB-CGELASTIC', 'CSUB-ELASTIC', 'CSUB-INELASTIC', 'CSUB-WATERCOMP'))[source]#
Return the term entries that are storage terms (by budget text).
Parameters:
Return type:

ndarray

gwsource.tracking.step_flows(step, storage_texts=('STO-SS', 'STO-SY', 'STORAGE', 'CSUB-CGELASTIC', 'CSUB-ELASTIC', 'CSUB-INELASTIC', 'CSUB-WATERCOMP'))[source]#
Return the edges, non-storage source and sink totals, storage and residual per node.
Parameters:
Return type:

StepFlows

gwsource.tracking.sfr_nodes(step)[source]#
Return a node mask of the SFR reaches of a network step.
Parameters:

step (FlowStep)

Return type:

ndarray

gwsource.tracking.gwt_volumes(step, flows, porosity, delt, previous=None)[source]#
Return V^new = porosity x saturated volume (or package volume), V^old = V^new + dt S.

SFR reaches whose storage flow is exactly 0 (SFR without kinematic-wave storage, or a
steady-state period) take the previous step's V^new as V^old instead (0 before the first
step), as MF6 SFT does (sft_get_volumes); previous is that step's V^new.
Parameters:
Return type:

StepVolumes

gwsource.tracking.carried_volumes(old, flows, delt)[source]#
Return flow-mode volumes: V^new = V^old - dt S (asserting that no volume is negative).
Parameters:
Return type:

StepVolumes

gwsource.tracking.substep_volumes(volumes, index, n_sub)[source]#
Return the volumes of sub-step index (0-based) of n_sub, linear between V^old and V^new.

The first sub-step starts at V^old and the last ends at V^new exactly, so n_sub = 1 returns
the step's own volumes.
Parameters:
Return type:

StepVolumes

gwsource.tracking.empty_rows(flows, volumes)[source]#
Return the nodes with no water and no flow (inactive or dry): identity rows, value 0.
Parameters:
Return type:

ndarray

gwsource.tracking.step_matrix(flows, volumes, delt, correction)[source]#
Return A = diag(V^new + dt (O + a r)) - dt Q^T of the forward step equation (1).
Parameters:
Return type:

csr_matrix

gwsource.tracking.group_flows(step, groups, tracked, sign)[source]#
Return (nodes, groups) totals of |q| over each group's tracked terms of the given sign.
Parameters:
Return type:

ndarray

gwsource.tracking.by_package(step, term_values, tracked)[source]#
Return {package name: sum of term_values} over the packages with tracked term entries.
Parameters:
Return type:

dict[str, float]

gwsource.tracking.check_time(totim, delt, previous, kstpkper)[source]#
Assert that a step starts where the previous one ended (every step must be saved).
Parameters:
Return type:

None

gwsource.tracking.initial_columns(step, initial_water)[source]#
Return {name: c^0} of the tagged-initial-water entries ('all' or node mask/indices).
Parameters:
Return type:

dict[str, ndarray]

gwsource.tracking.step_volumes(step, flows, config, delt, previous)[source]#
Return the step's volumes; previous is the previous step's V^new (None at the start).

gwt mode uses previous only for the SFR rule of gwt_volumes(); flow mode carries it, and
starts from porosity x saturated volume + dt S of the first step.
Parameters:
Return type:

StepVolumes

gwsource.tracking.water_balance_error(flows, volumes, delt, correction)[source]#
Return max |V^new - V^old - dt (I - O - a r)| over the nodes (non-storage flows).
Parameters:
Return type:

float

gwsource.tracking.track_forward(reader, sources, config, sinks=None)[source]#
Carry tagged water from source groups (and tagged initial water) through every step.

reader: StepReader (gwsource.open_model) or NetworkReader (gwsource.open_simulation).
sources, sinks: dict {group: select() kwargs} or a function step -> {group: Selection},
re-selected on every step; sinks only report the tagged outflow through each sink group.
Parameters:
Return type:

ForwardTracking

gwsource.tracking.flow_mode_end_volume(reader, steps, config)[source]#
Return V^new of the last step in flow mode (a forward pass over the storage terms).
Parameters:
Return type:

ndarray

gwsource.tracking.track_backward(reader, destinations, config, sources=None)[source]#
Solve the adjoint backward in time from the end of the simulation for each destination.

sources (optional): source groups whose share of each destination's outflow is reported
per step (volume_from_group); config.initial_water entries give initial_volume_from_group.
Parameters:
Return type:

BackwardTracking

gwsource.tracking.forward_balance_error(result, group)[source]#
Return max |error| of a group's tagged-water balance relative to its largest term.
Parameters:
Return type:

float

gwsource.tracking.backward_closure(result, group)[source]#
Return sum(volume from) + initial volume from = destination outflow over the simulation.

Exact with the flow-error correction (residual_term 0); without it the residual term is
sum_n dt r^n . w^n.
Parameters:
Return type:

Check

gwsource.tracking.tracking_duality(fwd, source, bwd, destination)[source]#
Return the tagged water from source to destination, forward vs adjoint (identity (3)).

source is a forward group (a source group or an initial_water entry); the forward run must
report the destination as a sink group, and the backward run the source as a source group
or initial_water entry. Check.lhs is the adjoint side, Check.rhs the forward side.
Parameters:
Return type:

Check