fix(prt): rework vertical tracking approach (#2066)

Address #2014. This fixes a crash, avoids a potential infinite loop condition, and reworks vertical tracking behavior more generally. The aim is to produce results matching MF6.5 and MP7 by default (barring a few edge cases), while giving more control over how particles behave in dry conditions. Also some miscellaneous reorganization and cleanup.

See the included dev notes document and/or mf6io for details on this change.

autotest/test_prt_disv1.py → autotest/test_prt_newton_disv1.py

@@ -6,13 +6,7 @@
 refined cells, instead of the new ternary method
 which applies more generally to polygonal cells.
 
-The simulation includes a single stress period
-with multiple time steps. This serves to test
-whether PRT properly solves trajectories over
-"internal" time steps, i.e. within the step's
-slice of simulation time, as well as extending
-tracking to termination or particle stop times
-during the simulation's final time step.
+The simulation uses the Newton formulation.
 
 Several cases are provided:
     - default: No user-specified tracking times, MP7 in pathline mode.

autotest/test_prt_voronoi1.py

@@ -10,7 +10,6 @@
 """
 
 from pathlib import Path
-from platform import processor, system
 
 import flopy
 import matplotlib as mpl
@@ -267,6 +266,7 @@ def plot_output(idx, test):
     ax = plt.subplot(1, 1, 1, aspect="equal")
     pmv = flopy.plot.PlotMapView(model=gwf, ax=ax)
     pmv.plot_grid(alpha=0.25)
+
     pmv.plot_ibound(alpha=0.5)
     headmesh = pmv.plot_array(head, alpha=0.25)
     cv = pmv.contour_array(head, levels=np.linspace(0, 1, 9), colors="black")
@@ -331,7 +331,7 @@ def plot_output(idx, test):
     plt.savefig(prt_ws / f"{name}.png")
 
 
-def check_output(idx, test, snapshot):
+def check_output(idx, test):
     name = test.name
     prt_ws = test.workspace / "prt"
     prt_name = get_model_name(name, "prt")
@@ -349,22 +349,20 @@ def check_output(idx, test, snapshot):
     pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False)
     endpts = pls[pls.ireason == 3]  # termination
 
-    # compare pathlines with snapshot. particles shouldn't
-    # have moved vertically. round for cross-platform error.
-    # skip macos-14 in CI because grid is slightly different
-    if not (system() == "Darwin" and processor() == "arm"):
-        assert snapshot == endpts.drop("name", axis=1).round(1).to_records(index=False)
+    assert np.allclose(endpts.z, 0.5)
+    assert np.isclose(endpts.y.min(), 1, atol=4)
+    assert np.isclose(endpts.y.max(), 996, atol=4)
 
 
 @requires_pkg("syrupy")
 @pytest.mark.slow
 @pytest.mark.parametrize("idx, name", enumerate(cases))
-def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot, plot):
+def test_mf6model(idx, name, function_tmpdir, targets, benchmark, plot):
     test = TestFramework(
         name=name,
         workspace=function_tmpdir,
         build=lambda t: build_models(idx, t),
-        check=lambda t: check_output(idx, t, array_snapshot),
+        check=lambda t: check_output(idx, t),
         plot=lambda t: plot_output(idx, t) if plot else None,
         targets=targets,
         compare=None,

doc/ReleaseNotes/develop.tex

@@ -51,8 +51,9 @@
 		\item The number of characters used to represent integers and floating point numbers in MODFLOW input files was restricted to 30.  The program was modified to accept any number of characters provided the number is valid.  This may be useful for parameter estimation programs that use character substitution to create new input files.
 		\item A string to character array conversion function in the BMI interface could fail on Apple silicon macOS with recent versions of GNU Fortran, producing array extent errors at runtime. This has been fixed by properly specifying the intent of a dummy argument in the relevant function.
 		\item The UZF Package will facilitate UZF objects with areas less than the area of the host cell.  However, within the GWE model type, the UZE package will not work properly when the area of UZF objects is not equal to the area of the host cell.  New code was added to ensure that the area of each UZF object is equal to that of the host grid cell.  When this condition is not satisfied, the new code will stop the simulation with an error message indicating which cell is in violation.
-	%	\item xxx
-
+		\item With a flow model using the Newton formulation, the PRT model could crash upon a particle's entry into a dry cell. This has been fixed.
+		\item With a flow model using the Newton formulation, the PRT model could enter an endless loop upon a particle's entry to a dry cell if that cell contains a boundary package (e.g. a pumping well). Where the particle should be captured and terminate, it would instead be passed back and forth between the cell bottom and the top of the cell below. To avoid this, particles are forbidden from backtracking (reentering the previous cell) within the same time step.
+		\item The PRT model now allows more control over vertical particle motion in dry conditions. In addition to the existing DRAPE option, which controls release-time behavior, the PRP package now provides a DRY\_TRACKING\_METHOD option which configures how dry particles (particles in dry cells, or above the water table in partially saturated cells) behave at tracking time. This option is relevant only when the Newton formulation is used, in which case dry cells remain active; otherwise, dry cells are inactive and particles will terminate. See the MF6IO document for a detailed explanation of DRY\_TRACKING\_METHOD.
 	\end{itemize}
 
 	%\underline{INTERNAL FLOW PACKAGES}

doc/mf6io/mf6ivar/dfn/prt-prp.dfn

@@ -167,6 +167,15 @@ optional true
 longname drape
 description is a text keyword to indicate that if a particle's release point is in a cell that happens to be inactive at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is inactive at release time.
 
+block options
+name dry_tracking_method
+type string
+valid drop stop stay
+reader urword
+optional true
+longname what to do in dry-but-active cells
+description is a string indicating how particles should behave in dry-but-active cells (as can occur with the Newton formulation).  The value can be ``DROP'', ``STOP'', or ``STAY''.  The default is ``DROP'', which passes particles vertically and instantaneously to the water table. ``STOP'' causes particles to terminate. ``STAY'' causes particles to remain stationary but active.
+
 block options
 name dev_forceternary
 type keyword

doc/mf6io/prt/prp.tex

@@ -2,6 +2,8 @@
 
 The PRP Package offers multiple ways to specify particle release times.  Particle release times may either be provided explicitly, relative to the simulation start time, or configured relative to the time discretization of stress periods via period block settings.  When multiple ways of specifying release times are used together, the resulting set of release times is the union of the times specified by each method.
 
+This package contains options to configure how particles behave in dry conditions. Each of these is described in detail below. Their interaction is summarized in the PRT chapter introduction.
+
 \vspace{5mm}
 \subsubsection{Structure of Blocks}
 \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/prt-prp-options.dat}

doc/mf6io/prt/prt.tex

@@ -21,6 +21,42 @@ \subsection{Specifying Cell Face Flows using IFLOWFACE}
 \subsection{Particle Mass Budget}
 A summary of all inflow (sources) and outflow (sinks) of particle mass is called a mass budget.  The particle mass budget is printed to the PRT Model Listing File for selected time steps.  In the current implementation, each particle is assigned unit mass, and the numerical value of the flow can be interpreted as particles per time.
 
+\subsection{Vertical Tracking}
+
+When a particle is in the flow field, vertical motion can be solved in the same way as lateral motion. Special handling is necessary above the water table.
+
+A dry cell is either an inactive cell or an active but dry cell, as can occur with the Newton formulation.
+
+Normally, an inactive cell might be dry or explicitly disabled (idomain).  With Newton, dry cells remain active.
+
+Release-time and tracking-time considerations are described (and implemented) separately.
+
+\subsubsection{Release}
+
+At release time, PRT decides whether to release each particle or to terminate it unreleased.
+
+If the release cell is active, the particle will be released at the user-specified location.
+
+If the release cell is inactive, behavior is determined by the DRAPE option. If the DRAPE option is enabled, the particle will be released from the top-most active cell beneath it, if any. If there is no active cell underneath the particle in any layer, or if DRAPE is not enabled, the particle will terminate unreleased (with status code 8).
+
+Since under the Newton formulation dry cells can remain active, the DRAPE option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured using the DRY\_TRACKING\_METHOD option discussed below.
+
+\subsubsection{Tracking}
+
+With the Newton formulation, particles can be released into dry-but-active cells.
+
+Particle trajectories are solved over the same time discretization used by the flow model. A particle may be immersed in the flow field in one time step, and find that the water table has dropped below it in the next.
+
+A particle which finds itself in an inactive cell will terminate with status code 7. This is consistent with MODPATH 7's behavior.
+
+A particle in a dry but active cell, or above the water table in a partially saturated cell, need not terminate. We call such a particle dry. The PRP package provides an option DRY\_TRACKING\_METHOD determining how dry particles should behave. Supported values are DROP (default), STOP, and STAY.
+
+If DROP is selected, or if a DRY\_TRACKING\_METHOD is unspecified, a dry particle is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. If the vertical column containing the particle is entirely dry, the particle will terminate upon reaching the bottom of the model grid.
+
+If STOP is selected, dry particles will be terminated.
+
+If STAY is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue, b) the particle terminates due to reaching its STOPTIME or STOPTRAVELTIME, or c) the simulation ends.
+
 \subsection{Particle Track Output}
 
 The PRT Model supports both binary and CSV particle track output files. A particle track CSV file contains the output data in tabular format. The first line of the CSV file contains column names. Each subsequent line in the file contains a row of data for a single particle track record, with the following fields:

src/Model/ModelUtilities/TrackFile.f90

@@ -35,7 +35,7 @@ module TrackFileModule
   !!     1: particle transitioned between cells
   !!     2: current time step ended****
   !!     3: particle terminated
-  !!     4: particle exited weak sink
+  !!     4: particle in weak sink
   !!     5: user-specified tracking time
   !!
   !! Each record has an "istatus" property, which is the tracking status;

src/Model/ParticleTracking/prt-prp.f90

@@ -49,6 +49,7 @@ module PrtPrpModule
     integer(I4B), pointer :: istopweaksink => null() !< weak sink option: 0 = no stop, 1 = stop
     integer(I4B), pointer :: istopzone => null() !< optional stop zone number: 0 = no stop zone
     integer(I4B), pointer :: idrape => null() !< drape option: 0 = do not drape, 1 = drape to topmost active cell
+    integer(I4B), pointer :: idrymeth => null() !< dry tracking method: 0 = drop, 1 = stop, 2 = stay
     integer(I4B), pointer :: itrkout => null() !< binary track file
     integer(I4B), pointer :: itrkhdr => null() !< track header file
     integer(I4B), pointer :: itrkcsv => null() !< CSV track file
@@ -89,6 +90,7 @@ module PrtPrpModule
     procedure :: release
     procedure :: log_release
     procedure :: validate_release_point
+    procedure :: initialize_particle
     procedure, public :: bnd_obs_supported => prp_obs_supported
     procedure, public :: bnd_df_obs => prp_df_obs
   end type PrtPrpType
@@ -158,6 +160,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%istopweaksink)
     call mem_deallocate(this%istopzone)
     call mem_deallocate(this%idrape)
+    call mem_deallocate(this%idrymeth)
     call mem_deallocate(this%nreleasepoints)
     call mem_deallocate(this%nreleasetimes)
     call mem_deallocate(this%nparticles)
@@ -247,6 +250,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
     call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
     call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath)
+    call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath)
     call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
     call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
     call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
@@ -270,6 +274,7 @@ subroutine prp_allocate_scalars(this)
     this%istopweaksink = 0
     this%istopzone = 0
     this%idrape = 0
+    this%idrymeth = 0
     this%nreleasepoints = 0
     this%nreleasetimes = 0
     this%nparticles = 0
@@ -413,47 +418,52 @@ subroutine release(this, ip, trelease)
     integer(I4B), intent(in) :: ip !< particle index
     real(DP), intent(in) :: trelease !< release time
     ! local
-    integer(I4B) :: irow, icol, ilay, icpl
-    integer(I4B) :: ic, icu
     integer(I4B) :: np
-    real(DP) :: x, y, z
-    real(DP) :: top, bot, hds
     type(ParticleType), pointer :: particle
 
+    call this%initialize_particle(particle, ip, trelease)
+
     ! Increment cumulative particle count
     np = this%nparticles + 1
     this%nparticles = np
 
-    ! Get reduced and user node numbers
-    ic = this%rptnode(ip)
-    icu = this%dis%get_nodeuser(ic)
+    ! Save the particle to the store
+    call this%particles%save_particle(particle, np)
+    deallocate (particle)
 
-    ! Load coordinates and transform if needed
-    x = this%rptx(ip)
-    y = this%rpty(ip)
-    if (this%ilocalz > 0) then
-      top = this%fmi%dis%top(ic)
-      bot = this%fmi%dis%bot(ic)
-      hds = this%fmi%gwfhead(ic)
-      z = bot + this%rptz(ip) * (hds - bot)
-    else
-      z = this%rptz(ip)
-    end if
+    ! Accumulate mass for release point
+    this%rptm(ip) = this%rptm(ip) + DONE
 
-    ! Make sure release point is in the grid/cell
-    call this%validate_release_point(ic, x, y, z)
+  end subroutine release
+
+  subroutine initialize_particle(this, particle, ip, trelease)
+    class(PrtPrpType), intent(inout) :: this !< this instance
+    type(ParticleType), pointer, intent(inout) :: particle !< the particle
+    integer(I4B), intent(in) :: ip !< particle index
+    real(DP), intent(in) :: trelease !< release time
+    ! local
+    integer(I4B) :: irow, icol, ilay, icpl
+    integer(I4B) :: ic, icu, ic_old
+    real(DP) :: x, y, z
+    real(DP) :: top, bot, hds
+
+    ic = this%rptnode(ip)
+    icu = this%dis%get_nodeuser(ic)
 
-    ! Initialize the particle
     call create_particle(particle)
+
     if (size(this%boundname) /= 0) then
       particle%name = this%boundname(ip)
     else
       particle%name = ''
     end if
+
     particle%irpt = ip
     particle%istopweaksink = this%istopweaksink
     particle%istopzone = this%istopzone
+    particle%idrymeth = this%idrymeth
     particle%icu = icu
+
     select type (dis => this%dis)
     type is (DisType)
       call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
@@ -462,28 +472,49 @@ subroutine release(this, ip, trelease)
     end select
     particle%ilay = ilay
     particle%izone = this%rptzone(ic)
+
     particle%istatus = 0
-    ! Handle inactive cells
+    ! If the cell is inactive, either drape the particle
+    ! to the top-most active cell beneath it if drape is
+    ! enabled, or else terminate permanently unreleased.
     if (this%ibound(ic) == 0) then
-      ! If drape option activated, release in highest active
-      ! cell vertically below release point.
-      if (this%idrape /= 0) &
+      ic_old = ic
+      if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-      ! If returned cell is inactive, do not release particle
-      if (this%ibound(ic) == 0) &
-        particle%istatus = 8 ! permanently unreleased
+        if (ic == ic_old .or. this%ibound(ic) == 0) &
+          particle%istatus = 8
+      else
+        particle%istatus = 8
+      end if
+    end if
+
+    ! Load coordinates and transform if needed
+    x = this%rptx(ip)
+    y = this%rpty(ip)
+    if (this%ilocalz > 0) then
+      top = this%fmi%dis%top(ic)
+      bot = this%fmi%dis%bot(ic)
+      hds = this%fmi%gwfhead(ic)
+      z = bot + this%rptz(ip) * (hds - bot)
+    else
+      z = this%rptz(ip)
     end if
+
+    call this%validate_release_point(ic, x, y, z)
+
     particle%x = x
     particle%y = y
     particle%z = z
     particle%trelease = trelease
-    ! Set stopping time to earlier of times specified by STOPTIME and STOPTRAVELTIME
+
+    ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME
     if (this%stoptraveltime == huge(1d0)) then
       particle%tstop = this%stoptime
     else
       particle%tstop = particle%trelease + this%stoptraveltime
       if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
     end if
+
     particle%ttrack = particle%trelease
     particle%idomain(1) = 0
     particle%iboundary(1) = 0
@@ -495,15 +526,7 @@ subroutine release(this, ip, trelease)
     particle%iexmeth = this%iexmeth
     particle%iextend = this%iextend
     particle%extol = this%extol
-
-    ! Store the particle's state and deallocate it
-    call this%particles%save_particle(particle, np)
-    deallocate (particle)
-
-    ! Accumulate particle mass for release point
-    this%rptm(ip) = this%rptm(ip) + DONE
-
-  end subroutine release
+  end subroutine initialize_particle
 
   !> @ brief Read and prepare period data for particle input
   subroutine prp_rp(this)
@@ -691,6 +714,25 @@ subroutine prp_options(this, option, found)
     case ('DRAPE')
       this%idrape = 1
       found = .true.
+    case ('DRY_TRACKING_METHOD')
+      call this%parser%GetStringCaps(keyword)
+      select case (keyword)
+      case ('DROP')
+        this%idrymeth = 0
+      case ('STOP')
+        this%idrymeth = 1
+      case ('STAY')
+        this%idrymeth = 2
+      case default
+        write (errmsg, '(a, a)') &
+          'Unknown dry tracking method: ', trim(keyword)
+        call store_error(errmsg)
+        write (errmsg, '(a, a)') &
+          'DRY must be "DROP", "STOP" or "STAY"'
+        call store_error(errmsg)
+        call this%parser%StoreErrorUnit()
+      end select
+      found = .true.
     case ('TRACK')
       call this%parser%GetStringCaps(keyword)
       if (keyword == 'FILEOUT') then
@@ -891,7 +933,7 @@ subroutine prp_read_packagedata(this)
 
         if (n < 1 .or. n > this%nreleasepoints) then
           write (errmsg, '(a,1x,i0,a)') &
-            'Release point number must be greater than 0 and less than ', &
+            'Release point number must be greater than 0 and less than', &
             'or equal to', this%nreleasepoints, '.'
           call store_error(errmsg)
           cycle