diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90
index e0441cac2e..e686261fdf 100644
--- a/src/diagnostics/MOM_obsolete_params.F90
+++ b/src/diagnostics/MOM_obsolete_params.F90
@@ -71,6 +71,7 @@ subroutine find_obsolete_params(param_file)
 
   call obsolete_real(param_file, "VSTAR_SCALE_COEF")
   call obsolete_real(param_file, "ZSTAR_RIGID_SURFACE_THRESHOLD")
+  call obsolete_logical(param_file, "HENYEY_IGW_BACKGROUND_NEW")
 
   ! Test for inconsistent parameter settings.
   split = .true. ; test_logic = .false.
diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90
index 9eabd8dc4d..179dfa470c 100644
--- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90
+++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90
@@ -280,7 +280,7 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL
   call get_param(param_file, mdl, "HENYEY_IGW_BACKGROUND_NEW", CS%Henyey_IGW_background_new, &
                  "If true, use a better latitude-dependent scaling for the "//&
                  "background diffusivity, as described in "//&
-                 "Harrison & Hallberg, JPO 2008.", default=.false.)
+                 "Harrison & Hallberg, JPO 2008. This option is obsolete.", default=.false.)
   if (CS%Henyey_IGW_background_new) call check_bkgnd_scheme(CS, "HENYEY_IGW_BACKGROUND_NEW")
 
   if (CS%Kd>0.0 .and. (trim(CS%bkgnd_scheme_str)=="BRYAN_LEWIS_DIFFUSIVITY" .or.&
@@ -316,7 +316,7 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL
 
   CS%Kd_via_Kdml_bug = .false.
   if ((CS%Kd /= CS%Kd_tot_ml) .and. .not.(CS%Kd_tanh_lat_fn .or. CS%physical_OBL_scheme .or. &
-                                     CS%Henyey_IGW_background .or. CS%Henyey_IGW_background_new .or. &
+                                     CS%Henyey_IGW_background .or. &
                                      CS%horiz_varying_background .or. CS%Bryan_Lewis_diffusivity)) then
     call get_param(param_file, mdl, "KD_BACKGROUND_VIA_KDML_BUG", CS%Kd_via_Kdml_bug, &
                  "If true and KDML /= KD and several other conditions apply, the background "//&
@@ -453,25 +453,6 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G,
       Kd_lay(i,k) = Kd_int(i,1)
     enddo ; enddo
 
-  elseif (CS%Henyey_IGW_background_new) then
-    I_x30 = 2.0 / invcosh(CS%N0_2Omega*2.0) ! This is evaluated at 30 deg.
-    I_2Omega = 0.5 / CS%omega
-    do k=1,nz ; do i=is,ie
-      abs_sinlat = max(min_sinlat, abs(sin(G%geoLatT(i,j)*deg_to_rad)))
-      N_2Omega = max(abs_sinlat, sqrt(N2_lay(i,k))*I_2Omega)
-      N02_N2 = (CS%N0_2Omega/N_2Omega)**2
-      Kd_lay(i,k) = max(CS%Kd_min, CS%Kd * &
-           ((abs_sinlat * invcosh(N_2Omega/abs_sinlat)) * I_x30)*N02_N2)
-    enddo ; enddo
-    ! Update Kd_int and Kv_bkgnd, based on Kd_lay.  These might be just used for diagnostic purposes.
-    do i=is,ie
-      Kd_int(i,1) = 0.0; Kv_bkgnd(i,1) = 0.0
-      Kd_int(i,nz+1) = 0.0; Kv_bkgnd(i,nz+1) = 0.0
-    enddo
-    do K=2,nz ; do i=is,ie
-      Kd_int(i,K) = 0.5*(Kd_lay(i,k-1) + Kd_lay(i,k))
-      Kv_bkgnd(i,K) = Kd_int(i,K) * CS%prandtl_bkgnd
-    enddo ; enddo
   else
     ! Set a potentially spatially varying surface value of diffusivity.
     if (CS%Henyey_IGW_background) then
diff --git a/src/parameterizations/vertical/_V_diffusivity.dox b/src/parameterizations/vertical/_V_diffusivity.dox
index f3b7ed5962..df1ce50e27 100644
--- a/src/parameterizations/vertical/_V_diffusivity.dox
+++ b/src/parameterizations/vertical/_V_diffusivity.dox
@@ -253,26 +253,6 @@ in \cite harrison2008, but that isn't what is in the MOM6 code. Instead, the sur
 value is propagated down, with the assumption that the tidal mixing parameterization
 will provide the deep mixing: \ref section_Internal_Tidal_Mixing.
 
-There is also a "new" Henyey version, taking into account the effect of stratification on
-TKE dissipation,
-
-\todo Harrison (personal communication) recommends that this option be made obsolete and
-eventually removed.
-
-\f[
-   \epsilon = \epsilon_0 \frac{f}{f_0} \frac{\mbox{acosh} (N/f)}{\mbox{acosh} (N_0 / f_0)}
-\f]
-
-where \f$N_0\f$ and \f$f_0\f$ are the reference buoyancy frequency and inertial frequencies, respectively
-and \f$\epsilon_0\f$ is the reference dissipation at \f$(N_0, f_0)\f$. In the previous version, \f$N =
-N_0\f$. Additionally, the relationship between diapycnal diffusivities and stratification is included:
-
-\f[
-   \kappa = \frac{\epsilon}{N^2}
-\f]
-This approach assumes that work done against gravity is uniformly distributed throughout the water column.
-The original version concentrates buoyancy work in regions of strong stratification.
-
 \subsection subsection_danabasoglu_back Danabasoglu background mixing
 
 The shape of the \cite danabasoglu2012 background mixing has a uniform background value, with a dip