diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst
index e64fad65f8..cb2747e520 100644
--- a/doc/recipe/preprocessor.rst
+++ b/doc/recipe/preprocessor.rst
@@ -2638,8 +2638,11 @@ This module contains the following preprocessor functions:
 
 * ``bias``: Calculate absolute or relative biases with respect to a reference
   dataset.
-* ``distance_metric``: Calculate absolute or relative biases with respect to a
-  reference dataset.
+* ``distance_metric``: Calculate distance metrics with respect to a reference
+  dataset.
+* ``t_test``: Calculate t-test (null hypothesis: 2 independent samples have
+  identical mean) with respect to a reference dataset and attach the
+  corresponding p-value as ancillary variable to the data.
 
 ``bias``
 --------
@@ -2937,8 +2940,117 @@ See also :func:`esmvalcore.preprocessor.distance_metric`.
 .. _Weighted Earth mover's distance: https://pythonot.github.io/
   quickstart.html#computing-wasserstein-distance
 
-.. _Other:
+.. _t_test:
+
+``t_test``
+----------
+
+This function performs a :func:`t-test <scipy.stats.ttest_ind>` (null
+hypothesis: 2 independent samples have identical mean) with respect to a
+reference dataset and attaches the corresponding p-value as ancillary variable
+to the data.
+
+With `Iris <https://scitools-iris.readthedocs.io/>`__, the p-value can be accessed with the following code:
+
+.. code-block:: python
+
+   import iris
+   cube = iris.load_cube("esmvaltool_output/recipe_t_test_20260721_095242/preproc/t-test/tas/CMIP6_BCC-ESM1_Amon_historical_r1i1p1f1_tas_gn_2000-2005.nc")
+   p_value = cube.ancillary_variable("p-value")
+   print(p_value)
+
+The returned p-value looks like this:
+
+.. code-block::
+
+  AncillaryVariable :  p-value / (1)
+      data: [
+          [0.21600565   , 0.18489039   , ..., 0.3673955    , 0.27788544   ],
+          [0.3407415    , 0.3905643    , ..., --           , 0.58700126   ],
+          ...,
+          [--           , --           , ..., --           , --           ],
+          [--           , --           , ..., --           , --           ]]
+      shape: (18, 36)
+      dtype: float32
+      long_name: 'p-value'
+      var_name: 'pvalue'
+
+For this preprocessor, exactly one input dataset needs to be declared as
+``reference_for_t_test: true`` in the recipe, e.g.,
+
+.. code-block:: yaml
+
+  datasets:
+    - {dataset: CanESM5, project: CMIP6, ensemble: r1i1p1f1, grid: gn}
+    - {dataset: CESM2,   project: CMIP6, ensemble: r1i1p1f1, grid: gn}
+    - {dataset: MIROC6,  project: CMIP6, ensemble: r1i1p1f1, grid: gn}
+    - {dataset: ERA-Interim, project: OBS6, tier: 3, type: reanaly, version: 1,
+       reference_for_t_test: true}
+
+In the example above, ERA-Interim is used as reference dataset for the t-test
+calculation.
+
+It is also possible to use the output from the :ref:`multi-model statistics` or
+:ref:`ensemble statistics` preprocessor as reference dataset.
+In this case, make sure to use ``reference_for_t_test: true`` for each dataset
+that will be used to create the reference dataset and use the option
+``keep_input_datasets: false`` for the multi-dataset preprocessor.
+For example:
+
+.. code-block:: yaml
+
+  datasets:
+    - {dataset: CanESM5, group: ref, reference_for_t_test: true}
+    - {dataset: CESM2,   group: ref, reference_for_t_test: true}
+    - {dataset: MIROC6,  group: notref}
+
+  preprocessors:
+    calculate_t_test:
+      custom_order: true
+      multi_model_statistics:
+        statistics: [mean]
+        span: overlap
+        groupby: [group]
+        keep_input_datasets: false
+      t_test:
+
+Here, the t-test for MIROC6 is calculated relative to the the multi-model mean
+from the models CanESM5 and CESM2.
+
+All datasets need to have the same shape.
+To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` might
+be helpful.
+
+The ``t_test`` preprocessor supports the following arguments in the recipe:
 
+* ``coords`` (:obj:`list` of :obj:`str`, default: ``None``): Coordinates over
+  which the t-test is calculated.
+  If ``None``, calculate the metric over all coordinates, which results in a
+  scalar cube.
+* Other parameters are passed to :func:`scipy.stats.ttest_ind`.
+
+Example:
+
+.. code-block:: yaml
+
+    preprocessors:
+      welchs_t_test:
+        t_test:
+          coords: [time]
+          equal_var: false
+          nan_policy: omit
+
+This will perform a `Welch's t-test
+<https://en.wikipedia.org/wiki/Welch%27s_t-test>`__, which does not assume
+equal variances between the samples.
+In addition, all masked values will be omitted in the t-test calculation (the
+default behavior is to return masked values if at least one input value is
+masked).
+
+See also :func:`esmvalcore.preprocessor.t_test`.
+
+
+.. _Other:
 
 Other
 =====
