diff --git a/autotest/test_gwf_libmf6_rch01.py b/autotest/test_gwf_libmf6_rch01.py
index bb8252e8cff..76b54f3cdce 100644
--- a/autotest/test_gwf_libmf6_rch01.py
+++ b/autotest/test_gwf_libmf6_rch01.py
@@ -177,7 +177,7 @@ def api_func(exe, idx, model_ws=None):
     max_iter = mf6.get_value(mxit_tag)
 
     # get copy of recharge array
-    rch_tag = mf6.get_var_address("BOUND", name, rch_pname)
+    rch_tag = mf6.get_var_address("RECHARGE", name, rch_pname)
     new_recharge = mf6.get_value(rch_tag)
 
     # model time loop
@@ -193,7 +193,7 @@ def api_func(exe, idx, model_ws=None):
         mf6.prepare_solve()
 
         # update recharge
-        new_recharge[:, 0] = rch_spd[idx] * area
+        new_recharge[:] = rch_spd[idx]
         mf6.set_value(rch_tag, new_recharge)
 
         while kiter < max_iter:
diff --git a/autotest/test_gwf_libmf6_rch02.py b/autotest/test_gwf_libmf6_rch02.py
index 1e7a5797a3e..6c644d5f8f3 100644
--- a/autotest/test_gwf_libmf6_rch02.py
+++ b/autotest/test_gwf_libmf6_rch02.py
@@ -172,7 +172,7 @@ def run_perturbation(mf6, max_iter, recharge, tag, rch):
     kiter = 0
     while kiter < max_iter:
         # update recharge
-        recharge[:, 0] = rch * area
+        recharge[:] = rch
         mf6.set_value(tag, recharge)
         # solve with updated well rate
         has_converged = mf6.solve()
@@ -222,7 +222,7 @@ def api_func(exe, idx, model_ws=None):
     max_iter = mf6.get_value(mxit_tag)
 
     # get copy of recharge array
-    rch_tag = mf6.get_var_address("BOUND", name, rch_pname)
+    rch_tag = mf6.get_var_address("RECHARGE", name, rch_pname)
     new_recharge = mf6.get_value(rch_tag).copy()
 
     # determine initial recharge value
diff --git a/doc/mf6io/mf6ivar/dfn/gwf-rch.dfn b/doc/mf6io/mf6ivar/dfn/gwf-rch.dfn
index 0a36ae078be..dbb47422bc0 100644
--- a/doc/mf6io/mf6ivar/dfn/gwf-rch.dfn
+++ b/doc/mf6io/mf6ivar/dfn/gwf-rch.dfn
@@ -1,5 +1,6 @@
 # --------------------- gwf rch options ---------------------
 # flopy multi-package
+# modflow6 aux-sfac-param recharge
 
 block options
 name fixed_cell
@@ -44,6 +45,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'recharge'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -52,6 +54,7 @@ reader urword
 optional true
 longname print recharge rates to listing file
 description REPLACE print_flows {'{#1}': 'recharge'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -60,6 +63,7 @@ reader urword
 optional true
 longname save recharge to budget file
 description REPLACE save_flows {'{#1}': 'recharge'}
+mf6internal ipakcb
 
 block options
 name ts_filerecord
@@ -169,6 +173,7 @@ shape (maxbound)
 reader urword
 longname
 description
+mf6internal spd
 
 block period
 name cellid
@@ -202,6 +207,7 @@ optional true
 time_series true
 longname auxiliary variables
 description REPLACE aux {'{#1}': 'recharge'}
+mf6internal auxvar
 
 block period
 name boundname
diff --git a/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn b/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn
index bc5874ef599..e8598e47da1 100644
--- a/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn
+++ b/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn
@@ -1,5 +1,6 @@
 # --------------------- gwf rcha options ---------------------
 # flopy multi-package
+# modflow6 aux-sfac-param recharge
 
 block options
 name readasarrays
@@ -45,6 +46,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'recharge'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -53,6 +55,7 @@ reader urword
 optional true
 longname print recharge rates to listing file
 description REPLACE print_flows {'{#1}': 'recharge'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -61,6 +64,7 @@ reader urword
 optional true
 longname save CHD flows to budget file
 description REPLACE save_flows {'{#1}': 'recharge'}
+mf6internal ipakcb
 
 block options
 name tas_filerecord
@@ -168,6 +172,7 @@ name recharge
 type double precision
 shape (ncol*nrow; ncpl)
 reader readarray
+time_series true
 longname recharge rate
 description is the recharge flux rate ($LT^{-1}$).  This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the "Using Time-Array Series in a Package" section).
 default_value 1.e-3
@@ -177,6 +182,8 @@ name aux
 type double precision
 shape (ncol*nrow; ncpl)
 reader readarray
+time_series true
 optional true
 longname auxiliary variable iaux
 description is an array of values for auxiliary variable aux(iaux), where iaux is a value from 1 to naux, and aux(iaux) must be listed as part of the auxiliary variables.  A separate array can be specified for each auxiliary variable.  If an array is not specified for an auxiliary variable, then a value of zero is assigned.  If the value specified here for the auxiliary variable is the same as auxmultname, then the recharge array will be multiplied by this array.
+mf6internal auxvar
diff --git a/make/makefile b/make/makefile
index 97658431f17..17e860d7f56 100644
--- a/make/makefile
+++ b/make/makefile
@@ -104,6 +104,8 @@ $(OBJDIR)/gwt1disu1idm.o \
 $(OBJDIR)/gwt1dis1idm.o \
 $(OBJDIR)/gwf3wel8idm.o \
 $(OBJDIR)/gwf3riv8idm.o \
+$(OBJDIR)/gwf3rch8idm.o \
+$(OBJDIR)/gwf3rcha8idm.o \
 $(OBJDIR)/gwf3npf8idm.o \
 $(OBJDIR)/gwf3idm.o \
 $(OBJDIR)/gwf3ghb8idm.o \
diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
index 09cc9fa9b92..c7d0fdaa754 100644
--- a/msvs/mf6core.vfproj
+++ b/msvs/mf6core.vfproj
@@ -146,6 +146,8 @@
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3obs8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3oc8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8idm.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rcha8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sfr8.f90"/>
diff --git a/src/Model/GroundWaterFlow/gwf3.f90 b/src/Model/GroundWaterFlow/gwf3.f90
index 90df71a9f31..8000150e376 100644
--- a/src/Model/GroundWaterFlow/gwf3.f90
+++ b/src/Model/GroundWaterFlow/gwf3.f90
@@ -1305,7 +1305,8 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
       call ghb_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
                       pakname, mempath)
     case ('RCH6')
-      call rch_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
+      call rch_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
+                      pakname, mempath)
     case ('EVT6')
       call evt_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
     case ('MAW6')
diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90
index 1380a9503b5..c734751439f 100644
--- a/src/Model/GroundWaterFlow/gwf3dis8.f90
+++ b/src/Model/GroundWaterFlow/gwf3dis8.f90
@@ -36,6 +36,7 @@ module GwfDisModule
     procedure, public :: read_layer_array
     procedure, public :: record_srcdst_list_header
     procedure, public :: nlarray_to_nodelist
+    procedure, public :: nlarray_to_nodelist2
     ! -- helper functions
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx3
@@ -1747,4 +1748,91 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, &
     ! -- return
   end subroutine nlarray_to_nodelist
 
