diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst
index de7d1156c..ccc309c23 100644
--- a/docs/api/changelog.rst
+++ b/docs/api/changelog.rst
@@ -58,6 +58,9 @@ Fixed
 - Fixed :func:`imod.prepare.spatial.polygonize` for polygons with holes.
 - :func:`imod.formats.prj.open_projectfile_data` now drops empty wells from the
   dataset, and logs a warning about it.
+- :meth:`imod.mf6.NodePropertyFlow.regrid_like` now regrids ``k33`` using the
+  correct method, namely ``mean`` instead of ``harmonic_mean``. As this is the
+  appropriate method for horizontal regridding of ``k33``.
 
 Changed
 ~~~~~~~
diff --git a/examples/user-guide/08-regridding.py b/examples/user-guide/08-regridding.py
index 9d932c541..15d1e33ef 100644
--- a/examples/user-guide/08-regridding.py
+++ b/examples/user-guide/08-regridding.py
@@ -372,22 +372,25 @@ def plot_histograms_side_by_side(array_original, array_regridded, title):
 # A note on regridding conductivity
 # ---------------------------------
 #
-# In the npf package, it is possible to use for definining the conductivity
-# tensor:
-#
-# - 1 array (K)
-# - 2 arrays (K and K22)
-# - 3 arrays (K, K22, K33)
-#
-# If 1 array is given the tensor is called isotropic. Defining only K gives the
-# same behavior as specifying K, K22 and K33 with the same value. When
-# regridding, K33 has a default method different from that of K and K22, but it
-# can only be applied if K33 exists in the source model in the first place. So
-# it is recommended to introduce K33 as a separate array in the source model
-# even if it is isotropic. Also note that default regridding methods were chosen
-# assuming that K and K22 are roughly horizontal and K33 roughly vertical. But
-# this may not be the case if the input arrays angle2 and/or angle3 have large
-# values.
+# By default, K and K22 are regrid with a different method than K33. Namely K
+# and K22 are regridded by default with a geometric mean (Bierkens, 1998),
+# whereas K33 is regrid with an arithmetic mean. Note that default regridding
+# methods were chosen assuming that K and K22 are roughly horizontal and K33
+# roughly vertical. But this may not be the case if the input arrays angle2
+# and/or angle3 have large values.
+#
+# Furthermore it is possible to only provide one array for the hydraulic
+# conductivity, K. In that case, the same k values are used in all directions
+# (K, K22, and K33), and the model is called isotropic. It is recommended to
+# introduce K33 as a separate array in the source model even if it is isotropic
+# when regridding simulations, as K33 is regridded with a different method.
+#
+# Reference:
+#
+# Bierkens, M., van der Gaast, J. Upscaling hydraulic conductivity: theory and
+# examples from geohydrological studies. *Nutrient Cycling in Agroecosystems* **50**,
+# 193–207 (1998). https://doi.org/10.1023/A:1009740328153
+#
 #
 # Regridding boundary conditions
 # ------------------------------
diff --git a/imod/mf6/regrid/regrid_schemes.py b/imod/mf6/regrid/regrid_schemes.py
index 496db71b8..0af7267af 100644
--- a/imod/mf6/regrid/regrid_schemes.py
+++ b/imod/mf6/regrid/regrid_schemes.py
@@ -289,7 +289,7 @@ class NodePropertyFlowRegridMethod(DataclassType):
     icelltype: tuple, defaults (RegridderType.OVERLAP, "mean")
     k: tuple, defaults ( RegridderType.OVERLAP,"geometric_mean")
     k22: tuple, defaults (RegridderType.OVERLAP,"geometric_mean")
-    k33: tuple, defaults (RegridderType.OVERLAP,"harmonic_mean")
+    k33: tuple, defaults (RegridderType.OVERLAP,"mean")
     angle1: tuple, defaults (RegridderType.OVERLAP, "mean")
     angle2: tuple, defaults (RegridderType.OVERLAP, "mean")
     angle3: tuple, defaults (RegridderType.OVERLAP, "mean")
@@ -319,7 +319,7 @@ class NodePropertyFlowRegridMethod(DataclassType):
     )  # horizontal if angle2 = 0 & angle3 = 0
     k33: RegridVarType = (
         RegridderType.OVERLAP,
-        "harmonic_mean",
+        "mean",
     )  # vertical if angle2 = 0 & angle3 = 0
     angle1: RegridVarType = (RegridderType.OVERLAP, "mean")
     angle2: RegridVarType = (RegridderType.OVERLAP, "mean")
diff --git a/pixi.lock b/pixi.lock
index 39e7316da..03d4b3985 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -14752,7 +14752,7 @@ packages:
 - pypi: ./
   name: imod
   version: 1.0.1.dev1
-  sha256: b242a1acd646ffa5410ebc6470e9fa558e81e42073304819f54be9e77949a462
+  sha256: d27105972c53a243e6793668ecb915b0b24268a493441e035e214e01a8a8e389
   requires_dist:
   - affine
   - cftime>=1