diff --git a/esmvalcore/_recipe/check.py b/esmvalcore/_recipe/check.py
index bb935cbc0a..436c30b311 100644
--- a/esmvalcore/_recipe/check.py
+++ b/esmvalcore/_recipe/check.py
@@ -627,6 +627,13 @@ def _check_ref_attributes(products: set, *, step: str, attr_name: str) -> None:
 )
 
 
+reference_for_t_test_preproc = partial(
+    _check_ref_attributes,
+    step="t_test",
+    attr_name="reference_for_t_test",
+)
+
+
 def statistics_preprocessors(settings: dict) -> None:
     """Check options of statistics preprocessors."""
     mm_stats = (
diff --git a/esmvalcore/_recipe/recipe.py b/esmvalcore/_recipe/recipe.py
index f15e61e92a..c974eda439 100644
--- a/esmvalcore/_recipe/recipe.py
+++ b/esmvalcore/_recipe/recipe.py
@@ -704,6 +704,7 @@ def _get_preprocessor_products(
 
     check.reference_for_bias_preproc(products)
     check.reference_for_distance_metric_preproc(products)
+    check.reference_for_t_test_preproc(products)
 
     _configure_multi_product_preprocessor(
         products=products,
diff --git a/esmvalcore/config/configurations/defaults/extra_facets_emac.yml b/esmvalcore/config/configurations/defaults/extra_facets_emac.yml
index 169a69860f..8ac3d9d05d 100644
--- a/esmvalcore/config/configurations/defaults/extra_facets_emac.yml
+++ b/esmvalcore/config/configurations/defaults/extra_facets_emac.yml
@@ -257,9 +257,12 @@ projects:
           ta:
             raw_name: [tm1_cav, tm1_ave, tm1]
             channel: Amon
+          # tro3:
+          #   raw_name: [O3_cav, O3_ave, O3]
+          #   channel: Amon
           tro3:
-            raw_name: [O3_cav, O3_ave, O3]
-            channel: Amon
+            raw_name: [O3_p39_cav, O3_cav, O3_ave, O3]
+            channel: AERmon
           ua:
             raw_name: [um1_cav, um1_ave, um1]
             channel: Amon
diff --git a/esmvalcore/preprocessor/__init__.py b/esmvalcore/preprocessor/__init__.py
index 3c17265530..7dac72fe4d 100644
--- a/esmvalcore/preprocessor/__init__.py
+++ b/esmvalcore/preprocessor/__init__.py
@@ -31,7 +31,11 @@
     meridional_statistics,
     zonal_statistics,
 )
-from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric
+from esmvalcore.preprocessor._compare_with_refs import (
+    bias,
+    distance_metric,
+    t_test,
+)
 from esmvalcore.preprocessor._concatenate import concatenate
 from esmvalcore.preprocessor._cycles import amplitude
 from esmvalcore.preprocessor._dask_progress import _compute_with_progress
@@ -218,6 +222,7 @@
     # Multi model statistics
     "multi_model_statistics",
     # Comparison with reference datasets
+    "t_test",
     "bias",
     "distance_metric",
     # Remove supplementary variables from cube
@@ -261,6 +266,7 @@
     "multi_model_statistics",
     "mask_multimodel",
     "mask_fillvalues",
+    "t_test",
 }
 
 
diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py
index 4e6fc6c180..b84350e81b 100644
--- a/esmvalcore/preprocessor/_compare_with_refs.py
+++ b/esmvalcore/preprocessor/_compare_with_refs.py
@@ -3,18 +3,20 @@
 from __future__ import annotations
 
 import logging
+import string
 from functools import partial
 from typing import TYPE_CHECKING, Any, Literal
 
 import dask
 import dask.array as da
+import dask.array.stats
 import iris.analysis
 import iris.analysis.stats
 import numpy as np
 from iris.common.metadata import CubeMetadata
-from iris.coords import CellMethod
+from iris.coords import AncillaryVariable, CellMethod
 from iris.cube import Cube, CubeList
-from scipy.stats import wasserstein_distance
+from scipy.stats import ttest_ind, wasserstein_distance
 
 from esmvalcore.iris_helpers import (
     ignore_iris_vague_metadata_warnings,
@@ -200,6 +202,16 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType) -> Cube:
     """Calculate bias for a single cube relative to a reference cube."""
     cube_metadata = cube.metadata
 
+    # Save ancillary variables and cell measures as they get removed in the
+    # bias calculation
+    ancillary_variables_and_dims = [
+        (a, cube.ancillary_variable_dims(a))
+        for a in cube.ancillary_variables()
+    ]
+    cell_measures_and_dims = [
+        (c, cube.cell_measure_dims(c)) for c in cube.cell_measures()
+    ]
+
     if bias_type == "absolute":
         cube = cube - reference
         new_units = cube.units
@@ -216,6 +228,12 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType) -> Cube:
     cube.metadata = cube_metadata
     cube.units = new_units
 
