Feature Request: Add a highly visible warning when users specify IFACE values that will be ignored by MP or PRT (ie. when iface is not along perimeter of model).

[pfjuckem-usgs]

Abrams et al. (2012) document proper methods for performing particle tracking with weak-sink cells associated with surface water boundary packages. Specifically, when forward tracking, particles should be set to "stop at weak sinks" (WeakSinkOption=2 and IFACE should be set to 6 (top of the cell). This works properly with the provisional MP7 code (mp7_win64_20231016_86b38df.exe) as long as the boundary package cell is along the boundary of the model domain (for IFACE=6, or surface water packages, this means being placed in layer 1 where the top of the model is the boundary of the model domain; for CHD or GHB or other packages where IFACE <6 may be desired, the IFACE designation must align with the perimeter of the domain. For example, MP7 would ignore an IFACE=1 setting in a GHB cell along the right side of the model domain). This is documented on page 19 of the MP7 manual (Pollock, 2016), so the issue we're submitting does not identify a bug. However, this highly nuanced characteristic is not abundantly clear in the documentation, and we anticipate that very few users fully understand and appreciate these nuances. Indeed, we have encountered several MODFLOW models in which users have assigned surface water boundary cells to layers > 1 and have NOT inactivated the cells above. This may not have been deemed problematic for their flow-orientated objectives with their MODFLOW models, but it creates a large problem when using the model with MODPATH for particle tracking. Thus, we request that a very eye-catching warning be printed to the DOS prompt AND added to the MP7 output (MPLIST) file in a very visible location (the bottom of the file) whenever an IFACE value entered by a user does not match with the perimeter of the domain for that cell. We also encourage you to include greater discussion of this topic in future documentation (if not for MODPATH, then the PRT package for MF6).

Concerns about how this issue can influence results are illustrated below and with the files in the attached ZIP file. Example ex01b_mf2005 from Pollock (2016) was modified by (A) setting WeakSinkOption=2, and (B) adding two RIV cells into layer 2. Then MP7 was run for two instances. In the first instance (ex01b_mf2005_addriv.mpsim), IBOUND remained unchanged. In the second instance (ex01b_mf2005_rivibound.mpsim), IBOUND was set to zero in the cells above these new RIV cells. Both instances were evaluated using the "ept_eval.py" code that's included in the ZIP file (along with the associated MP7 input files; MF2005 would need to be re-run to generate the binary head and flux files required by MP7 [or I can provide them if requested]).

The images illustrate the forward-tracked particle starting locations, with colors used to identify their final end point. Red particles discharge to the well. Blues, greens, and yellows correspond with the RIV cells, with the individual color matching the unique "node" number for each RIV cell. Yellow cells and particles are associated with the RIV cells in layer 2, so those are the particles of interest for this issue. The first image is for the instance when all cells are active ("addriv"). The second image is for the instance when cells above the RIV cells have been deactivated ("rivibound"), and represents the correct solution. The images show a smaller contributing area for the correct instance where overlying cells have been deactivated ("rivibound"). (Note that particles started in the inactive cells are not tracked by MP7, but they don't discharge to the layer-2 RIV cells anyway, so they're irrelevant, even though it's an eye-catching difference in the plots).

The ZIP file contains ASCII files ("...summary.txt") that summarize important results from the endpoint files. Of note is that the first instance with active cells above RIV cells ("addriv") has 77 particles that stop at RIV cells in layer 2, but almost half stop at IFACE=5 (cell bottom) and the average age of those particles is nearly 3x that of the particles that stop at IFACE=6. While it's not included in the summary, investigation of the *.endpoint7 file shows that both RIV cells have a mix of particles that terminate at IFACE 5 & 6 with Status values of 3 & 5. In comparison, for the instance where overlying cells were deactivated ("rivibound"), only 36 particles were captured by the layer-2 RIV cells (more were allowed to pass-through) because only particles that intersected IFACE=6 were stopped. Moreover, all particles had a status=2 (terminated at boundary face). Equally important is that the average age for instance 2 ("rivibound") is much older than the IFACE=6 cells from the first instance ("addriv" simulation), indicating that one can't simply post-process the endpoint file to select only IFACE=6 particles; MP7 handles the actual tracking differently depending upon the IBOUND status above the layer-2 RIV cells!

ex01b_mf2005_addriver.zip

Thank you for considering this request, and please reach out to Leon Kauffman or I with any questions or comments. Thanks! Paul

References:
Abrams, D., Haitjema, H., and Kauffman, L., 2012. “On Modeling Weak Sinks in MODPATH.” Groundwater 51 (4): 597–602. https://doi.org/10.1111/j.1745-6584.2012.00995.x.
Pollock, D. W., 2016. “User Guide for MODPATH Version 7 -- A Particle-Tracking Model for MODFLOW: U.S. Geological Survey Open-File Report.” U.S. Geological Survey. https://doi.org/https://doi.org/10.3133/ofr20161086.

[aprovost-usgs]

@pfjuckem-usgs Thanks for this suggestion, which follows on a discussion we had earlier. I very much appreciate your thorough explanation and documentation of the issue, which will help us understand the ramifications and decide how best to proceed. As we discussed, there are both practical and philosophical considerations as to when/how/where to issue warnings. We'll give your suggestion careful consideration in the context of both MP7 and MF6-PRT.