diff --git a/doc/mf6io/prt/prt.tex b/doc/mf6io/prt/prt.tex
index 6ead0ffc30d..39ed9bfb0a7 100644
--- a/doc/mf6io/prt/prt.tex
+++ b/doc/mf6io/prt/prt.tex
@@ -115,6 +115,7 @@ \subsection{Particle Track Output}
 \item \texttt{7}: particle terminated in an inactive cell
 \item \texttt{8}: particle terminated immediately upon release into a dry cell
 \item \texttt{9}: particle terminated in a subcell with no exit face
+\item \texttt{10}: particle terminated at stop time or end of simulation
 \end{description}
 
 The IREASON field indicates the reason the particle track record was saved:
diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90
index d0d23e5ab71..2ef386b86e0 100644
--- a/src/Model/ParticleTracking/prt-prp.f90
+++ b/src/Model/ParticleTracking/prt-prp.f90
@@ -445,7 +445,7 @@ subroutine release(this, ip, trelease)
   end subroutine release
 
   subroutine initialize_particle(this, particle, ip, trelease)
-    use ParticleModule, only: TERM_UNRELEASED
+    use ParticleModule, only: TERM_UNRELEASED, LVL_MODEL, LVL_CELL, LVL_SUBCELL
     class(PrtPrpType), intent(inout) :: this !< this instance
     type(ParticleType), pointer, intent(inout) :: particle !< the particle
     integer(I4B), intent(in) :: ip !< particle index
@@ -528,12 +528,12 @@ subroutine initialize_particle(this, particle, ip, trelease)
     end if
 
     particle%ttrack = particle%trelease
-    particle%idomain(1) = 0
-    particle%iboundary(1) = 0
-    particle%idomain(2) = ic
-    particle%iboundary(2) = 0
-    particle%idomain(3) = 0
-    particle%iboundary(3) = 0
+    particle%itrdomain(LVL_MODEL) = 0
+    particle%iboundary(LVL_MODEL) = 0
+    particle%itrdomain(LVL_CELL) = ic
+    particle%iboundary(LVL_CELL) = 0
+    particle%itrdomain(LVL_SUBCELL) = 0
+    particle%iboundary(LVL_SUBCELL) = 0
     particle%ifrctrn = this%ifrctrn
     particle%iexmeth = this%iexmeth
     particle%iextend = this%iextend
diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90
index 9dfe1036f6a..ef87ff376b1 100644
--- a/src/Model/ParticleTracking/prt.f90
+++ b/src/Model/ParticleTracking/prt.f90
@@ -405,6 +405,7 @@ end subroutine prt_cq
 
   !> @brief Calculate particle mass storage
   subroutine prt_cq_sto(this)
+    use ParticleModule, only: LVL_CELL
     ! modules
     use TdisModule, only: delt
     use PrtPrpModule, only: PrtPrpType
@@ -436,7 +437,7 @@ subroutine prt_cq_sto(this)
           istatus = packobj%particles%istatus(np)
           ! this may need to change if istatus flags change
           if ((istatus > 0) .and. (istatus /= 8)) then
-            n = packobj%particles%idomain(np, 2)
+            n = packobj%particles%itrdomain(np, LVL_CELL)
             ! Each particle currently assigned unit mass
             this%masssto(n) = this%masssto(n) + DONE
           end if
diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
index 82564f80165..4318faf8c48 100644
--- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90
+++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
@@ -43,7 +43,7 @@ end subroutine deallocate
 
   !> @brief Pass particle vertically and instantaneously to the cell bottom
   subroutine apply_ptb(this, particle, tmax)
-    use ParticleModule, only: TERM_NO_EXITS
+    use ParticleModule, only: TERM_NO_EXITS, LVL_CELL
     ! dummy
     class(MethodCellPassToBotType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -57,7 +57,7 @@ subroutine apply_ptb(this, particle, tmax)
 
     ! Pass to bottom face
     particle%z = this%cell%defn%bot
-    particle%iboundary(2) = this%cell%defn%npolyverts + 2
+    particle%iboundary(LVL_CELL) = this%cell%defn%npolyverts + 2
 
     ! Terminate if in bottom layer
     select type (dis => this%fmi%dis)
diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90
index 06790e8f7f5..b90e915b5ba 100644
--- a/src/Solution/ParticleTracker/MethodCellPollock.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollock.f90
@@ -68,12 +68,13 @@ subroutine load_mcp(this, particle, next_level, submethod)
       trackctl=this%trackctl, &
       tracktimes=this%tracktimes)
     submethod => method_subcell_plck
