diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md
index a726ac597..338efbec1 100644
--- a/docs/user_guide/examples/explanation_kernelloop.md
+++ b/docs/user_guide/examples/explanation_kernelloop.md
@@ -103,11 +103,10 @@ First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]`
 ```{code-cell}
 :tags: [hide-output]
 npart = 10
-z = np.repeat(ds_fields.depth[0].values, npart)
 lons = np.repeat(32.2, npart)
 lats = np.linspace(-32.5, -30.5, npart)
 
-pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons)
+pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, y=lats, x=lons)
 output_file = parcels.ParticleFile(
     path="advection_then_wind.parquet", outputdt=np.timedelta64(6,'h')
 )
@@ -124,7 +123,7 @@ Then also run a simulation where we apply the Kernels in the reverse order as `[
 ```{code-cell}
 :tags: [hide-output]
 pset_reverse = parcels.ParticleSet(
-    fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons
+    fieldset, pclass=parcels.Particle, y=lats, x=lons
 )
 output_file_reverse = parcels.ParticleFile(
     path="wind_then_advection.parquet", outputdt=np.timedelta64(6,"h")
diff --git a/docs/user_guide/examples/tutorial_croco_3D.ipynb b/docs/user_guide/examples/tutorial_croco_3D.ipynb
index 3458e3df1..a393739f1 100644
--- a/docs/user_guide/examples/tutorial_croco_3D.ipynb
+++ b/docs/user_guide/examples/tutorial_croco_3D.ipynb
@@ -261,7 +261,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "When using Croco data, Parcels needs to convert the particles.depth to sigma-coordinates, before doing any interpolation. This is done with the `convert_z_to_sigma_croco()`  function. Interpolating onto a Field is then done like:\n",
+    "When using Croco data, Parcels needs to convert the particles.z to sigma-coordinates, before doing any interpolation. This is done with the `convert_z_to_sigma_croco()`  function. Interpolating onto a Field is then done like:\n",
     "\n",
     "```python\n",
     "def SampleTempCroco(particles, fieldset):\n",
diff --git a/docs/user_guide/examples/tutorial_delaystart.ipynb b/docs/user_guide/examples/tutorial_delaystart.ipynb
index d8dfa2ce3..3b57bb742 100644
--- a/docs/user_guide/examples/tutorial_delaystart.ipynb
+++ b/docs/user_guide/examples/tutorial_delaystart.ipynb
@@ -93,10 +93,9 @@
     "lat = np.linspace(-31.5, -30.5, npart, dtype=np.float32)\n",
     "# release every particle one hour later from the initial fieldset time\n",
     "time = ds_fields.time.values[0] + np.arange(0, npart) * np.timedelta64(1, \"h\")\n",
-    "z = np.repeat(ds_fields.depth.values[0], npart)\n",
     "\n",
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time, z=z\n",
+    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time\n",
     ")"
    ]
   },
@@ -209,7 +208,6 @@
     "lon_i = 32 * np.ones(npart)\n",
     "lat_i = np.linspace(-31.5, -30.5, npart, dtype=np.float32)\n",
     "time_i = np.repeat(ds_fields.time.values[0], npart)\n",
-    "z_i = np.repeat(ds_fields.depth.values[0], npart)\n",
     "\n",
     "# repeat release at frequency `repeatdt` for `nrepeat` different releases\n",
     "repeatdt = np.timedelta64(6, \"h\")\n",
@@ -222,10 +220,9 @@
     "    + np.arange(0, nrepeat)[:, np.newaxis] * repeatdt\n",
     ")\n",
     "print(f\"pset.time shape = {time.shape}\")\n",
-    "z = np.broadcast_to(z_i, (nrepeat, npart))\n",
     "\n",
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time, z=z\n",
+    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time\n",
     ")"
    ]
   },
@@ -400,7 +397,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Parcels:default (3.14.6)",
    "language": "python",
    "name": "python3"
   },
