diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90
index 15f1116190..10fe37da89 100644
--- a/src/parameterizations/vertical/MOM_geothermal.F90
+++ b/src/parameterizations/vertical/MOM_geothermal.F90
@@ -35,6 +35,10 @@ module MOM_geothermal
   type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                              !! regulate the timing of diagnostic output.
+  integer :: id_internal_heat_heat_tendency = -1   !< ID for diagnostic of heat tendency
+  integer :: id_internal_heat_temp_tendency = -1   !< ID for diagnostic of temperature tendency
+  integer :: id_internal_heat_h_tendency = -1      !< ID for diagnostic of thickness tendency
+
 end type geothermal_CS
 
 contains
@@ -100,6 +104,17 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
   real :: Irho_cp       ! inverse of heat capacity per unit layer volume
                         ! [degC H m2 J-1 ~> degC m3 J-1 or degC kg J-1]
 
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_old   ! Temperature of each layer
+                                                      ! before any heat is added,
+                                                      ! for diagnostics [degC]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_old   ! Thickness of each layer
+                                                      ! before any heat is added,
+                                                      ! for diagnostics [m or kg m-2]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d ! Scratch variable used to
+                                                      ! calculate change in heat
+                                                      ! due to geothermal
+  real :: Idt           ! inverse of the timestep [s-1]
+
   logical :: do_i(SZI_(G))
   integer :: i, j, k, is, ie, js, je, nz, k2, i2
   integer :: isj, iej, num_start, num_left, nkmb, k_tgt
@@ -119,6 +134,7 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
   Angstrom  = GV%Angstrom_H
   H_neglect = GV%H_subroundoff
   p_ref(:)  = tv%P_Ref
+  Idt       = 1.0 / dt
 
   if (.not.associated(tv%T)) call MOM_error(FATAL, "MOM geothermal: "//&
       "Geothermal heating can only be applied if T & S are state variables.")
@@ -175,6 +191,18 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
     do k=nz,1,-1
       do i=isj,iej ; if (do_i(i)) then
 
+        ! Save temperature and thickness before any changes are made (for diagnostic)
+        if (CS%id_internal_heat_h_tendency > 0 &
+             .or. CS%id_internal_heat_heat_tendency > 0 &
+             .or. CS%id_internal_heat_temp_tendency > 0 ) then
+          h_old(i,j,k) = h(i,j,k)
+        endif
+        if (CS%id_internal_heat_heat_tendency > 0 &
+             .or. CS%id_internal_heat_temp_tendency > 0) then
+          T_old(i,j,k) = tv%T(i,j,k)
+        endif
+
+
         if (h(i,j,k) > Angstrom) then
           if ((h(i,j,k)-Angstrom) >= h_geo_rem(i)) then
             h_heated = h_geo_rem(i)
@@ -294,6 +322,12 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
             ! endif
           endif
         endif
+
+        ! Calculate heat tendency due to addition and transfer of internal heat
+        if (CS%id_internal_heat_heat_tendency > 0) then
+          work_3d(i,j,k) = ((GV%H_to_kg_m2 * tv%C_p) * Idt) * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * T_old(i,j,k))
+        endif
+
       endif ; enddo
       if (num_left <= 0) exit
     enddo ! k-loop
@@ -304,6 +338,23 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo)
     enddo ; endif
   enddo ! j-loop
 
+  ! Post diagnostic of 3D tendencies (heat, temperature, and thickness) due to internal heat
+  if (CS%id_internal_heat_heat_tendency > 0) then
+    call post_data(CS%id_internal_heat_heat_tendency, work_3d, CS%diag, alt_h = h_old)
+  endif
+  if (CS%id_internal_heat_temp_tendency > 0) then
+    do j=js,je; do i=is,ie; do k=nz,1,-1
+      work_3d(i,j,k) = Idt * (tv%T(i,j,k) - T_old(i,j,k))
+    enddo; enddo; enddo
+    call post_data(CS%id_internal_heat_temp_tendency, work_3d, CS%diag, alt_h = h_old)
+  endif
+  if (CS%id_internal_heat_h_tendency > 0) then
+    do j=js,je; do i=is,ie; do k=nz,1,-1
+      work_3d(i,j,k) = Idt * (h(i,j,k) - h_old(i,j,k))
+    enddo; enddo; enddo
+    call post_data(CS%id_internal_heat_h_tendency, work_3d, CS%diag, alt_h = h_old)
+  endif
+
 !  do i=is,ie ; do j=js,je
 !    resid(i,j) = tv%internal_heat(i,j) - resid(i,j) - GV%H_to_kg_m2 * &
 !           (G%mask2dT(i,j) * (CS%geo_heat(i,j) * (dt*Irho_cp)))
@@ -392,6 +443,20 @@ subroutine geothermal_init(Time, G, param_file, diag, CS)
         x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
   if (id > 0) call post_data(id, CS%geo_heat, diag, .true.)
 
+  ! Diagnostic for tendencies due to internal heat (in 3d)
+  CS%id_internal_heat_heat_tendency=register_diag_field('ocean_model', &
+        'internal_heat_heat_tendency', diag%axesTL, Time,              &
+        'Heat tendency (in 3D) due to internal (geothermal) sources',  &
+        'W m-2', v_extensive = .true.)
+  CS%id_internal_heat_temp_tendency=register_diag_field('ocean_model', &
+        'internal_heat_temp_tendency', diag%axesTL, Time,              &
+        'Temperature tendency (in 3D) due to internal (geothermal) sources', &
+        'degC s-1', v_extensive = .true.)
+  CS%id_internal_heat_h_tendency=register_diag_field('ocean_model',    &
+        'internal_heat_h_tendency', diag%axesTL, Time,                &
+        'Thickness tendency (in 3D) due to internal (geothermal) sources', &
+        'm OR kg m-2', v_extensive = .true.)
+
 end subroutine geothermal_init
 
 !> Clean up and deallocate memory associated with the geothermal heating module.