Coverage for src/gwtransport/fronttracking/output.py: 89%

1"""

2Concentration Extraction from Front Tracking Solutions.

4This module computes outlet concentrations from front tracking wave solutions

5using exact analytical integration. All calculations maintain machine precision

6with no numerical dispersion.

8Functions

9---------

10concentration_at_point(v, t, waves, sorption)

11 Compute concentration at any point (v, t) in domain

12compute_breakthrough_curve(t_array, v_outlet, waves, sorption)

13 Compute concentration at outlet over time array

14compute_bin_averaged_concentration_exact(t_edges, v_outlet, waves, sorption)

15 Compute bin-averaged concentrations using exact analytical integration

16compute_domain_mass(t, v_outlet, waves, sorption)

17 Compute total mass in domain [0, v_outlet] at time t using exact analytical integration

18compute_cumulative_inlet_mass(t, cin, flow, tedges_days)

19 Compute cumulative mass entering domain from t=0 to t

20compute_cumulative_outlet_mass(t, v_outlet, waves, sorption, flow, tedges_days)

21 Compute cumulative mass exiting domain from t=0 to t

22find_last_rarefaction_start_time(v_outlet, waves)

23 Find time when last rarefaction head reaches outlet

24integrate_rarefaction_total_mass(raref, v_outlet, t_start, sorption, flow)

25 Compute total mass exiting through rarefaction to infinity

26compute_total_outlet_mass(v_outlet, waves, sorption, flow, tedges_days)

27 Compute total integrated outlet mass until all mass has exited

29This file is part of gwtransport which is released under AGPL-3.0 license.

30See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.

31"""

33from operator import itemgetter

35import mpmath as mp

36import numpy as np

37import numpy.typing as npt

38from scipy.special import beta as beta_func

39from scipy.special import betainc

41from gwtransport.fronttracking.events import find_outlet_crossing

42from gwtransport.fronttracking.math import ConstantRetardation, FreundlichSorption

43from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave, Wave

45# Numerical tolerance constants

46EPSILON_VELOCITY = 1e-15 # Tolerance for checking if velocity is effectively zero

47EPSILON_BETA = 1e-15 # Tolerance for checking if beta is effectively zero (linear case)

48EPSILON_TIME = 1e-15 # Tolerance for negligible time segments

49EPSILON_TIME_MATCH = 1e-6 # Tolerance for matching arrival times (for rarefaction identification)

50EPSILON_CONCENTRATION = 1e-10 # Tolerance for checking if concentration is effectively zero

53def concentration_at_point(

54 v: float,

55 t: float,

56 waves: list[Wave],

57 sorption: FreundlichSorption | ConstantRetardation, # noqa: ARG001

58) -> float:

59 """

60 Compute concentration at point (v, t) with exact analytical value.

62 Searches through all waves to find which wave controls the concentration

63 at the given point in space-time. Uses exact analytical solutions for

64 characteristics, shocks, and rarefaction fans.

66 Parameters

67 ----------

68 v : float

69 Position [m³].

70 t : float

71 Time [days].

72 waves : list of Wave

73 All waves in the simulation (active and inactive).

74 sorption : FreundlichSorption or ConstantRetardation

75 Sorption model.

77 Returns

78 -------

79 concentration : float

80 Concentration at point (v, t) [mass/volume].

82 Notes

83 -----

84 **Wave Priority**:

85 The algorithm checks waves in this order:

86 1. Rarefaction waves (if point is inside rarefaction fan)

87 2. Shocks (discontinuities)

88 3. Rarefaction tails (if rarefaction tail has passed point)

89 4. Characteristics (smooth regions)

91 If no active wave controls the point, returns 0.0 (initial condition).

93 **Rarefaction Tail Behavior**:

94 After a rarefaction tail passes a query point, that point maintains the

95 tail concentration as a plateau. This ensures proper concentration behavior

96 after rarefaction waves pass through.

98 **Machine Precision**:

99 All position and concentration calculations use exact analytical formulas.

100 Numerical tolerance is only used for equality checks (v == v_shock).

101

102 Examples

103 --------

104 ::

105

106 sorption = FreundlichSorption(

107 k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3

108 )

109 # After running simulation with waves...

110 c = concentration_at_point(v=250.0, t=5.0, waves=all_waves, sorption=sorption)

111 c >= 0.0

112 """

113 # Check rarefactions first (they have spatial extent and override other waves)

114 for wave in waves:

115 if isinstance(wave, RarefactionWave) and wave.is_active:

116 c = wave.concentration_at_point(v, t)

117 if c is not None:

118 return c

119

120 # Track the most recent wave to affect position v

121 # We need to compare crossing times of shocks and rarefaction tails

122 latest_wave_time = -np.inf

123 latest_wave_c = None

124

125 # Check shocks - track which shocks control this position

126 for wave in waves:

127 if isinstance(wave, ShockWave) and wave.is_active:

128 v_shock = wave.position_at_time(t)

129 if v_shock is not None:

130 # Tolerance for exact shock position

131 tol = 1e-15

132

133 if abs(v - v_shock) < tol:

134 # Exactly at shock position

135 return 0.5 * (wave.c_left + wave.c_right)

136

137 # Determine if shock has crossed position v and when

138 if wave.velocity is not None and abs(wave.velocity) > EPSILON_VELOCITY:

139 t_cross = wave.t_start + (v - wave.v_start) / wave.velocity

140

141 if t_cross <= t:

142 # Shock has crossed position v by time t

143 # After crossing, point sees c_left (concentration behind shock)

144 if t_cross > latest_wave_time:

145 latest_wave_time = t_cross

146 latest_wave_c = wave.c_left

147 elif v > v_shock and wave.t_start > latest_wave_time:

148 # Point is ahead of shock (shock hasn't reached it yet)

149 # Check if this is the closest shock ahead of us

150 # In this case, we see c_right unless overridden by another wave

151 # We track this with a negative time (shock formation time) to indicate

152 # it's a "passive" state (not actively changed)

153 latest_wave_time = wave.t_start

154 latest_wave_c = wave.c_right

155

156 # Check rarefaction tails - they can override previous waves

157 for wave in waves:

158 if isinstance(wave, RarefactionWave) and wave.is_active:

159 v_tail = wave.tail_position_at_time(t)

160 if v_tail is not None and v_tail > v + 1e-15:

161 # Rarefaction tail has passed position v

162 # Find when it passed: v_start + tail_velocity * (t_pass - t_start) = v

163 tail_vel = wave.tail_velocity()

164 if tail_vel > EPSILON_VELOCITY:

165 t_pass = wave.t_start + (v - wave.v_start) / tail_vel

166 if t_pass <= t and t_pass > latest_wave_time:

167 latest_wave_time = t_pass

168 latest_wave_c = wave.c_tail

169

170 if latest_wave_c is not None:

171 return latest_wave_c

172

173 # Check characteristics

174 # Find the most recent characteristic that has reached position v by time t

175 latest_c = 0.0

176 latest_time = -np.inf

177

178 for wave in waves:

179 if isinstance(wave, CharacteristicWave) and wave.is_active:

180 # Check if this characteristic has reached position v by time t

181 v_char_at_t = wave.position_at_time(t)

182

183 if v_char_at_t is not None and v_char_at_t >= v - 1e-15:

184 # This characteristic has passed through v

185 # Find when it passed: v_start + vel*(t_pass - t_start) = v