-    particle%idomain(next_level) = 1
+    particle%itrdomain(next_level) = 1
   end subroutine load_mcp
 
   !> @brief Having exited the lone subcell, pass the particle to the cell face
   !! In this case the lone subcell is the cell.
   subroutine pass_mcp(this, particle)
+    use ParticleModule, only: LVL_CELL, LVL_SUBCELL
     ! dummy
     class(MethodCellPollockType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -81,7 +82,7 @@ subroutine pass_mcp(this, particle)
     integer(I4B) :: exitface
     integer(I4B) :: entryface
 
-    exitface = particle%iboundary(3)
+    exitface = particle%iboundary(LVL_SUBCELL)
     ! Map subcell exit face to cell face
     select case (exitface) ! note: exitFace uses Dave's iface convention
     case (0)
@@ -100,7 +101,7 @@ subroutine pass_mcp(this, particle)
       entryface = 7
     end select
     if (entryface .eq. -1) then
-      particle%iboundary(2) = 0
+      particle%iboundary(LVL_CELL) = 0
     else
       if ((entryface .ge. 1) .and. (entryface .le. 4)) then
         ! Account for local cell rotation
@@ -110,7 +111,7 @@ subroutine pass_mcp(this, particle)
         end select
         if (entryface .gt. 4) entryface = entryface - 4
       end if
-      particle%iboundary(2) = entryface
+      particle%iboundary(LVL_CELL) = entryface
     end if
   end subroutine pass_mcp
 
diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
index 5bd749e192d..85c1e0c9c32 100644
--- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
@@ -71,6 +71,7 @@ end subroutine load_mcpq
 
   !> @brief Pass particle to next subcell if there is one, or to the cell face
   subroutine pass_mcpq(this, particle)
+    use ParticleModule, only: LVL_CELL, LVL_SUBCELL
     ! dummy
     class(MethodCellPollockQuadType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -79,8 +80,8 @@ subroutine pass_mcpq(this, particle)
 
     select type (cell => this%cell)
     type is (CellRectQuadType)
-      exitFace = particle%iboundary(3)
-      isc = particle%idomain(3)
+      exitFace = particle%iboundary(LVL_SUBCELL)
+      isc = particle%itrdomain(LVL_SUBCELL)
       npolyverts = cell%defn%npolyverts
 
       ! exitFace uses MODPATH 7 iface convention here
@@ -92,13 +93,13 @@ subroutine pass_mcpq(this, particle)
         select case (isc)
         case (1)
           ! W face, subcell 1 --> E face, subcell 4  (cell interior)
-          particle%idomain(3) = 4
-          particle%iboundary(3) = 2
+          particle%itrdomain(LVL_SUBCELL) = 4
+          particle%iboundary(LVL_SUBCELL) = 2
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (2)
           ! W face, subcell 2 --> E face, subcell 3 (cell interior)
-          particle%idomain(3) = 3
-          particle%iboundary(3) = 2
+          particle%itrdomain(LVL_SUBCELL) = 3
+          particle%iboundary(LVL_SUBCELL) = 2
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (3)
           ! W face, subcell 3 (cell face)
@@ -121,21 +122,21 @@ subroutine pass_mcpq(this, particle)
           infaceoff = -1
         case (3)
           ! E face, subcell 3 --> W face, subcell 2 (cell interior)
-          particle%idomain(3) = 2
-          particle%iboundary(3) = 1
+          particle%itrdomain(LVL_SUBCELL) = 2
+          particle%iboundary(LVL_SUBCELL) = 1
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (4)
           ! E face, subcell 4 --> W face subcell 1 (cell interior)
-          particle%idomain(3) = 1
-          particle%iboundary(3) = 1
+          particle%itrdomain(LVL_SUBCELL) = 1
+          particle%iboundary(LVL_SUBCELL) = 1
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         end select
       case (3)
         select case (isc)
         case (1)
           ! S face, subcell 1 --> N face, subcell 2 (cell interior)
-          particle%idomain(3) = 2
-          particle%iboundary(3) = 4
+          particle%itrdomain(LVL_SUBCELL) = 2
+          particle%iboundary(LVL_SUBCELL) = 4
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (2)
           ! S face, subcell 2 (cell face)
@@ -147,8 +148,8 @@ subroutine pass_mcpq(this, particle)
           infaceoff = -1
         case (4)
           ! S face, subcell 4 --> N face, subcell 3 (cell interior)
-          particle%idomain(3) = 3
-          particle%iboundary(3) = 4
+          particle%itrdomain(LVL_SUBCELL) = 3
+          particle%iboundary(LVL_SUBCELL) = 4
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         end select
       case (4)
@@ -159,13 +160,13 @@ subroutine pass_mcpq(this, particle)
           infaceoff = -1
         case (2)
           ! N face, subcell 2 --> S face, subcell 1 (cell interior)
-          particle%idomain(3) = 1
-          particle%iboundary(3) = 3
+          particle%itrdomain(LVL_SUBCELL) = 1
+          particle%iboundary(LVL_SUBCELL) = 3
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (3)
           ! N face, subcell 3 --> S face, subcell 4 (cell interior)
-          particle%idomain(3) = 4
-          particle%iboundary(3) = 3
+          particle%itrdomain(LVL_SUBCELL) = 4
+          particle%iboundary(LVL_SUBCELL) = 3
           inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
         case (4)
           ! N face, subcell 4 (cell face)
@@ -181,9 +182,9 @@ subroutine pass_mcpq(this, particle)
       end select
 
       if (inface .eq. -1) then
-        particle%iboundary(2) = 0
+        particle%iboundary(LVL_CELL) = 0
       else if (inface .eq. 0) then
-        particle%iboundary(2) = 0
+        particle%iboundary(LVL_CELL) = 0
       else
         if ((inface .ge. 1) .and. (inface .le. 4)) then
           ! Account for local cell rotation
@@ -192,7 +193,7 @@ subroutine pass_mcpq(this, particle)
           inface = cell%irectvert(inface) + infaceoff
           if (inface .lt. 1) inface = inface + npolyverts
         end if
-        particle%iboundary(2) = inface
+        particle%iboundary(LVL_CELL) = inface
       end if
     end select
   end subroutine pass_mcpq
@@ -234,6 +235,7 @@ end subroutine apply_mcpq
 
   !> @brief Load the rectangular subcell from the rectangular cell
   subroutine load_subcell(this, particle, subcell)
+    use ParticleModule, only: LVL_SUBCELL
     ! dummy
     class(MethodCellPollockQuadType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -251,7 +253,7 @@ subroutine load_subcell(this, particle, subcell)
       factor = factor / cell%defn%porosity
       npolyverts = cell%defn%npolyverts
 
-      isc = particle%idomain(3)
+      isc = particle%itrdomain(LVL_SUBCELL)
       ! Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D,
       ! respectively
 
@@ -277,7 +279,7 @@ subroutine load_subcell(this, particle, subcell)
         end if
 
         subcell%isubcell = isc
-        particle%idomain(3) = isc
+        particle%itrdomain(LVL_SUBCELL) = isc
       end if
       dx = 5d-1 * dx
       dy = 5d-1 * dy
diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90
index 9d653250dae..c523e21ac33 100644
--- a/src/Solution/ParticleTracker/MethodCellTernary.f90
+++ b/src/Solution/ParticleTracker/MethodCellTernary.f90
@@ -99,6 +99,7 @@ end subroutine load_mct
 
   !> @brief Pass particle to next subcell if there is one, or to the cell face
   subroutine pass_mct(this, particle)
+    use ParticleModule, only: LVL_CELL, LVL_SUBCELL
     ! dummy
     class(MethodCellTernaryType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -107,8 +108,8 @@ subroutine pass_mct(this, particle)
     integer(I4B) :: exitFace
     integer(I4B) :: inface
 
-    exitFace = particle%iboundary(3)
-    isc = particle%idomain(3)
+    exitFace = particle%iboundary(LVL_SUBCELL)
+    isc = particle%itrdomain(LVL_SUBCELL)
 
     select case (exitFace)
     case (0)
@@ -122,15 +123,15 @@ subroutine pass_mct(this, particle)
       ! Subcell face --> next subcell in "cycle" (cell interior)
       isc = isc + 1
       if (isc .gt. this%nverts) isc = 1
-      particle%idomain(3) = isc
-      particle%iboundary(3) = 3
+      particle%itrdomain(LVL_SUBCELL) = isc
+      particle%iboundary(LVL_SUBCELL) = 3
       inface = 0
     case (3)
       ! Subcell face --> preceding subcell in "cycle" (cell interior)
       isc = isc - 1
       if (isc .lt. 1) isc = this%nverts
-      particle%idomain(3) = isc
-      particle%iboundary(3) = 2
+      particle%itrdomain(LVL_SUBCELL) = isc
+      particle%iboundary(LVL_SUBCELL) = 2
       inface = 0
     case (4)
       ! Subcell bottom (cell bottom)
@@ -141,11 +142,11 @@ subroutine pass_mct(this, particle)
     end select
 
     if (inface .eq. -1) then
-      particle%iboundary(2) = 0
+      particle%iboundary(LVL_CELL) = 0
     else if (inface .eq. 0) then
-      particle%iboundary(2) = 0
+      particle%iboundary(LVL_CELL) = 0
     else
-      particle%iboundary(2) = inface
+      particle%iboundary(LVL_CELL) = inface
     end if
   end subroutine pass_mct
 
@@ -248,6 +249,7 @@ end subroutine apply_mct
 
   !> @brief Loads a triangular subcell from the polygonal cell
   subroutine load_subcell(this, particle, subcell)
+    use ParticleModule, only: LVL_SUBCELL
     ! dummy
     class(MethodCellTernaryType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -280,7 +282,7 @@ subroutine load_subcell(this, particle, subcell)
     type is (CellPolyType)
       ic = cell%defn%icell
       subcell%icell = ic
-      isc = particle%idomain(3)
+      isc = particle%itrdomain(LVL_SUBCELL)
       if (isc .le. 0) then
         xi = particle%x
         yi = particle%y
@@ -315,7 +317,7 @@ subroutine load_subcell(this, particle, subcell)
 
           call pstop(1)
         else
-          particle%idomain(3) = isc
+          particle%itrdomain(LVL_SUBCELL) = isc
         end if
       end if
       subcell%isubcell = isc
diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90
index 2d3688df8ad..866dac54e05 100644
--- a/src/Solution/ParticleTracker/MethodDis.f90
+++ b/src/Solution/ParticleTracker/MethodDis.f90
@@ -144,7 +144,7 @@ subroutine load_dis(this, particle, next_level, submethod)
 
     select type (cell => this%cell)
     type is (CellRectType)
-      ic = particle%idomain(next_level)
+      ic = particle%itrdomain(next_level)
       call this%load_celldefn(ic, cell%defn)
       call this%load_cell(ic, cell)
       if (this%fmi%ibdgwfsat0(ic) == 0) then
@@ -168,7 +168,7 @@ end subroutine load_dis
   !> @brief Load cell properties into the particle, including
   ! the z coordinate, entry face, and node and layer numbers.
   subroutine load_particle(this, cell, particle)
-    use ParticleModule, only: TERM_BOUNDARY
+    use ParticleModule, only: TERM_BOUNDARY, LVL_CELL
     ! dummy
     class(MethodDisType), intent(inout) :: this
     type(CellRectType), pointer, intent(inout) :: cell
@@ -194,7 +194,7 @@ subroutine load_particle(this, cell, particle)
     select type (dis => this%fmi%dis)
     type is (DisType)
       ! compute reduced/user node numbers and layer
-      inface = particle%iboundary(2)
+      inface = particle%iboundary(LVL_CELL)
       inbr = cell%defn%facenbr(inface)
       idiag = this%fmi%dis%con%ia(cell%defn%icell)
       ipos = idiag + inbr
@@ -210,18 +210,18 @@ subroutine load_particle(this, cell, particle)
       ! in the previous cell.
       if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
         particle%advancing = .false.
-        particle%idomain(2) = particle%icp
+        particle%itrdomain(LVL_CELL) = particle%icp
         particle%istatus = TERM_BOUNDARY
         particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
       else
-        particle%icp = particle%idomain(2)
+        particle%icp = particle%itrdomain(LVL_CELL)
         particle%izp = particle%izone
       end if
 
       ! update node numbers and layer
-      particle%idomain(2) = ic
+      particle%itrdomain(LVL_CELL) = ic
       particle%icu = icu
       particle%ilay = ilay
 
@@ -239,7 +239,7 @@ subroutine load_particle(this, cell, particle)
       else if (inface .eq. 7) then
         inface = 6
       end if
-      particle%iboundary(2) = inface
+      particle%iboundary(LVL_CELL) = inface
 
       ! map z between cells
       z = particle%z
@@ -259,6 +259,7 @@ end subroutine load_particle
   !> @brief Update cell-cell flows of particle mass.
   !! Every particle is currently assigned unit mass.
   subroutine update_flowja(this, cell, particle)
+    use ParticleModule, only: LVL_CELL
     ! dummy
     class(MethodDisType), intent(inout) :: this
     type(CellRectType), pointer, intent(inout) :: cell
@@ -269,7 +270,7 @@ subroutine update_flowja(this, cell, particle)
     integer(I4B) :: inface
     integer(I4B) :: ipos
 
-    inface = particle%iboundary(2)
+    inface = particle%iboundary(LVL_CELL)
     inbr = cell%defn%facenbr(inface)
     idiag = this%fmi%dis%con%ia(cell%defn%icell)
     ipos = idiag + inbr
@@ -284,7 +285,7 @@ end subroutine update_flowja
 
   !> @brief Pass a particle to the next cell, if there is one
   subroutine pass_dis(this, particle)
-    use ParticleModule, only: TERM_BOUNDARY
+    use ParticleModule, only: TERM_BOUNDARY, LVL_CELL
     ! dummy
     class(MethodDisType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -297,7 +298,7 @@ subroutine pass_dis(this, particle)
       ! If the entry face has no neighbors it's a
       ! boundary face, so terminate the particle.
       ! todo AMP: reconsider when multiple models supported
-      if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
+      if (cell%defn%facenbr(particle%iboundary(LVL_CELL)) .eq. 0) then
         particle%istatus = TERM_BOUNDARY
         particle%advancing = .false.
         call this%save(particle, reason=3)
diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90
index 8f178ed2272..f77f4512ea1 100644
--- a/src/Solution/ParticleTracker/MethodDisv.f90
+++ b/src/Solution/ParticleTracker/MethodDisv.f90
@@ -85,7 +85,7 @@ subroutine load_disv(this, particle, next_level, submethod)
 
     select type (cell => this%cell)
     type is (CellPolyType)
-      ic = particle%idomain(next_level)
+      ic = particle%itrdomain(next_level)
       call this%load_cell_defn(ic, cell%defn)
       if (this%fmi%ibdgwfsat0(ic) == 0) then
         ! Cell is active but dry, so select and initialize pass-to-bottom
@@ -141,7 +141,7 @@ end subroutine load_disv
   subroutine load_particle(this, cell, particle)
     ! modules
     use DisvModule, only: DisvType
-    use ParticleModule, only: TERM_BOUNDARY
+    use ParticleModule, only: TERM_BOUNDARY, LVL_CELL, LVL_SUBCELL
     ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(CellPolyType), pointer, intent(inout) :: cell
@@ -159,7 +159,7 @@ subroutine load_particle(this, cell, particle)
 
     select type (dis => this%fmi%dis)
     type is (DisvType)
-      inface = particle%iboundary(2)
+      inface = particle%iboundary(LVL_CELL)
       idiag = dis%con%ia(cell%defn%icell)
       inbr = cell%defn%facenbr(inface)
       ipos = idiag + inbr
@@ -174,32 +174,33 @@ subroutine load_particle(this, cell, particle)
       ! in the previous cell.
       if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
         particle%advancing = .false.
-        particle%idomain(2) = particle%icp
+        particle%itrdomain(LVL_CELL) = particle%icp
         particle%istatus = TERM_BOUNDARY
         particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
       else
-        particle%icp = particle%idomain(2)
+        particle%icp = particle%itrdomain(LVL_CELL)
         particle%izp = particle%izone
       end if
 
-      particle%idomain(2) = ic
+      particle%itrdomain(LVL_CELL) = ic
       particle%icu = icu
       particle%ilay = ilay
 
       z = particle%z
       call this%map_neighbor(cell%defn, inface, z)
 
-      particle%iboundary(2) = inface
-      particle%idomain(3:) = 0
-      particle%iboundary(3:) = 0
+      particle%iboundary(LVL_CELL) = inface
+      particle%itrdomain(LVL_SUBCELL:) = 0
+      particle%iboundary(LVL_SUBCELL:) = 0
       particle%z = z
     end select
 
   end subroutine
 
   subroutine update_flowja(this, cell, particle)
+    use ParticleModule, only: LVL_CELL
     ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(CellPolyType), pointer, intent(inout) :: cell
@@ -210,7 +211,7 @@ subroutine update_flowja(this, cell, particle)
     integer(I4B) :: ipos
 
     idiag = this%fmi%dis%con%ia(cell%defn%icell)
-    inbr = cell%defn%facenbr(particle%iboundary(2))
+    inbr = cell%defn%facenbr(particle%iboundary(LVL_CELL))
     ipos = idiag + inbr
 
     ! leaving old cell
@@ -224,7 +225,7 @@ end subroutine update_flowja
 
   !> @brief Pass a particle to the next cell, if there is one
   subroutine pass_disv(this, particle)
-    use ParticleModule, only: TERM_BOUNDARY
+    use ParticleModule, only: TERM_BOUNDARY, LVL_CELL
     ! dummy
     class(MethodDisvType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -237,7 +238,7 @@ subroutine pass_disv(this, particle)
       ! If the entry face has no neighbors it's a
       ! boundary face, so terminate the particle.
       ! todo AMP: reconsider when multiple models supported
-      if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
+      if (cell%defn%facenbr(particle%iboundary(LVL_CELL)) .eq. 0) then
         particle%istatus = TERM_BOUNDARY
         particle%advancing = .false.
         call this%save(particle, reason=3)
diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90
index c002a567b5a..75c3178ac58 100644
--- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90
+++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90
@@ -84,7 +84,7 @@ end subroutine apply_msp
   !! this context and for any modifications or errors.
   !<
   subroutine track_subcell(this, subcell, particle, tmax)
-    use ParticleModule, only: ACTIVE, TERM_NO_EXITS_SUB
+    use ParticleModule, only: ACTIVE, TERM_NO_EXITS_SUB, LVL_SUBCELL
     ! dummy
     class(MethodSubcellPollockType), intent(inout) :: this
     class(SubcellRectType), intent(in) :: subcell
@@ -260,7 +260,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     particle%y = y * subcell%dy
     particle%z = z * subcell%dz
     particle%ttrack = t
-    particle%iboundary(3) = exitFace
+    particle%iboundary(LVL_SUBCELL) = exitFace
 
     ! Save particle track record
     if (reason >= 0) call this%save(particle, reason=reason)
diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90
index 98014080b0f..1b882338ec4 100644
--- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90
+++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90
@@ -63,7 +63,7 @@ end subroutine apply_mst
 
   !> @brief Track a particle across a triangular subcell.
   subroutine track_subcell(this, subcell, particle, tmax)
-    use ParticleModule, only: ACTIVE, TERM_NO_EXITS_SUB
+    use ParticleModule, only: ACTIVE, TERM_NO_EXITS_SUB, LVL_SUBCELL
     ! dummy
     class(MethodSubcellTernaryType), intent(inout) :: this
     class(SubcellTriType), intent(in) :: subcell
@@ -202,7 +202,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     vziodz = vzi / dz
 
     ! If possible, track the particle across the subcell.
-    itrifaceenter = particle%iboundary(3) - 1
+    itrifaceenter = particle%iboundary(LVL_SUBCELL) - 1
     if (itrifaceenter == -1) itrifaceenter = 999
     call traverse_triangle(isolv, tol, &
                            dtexitxy, alpexit, betexit, &
@@ -292,7 +292,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     particle%y = y
     particle%z = z
     particle%ttrack = t
-    particle%iboundary(3) = exitFace
+    particle%iboundary(LVL_SUBCELL) = exitFace
 
     call this%save(particle, reason=reason)
   end subroutine track_subcell
diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90
index db671f4d639..0d39228179b 100644
--- a/src/Solution/ParticleTracker/Particle.f90
+++ b/src/Solution/ParticleTracker/Particle.f90
@@ -14,6 +14,18 @@ module ParticleModule
   !! more methods as appropriate for finer-grained levels.
   integer, parameter :: MAX_LEVEL = 4
 
+  !> @brief Particle tracking level.
+  !!
+  !! Identifies the domain over which a tracking method is responsible for
+  !! moving a particle. Tracking methods at lower levels (coarser-grained)
+  !! defer to higher-level (finer-grained) methods within their subdomains.
+  !<
+  enum, bind(C)
+    enumerator :: LVL_MODEL = 1
+    enumerator :: LVL_CELL = 2
+    enumerator :: LVL_SUBCELL = 3
+  end enum
+
   !> @brief Particle status enumeration.
   !!
   !! Particles begin in status 1 (active) at release time. Status may only
@@ -88,7 +100,7 @@ module ParticleModule
     integer(I4B), public :: istopzone !< stop zone number
     integer(I4B), public :: idrymeth !< dry tracking method
     ! state
-    integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy ! TODO: rename to itdomain? idomain
+    integer(I4B), allocatable, public :: itrdomain(:) !< tracking domain
     integer(I4B), allocatable, public :: iboundary(:) !< tracking domain boundaries
     integer(I4B), public :: icp !< previous cell number (reduced)
     integer(I4B), public :: icu !< user cell number
@@ -132,7 +144,7 @@ module ParticleModule
     integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number
     integer(I4B), dimension(:), pointer, public, contiguous :: idrymeth !< stop in dry cells
     ! state
-    integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy
+    integer(I4B), dimension(:, :), pointer, public, contiguous :: itrdomain !< array of indices for domains in the tracking domain hierarchy
     integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries
     integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user)
     integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer
@@ -163,7 +175,7 @@ module ParticleModule
   subroutine create_particle(particle)
     type(ParticleType), pointer :: particle !< particle
     allocate (particle)
-    allocate (particle%idomain(MAX_LEVEL))
+    allocate (particle%itrdomain(MAX_LEVEL))
     allocate (particle%iboundary(MAX_LEVEL))
   end subroutine create_particle
 
@@ -196,7 +208,7 @@ subroutine create_particle_store(store, np, mempath)
     call mem_allocate(store%iexmeth, np, 'PLIEXMETH', mempath)
     call mem_allocate(store%extol, np, 'PLEXTOL', mempath)
     call mem_allocate(store%extend, np, 'PLIEXTEND', mempath)
-    call mem_allocate(store%idomain, np, MAX_LEVEL, 'PLIDOMAIN', mempath)
+    call mem_allocate(store%itrdomain, np, MAX_LEVEL, 'PLITRDOMAIN', mempath)
     call mem_allocate(store%iboundary, np, MAX_LEVEL, 'PLIBOUNDARY', mempath)
   end subroutine create_particle_store
 
@@ -227,7 +239,7 @@ subroutine destroy(this, mempath)
     call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath)
     call mem_deallocate(this%extol, 'PLEXTOL', mempath)
     call mem_deallocate(this%extend, 'PLIEXTEND', mempath)
-    call mem_deallocate(this%idomain, 'PLIDOMAIN', mempath)
+    call mem_deallocate(this%itrdomain, 'PLITRDOMAIN', mempath)
     call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath)
   end subroutine destroy
 
@@ -261,7 +273,7 @@ subroutine resize(this, np, mempath)
     call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath)
     call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
     call mem_reallocate(this%extend, np, 'PLIEXTEND', mempath)
-    call mem_reallocate(this%idomain, np, MAX_LEVEL, 'PLIDOMAIN', mempath)
+    call mem_reallocate(this%itrdomain, np, MAX_LEVEL, 'PLITRDOMAIN', mempath)
     call mem_reallocate(this%iboundary, np, MAX_LEVEL, 'PLIBOUNDARY', mempath)
   end subroutine resize
 
@@ -299,9 +311,9 @@ subroutine get(this, particle, imdl, iprp, ip)
     particle%tstop = this%tstop(ip)
     particle%ttrack = this%ttrack(ip)
     particle%advancing = .true.
-    particle%idomain(1:MAX_LEVEL) = &
-      this%idomain(ip, 1:MAX_LEVEL)
-    particle%idomain(1) = imdl
+    particle%itrdomain(1:MAX_LEVEL) = &
+      this%itrdomain(ip, 1:MAX_LEVEL)
+    particle%itrdomain(1) = imdl
     particle%iboundary(1:MAX_LEVEL) = &
       this%iboundary(ip, 1:MAX_LEVEL)
     particle%ifrctrn = this%ifrctrn(ip)
@@ -334,10 +346,10 @@ subroutine put(this, particle, ip)
     this%trelease(ip) = particle%trelease
     this%tstop(ip) = particle%tstop
     this%ttrack(ip) = particle%ttrack
-    this%idomain( &
+    this%itrdomain( &
       ip, &
       1:MAX_LEVEL) = &
-      particle%idomain(1:MAX_LEVEL)
+      particle%itrdomain(1:MAX_LEVEL)
     this%iboundary( &
       ip, &
       1:MAX_LEVEL) = &