diff --git a/autotest/test_gwf_libmf6_riv01.py b/autotest/test_gwf_libmf6_riv01.py
index 0fc8dd7e598..0ace1a457f4 100644
--- a/autotest/test_gwf_libmf6_riv01.py
+++ b/autotest/test_gwf_libmf6_riv01.py
@@ -178,8 +178,12 @@ def api_func(exe, idx, model_ws=None):
     end_time = mf6.get_end_time()
 
     # get copy of (multi-dim) array with river parameters
-    riv_tag = mf6.get_var_address("BOUND", name, riv_packname)
-    new_spd = mf6.get_value(riv_tag)
+    stage_tag = mf6.get_var_address("STAGE", name, riv_packname)
+    cond_tag = mf6.get_var_address("COND", name, riv_packname)
+    rbot_tag = mf6.get_var_address("RBOT", name, riv_packname)
+    new_stage = mf6.get_value(stage_tag)
+    new_cond = mf6.get_value(cond_tag)
+    new_rbot = mf6.get_value(rbot_tag)
 
     # model time loop
     idx = 0
@@ -192,16 +196,18 @@ def api_func(exe, idx, model_ws=None):
         mf6.prepare_time_step(dt)
 
         # set the RIV data through the BMI
+        # change cond and rbot data
+        new_cond[:] = [riv_cond]
+        new_rbot[:] = [riv_bot]
+        mf6.set_value(cond_tag, new_cond)
+        mf6.set_value(rbot_tag, new_rbot)
+        # change stage data
         if current_time < 5:
-            # set columns of BOUND data (we're setting entire columns of the
-            # 2D array for convenience, setting only the value for the active
-            # stress period should work too)
-            new_spd[:] = [riv_stage, riv_cond, riv_bot]
-            mf6.set_value(riv_tag, new_spd)
+            new_stage[:] = [riv_stage]
+            mf6.set_value(stage_tag, new_stage)
         else:
-            # change only stage data
-            new_spd[:] = [riv_stage2, riv_cond, riv_bot]
-            mf6.set_value(riv_tag, new_spd)
+            new_stage[:] = [riv_stage2]
+            mf6.set_value(stage_tag, new_stage)
 
         kiter = 0
         mf6.prepare_solve()
diff --git a/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn b/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn
index a57b653cc0c..385f3f7b793 100644
--- a/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn
+++ b/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn
@@ -35,6 +35,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'river'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -43,6 +44,7 @@ reader urword
 optional true
 longname print calculated flows to listing file
 description REPLACE print_flows {'{#1}': 'river'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -51,6 +53,7 @@ reader urword
 optional true
 longname save CHD flows to budget file
 description REPLACE save_flows {'{#1}': 'river'}
+mf6internal ipakcb
 
 block options
 name ts_filerecord
@@ -169,6 +172,7 @@ shape (maxbound)
 reader urword
 longname
 description
+mf6internal spd
 
 block period
 name cellid
@@ -224,6 +228,7 @@ optional true
 time_series true
 longname auxiliary variables
 description REPLACE aux {'{#1}': 'river'}
+mf6internal auxvar
 
 block period
 name boundname
diff --git a/make/makefile b/make/makefile
index 64667ddf656..bf255d76c63 100644
--- a/make/makefile
+++ b/make/makefile
@@ -100,6 +100,7 @@ $(OBJDIR)/gwt1disv1idm.o \
 $(OBJDIR)/gwt1disu1idm.o \
 $(OBJDIR)/gwt1dis1idm.o \
 $(OBJDIR)/gwf3wel8idm.o \
+$(OBJDIR)/gwf3riv8idm.o \
 $(OBJDIR)/gwf3npf8idm.o \
 $(OBJDIR)/gwf3idm.o \
 $(OBJDIR)/gwf3drn8idm.o \
diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
index cf6c91bc6eb..557e437198f 100644
--- a/msvs/mf6core.vfproj
+++ b/msvs/mf6core.vfproj
@@ -146,6 +146,7 @@
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3oc8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sfr8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sto8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3tvbase8.f90"/>
diff --git a/src/Model/GroundWaterFlow/gwf3.f90 b/src/Model/GroundWaterFlow/gwf3.f90
index d955ae74b10..d57b12ab3bb 100644
--- a/src/Model/GroundWaterFlow/gwf3.f90
+++ b/src/Model/GroundWaterFlow/gwf3.f90
@@ -1299,7 +1299,8 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
       call drn_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
                       pakname, mempath)
     case ('RIV6')
