diff --git a/docs/user-guide/index.md b/docs/user-guide/index.md
index 24eee0c8..d7b25e4f 100644
--- a/docs/user-guide/index.md
+++ b/docs/user-guide/index.md
@@ -18,4 +18,5 @@ documentation/full_sensor_list.md
 documentation/copernicus_products.md
 documentation/pre_download_data.md
 documentation/example_copernicus_download.ipynb
+documentation/full_sensor_list.md
 ```
diff --git a/pixi.toml b/pixi.toml
index ba024968..6cd97d08 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -18,13 +18,12 @@ setuptools = "*"
 setuptools_scm = "*"
 
 [package.run-dependencies] # Keep in sync with `pyproject.toml` and feedstock recipe
-python = ">=3.10"
+python = "3.11.*"
 click = "*"
-parcels = ">3.1.0"
 pyproj = ">=3,<4"
 sortedcontainers = "==2.4.0"
 opensimplex = "==0.4.5"
-numpy = ">=1,<2"
+numpy = ">=2.1.0"
 pydantic = ">=2,<3"
 pyyaml = "*"
 copernicusmarine = ">=2.2.2"
@@ -33,15 +32,25 @@ textual = "*"
 
 [dependencies]
 virtualship = { path = "." }
+# Pre-install as conda packages to avoid PyPI source builds
+netcdf4 = "*"
+numpy = ">=2.1.0"
+dask = "*"
+zarr = ">=3"
+ipdb = ">=0.13.13,<0.14"
 
-[feature.py310.dependencies]
-python = "3.10.*"
+[pypi-dependencies]
+parcels = { git = "https://github.com/Parcels-code/Parcels", branch = "main" }
 
-[feature.py311.dependencies]
-python = "3.11.*"
+# Commented out whilst parcels v4 alpha only supports Python 3.11 (?)
+# [feature.py310.dependencies]
+# python = "3.10.*"
+
+# [feature.py311.dependencies]
+# python = "3.11.*"
 
-[feature.py312.dependencies]
-python = "3.12.*"
+# [feature.py312.dependencies]
+# python = "3.12.*"
 
 [feature.test.dependencies]
 pytest = "*"
@@ -98,11 +107,8 @@ lxml = "*"
 typing = "mypy src/virtualship --install-types"
 
 [environments]
-default = { features = ["test", "notebooks", "typing", "pre-commit", "analysis"] }
+default = { features = ["test", "notebooks", "typing", "pre-commit", "analysis"] }
 test-latest = { features = ["test"], solve-group = "test" }
-test-py310 = { features = ["test", "py310"] }
-test-py311 = { features = ["test", "py311"] }
-test-py312 = { features = ["test", "py312"] }
 test-notebooks = { features = ["test", "notebooks"], solve-group = "test" }
 analysis = { features = ["analysis"], solve-group = "analysis" }
 docs = { features = ["docs"], solve-group = "docs" }
diff --git a/pyproject.toml b/pyproject.toml
index 7f9a2108..fd0d612b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ name = "virtualship"
 description = "Code for the Virtual Ship Classroom, where Marine Scientists can combine Copernicus Marine Data with an OceanParcels ship to go on a virtual expedition."
 readme = "README.md"
 dynamic = ["version"]
-authors = [{ name = "oceanparcels.org team" }]
+authors = [{ name = "parcels-code.org team" }]
 requires-python = ">=3.10"
 license = { file = "LICENSE" }
 classifiers = [
@@ -26,11 +26,11 @@ classifiers = [
 ]
 dependencies = [
     "click",
-    "parcels >3.1.0",
+    "parcels >=4.0.0alpha",
     "pyproj >= 3, < 4",
     "sortedcontainers == 2.4.0",
     "opensimplex == 0.4.5",
-    "numpy >=1, < 2",
+    "numpy >=2.1.0",
     "pydantic >=2, <3",
     "PyYAML",
     "copernicusmarine >= 2.2.2",
@@ -40,7 +40,7 @@ dependencies = [
 ]
 
 [project.urls]
-Homepage = "https://oceanparcels.org/" # TODO: Update this to just be repo?
+Homepage = "https://virtualship.parcels-code.org/"
 Repository = "https://github.com/OceanParcels/virtualship"
 Documentation = "https://virtualship.readthedocs.io/"
 "Bug Tracker" = "https://github.com/OceanParcels/virtualship/issues"
@@ -69,7 +69,8 @@ filterwarnings = [
   "error",
   "default::DeprecationWarning",
   "error::DeprecationWarning:virtualship",
-  "ignore:ParticleSet is empty.*:RuntimeWarning" # TODO: Probably should be ignored in the source code
+  "ignore:ParticleSet is empty.*:RuntimeWarning", # TODO: Probably should be ignored in the source code
+  "ignore:This is an alpha version of Parcels v4.*:UserWarning" # TODO: necessary whilst Parcels v4 is still alpha
 ]
 log_cli_level = "INFO"
 testpaths = [
diff --git a/src/virtualship/cli/_run.py b/src/virtualship/cli/_run.py
index f2622be3..10afc5c5 100644
--- a/src/virtualship/cli/_run.py
+++ b/src/virtualship/cli/_run.py
@@ -35,11 +35,9 @@
     get_instrument_class,
 )
 
-# parcels logger (suppress INFO messages to prevent log being flooded)
-external_logger = logging.getLogger("parcels.tools.loggers")
-external_logger.setLevel(logging.WARNING)
-
-# copernicusmarine logger (suppress INFO messages to prevent log being flooded)
+# suppress INFO messages from copernicusmarine and parcels loggers; prevent log flooding
+parcels_logger = logging.getLogger("parcels._logger")
+parcels_logger.setLevel(logging.WARNING)
 logging.getLogger("copernicusmarine").setLevel("ERROR")
 
 
@@ -202,6 +200,7 @@ def _run(
             )
 
             # execute simulation
+            # TODO: outpath will be Parquet with v4...
             instrument.execute(
                 measurements=measurements,
                 out_path=expedition_dir.joinpath(RESULTS, f"{itype.name.lower()}.zarr"),
diff --git a/src/virtualship/instruments/adcp.py b/src/virtualship/instruments/adcp.py
index b2da6582..626ca354 100644
--- a/src/virtualship/instruments/adcp.py
+++ b/src/virtualship/instruments/adcp.py
@@ -3,7 +3,7 @@
 from typing import ClassVar
 
 import numpy as np
-from parcels import ParticleSet, ScipyParticle
+from parcels import ParticleSet
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
@@ -35,9 +35,13 @@ class ADCP:
 # =====================================================
 
 
-def _sample_velocity(particle, fieldset, time):
-    particle.U, particle.V = fieldset.UV.eval(
-        time, particle.depth, particle.lat, particle.lon, applyConversion=False
+def _sample_velocity(particles, fieldset):
+    particles.U, particles.V = fieldset.UV.eval(
+        particles.time,
+        particles.z,
+        particles.lat,
+        particles.lon,
+        applyConversion=False,
     )
 
 
@@ -96,7 +100,7 @@ def simulate(self, measurements, out_path) -> None:
         # build dynamic particle class from the active sensors
         adcp_config = self.expedition.instruments_config.adcp_config
         _ADCPParticle = build_particle_class_from_sensors(
-            adcp_config.sensors, _ADCP_NONSENSOR_VARIABLES, ScipyParticle
+            adcp_config.sensors, _ADCP_NONSENSOR_VARIABLES
         )
 
         bins = np.linspace(MAX_DEPTH, MIN_DEPTH, NUM_BINS)
diff --git a/src/virtualship/instruments/argo_float.py b/src/virtualship/instruments/argo_float.py
index 8c90cfb2..4eae6350 100644
--- a/src/virtualship/instruments/argo_float.py
+++ b/src/virtualship/instruments/argo_float.py
@@ -1,11 +1,11 @@
-import math
 from collections.abc import Callable
 from dataclasses import dataclass
 from datetime import timedelta
 from typing import ClassVar
 
 import numpy as np
-from parcels import AdvectionRK4, JITParticle, ParticleSet, StatusCode, Variable
+from parcels import ParticleSet, StatusCode, Variable
+from parcels.kernels import AdvectionRK2
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
@@ -53,103 +53,125 @@ class ArgoFloat:
 # SECTION: Kernels
 # =====================================================
 
-
-def _argo_float_vertical_movement(particle, fieldset, time):
-    if particle.cycle_phase == 0:
-        # Phase 0: Sinking with vertical_speed until depth is drift_depth
-        particle_ddepth += (  # noqa
-            particle.vertical_speed * particle.dt
+# TODO: need to add back in the shallow bathymetry checks (to phases 0 and 2?!)
+# TODO: can this be refactored as well to a helper function?
+
+
+def _argo_float_vertical_movement(particles, fieldset):
+    # Split particles based on their current cycle_phase
+    ptcls0 = particles[particles.cycle_phase == 0]
+    ptcls1 = particles[particles.cycle_phase == 1]
+    ptcls2 = particles[particles.cycle_phase == 2]
+    ptcls3 = particles[particles.cycle_phase == 3]
+    ptcls4 = particles[particles.cycle_phase == 4]
+
+    # Phase 0: Sinking with vertical_speed until depth is driftdepth
+    ptcls0.dz += particles.vertical_speed * ptcls0.dt
+    loc_bathy = fieldset.bathymetry.eval(ptcls0.time, ptcls0.z, ptcls0.lat, ptcls0.lon)
+    driftdepth_mask = ptcls0.z + ptcls0.dz >= particles.drift_depth
+    bathy_mask = ptcls0.z + ptcls0.dz >= loc_bathy
+    next_phase = np.logical_and(
+        driftdepth_mask, bathy_mask
+    )  # combined mask; not at drift depth yet and not hitting bathymetry
+    ptcls0.cycle_phase[next_phase] = 1
+    ptcls0.dz[next_phase] = (
+        particles.drift_depth - ptcls0.z[next_phase]
+    )  # avoid overshoot
+
+    # Phase 0.5: Check for grounding at bathymetry and raise if necessary
+    ptcls0.grounded[~bathy_mask] = 1
+    if np.any(~bathy_mask):
+        print(
+            "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to drift depth. Raising by 50m above bathymetry and continuing cycle."
         )
-
-        # bathymetry at particle location
-        loc_bathy = fieldset.bathymetry.eval(
-            time, particle.depth, particle.lat, particle.lon
+        ptcls0.dz[~bathy_mask] = (
+            loc_bathy[~bathy_mask] - ptcls0.z[~bathy_mask] + 50.0
+        )  # raise to 50m above bathymetry
+        ptcls0.cycle_phase[~bathy_mask] = 1
+
+    # Phase 1: Drifting at depth for drifttime seconds
+    ptcls1.drift_age += ptcls1.dt
+    next_phase = ptcls1.drift_age >= particles.drift_days * 86400  # [seconds]
+    ptcls1.cycle_phase[next_phase] = 2
+    ptcls1.drift_age[next_phase] = 0  # reset drift_age for next cycle
+
+    # Phase 2: Sinking further to maxdepth
+    ptcls2.dz += particles.vertical_speed * ptcls2.dt
+    loc_bathy = fieldset.bathymetry.eval(ptcls2.time, ptcls2.z, ptcls2.lat, ptcls2.lon)
+    maxdepth_mask = ptcls2.z + ptcls2.dz >= particles.max_depth
+    bathy_mask = ptcls2.z + ptcls2.dz >= loc_bathy
+    next_phase = np.logical_and(
+        maxdepth_mask, bathy_mask
+    )  # combined mask; not at max depth yet and not hitting bathymetry
+    ptcls2.cycle_phase[next_phase] = 3
+    ptcls2.dz[next_phase] = (
+        particles.max_depth - ptcls2.z[next_phase]
+    )  # avoid overshoot
+
+    # Phase 2.5: Check for grounding at bathymetry and raise if necessary
+    ptcls2.grounded[~bathy_mask] = 1
+    if np.any(~bathy_mask):
+        print(
+            "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to max depth. Raising by 50m above bathymetry and continuing cycle."
         )
-        if particle.depth + particle_ddepth <= loc_bathy:
-            particle_ddepth = loc_bathy - particle.depth + 50.0  # 50m above bathy
-            particle.cycle_phase = 1
-            particle.grounded = 1
-            print(
-                "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to drift depth. Raising by 50m above bathymetry and continuing cycle."
-            )
+        ptcls2.dz[~bathy_mask] = (
+            loc_bathy[~bathy_mask] - ptcls2.z[~bathy_mask] + 50.0
+        )  # raise to 50m above bathymetry
+        ptcls2.cycle_phase[~bathy_mask] = 3
 
-        elif particle.depth + particle_ddepth <= particle.drift_depth:
-            particle_ddepth = particle.drift_depth - particle.depth
-            particle.cycle_phase = 1
-
-    elif particle.cycle_phase == 1:
-        # Phase 1: Drifting at depth for drifttime seconds
-        particle.drift_age += particle.dt
-        if particle.drift_age >= particle.drift_days * 86400:
-            particle.drift_age = 0  # reset drift_age for next cycle
-            particle.cycle_phase = 2
-
-    elif particle.cycle_phase == 2:
-        # Phase 2: Sinking further to max_depth
-        particle_ddepth += particle.vertical_speed * particle.dt
-        loc_bathy = fieldset.bathymetry.eval(
-            time, particle.depth, particle.lat, particle.lon
-        )
-        if particle.depth + particle_ddepth <= loc_bathy:
-            particle_ddepth = loc_bathy - particle.depth + 50.0  # 50m above bathy
-            particle.cycle_phase = 3
-            particle.grounded = 1
-            print(
-                "Shallow bathymetry warning: Argo float grounded at bathymetry depth during sinking to max depth. Raising by 50m above bathymetry and continuing cycle."
-            )
-        elif particle.depth + particle_ddepth <= particle.max_depth:
-            particle_ddepth = particle.max_depth - particle.depth
-            particle.cycle_phase = 3
-
-    elif particle.cycle_phase == 3:
-        # Phase 3: Rising with vertical_speed until at surface
-        particle_ddepth -= particle.vertical_speed * particle.dt
-        particle.cycle_age += (
-            particle.dt
-        )  # solve issue of not updating cycle_age during ascent
-        particle.grounded = 0
-        if particle.depth + particle_ddepth >= particle.min_depth:
-            particle_ddepth = particle.min_depth - particle.depth
-            particle.cycle_phase = 4
+    # Phase 3: Rising with vertical_speed until at surface
+    ptcls3.dz -= particles.vertical_speed * ptcls3.dt
+    ptcls3.temp = fieldset.thetao[ptcls3.time, ptcls3.z, ptcls3.lat, ptcls3.lon]
+    next_phase = ptcls3.z + ptcls3.dz <= particles.min_depth
+    ptcls3.cycle_phase[next_phase] = 4
+    ptcls3.dz[next_phase] = (
+        particles.min_depth - ptcls3.z[next_phase]
+    )  # avoid overshoot
 
-    elif particle.cycle_phase == 4:
-        # Phase 4: Transmitting at surface until cycletime is reached
-        if particle.cycle_age > particle.cycle_days * 86400:
-            particle.cycle_phase = 0
-            particle.cycle_age = 0
+    # Phase 4: Transmitting at surface until cycletime is reached
+    next_phase = ptcls4.cycle_age >= particles.cycle_days * 86400
+    ptcls4.cycle_phase[next_phase] = 0
+    ptcls4.cycle_age[next_phase] = 0  # reset cycle_age for next cycle
+    ptcls4.temp = np.nan  # no temperature measurement when at surface
 
-    if particle.state == StatusCode.Evaluate:
-        particle.cycle_age += particle.dt  # update cycle_age
+    particles.cycle_age += particles.dt  # update cycle_age
 
 
-def _keep_at_surface(particle, fieldset, time):
-    # Prevent error when float reaches surface
-    if particle.state == StatusCode.ErrorThroughSurface:
-        particle.depth = particle.min_depth
-        particle.state = StatusCode.Success
+def _keep_at_surface(particles, fieldset):
+    through_surface = particles.state == StatusCode.ErrorThroughSurface
+    particles.z[through_surface] = particles.min_depth[through_surface]
+    particles.state[through_surface] = StatusCode.Success
 
 
-def _check_error(particle, fieldset, time):
-    if particle.state >= 50:  # This captures all Errors
-        particle.delete()
+def _check_error(particles, fieldset):
+    errors = particles.state >= 50  # captures all Errors
+    particles.state[errors] = StatusCode.Delete
 
 
-def _argo_sample_temperature(particle, fieldset, time):
+def _argo_sample_temperature(particles, fieldset):
     # Phase 3: ascending — sample temperature; NaN otherwise
-    if particle.cycle_phase == 3 and particle.depth < particle.min_depth:
-        particle.temperature = fieldset.T[
-            time, particle.depth, particle.lat, particle.lon
-        ]
-    else:
-        particle.temperature = math.nan
-
-
-def _argo_sample_salinity(particle, fieldset, time):
+    phase_mask = particles.cycle_phase == 3
+    depth_mask = particles.depth < particles.min_depth
+    sampling_particles = particles[np.logical_and(phase_mask, depth_mask)]
+    sampling_particles.temperature = fieldset.T[
+        sampling_particles.time,
+        sampling_particles.depth,
+        sampling_particles.lat,
+        sampling_particles.lon,
+    ]
+
+
+def _argo_sample_salinity(particles, fieldset):
     # Phase 3: ascending — sample salinity; NaN otherwise
-    if particle.cycle_phase == 3 and particle.depth < particle.min_depth:
-        particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon]
-    else:
-        particle.salinity = math.nan
+    phase_mask = particles.cycle_phase == 3
+    depth_mask = particles.depth < particles.min_depth
+    sampling_particles = particles[np.logical_and(phase_mask, depth_mask)]
+    sampling_particles.salinity = fieldset.S[
+        sampling_particles.time,
+        sampling_particles.depth,
+        sampling_particles.lat,
+        sampling_particles.lon,
+    ]
 
 
 # =====================================================
@@ -229,9 +251,7 @@ def simulate(self, measurements, out_path) -> None:
         # build dynamic particle class from the active sensors
         argo_float_config = self.expedition.instruments_config.argo_float_config
         _ArgoParticle = build_particle_class_from_sensors(
-            argo_float_config.sensors,
-            _ARGO_NONSENSOR_VARIABLES,
-            JITParticle,
+            argo_float_config.sensors, _ARGO_NONSENSOR_VARIABLES
         )
 
         # define parcel particles
@@ -272,7 +292,7 @@ def simulate(self, measurements, out_path) -> None:
             [
                 _argo_float_vertical_movement,
                 *sampling_kernels,
-                AdvectionRK4,
+                AdvectionRK2,
                 _keep_at_surface,
                 _check_error,
             ],
diff --git a/src/virtualship/instruments/base.py b/src/virtualship/instruments/base.py
index d4e078e6..f2ad8a47 100644
--- a/src/virtualship/instruments/base.py
+++ b/src/virtualship/instruments/base.py
@@ -8,8 +8,8 @@
 from typing import TYPE_CHECKING, ClassVar
 
 import copernicusmarine
+import parcels
 import xarray as xr
-from parcels import FieldSet
 from yaspin import yaspin
 
 from virtualship.errors import CopernicusCatalogueError
@@ -86,7 +86,7 @@ def __init__(
         self.min_lat, self.max_lat = min(wp_lats), max(wp_lats)
         self.min_lon, self.max_lon = min(wp_lons), max(wp_lons)
 
-    def load_input_data(self) -> FieldSet:
+    def load_input_data(self) -> parcels.FieldSet:
         """Load and return the input data as a FieldSet for the instrument."""
         try:
             fieldset = self._generate_fieldset()
@@ -97,22 +97,14 @@ def load_input_data(self) -> FieldSet:
 
         # interpolation methods
         for var in (v for v in self.variables if v not in ("U", "V")):
-            getattr(fieldset, var).interp_method = "linear_invdist_land_tracer"
-
-        # depth negative
-        for g in fieldset.gridset.grids:
-            g.negate_depth()
+            getattr(fieldset, var).interp_method = parcels.interpolators.XLinear
 
         # bathymetry data
         if self.add_bathymetry:
-            bathymetry_field = _get_bathy_data(
-                self.min_lat,
-                self.max_lat,
-                self.min_lon,
-                self.max_lon,
-                from_data=self.from_data,
-            ).bathymetry
-            bathymetry_field.data = -bathymetry_field.data
+            bathymetry_field = _get_bathy_data(from_data=self.from_data).bathymetry
+            bathymetry_field.data = (
+                -bathymetry_field.data
+            )  # TODO: how does v4 handle? positive up or down?
             fieldset.add_field(bathymetry_field)
 
         return fieldset
@@ -183,11 +175,11 @@ def _get_copernicus_ds(
             coordinates_selection_method="outside",
         )
 
-    def _generate_fieldset(self) -> FieldSet:
+    def _generate_fieldset(self) -> parcels.FieldSet:
         """
         Create and combine FieldSets for each variable, supporting both local and Copernicus Marine data sources.
 
-        Per variable avoids issues when using copernicusmarine and creating directly one FieldSet of ds's sourced from different Copernicus Marine product IDs, which is often the case for BGC variables.
+        N.B. Per variable avoids issues when using copernicusmarine and creating directly one FieldSet of ds's sourced from different Copernicus Marine product IDs (which can also have different temporal resolutions), which is often the case for BGC variables.
         """
         fieldsets_list = []
         keys = list(self.variables.keys())
@@ -217,12 +209,11 @@ def _generate_fieldset(self) -> FieldSet:
                     [data_dir.joinpath(f) for f in files]
                 )  # using: ds --> .from_xarray_dataset seems more robust than .from_netcdf for handling different temporal resolutions for different variables ...
 
-                fs = FieldSet.from_xarray_dataset(
-                    ds,
-                    variables={key: full_var_name},
-                    dimensions=self.dimensions,
-                    mesh="spherical",
-                )
+                # TODO: do docs on pre-downloading data need to be updated for these changes? Anything about conventions etc.?
+                fields = {key: ds[full_var_name]}
+                ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
+                fs = parcels.FieldSet.from_sgrid_conventions(ds_fset)
+
             else:  # stream via Copernicus Marine Service
                 physical = var in COPERNICUSMARINE_PHYS_VARIABLES
                 ds = self._get_copernicus_ds(
@@ -230,9 +221,11 @@ def _generate_fieldset(self) -> FieldSet:
                     physical=physical,
                     var=var,
                 )
-                fs = FieldSet.from_xarray_dataset(
-                    ds, {key: var}, self.dimensions, mesh="spherical"
-                )
+
+                fields = {key: ds[var]}
+                ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
+                fs = parcels.FieldSet.from_sgrid_conventions(ds_fset)
+
             fieldsets_list.append(fs)
 
         base_fieldset = fieldsets_list[0]
diff --git a/src/virtualship/instruments/ctd.py b/src/virtualship/instruments/ctd.py
index 583a099c..b48677ad 100644
--- a/src/virtualship/instruments/ctd.py
+++ b/src/virtualship/instruments/ctd.py
@@ -4,13 +4,13 @@
 from typing import TYPE_CHECKING, ClassVar
 
 import numpy as np
-from parcels import JITParticle, ParticleSet, Variable
+from parcels import ParticleSet, Variable
+from parcels._core.statuscodes import StatusCode
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
 from virtualship.utils import (
-    add_dummy_UV,
     build_particle_class_from_sensors,
     register_instrument,
 )
@@ -52,60 +52,90 @@ class CTD:
 ## physical variables
 
 
-def _sample_temperature(particle, fieldset, time):
-    particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
+def _sample_temperature(particles, fieldset):
+    particles.temperature = fieldset.T[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_salinity(particle, fieldset, time):
-    particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon]
+def _sample_salinity(particles, fieldset):
+    particles.salinity = fieldset.S[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
 ## bgc variables
 
 
-def _sample_o2(particle, fieldset, time):
-    particle.o2 = fieldset.o2[time, particle.depth, particle.lat, particle.lon]
+def _sample_o2(particles, fieldset):
+    particles.o2 = fieldset.o2[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_chlorophyll(particle, fieldset, time):
-    particle.chl = fieldset.chl[time, particle.depth, particle.lat, particle.lon]
+def _sample_chlorophyll(particles, fieldset):
+    particles.chl = fieldset.chl[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_nitrate(particle, fieldset, time):
-    particle.no3 = fieldset.no3[time, particle.depth, particle.lat, particle.lon]
+def _sample_nitrate(particles, fieldset):
+    particles.no3 = fieldset.no3[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_phosphate(particle, fieldset, time):
-    particle.po4 = fieldset.po4[time, particle.depth, particle.lat, particle.lon]
+def _sample_phosphate(particles, fieldset):
+    particles.po4 = fieldset.po4[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_ph(particle, fieldset, time):
-    particle.ph = fieldset.ph[time, particle.depth, particle.lat, particle.lon]
+def _sample_ph(particles, fieldset):
+    particles.ph = fieldset.ph[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_phytoplankton(particle, fieldset, time):
-    particle.phyc = fieldset.phyc[time, particle.depth, particle.lat, particle.lon]
+def _sample_phytoplankton(particles, fieldset):
+    particles.phyc = fieldset.phyc[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _sample_primary_production(particle, fieldset, time):
-    particle.nppv = fieldset.nppv[time, particle.depth, particle.lat, particle.lon]
+def _sample_primary_production(particles, fieldset):
+    particles.nppv = fieldset.nppv[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
 ## cast
 
 
-def _ctd_cast(particle, fieldset, time):
+def _ctd_cast(particles, fieldset):
+    particles_lowering = particles[particles.raising == 0]
+    particles_raising = particles[particles.raising == 1]
+
+    # TODO: change to boolean masking, like with Argo Floats?
+    # TODO: different handling of positive down for z now?! Doing positive down now... think kernels need adjusting...
+    # TODO: need to check on all other instrument kernels as well...
+    # TODO: plus how the configs are inputted in e.g. expedition.yaml
+
     # lowering
-    if particle.raising == 0:
-        particle_ddepth = -particle.winch_speed * particle.dt
-        if particle.depth + particle_ddepth < particle.max_depth:
-            particle.raising = 1
-            particle_ddepth = -particle_ddepth
+    particles_lowering.dz = -particles_lowering.winch_speed * particles_lowering.dt
+    particles_lowering.raising = np.where(
+        particles_lowering.z + particles_lowering.dz < particles_lowering.max_depth,
+        1,
+        particles_lowering.raising,
+    )
+
     # raising
-    else:
-        particle_ddepth = particle.winch_speed * particle.dt
-        if particle.depth + particle_ddepth > particle.min_depth:
-            particle.delete()
+    particles_raising.dz = particles_raising.winch_speed * particles_raising.dt
+    particles_raising.state = np.where(
+        particles_raising.z + particles_raising.dz > particles_raising.min_depth,
+        StatusCode.Delete,
+        particles_raising.state,
+    )
 
 
 # =====================================================
@@ -162,9 +192,6 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # add dummy U
-        add_dummy_UV(fieldset)  # TODO: parcels v3 bodge; remove when parcels v4 is used
-
         # use first active field for time reference
         _time_ref_key = next(iter(self.variables))
         _time_ref_field = getattr(fieldset, _time_ref_key)
@@ -208,7 +235,7 @@ def simulate(self, measurements, out_path) -> None:
         # build dynamic particle class from the active sensors
         ctd_config = self.expedition.instruments_config.ctd_config
         _CTDParticle = build_particle_class_from_sensors(
-            ctd_config.sensors, _CTD_NONSENSOR_VARIABLES, JITParticle
+            ctd_config.sensors, _CTD_NONSENSOR_VARIABLES
         )
 
         # define parcel particles
diff --git a/src/virtualship/instruments/drifter.py b/src/virtualship/instruments/drifter.py
index 379334b3..79a1e34b 100644
--- a/src/virtualship/instruments/drifter.py
+++ b/src/virtualship/instruments/drifter.py
@@ -4,7 +4,9 @@
 from typing import ClassVar
 
 import numpy as np
-from parcels import AdvectionRK4, JITParticle, ParticleSet, Variable
+from parcels import ParticleSet, Variable
+from parcels._core.statuscodes import StatusCode
+from parcels.kernels import AdvectionRK2
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
@@ -46,15 +48,21 @@ class Drifter:
 # =====================================================
 
 
-def _sample_temperature(particle, fieldset, time):
-    particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
+def _sample_temperature(particles, fieldset):
+    particles.temperature = fieldset.T[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _check_lifetime(particle, fieldset, time):
-    if particle.has_lifetime == 1:
-        particle.age += particle.dt
-        if particle.age >= particle.lifetime:
-            particle.delete()
+def _check_lifetime(particles, fieldset):
+    particles_wlifetime = particles[particles.has_lifetime == 1]
+
+    particles_wlifetime.age += particles_wlifetime.dt
+    particles_wlifetime.state = np.where(
+        particles_wlifetime.age >= particles_wlifetime.lifetime,
+        StatusCode.Delete,
+        particles_wlifetime.state,
+    )
 
 
 # =====================================================
@@ -123,7 +131,7 @@ def simulate(self, measurements, out_path) -> None:
         # build dynamic particle class from the active sensors
         drifter_config = self.expedition.instruments_config.drifter_config
         _DrifterParticle = build_particle_class_from_sensors(
-            drifter_config.sensors, _DRIFTER_NONSENSOR_VARIABLES, JITParticle
+            drifter_config.sensors, _DRIFTER_NONSENSOR_VARIABLES
         )
 
         # define parcel particles
@@ -169,7 +177,7 @@ def simulate(self, measurements, out_path) -> None:
 
         # execute simulation
         drifter_particleset.execute(
-            [AdvectionRK4, *sampling_kernels, _check_lifetime],
+            [AdvectionRK2, *sampling_kernels, _check_lifetime],
             endtime=endtime,
             dt=DT,
             output_file=out_file,
diff --git a/src/virtualship/instruments/ship_underwater_st.py b/src/virtualship/instruments/ship_underwater_st.py
index 6a564cc0..d844b7a2 100644
--- a/src/virtualship/instruments/ship_underwater_st.py
+++ b/src/virtualship/instruments/ship_underwater_st.py
@@ -3,13 +3,12 @@
 from typing import ClassVar
 
 import numpy as np
-from parcels import ParticleSet, ScipyParticle
+from parcels import ParticleSet
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
 from virtualship.utils import (
-    add_dummy_UV,
     build_particle_class_from_sensors,
     register_instrument,
 )
@@ -40,13 +39,13 @@ class Underwater_ST:
 
 
 # define function sampling Salinity
-def _sample_salinity(particle, fieldset, time):
-    particle.salinity = fieldset.S[time, particle.depth, particle.lat, particle.lon]
+def _sample_salinity(particles, fieldset):
+    particles.S = fieldset.S[particles.time, particles.z, particles.lat, particles.lon]
 
 
 # define function sampling Temperature
-def _sample_temperature(particle, fieldset, time):
-    particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
+def _sample_temperature(particles, fieldset):
+    particles.T = fieldset.T[particles.time, particles.z, particles.lat, particles.lon]
 
 
 # =====================================================
@@ -95,13 +94,10 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # add dummy U
-        add_dummy_UV(fieldset)  # TODO: parcels v3 bodge; remove when parcels v4 is used
-
         # build dynamic particle class from the active sensors
         st_config = self.expedition.instruments_config.ship_underwater_st_config
         _ShipSTParticle = build_particle_class_from_sensors(
-            st_config.sensors, _ST_NONSENSOR_VARIABLES, ScipyParticle
+            st_config.sensors, _ST_NONSENSOR_VARIABLES
         )
 
         particleset = ParticleSet.from_list(
diff --git a/src/virtualship/instruments/xbt.py b/src/virtualship/instruments/xbt.py
index 051bf1fa..d02afb7a 100644
--- a/src/virtualship/instruments/xbt.py
+++ b/src/virtualship/instruments/xbt.py
@@ -4,14 +4,14 @@
 from typing import ClassVar
 
 import numpy as np
-from parcels import JITParticle, ParticleSet, Variable
+from parcels import ParticleSet, Variable
+from parcels._core.statuscodes import StatusCode
 
 from virtualship.instruments.base import Instrument
 from virtualship.instruments.sensors import SensorType
 from virtualship.instruments.types import InstrumentType
 from virtualship.models.spacetime import Spacetime
 from virtualship.utils import (
-    add_dummy_UV,
     build_particle_class_from_sensors,
     register_instrument,
 )
@@ -50,26 +50,32 @@ class XBT:
 # =====================================================
 
 
-def _sample_temperature(particle, fieldset, time):
-    particle.temperature = fieldset.T[time, particle.depth, particle.lat, particle.lon]
+def _sample_temperature(particles, fieldset):
+    particles.temperature = fieldset.T[
+        particles.time, particles.z, particles.lat, particles.lon
+    ]
 
 
-def _xbt_cast(particle, fieldset, time):
-    particle_ddepth = -particle.fall_speed * particle.dt
+def _xbt_cast(particles, fieldset):
+    particles.dz = -particles.fall_speed * particles.dt
 
     # update the fall speed from the quadractic fall-rate equation
     # check https://doi.org/10.5194/os-7-231-2011
-    particle.fall_speed = (
-        particle.fall_speed - 2 * particle.deceleration_coefficient * particle.dt
+    particles.fall_speed = (
+        particles.fall_speed - 2 * particles.deceleration_coefficient * particles.dt
     )
 
     # delete particle if depth is exactly max_depth
-    if particle.depth == particle.max_depth:
-        particle.delete()
+    particles.state = np.where(
+        particles.z == particles.max_depth, StatusCode.Delete, particles.state
+    )
 
     # set particle depth to max depth if it's too deep
-    if particle.depth + particle_ddepth < particle.max_depth:
-        particle_ddepth = particle.max_depth - particle.depth
+    particles.dz = np.where(
+        particles.z + particles.dz < particles.max_depth,
+        particles.max_depth - particles.z,
+        particles.z,
+    )
 
 
 # =====================================================
@@ -117,9 +123,6 @@ def simulate(self, measurements, out_path) -> None:
 
         fieldset = self.load_input_data()
 
-        # add dummy U
-        add_dummy_UV(fieldset)  # TODO: parcels v3 bodge; remove when parcels v4 is used
-
         # use first active field for time reference
         _time_ref_key = next(iter(self.variables))
         _time_ref_field = getattr(fieldset, _time_ref_key)
@@ -166,7 +169,7 @@ def simulate(self, measurements, out_path) -> None:
         # build dynamic particle class from the active sensors
         xbt_config = self.expedition.instruments_config.xbt_config
         _XBTParticle = build_particle_class_from_sensors(
-            xbt_config.sensors, _XBT_NONSENSOR_VARIABLES, JITParticle
+            xbt_config.sensors, _XBT_NONSENSOR_VARIABLES
         )
 
         # define xbt particles
diff --git a/src/virtualship/models/expedition.py b/src/virtualship/models/expedition.py
index eef23c76..d212a686 100644
--- a/src/virtualship/models/expedition.py
+++ b/src/virtualship/models/expedition.py
@@ -17,7 +17,6 @@
     _calc_sail_time,
     _calc_wp_stationkeeping_time,
     _get_bathy_data,
-    _get_waypoint_latlons,
     _validate_numeric_to_timedelta,
     get_supported_sensors,
     register_instrument_config,
@@ -131,14 +130,7 @@ def verify(
         land_waypoints = []
         if not ignore_land_test:
             try:
-                wp_lats, wp_lons = _get_waypoint_latlons(self.waypoints)
-                bathymetry_field = _get_bathy_data(
-                    min(wp_lats),
-                    max(wp_lats),
-                    min(wp_lons),
-                    max(wp_lons),
-                    from_data=from_data,
-                ).bathymetry
+                bathymetry_field = _get_bathy_data(from_data=from_data).bathymetry
             except Exception as e:
                 raise ScheduleError(
                     f"Problem loading bathymetry data (used to verify waypoints are in water) directly via copernicusmarine. \n\n original message: {e}"
@@ -147,7 +139,7 @@ def verify(
             for wp_i, wp in enumerate(self.waypoints):
                 try:
                     value = bathymetry_field.eval(
-                        0,  # time
+                        np.float64(0.0),  # time
                         0,  # depth (surface)
                         wp.location.lat,
                         wp.location.lon,
diff --git a/src/virtualship/utils.py b/src/virtualship/utils.py
index 7ad275cd..78657de5 100644
--- a/src/virtualship/utils.py
+++ b/src/virtualship/utils.py
@@ -13,9 +13,10 @@
 
 import copernicusmarine
 import numpy as np
+import parcels
 import pyproj
 import xarray as xr
-from parcels import FieldSet, Variable
+from parcels import FieldSet, Particle, Variable
 
 from virtualship.errors import CopernicusCatalogueError
 
@@ -340,29 +341,6 @@ def _get_expedition(expedition_dir: Path) -> Expedition:
         ) from e
 
 
-def add_dummy_UV(fieldset: FieldSet):
-    """Add a dummy U and V field to a FieldSet to satisfy parcels FieldSet completeness checks."""
-    if "U" not in fieldset.__dict__.keys():
-        for uv_var in ["U", "V"]:
-            dummy_field = getattr(
-                FieldSet.from_data(
-                    {"U": 0, "V": 0}, {"lon": 0, "lat": 0}, mesh="spherical"
-                ),
-                uv_var,
-            )
-            fieldset.add_field(dummy_field)
-            try:
-                fieldset.time_origin = (
-                    fieldset.T.grid.time_origin
-                    if "T" in fieldset.__dict__.keys()
-                    else fieldset.o2.grid.time_origin
-                )
-            except Exception:
-                raise ValueError(
-                    "Cannot determine time_origin for dummy UV fields. Assert T or o2 exists in fieldset."
-                ) from None
-
-
 def _select_product_id(
     physical: bool,
     schedule_start,
@@ -448,43 +426,36 @@ def _start_end_in_product_timerange(
     )
 
 
-def _get_bathy_data(
-    min_lat: float,
-    max_lat: float,
-    min_lon: float,
-    max_lon: float,
-    from_data: Path | None = None,
-) -> FieldSet:
+def _get_bathy_data(from_data: Path | None = None) -> FieldSet:
     """Bathymetry data from local or 'streamed' directly from Copernicus Marine."""
+    VAR = "deptho"
     if from_data is not None:  # load from local data
-        var = "deptho"
         bathy_dir = from_data.joinpath("bathymetry")
         try:
-            filename, _ = _find_nc_file_with_variable(bathy_dir, var)
+            filename, _ = _find_nc_file_with_variable(bathy_dir, VAR)
         except Exception as e:
             raise RuntimeError(
-                f"\n\n❗️ Could not find bathymetry variable '{var}' in data directory '{from_data}/bathymetry/'.\n\n❗️ Is the pre-downloaded data directory structure compliant with VirtualShip expectations?\n\n❗️ See the docs for more information on expectations: https://virtualship.readthedocs.io/en/latest/user-guide/index.html#documentation\n"
+                f"\n\n❗️ Could not find bathymetry variable '{VAR}' in data directory '{from_data}/bathymetry/'.\n\n❗️ Is the pre-downloaded data directory structure compliant with VirtualShip expectations?\n\n❗️ See the docs for more information on expectations: https://virtualship.readthedocs.io/en/latest/user-guide/index.html#documentation\n"
             ) from e
         ds_bathymetry = xr.open_dataset(bathy_dir.joinpath(filename))
-        bathymetry_variables = {"bathymetry": "deptho"}
-        bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"}
-        return FieldSet.from_xarray_dataset(
-            ds_bathymetry, bathymetry_variables, bathymetry_dimensions
-        )
 
     else:  # stream via Copernicus Marine Service
         ds_bathymetry = copernicusmarine.open_dataset(
             dataset_id=BATHYMETRY_ID,
-            variables=["deptho"],
+            variables=[VAR],
             coordinates_selection_method="outside",
         )
-        bathymetry_variables = {"bathymetry": "deptho"}
-        bathymetry_dimensions = {"lon": "longitude", "lat": "latitude"}
 
-    return FieldSet.from_xarray_dataset(
-        ds_bathymetry, bathymetry_variables, bathymetry_dimensions
+    ds_bathymetry = ds_bathymetry.expand_dims(
+        {"depth": 1}
+    )  # TODO: bodge whilst parcels v4 does not support 2D fields and seeks depth dim; change when parcels v4 released
+
+    ds_fset = parcels.convert.copernicusmarine_to_sgrid(
+        fields={"bathymetry": ds_bathymetry[VAR]}
     )
 
+    return FieldSet.from_sgrid_conventions(ds_fset)
+
 
 def expedition_cost(schedule_results: ScheduleOk, time_past: timedelta) -> float:
     """
@@ -677,14 +648,13 @@ def _make_hash(s: str, length: int) -> str:
 def build_particle_class_from_sensors(
     sensors: list[SensorConfig],
     nonsensor_variables: list[Variable],
-    particle_class: type,  # generic type annotation needed for v3 particle class behaviour # TODO: Update with Parcels v4
 ) -> type:
-    """Build a Particle class (JITParticle or ScipyParticle) from nonsensor variables and active sensors."""
+    """Build a Particle class from nonsensor variables and active sensors."""
     sensor_variables = [
         variable for sc in sensors if sc.enabled for variable in sc.meta.particle_vars
     ]
 
-    return particle_class.add_variables(nonsensor_variables + sensor_variables)
+    return Particle.add_variables(nonsensor_variables + sensor_variables)
 
 
 # =====================================================
diff --git a/tests/instruments/test_base.py b/tests/instruments/test_base.py
index bbcfea44..a17f95bf 100644
--- a/tests/instruments/test_base.py
+++ b/tests/instruments/test_base.py
@@ -33,7 +33,6 @@ def test_load_input_data(mock_copernicusmarine, mock_select_product_id, mock_Fie
     mock_fieldset = MagicMock()
     mock_FieldSet.from_netcdf.return_value = mock_fieldset
     mock_FieldSet.from_xarray_dataset.return_value = mock_fieldset
-    mock_fieldset.gridset.grids = [MagicMock(negate_depth=MagicMock())]
     mock_fieldset.__getitem__.side_effect = lambda k: MagicMock()
     mock_copernicusmarine.open_dataset.return_value = MagicMock()
     # Create a mock waypoint with latitude and longitude