diff --git a/docs/user_guide/examples/tutorial_diffusion.ipynb b/docs/user_guide/examples/tutorial_diffusion.ipynb
index b0314fbfb..d349dcd49 100644
--- a/docs/user_guide/examples/tutorial_diffusion.ipynb
+++ b/docs/user_guide/examples/tutorial_diffusion.ipynb
@@ -554,7 +554,7 @@
     "time = np.repeat(fieldset.time_interval.left, 13)\n",
     "lon = [32] * len(time)\n",
     "lat = [-30.1] * len(time)\n",
-    "z = [ds_fields.depth[0]] * len(time)\n",
+    "z = [0] * len(time)\n",
     "pset = parcels.ParticleSet(\n",
     "    fieldset=fieldset,\n",
     "    pclass=parcels.Particle,\n",
diff --git a/docs/user_guide/examples/tutorial_dt_integrators.ipynb b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
index 9ca44604d..8dba6b1f9 100644
--- a/docs/user_guide/examples/tutorial_dt_integrators.ipynb
+++ b/docs/user_guide/examples/tutorial_dt_integrators.ipynb
@@ -198,8 +198,7 @@
     "initial_release_lats = np.linspace(-32.5, -32, npart, dtype=np.float32)\n",
     "initial_release_times = np.repeat(\n",
     "    ds_fields.time.values[0], npart\n",
-    ")  # release all particles at the start time of the fieldset\n",
-    "initial_release_zs = np.repeat(ds_fields.depth.values[0], npart)"
+    ")  # release all particles at the start time of the fieldset"
    ]
   },
   {
@@ -239,7 +238,6 @@
     "        fieldset=fieldset,\n",
     "        pclass=parcels.Particle,\n",
     "        t=initial_release_times,\n",
-    "        z=initial_release_zs,\n",
     "        y=initial_release_lats,\n",
     "        x=initial_release_lons,\n",
     "    )\n",
@@ -580,7 +578,6 @@
     "            fieldset=fieldset,\n",
     "            pclass=parcels.Particle,\n",
     "            t=initial_release_times,\n",
-    "            z=initial_release_zs,\n",
     "            y=initial_release_lats,\n",
     "            x=initial_release_lons,\n",
     "        )\n",
@@ -812,7 +809,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Parcels:docs (3.14.6)",
+   "display_name": "Parcels:test (3.14.6)",
    "language": "python",
    "name": "python3"
   },
diff --git a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb
index e536de38b..c023e68e8 100644
--- a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb
+++ b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb
@@ -125,11 +125,10 @@
     "fieldset = fieldset.to_windowed_arrays()\n",
     "\n",
     "npart = 10\n",
-    "z = np.repeat(ds_fields.depth[0].values, npart)\n",
     "lons = np.repeat(32.2, npart)\n",
     "lats = np.linspace(-32.5, -30.5, npart)\n",
     "\n",
-    "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons)\n",
+    "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, y=lats, x=lons)\n",
     "output_file = parcels.ParticleFile(\n",
     "    path=\"summed_advection_wind.parquet\",\n",
     "    outputdt=np.timedelta64(6, \"h\"),\n",
@@ -168,7 +167,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Parcels:docs (3.14.6)",
+   "display_name": "Parcels:test (3.14.6)",
    "language": "python",
    "name": "python3"
   },
diff --git a/docs/user_guide/examples/tutorial_sampling.ipynb b/docs/user_guide/examples/tutorial_sampling.ipynb
index f81e5ec31..a2db9f2ad 100644
--- a/docs/user_guide/examples/tutorial_sampling.ipynb
+++ b/docs/user_guide/examples/tutorial_sampling.ipynb
@@ -90,7 +90,7 @@
     "time = np.repeat(\n",
     "    ds_fields.time.values[0], npart\n",
     ")  # release all particles at the start time of the fieldset\n",
-    "z = np.repeat(ds_fields.depth.values[0], npart)\n",
+    "z = np.repeat(0.5, npart)  # upper level is at 0.5m depth\n",
     "\n",
     "# Plot temperature field and initial particle locations\n",
     "plt.figure()\n",