-      call riv_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
+      call riv_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
+                      pakname, mempath)
     case ('GHB6')
       call ghb_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
     case ('RCH6')
diff --git a/src/Model/GroundWaterFlow/gwf3buy8.f90 b/src/Model/GroundWaterFlow/gwf3buy8.f90
index 2d2e26ef628..b30c789a9bc 100644
--- a/src/Model/GroundWaterFlow/gwf3buy8.f90
+++ b/src/Model/GroundWaterFlow/gwf3buy8.f90
@@ -648,6 +648,7 @@ subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, &
 ! ------------------------------------------------------------------------------
     ! -- modules
     use BndModule, only: BndType
+    use RivModule, only: RivType
     class(BndType), pointer :: packobj
     ! -- dummy
     real(DP), intent(in), dimension(:) :: hnew
@@ -674,39 +675,42 @@ subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, &
 ! ------------------------------------------------------------------------------
     !
     ! -- Process density terms for each RIV
-    do n = 1, packobj%nbound
-      node = packobj%nodelist(n)
-      if (packobj%ibound(node) <= 0) cycle
-      !
-      ! -- density
-      denseriv = get_bnd_density(n, locdense, locconc, denseref, &
-                                 drhodc, crhoref, ctemp, packobj%auxvar)
-      !
-      ! -- elevation
-      elevriv = elev(node)
-      if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
-      !
-      ! -- boundary head and conductance
-      hriv = packobj%bound(1, n)
-      cond = packobj%bound(2, n)
-      rbot = packobj%bound(3, n)
-      !
-      ! -- calculate and add terms depending on whether head is above rbot
-      if (hnew(node) > rbot) then
+    select type (packobj)
+    type is (RivType)
+      do n = 1, packobj%nbound
+        node = packobj%nodelist(n)
+        if (packobj%ibound(node) <= 0) cycle
         !
-        ! --calculate HCOF and RHS terms, similar to GHB in this case
-        call calc_ghb_hcof_rhs_terms(denseref, denseriv, dense(node), &
-                                     elevriv, elev(node), hriv, hnew(node), &
-                                     cond, iform, rhsterm, hcofterm)
-      else
-        hcofterm = DZERO
-        rhsterm = cond * (denseriv / denseref - DONE) * (hriv - rbot)
-      end if
-      !
-      ! -- Add terms to package hcof and rhs accumulators
-      packobj%hcof(n) = packobj%hcof(n) + hcofterm
-      packobj%rhs(n) = packobj%rhs(n) - rhsterm
-    end do
+        ! -- density
+        denseriv = get_bnd_density(n, locdense, locconc, denseref, &
+                                   drhodc, crhoref, ctemp, packobj%auxvar)
+        !
+        ! -- elevation
+        elevriv = elev(node)
+        if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
+        !
+        ! -- boundary head and conductance
+        hriv = packobj%stage(n)
+        cond = packobj%cond(n)
+        rbot = packobj%rbot(n)
+        !
+        ! -- calculate and add terms depending on whether head is above rbot
+        if (hnew(node) > rbot) then
+          !
+          ! --calculate HCOF and RHS terms, similar to GHB in this case
+          call calc_ghb_hcof_rhs_terms(denseref, denseriv, dense(node), &
+                                       elevriv, elev(node), hriv, hnew(node), &
+                                       cond, iform, rhsterm, hcofterm)
+        else
+          hcofterm = DZERO
+          rhsterm = cond * (denseriv / denseref - DONE) * (hriv - rbot)
+        end if
+        !
+        ! -- Add terms to package hcof and rhs accumulators
+        packobj%hcof(n) = packobj%hcof(n) + hcofterm
+        packobj%rhs(n) = packobj%rhs(n) - rhsterm
+      end do
+    end select
     !
     ! -- Return
     return
