fix(prt): fix drape behavior (MODFLOW-ORG#2824)

Fix the PRT PRP package's DRAPE option's behavior. A draped particle's release was previously reported in the original release position, followed by a "dropped" (ireason=6) event due to the particle reaching the water table at tracking time as a result of the default DRY_TRACKING_METHOD setting DROP. This is inconsistent with the intent of the DRAPE option: DRAPE is applied before release to control a particle's release position. Draped particles will now be released in the highest active cell.

While this fix produces essentially identical pathlines (ignoring the release event and subsequent dropped event) as the previous logic, it is worth pointing out that PRT's DRAPE semantics diverge from MP7 in an important way: MP7 unconditionally calculates a particle release position's z coordinate with the given local z value, regardless whether the particle was draped. This means a draped particle may start beneath the water table if local z < 1. Whereas, PRT did previously and will continue to start draped particles at the water table in convertible cells. Both PRT and MP7 start draped particles at the geometric top of confined cells.

Author: wpbonelli

autotest/__snapshots__/test_prt_dry/test_mf6model[1-prtdry_drape].npy

@@ -276,7 +276,8 @@ def check_output(idx, test):
         )
 
     if drape:
-        assert mf6_pls.shape[0] == 36
+        assert mf6_pls.shape[0] == 27
+        assert mf6_pls[mf6_pls.ireason == 6].empty
     else:
         # expect no movement without drape
         assert mf6_pls.shape[0] == 9

autotest/test_prt_drape.py

@@ -68,3 +68,8 @@ description = "Fix an indexing issue in structured NetCDF export where READASARR
 section = "fixes"
 subsection = "basic"
 description = "Stop with an informative error if a discretization package is missing after loading the model input file's PACKAGES block."
+
+[[items]]
+section = "fixes"
+subsection = "model"
+description = "Fix the PRT PRP package's DRAPE option's behavior. A draped particle's release was previously reported in the original release position, followed by a DROPPED (ireason=6) event due to the particle reaching the water table at tracking time as a result of the default DRY\\_TRACKING\\_METHOD setting DROP. This is inconsistent with the intent of the DRAPE option: DRAPE is applied before release to control a particle's release position. Draped particles will now be released in the highest active cell, at the water table if the cell is convertible or at the geometric top if the cell is confined, and no DROPPED event will be reported."

doc/ReleaseNotes/develop.toml

@@ -605,6 +605,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     integer(I4B), intent(in) :: ip !< particle index
     real(DP), intent(in) :: trelease !< release time
     ! local
+    logical(LGP) :: draped
     integer(I4B) :: irow, icol, ilay, icpl
     integer(I4B) :: ic, icu, ic_old
     real(DP) :: x, y, z
@@ -617,7 +618,6 @@ subroutine initialize_particle(this, particle, ip, trelease)
       &Make sure a GWFGRID entry is configured in the PRT FMI package.')"
 
     ic = this%rptnode(ip)
-    icu = this%dis%get_nodeuser(ic)
 
     call create_particle(particle)
 
@@ -631,25 +631,18 @@ subroutine initialize_particle(this, particle, ip, trelease)
     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)
-    type is (DisvType)
-      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
-    end select
-    particle%ilay = ilay
-    particle%izone = this%rptzone(ic)
     particle%istatus = 0 ! status 0 until tracking starts
+
     ! 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.
+    draped = .false.
     if (this%ibound(ic) == 0) then
       ic_old = ic
       if (this%drape) then
         call this%dis%highest_active(ic, this%ibound)
-        if (ic == ic_old .or. this%ibound(ic) == 0) then
+        draped = ic /= ic_old
+        if (.not. draped .and. this%ibound(ic) == 0) then
           ! negative unreleased status signals to the
           ! tracking method that we haven't yet saved
           ! a termination record, it needs to do so.
@@ -660,24 +653,35 @@ subroutine initialize_particle(this, particle, ip, trelease)
       end if
     end if
 
-    ! load coordinates
-    x = this%rptx(ip)
-    y = this%rpty(ip)
-    if (this%localz) then
-      ! make sure FMI has cell type array. we need
-      ! it to distinguish convertible and confined
-      ! cells if release z coordinates are local
-      if (this%fmi%igwfceltyp /= 1) then
-        write (errmsg, fmticterr) trim(this%text)
-        call store_error(errmsg, terminate=.TRUE.)
-      end if
+    icu = this%dis%get_nodeuser(ic)
+    particle%icu = icu
+    select type (dis => this%dis)
+    type is (DisType)
+      call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
+    type is (DisvType)
+      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
+    end select
+    particle%ilay = ilay
+    particle%izone = this%rptzone(ic)
+
+    ! if the particle was draped, override the release z coord and
+    ! set it to the saturated top of the cell. this puts a draped
+    ! a draped particle at the water table for a convertible cell
+    ! or at the geometric cell top for a confined cell. if it was
+    ! not draped and localz is enabled, calculate a model z coord
+    ! using the geometric cell top if the cell is confined or the
+    ! water table as the effective top if the cell is convertible.
+    if (draped) then
+      z = this%fmi%dis%bot(ic) + &
+          this%fmi%gwfsat(ic) * &
+          (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
+    else if (this%localz) then
+      ! TODO: is this sufficient instead of the below??
+      ! z = this%fmi%dis%bot(ic) + &
+      !     this%rptz(ip) * &
+      !     this%fmi%gwfsat(ic) * &
+      !     (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
 
-      ! calculate model z coord from local z coord.
-      ! if cell is confined (icelltype == 0) use the
-      ! actual cell height (geometric top - bottom).
-      ! otherwise use head as cell top, clamping to
-      ! the cell bottom if head is below the bottom
-      ! and to geometric cell top if head is above top
       top = this%fmi%dis%top(ic)
       bot = this%fmi%dis%bot(ic)
       if (this%fmi%gwfceltyp(icu) /= 0) then
@@ -690,6 +694,9 @@ subroutine initialize_particle(this, particle, ip, trelease)
       z = this%rptz(ip)
     end if
 
+    x = this%rptx(ip)
+    y = this%rpty(ip)
+
     if (this%ichkmeth > 0) &
       call this%validate_release_point(ic, x, y, z)