diff --git a/docs/user_guide/examples/tutorial_statuscodes.md b/docs/user_guide/examples/tutorial_statuscodes.md
index 461992b7e..708c04c65 100644
--- a/docs/user_guide/examples/tutorial_statuscodes.md
+++ b/docs/user_guide/examples/tutorial_statuscodes.md
@@ -49,7 +49,7 @@ def KeepInOcean(particles, fieldset):
     through_surface = particles.state == parcels.StatusCode.ErrorThroughSurface
 
     # move particles to surface
-    particles[through_surface].dz = fieldset.W.grid.depth[0] - particles[through_surface].z
+    particles[through_surface].dz = fieldset.surface - particles[through_surface].z
 
     # change state from error to evaluate
     particles[through_surface].state = parcels.StatusCode.Evaluate
@@ -71,6 +71,7 @@ dx, dy = 1.0 / len(ds.XG), 1.0 / len(ds.YG)
 ds["W"] = ds["U"] - 0.1 # 0.1 m/s towards the surface
 
 fieldset = parcels.FieldSet.from_sgrid_conventions(ds, mesh="flat")
+fieldset.add_context("surface", 0)  # surface is at z=0
 ```
 
 If we advect particles with the `AdvectionRK2_3D` kernel, Parcels will raise a `FieldOutOfBoundSurfaceError`:
diff --git a/docs/user_guide/getting_started/tutorial_output.ipynb b/docs/user_guide/getting_started/tutorial_output.ipynb
index e6444473d..9cfb0a91b 100644
--- a/docs/user_guide/getting_started/tutorial_output.ipynb
+++ b/docs/user_guide/getting_started/tutorial_output.ipynb
@@ -57,7 +57,7 @@
     "# Load the CopernicusMarine data in the Agulhas region from the example_datasets\n",
     "ds_fields = parcels.tutorial.open_dataset(\n",
     "    \"CopernicusMarine_data_for_Argo_tutorial/data\"\n",
-    ")\n",
+    ").isel(depth=0)  # select the first depth level (surface)\n",
     "\n",
     "# Convert to SGRID-compliant dataset and create FieldSet\n",
     "fields = {\"U\": ds_fields[\"uo\"], \"V\": ds_fields[\"vo\"]}\n",
@@ -78,11 +78,10 @@
     "npart = 10  # number of particles to be released\n",
     "lon = 32 * np.ones(npart)\n",
     "lat = np.linspace(-32.5, -30.5, npart, dtype=np.float32)\n",
-    "z = np.repeat(ds_fields.depth.values[0], npart)\n",
     "time = ds_fields.time.values[0] + np.arange(0, npart) * np.timedelta64(2, \"h\")\n",
     "\n",
     "pset = parcels.ParticleSet(\n",
-    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, z=z, t=time\n",
+    "    fieldset=fieldset, pclass=parcels.Particle, x=lon, y=lat, t=time\n",
     ")\n",
     "\n",
     "output_file = parcels.ParticleFile(\"output.parquet\", outputdt=np.timedelta64(2, \"h\"))"
diff --git a/docs/user_guide/getting_started/tutorial_quickstart.md b/docs/user_guide/getting_started/tutorial_quickstart.md
index c647173a4..896d27f09 100644
--- a/docs/user_guide/getting_started/tutorial_quickstart.md
+++ b/docs/user_guide/getting_started/tutorial_quickstart.md
@@ -89,7 +89,7 @@ npart = 10  # number of particles to be released
 lat = np.linspace(-32.5, -30.5, npart)
 lon = np.repeat(32, npart)
 time = np.repeat(ds_fields.time.values[0], npart) # at initial time of input data
-z = np.repeat(ds_fields.depth.values[0], npart) # at the first depth (surface)
+z = np.repeat(0.5, npart) # at 0.5 meter depth (first level of the input data)
 
 pset = parcels.ParticleSet(
     fieldset=fieldset, pclass=parcels.Particle, t=time, z=z, y=lat, x=lon