diff --git a/src/Model/GroundWaterFlow/gwf3riv8.f90 b/src/Model/GroundWaterFlow/gwf3riv8.f90
index 8feffb43f43..7033f91c3ac 100644
--- a/src/Model/GroundWaterFlow/gwf3riv8.f90
+++ b/src/Model/GroundWaterFlow/gwf3riv8.f90
@@ -1,11 +1,12 @@
 module rivmodule
   use KindModule, only: DP, I4B
   use ConstantsModule, only: DZERO, LENFTYPE, LENPACKAGENAME
+  use SimVariablesModule, only: errmsg
+  use SimModule, only: count_errors, store_error, store_error_filename
   use MemoryHelperModule, only: create_mem_path
   use BndModule, only: BndType
+  use BndExtModule, only: BndExtType
   use ObsModule, only: DefaultObsIdProcessor
-  use TimeSeriesLinkModule, only: TimeSeriesLinkType, &
-                                  GetTimeSeriesLinkFromList
   use MatrixBaseModule
   !
   implicit none
@@ -17,23 +18,32 @@ module rivmodule
   character(len=LENFTYPE) :: ftype = 'RIV'
   character(len=LENPACKAGENAME) :: text = '             RIV'
   !
-  type, extends(BndType) :: RivType
+  type, extends(BndExtType) :: RivType
+    real(DP), dimension(:), pointer, contiguous :: stage => null() !< RIV head
+    real(DP), dimension(:), pointer, contiguous :: cond => null() !< RIV bed hydraulic conductance
+    real(DP), dimension(:), pointer, contiguous :: rbot => null() !< RIV bed bottom elevation
   contains
-    procedure :: bnd_options => riv_options
+    procedure :: allocate_arrays => riv_allocate_arrays
+    procedure :: source_options => riv_options
+    procedure :: log_riv_options
+    procedure :: bnd_rp => riv_rp
     procedure :: bnd_ck => riv_ck
     procedure :: bnd_cf => riv_cf
     procedure :: bnd_fc => riv_fc
+    procedure :: bnd_da => riv_da
     procedure :: define_listlabel
+    procedure :: bound_value => riv_bound_value
+    procedure :: cond_mult
     ! -- methods for observations
     procedure, public :: bnd_obs_supported => riv_obs_supported
     procedure, public :: bnd_df_obs => riv_df_obs
-    ! -- method for time series
-    procedure, public :: bnd_rp_ts => riv_rp_ts
+    procedure, public :: riv_store_user_cond
   end type RivType
 
 contains
 
-  subroutine riv_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
+  subroutine riv_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, &
+                        mempath)
 ! ******************************************************************************
 ! riv_create -- Create a New Riv Package
 ! Subroutine: (1) create new-style package
@@ -50,6 +60,7 @@ subroutine riv_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(RivType), pointer :: rivobj
 ! ------------------------------------------------------------------------------
@@ -59,7 +70,7 @@ subroutine riv_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
     packobj => rivobj
     !
     ! -- 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
@@ -80,38 +91,157 @@ subroutine riv_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
     return
   end subroutine riv_create
 
-  subroutine riv_options(this, option, found)
+  subroutine riv_da(this)
+! ******************************************************************************
+! riv_da -- deallocate
 ! ******************************************************************************
-! riv_options -- set options specific to RivType
 !
