diff --git a/esmvalcore/_recipe.py b/esmvalcore/_recipe.py
index f362b74b65..35d429a336 100644
--- a/esmvalcore/_recipe.py
+++ b/esmvalcore/_recipe.py
@@ -754,11 +754,16 @@ def _update_multiproduct(input_products, order, preproc_dir, step):
         common_attributes = _get_common_attributes(products)
 
         for statistic in settings.get('statistics'):
-            common_attributes[step] = _get_tag(step, identifier, statistic)
-            filename = get_multiproduct_filename(common_attributes,
+            statistic_attributes = dict(common_attributes)
+            statistic_attributes[step] = _get_tag(step, identifier, statistic)
+            statistic_attributes.setdefault('alias',
+                                            statistic_attributes[step])
+            statistic_attributes.setdefault('dataset',
+                                            statistic_attributes[step])
+            filename = get_multiproduct_filename(statistic_attributes,
                                                  preproc_dir)
-            common_attributes['filename'] = filename
-            statistic_product = PreprocessorFile(common_attributes,
+            statistic_attributes['filename'] = filename
+            statistic_product = PreprocessorFile(statistic_attributes,
                                                  downstream_settings)
             output_products.add(statistic_product)
             relevant_settings['output_products'][identifier][
diff --git a/tests/__init__.py b/tests/__init__.py
index 41d3e3f830..03f1dff082 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -73,12 +73,17 @@ def assert_array_equal(self, a, b):
 class PreprocessorFile(mock.Mock):
     """Mocked PreprocessorFile."""
 
-    def __init__(self, cubes, filename, attributes, **kwargs):
+    def __init__(self, cubes, filename, attributes, settings=None, **kwargs):
         """Initialize with cubes."""
         super().__init__(spec=PreprocessorFileBase, **kwargs)
         self.cubes = cubes
         self.filename = filename
         self.attributes = attributes
-        self.settings = {}
+        if settings is None:
+            self.settings = {}
+        else:
+            self.settings = settings
         self.mock_ancestors = set()
         self.wasderivedfrom = mock.Mock(side_effect=self.mock_ancestors.add)
+
+    group = PreprocessorFileBase.group
diff --git a/tests/unit/test_recipe.py b/tests/unit/test_recipe.py
index b48a731958..844cdcdae5 100644
--- a/tests/unit/test_recipe.py
+++ b/tests/unit/test_recipe.py
@@ -356,6 +356,127 @@ def test_multi_model_filename():
     assert attributes['timerange'] == '1989/1992'
 
 
+def test_update_multiproduct_multi_model_statistics():
+    """Test ``_update_multiproduct``."""
+    settings = {'multi_model_statistics': {'statistics': ['mean', 'std_dev']}}
+    common_attributes = {
+        'project': 'CMIP6',
+        'timerange': '2000/2000',
+        'diagnostic': 'd',
+        'variable_group': 'var',
+    }
+    cube = iris.cube.Cube(np.array([1]))
+    products = [
+        PreprocessorFile(cube, 'A',
+                         attributes={'dataset': 'a', **common_attributes},
+                         settings=settings),
+        PreprocessorFile(cube, 'B',
+                         attributes={'dataset': 'b', **common_attributes},
+                         settings=settings),
+        PreprocessorFile(cube, 'C',
+                         attributes={'dataset': 'c', **common_attributes},
+                         settings=settings),
+        PreprocessorFile(cube, 'D',
+                         attributes={'dataset': 'd', **common_attributes},
+                         settings=settings),
+    ]
+    order = ('load', 'multi_model_statistics', 'save')
+    preproc_dir = '/preproc'
+    step = 'multi_model_statistics'
+    output, settings = _recipe._update_multiproduct(
+        products, order, preproc_dir, step)
+
+    assert len(output) == 2
+
+    filenames = [p.filename for p in output]
+    assert '/preproc/d/var/CMIP6_MultiModelMean_2000-2000.nc' in filenames
+    assert '/preproc/d/var/CMIP6_MultiModelStd_Dev_2000-2000.nc' in filenames
+
+    for product in output:
+        for attr in common_attributes:
+            assert attr in product.attributes
+            assert product.attributes[attr] == common_attributes[attr]
+            assert 'alias' in product.attributes
+            assert 'dataset' in product.attributes
+            assert 'multi_model_statistics' in product.attributes
+        if 'MultiModelStd_Dev' in product.filename:
+            assert product.attributes['alias'] == 'MultiModelStd_Dev'
+            assert product.attributes['dataset'] == 'MultiModelStd_Dev'
+            assert (product.attributes['multi_model_statistics'] ==
+                    'MultiModelStd_Dev')
+        elif 'MultiModelMean' in product.filename:
+            assert product.attributes['alias'] == 'MultiModelMean'
+            assert product.attributes['dataset'] == 'MultiModelMean'
+            assert (product.attributes['multi_model_statistics'] ==
+                    'MultiModelMean')
+
+    assert len(settings) == 1
+    output_products = settings['output_products']
+    assert len(output_products) == 1
+    stats = output_products['']
+    assert len(stats) == 2
+    assert 'mean' in stats
+    assert 'std_dev' in stats
+    assert 'MultiModelMean' in stats['mean'].filename
+    assert 'MultiModelStd_Dev' in stats['std_dev'].filename
+
+
+def test_update_multiproduct_ensemble_statistics():
+    """Test ``_update_multiproduct``."""
+    settings = {'ensemble_statistics': {'statistics': ['median']}}
+    common_attributes = {
+        'dataset': 'CanESM2',
+        'project': 'CMIP6',
+        'timerange': '2000/2000',
+        'diagnostic': 'd',
+        'variable_group': 'var',
+    }
+    cube = iris.cube.Cube(np.array([1]))
+    products = [
+        PreprocessorFile(cube, 'A',
+                         attributes=common_attributes,
+                         settings=settings),
+        PreprocessorFile(cube, 'B',
+                         attributes=common_attributes,
+                         settings=settings),
+        PreprocessorFile(cube, 'C',
+                         attributes=common_attributes,
+                         settings=settings),
+        PreprocessorFile(cube, 'D',
+                         attributes=common_attributes,
+                         settings=settings),
+    ]
+    order = ('load', 'ensemble_statistics', 'save')
+    preproc_dir = '/preproc'
+    step = 'ensemble_statistics'
+    output, settings = _recipe._update_multiproduct(
+        products, order, preproc_dir, step)
+
+    assert len(output) == 1
+    product = list(output)[0]
+    assert (product.filename ==
+            '/preproc/d/var/CMIP6_CanESM2_EnsembleMedian_2000-2000.nc')
+
+    for attr in common_attributes:
+        assert attr in product.attributes
+        assert product.attributes[attr] == common_attributes[attr]
+        assert 'alias' in product.attributes
+        assert product.attributes['alias'] == 'EnsembleMedian'
+        assert 'dataset' in product.attributes
+        assert product.attributes['dataset'] == 'CanESM2'
+        assert 'ensemble_statistics' in product.attributes
+        assert product.attributes['ensemble_statistics'] == 'EnsembleMedian'
+
+    assert len(settings) == 1
+    output_products = settings['output_products']
+    assert len(output_products) == 1
+    stats = output_products['CMIP6_CanESM2']
+    assert len(stats) == 1
+    assert 'median' in stats
+    assert (stats['median'].filename ==
+            '/preproc/d/var/CMIP6_CanESM2_EnsembleMedian_2000-2000.nc')
+
+
 def test_update_multiproduct_no_product():
     cube = iris.cube.Cube(np.array([1]))
     products = [