186 vel = wave.velocity()

187

188 if vel > EPSILON_VELOCITY:

189 t_pass = wave.t_start + (v - wave.v_start) / vel

190

191 if t_pass <= t and t_pass > latest_time:

192 latest_time = t_pass

193 latest_c = wave.concentration

194

195 return latest_c

196

197

198def compute_breakthrough_curve(

199 t_array: npt.NDArray[np.floating],

200 v_outlet: float,

201 waves: list[Wave],

202 sorption: FreundlichSorption | ConstantRetardation,

203) -> npt.NDArray[np.floating]:

204 """

205 Compute concentration at outlet over time array.

206

207 This is the breakthrough curve: C(v_outlet, t) for all t in t_array.

208

209 Parameters

210 ----------

211 t_array : array-like

212 Time points [days]. Must be sorted in ascending order.

213 v_outlet : float

214 Outlet position [m³].

215 waves : list of Wave

216 All waves in the simulation.

217 sorption : FreundlichSorption or ConstantRetardation

218 Sorption model.

219

220 Returns

221 -------

222 c_out : numpy.ndarray

223 Array of concentrations matching t_array [mass/volume].

224

225 Examples

226 --------

227 ::

228

229 t_array = np.linspace(0, 100, 1000)

230 c_out = compute_breakthrough_curve(

231 t_array, v_outlet=500.0, waves=all_waves, sorption=sorption

232 )

233 len(c_out) == len(t_array)

234

235 See Also

236 --------

237 concentration_at_point : Point-wise concentration

238 compute_bin_averaged_concentration_exact : Bin-averaged concentrations

239 """

240 t_array = np.asarray(t_array, dtype=float)

241 c_out = np.zeros(len(t_array))

242

243 for i, t in enumerate(t_array):

244 c_out[i] = concentration_at_point(v_outlet, t, waves, sorption)

245

246 return c_out

247

248

249def identify_outlet_segments(

250 t_start: float, t_end: float, v_outlet: float, waves: list[Wave], sorption: FreundlichSorption | ConstantRetardation

251) -> list[dict]:

252 """

253 Identify which waves control outlet concentration in time interval [t_start, t_end].

254

255 Finds all wave crossing events at the outlet and constructs segments where

256 concentration is constant or varying (rarefaction).

257

258 Parameters

259 ----------

260 t_start : float

261 Start of time interval [days].

262 t_end : float

263 End of time interval [days].

264 v_outlet : float

265 Outlet position [m³].

266 waves : list of Wave

267 All waves in the simulation.

268 sorption : FreundlichSorption or ConstantRetardation

269 Sorption model.

270

271 Returns

272 -------

273 segments : list of dict

274 List of segment dictionaries, each containing:

275

276 - 't_start' : float

277 Segment start time

278 - 't_end' : float

279 Segment end time

280 - 'type' : str

281 'constant' or 'rarefaction'

282 - 'concentration' : float

283 For constant segments

284 - 'wave' : Wave

285 For rarefaction segments

286 - 'c_start' : float

287 Concentration at segment start

288 - 'c_end' : float

289 Concentration at segment end

290

291 Notes

292 -----

293 Segments are constructed by:

294 1. Finding all wave crossing events at outlet in [t_start, t_end]

295 2. Sorting events chronologically

296 3. Creating constant-concentration segments between events

297 4. Handling rarefaction fans with time-varying concentration

298

299 The segments completely partition the time interval [t_start, t_end].

300 """

301 # Find all waves that cross outlet in this time range

302 outlet_events = []

303

304 # Track rarefactions that already contain the outlet at t_start

305 # These need to be handled separately since they don't generate crossing events

306 active_rarefactions_at_start = []

307

308 for wave in waves:

309 if not wave.is_active:

310 continue

311

312 # For rarefactions, detect both head and tail crossings

313 if isinstance(wave, RarefactionWave):

314 # Check if outlet is already inside this rarefaction at t_start

315 if wave.contains_point(v_outlet, t_start):

316 active_rarefactions_at_start.append(wave)

317 # Don't add crossing events for this wave since we're already inside it

318 # But we still need to detect when the tail crosses during [t_start, t_end]

319 v_tail = wave.tail_position_at_time(t_start)

320 if v_tail is not None and v_tail < v_outlet:

321 vel_tail = wave.tail_velocity()

322 if vel_tail > 0:

323 dt = (v_outlet - v_tail) / vel_tail

324 t_cross = t_start + dt

325 if t_start < t_cross <= t_end:

326 outlet_events.append({

327 "time": t_cross,

328 "wave": wave,

329 "boundary": "tail",

330 "c_after": wave.c_tail,

331 })

332 continue

333

334 # Head crossing

335 v_head = wave.head_position_at_time(t_start)

336 if v_head is not None and v_head < v_outlet:

337 vel_head = wave.head_velocity()

338 if vel_head > 0:

339 dt = (v_outlet - v_head) / vel_head

340 t_cross = t_start + dt

341 if t_start <= t_cross <= t_end:

342 outlet_events.append({

343 "time": t_cross,

344 "wave": wave,

345 "boundary": "head",

346 "c_after": wave.c_head,

347 })

348

349 # Tail crossing

350 v_tail = wave.tail_position_at_time(t_start)

351 if v_tail is not None and v_tail < v_outlet:

352 vel_tail = wave.tail_velocity()

353 if vel_tail > 0:

354 dt = (v_outlet - v_tail) / vel_tail

355 t_cross = t_start + dt

356 if t_start <= t_cross <= t_end:

357 outlet_events.append({

358 "time": t_cross,

359 "wave": wave,

360 "boundary": "tail",

361 "c_after": wave.c_tail,

362 })

363 else:

364 # Characteristics and shocks

365 t_cross = find_outlet_crossing(wave, v_outlet, t_start)

366

367 if t_cross is not None and t_start <= t_cross <= t_end:

368 if isinstance(wave, CharacteristicWave):

369 c_after = wave.concentration

370 elif isinstance(wave, ShockWave):

371 # After shock passes outlet, outlet sees left (upstream) state

372 c_after = wave.c_left

373 else:

374 c_after = 0.0

375

376 outlet_events.append({"time": t_cross, "wave": wave, "boundary": None, "c_after": c_after})

377

378 # Sort events by time

379 outlet_events.sort(key=itemgetter("time"))

380

381 # Create segments between events

382 segments = []

383 current_t = t_start

384 current_c = concentration_at_point(v_outlet, t_start, waves, sorption)

385

386 # Handle case where we start inside a rarefaction

387 if active_rarefactions_at_start:

388 # Should only be one rarefaction containing the outlet at t_start

389 raref = active_rarefactions_at_start[0]

390

391 # Find when tail crosses (if it does)

392 tail_cross_time = None

393 for event in outlet_events:

394 if event["wave"] is raref and event["boundary"] == "tail" and event["time"] > t_start:

395 tail_cross_time = event["time"]

396 break

397

398 # Create rarefaction segment from t_start

399 raref_end = min(tail_cross_time or t_end, t_end)

400

401 c_start = concentration_at_point(v_outlet, t_start, waves, sorption)

402 c_end = None

403 if tail_cross_time and tail_cross_time <= t_end:

404 c_end = raref.c_tail

405

406 segments.append({

407 "t_start": t_start,

408 "t_end": raref_end,

409 "type": "rarefaction",

410 "wave": raref,

411 "c_start": c_start,

412 "c_end": c_end,

413 })

