Solution you refer to might work on very dense stream network. Obviously it is not your case.
I assume that what you call HAND is height above elevation of nearest drainage points downstream, i.e. first point downhill.
Use raster calculator to define your ‘streams’ from flow accumulation raster and assign elevation to it from Filled raster, e.g.:
Con("facc"=>350000,"filled")
Call output raster “flood”.
Replace 3 lines in below script pointing to flow direction grid, “flood” and output grid name (at the bottom of the script)
Script assigns elevation of drainage point to cells within point catchment, and saves result to output raster. In order to obtain HAND you’ll need to subtract output raster from filled DEM.

import arcpy, os, traceback, sys, numpy
from arcpy import env
from arcpy.sa import *
fdir=r'D:\Scratch\fdir'
demG=r'D:\Scratch\flood'
d=arcpy.Describe(demG)
##SR=d.spatialReference
try:
def showPyMessage():
arcpy.AddMessage(str(time.ctime()) + " - " + message)
dirArray = arcpy.RasterToNumPyArray(fdir,"","","",-9999)
dem = arcpy.RasterToNumPyArray(demG,"","","",-9999)
nRows,nCols=dirArray.shape
blankArray=arcpy.RasterToNumPyArray(demG,"","","",-9999)
cellsTotal=nCols*nRows
d=arcpy.Describe(fdir)
origin=d.extent.lowerLeft
cSize=arcpy.Raster(fdir).meanCellHeight
## directions to find cell neighbour
fDirs=(1,2,4,8,16,32,64,128)
dCol=(1, 1, 0, -1, -1,-1, 0,1)
dRow=(0, -1, -1, -1, 0, 1, 1,1)
## flipped
dRow=(0, 1, 1, 1, 0, -1, -1,-1)
## main loop
arcpy.SetProgressor("step", "", 0, nRows)
for nRow in range (nRows):
for nCol in range (nCols):
S=dirArray[nRow,nCol]
if S in (-1,-9999): continue
nR,nC=nRow,nCol
cells=[(nR,nC)]
z=dem[nR,nC]
while True:
direction=dirArray[nR,nC]
i=fDirs.index(direction)
dX=dCol[i];nC+=dX
## if nC in (0,nCols): break
if nC not in range(nCols): break
dY=dRow[i];nR+=dY
## if nR in (0,nRows): break
if nR not in range(nRows): break
cells.append((nR,nC))
S=dirArray[nR,nC]
z=dem[nR,nC]
if S<-9998 or z>-9998:
break
if S==-1:
z=blankArray[nR,nC]
break
for nR,nC in cells:
blankArray[nR,nC]=z
dirArray[nR,nC]=-1
arcpy.SetProgressorPosition()
myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
oneGrid=Con(myRaster<>-9999,myRaster)
oneGrid.save(r'D:\Aerials\nztm\toriver')
del dirArray,blankArray
except:
message = "\n*** PYTHON ERRORS *** "; showPyMessage()
message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()
Note: It takes about 5 minutes to process 3000*3000 cells inputs.