-! riv_options overrides BndType%bnd_options
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use MemoryManagerModule, only: mem_deallocate
+    ! -- dummy
+    class(RivType) :: this
+! ------------------------------------------------------------------------------
+    !
+    ! -- Deallocate parent package
+    call this%BndExtType%bnd_da()
+    !
+    ! -- arrays
+    call mem_deallocate(this%stage, 'STAGE', this%memoryPath)
+    call mem_deallocate(this%cond, 'COND', this%memoryPath)
+    call mem_deallocate(this%rbot, 'RBOT', this%memoryPath)
+    !
+    ! -- return
+    return
+  end subroutine riv_da
+
+  subroutine riv_options(this)
+! ******************************************************************************
+! riv_options -- set options specific to RivType
 ! ******************************************************************************
 !
 !    SPECIFICATIONS:
 ! ------------------------------------------------------------------------------
-    use InputOutputModule, only: urword
+    use MemoryManagerExtModule, only: mem_set_value
+    use CharacterStringModule, only: CharacterStringType
+    use GwfRivInputModule, only: GwfRivParamFoundType
     ! -- dummy
     class(RivType), intent(inout) :: this
-    character(len=*), intent(inout) :: option
-    logical, intent(inout) :: found
     ! -- local
+    type(GwfRivParamFoundType) :: found
 ! ------------------------------------------------------------------------------
     !
-    select case (option)
-    case ('MOVER')
-      this%imover = 1
-      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
-      found = .true.
-    case default
-      !
-      ! -- No options found
-      found = .false.
-    end select
+    ! -- source base class options
+    call this%BndExtType%source_options()
+    !
+    ! -- source options from input context
+    call mem_set_value(this%imover, 'MOVER', this%input_mempath, found%mover)
+    !
+    ! -- log riv specific options
+    call this%log_riv_options(found)
     !
     ! -- return
     return
   end subroutine riv_options
 
+  subroutine log_riv_options(this, found)
+! ******************************************************************************
+! log_riv_options -- log options specific to RivType
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use GwfRivInputModule, only: GwfRivParamFoundType
+    ! -- dummy variables
+    class(RivType), intent(inout) :: this !< BndExtType object
+    type(GwfRivParamFoundType), intent(in) :: found
+    ! -- local variables
+    ! -- format
+    !
+    ! -- log found options
+    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%text)) &
+      //' OPTIONS'
+    !
+    if (found%mover) then
+      write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
+    end if
+    !
+    ! -- close logging block
+    write (this%iout, '(1x,a)') &
+      'END OF '//trim(adjustl(this%text))//' OPTIONS'
+    !
+    ! -- return
+    return
+  end subroutine log_riv_options
+
+  subroutine riv_allocate_arrays(this, nodelist, auxvar)
+! ******************************************************************************
+! riv_allocate_arrays -- allocate arrays
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use MemoryManagerModule, only: mem_allocate, mem_setptr, mem_checkin
+    ! -- dummy
+    class(RivType) :: this
+    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
+    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
+    ! -- local
+! ------------------------------------------------------------------------------
+    !
+    ! -- call base type allocate arrays
+    call this%BndExtType%allocate_arrays(nodelist, auxvar)
+    !
+    ! -- set riv input context pointers
+    call mem_setptr(this%stage, 'STAGE', this%input_mempath)
+    call mem_setptr(this%cond, 'COND', this%input_mempath)
+    call mem_setptr(this%rbot, 'RBOT', this%input_mempath)
+    !
+    ! --checkin riv input context pointers
+    call mem_checkin(this%stage, 'STAGE', this%memoryPath, &
+                     'STAGE', this%input_mempath)
+    call mem_checkin(this%cond, 'COND', this%memoryPath, &
+                     'COND', this%input_mempath)
+    call mem_checkin(this%rbot, 'RBOT', this%memoryPath, &
+                     'RBOT', this%input_mempath)
+    !
+    ! -- return
+    return
+  end subroutine riv_allocate_arrays
+
+  subroutine riv_rp(this)
+! ******************************************************************************
+! riv_rp -- Read and prepare
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    use TdisModule, only: kper
+    ! -- dummy
+    class(RivType), intent(inout) :: this
+    ! -- local
+! ------------------------------------------------------------------------------
+    if (this%iper /= kper) return
+    !
+    ! -- Call the parent class read and prepare
+    call this%BndExtType%bnd_rp()
+    !
+    ! -- store user cond
+    if (this%ivsc == 1) then
+      call this%riv_store_user_cond()
+    end if
+    !
+    ! -- Write the list to iout if requested
+    if (this%iprpak /= 0) then
+      call this%write_list()
+    end if
+    !
+    ! -- return
+    return
+  end subroutine riv_rp
+
   subroutine riv_ck(this)
 ! ******************************************************************************
 ! riv_ck -- Check river boundary condition data