414

415 current_t = raref_end

416 current_c = (

417 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < t_end else current_c

418 )

419

420 for event in outlet_events:

421 # Check if we're entering a rarefaction fan

422 if isinstance(event["wave"], RarefactionWave) and event["boundary"] == "head":

423 # Before rarefaction head: constant segment

424 if event["time"] > current_t:

425 segments.append({

426 "t_start": current_t,

427 "t_end": event["time"],

428 "type": "constant",

429 "concentration": current_c,

430 "c_start": current_c,

431 "c_end": current_c,

432 })

433

434 # Find when tail crosses (if it does)

435 raref = event["wave"]

436 tail_cross_time = None

437

438 for later_event in outlet_events:

439 if (

440 later_event["wave"] is raref

441 and later_event["boundary"] == "tail"

442 and later_event["time"] > event["time"]

443 ):

444 tail_cross_time = later_event["time"]

445 break

446

447 # Rarefaction segment

448 raref_end = min(tail_cross_time or t_end, t_end)

449

450 segments.append({

451 "t_start": event["time"],

452 "t_end": raref_end,

453 "type": "rarefaction",

454 "wave": raref,

455 "c_start": raref.c_head,

456 "c_end": raref.c_tail if tail_cross_time and tail_cross_time <= t_end else None,

457 })

458

459 current_t = raref_end

460 current_c = (

461 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < t_end else current_c

462 )

463 else:

464 # Regular event (characteristic or shock crossing)

465 # Segment before event

466 if event["time"] > current_t:

467 segments.append({

468 "t_start": current_t,

469 "t_end": event["time"],

470 "type": "constant",

471 "concentration": current_c,

472 "c_start": current_c,

473 "c_end": current_c,

474 })

475

476 current_t = event["time"]

477 current_c = event["c_after"]

478

479 # Final segment

480 if t_end > current_t:

481 segments.append({

482 "t_start": current_t,

483 "t_end": t_end,

484 "type": "constant",

485 "concentration": current_c,

486 "c_start": current_c,

487 "c_end": current_c,

488 })

489

490 return segments

491

492

493def integrate_rarefaction_exact(

494 raref: RarefactionWave, v_outlet: float, t_start: float, t_end: float, sorption: FreundlichSorption

495) -> float:

496 """

497 Exact analytical integral of rarefaction concentration over time at fixed position.

498

499 For Freundlich sorption, the concentration within a rarefaction fan varies as:

500 R(C) = flow*(t - t_origin)/(v_outlet - v_origin)

501

502 This function computes the exact integral:

503 ∫_{t_start}^{t_end} C(t) dt

504

505 where C(t) is obtained by inverting the retardation relation.

506

507 Parameters

508 ----------

509 raref : RarefactionWave

510 Rarefaction wave controlling the outlet.

511 v_outlet : float

512 Outlet position [m³].

513 t_start : float

514 Integration start time [days]. Can be -np.inf.

515 t_end : float

516 Integration end time [days]. Can be np.inf.

517 sorption : FreundlichSorption

518 Freundlich sorption model.

519

520 Returns

521 -------

522 integral : float

523 Exact integral value [concentration * time].

524

525 Notes

526 -----

527 **Derivation**:

528

529 For Freundlich: R(C) = 1 + alpha*C^beta where:

530 alpha = rho_b*k_f/(n_por*n)

531 beta = 1/n - 1

532

533 At outlet: R = kappa*t + mu where:

534 kappa = flow/(v_outlet - v_origin)

535 mu = -flow*t_origin/(v_outlet - v_origin)

536

537 Inverting: C = [(R-1)/alpha]^(1/beta) = [(kappa*t + mu - 1)/alpha]^(1/beta)

538

539 Integral:

540 ∫ C dt = (1/alpha^(1/beta)) * ∫ (kappa*t + mu - 1)^(1/beta) dt

541 = (1/alpha^(1/beta)) * (1/kappa) * [(kappa*t + mu - 1)^(1/beta + 1) / (1/beta + 1)]

542

543 evaluated from t_start to t_end.

544

545 **Special Cases**:

546 - If (kappa*t + mu - 1) <= 0, concentration is 0 (unphysical region)

547 - For beta = 0 (n = 1), use ConstantRetardation instead

548 - For t_end = +∞ with exponent < 0 (n>1), integral converges to 0

549 - For t_start = -∞, antiderivative evaluates to 0

550

551 Examples

552 --------

553 ::

554

555 sorption = FreundlichSorption(

556 k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3

557 )

558 raref = RarefactionWave(

559 t_start=0.0,

560 v_start=0.0,

561 flow=100.0,

562 c_head=10.0,

563 c_tail=2.0,

564 sorption=sorption,

565 )

566 integral = integrate_rarefaction_exact(

567 raref, v_outlet=500.0, t_start=5.0, t_end=15.0, sorption=sorption

568 )

569 integral > 0

570 # Integration to infinity

571 integral_inf = integrate_rarefaction_exact(

572 raref, v_outlet=500.0, t_start=5.0, t_end=np.inf, sorption=sorption

573 )

574 integral_inf > integral

575 """

576 # Extract parameters

577 t_origin = raref.t_start

578 v_origin = raref.v_start

579 flow = raref.flow

580

581 # Coefficients in R = kappa*t + μ

582 kappa = flow / (v_outlet - v_origin)

583 mu = -flow * t_origin / (v_outlet - v_origin)

584

585 # Freundlich parameters: R = 1 + alpha*C^beta

586 alpha = sorption.bulk_density * sorption.k_f / (sorption.porosity * sorption.n)

587 beta = 1.0 / sorption.n - 1.0

588

589 # For integration, we need exponent = 1/beta + 1

590 if abs(beta) < EPSILON_BETA:

591 # Linear case - shouldn't happen with Freundlich

592 # Fall back to numerical integration or raise error

593 msg = "integrate_rarefaction_exact requires nonlinear sorption (n != 1)"

594 raise ValueError(msg)

595

596 exponent = 1.0 / beta + 1.0

597

598 # Coefficient for antiderivative

599 # F(t) = (1/(alpha^(1/beta)*kappa*exponent)) * (kappa*t + mu - 1)^exponent

600 coeff = 1.0 / (alpha ** (1.0 / beta) * kappa * exponent)

601

602 def antiderivative(t: float) -> float:

603 """Evaluate antiderivative at time t."""

604 # Handle infinite bounds

605 if np.isinf(t):

606 if t > 0:

607 # t = +∞

608 if exponent < 0:

609 # For n > 1 (higher C travels faster, beta < 0, exponent < 0):

610 # As t → +∞, base → +∞, so base^exponent → 0

611 return 0.0

612 # For n < 1 (lower C travels faster, beta > 0, exponent > 0):

613 # As t → +∞, base → +∞, so base^exponent → +∞

614 # This should not happen for physical rarefactions to C=0

615 msg = f"Integral diverges at t=+∞ with exponent={exponent} > 0"

616 raise ValueError(msg)

617 # t = -∞: always returns 0

618 return 0.0

619

620 base = kappa * t + mu - 1.0

621

622 # Check if we're in physical region (R > 1)

623 if base <= 0:

624 return 0.0

625

626 return coeff * base**exponent

627

628 # Evaluate definite integral

629 return antiderivative(t_end) - antiderivative(t_start)

