diff --git a/esmvalcore/preprocessor/_time.py b/esmvalcore/preprocessor/_time.py
index 3338f0b519..217fdd443f 100644
--- a/esmvalcore/preprocessor/_time.py
+++ b/esmvalcore/preprocessor/_time.py
@@ -3,18 +3,19 @@
 Allows for selecting data subsets using certain time bounds;
 constructing seasonal and area averages.
 """
+import copy
 import datetime
 import logging
 from warnings import filterwarnings
 
 import dask.array as da
 import iris
-import iris.cube
 import iris.coord_categorisation
+import iris.cube
 import iris.exceptions
 import iris.util
-from iris.time import PartialDateTime
 import numpy as np
+from iris.time import PartialDateTime
 
 from ._shared import get_iris_analysis_operation, operator_accept_weights
 
@@ -460,7 +461,9 @@ def anomalies(cube, period, reference=None, standardize=False):
         reference_cube = extract_time(cube, **reference)
     reference = climate_statistics(reference_cube, period=period)
     if period in ['full']:
+        metadata = copy.deepcopy(cube.metadata)
         cube = cube - reference
+        cube.metadata = metadata
         if standardize:
             cube_stddev = climate_statistics(
                 cube, operator='std_dev', period=period)
diff --git a/tests/unit/preprocessor/_time/test_time.py b/tests/unit/preprocessor/_time/test_time.py
index 159eed1463..f9fd445bc8 100644
--- a/tests/unit/preprocessor/_time/test_time.py
+++ b/tests/unit/preprocessor/_time/test_time.py
@@ -1,25 +1,26 @@
 """Unit tests for the :func:`esmvalcore.preprocessor._time` module."""
 
+import copy
 import unittest
-import pytest
-import tests
-
-import numpy as np
-from numpy.testing import assert_array_equal, assert_array_almost_equal
 
-from cf_units import Unit
 import iris
 import iris.coord_categorisation
 import iris.coords
+import numpy as np
+import pytest
+from cf_units import Unit
 from iris.cube import Cube
+from numpy.testing import assert_array_almost_equal, assert_array_equal
 
-from esmvalcore.preprocessor._time import (
-    extract_month, extract_season, extract_time,
-    regrid_time,
-    decadal_statistics, annual_statistics, seasonal_statistics,
-    monthly_statistics, daily_statistics, timeseries_filter,
-    climate_statistics, anomalies
-)
+import tests
+from esmvalcore.preprocessor._time import (annual_statistics, anomalies,
+                                           climate_statistics,
+                                           daily_statistics,
+                                           decadal_statistics, extract_month,
+                                           extract_season, extract_time,
+                                           monthly_statistics, regrid_time,
+                                           seasonal_statistics,
+                                           timeseries_filter)
 
 
 def _create_sample_cube():
@@ -1022,6 +1023,20 @@ def test_standardized_anomalies(period, standardize=True):
             )
 
 
+@pytest.mark.parametrize('period, reference', PARAMETERS)
+def test_anomalies_preserve_metadata(period, reference, standardize=False):
+    cube = make_map_data(number_years=2)
+    cube.var_name = "si"
+    cube.units = "m"
+    metadata = copy.deepcopy(cube.metadata)
+    result = anomalies(cube, period, reference, standardize=standardize)
+    assert result.metadata == metadata
+    for coord_cube, coord_res in zip(cube.coords(), result.coords()):
+        if coord_cube.has_bounds() and coord_res.has_bounds():
+            assert_array_equal(coord_cube.bounds, coord_res.bounds)
+        assert coord_cube == coord_res
+
+
 @pytest.mark.parametrize('period, reference', PARAMETERS)
 def test_anomalies(period, reference, standardize=False):
     cube = make_map_data(number_years=2)