+    # Reattach ancillary variables and cell measures
+    for ancillary_variable, dims in ancillary_variables_and_dims:
+        cube.add_ancillary_variable(ancillary_variable, dims)
+    for cell_measure, dims in cell_measures_and_dims:
+        cube.add_cell_measure(cell_measure, dims)
+
     return cube
 
 
@@ -597,3 +615,200 @@ def _get_emd(
     if np.ma.is_masked(arr) or np.ma.is_masked(ref_arr):
         return np.ma.masked  # this is safe because PMFs will be masked arrays
     return wasserstein_distance(bin_centers, bin_centers, arr, ref_arr)
+
+
+def t_test(
+    products: set[PreprocessorFile] | Iterable[Cube],
+    reference: Cube | None = None,
+    coords: Iterable[Coord] | Iterable[str] | None = None,
+    **kwargs: Any,
+) -> set[PreprocessorFile] | CubeList:
+    """Perform t-test between dataset and reference dataset.
+
+    This is a test for the null hypothesis that 2 independent samples (dataset
+    and reference dataset) have identical average (expected) values.
+
+    All input datasets need to have the same shape. To ensure this, the
+    preprocessors :func:`esmvalcore.preprocessor.regrid` might be helpful.
+
+    Handles lazy data and masked data.
+
+    Notes
+    -----
+    The reference dataset can be specified with the ``reference`` argument. If
+    ``reference`` is ``None``, exactly one input dataset in the ``products``
+    set needs to have the facet ``reference_for_t_test: true`` defined in the
+    recipe. Please do **not** specify the option ``reference`` when using this
+    preprocessor function in a recipe.
+
+    Parameters
+    ----------
+    products:
+        Input datasets/cubes for which the t-test is performed relative to a
+        reference dataset/cube.
+    reference:
+        Cube which is used as reference for the t-test calculation. If
+        ``None``, `products` needs to be a :obj:`set` of
+        :class:`~esmvalcore.preprocessor.PreprocessorFile` objects and exactly
+        one dataset in ``products`` needs the facet ``reference_for_t_test:
+        true``. Do not specify this argument in a recipe.
+    coords:
+        Coordinates over which the t-test is calculated. If ``None``, calculate
+        the t-test over all coordinates, which results in a scalar cube.
+    **kwargs:
+        Additional keyword arguments passed to :func:`scipy.stats.ttest_ind`.
+
+    Returns
+    -------
+    :
+        Output datasets/cubes. Will be a :obj:`set` of
+        :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if
+        ``products`` is also one, a :class:`~iris.cube.CubeList` otherwise. The
+        cubes contain an attached :class:`~iris.coords.AncillaryVariable` which
+        is the p-value w.r.t. the chosen t-test. It can be accessed with
+        ``cube.ancillary_variable("p-value")``.
+
+    Raises
+    ------
+    ValueError
+        Not exactly one input datasets contains the facet
+        ``reference_for_t_test: true`` if ``reference=None``;
+        ``reference=None`` and the input products are given as iterable of
+        :class:`~iris.cube.Cube` objects.
+
+    """
+    ref_product = None
+    all_cubes_given = all(isinstance(p, Cube) for p in products)
+
+    # Get reference cube if not explicitly given
+    if reference is None:
+        if all_cubes_given:
+            msg = (
+                "A list of Cubes is given to this preprocessor; please "
+                "specify a `reference`"
+            )
+            raise ValueError(msg)
+        (reference, ref_product) = _get_ref(products, "reference_for_t_test")
+    else:
+        ref_product = None
+
+    # If input is an Iterable of Cube objects, calculate t-test for each element
+    if all_cubes_given:
+        cubes = [
+            _calculate_t_test(c, reference, coords=coords, **kwargs)
+            for c in products
+        ]
+        return CubeList(cubes)
+
+    # Otherwise, iterate over all input products, calculate t-test and adapt
+    # metadata and provenance information accordingly
+    output_products = set()
+    for product in products:
+        if product == ref_product:
+            continue
+        cube = concatenate(product.cubes)
+
+        # Calculate t-test
+        cube = _calculate_t_test(cube, reference, coords=coords, **kwargs)
+
+        # Adapt metadata and provenance information
+        if ref_product is not None:
+            product.wasderivedfrom(ref_product)
+
+        product.cubes = CubeList([cube])
+        output_products.add(product)
+
+    # Add reference dataset to output
+    if ref_product is not None:
+        output_products.add(ref_product)
+
+    return output_products
+
+
+def _calculate_t_test(
+    cube: Cube,
+    reference: Cube,
+    *,
+    coords: Iterable[Coord] | Iterable[str] | None = None,
+    **kwargs: Any,
+) -> Cube:
+    """Calculate the t-test and attach the p-value as ancillary variable to cube."""
+    cube = cube.copy()  # do not modify input cube
+
+    # Ensure that data is not chunked along desired coords
+    coords = get_all_coords(cube, coords)
+    cube = rechunk_cube(cube, coords)
+    reference = rechunk_cube(reference, coords)
+
+    axes = get_all_coord_dims(cube, coords)
+    n_axes = len(axes)
+
+    if cube.has_lazy_data() and reference.has_lazy_data():
+        # da.apply_gufunc transposes the input array so that the axes given by
+        # the `axes` argument to da.apply_gufunc are the rightmost dimensions.
+        # Thus, we need to use `along_axes=(ndim-n_axes, ..., ndim-2, ndim-1)`
+        # for _get_pvalue_from_ttest_ind here.
+        axes_in_chunk = tuple(range(cube.ndim - n_axes, cube.ndim))
+
+        # The call signature depends also on the number of axes in `axes`, and
+        # will be (a,b,...)->(nbins) where a,b,... are the data dimensions that
+        # are collapsed, and nbins the number of bin centers
+        in_signature = f"({','.join(list(string.ascii_lowercase)[:n_axes])})"
+        p_value_arr = da.apply_gufunc(
+            _get_pvalue_from_ttest_ind,
+            f"{in_signature},{in_signature}->()",
+            cube.lazy_data(),
+            reference.lazy_data(),
+            axes=[axes, axes, ()],
+            output_dtypes=cube.dtype,
+            along_axes=axes_in_chunk,
+            **kwargs,
+        )
+    else:
+        # Avoid realizing cube.data
+        cube_data = (
+            cube.lazy_data().compute() if cube.has_lazy_data() else cube.data
+        )
+        ref_data = (
+            reference.lazy_data().compute()
+            if reference.has_lazy_data()
+            else reference.data
+        )
+        p_value_arr = _get_pvalue_from_ttest_ind(
+            cube_data,
+            ref_data,
+            along_axes=axes,
+            **kwargs,
+        )
+
+    remaining_axes = tuple(sorted(set(range(cube.ndim)) - set(axes)))
+    cube.add_ancillary_variable(
+        AncillaryVariable(
+            p_value_arr,
+            long_name="p-value",
+            var_name="pvalue",
+            units="1",
+        ),
+        remaining_axes,
+    )
+
+    return cube
+
+
+@preserve_float_dtype
+def _get_pvalue_from_ttest_ind(
+    arr: np.ndarray,
+    ref_arr: np.ndarray,
+    *,
+    along_axes: int | Iterable[int] | None,
+    **kwargs: Any,
+) -> np.ndarray:
+    """Calculate :func:`scipy.stats.ttest_ind` and return p-value."""
+    # To support masked arrays, first convert masked values to NaNs;
+    # afterwards, convert NaNs back to masked values
+    arr = np.ma.filled(arr, np.nan)
+    ref_arr = np.ma.filled(ref_arr, np.nan)
+
+    t_test_result = ttest_ind(arr, ref_arr, axis=along_axes, **kwargs)
+
+    return np.ma.masked_invalid(t_test_result.pvalue)
diff --git a/tests/integration/recipe/test_recipe.py b/tests/integration/recipe/test_recipe.py
index 8dbe034e61..96dae3fe9a 100644
--- a/tests/integration/recipe/test_recipe.py
+++ b/tests/integration/recipe/test_recipe.py
@@ -3998,3 +3998,71 @@ def test_align_metadata_missing_arg(tmp_path, patched_datafinder, session):
     msg = "Missing required argument"
     with pytest.raises(ValueError, match=msg):
         get_recipe(tmp_path, content, session)
+
+
+def test_t_test_no_ref(tmp_path, patched_datafinder, session):
+    content = dedent("""
+        preprocessors:
+          test_t_test:
+            t_test:
+
+        diagnostics:
+          diagnostic_name:
+            variables:
+              ta:
+                preprocessor: test_t_test
+                project: CMIP6
+                mip: Amon
+                exp: historical
+                timerange: 20000101/20001231
+                ensemble: r1i1p1f1
+                grid: gn
+                additional_datasets:
+                  - {dataset: CanESM5}
+                  - {dataset: CESM2}
+
+            scripts: null
+        """)
+    msg = (
+        "Expected exactly 1 dataset with 'reference_for_t_test: true' in "
+        "products"
+    )
+    with pytest.raises(RecipeError) as exc:
+        get_recipe(tmp_path, content, session)
+    assert str(exc.value) == INITIALIZATION_ERROR_MSG
+    assert msg in exc.value.failed_tasks[0].message
+    assert "found 0" in exc.value.failed_tasks[0].message
+
+
+def test_t_test_two_refs(tmp_path, patched_datafinder, session):
+    content = dedent("""
+        preprocessors:
+          test_t_test:
+            t_test:
+
+        diagnostics:
+          diagnostic_name:
+            variables:
+              ta:
+                preprocessor: test_t_test
+                project: CMIP6
+                mip: Amon
+                exp: historical
+                timerange: 20000101/20001231
+                ensemble: r1i1p1f1
+                grid: gn
+                additional_datasets:
+                  - {dataset: CanESM5, reference_for_t_test: true}
+                  - {dataset: CESM2, reference_for_t_test: true}
+
+            scripts: null
+        """)
+    msg = (
+        "Expected exactly 1 dataset with 'reference_for_t_test: true' in "
+        "products"
+    )
+    with pytest.raises(RecipeError) as exc:
+        get_recipe(tmp_path, content, session)
+    assert str(exc.value) == INITIALIZATION_ERROR_MSG
+    assert msg in exc.value.failed_tasks[0].message
+    assert "found 2" in exc.value.failed_tasks[0].message
diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py
index d9de9c7d4b..c5ee399633 100644
--- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py
+++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py
@@ -1,17 +1,25 @@
 """Unit tests for :mod:`esmvalcore.preprocessor._compare_with_refs`."""
 
+from __future__ import annotations
+
 import contextlib
+import re
+from typing import Any
 
 import dask.array as da
 import iris
 import numpy as np
 import pytest
 from cf_units import Unit
-from iris.coords import CellMeasure, CellMethod
+from iris.coords import AncillaryVariable, AuxCoord, CellMeasure, CellMethod
 from iris.cube import Cube, CubeList
 from iris.exceptions import CoordinateNotFoundError
 
-from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric
+from esmvalcore.preprocessor._compare_with_refs import (
+    bias,
+    distance_metric,
+    t_test,
+)
 from tests import PreprocessorFile
 
 
@@ -343,6 +351,41 @@ def test_bias_products_and_ref_cube(
     assert product_a.mock_ancestors == set()
 
 
+@pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS)
+def test_bias_cubes_preserve_metadata(
+    regular_cubes,
+    ref_cubes,
+    bias_type,
+    data,
+    units,
+):
+    """Test calculation of bias with cubes."""
+    cube = regular_cubes[0]
+    ancillary_var = AncillaryVariable([0, 1], var_name="a")
+    cell_measure = CellMeasure([0, 1], var_name="c")
+    cube.add_ancillary_variable(ancillary_var, 0)
+    cube.add_cell_measure(cell_measure, 1)
+    ref_cube = ref_cubes[0]
+
+    out_cubes = bias(regular_cubes, ref_cube, bias_type=bias_type)
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == units
+    assert out_cube.dim_coords == cube.dim_coords
+    assert out_cube.aux_coords == cube.aux_coords
+    assert out_cube.ancillary_variables() == [ancillary_var]
+    assert out_cube.ancillary_variable_dims(ancillary_var) == (0,)
+    assert out_cube.cell_measures() == [cell_measure]
+    assert out_cube.cell_measure_dims(cell_measure) == (1,)
+
+
 def test_no_reference_for_bias(regular_cubes, ref_cubes):
     """Test fail when no reference_for_bias is given."""
     products = {
@@ -382,8 +425,8 @@ def test_invalid_bias_type(regular_cubes, ref_cubes):
         bias(products, bias_type="invalid_bias_type")
 
 
-def test_reference_none_cubes(regular_cubes):
-    """Test reference=None with with cubes."""
+def test_bias_reference_none_cubes(regular_cubes):
+    """Test bias with reference=None and cubes given."""
     msg = (
         "A list of Cubes is given to this preprocessor; please specify a "
         "`reference`"
@@ -932,3 +975,381 @@ def test_distance_metric_no_lon_for_area_weights(regular_cubes, metric, error):
             reference=ref_cube,
             coords=["time", "latitude"],
         )
+
+
+def assert_correct_p_value(
+    cube: Cube,
+    data: Any,  # noqa: ANN401
+    *,
+    dims: tuple[int, ...],
+    is_lazy: bool = False,
+) -> None:
+    assert len(cube.ancillary_variables()) == 1
+    p_value = cube.ancillary_variables()[0]
+    assert p_value.standard_name is None
+    assert p_value.long_name == "p-value"
+    assert p_value.var_name == "pvalue"
+    assert p_value.units == "1"
+    assert p_value.attributes == {}
+    assert p_value.has_lazy_data() is is_lazy
+    assert p_value.dtype == np.float32
+    print("p-value:", p_value.data)
+    if np.ma.is_masked(data):
+        np.testing.assert_equal(p_value.data.mask, data.mask)
+    np.testing.assert_allclose(p_value.data, data)
+    assert cube.ancillary_variable_dims(p_value) == dims
+
+
+@pytest.mark.parametrize("use_reference_product", [True, False])
+def test_t_test_products(regular_cubes, ref_cubes, use_reference_product):  # noqa: PLR0915
+    products = {
+        PreprocessorFile(regular_cubes, "A", {"dataset": "a"}),
+        PreprocessorFile(regular_cubes, "B", {"dataset": "b"}),
+    }
+
+    if use_reference_product:
+        ref_product = PreprocessorFile(
+            ref_cubes,
+            "REF",
+            {"reference_for_t_test": True},
+        )
+        products.add(ref_product)
+        expected_ancestors = {ref_product}
+        reference = None
+    else:
+        expected_ancestors = set()
+        reference = ref_cubes[0]
+
+    out_products = t_test(products, reference=reference, coords=["time"])
+
+    assert isinstance(out_products, set)
+    out_dict = products_set_to_dict(out_products)
+    assert len(out_dict) == len(products)
+
+    product_a = out_dict["A"]
+    assert product_a.filename == "A"
+    assert product_a.attributes == {"dataset": "a"}
+    assert len(product_a.cubes) == 1
+    out_cube = product_a.cubes[0]
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, regular_cubes[0].data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(
+        out_cube,
+        [[1.0, 0.6666667], [0.42264974, 0.46547753]],
+        dims=(1, 2),
+    )
+    assert product_a.wasderivedfrom.call_count == len(expected_ancestors)
+    assert product_a.mock_ancestors == expected_ancestors
+
+    product_b = out_dict["B"]
+    assert product_b.filename == "B"
+    assert product_b.attributes == {"dataset": "b"}
+    assert len(product_b.cubes) == 1
+    out_cube = product_b.cubes[0]
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, regular_cubes[0].data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(
+        out_cube,
+        [[1.0, 0.6666667], [0.42264974, 0.46547753]],
+        dims=(1, 2),
+    )
+    assert product_b.wasderivedfrom.call_count == len(expected_ancestors)
+    assert product_b.mock_ancestors == expected_ancestors
+
+    if use_reference_product:
+        product_ref = out_dict["REF"]
+        assert product_ref.filename == "REF"
+        assert product_ref.attributes == {"reference_for_t_test": True}
+        assert len(product_ref.cubes) == 1
+        out_cube = product_ref.cubes[0]
+        assert out_cube.dtype == np.float32
+        assert_allclose(out_cube.data, ref_cubes[0].data)
+        assert out_cube.var_name == "tas"
+        assert out_cube.standard_name == "air_temperature"
+        assert out_cube.units == "K"
+        assert out_cube.dim_coords == ref_cubes[0].dim_coords
+        assert out_cube.aux_coords == ref_cubes[0].aux_coords
+        assert out_cube.ancillary_variables() == []
+        product_ref.wasderivedfrom.assert_not_called()
+        assert product_ref.mock_ancestors == set()
+
+
+@pytest.mark.parametrize("coords_as_str", [True, False])
+@pytest.mark.parametrize("lazy", [True, False])
+def test_t_test_cubes_1d(regular_cubes, ref_cubes, lazy, coords_as_str):
+    cube = regular_cubes[0]
+    ref_cube = ref_cubes[0]
+    if lazy:
+        cube.data = cube.lazy_data().rechunk((1, 1, 1))
+        ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1))
+    coords = ["time"] if coords_as_str else [cube.coord("time")]
+
+    out_cubes = t_test([cube], ref_cube, coords=coords)
+
+    assert cube.has_lazy_data() is lazy
+    assert ref_cube.has_lazy_data() is lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(
+        out_cube,
+        [[1.0, 0.6666667], [0.42264974, 0.46547753]],
+        dims=(1, 2),
+        is_lazy=lazy,
+    )
+
+
+@pytest.mark.parametrize(
+    ("dims", "p_value", "out_dims"),
+    [
+        ((0, 1), [0.4679941236972809, 0.3202061057090759], (2,)),
+        ((0, 2), [0.6890520453453064, 0.17230847477912903], (1,)),
+    ],
+)
+@pytest.mark.parametrize("lazy", [True, False])
+def test_t_test_cubes_2d(
+    regular_cubes,
+    ref_cubes,
+    lazy,
+    dims,
+    p_value,
+    out_dims,
+):
+    cube = regular_cubes[0]
+    ref_cube = ref_cubes[0]
+    x_coord = AuxCoord([[0, 0], [0, 0]], var_name="x")
+    cube.add_aux_coord(x_coord, dims)
+    ref_cube.add_aux_coord(x_coord, dims)
+    if lazy:
+        cube.data = cube.lazy_data().rechunk((1, 1, 1))
+        ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1))
+
+    out_cubes = t_test([cube], reference=ref_cube, coords=["x"])
+
+    assert cube.has_lazy_data() is lazy
+    assert ref_cube.has_lazy_data() is lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(out_cube, p_value, dims=out_dims, is_lazy=lazy)
+
+
+@pytest.mark.parametrize(
+    ("coords", "p_value", "out_dims"),
+    [
+        (("time", "latitude"), [0.4679941236972809, 0.3202061057090759], (2,)),
+        (
+            ("time", "longitude"),
+            [0.6890520453453064, 0.17230847477912903],
+            (1,),
+        ),
+    ],
+)
+@pytest.mark.parametrize("lazy", [True, False])
+def test_t_test_cubes_multiple_coords(
+    regular_cubes,
+    ref_cubes,
+    lazy,
+    coords,
+    p_value,
+    out_dims,
+):
+    cube = regular_cubes[0]
+    ref_cube = ref_cubes[0]
+    if lazy:
+        cube.data = cube.lazy_data()
+        ref_cube.data = ref_cube.lazy_data()
+
+    out_cubes = t_test([cube], reference=ref_cube, coords=coords)
+
+    assert cube.has_lazy_data() is lazy
+    assert ref_cube.has_lazy_data() is lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(out_cube, p_value, dims=out_dims, is_lazy=lazy)
+
+
+@pytest.mark.parametrize("lazy", [True, False])
+def test_t_test_cubes_scalar_result(regular_cubes, ref_cubes, lazy):
+    cube = regular_cubes[0]
+    ref_cube = ref_cubes[0]
+    if lazy:
+        cube.data = cube.lazy_data().rechunk((1, 1, 1))
+        ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1))
+
+    out_cubes = t_test(
+        [cube],
+        ref_cube,
+        equal_var=False,
+    )
+
+    assert cube.has_lazy_data() is lazy
+    assert ref_cube.has_lazy_data() is lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == cube.dim_coords
+    assert out_cube.aux_coords == cube.aux_coords
+    assert_correct_p_value(
+        out_cube,
+        [0.20223067700862885],
+        dims=(),
+        is_lazy=lazy,
+    )
+
+
+@pytest.mark.parametrize("ref_lazy", [True, False])
+def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy):
+    cube = regular_cubes[0]
+    ref_cube = ref_cubes[0]
+    if ref_lazy:
+        ref_cube.data = ref_cube.lazy_data()
+    else:
+        cube.data = cube.lazy_data()
+
+    out_cubes = t_test([cube], ref_cube, coords=["time"])
+
+    assert cube.has_lazy_data() is not ref_lazy
+    assert ref_cube.has_lazy_data() is ref_lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is not ref_lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(
+        out_cube,
+        [[1.0, 0.6666667], [0.42264974, 0.46547753]],
+        dims=(1, 2),
+    )
+
+
+@pytest.mark.parametrize("lazy", [True, False])
+def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy):
+    cube = regular_cubes[0]
+    cube.data = np.ma.masked_inside(cube.data, 1.5, 3.5)
+    cube.data = np.ma.masked_greater(cube.data, 6.5)
+    ref_cube = ref_cubes[0]
+    ref_cube.data = np.ma.masked_invalid(
+        np.array(
+            [[[2.0, np.nan], [2.0, 2.0]], [[2.0, 2.0], [2.0, 2.0]]],
+            dtype=np.float32,
+        ),
+    )
+    if lazy:
+        cube.data = cube.lazy_data()
+        ref_cube.data = ref_cube.lazy_data()
+
+    out_cubes = t_test([cube], ref_cube, coords=["time"])
+
+    assert cube.has_lazy_data() is lazy
+    assert ref_cube.has_lazy_data() is lazy
+
+    assert isinstance(out_cubes, CubeList)
+    assert len(out_cubes) == 1
+    out_cube = out_cubes[0]
+
+    assert out_cube.has_lazy_data() is lazy
+    assert out_cube.dtype == np.float32
+    assert_allclose(out_cube.data, cube.data)
+    assert out_cube.var_name == "tas"
+    assert out_cube.standard_name == "air_temperature"
+    assert out_cube.units == "K"
+    assert out_cube.dim_coords == regular_cubes[0].dim_coords
+    assert out_cube.aux_coords == regular_cubes[0].aux_coords
+    assert_correct_p_value(
+        out_cube,
+        np.ma.masked_invalid([[1.0, np.nan], [np.nan, np.nan]]),
+        dims=(1, 2),
+        is_lazy=lazy,
+    )
+
+
+def test_t_test_reference_none_cubes_fail(regular_cubes):
+    """Test t_test with reference=None and cubes given."""
+    msg = (
+        r"A list of Cubes is given to this preprocessor; please specify a "
+        r"`reference`"
+    )
+    with pytest.raises(ValueError, match=re.escape(msg)):
+        t_test(regular_cubes)
+
+
+def test_no_reference_for_t_test_fail(regular_cubes, ref_cubes):
+    """Test fail when no reference_for_t_test is given."""
+    products = {
+        PreprocessorFile(regular_cubes, "A", {}),
+        PreprocessorFile(regular_cubes, "B", {}),
+        PreprocessorFile(ref_cubes, "REF", {}),
+    }
+    msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 0"
+    with pytest.raises(ValueError, match=re.escape(msg)):
+        t_test(products)
+
+
+def test_two_references_for_t_test_fail(regular_cubes, ref_cubes):
+    """Test fail when two reference_for_t_test products are given."""
+    products = {
+        PreprocessorFile(regular_cubes, "A", {"reference_for_t_test": False}),
+        PreprocessorFile(ref_cubes, "REF1", {"reference_for_t_test": True}),
+        PreprocessorFile(ref_cubes, "REF2", {"reference_for_t_test": True}),
+    }
+    msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 2"
+    with pytest.raises(ValueError, match=re.escape(msg)):
+        t_test(products)