630

631

632def compute_bin_averaged_concentration_exact(

633 t_edges: npt.NDArray[np.floating],

634 v_outlet: float,

635 waves: list[Wave],

636 sorption: FreundlichSorption | ConstantRetardation,

637) -> npt.NDArray[np.floating]:

638 """

639 Compute bin-averaged concentration using EXACT analytical integration.

640

641 For each time bin [t_i, t_{i+1}], computes:

642 C_avg = (1/(t_{i+1} - t_i)) * ∫_{t_i}^{t_{i+1}} C(v_outlet, t) dt

643

644 This is the critical function for maintaining machine precision in output.

645 All integrations use exact analytical formulas with no numerical quadrature.

646

647 Parameters

648 ----------

649 t_edges : array-like

650 Time bin edges [days]. Length N+1 for N bins.

651 v_outlet : float

652 Outlet position [m³].

653 waves : list of Wave

654 All waves from front tracking simulation.

655 sorption : FreundlichSorption or ConstantRetardation

656 Sorption model.

657

658 Returns

659 -------

660 c_avg : numpy.ndarray

661 Bin-averaged concentrations [mass/volume]. Length N.

662

663 Notes

664 -----

665 **Algorithm**:

666

667 1. For each bin [t_i, t_{i+1}]:

668

669 a. Identify which wave segments control outlet during this period

670 b. For each segment, compute: Constant C gives integral = C * Δt,

671 Rarefaction C(t) uses exact analytical integral formula

672 c. Sum segment integrals and divide by bin width

673

674 **Machine Precision**:

675

676 - Constant segments: exact to floating-point precision

677 - Rarefaction segments: uses closed-form antiderivative

678 - No numerical quadrature or interpolation

679 - Maintains mass balance to ~1e-14 relative error

680

681 **Rarefaction Integration**:

682

683 For Freundlich sorption, rarefaction concentration at outlet varies as::

684

685 C(t) = [(kappa*t + mu - 1)/alpha]^(1/beta)

686

687 The exact integral is::

688

689 ∫ C dt = (1/(alpha^(1/beta)*kappa*exponent)) * (kappa*t + mu - 1)^exponent

690

691 where exponent = 1/beta + 1.

692

693 Examples

694 --------

695 ::

696

697 # After running front tracking simulation

698 t_edges = np.array([0.0, 10.0, 20.0, 30.0])

699 c_avg = compute_bin_averaged_concentration_exact(

700 t_edges=t_edges,

701 v_outlet=500.0,

702 waves=tracker.state.waves,

703 sorption=sorption,

704 )

705 len(c_avg) == len(t_edges) - 1

706 np.all(c_avg >= 0)

707

708 See Also

709 --------

710 concentration_at_point : Point-wise concentration

711 compute_breakthrough_curve : Breakthrough curve

712 integrate_rarefaction_exact : Exact rarefaction integration

713 """

714 t_edges = np.asarray(t_edges, dtype=float)

715 n_bins = len(t_edges) - 1

716 c_avg = np.zeros(n_bins)

717

718 for i in range(n_bins):

719 t_start = t_edges[i]

720 t_end = t_edges[i + 1]

721 dt = t_end - t_start

722

723 if dt <= 0:

724 msg = f"Invalid time bin: t_edges[{i}]={t_start} >= t_edges[{i + 1}]={t_end}"

725 raise ValueError(msg)

726

727 # Identify wave segments controlling outlet in this time bin

728 segments = identify_outlet_segments(t_start, t_end, v_outlet, waves, sorption)

729

730 # Integrate each segment

731 total_integral = 0.0

732

733 for seg in segments:

734 seg_t_start = max(seg["t_start"], t_start)

735 seg_t_end = min(seg["t_end"], t_end)

736 seg_dt = seg_t_end - seg_t_start

737

738 if seg_dt <= EPSILON_TIME: # Skip negligible segments

739 continue

740

741 if seg["type"] == "constant":

742 # C is constant over segment - exact integral

743 integral = seg["concentration"] * seg_dt

744

745 elif seg["type"] == "rarefaction":

746 # C(t) given by self-similar solution - use exact analytical integral

747 if isinstance(sorption, FreundlichSorption):

748 raref = seg["wave"]

749 integral = integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)

750 else:

751 # ConstantRetardation - rarefactions shouldn't form, use constant approximation

752 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)

753 integral = c_mid * seg_dt

754 else:

755 msg = f"Unknown segment type: {seg['type']}"

756 raise ValueError(msg)

757

758 total_integral += integral

759

760 c_avg[i] = total_integral / dt

761

762 return c_avg

763

764

765def compute_domain_mass(

766 t: float,

767 v_outlet: float,

768 waves: list[Wave],

769 sorption: FreundlichSorption | ConstantRetardation,

770) -> float:

771 """

772 Compute total mass in domain [0, v_outlet] at time t using exact analytical integration.

773

774 Implements runtime mass balance verification as described in High Priority #3

775 of FRONT_TRACKING_REBUILD_PLAN.md. Integrates concentration over space:

776

777 M(t) = ∫₀^v_outlet C(v, t) dv

778

779 using exact analytical formulas for each wave type.

780

781 Parameters

782 ----------

783 t : float

784 Time at which to compute domain mass [days].

785 v_outlet : float

786 Outlet position (domain extent) [m³].

787 waves : list of Wave

788 All waves in the simulation.

789 sorption : FreundlichSorption or ConstantRetardation

790 Sorption model.

791

792 Returns

793 -------

794 mass : float

795 Total mass in domain [mass]. Computed to machine precision (~1e-14).

796

797 Notes

798 -----

799 **Algorithm**:

800

801 1. Partition spatial domain [0, v_outlet] into segments where concentration

802 is controlled by a single wave or is constant.

803 2. For each segment, compute mass analytically:

804 - Constant C: mass = C * Δv

805 - Rarefaction C(v): use exact spatial integration formula

806 3. Sum all segment masses.

807

808 **Wave Priority** (same as concentration_at_point):

809

810 1. Rarefactions (if position is inside rarefaction fan)

811 2. Shocks and rarefaction tails (most recent to pass)

812 3. Characteristics (smooth regions)

813

814 **Rarefaction Spatial Integration**:

815

816 For a rarefaction fan at fixed time t, concentration varies with position v

817 according to the self-similar solution:

818

819 R(C) = flow*(t - t_origin)/(v - v_origin)

820

821 The spatial integral ∫ C dv is computed analytically using the inverse

822 retardation relation.

823

824 **Integration Precision**:

825

826 - Constant concentration regions: Exact analytical (C_total * dv)

827 - Rarefaction regions: High-precision trapezoidal quadrature (10+ points)

828 - Overall accuracy: ~1e-10 to 1e-12 relative error

829 - Sufficient for runtime verification; primary outputs remain exact

830

831 Examples

832 --------

833 ::

834

835 # After running simulation to time t=10.0

836 mass = compute_domain_mass(

837 t=10.0, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption

838 )

839 mass >= 0.0

840

841 See Also

842 --------

843 compute_cumulative_inlet_mass : Cumulative inlet mass

844 compute_cumulative_outlet_mass : Cumulative outlet mass

845 concentration_at_point : Point-wise concentration

846 """

847 # Partition spatial domain into segments based on wave structure

848 # We'll evaluate concentration at many points and identify constant regions

849

850 # Collect all wave positions at time t