@@ -147,8 +277,8 @@ subroutine riv_ck(this)
     do i = 1, this%nbound
       node = this%nodelist(i)
       bt = this%dis%bot(node)
-      stage = this%bound(1, i)
-      rbot = this%bound(3, i)
+      stage = this%stage(i)
+      rbot = this%rbot(i)
       ! -- accumulate errors
       if (rbot < bt .and. this%icelltype(node) /= 0) then
         write (errmsg, fmt=fmtriverr) i, rbot, bt
@@ -209,9 +339,9 @@ subroutine riv_cf(this, reset_mover)
         this%rhs(i) = DZERO
         cycle
       end if
-      hriv = this%bound(1, i)
-      criv = this%bound(2, i)
-      rbot = this%bound(3, i)
+      hriv = this%stage(i)
+      criv = this%cond_mult(i)
+      rbot = this%rbot(i)
       if (this%xnew(node) <= rbot) then
         this%rhs(i) = -criv * (hriv - rbot)
         this%hcof(i) = DZERO
@@ -257,9 +387,9 @@ subroutine riv_fc(this, rhs, ia, idxglo, matrix_sln)
       !
       ! -- If mover is active and this river cell is discharging,
       !    store available water (as positive value).
-      stage = this%bound(1, i)
+      stage = this%stage(i)
       if (this%imover == 1 .and. this%xnew(n) > stage) then
-        cond = this%bound(2, i)
+        cond = this%cond_mult(i)
         qriv = cond * (this%xnew(n) - stage)
         call this%pakmvrobj%accumulate_qformvr(i, qriv)
       end if
@@ -348,35 +478,73 @@ subroutine riv_df_obs(this)
     return
   end subroutine riv_df_obs
 
-  ! -- Procedure related to time series
-
-  subroutine riv_rp_ts(this)
-    ! -- Assign tsLink%Text appropriately for
-    !    all time series in use by package.
-    !    In RIV package variables STAGE, COND, and RBOT
-    !    can be controlled by time series.
-    ! -- dummy
-    class(RivType), intent(inout) :: this
-    ! -- 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 = 'STAGE'
-        case (2)
-          tslink%Text = 'COND'
-        case (3)
-          tslink%Text = 'RBOT'
-        end select
-      end if
+  subroutine riv_store_user_cond(this)
+    ! -- modules
+    ! -- dummy variables
+    class(RivType), intent(inout) :: this !< BndExtType object
+    ! -- local variables
+    integer(I4B) :: n
+    !
+    ! -- store backup copy of conductance values
+    do n = 1, this%nbound
+      this%condinput(n) = this%cond_mult(n)
     end do
     !
