diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py
index 988636f13..b2ad6cb96 100644
--- a/src/parcels/_core/field.py
+++ b/src/parcels/_core/field.py
@@ -23,6 +23,7 @@
 from parcels._core.utils.time import TimeInterval
 from parcels._core.uxgrid import UxGrid
 from parcels._core.xgrid import XGrid, _transpose_xfield_data_to_tzyx, assert_all_field_dims_have_axis
+from parcels._logger import logger
 from parcels._python import assert_same_function_signature
 from parcels._reprs import default_repr
 from parcels._typing import VectorType
@@ -137,6 +138,14 @@ def __init__(
         assert_same_function_signature(interp_method, ref=ZeroInterpolator, context="Interpolation")
         self._interp_method = interp_method
 
+        # Setting the land value if present
+        print("DATA ATTRS:", data.attrs)
+        if "_FillValue" in data.attrs:
+            self._land_value = data.attrs["_FillValue"]
+        else:
+            logger.info(f"No _FillValue attribute found for field {self.name!r}, defaulting land_value to 0.0")
+            self._land_value = 0.0
+
         self.igrid = -1  # Default the grid index to -1
 
         if self.grid._mesh == "flat" or (self.name not in _unitconverters_map.keys()):
@@ -193,6 +202,15 @@ def interp_method(self, method: Callable):
         assert_same_function_signature(method, ref=ZeroInterpolator, context="Interpolation")
         self._interp_method = method
 
+    @property
+    def land_value(self):
+        """The value in the field data that represents land points. Can be used to mask land in Interpolators"""
+        return self._land_value
+
+    @land_value.setter
+    def land_value(self, value):
+        self._land_value = value
+
     def _check_velocitysampling(self):
         if self.name in ["U", "V", "W"]:
             warnings.warn(
diff --git a/src/parcels/_datasets/structured/circulation_models.py b/src/parcels/_datasets/structured/circulation_models.py
index 7933db367..d154a0d1b 100644
--- a/src/parcels/_datasets/structured/circulation_models.py
+++ b/src/parcels/_datasets/structured/circulation_models.py
@@ -24,6 +24,7 @@ def _copernicusmarine():
                     "long_name": "Eastward velocity",
                     "standard_name": "eastward_sea_water_velocity",
                     "valid_min": -5.0,
+                    "_FillValue": -999.99,
                 },
             ),
             "vo": (
@@ -36,6 +37,7 @@ def _copernicusmarine():
                     "long_name": "Northward velocity",
                     "standard_name": "northward_sea_water_velocity",
                     "valid_min": -5.0,
+                    "_FillValue": -999.99,
                 },
             ),
         },
@@ -103,6 +105,7 @@ def _copernicusmarine_waves():
                     "cell_methods": "time:point area:mean",
                     "missing_value": -32767,
                     "type_of_analysis": "spectral analysis",
+                    "_FillValue": -999.99,
                 },
             ),
             "VSDY": (
@@ -116,6 +119,7 @@ def _copernicusmarine_waves():
                     "cell_methods": "time:point area:mean",
                     "missing_value": -32767,
                     "type_of_analysis": "spectral analysis",
+                    "_FillValue": -999.99,
                 },
             ),
         },
diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py
index 8e23d4952..731b9a803 100644
--- a/src/parcels/interpolators.py
+++ b/src/parcels/interpolators.py
@@ -425,7 +425,9 @@ def _Spatialslip(
     def is_land(ti: int, zi: int, yi: int, xi: int):
         uval = corner_dataU[ti, zi, yi, xi, :]
         vval = corner_dataV[ti, zi, yi, xi, :]
-        return np.where(np.isclose(uval, 0.0) & np.isclose(vval, 0.0), True, False)
+        return np.where(
+            np.isclose(uval, vectorfield.U._land_value) & np.isclose(vval, vectorfield.V._land_value), True, False
+        )
 
     f_u = np.ones_like(xsi)
     f_v = np.ones_like(eta)
@@ -604,7 +606,7 @@ def XLinearInvdistLandTracer(
 
     corner_data = _get_corner_data_Agrid(field.data, ti, zi, yi, xi, lenT, lenZ, len(xsi), axis_dim)
 
-    land_mask = np.isclose(corner_data, 0.0)
+    land_mask = np.isclose(corner_data, field._land_value)
     nb_land = np.sum(land_mask, axis=(0, 1, 2, 3))
 
     if np.any(nb_land):
diff --git a/tests/test_field.py b/tests/test_field.py
index 7795c953d..b10e1951b 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -159,6 +159,20 @@ def invalid_interpolator_wrong_signature(particle_positions, grid_positions, inv
         )
 
 
+@pytest.mark.parametrize("fill_value", [-999.99, 0.0, 42, None])
+def test_field_land_value(fill_value):
+    ds = datasets_structured["ds_2d_left"].copy()
+    if fill_value is not None:
+        ds["data_g"].attrs["_FillValue"] = fill_value
+    grid = XGrid.from_dataset(ds, mesh="flat")
+
+    field = Field(name="test_field", data=ds["data_g"], grid=grid, interp_method=XLinear)
+    if fill_value is None:
+        assert field.land_value == 0.0
+    else:
+        assert field.land_value == fill_value
+
+
 def test_vectorfield_invalid_interpolator():
     ds = datasets_structured["ds_2d_left"]
     grid = XGrid.from_dataset(ds, mesh="flat")