851 wave_positions = []

852

853 for wave in waves:

854 if not wave.is_active:

855 continue

856

857 if isinstance(wave, (CharacteristicWave, ShockWave)):

858 v_pos = wave.position_at_time(t)

859 if v_pos is not None and 0 <= v_pos <= v_outlet:

860 wave_positions.append(v_pos)

861

862 elif isinstance(wave, RarefactionWave):

863 v_head = wave.head_position_at_time(t)

864 v_tail = wave.tail_position_at_time(t)

865

866 if v_head is not None and 0 <= v_head <= v_outlet:

867 wave_positions.append(v_head)

868 if v_tail is not None and 0 <= v_tail <= v_outlet:

869 wave_positions.append(v_tail)

870

871 # Add domain boundaries

872 wave_positions.extend([0.0, v_outlet])

873

874 # Sort and remove duplicates

875 wave_positions = sorted(set(wave_positions))

876

877 # Remove positions outside domain

878 wave_positions = [v for v in wave_positions if 0 <= v <= v_outlet]

879

880 # Compute mass in each segment using refined integration

881 total_mass = 0.0

882

883 for i in range(len(wave_positions) - 1):

884 v_start = wave_positions[i]

885 v_end = wave_positions[i + 1]

886 dv = v_end - v_start

887

888 if dv < EPSILON_VELOCITY:

889 continue

890

891 # Check if this segment is inside a rarefaction fan

892 # Sample at midpoint to check

893 v_mid = 0.5 * (v_start + v_end)

894 inside_rarefaction = False

895 raref_wave = None

896

897 for wave in waves:

898 if isinstance(wave, RarefactionWave) and wave.is_active and wave.contains_point(v_mid, t):

899 inside_rarefaction = True

900 raref_wave = wave

901 break

902

903 if inside_rarefaction and raref_wave is not None:

904 # Rarefaction: concentration varies with position

905 # Use EXACT analytical spatial integration

906 mass_segment = _integrate_rarefaction_spatial_exact(raref_wave, v_start, v_end, t, sorption)

907 else:

908 # Constant concentration region - exact integration

909 v_mid = 0.5 * (v_start + v_end)

910 c = concentration_at_point(v_mid, t, waves, sorption)

911 c_total = sorption.total_concentration(c)

912 mass_segment = c_total * dv

913

914 total_mass += mass_segment

915

916 return total_mass

917

918

919def _integrate_rarefaction_spatial_exact(

920 raref: RarefactionWave,

921 v_start: float,

922 v_end: float,

923 t: float,

924 sorption: FreundlichSorption | ConstantRetardation,

925) -> float:

926 """

927 EXACT analytical spatial integral of rarefaction total concentration at fixed time.

928

929 Computes ∫_{v_start}^{v_end} C_total(v) dv analytically using closed-form

930 antiderivatives for specific Freundlich n values. This maintains machine precision

931 for the mass balance diagnostic.

932

933 For rarefaction at time t: R(C) = kappa/(v - v₀) where kappa = flow*(t - t₀)

934

935 For Freundlich: R = 1 + alpha*C^beta where alpha = rho_b*k_f/(n_por*n), beta = 1/n - 1

936 Total concentration: C_total = C + (rho_b/n_por)*k_f*C^(1/n)

937

938 **Exact formulas implemented**:

939 - n = 2 (beta = -0.5): Closed-form using polynomial expansion (optimized path)

940 - n = 1: No rarefactions (constant retardation)

941 - All other n > 0: Unified formula using generalized incomplete beta function via mpmath

942

943 Parameters

944 ----------

945 raref : RarefactionWave

946 v_start, v_end : float

947 Integration bounds [m³]

948 t : float

949 Time [days]

950 sorption : FreundlichSorption or ConstantRetardation

951

952 Returns

953 -------

954 mass : float

955 Exact mass in segment to machine precision

956

957 Notes

958 -----

959 **Derivation for n=2** (beta = -0.5):

960

961 C(v) = [(kappa/(v-v₀) - 1) / alpha]^(-2) = [alpha*(v-v₀)]² / (kappa - (v-v₀))²

962

963 Let u = v - v₀, then ∫ C du requires integrating u²/(kappa-u)².

964 Using substitution w = kappa - u:

965

966 ∫ u²/(kappa-u)² du = alpha²[kappa²/(kappa-u) + 2kappa*ln|kappa-u| - (kappa-u)]

967

968 For C_total, we integrate both C and (rho_b/n_por)*k_f*C^0.5 terms.

969

970 **Unified Formula for All n > 0**:

971

972 The spatial integral has the structure:

973

974 ∫ C_total(u) du = ∫ C(u) du + ∫ (rho_b/n_por)*k_f*C(u)^(1/n) du

975

976 where C(u) = [(kappa-u)/(alpha*u)]^(1/beta) with beta = 1/n - 1.

977

978 Both integrals reduce to power-law forms:

979

980 ∫ u^p (kappa-u)^q du

981

982 which can be expressed using the generalized incomplete beta function:

983

984 ∫_{u₁}^{u₂} u^p (kappa-u)^q du = kappa^(p+q+1) B(u₁/kappa, u₂/kappa; p+1, q+1)

985

986 where B(x₁, x₂; a, b) = ∫_{x₁}^{x₂} t^(a-1) (1-t)^(b-1) dt.

987

988 For many n values, the parameters a and b can be negative or zero, which

989 requires analytic continuation beyond the standard incomplete beta function.

990 The mpmath library provides this through mpmath.betainc(), which handles

991 all parameter values and achieves machine precision (~1e-15 relative error).

992

993 **Mathematical Properties**:

994 - Continuous for all n > 0 (no discontinuities)

995 - Achieves machine precision (~1e-14 to 1e-15 relative error)

996 - No numerical quadrature (uses special function evaluation)

997 - Single unified formula (no conditional logic based on n)

998 """

999 if isinstance(sorption, ConstantRetardation):

1000 # Constant retardation: no rarefactions

1001 v_mid = 0.5 * (v_start + v_end)

1002 c = raref.concentration_at_point(v_mid, t) or 0.0

1003 c_total = sorption.total_concentration(c)

1004 return c_total * (v_end - v_start)

1005

1006 # Extract parameters

1007 t_origin = raref.t_start

1008 v_origin = raref.v_start

1009 flow = raref.flow

1010

1011 if t <= t_origin:

1012 return 0.0

1013

1014 kappa = flow * (t - t_origin) # kappa

1015 alpha = sorption.bulk_density * sorption.k_f / (sorption.porosity * sorption.n) # alpha

1016 rho_b = sorption.bulk_density

1017 n_por = sorption.porosity

1018 k_f = sorption.k_f

1019 n = sorption.n

1020

1021 # Transform to u coordinates

1022 u_start = v_start - v_origin

1023 u_end = v_end - v_origin

1024

1025 if u_start <= 0 or u_end <= 0:

1026 return 0.0

1027

1028 # Unified formula for all n > 0 using generalized incomplete beta function

1029 beta = 1 / n - 1 # beta parameter

1030

1031 # Transform to [0,1] domain: t = u/kappa (use regular floats)

1032 t_start = u_start / kappa

1033 t_end = u_end / kappa

1034

1035 # Integral 1: Dissolved concentration

1036 # ∫ C(u) du = ∫ [(kappa-u)/(alpha*u)]^(1/beta) du