+    ! -- return
+    return
+  end subroutine riv_store_user_cond
+
+  function cond_mult(this, row) result(cond)
+    ! -- modules
+    use ConstantsModule, only: DZERO
+    ! -- dummy variables
+    class(RivType), intent(inout) :: this !< BndExtType object
+    integer(I4B), intent(in) :: row
+    ! -- result
+    real(DP) :: cond
+    !
+    if (this%iauxmultcol > 0) then
+      cond = this%cond(row) * this%auxvar(this%iauxmultcol, row)
+    else
+      cond = this%cond(row)
+    end if
+    !
+    ! -- return
+    return
+  end function cond_mult
+
+! ******************************************************************************
+! riv_bound_value -- return requested boundary value
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+  function riv_bound_value(this, col, row) result(bndval)
+    ! -- modules
+    use ConstantsModule, only: DZERO
+    ! -- dummy variables
+    class(RivType), intent(inout) :: this !< BndExtType object
+    integer(I4B), intent(in) :: col
+    integer(I4B), intent(in) :: row
+    ! -- result
+    real(DP) :: bndval
+    !
+    select case (col)
+    case (1)
+      bndval = this%stage(row)
+    case (2)
+      bndval = this%cond_mult(row)
+    case (3)
+      bndval = this%rbot(row)
+    case default
+      errmsg = 'Programming error. RIV bound value requested column '&
+               &'outside range of ncolbnd (3).'
+      call store_error(errmsg)
+      call store_error_filename(this%input_fname)
+    end select
+    !
+    ! -- return
     return
-  end subroutine riv_rp_ts
+  end function riv_bound_value
 
 end module rivmodule
