diff --git a/src/earthkit/data/utils/dates.py b/src/earthkit/data/utils/dates.py
index baa8cecf..e3dcfe4b 100644
--- a/src/earthkit/data/utils/dates.py
+++ b/src/earthkit/data/utils/dates.py
@@ -67,6 +67,22 @@ def to_datetime(dt):
     if hasattr(dt, "isoformat"):
         return datetime.datetime.fromisoformat(dt.isoformat())
 
+    # if dt is NaT-like, then return None
+    if dt is None:
+        return dt
+    try:
+        # covers nan, np.datetime64('NaT') and np.timedelta64('NaT')
+        if np.isnan(dt):
+            return None
+    except Exception:
+        pass
+    try:
+        # covers pd.NaT
+        if np.isnan(dt.second):
+            return None
+    except Exception:
+        pass
+
     dt = from_object(dt)
 
     return to_datetime(dt.to_datetime())
diff --git a/tests/xr_engine/test_xr_engine_add_valid_time_coord.py b/tests/xr_engine/test_xr_engine_add_valid_time_coord.py
index 5438e12c..4c551a1c 100644
--- a/tests/xr_engine/test_xr_engine_add_valid_time_coord.py
+++ b/tests/xr_engine/test_xr_engine_add_valid_time_coord.py
@@ -243,6 +243,34 @@ def test_fixed_dims_different_order(pl_fl, lazy_load, allow_holes):
     assert ds.coords["valid_time"].shape == (2, 4)
 
 
+# -------------------------------------------------------------------------
+# valid_time aux coord should work even if there is no field for a given
+# (forecast_reference_time, step) pair (i.e. holes in the data) and thus
+# the corresponding valid_time value is NaT
+# -------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("lazy_load", [True, False])
+@pytest.mark.parametrize("dim_name_from_role_name", [True, False])
+def test_NaT(pl_fl, lazy_load, dim_name_from_role_name):
+    fl_with_holes = pl_fl.to_fieldlist().sel({"vertical.level": 500})[2:]
+    ds = fl_with_holes.to_xarray(
+        profile="earthkit",
+        add_valid_time_coord=True,
+        dim_name_from_role_name=dim_name_from_role_name,
+        lazy_load=lazy_load,
+        allow_holes=True,
+    )
+
+    assert "valid_time" in ds.coords
+    assert "valid_time" not in ds.sizes
+    assert ds.coords["valid_time"].dims == ("forecast_reference_time", "step")
+    assert ds.coords["valid_time"].shape == (4, 2)
+    _valid_time = np.array(VALID_TIME_FRT_STEP)
+    _valid_time[0, 0] = np.datetime64("NaT")
+    np.testing.assert_array_equal(ds.coords["valid_time"].values, _valid_time)
+
+
 # -------------------------------------------------------------------------
 # Edge case: add_valid_time_coord=False should not add it
 # -------------------------------------------------------------------------