+  subroutine nlarray_to_nodelist2(this, darray, nodelist, maxbnd, nbound, aname)
+! ******************************************************************************
+! nlarray_to_nodelist -- Read an integer array into nodelist. For structured
+!                        model, integer array is layer number; for unstructured
+!                        model, integer array is node number.
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use InputOutputModule, only: get_node
+    use ConstantsModule, only: LINELENGTH
+    ! -- dummy
+    class(GwfDisType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    ! -- local
+    integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
+! ------------------------------------------------------------------------------
+    !
+    ! -- set variables
+    nlay = this%mshape(1)
+    nrow = this%mshape(2)
+    ncol = this%mshape(3)
+    !
+    if (this%ndim > 1) then
+      !
+      nval = ncol * nrow
+      !
+      ! -- Copy array into nodelist
+      ipos = 1
+      ierr = 0
+      do ir = 1, nrow
+        do ic = 1, ncol
+          nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
+          il = darray(nodeu)
+          if (il < 1 .or. il > nlay) then
+            write (errmsg, '(a,1x,i0)') 'Invalid layer number:', il
+            call store_error(errmsg, terminate=.TRUE.)
+          end if
+          nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
+          noder = this%get_nodenumber(nodeu, 0)
+          if (ipos > maxbnd) then
+            ierr = ipos
+          else
+            nodelist(ipos) = noder
+          end if
+          ipos = ipos + 1
+        end do
+      end do
+      !
+      ! -- Check for errors
+      nbound = ipos - 1
+      if (ierr > 0) then
+        write (errmsg, '(a,1x,i0)') &
+          'MAXBOUND dimension is too small.'// &
+          'INCREASE MAXBOUND TO:', ierr
+        call store_error(errmsg, terminate=.TRUE.)
+      end if
+      !
+      ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
+      if (nbound < maxbnd) then
+        do ipos = nbound + 1, maxbnd
+          nodelist(ipos) = 0
+        end do
+      end if
+      !
+    else
+      !
+      ! -- For unstructured, read nodelist directly, then check node numbers
+      nodelist = darray
+      do noder = 1, maxbnd
+        if (noder < 1 .or. noder > this%nodes) then
+          write (errmsg, '(a,1x,i0)') 'Invalid node number:', noder
+          call store_error(errmsg, terminate=.TRUE.)
+        end if
+      end do
+      nbound = maxbnd
+      !
+    end if
+    !
+    ! -- return
+  end subroutine nlarray_to_nodelist2
+
 end module GwfDisModule
diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90
index 06d07c8d78f..a62bb66b3de 100644
--- a/src/Model/GroundWaterFlow/gwf3disv8.f90
+++ b/src/Model/GroundWaterFlow/gwf3disv8.f90
@@ -38,6 +38,7 @@ module GwfDisvModule
     procedure, public :: read_layer_array
     procedure, public :: record_srcdst_list_header
     procedure, public :: nlarray_to_nodelist
+    procedure, public :: nlarray_to_nodelist2
     ! -- helper functions
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx2
@@ -1981,4 +1982,74 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, &
     ! -- return
   end subroutine nlarray_to_nodelist
 
+  subroutine nlarray_to_nodelist2(this, darray, nodelist, maxbnd, nbound, aname)
+! ******************************************************************************
+! nlarray_to_nodelist -- Read an integer array into nodelist. For structured
+!                        model, integer array is layer number; for unstructured
+!                        model, integer array is node number.
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use InputOutputModule, only: get_node
+    ! -- dummy
+    class(GwfDisvType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    ! -- local
+    integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
+! ------------------------------------------------------------------------------
+    !
+    ! -- set variables
+    nlay = this%mshape(1)
+    nrow = 1
+    ncol = this%mshape(2)
+    !
+    nval = ncol * nrow
+    !
+    ! -- Copy array into nodelist
+    ipos = 1
+    ierr = 0
+    do ir = 1, nrow
+      do ic = 1, ncol
+        nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
+        il = darray(nodeu)
+        if (il < 1 .or. il > nlay) then
+          write (errmsg, '(a,i0,a)') &
+            'Invalid layer number (', il, ').'
+          call store_error(errmsg, terminate=.TRUE.)
+        end if
+        nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
+        noder = this%get_nodenumber(nodeu, 0)
+        if (ipos > maxbnd) then
+          ierr = ipos
+        else
+          nodelist(ipos) = noder
+        end if
+        ipos = ipos + 1
+      end do
+    end do
+    !
+    ! -- Check for errors
+    nbound = ipos - 1
+    if (ierr > 0) then
+      write (errmsg, '(a,i0,a)') &
+        'MAXBOUND dimension is too small. Increase MAXBOUND to ', ierr, '.'
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
+    !
+    ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
+    if (nbound < maxbnd) then
+      do ipos = nbound + 1, maxbnd
+        nodelist(ipos) = 0
+      end do
+    end if
+    !
+    ! -- return
+  end subroutine nlarray_to_nodelist2
+
 end module GwfDisvModule
diff --git a/src/Model/GroundWaterFlow/gwf3rch8.f90 b/src/Model/GroundWaterFlow/gwf3rch8.f90
index d0add5b9f79..539ba60b092 100644
--- a/src/Model/GroundWaterFlow/gwf3rch8.f90
+++ b/src/Model/GroundWaterFlow/gwf3rch8.f90
@@ -1,17 +1,17 @@
 module RchModule
   !
-  use KindModule, only: DP, I4B
+  use KindModule, only: DP, I4B, LGP
   use ConstantsModule, only: DZERO, LENFTYPE, LENPACKAGENAME, MAXCHARLEN, &
-                             IWETLAKE
+                             IWETLAKE, LINELENGTH
   use MemoryHelperModule, only: create_mem_path
   use BndModule, only: BndType
-  use SimModule, only: store_error, store_error_unit
+  use BndExtModule, only: BndExtType
+  use SimModule, only: store_error, store_error_filename
   use SimVariablesModule, only: errmsg
   use ObsModule, only: DefaultObsIdProcessor
   use TimeArraySeriesLinkModule, only: TimeArraySeriesLinkType
-  use TimeSeriesLinkModule, only: TimeSeriesLinkType, &
-                                  GetTimeSeriesLinkFromList
   use BlockParserModule, only: BlockParserType
+  use CharacterStringModule, only: CharacterStringType
   use MatrixBaseModule
   !
   implicit none
@@ -23,25 +23,25 @@ module RchModule
   character(len=LENPACKAGENAME) :: text = '             RCH'
   character(len=LENPACKAGENAME) :: texta = '            RCHA'
   !
-  type, extends(BndType) :: RchType
-    integer(I4B), pointer :: inirch => NULL()
+  type, extends(BndExtType) :: RchType
+    real(DP), dimension(:), pointer, contiguous :: recharge => null() !< boundary recharge array
     integer(I4B), dimension(:), pointer, contiguous :: nodesontop => NULL() ! User provided cell numbers; nodelist is cells where recharge is applied)
-    logical, private :: fixed_cell = .false.
-    logical, private :: read_as_arrays = .false.
+    logical, pointer, private :: fixed_cell
+    logical, pointer, private :: read_as_arrays
   contains
     procedure :: rch_allocate_scalars
-    procedure :: bnd_options => rch_options
-    procedure :: read_dimensions => rch_read_dimensions
+    procedure :: allocate_arrays => rch_allocate_arrays
+    procedure :: source_options => rch_source_options
+    procedure :: source_dimensions => rch_source_dimensions
+    procedure :: log_rch_options
     procedure :: read_initial_attr => rch_read_initial_attr
     procedure :: bnd_rp => rch_rp
-    procedure :: set_nodesontop
     procedure :: bnd_cf => rch_cf
     procedure :: bnd_fc => rch_fc
     procedure :: bnd_da => rch_da
+    procedure :: set_nodesontop
     procedure :: define_listlabel => rch_define_listlabel
-    procedure, public :: bnd_rp_ts => rch_rp_ts
-    procedure, private :: rch_rp_array
-    procedure, private :: rch_rp_list
+    procedure :: bound_value => rch_bound_value
     procedure, private :: default_nodelist
     ! -- for observations
     procedure, public :: bnd_obs_supported => rch_obs_supported
@@ -50,7 +50,8 @@ module RchModule
 
 contains
 
-  subroutine rch_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
+  subroutine rch_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
+                        mempath)
 ! ******************************************************************************
 ! rch_create -- Create a New Recharge Package
 ! Subroutine: (1) create new-style package
@@ -67,6 +68,7 @@ subroutine rch_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
     integer(I4B), intent(in) :: iout
     character(len=*), intent(in) :: namemodel
     character(len=*), intent(in) :: pakname
+    character(len=*), intent(in) :: mempath
     ! -- local
     type(rchtype), pointer :: rchobj
 ! ------------------------------------------------------------------------------
