diff --git a/esmvalcore/preprocessor/_regrid.py b/esmvalcore/preprocessor/_regrid.py
index b826f7bd3c..67e3e48c90 100644
--- a/esmvalcore/preprocessor/_regrid.py
+++ b/esmvalcore/preprocessor/_regrid.py
@@ -986,6 +986,24 @@ def regrid(
     # Rechunk and actually perform the regridding
     cube = _rechunk(cube, target_grid_cube)
     result = regridder(cube)
+    for ancillary_var in cube.ancillary_variables():
+        ancillary_dims = cube.ancillary_variable_dims(ancillary_var)
+        ancillary_slice = tuple(
+            slice(None) if i in ancillary_dims else 0 for i in range(cube.ndim)
+        )
+        ancillary_cube = cube[ancillary_slice].copy(ancillary_var.core_data())
+        ancillary_cube = _rechunk(ancillary_cube, target_grid_cube)
+        ancillary_result = regridder(ancillary_cube)
+        if result.has_lazy_data() and ancillary_result.has_lazy_data():
+            # Keep the chunks of the ancillary variable aligned with the
+            # regridded data variable.
+            ancillary_result.data = ancillary_result.lazy_data().rechunk(
+                result[ancillary_slice].lazy_data().chunks,
+            )
+        result.add_ancillary_variable(
+            ancillary_var.copy(ancillary_result.core_data()),
+            ancillary_dims,
+        )
     # Iris only supports regridding of 1D coordinates and iris-esmf-regrid
     # uses only the DimCoords if both grid_latitude/grid_longitude or
     # projection_x_coordinate/projection_y_coordinate DimCoords and latitude
diff --git a/tests/integration/preprocessor/_regrid/test_regrid.py b/tests/integration/preprocessor/_regrid/test_regrid.py
index e72c8b0383..0d946a782f 100644
--- a/tests/integration/preprocessor/_regrid/test_regrid.py
+++ b/tests/integration/preprocessor/_regrid/test_regrid.py
@@ -1,6 +1,11 @@
 """Integration tests for :func:`esmvalcore.preprocessor.regrid`."""
 
+import dask.array as da
 import iris
+import iris.coord_systems
+import iris.coords
+import iris.cube
+import iris.fileformats.pp
 import numpy as np
 import pytest
 from numpy import ma
@@ -364,6 +369,26 @@ def test_regrid__linear_with_mask(self, cache_weights):
         expected[:, 1, 1] = np.array([1.5, 5.5, 9.5])
         assert_array_equal(result.data, expected)
 
+    def test_regrid__linear_with_ancillary(self) -> None:
+        """Test that ancillary coordinates are also regridded."""
+        cube = self.cube.copy()
+        cube.data = cube.lazy_data()
+        cube.add_ancillary_variable(
+            iris.coords.AncillaryVariable(
+                da.arange(2, 6).astype(np.float32).reshape(2, 2),
+                var_name="ancillary",
+            ),
+            (1, 2),
+        )
+        result = regrid(cube, self.grid_for_linear, "linear")
+        ancillary_result = result.ancillary_variable("ancillary")
+        assert isinstance(ancillary_result, iris.coords.AncillaryVariable)
+        assert ancillary_result.has_lazy_data()
+        assert_array_equal(
+            ancillary_result.data,
+            np.array([3.5], dtype=np.float32).reshape(1, 1),
+        )
+
     @pytest.mark.parametrize("cache_weights", [True, False])
     def test_regrid__nearest(self, cache_weights):
         data = np.empty((1, 1))