1037 # = (1/alpha)^(1/beta) ∫ u^(-1/beta) (kappa-u)^(1/beta) du

1038 # = (1/alpha)^(1/beta) * kappa * B(t_start, t_end; 1-1/beta, 1+1/beta)

1039 #

1040 # where B(x1, x2; a, b) = ∫_{x1}^{x2} t^(a-1) (1-t)^(b-1) dt

1041

1042 # Beta function parameters (use regular floats)

1043 a_diss = 1 - 1 / beta

1044 b_diss = 1 + 1 / beta

1045

1046 # Use scipy.special.betainc when parameters are positive (faster),

1047 # fall back to mpmath for negative parameters (requires analytic continuation)

1048 if a_diss > 0 and b_diss > 0:

1049 # scipy.special.betainc is regularized, so multiply by beta function

1050 beta_diss = betainc(a_diss, b_diss, t_end) - betainc(a_diss, b_diss, t_start)

1051 beta_diss *= beta_func(a_diss, b_diss)

1052 else:

1053 # Use mpmath for negative parameters (analytic continuation)

1054 beta_diss = float(mp.betainc(a_diss, b_diss, t_start, t_end, regularized=False))

1055

1056 # Compute coefficient using regular float arithmetic

1057 coeff_diss = (1 / alpha) ** (1 / beta)

1058 mass_dissolved = coeff_diss * kappa * beta_diss

1059

1060 # Integral 2: Sorbed concentration

1061 # ∫ (rho_b/n_por)*k_f*C^(1/n) du

1062 #

1063 # where C^(1/n) = [(kappa-u)/(alpha*u)]^(1/(1-n))

1064 #

1065 # This gives: ∫ u^(-1/(1-n)) (kappa-u)^(1/(1-n)) du

1066 #

1067 # Using same transformation: = kappa * B(t_start, t_end; 1-1/(1-n), 1+1/(1-n))

1068

1069 # Beta function parameters (use regular floats)

1070 exponent_sorb = 1 / (1 - n)

1071 a_sorb = 1 - exponent_sorb

1072 b_sorb = 1 + exponent_sorb

1073

1074 # Use scipy when possible, mpmath for negative parameters

1075 if a_sorb > 0 and b_sorb > 0:

1076 beta_sorb = betainc(a_sorb, b_sorb, t_end) - betainc(a_sorb, b_sorb, t_start)

1077 beta_sorb *= beta_func(a_sorb, b_sorb)

1078 else:

1079 beta_sorb = float(mp.betainc(a_sorb, b_sorb, t_start, t_end, regularized=False))

1080

1081 # Compute coefficient using regular float arithmetic

1082 coeff_sorb = (rho_b / n_por) * k_f / (alpha**exponent_sorb)

1083 mass_sorbed = coeff_sorb * kappa * beta_sorb

1084

1085 return mass_dissolved + mass_sorbed

1086

1087

1088def compute_cumulative_inlet_mass(

1089 t: float,

1090 cin: npt.NDArray[np.floating],

1091 flow: npt.NDArray[np.floating],

1092 tedges_days: npt.NDArray[np.floating],

1093) -> float:

1094 """

1095 Compute cumulative mass entering domain from t=0 to t.

1096

1097 Integrates inlet mass flux over time:

1098 M_in(t) = ∫₀^t cin(τ) * flow(τ) dτ

1099

1100 using exact analytical integration of piecewise-constant functions.

1101

1102 Parameters

1103 ----------

1104 t : float

1105 Time up to which to integrate [days].

1106 cin : array-like

1107 Inlet concentration time series [mass/volume].

1108 Piecewise constant within bins defined by tedges_days.

1109 flow : array-like

1110 Flow rate time series [m³/day].

1111 Piecewise constant within bins defined by tedges_days.

1112 tedges_days : numpy.ndarray

1113 Time bin edges [days]. Length len(cin) + 1.

1114

1115 Returns

1116 -------

1117 mass_in : float

1118 Cumulative mass entered [mass].

1119

1120 Notes

1121 -----

1122 For piecewise-constant cin and flow:

1123 M_in = Σ cin[i] * flow[i] * dt[i]

1124

1125 where the sum is over all bins from tedges_days[0] to t.

1126 Partial bins are handled exactly.

1127

1128 Examples

1129 --------

1130 ::

1131

1132 mass_in = compute_cumulative_inlet_mass(

1133 t=50.0, cin=cin, flow=flow, tedges_days=tedges_days

1134 )

1135 mass_in >= 0.0

1136 """

1137 tedges_arr = np.asarray(tedges_days, dtype=float)

1138 cin_arr = np.asarray(cin, dtype=float)

1139 flow_arr = np.asarray(flow, dtype=float)

1140

1141 # Find which bin t falls into

1142 if t <= tedges_arr[0]:

1143 return 0.0

1144

1145 if t >= tedges_arr[-1]:

1146 # Integrate all bins

1147 dt = np.diff(tedges_arr)

1148 return float(np.sum(cin_arr * flow_arr * dt))

1149

1150 # Find bin containing t

1151 bin_idx = np.searchsorted(tedges_arr, t, side="right") - 1

1152

1153 # Mass flux across inlet boundary = Q * C_in (aqueous concentration)

1154 # This is correct for sorbing solutes: only dissolved mass flows with water

1155 # Integrate complete bins before t

1156 if bin_idx > 0:

1157 dt_complete = np.diff(tedges_arr[: bin_idx + 1])

1158 mass_complete = np.sum(cin_arr[:bin_idx] * flow_arr[:bin_idx] * dt_complete)

1159 else:

1160 mass_complete = 0.0

1161

1162 # Add partial bin

1163 if bin_idx >= 0 and bin_idx < len(cin_arr):

1164 dt_partial = t - tedges_arr[bin_idx]

1165 mass_partial = cin_arr[bin_idx] * flow_arr[bin_idx] * dt_partial

1166 else:

1167 mass_partial = 0.0

1168

1169 return float(mass_complete + mass_partial)

1170

1171

1172def find_last_rarefaction_start_time(

1173 v_outlet: float,

1174 waves: list[Wave],

1175) -> float:

1176 """

1177 Find the time when the last rarefaction head reaches the outlet.

1178

1179 For rarefactions, we integrate analytically so we only need to know

1180 when the rarefaction starts at the outlet (head arrival).

1181

1182 Parameters

1183 ----------

1184 v_outlet : float

1185 Outlet position [m³].

1186 waves : list of Wave

1187 All waves in the simulation.

1188

1189 Returns

1190 -------

1191 t_last : float

1192 Time when last rarefaction head reaches outlet [days].

1193 For non-rarefaction waves, uses their arrival time.

1194 Returns 0.0 if no waves reach the outlet.

1195

1196 Notes

1197 -----

1198 This function finds when the last wave structure *starts* at the outlet.

1199 For rarefactions, this is the head arrival time. The tail may arrive

1200 much later (or at infinite time for rarefactions to C=0), but the

1201 total mass in the rarefaction is computed analytically.

1202

1203 Examples

1204 --------

1205 ::

1206

1207 t_last = find_last_rarefaction_start_time(v_outlet=500.0, waves=all_waves)

1208 t_last >= 0.0

1209 """

1210 t_last = 0.0

1211

1212 for wave in waves:

1213 if not wave.is_active:

1214 continue

1215

1216 if isinstance(wave, RarefactionWave):