diff --git a/src/Model/GroundWaterFlow/gwf3riv8idm.f90 b/src/Model/GroundWaterFlow/gwf3riv8idm.f90
new file mode 100644
index 00000000000..5b2b81ad9de
--- /dev/null
+++ b/src/Model/GroundWaterFlow/gwf3riv8idm.f90
@@ -0,0 +1,471 @@
+! ** Do Not Modify! MODFLOW 6 system generated file. **
+module GwfRivInputModule
+  use ConstantsModule, only: LENVARNAME
+  use InputDefinitionModule, only: InputParamDefinitionType, &
+                                   InputBlockDefinitionType
+  private
+  public gwf_riv_param_definitions
+  public gwf_riv_aggregate_definitions
+  public gwf_riv_block_definitions
+  public GwfRivParamFoundType
+  public gwf_riv_multi_package
+  public gwf_riv_aux_sfac_param
+
+  type GwfRivParamFoundType
+    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 :: mover = .false.
+    logical :: maxbound = .false.
+    logical :: cellid = .false.
+    logical :: stage = .false.
+    logical :: cond = .false.
+    logical :: rbot = .false.
+    logical :: auxvar = .false.
+    logical :: boundname = .false.
+  end type GwfRivParamFoundType
+
+  logical :: gwf_riv_multi_package = .true.
+
+  character(len=LENVARNAME) :: gwf_riv_aux_sfac_param = ''
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_auxiliary = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_auxmultname = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_boundnames = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_iprpak = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_iprflow = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_ipakcb = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_ts_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_ts6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_filein = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_ts6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_obs_filerecord = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_obs6 = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_obs6_filename = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_mover = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! subcomponent
+    'OPTIONS', & ! block
+    'MOVER', & ! tag name
+    'MOVER', & ! fortran variable
+    'KEYWORD', & ! type
+    '', & ! shape
+    .false., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_maxbound = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_cellid = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_stage = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! subcomponent
+    'PERIOD', & ! block
+    'STAGE', & ! tag name
+    'STAGE', & ! fortran variable
+    'DOUBLE', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_cond = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! subcomponent
+    'PERIOD', & ! block
+    'COND', & ! tag name
+    'COND', & ! fortran variable
+    'DOUBLE', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_rbot = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! subcomponent
+    'PERIOD', & ! block
+    'RBOT', & ! tag name
+    'RBOT', & ! fortran variable
+    'DOUBLE', & ! type
+    '', & ! shape
+    .true., & ! required
+    .true., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .true. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_auxvar = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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 :: &
+    gwfriv_boundname = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! 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_riv_param_definitions(*) = &
+    [ &
+    gwfriv_auxiliary, &
+    gwfriv_auxmultname, &
+    gwfriv_boundnames, &
+    gwfriv_iprpak, &
+    gwfriv_iprflow, &
+    gwfriv_ipakcb, &
+    gwfriv_ts_filerecord, &
+    gwfriv_ts6, &
+    gwfriv_filein, &
+    gwfriv_ts6_filename, &
+    gwfriv_obs_filerecord, &
+    gwfriv_obs6, &
+    gwfriv_obs6_filename, &
+    gwfriv_mover, &
+    gwfriv_maxbound, &
+    gwfriv_cellid, &
+    gwfriv_stage, &
+    gwfriv_cond, &
+    gwfriv_rbot, &
+    gwfriv_auxvar, &
+    gwfriv_boundname &
+    ]
+
+  type(InputParamDefinitionType), parameter :: &
+    gwfriv_spd = InputParamDefinitionType &
+    ( &
+    'GWF', & ! component
+    'RIV', & ! subcomponent
+    'PERIOD', & ! block
+    'STRESS_PERIOD_DATA', & ! tag name
+    'SPD', & ! fortran variable
+    'RECARRAY CELLID STAGE COND RBOT AUX BOUNDNAME', & ! type
+    'MAXBOUND', & ! shape
+    .true., & ! required
+    .false., & ! multi-record
+    .false., & ! preserve case
+    .false., & ! layered
+    .false. & ! timeseries
+    )
+
+  type(InputParamDefinitionType), parameter :: &
+    gwf_riv_aggregate_definitions(*) = &
+    [ &
+    gwfriv_spd &
+    ]
+
+  type(InputBlockDefinitionType), parameter :: &
+    gwf_riv_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 GwfRivInputModule
diff --git a/src/Model/GroundWaterFlow/gwf3vsc8.f90 b/src/Model/GroundWaterFlow/gwf3vsc8.f90
index ef8f55d4026..4b164d36f8c 100644
--- a/src/Model/GroundWaterFlow/gwf3vsc8.f90
+++ b/src/Model/GroundWaterFlow/gwf3vsc8.f90
@@ -520,6 +520,7 @@ subroutine vsc_ad_standard_bnd(packobj, hnew, visc, viscref, locelev, &
     ! -- modules
     use BndModule, only: BndType
     use DrnModule, only: DrnType
+    use RivModule, only: RivType
     class(BndType), pointer :: packobj
     ! -- dummy
     real(DP), intent(in), dimension(:) :: hnew
@@ -559,6 +560,12 @@ subroutine vsc_ad_standard_bnd(packobj, hnew, visc, viscref, locelev, &
           packobj%cond(n) = update_bnd_cond(viscbnd, viscref, &
                                             packobj%condinput(n))
         end select
+      case ('RIV')
+        select type (packobj)
+        type is (RivType)
+          packobj%cond(n) = update_bnd_cond(viscbnd, viscref, &
+                                            packobj%condinput(n))
+        end select
       case default
         packobj%bound(2, n) = update_bnd_cond(viscbnd, viscref, &
                                               packobj%condinput(n))
diff --git a/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90 b/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
index c9c15d09b60..a8477edb9ca 100644
--- a/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
+++ b/src/Utilities/Idm/selector/IdmGwfDfnSelector.f90
@@ -11,6 +11,7 @@ module IdmGwfDfnSelectorModule
   use GwfDisvInputModule
   use GwfDrnInputModule
   use GwfNpfInputModule
+  use GwfRivInputModule
   use GwfWelInputModule
   use GwfNamInputModule
 
@@ -116,6 +117,8 @@ module IdmGwfDfnSelectorModule
     logical :: angle2 = .false.
     logical :: angle3 = .false.
     logical :: wetdry = .false.
+    logical :: stage = .false.
+    logical :: rbot = .false.
     logical :: flowred = .false.
     logical :: afrcsv_rec = .false.
     logical :: afrcsv = .false.
@@ -165,6 +168,8 @@ function gwf_param_definitions(subcomponent) result(input_definition)
       call set_param_pointer(input_definition, gwf_drn_param_definitions)
     case ('NPF')
       call set_param_pointer(input_definition, gwf_npf_param_definitions)
+    case ('RIV')
+      call set_param_pointer(input_definition, gwf_riv_param_definitions)
     case ('WEL')
       call set_param_pointer(input_definition, gwf_wel_param_definitions)
     case ('NAM')
@@ -191,6 +196,8 @@ function gwf_aggregate_definitions(subcomponent) result(input_definition)
       call set_param_pointer(input_definition, gwf_drn_aggregate_definitions)
     case ('NPF')
       call set_param_pointer(input_definition, gwf_npf_aggregate_definitions)
+    case ('RIV')
+      call set_param_pointer(input_definition, gwf_riv_aggregate_definitions)
     case ('WEL')
       call set_param_pointer(input_definition, gwf_wel_aggregate_definitions)
     case ('NAM')
@@ -217,6 +224,8 @@ function gwf_block_definitions(subcomponent) result(input_definition)
       call set_block_pointer(input_definition, gwf_drn_block_definitions)
     case ('NPF')
       call set_block_pointer(input_definition, gwf_npf_block_definitions)
+    case ('RIV')
+      call set_block_pointer(input_definition, gwf_riv_block_definitions)
     case ('WEL')
       call set_block_pointer(input_definition, gwf_wel_block_definitions)
     case ('NAM')
@@ -242,6 +251,8 @@ function gwf_idm_multi_package(subcomponent) result(multi_package)
       multi_package = gwf_drn_multi_package
     case ('NPF')
       multi_package = gwf_npf_multi_package
+    case ('RIV')
+      multi_package = gwf_riv_multi_package
     case ('WEL')
       multi_package = gwf_wel_multi_package
     case ('NAM')
@@ -270,6 +281,8 @@ function gwf_idm_sfac_param(subcomponent) result(sfac_param)
       sfac_param = gwf_drn_aux_sfac_param
     case ('NPF')
       sfac_param = gwf_npf_aux_sfac_param
+    case ('RIV')
+      sfac_param = gwf_riv_aux_sfac_param
     case ('WEL')
       sfac_param = gwf_wel_aux_sfac_param
     case ('NAM')
@@ -299,6 +312,8 @@ function gwf_idm_integrated(subcomponent) result(integrated)
       integrated = .true.
     case ('NPF')
       integrated = .true.
+    case ('RIV')
+      integrated = .true.
     case ('WEL')
       integrated = .true.
     case ('NAM')
diff --git a/src/meson.build b/src/meson.build
index 83fc181b801..79f2da7f3d6 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -74,6 +74,7 @@ modflow_sources = files(
     'Model' / 'GroundWaterFlow' / 'gwf3oc8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3rch8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3riv8.f90',
+    'Model' / 'GroundWaterFlow' / 'gwf3riv8idm.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3sfr8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3sto8.f90',
     'Model' / 'GroundWaterFlow' / 'gwf3tvbase8.f90',
diff --git a/utils/idmloader/scripts/dfn2f90.py b/utils/idmloader/scripts/dfn2f90.py
index 1b8d67f3264..99e9c350619 100644
--- a/utils/idmloader/scripts/dfn2f90.py
+++ b/utils/idmloader/scripts/dfn2f90.py
@@ -992,6 +992,10 @@ 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-riv.dfn"),
+            Path("../../../src/Model/GroundWaterFlow", "gwf3riv8idm.f90"),
+        ],
         [
             Path("../../../doc/mf6io/mf6ivar/dfn", "gwf-wel.dfn"),
             Path("../../../src/Model/GroundWaterFlow", "gwf3wel8idm.f90"),