@@ -76,7 +78,7 @@ subroutine rch_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
     packobj => rchobj
     !
     ! -- create name and memory path
-    call packobj%set_names(ibcnum, namemodel, pakname, ftype)
+    call packobj%set_names(ibcnum, namemodel, pakname, ftype, mempath)
     packobj%text = text
     !
     ! -- allocate scalars
@@ -92,10 +94,6 @@ subroutine rch_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
     packobj%ncolbnd = 1
     packobj%iscloc = 1 ! sfac applies to recharge rate
     packobj%ictMemPath = create_mem_path(namemodel, 'NPF')
-    ! indxconvertflux is Column index of bound that will be multiplied by
-    ! cell area to convert flux rates to flow rates
-    packobj%indxconvertflux = 1
-    packobj%AllowTimeArraySeries = .true.
     !
     ! -- return
     return
@@ -103,104 +101,154 @@ end subroutine rch_create
 
   subroutine rch_allocate_scalars(this)
 ! ******************************************************************************
-! allocate_scalars -- allocate scalar members
+! rch_allocate_scalars -- allocate scalar members
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
 ! ------------------------------------------------------------------------------
     ! -- modules
-    use MemoryManagerModule, only: mem_allocate
     ! -- dummy
     class(RchType), intent(inout) :: this
 ! ------------------------------------------------------------------------------
     !
-    ! -- call standard BndType allocate scalars
-    call this%BndType%allocate_scalars()
+    ! -- allocate base scalars
+    call this%BndExtType%allocate_scalars()
     !
-    ! -- allocate the object and assign values to object variables
-    call mem_allocate(this%inirch, 'INIRCH', this%memoryPath)
+    ! -- allocate internal members
+    allocate (this%fixed_cell)
+    allocate (this%read_as_arrays)
     !
     ! -- Set values
-    this%inirch = 0
     this%fixed_cell = .false.
+    this%read_as_arrays = .false.
     !
     ! -- return
     return
   end subroutine rch_allocate_scalars
 
-  subroutine rch_options(this, option, found)
+  subroutine rch_allocate_arrays(this, nodelist, auxvar)
+! ******************************************************************************
+! rch_allocate_arrays -- allocate arrays
 ! ******************************************************************************
-! rch_options -- set options specific to RchType
 !
-! rch_options overrides BndType%bnd_options
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use MemoryManagerModule, only: mem_setptr, mem_checkin
+    ! -- dummy
+    class(RchType) :: this
+    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
+    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
+    ! -- local
+! ------------------------------------------------------------------------------
+    !
+    ! -- allocate base arrays
+    call this%BndExtType%allocate_arrays(nodelist, auxvar)
+    !
+    ! -- set recharge input context pointer
+    call mem_setptr(this%recharge, 'RECHARGE', this%input_mempath)
+    !
+    ! -- checkin recharge input context pointer
+    call mem_checkin(this%recharge, 'RECHARGE', this%memoryPath, &
+                     'RECHARGE', this%input_mempath)
+    !
+    ! -- return
+    return
+  end subroutine rch_allocate_arrays
+
+  subroutine rch_source_options(this)
+! ******************************************************************************
+! rch_source_options -- source options specific to RchType
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
 ! ------------------------------------------------------------------------------
-    use ConstantsModule, only: DZERO
-    use SimModule, only: store_error
+    use MemoryManagerExtModule, only: mem_set_value
+    use IdmGwfDfnSelectorModule, only: GwfParamFoundType
+    implicit none
+    ! -- dummy
+    class(RchType), intent(inout) :: this
+    ! -- local
+    type(GwfParamFoundType) :: found
+! ------------------------------------------------------------------------------
+    !
+    ! -- source common bound options
+    call this%BndExtType%source_options()
+    !
+    ! -- update defaults with idm sourced values
+    call mem_set_value(this%fixed_cell, 'FIXED_CELL', this%input_mempath, &
+                       found%fixed_cell)
+    call mem_set_value(this%read_as_arrays, 'READASARRAYS', this%input_mempath, &
+                       found%readasarrays)
+    !
+    if (found%readasarrays) then
+      if (this%dis%supports_layers()) then
+        this%text = texta
+      else
+        errmsg = 'READASARRAYS option is not compatible with selected'// &
+                 ' discretization type.'
+        call store_error(errmsg)
+        call store_error_filename(this%input_fname)
+      end if
+    end if
+    !
+    ! -- log rch params
+    call this%log_rch_options(found)
+    !
+    ! -- return
+    return
+  end subroutine rch_source_options
+
+  subroutine log_rch_options(this, found)
+! ******************************************************************************
+! log_rch_options -- log options specific to RchType
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use IdmGwfDfnSelectorModule, only: GwfParamFoundType
     implicit none
     ! -- dummy
     class(RchType), intent(inout) :: this
-    character(len=*), intent(inout) :: option
-    logical, intent(inout) :: found
+    type(GwfParamFoundType), intent(in) :: found
     ! -- local
-    character(len=MAXCHARLEN) :: ermsg
     ! -- formats
-    character(len=*), parameter :: fmtihact = &
-      &"(4x, 'RECHARGE WILL BE APPLIED TO HIGHEST ACTIVE CELL.')"
     character(len=*), parameter :: fmtfixedcell = &
       &"(4x, 'RECHARGE WILL BE APPLIED TO SPECIFIED CELL.')"
     character(len=*), parameter :: fmtreadasarrays = &
       &"(4x, 'RECHARGE INPUT WILL BE READ AS ARRAY(S).')"
 ! ------------------------------------------------------------------------------
     !
-    ! -- Check for FIXED_CELL and READASARRAYS
-    select case (option)
-    case ('FIXED_CELL')
-      this%fixed_cell = .true.
+    ! -- log found options
+    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
+      //' OPTIONS'
+    !
+    if (found%fixed_cell) then
       write (this%iout, fmtfixedcell)
-      found = .true.
-    case ('READASARRAYS')
-      if (this%dis%supports_layers()) then
-        this%read_as_arrays = .true.
-        this%text = texta
-      else
-        ermsg = 'READASARRAYS option is not compatible with selected'// &
-                ' discretization type.'
-        call store_error(ermsg)
-        call this%parser%StoreErrorUnit()
-      end if
-      !
-      ! -- Write option
+    end if
+    !
+    if (found%readasarrays) then
       write (this%iout, fmtreadasarrays)
-      !
-      found = .true.
-    case default
-      !
-      ! -- No options found
-      found = .false.
-    end select
+    end if
+    !
+    ! -- close logging block
+    write (this%iout, '(1x,a)') &
+      'END OF '//trim(adjustl(this%text))//' OPTIONS'
     !
     ! -- return
     return
-  end subroutine rch_options
+  end subroutine log_rch_options
 
-  subroutine rch_read_dimensions(this)
+  subroutine rch_source_dimensions(this)
 ! ******************************************************************************
-! bnd_read_dimensions -- Read the dimensions for this package
+! rch_source_dimensions -- Source the dimensions for this package
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
 ! ------------------------------------------------------------------------------
-    use ConstantsModule, only: LINELENGTH
-    use SimModule, only: store_error, store_error_unit
     ! -- dummy
     class(RchType), intent(inout) :: this
     ! -- local
-    character(len=LINELENGTH) :: keyword
-    integer(I4B) :: ierr
-    logical :: isfound, endOfBlock
     ! -- format
 ! ------------------------------------------------------------------------------
     !
@@ -209,45 +257,19 @@ subroutine rch_read_dimensions(this)
     !   (2) READASARRAYS option has been specified.
     if (this%read_as_arrays) then
       this%maxbound = this%dis%get_ncpl()
-    else
-      ! -- get dimensions block
-      call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
-                                supportOpenClose=.true.)
       !