1217 # For rarefaction, use head arrival (when rarefaction starts)

1218 head_vel = wave.head_velocity()

1219 if head_vel > EPSILON_VELOCITY:

1220 t_arrival = wave.t_start + (v_outlet - wave.v_start) / head_vel

1221 t_last = max(t_last, t_arrival)

1222 elif isinstance(wave, (CharacteristicWave, ShockWave)):

1223 # For characteristics and shocks, compute arrival time

1224 vel = wave.velocity if isinstance(wave, ShockWave) else wave.velocity()

1225 if vel is not None and vel > EPSILON_VELOCITY:

1226 t_arrival = wave.t_start + (v_outlet - wave.v_start) / vel

1227 t_last = max(t_last, t_arrival)

1228

1229 return t_last

1230

1231

1232def compute_cumulative_outlet_mass(

1233 t: float,

1234 v_outlet: float,

1235 waves: list[Wave],

1236 sorption: FreundlichSorption | ConstantRetardation,

1237 flow: npt.NDArray[np.floating],

1238 tedges_days: npt.NDArray[np.floating],

1239) -> float:

1240 """

1241 Compute cumulative mass exiting domain from t=0 to t.

1242

1243 Integrates outlet mass flux over time:

1244 M_out(t) = ∫₀^t cout(τ) * flow(τ) dτ

1245

1246 using exact analytical integration. Outlet concentration cout(τ) is obtained

1247 from the wave solution, and flow(τ) is piecewise constant.

1248

1249 Parameters

1250 ----------

1251 t : float

1252 Time up to which to integrate [days].

1253 v_outlet : float

1254 Outlet position [m³].

1255 waves : list of Wave

1256 All waves in the simulation.

1257 sorption : FreundlichSorption or ConstantRetardation

1258 Sorption model.

1259 flow : array-like

1260 Flow rate time series [m³/day].

1261 Piecewise constant within bins defined by tedges_days.

1262 tedges_days : numpy.ndarray

1263 Time bin edges [days]. Length len(flow) + 1.

1264

1265 Returns

1266 -------

1267 mass_out : float

1268 Cumulative mass exited [mass].

1269

1270 Notes

1271 -----

1272 The outlet concentration is obtained from wave solution via

1273 concentration_at_point(v_outlet, τ, waves, sorption).

1274

1275 For each flow bin [t_i, t_{i+1}], the mass flux integral is computed

1276 exactly using identify_outlet_segments and analytical integration

1277 (constant segments and rarefaction segments).

1278

1279 Examples

1280 --------

1281 ::

1282

1283 mass_out = compute_cumulative_outlet_mass(

1284 t=50.0,

1285 v_outlet=500.0,

1286 waves=tracker.state.waves,

1287 sorption=sorption,

1288 flow=flow,

1289 tedges_days=tedges_days,

1290 )

1291 mass_out >= 0.0

1292 """

1293 tedges_arr = np.asarray(tedges_days, dtype=float)

1294 flow_arr = np.asarray(flow, dtype=float)

1295

1296 if t <= tedges_arr[0]:

1297 return 0.0

1298

1299 # Integrate bin by bin through all flow bins, then continue to t if needed

1300 total_mass = 0.0

1301

1302 # Process bins within tedges range

1303 n_flow_bins = len(flow_arr)

1304

1305 for i in range(n_flow_bins):

1306 t_bin_start = tedges_arr[i]

1307 t_bin_end = tedges_arr[i + 1]

1308

1309 # Skip bins entirely before t

1310 if t_bin_end <= tedges_arr[0]:

1311 continue

1312

1313 # Clip to [tedges[0], t]

1314 t_bin_start = max(t_bin_start, tedges_arr[0])

1315 t_bin_end = min(t_bin_end, t)

1316

1317 flow_i = flow_arr[i]

1318 dt_i = t_bin_end - t_bin_start

1319

1320 if dt_i <= 0:

1321 continue

1322

1323 # Compute ∫_{t_bin_start}^{t_bin_end} cout(τ) dτ using exact integration

1324 segments = identify_outlet_segments(t_bin_start, t_bin_end, v_outlet, waves, sorption)

1325

1326 integral_c_dt = 0.0

1327

1328 for seg in segments:

1329 seg_t_start = max(seg["t_start"], t_bin_start)

1330 seg_t_end = min(seg["t_end"], t_bin_end)

1331 seg_dt = seg_t_end - seg_t_start

1332

1333 if seg_dt <= EPSILON_TIME:

1334 continue

1335

1336 if seg["type"] == "constant":

1337 integral_c_dt += seg["concentration"] * seg_dt

1338 elif seg["type"] == "rarefaction":

1339 if isinstance(sorption, FreundlichSorption):

1340 raref = seg["wave"]

1341 integral_c_dt += integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)

1342 else:

1343 # ConstantRetardation - use midpoint

1344 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)

1345 integral_c_dt += c_mid * seg_dt

1346

1347 # Mass for this bin = flow * ∫ cout dt

1348 total_mass += flow_i * integral_c_dt

1349

1350 return float(total_mass)

1351

1352

1353def integrate_rarefaction_total_mass(

1354 raref: RarefactionWave,

1355 v_outlet: float,

1356 t_start: float,

1357 sorption: FreundlichSorption | ConstantRetardation,

1358 flow: float,

1359) -> float:

1360 """

1361 Compute total mass exiting through a rarefaction.

1362

1363 For a rarefaction that starts at the outlet at time t_start, compute the

1364 total mass that will exit through the rarefaction. Integration endpoint

1365 depends on the rarefaction tail concentration:

1366

1367 - If c_tail ≈ 0: Integrate to infinity (tail extends infinitely)

1368 - If c_tail > 0: Integrate to finite tail arrival time

1369

1370 Parameters

1371 ----------

1372 raref : RarefactionWave

1373 Rarefaction wave.

1374 v_outlet : float

1375 Outlet position [m³].

1376 t_start : float

1377 Time when rarefaction head reaches outlet [days].

1378 sorption : FreundlichSorption or ConstantRetardation

1379 Sorption model.

1380 flow : float

1381 Flow rate [m³/day] (assumed constant).

1382

1383 Returns

1384 -------

1385 total_mass : float

1386 Total mass that exits through rarefaction [mass].

1387

1388 Notes

1389 -----

1390 Uses the exact analytical integral:

1391 M_total = ∫_{t_start}^{t_end} Q * C(t) dt

1392

1393 where C(t) is the concentration at the outlet from the rarefaction wave.

1394

1395 For n > 1 (favorable sorption), rarefactions typically decrease to C=0,

1396 so t_end = ∞ and the integral converges.

1397

1398 For n < 1 (unfavorable sorption), rarefactions typically increase from

1399 low C to high C, so c_tail > 0 and the tail arrives at finite time

1400 t_end = t_start + (v_outlet - v_start) / tail_velocity.

1401

1402 Examples

1403 --------

1404 ::

1405

1406 mass = integrate_rarefaction_total_mass(

1407 raref=raref_wave,

1408 v_outlet=500.0,

1409 t_start=40.0,

1410 sorption=sorption,

1411 flow=100.0,

1412 )

1413 mass >= 0.0

1414 """

1415 if isinstance(sorption, ConstantRetardation):

1416 # No rarefactions with constant retardation

1417 return 0.0

1418

1419 # Determine integration endpoint based on c_tail

1420 # For rarefactions with c_tail ≈ 0, the tail extends to infinity

