diff --git a/docs/user_guide/examples/tutorial_sampling.ipynb b/docs/user_guide/examples/tutorial_sampling.ipynb
index 13ec95396..53bbc82a5 100644
--- a/docs/user_guide/examples/tutorial_sampling.ipynb
+++ b/docs/user_guide/examples/tutorial_sampling.ipynb
@@ -272,6 +272,15 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "```{note}\n",
+    "Note that you need to explicitly set the `t` values in the ParticleSet; otherwise you will get a `ValueError: Time values cannot be NaN.` error. This is because (if not specified) the time values are NaN until they get set on `pset.execute()`, when the sign of `dt` is known.\n",
+    "```"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py
index 0cff9adb2..888cdcca3 100644
--- a/src/parcels/_core/field.py
+++ b/src/parcels/_core/field.py
@@ -393,6 +393,9 @@ def _assert_same_time_interval(fields: Sequence[Field]) -> None:
 
 def _get_positions(field: Field, t, z, y, x, particles, _ei) -> tuple[dict, dict]:
     """Initialize and populate particle_positions and grid_positions dictionaries"""
+    if np.any(np.isnan(t)):
+        nan_indices = np.where(np.isnan(t))[0]
+        raise ValueError(f"Time values for particles with indices {nan_indices} cannot be NaN.")
     particle_positions = {"t": t, "z": z, "y": y, "x": x}
     grid_positions = {}
     grid_positions.update(_search_time_index(field, t))
diff --git a/tests/test_particleset.py b/tests/test_particleset.py
index 4e3ce9dc5..d0721b7ed 100644
--- a/tests/test_particleset.py
+++ b/tests/test_particleset.py
@@ -211,3 +211,22 @@ def test_pset_default_z_closest_to_zero(depths):
     pset = ParticleSet(fieldset, x=[0], y=[0])
     expected_z = depths[np.argmin(np.abs(depths))]
     assert np.isclose(pset.z[0], expected_z)
+
+
+@pytest.mark.parametrize("npart", [1, 10])
+@pytest.mark.parametrize("witht", [True, False])
+def test_sampling_pset(fieldset, npart, witht):
+    # Test that inital value of a field gets sampled
+    fieldset.U.data[:] = 2.0
+
+    x = np.zeros(npart)
+    y = np.zeros(npart)
+    if witht:
+        t = npart * [np.timedelta64(0, "s")]
+        pset = ParticleSet(fieldset, x=x, y=y, t=t)
+        pset.sample, _ = fieldset.UV[pset]
+        np.testing.assert_allclose(pset.sample, 2.0, rtol=1e-12)
+    else:
+        with pytest.raises(ValueError, match="Time values for particles with indices .* cannot be NaN."):
+            pset = ParticleSet(fieldset, x=x, y=y)
+            pset.sample, _ = fieldset.UV[pset]