-      ! -- parse dimensions block if detected
-      if (isfound) then
-        write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text))// &
-          ' DIMENSIONS'
-        do
-          call this%parser%GetNextLine(endOfBlock)
-          if (endOfBlock) exit
-          call this%parser%GetStringCaps(keyword)
-          select case (keyword)
-          case ('MAXBOUND')
-            this%maxbound = this%parser%GetInteger()
-            write (this%iout, '(4x,a,i7)') 'MAXBOUND = ', this%maxbound
-          case default
-            write (errmsg, '(a,a)') &
-              'Unknown '//trim(this%text)//' DIMENSION: ', trim(keyword)
-            call store_error(errmsg)
-            call this%parser%StoreErrorUnit()
-          end select
-        end do
-        !
-        write (this%iout, '(1x,a)') &
-          'END OF '//trim(adjustl(this%text))//' DIMENSIONS'
-      else
-        call store_error('Required DIMENSIONS block not found.')
-        call this%parser%StoreErrorUnit()
+      ! -- verify dimensions were set
+      if (this%maxbound <= 0) then
+        write (errmsg, '(a)') &
+          'MAXBOUND must be an integer greater than zero.'
+        call store_error(errmsg)
+        call store_error_filename(this%input_fname)
       end if
-    end if
-    !
-    ! -- verify dimensions were set
-    if (this%maxbound <= 0) then
-      write (errmsg, '(a)') &
-        'MAXBOUND must be an integer greater than zero.'
-      call store_error(errmsg)
-      call this%parser%StoreErrorUnit()
+      !
+    else
+      !
+      ! -- source maxbound
+      call this%BndExtType%source_dimensions()
     end if
     !
     ! -- Call define_listlabel to construct the list label that is written
@@ -256,12 +278,11 @@ subroutine rch_read_dimensions(this)
     !
     ! -- return
     return
-  end subroutine rch_read_dimensions
+  end subroutine rch_source_dimensions
 
   subroutine rch_read_initial_attr(this)
 ! ******************************************************************************
 ! rch_read_initial_attr -- Part of allocate and read
-! If READASARRAYS has been specified, assign default IRCH = 1
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
@@ -286,334 +307,38 @@ subroutine rch_rp(this)
 !
 !    SPECIFICATIONS:
 ! ------------------------------------------------------------------------------
-    use ConstantsModule, only: LINELENGTH
-    use TdisModule, only: kper, nper
-    use SimModule, only: store_error
+    use TdisModule, only: kper
     implicit none
     ! -- dummy
     class(RchType), intent(inout) :: this
     ! -- local
-    integer(I4B) :: ierr
-    integer(I4B) :: node, n
-    integer(I4B) :: inirch, inrech
-    logical :: isfound
-    logical :: supportopenclose
-    character(len=LINELENGTH) :: line
-    ! -- formats
-    character(len=*), parameter :: fmtblkerr = &
-      &"('Looking for BEGIN PERIOD iper.  Found ', a, ' instead.')"
-    character(len=*), parameter :: fmtlsp = &
-      &"(1X,/1X,'REUSING ',A,'S FROM LAST STRESS PERIOD')"
-    character(len=*), parameter :: fmtnbd = &
-      "(1X,/1X,'THE NUMBER OF ACTIVE ',A,'S (',I6, &
-      &') IS GREATER THAN MAXIMUM(',I6,')')"
-    character(len=*), parameter :: fmtdimlayered = &
-      "('When READASARRAYS is specified for the selected discretization &
-      &package, DIMENSIONS block must be omitted.')"
 ! ------------------------------------------------------------------------------
     !
-    if (this%inunit == 0) return
+    if (this%iper /= kper) return
     !
-    ! -- Set ionper to the stress period number for which a new block of data
-    !    will be read.
-    if (this%ionper < kper) then
-      !
-      ! -- get period block
-      supportopenclose = .not. this%read_as_arrays
-      ! When reading a list, OPEN/CLOSE is handled by list reader,
-      ! so supportOpenClose needs to be false in call the GetBlock.
-      ! When reading as arrays, set supportOpenClose as desired.
-      call this%parser%GetBlock('PERIOD', isfound, ierr, &
-                                blockRequired=.false.)
-      if (isfound) then
-        !
-        ! -- read ionper and check for increasing period numbers
-        call this%read_check_ionper()
-      else
-        !
-        ! -- PERIOD block not found
-        if (ierr < 0) then
-          ! -- End of file found; data applies for remainder of simulation.
-          this%ionper = nper + 1
-        else
-          ! -- Found invalid block
-          call this%parser%GetCurrentLine(line)
-          write (errmsg, fmtblkerr) adjustl(trim(line))
-          call store_error(errmsg)
-          if (this%read_as_arrays) then
-            write (errmsg, fmtdimlayered)
-            call store_error(errmsg)
-          end if
-          call this%parser%StoreErrorUnit()
-        end if
-      end if
-    end if
-    !
-    ! -- Read data if ionper == kper
-    inrech = 0
-    inirch = 0
-    if (this%ionper == kper) then
+    if (this%read_as_arrays) then
       !
-      ! -- Remove all time-series links associated with this package
-      call this%TsManager%Reset(this%packName)
-      call this%TasManager%Reset(this%packName)
-      !
-      if (.not. this%read_as_arrays) then
-        ! -- Read RECHARGE and other input as a list
-        call this%rch_rp_list(inrech)
-        call this%bnd_rp_ts()
-      else
-        ! -- Read RECHARGE, IRCH, and AUX variables as arrays
-        call this%rch_rp_array(line, inrech)
-      end if
+      ! -- update nodelist based on IRCH input
+      call nodelist_update(this%nodelist, this%nbound, this%maxbound, &
+                           this%dis, this%input_mempath)
       !
     else