1421 # For rarefactions with c_tail > 0, the tail arrives at finite time

1422 if raref.c_tail < EPSILON_CONCENTRATION:

1423 # Rarefaction tail goes to C≈0, extends to infinite time

1424 # This is typical for n > 1 rarefactions (concentration decreases)

1425 t_end = np.inf

1426 else:

1427 # Rarefaction tail has finite concentration, arrives at finite time

1428 # This is typical for n < 1 rarefactions (concentration increases)

1429 tail_velocity = raref.tail_velocity()

1430 if tail_velocity < EPSILON_VELOCITY:

1431 # Tail velocity is effectively zero, extends to infinity

1432 t_end = np.inf

1433 else:

1434 # Compute finite tail arrival time

1435 t_end = raref.t_start + (v_outlet - raref.v_start) / tail_velocity

1436

1437 # Integrate from t_start to t_end

1438 integral_c_dt = integrate_rarefaction_exact(raref, v_outlet, t_start, t_end, sorption)

1439

1440 return flow * integral_c_dt

1441

1442

1443def compute_total_outlet_mass(

1444 v_outlet: float,

1445 waves: list[Wave],

1446 sorption: FreundlichSorption | ConstantRetardation,

1447 flow: npt.NDArray[np.floating],

1448 tedges_days: npt.NDArray[np.floating],

1449) -> tuple[float, float]:

1450 """

1451 Compute total integrated outlet mass until all mass has exited.

1452

1453 Automatically determines when the last wave passes the outlet and

1454 integrates the outlet mass flux until that time, regardless of the

1455 provided tedges extent.

1456

1457 Parameters

1458 ----------

1459 v_outlet : float

1460 Outlet position [m³].

1461 waves : list of Wave

1462 All waves in the simulation.

1463 sorption : FreundlichSorption or ConstantRetardation

1464 Sorption model.

1465 flow : array-like

1466 Flow rate time series [m³/day].

1467 Piecewise constant within bins defined by tedges_days.

1468 tedges_days : numpy.ndarray

1469 Time bin edges [days]. Length len(flow) + 1.

1470

1471 Returns

1472 -------

1473 total_mass_out : float

1474 Total mass that exits through outlet [mass].

1475 t_integration_end : float

1476 Time until which integration was performed [days].

1477 This is the time when the last wave passes the outlet.

1478

1479 Notes

1480 -----

1481 This function:

1482 1. Finds when the last rarefaction *starts* at the outlet (head arrival)

1483 2. Integrates outlet mass flux until that time

1484 3. Adds analytical integral of rarefaction mass from start to infinity

1485

1486 For rarefactions to C=0, the tail has infinite arrival time but the

1487 total mass is finite and computed analytically.

1488

1489 Examples

1490 --------

1491 ::

1492

1493 total_mass, t_end = compute_total_outlet_mass(

1494 v_outlet=500.0,

1495 waves=tracker.state.waves,

1496 sorption=sorption,

1497 flow=flow,

1498 tedges_days=tedges_days,

1499 )

1500 total_mass >= 0.0

1501 t_end >= tedges_days[0]

1502

1503 See Also

1504 --------

1505 compute_cumulative_outlet_mass : Cumulative outlet mass up to time t

1506 find_last_rarefaction_start_time : Find when last rarefaction starts

1507 integrate_rarefaction_total_mass : Total mass in rarefaction to infinity

1508 """

1509 # Find when the last rarefaction starts at the outlet

1510 t_last_raref_start = find_last_rarefaction_start_time(v_outlet, waves)

1511

1512 tedges_arr = np.asarray(tedges_days, dtype=float)

1513 flow_arr = np.asarray(flow, dtype=float)

1514

1515 # Integrate up to when last rarefaction starts

1516 if t_last_raref_start <= tedges_arr[-1]:

1517 # All waves start within provided time range

1518 mass_up_to_raref_start = compute_cumulative_outlet_mass(

1519 t_last_raref_start, v_outlet, waves, sorption, flow_arr, tedges_arr

1520 )

1521 flow_at_raref_start = flow_arr[-1] # Use last flow value

1522 else:

1523 # Need to extend beyond tedges to reach rarefaction start

1524 # First, compute mass up to tedges[-1]

1525 mass_within_tedges = compute_cumulative_outlet_mass(

1526 tedges_arr[-1], v_outlet, waves, sorption, flow_arr, tedges_arr

1527 )

1528

1529 # Then, integrate from tedges[-1] to t_last_raref_start

1530 flow_extended = flow_arr[-1]

1531 t_start_extended = tedges_arr[-1]

1532 t_end_extended = t_last_raref_start

1533

1534 # Get outlet segments for extended period

1535 segments = identify_outlet_segments(t_start_extended, t_end_extended, v_outlet, waves, sorption)

1536

1537 integral_c_dt = 0.0

1538

1539 for seg in segments:

1540 seg_t_start = max(seg["t_start"], t_start_extended)

1541 seg_t_end = min(seg["t_end"], t_end_extended)

1542 seg_dt = seg_t_end - seg_t_start

1543

1544 if seg_dt <= EPSILON_TIME:

1545 continue

1546

1547 if seg["type"] == "constant":

1548 integral_c_dt += seg["concentration"] * seg_dt

1549 elif seg["type"] == "rarefaction":

1550 if isinstance(sorption, FreundlichSorption):

1551 raref = seg["wave"]

1552 integral_c_dt += integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)

1553 else:

1554 # ConstantRetardation - use midpoint

1555 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)

1556 integral_c_dt += c_mid * seg_dt

1557

1558 mass_up_to_raref_start = mass_within_tedges + flow_extended * integral_c_dt

1559 flow_at_raref_start = flow_extended

1560

1561 # Find rarefactions that are active at the outlet after t_last_raref_start

1562 # and add their total integrated mass

1563 total_raref_mass = 0.0

1564

1565 for wave in waves:

1566 if not wave.is_active:

1567 continue

1568

1569 if isinstance(wave, RarefactionWave):

1570 # Check if this rarefaction reaches the outlet

1571 head_vel = wave.head_velocity()

1572 if head_vel > EPSILON_VELOCITY and wave.v_start < v_outlet:

1573 t_raref_start_at_outlet = wave.t_start + (v_outlet - wave.v_start) / head_vel

1574

1575 # If this rarefaction starts at or after t_last_raref_start, include its total mass

1576 # (with small tolerance for numerical precision)

1577 if abs(t_raref_start_at_outlet - t_last_raref_start) < EPSILON_TIME_MATCH:

1578 # This is the last rarefaction - integrate to infinity

1579 raref_mass = integrate_rarefaction_total_mass(

1580 raref=wave,

1581 v_outlet=v_outlet,

1582 t_start=t_raref_start_at_outlet,

1583 sorption=sorption,

1584 flow=flow_at_raref_start,

1585 )

1586 total_raref_mass += raref_mass

1587

1588 # For rarefactions with finite tails (c_tail > 0, typical for n < 1),

1589 # all mass is already accounted for in the rarefaction integration

1590 # from head to tail. No additional mass needs to be integrated after

1591 # the tail - if there were more waves, they would be in the wave list.

1592 total_mass = mass_up_to_raref_start + total_raref_mass

1593

1594 return float(total_mass), t_last_raref_start

Coverage for src / gwtransport / fronttracking / output.py: 89%

417 statements