The grid artifacts appear to be the meshblocks for mpi parallelization when running the simulation. This doesn’t occur for any other field.
I am working with an Athena++ athdf dataset and using yt version 3.7.dev0.
Here is the script I used to add the desired field:

import yt
from yt.fields.derived_field import ValidateSpatial
import yt.units as u

kappa0 = 9.68e4

ds = yt.load('../pleiades/grf/grf_pr_2_3_nm5_3/TI.uov.00200.athdf',\
             units_override=units_override,unit_system="cgs")

grad_fields = ds.add_gradient_fields(("gas","temperature"))

def _therm_cond(field,data):
    kappa = kappa0*u.erg/(u.cm*u.K*u.s)
    sl_left = slice(None, -2, None)
    sl_right = slice(2, None, None)
    div_fac = 2.0
    ds = div_fac * data['dx'].flat[0]
    f  = data["temperature_gradient_x"][sl_right,1:-1,1:-1]/ds
    f -= data["temperature_gradient_x"][sl_left ,1:-1,1:-1]/ds
    if data.ds.dimensionality > 1:
        ds = div_fac * data['dy'].flat[0]
        f += data["temperature_gradient_y"][1:-1,sl_right,1:-1]/ds
        f -= data["temperature_gradient_y"][1:-1,sl_left ,1:-1]/ds
    if data.ds.dimensionality > 2:
        ds = div_fac * data['dz'].flat[0]
        f += data["temperature_gradient_z"][1:-1,1:-1,sl_right]/ds
        f -= data["temperature_gradient_z"][1:-1,1:-1,sl_left ]/ds
    new_field = np.zeros(data["temperature_gradient_x"].shape, dtype='float64')
    new_field = data.ds.arr(new_field,'K/cm**2')
    new_field[1:-1,1:-1,1:-1] = f
    return kappa*new_field/data['density']

ds.add_field(('gas','therm_cond'), function=_therm_cond, units="erg/(g*s)", sampling_type="cell", take_log=False,
             display_name=r'$\nabla \cdot (\kappa\nabla T)/\rho$', force_override=True, validators = [ValidateSpatial(1,
                           ['temperature_gradient_x','temperature_gradient_y','temperature_gradient_z'])])

Has anyone come across this issue before? Is this a known bug or am I just implementing the divergence field incorrectly? Any help is greatly appreciated.

Thanks,
Michael