-      write (this%iout, fmtlsp) trim(this%filtyp)
-    end if
-    !
-    ! -- If recharge was read, then multiply by cell area.  If inrech = 2, then
-    !    recharge is begin managed as a time series, and the time series object
-    !    will multiply the recharge rate by the cell area.
-    if (inrech == 1) then
-      do n = 1, this%nbound
-        node = this%nodelist(n)
-        if (node > 0) then
-          this%bound(1, n) = this%bound(1, n) * this%dis%get_area(node)
-        end if
-      end do
-    end if
-    !
-    ! -- return
-    return
-  end subroutine rch_rp
-
-  subroutine rch_rp_array(this, line, inrech)
-! ******************************************************************************
-! rch_rp_array -- Read and Prepare Recharge as arrays
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    use ConstantsModule, only: LENTIMESERIESNAME, LINELENGTH
-    use SimModule, only: store_error
-    use ArrayHandlersModule, only: ifind
-    implicit none
-    ! -- dummy
-    class(RchType), intent(inout) :: this
-    character(len=LINELENGTH), intent(inout) :: line
-    integer(I4B), intent(inout) :: inrech
-    ! -- local
-    integer(I4B) :: n
-    integer(I4B) :: ipos
-    integer(I4B) :: jcol, jauxcol, lpos, ivarsread
-    character(len=LENTIMESERIESNAME) :: tasName
-    character(len=24), dimension(2) :: aname
-    character(len=LINELENGTH) :: keyword, atemp
-    logical :: found, endOfBlock
-    logical :: convertFlux
-    !
-    ! -- these time array series pointers need to be non-contiguous
-    !    beacuse a slice of bound is passed
-    real(DP), dimension(:), pointer :: bndArrayPtr => null()
-    real(DP), dimension(:), pointer :: auxArrayPtr => null()
-    real(DP), dimension(:), pointer :: auxMultArray => null()
-    type(TimeArraySeriesLinkType), pointer :: tasLink => null()
-    ! -- formats
-    character(len=*), parameter :: fmtrchauxmult = &
-      "(4x, 'THE RECHARGE ARRAY IS BEING MULTIPLED BY THE AUXILIARY ARRAY WITH &
-        &THE NAME: ', A)"
-    ! -- data
-    data aname(1)/'     LAYER OR NODE INDEX'/
-    data aname(2)/'                RECHARGE'/
-    !
-! ------------------------------------------------------------------------------
-    !
-    ! -- Initialize
-    jauxcol = 0
-    ivarsread = 0
-    !
-    ! -- Read RECHARGE, IRCH, and AUX variables as arrays
-    do
-      call this%parser%GetNextLine(endOfBlock)
-      if (endOfBlock) exit
-      call this%parser%GetStringCaps(keyword)
-      !
-      ! -- Parse the keywords
-      select case (keyword)
-      case ('RECHARGE')
-        !
-        ! -- Look for keyword TIMEARRAYSERIES and time-array series
-        !    name on line, following RECHARGE
-        call this%parser%GetStringCaps(keyword)
-        if (keyword == 'TIMEARRAYSERIES') then
-          ! -- Get time-array series name
-          call this%parser%GetStringCaps(tasName)
-          jcol = 1 ! for recharge rate
-          bndArrayPtr => this%bound(jcol, :)
-          ! Make a time-array-series link and add it to the list of links
-          ! contained in the TimeArraySeriesManagerType object.
-          convertflux = .true.
-          call this%TasManager%MakeTasLink(this%packName, bndArrayPtr, &
-                                           this%iprpak, tasName, 'RECHARGE', &
-                                           convertFlux, this%nodelist, &
-                                           this%parser%iuactive)
-          lpos = this%TasManager%CountLinks()
-          tasLink => this%TasManager%GetLink(lpos)
-          inrech = 2
-        else
-          !
-          ! -- Read the recharge array, then indicate
-          !    that recharge was read by setting inrech
-          call this%dis%read_layer_array(this%nodelist, this%bound, &
-                                         this%ncolbnd, this%maxbound, 1, &
-                                         aname(2), this%parser%iuactive, &
-                                         this%iout)
-          inrech = 1
-        end if
-        !
-      case ('IRCH')
-        !
-        ! -- Check to see if other variables have already been read.  If so,
-        !    then terminate with an error that IRCH must be read first.
-        if (ivarsread > 0) then
-          call store_error('IRCH IS NOT FIRST VARIABLE IN &
-            &PERIOD BLOCK OR IT IS SPECIFIED MORE THAN ONCE.')
-          call this%parser%StoreErrorUnit()
-        end if
-        !
-        ! -- Read the IRCH array
-        call this%dis%nlarray_to_nodelist(this%nodelist, this%maxbound, &
-                                          this%nbound, aname(1), &
-                                          this%parser%iuactive, this%iout)
-        !
-        ! -- set flag to indicate that irch array has been read
-        this%inirch = 1
-        !
-        ! -- if fixed_cell option not set, then need to store nodelist
-        !    in the nodesontop array
-        if (.not. this%fixed_cell) call this%set_nodesontop()
-        !
-      case default
-        !
-        ! -- Check for auxname, and if found, then read into auxvar array
-        found = .false.
-        ipos = ifind(this%auxname, keyword)
-        if (ipos > 0) then
-          found = .true.
-          atemp = keyword
-          !
-          ! -- Look for keyword TIMEARRAYSERIES and time-array series
-          !    name on line, following auxname
-          call this%parser%GetStringCaps(keyword)
-          if (keyword == 'TIMEARRAYSERIES') then
-            ! -- Get time-array series name
-            call this%parser%GetStringCaps(tasName)
-            jauxcol = jauxcol + 1
-            auxArrayPtr => this%auxvar(jauxcol, :)
-            ! Make a time-array-series link and add it to the list of links
-            ! contained in the TimeArraySeriesManagerType object.
-            convertflux = .false.
-            call this%TasManager%MakeTasLink(this%packName, auxArrayPtr, &
-                                             this%iprpak, tasName, &
-                                             this%auxname(ipos), convertFlux, &
-                                             this%nodelist, &
-                                             this%parser%iuactive)
-          else
-            !
-            ! -- Read the aux variable array
-            call this%dis%read_layer_array(this%nodelist, this%auxvar, &
-                                           this%naux, this%maxbound, ipos, &
-                                           atemp, this%parser%iuactive, this%iout)
-          end if
-        end if
-        !
-        ! -- Nothing found
-        if (.not. found) then
-          call this%parser%GetCurrentLine(line)
-          errmsg = 'LOOKING FOR VALID VARIABLE NAME.  FOUND: '//trim(line)
-          call store_error(errmsg)
-          call this%parser%StoreErrorUnit()
-        end if
-        !
-        ! -- If this aux variable has been designated as a multiplier array
-        !    by presence of AUXMULTNAME, set local pointer appropriately.
-        if (this%iauxmultcol > 0 .and. this%iauxmultcol == ipos) then
-          auxMultArray => this%auxvar(this%iauxmultcol, :)
-        end if
-      end select
       !
-      ! -- Increment the number of variables read
-      ivarsread = ivarsread + 1
+      call this%BndExtType%bnd_rp()
       !
-    end do
-    !
-    ! -- If the multiplier-array pointer has been assigned and
-    !    stress is controlled by a time-array series, assign
-    !    multiplier-array pointer in time-array series link.
-    if (associated(auxMultArray)) then
-      if (associated(tasLink)) then
-        tasLink%RMultArray => auxMultArray
-      end if
     end if
     !
-    ! -- If recharge was read and auxmultcol was specified, then multiply
-    !    the recharge rate by the multplier column
-    if (inrech == 1 .and. this%iauxmultcol > 0) then
-      write (this%iout, fmtrchauxmult) this%auxname(this%iauxmultcol)
-      do n = 1, this%nbound
-        this%bound(this%iscloc, n) = this%bound(this%iscloc, n) * &
-                                     this%auxvar(this%iauxmultcol, n)
-      end do
-    end if
-    !
-    return
-  end subroutine rch_rp_array
-
-  subroutine rch_rp_list(this, inrech)
-! ******************************************************************************
-! rch_rp_list -- Read and Prepare Recharge as a list
-! ******************************************************************************
-!
-!    SPECIFICATIONS:
-! ------------------------------------------------------------------------------
-    implicit none
-    ! -- dummy
-    class(RchType), intent(inout) :: this
-    integer(I4B), intent(inout) :: inrech
-    ! -- local
-    integer(I4B) :: maxboundorig, nlist
-    !
-! ------------------------------------------------------------------------------
-    !
-    ! -- initialize
-    nlist = -1
-    maxboundorig = this%maxbound
-    !
-    ! -- read the list of recharge values; scale the recharge by auxmultcol
-    !    if it is specified.
-    call this%dis%read_list(this%parser%line_reader, &
-                            this%parser%iuactive, this%iout, this%iprpak, &
-                            nlist, this%inamedbound, this%iauxmultcol, &
-                            this%nodelist, this%bound, this%auxvar, &
-                            this%auxname, this%boundname, this%listlabel, &
-                            this%packName, this%tsManager, this%iscloc, &
-                            this%indxconvertflux)
-    this%nbound = nlist
-    if (this%maxbound > maxboundorig) then
-      ! -- The arrays that belong to BndType have been extended.
-      ! Now, RCH array nodesontop needs to be recreated.
-      if (associated(this%nodesontop)) then
-        deallocate (this%nodesontop)
-      end if
-    end if
+    ! -- copy nodelist to nodesontop if not fixed cell
     if (.not. this%fixed_cell) call this%set_nodesontop()
-    inrech = 1
     !
-    ! -- terminate the period block
-    call this%parser%terminateblock()
+    ! -- Write the list to iout if requested
+    if (this%iprpak /= 0) then
+      call this%write_list()
+    end if
     !
+    ! -- return
     return
-  end subroutine rch_rp_list
+  end subroutine rch_rp
 
   subroutine set_nodesontop(this)
 ! ******************************************************************************
@@ -689,7 +414,12 @@ subroutine rch_cf(this, reset_mover)
       !
       ! -- Set rhs and hcof
       this%hcof(i) = DZERO
-      this%rhs(i) = -this%bound(1, i)
+      if (this%iauxmultcol > 0) then
+        this%rhs(i) = -this%recharge(i) * this%dis%get_area(node) * &
+                      this%auxvar(this%iauxmultcol, i)
+      else
+        this%rhs(i) = -this%recharge(i) * this%dis%get_area(node)
+      end if
       if (this%ibound(node) <= 0) then
         this%rhs(i) = DZERO
         cycle
@@ -754,13 +484,15 @@ subroutine rch_da(this)
 ! ------------------------------------------------------------------------------
     !
     ! -- Deallocate parent package
-    call this%BndType%bnd_da()
+    call this%BndExtType%bnd_da()
     !
     ! -- scalars
-    call mem_deallocate(this%inirch)
+    deallocate (this%fixed_cell)
+    deallocate (this%read_as_arrays)
     !
     ! -- arrays
     if (associated(this%nodesontop)) deallocate (this%nodesontop)
+    call mem_deallocate(this%recharge, 'RECHARGE', this%memoryPath)
     !
     ! -- return
     return
@@ -810,8 +542,6 @@ subroutine default_nodelist(this)
 ! ------------------------------------------------------------------------------
     ! -- modules
     use InputOutputModule, only: get_node
-    use SimModule, only: store_error
-    use ConstantsModule, only: LINELENGTH
     ! -- dummy
     class(RchType) :: this
     ! -- local
@@ -841,8 +571,7 @@ subroutine default_nodelist(this)
       end do
     end do
     !
-    ! Set flag that indicates IRCH has been assigned, and assign nbound.
-    this%inirch = 1
+    ! -- Assign nbound
     this%nbound = ipos - 1
     !
     ! -- if fixed_cell option not set, then need to store nodelist
@@ -893,32 +622,86 @@ subroutine rch_df_obs(this)
     return
   end subroutine rch_df_obs
 
-  !
-  ! -- Procedure related to time series
-  subroutine rch_rp_ts(this)
-    ! -- Assign tsLink%Text appropriately for
-    !    all time series in use by package.
-    !    In RCH package only the RECHARGE variable
-    !    can be controlled by time series.
+! ******************************************************************************
+! rch_bound_value -- return requested boundary value
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+  function rch_bound_value(this, col, row) result(bndval)
+    ! -- modules
+    use ConstantsModule, only: DZERO
+    ! -- dummy variables
+    class(RchType), intent(inout) :: this !< BndExtType object
+    integer(I4B), intent(in) :: col
+    integer(I4B), intent(in) :: row
+    ! -- result
+    real(DP) :: bndval
+    !
+    select case (col)
+    case (1)
+      if (this%iauxmultcol > 0) then
+        bndval = this%recharge(row) * this%auxvar(this%iauxmultcol, row)
+      else
+        bndval = this%recharge(row)
+      end if
+    case default
+      errmsg = 'Programming error. RCH bound value requested column '&
+               &'outside range of ncolbnd (1).'
+      call store_error(errmsg)
+      call store_error_filename(this%input_fname)
+    end select
+    !
+    ! -- return
+    return
+  end function rch_bound_value
+
+  !> @brief Update the nodelist based on IRCH input
+  !!
+  !! This is a module scoped routine to check for IRCH
+  !! input. If array input was provided, INIRCH and IRCH
+  !! will be allocated in the input context.  If the read
+  !! state variable INIRCH is set to 1 during this period
+  !! update, IRCH input was read and is used here to update
+  !! the nodelist.
+  !!
+  !<
+  subroutine nodelist_update(nodelist, nbound, maxbound, &
+                             dis, input_mempath)
+    ! -- modules
+    use MemoryManagerModule, only: mem_setptr
+    use BaseDisModule, only: DisBaseType
     ! -- dummy
-    class(RchType), intent(inout) :: this
+    integer(I4B), dimension(:), contiguous, &
+      pointer, intent(inout) :: nodelist
+    class(DisBaseType), pointer, intent(in) :: dis
+    character(len=*), intent(in) :: input_mempath
+    integer(I4B), intent(inout) :: nbound
+    integer(I4B), intent(in) :: maxbound
+    character(len=24) :: aname = '     LAYER OR NODE INDEX'
     ! -- local
-    integer(I4B) :: i, nlinks
-    type(TimeSeriesLinkType), pointer :: tslink => null()
-    !
-    nlinks = this%TsManager%boundtslinks%Count()
-    do i = 1, nlinks
-      tslink => GetTimeSeriesLinkFromList(this%TsManager%boundtslinks, i)
-      if (associated(tslink)) then
-        select case (tslink%JCol)
-        case (1)
-          tslink%Text = 'RECHARGE'
-        end select
-      end if
-    end do
+    integer(I4B), dimension(:), contiguous, &
+      pointer :: irch => null()
+    integer(I4B), pointer :: inirch => NULL()
     !
+    ! -- set pointer to input context INIRCH
+    call mem_setptr(inirch, 'INIRCH', input_mempath)
+    !
+    ! -- check INIRCH read state
+    if (inirch == 1) then
+      ! -- irch was read this period
+      !
+      ! -- set pointer to input context IRCH
+      call mem_setptr(irch, 'IRCH', input_mempath)
+      !
+      ! -- update nodelist
+      call dis%nlarray_to_nodelist2(irch, nodelist, &
+                                    maxbound, nbound, aname)
+    end if
+    !
+    ! -- return
     return
-  end subroutine rch_rp_ts
+  end subroutine nodelist_update
 
 end module RchModule
 
diff --git a/src/Model/GroundWaterFlow/gwf3rch8idm.f90 b/src/Model/GroundWaterFlow/gwf3rch8idm.f90
new file mode 100644
index 00000000000..89b0d51869e
--- /dev/null
+++ b/src/Model/GroundWaterFlow/gwf3rch8idm.f90
@@ -0,0 +1,433 @@
+! ** Do Not Modify! MODFLOW 6 system generated file. **
+module GwfRchInputModule
+  use ConstantsModule, only: LENVARNAME
+  use InputDefinitionModule, only: InputParamDefinitionType, &
+                                   InputBlockDefinitionType
+  private
+  public gwf_rch_param_definitions
+  public gwf_rch_aggregate_definitions
+  public gwf_rch_block_definitions
+  public GwfRchParamFoundType
+  public gwf_rch_multi_package
+  public gwf_rch_aux_sfac_param
+
+  type GwfRchParamFoundType
+    logical :: fixed_cell = .false.
+    logical :: auxiliary = .false.
+    logical :: auxmultname = .false.
+    logical :: boundnames = .false.
+    logical :: iprpak = .false.
+    logical :: iprflow = .false.
+    logical :: ipakcb = .false.
+    logical :: ts_filerecord = .false.
+    logical :: ts6 = .false.
+    logical :: filein = .false.
+    logical :: ts6_filename = .false.
+    logical :: obs_filerecord = .false.
+    logical :: obs6 = .false.
+    logical :: obs6_filename = .false.
+    logical :: maxbound = .false.
+    logical :: cellid = .false.
+    logical :: recharge = .false.
+    logical :: auxvar = .false.
+    logical :: boundname = .false.
+  end type GwfRchParamFoundType
+
+  logical :: gwf_rch_multi_package = .true.
+
+  character(len=LENVARNAME) :: gwf_rch_aux_sfac_param = 'RECHARGE'
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_fixed_cell = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'FIXED_CELL', & ! tag name
+    'FIXED_CELL', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_auxiliary = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'AUXILIARY', & ! tag name
+    'AUXILIARY', & ! fortran variable
+    'STRING', & ! type
+    'NAUX', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_auxmultname = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'AUXMULTNAME', & ! tag name
+    'AUXMULTNAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_boundnames = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'BOUNDNAMES', & ! tag name
+    'BOUNDNAMES', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_iprpak = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'PRINT_INPUT', & ! tag name
+    'IPRPAK', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_iprflow = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'PRINT_FLOWS', & ! tag name
+    'IPRFLOW', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_ipakcb = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'SAVE_FLOWS', & ! tag name
+    'IPAKCB', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_ts_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'TS_FILERECORD', & ! tag name
+    'TS_FILERECORD', & ! fortran variable
+    'RECORD TS6 FILEIN TS6_FILENAME', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_ts6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'TS6', & ! tag name
+    'TS6', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_filein = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'FILEIN', & ! tag name
+    'FILEIN', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_ts6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'TS6_FILENAME', & ! tag name
+    'TS6_FILENAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .true., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_obs_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS_FILERECORD', & ! tag name
+    'OBS_FILERECORD', & ! fortran variable
+    'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_obs6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS6', & ! tag name
+    'OBS6', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_obs6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS6_FILENAME', & ! tag name
+    'OBS6_FILENAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .true., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_maxbound = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'DIMENSIONS', & ! block
+    'MAXBOUND', & ! tag name
+    'MAXBOUND', & ! fortran variable
+    'INTEGER', & ! type
+    '', & ! shape
+    .true., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_cellid = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'PERIOD', & ! block
+    'CELLID', & ! tag name
+    'CELLID', & ! fortran variable
+    'INTEGER1D', & ! type
+    'NCELLDIM', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_recharge = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'PERIOD', & ! block
+    'RECHARGE', & ! tag name
+    'RECHARGE', & ! fortran variable
+    'DOUBLE', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_auxvar = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'PERIOD', & ! block
+    'AUX', & ! tag name
+    'AUXVAR', & ! fortran variable
+    'DOUBLE1D', & ! type
+    'NAUX', & ! shape
+    .false., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_boundname = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'PERIOD', & ! block
+    'BOUNDNAME', & ! tag name
+    'BOUNDNAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .false., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwf_rch_param_definitions(*) = &
+    [ &
+    gwfrch_fixed_cell, &
+    gwfrch_auxiliary, &
+    gwfrch_auxmultname, &
+    gwfrch_boundnames, &
+    gwfrch_iprpak, &
+    gwfrch_iprflow, &
+    gwfrch_ipakcb, &
+    gwfrch_ts_filerecord, &
+    gwfrch_ts6, &
+    gwfrch_filein, &
+    gwfrch_ts6_filename, &
+    gwfrch_obs_filerecord, &
+    gwfrch_obs6, &
+    gwfrch_obs6_filename, &
+    gwfrch_maxbound, &
+    gwfrch_cellid, &
+    gwfrch_recharge, &
+    gwfrch_auxvar, &
+    gwfrch_boundname &
+    ]
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrch_spd = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCH', & ! subcomponent
+    'PERIOD', & ! block
+    'STRESS_PERIOD_DATA', & ! tag name
+    'SPD', & ! fortran variable
+    'RECARRAY CELLID RECHARGE AUX BOUNDNAME', & ! type
+    'MAXBOUND', & ! shape
+    .true., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwf_rch_aggregate_definitions(*) = &
+    [ &
+    gwfrch_spd &
+    ]
+
+  type(InputBlockDefinitionType), parameter :: &
+    gwf_rch_block_definitions(*) = &
+    [ &
+    InputBlockDefinitionType( &
+    'OPTIONS', & ! blockname
+    .false., & ! required
+    .false., & ! aggregate
+    .false. & ! block_variable
+    ), &
+    InputBlockDefinitionType( &
+    'DIMENSIONS', & ! blockname
+    .true., & ! required
+    .false., & ! aggregate
+    .false. & ! block_variable
+    ), &
+    InputBlockDefinitionType( &
+    'PERIOD', & ! blockname
+    .true., & ! required
+    .true., & ! aggregate
+    .true. & ! block_variable
+    ) &
+    ]
+
+end module GwfRchInputModule
diff --git a/src/Model/GroundWaterFlow/gwf3rcha8idm.f90 b/src/Model/GroundWaterFlow/gwf3rcha8idm.f90
new file mode 100644
index 00000000000..ee044f0fcc8
--- /dev/null
+++ b/src/Model/GroundWaterFlow/gwf3rcha8idm.f90
@@ -0,0 +1,386 @@
+! ** Do Not Modify! MODFLOW 6 system generated file. **
+module GwfRchaInputModule
+  use ConstantsModule, only: LENVARNAME
+  use InputDefinitionModule, only: InputParamDefinitionType, &
+                                   InputBlockDefinitionType
+  private
+  public gwf_rcha_param_definitions
+  public gwf_rcha_aggregate_definitions
+  public gwf_rcha_block_definitions
+  public GwfRchaParamFoundType
+  public gwf_rcha_multi_package
+  public gwf_rcha_aux_sfac_param
+
+  type GwfRchaParamFoundType
+    logical :: readasarrays = .false.
+    logical :: fixed_cell = .false.
+    logical :: auxiliary = .false.
+    logical :: auxmultname = .false.
+    logical :: iprpak = .false.
+    logical :: iprflow = .false.
+    logical :: ipakcb = .false.
+    logical :: tas_filerecord = .false.
+    logical :: tas6 = .false.
+    logical :: filein = .false.
+    logical :: tas6_filename = .false.
+    logical :: obs_filerecord = .false.
+    logical :: obs6 = .false.
+    logical :: obs6_filename = .false.
+    logical :: irch = .false.
+    logical :: recharge = .false.
+    logical :: auxvar = .false.
+  end type GwfRchaParamFoundType
+
+  logical :: gwf_rcha_multi_package = .true.
+
+  character(len=LENVARNAME) :: gwf_rcha_aux_sfac_param = 'RECHARGE'
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_readasarrays = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'READASARRAYS', & ! tag name
+    'READASARRAYS', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_fixed_cell = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'FIXED_CELL', & ! tag name
+    'FIXED_CELL', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_auxiliary = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'AUXILIARY', & ! tag name
+    'AUXILIARY', & ! fortran variable
+    'STRING', & ! type
+    'NAUX', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_auxmultname = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'AUXMULTNAME', & ! tag name
+    'AUXMULTNAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_iprpak = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'PRINT_INPUT', & ! tag name
+    'IPRPAK', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_iprflow = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'PRINT_FLOWS', & ! tag name
+    'IPRFLOW', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_ipakcb = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'SAVE_FLOWS', & ! tag name
+    'IPAKCB', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_tas_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'TAS_FILERECORD', & ! tag name
+    'TAS_FILERECORD', & ! fortran variable
+    'RECORD TAS6 FILEIN TAS6_FILENAME', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_tas6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'TAS6', & ! tag name
+    'TAS6', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_filein = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'FILEIN', & ! tag name
+    'FILEIN', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_tas6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'TAS6_FILENAME', & ! tag name
+    'TAS6_FILENAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .true., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_obs_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS_FILERECORD', & ! tag name
+    'OBS_FILERECORD', & ! fortran variable
+    'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_obs6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS6', & ! tag name
+    'OBS6', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_obs6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'OPTIONS', & ! block
+    'OBS6_FILENAME', & ! tag name
+    'OBS6_FILENAME', & ! fortran variable
+    'STRING', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .true., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_irch = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'PERIOD', & ! block
+    'IRCH', & ! tag name
+    'IRCH', & ! fortran variable
+    'INTEGER1D', & ! type
+    'NCPL', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_recharge = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'PERIOD', & ! block
+    'RECHARGE', & ! tag name
+    'RECHARGE', & ! fortran variable
+    'DOUBLE1D', & ! type
+    'NCPL', & ! shape
+    .true., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfrcha_auxvar = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RCHA', & ! subcomponent
+    'PERIOD', & ! block
+    'AUX', & ! tag name
+    'AUXVAR', & ! fortran variable
+    'DOUBLE2D', & ! type
+    'NAUX NCPL', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwf_rcha_param_definitions(*) = &
+    [ &
+    gwfrcha_readasarrays, &
+    gwfrcha_fixed_cell, &
+    gwfrcha_auxiliary, &
+    gwfrcha_auxmultname, &
+    gwfrcha_iprpak, &
+    gwfrcha_iprflow, &
+    gwfrcha_ipakcb, &
+    gwfrcha_tas_filerecord, &
+    gwfrcha_tas6, &
+    gwfrcha_filein, &
+    gwfrcha_tas6_filename, &
+    gwfrcha_obs_filerecord, &
+    gwfrcha_obs6, &
+    gwfrcha_obs6_filename, &
+    gwfrcha_irch, &
+    gwfrcha_recharge, &
+    gwfrcha_auxvar &
+    ]
+
+  type(InputParamDefinitionType), parameter :: &
+    gwf_rcha_aggregate_definitions(*) = &
+    [ &
+    InputParamDefinitionType &
+    ( &
+    '', & ! component
+    '', & ! subcomponent
+    '', & ! block
+    '', & ! tag name
+    '', & ! fortran variable
+    '', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    ) &
+    ]
+
+  type(InputBlockDefinitionType), parameter :: &
+    gwf_rcha_block_definitions(*) = &
+    [ &
+    InputBlockDefinitionType( &
+    'OPTIONS', & ! blockname
+    .true., & ! required
+    .false., & ! aggregate
+    .false. & ! block_variable
+    ), &
+    InputBlockDefinitionType( &
+    'PERIOD', & ! blockname
+    .true., & ! required
+    .false., & ! aggregate
+    .true. & ! block_variable
+    ) &
+    ]
+
+end module GwfRchaInputModule
diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90
index ab19fd6afe5..78bea23b746 100644
--- a/src/Model/ModelUtilities/DiscretizationBase.f90
+++ b/src/Model/ModelUtilities/DiscretizationBase.f90
@@ -105,6 +105,7 @@ module BaseDisModule
     procedure, private :: record_srcdst_list_entry
     generic, public :: record_mf6_list_entry => record_srcdst_list_entry
     procedure, public :: nlarray_to_nodelist
+    procedure, public :: nlarray_to_nodelist2
     procedure, public :: highest_active
     procedure, public :: get_area
     procedure, public :: get_area_factor
@@ -1443,6 +1444,35 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, &
     return
   end subroutine nlarray_to_nodelist
 
+  subroutine nlarray_to_nodelist2(this, darray, nodelist, maxbnd, nbound, aname)
+! ******************************************************************************
+! nlarray_to_nodelist -- Read an integer array into nodelist. For structured
+!                        model, integer array is layer number; for unstructured
+!                        model, integer array is node number.
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use SimModule, only: store_error
+    use ConstantsModule, only: LINELENGTH
+    ! -- dummy
+    class(DisBaseType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    !
+    ! --
+    errmsg = 'Programmer error: nlarray_to_nodelist needs to be &
+            &overridden in any DIS type that extends DisBaseType'
+    call store_error(errmsg, terminate=.TRUE.)
+    !
+    ! -- return
+    return
+  end subroutine nlarray_to_nodelist2
+
   subroutine highest_active(this, n, ibound)
 ! ******************************************************************************
 ! highest_active -- Find the first highest active cell beneath cell n
diff --git a/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90 b/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
index 386ca2b8726..5bac4554e36 100644
--- a/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
+++ b/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
@@ -12,6 +12,8 @@ module IdmGwfDfnSelectorModule
   use GwfDrnInputModule
   use GwfGhbInputModule
   use GwfNpfInputModule
+  use GwfRchInputModule
+  use GwfRchaInputModule
   use GwfRivInputModule
   use GwfWelInputModule
   use GwfNamInputModule
@@ -119,6 +121,13 @@ module IdmGwfDfnSelectorModule
     logical :: angle2 = .false.
     logical :: angle3 = .false.
     logical :: wetdry = .false.
+    logical :: fixed_cell = .false.
+    logical :: recharge = .false.
+    logical :: readasarrays = .false.
+    logical :: tas_filerecord = .false.
+    logical :: tas6 = .false.
+    logical :: tas6_filename = .false.
+    logical :: irch = .false.
     logical :: stage = .false.
     logical :: rbot = .false.
     logical :: flowred = .false.
@@ -172,6 +181,10 @@ function gwf_param_definitions(subcomponent) result(input_definition)
       call set_param_pointer(input_definition, gwf_ghb_param_definitions)
     case ('NPF')
       call set_param_pointer(input_definition, gwf_npf_param_definitions)
+    case ('RCH')
+      call set_param_pointer(input_definition, gwf_rch_param_definitions)
+    case ('RCHA')
+      call set_param_pointer(input_definition, gwf_rcha_param_definitions)
     case ('RIV')
       call set_param_pointer(input_definition, gwf_riv_param_definitions)
     case ('WEL')
@@ -202,6 +215,10 @@ function gwf_aggregate_definitions(subcomponent) result(input_definition)
       call set_param_pointer(input_definition, gwf_ghb_aggregate_definitions)
     case ('NPF')
       call set_param_pointer(input_definition, gwf_npf_aggregate_definitions)
+    case ('RCH')
+      call set_param_pointer(input_definition, gwf_rch_aggregate_definitions)
+    case ('RCHA')
+      call set_param_pointer(input_definition, gwf_rcha_aggregate_definitions)
     case ('RIV')
       call set_param_pointer(input_definition, gwf_riv_aggregate_definitions)
     case ('WEL')
@@ -232,6 +249,10 @@ function gwf_block_definitions(subcomponent) result(input_definition)
       call set_block_pointer(input_definition, gwf_ghb_block_definitions)
     case ('NPF')
       call set_block_pointer(input_definition, gwf_npf_block_definitions)
+    case ('RCH')
+      call set_block_pointer(input_definition, gwf_rch_block_definitions)
+    case ('RCHA')
+      call set_block_pointer(input_definition, gwf_rcha_block_definitions)
     case ('RIV')
       call set_block_pointer(input_definition, gwf_riv_block_definitions)
     case ('WEL')
@@ -261,6 +282,10 @@ function gwf_idm_multi_package(subcomponent) result(multi_package)
       multi_package = gwf_ghb_multi_package
     case ('NPF')
       multi_package = gwf_npf_multi_package
+    case ('RCH')
+      multi_package = gwf_rch_multi_package
+    case ('RCHA')
+      multi_package = gwf_rcha_multi_package
     case ('RIV')
       multi_package = gwf_riv_multi_package
     case ('WEL')
@@ -293,6 +318,10 @@ function gwf_idm_sfac_param(subcomponent) result(sfac_param)
       sfac_param = gwf_ghb_aux_sfac_param
     case ('NPF')
       sfac_param = gwf_npf_aux_sfac_param
+    case ('RCH')
+      sfac_param = gwf_rch_aux_sfac_param
+    case ('RCHA')
+      sfac_param = gwf_rcha_aux_sfac_param
     case ('RIV')
       sfac_param = gwf_riv_aux_sfac_param
     case ('WEL')
@@ -326,6 +355,10 @@ function gwf_idm_integrated(subcomponent) result(integrated)
       integrated = .true.
     case ('NPF')
       integrated = .true.
+    case ('RCH')
+      integrated = .true.
+    case ('RCHA')
+      integrated = .true.
     case ('RIV')
       integrated = .true.
     case ('WEL')
diff --git a/src/meson.build b/src/meson.build
index c3533a7e600..731b30ad3fe 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -74,6 +74,8 @@ modflow_sources = files(
     'Model' / 'GroundWaterFlow' / 'gwf3obs8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3oc8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3rch8.f90',
+    'Model' / 'GroundWaterFlow' / 'gwf3rch8idm.f90',
+    'Model' / 'GroundWaterFlow' / 'gwf3rcha8idm.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3riv8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3riv8idm.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3sfr8.f90',
diff --git a/utils/idmloader/scripts/dfn2f90.py b/utils/idmloader/scripts/dfn2f90.py
index fefe869e865..cc27bf78010 100644
--- a/utils/idmloader/scripts/dfn2f90.py
+++ b/utils/idmloader/scripts/dfn2f90.py
@@ -996,6 +996,14 @@ def _write_master_component(self, fh=None):
             Path("../../../doc/mf6io/mf6ivar/dfn", "gwf-npf.dfn"),
             Path("../../../src/Model/GroundWaterFlow", "gwf3npf8idm.f90"),
         ],
+        [
+            Path("../../../doc/mf6io/mf6ivar/dfn", "gwf-rch.dfn"),
+            Path("../../../src/Model/GroundWaterFlow", "gwf3rch8idm.f90"),
+        ],
+        [
+            Path("../../../doc/mf6io/mf6ivar/dfn", "gwf-rcha.dfn"),
+            Path("../../../src/Model/GroundWaterFlow", "gwf3rcha8idm.f90"),
+        ],
         [
             Path("../../../doc/mf6io/mf6ivar/dfn", "gwf-riv.dfn"),
             Path("../../../src/Model/GroundWaterFlow", "gwf3riv8idm.f90"),