diff --git a/src/CSET/operators/__init__.py b/src/CSET/operators/__init__.py
index 01ccb8e57..bf55c0f66 100644
--- a/src/CSET/operators/__init__.py
+++ b/src/CSET/operators/__init__.py
@@ -82,6 +82,7 @@
     "transect",
     "wind",
     "write",
+    "time_utils",
 ]
 
 logger = logging.getLogger(__name__)
diff --git a/src/CSET/operators/aggregate.py b/src/CSET/operators/aggregate.py
index 451bb0eed..498b60242 100644
--- a/src/CSET/operators/aggregate.py
+++ b/src/CSET/operators/aggregate.py
@@ -15,6 +15,7 @@
 """Operators to aggregate across either 1 or 2 dimensions."""
 
 import logging
+from typing import NamedTuple
 
 import iris
 import iris.analysis
@@ -24,6 +25,8 @@
 import iris.util
 import isodate
 import numpy as np
+from iris.coords import AuxCoord, DimCoord
+from iris.cube import Cube, CubeList
 
 from CSET._common import iter_maybe
 from CSET.operators._utils import is_time_aggregatable
@@ -43,16 +46,16 @@ def _add_nref(cube: iris.cube.Cube):
 
 
 def time_aggregate(
-    cube: iris.cube.Cube,
+    cubes: iris.cube.Cube | iris.cube.CubeList,
     method: str,
     interval_iso: str,
     **kwargs,
-) -> iris.cube.Cube:
+) -> iris.cube.Cube | iris.cube.CubeList:
     """Aggregate cube by its time coordinate.
 
     Aggregates similar (stash) fields in a cube for the specified coordinate and
     using the method supplied. The aggregated cube will keep the coordinate and
-    add a further coordinate with the aggregated end time points.
+    add a coordinate with the aggregated end time points.
 
     Examples are: 1. Generating hourly or 6-hourly precipitation accumulations
     given an interval for the new time coordinate.
@@ -66,8 +69,8 @@ def time_aggregate(
 
     Arguments
     ---------
-    cube: iris.cube.Cube
-        Cube to aggregate and iterate over one dimension
+    cubes: iris.cube.Cube | iris.cube.CubeList
+        Cube or CubeList to aggregate and iterate over one dimension
     coordinate: str
         Coordinate to aggregate over i.e. 'time', 'longitude',
         'latitude','model_level_number'.
@@ -90,19 +93,33 @@ def time_aggregate(
     # Duration of ISO timedelta.
     timedelta = isodate.parse_duration(interval_iso)
 
+    if timedelta == "0":
+        return cubes
+
+    resampled_cubes = iris.cube.CubeList()
+
     # Convert interval format to whole hours.
     interval = int(timedelta.total_seconds() / 3600)
 
-    # Add time categorisation overwriting hourly increment via lambda coord.
-    # https://scitools-iris.readthedocs.io/en/latest/_modules/iris/coord_categorisation.html
-    iris.coord_categorisation.add_categorised_coord(
-        cube, "interval", "time", lambda coord, cell: cell // interval * interval
-    )
+    cubes = iter_maybe(cubes)
+
+    for cube in cubes:
+        # Add time categorisation overwriting hourly increment via lambda coord.
+        # https://scitools-iris.readthedocs.io/en/latest/_modules/iris/coord_categorisation.html
+        iris.coord_categorisation.add_categorised_coord(
+            cube, "interval", "time", lambda coord, cell: cell // interval * interval
+        )
 
-    # Aggregate cube using supplied method.
-    aggregated_cube = cube.aggregated_by("interval", getattr(iris.analysis, method))
-    aggregated_cube.remove_coord("interval")
-    return aggregated_cube
+        # Aggregate cube using supplied method.
+        aggregated_cube = cube.aggregated_by("interval", getattr(iris.analysis, method))
+        aggregated_cube.remove_coord("interval")
+
+        resampled_cubes.append(aggregated_cube)
+
+    if len(resampled_cubes) == 1:
+        return resampled_cubes[0]
+    else:
+        return resampled_cubes
 
 
 def ensure_aggregatable_across_cases(
@@ -208,6 +225,42 @@ def get_buckets(self) -> list[iris.cube.CubeList]:
     return aggregatable_cubes
 
 
+def combine_obs_across_forecasts(cubes: CubeList) -> Cube:
+    """
+    Combine observation cubes from multiple forecast_reference_times.
+
+    Input:
+        CubeList of cubes with dimensions
+
+            (time, station)
+
+    Output:
+        Cube with dimensions
+
+            (forecast_reference_time,
+             forecast_period,
+             station)
+
+    where
+
+        time
+
+    becomes a 2D auxiliary coordinate attached to
+
+        (forecast_reference_time, forecast_period)
+
+    Only stations present in every forecast are retained.
+    All station metadata coordinates are preserved.
+    """
+    if len(cubes) < 2:
+        raise ValueError("Need at least two cubes")
+
+    common_stations = _get_common_stations(cubes)
+    station_lookup = _build_station_lookup(cubes, common_stations)
+    fp_hours = _generate_forecast_period(cubes)
+    return _make_aggregated_obs_cube(cubes, station_lookup, common_stations, fp_hours)
+
+
 def add_hour_coordinate(
     cubes: iris.cube.Cube | iris.cube.CubeList,
 ) -> iris.cube.Cube | iris.cube.CubeList:
@@ -284,3 +337,221 @@ def rolling_window_time_aggregation(
         return new_cubelist[0]
     else:
         return new_cubelist
+
+
+def _get_common_stations(cubes: CubeList) -> list[str]:
+    # --------------------------------------------------------------
+    # Find stations common to all cubes with complete data
+    # --------------------------------------------------------------
+    if len(cubes) < 2:
+        raise ValueError(
+            f"Need at least two cubes to find common stations, but got {len(cubes)}"
+        )
+
+    valid_station_sets = []
+
+    for cb in cubes:
+        names = cb.coord("Station_Name").points
+        data = cb.data
+
+        # Handle masked and unmasked arrays
+        if np.ma.isMaskedArray(data):
+            mask = np.ma.getmaskarray(data)
+            station_valid = ~np.any(mask, axis=0) & np.all(
+                np.isfinite(data.filled(np.nan)), axis=0
+            )
+        else:
+            station_valid = np.all(np.isfinite(data), axis=0)
+        valid_station_sets.append(set(names[station_valid]))
+
+    common_stations = sorted(set.intersection(*valid_station_sets))
+
+    logger.info(
+        "Retaining %s stations with complete observations", len(common_stations)
+    )
+
+    if not common_stations:
+        raise ValueError(
+            "No stations with complete data in all forecast_reference_times"
+        )
+    return common_stations
+
+
+class StationLookup(NamedTuple):
+    """Station lookup data structure."""
+
+    subset_data: list[np.ma.MaskedArray | np.ndarray]
+    frt_points: list[float]
+    time_points: list[np.ndarray]
+
+
+def _build_station_lookup(cubes: CubeList, common_stations: list[str]) -> StationLookup:
+    # --------------------------------------------------------------
+    # Build station lookup for every cube
+    # --------------------------------------------------------------
+
+    subset_data = []
+    frt_points = []
+    time_points = []
+
+    for cube in cubes:
+        names = cube.coord("Station_Name").points
+        lookup = {name: idx for idx, name in enumerate(names)}
+        station_indices = [lookup[name] for name in common_stations]
+        subcube = cube[:, station_indices]
+        subset_data.append(subcube.data)
+        frt_points.append(cube.coord("forecast_reference_time").points[0])
+        time_points.append(cube.coord("time").points)
+
+    # --------------------------------------------------------------
+    # Check all cubes have same time axis length
+    # --------------------------------------------------------------
+
+    ntime = len(time_points[0])
+
+    for t in time_points[1:]:
+        if len(t) != ntime:
+            raise ValueError("Forecasts have different numbers of lead times")
+
+    return StationLookup(subset_data, frt_points, time_points)
+
+
+def _generate_forecast_period(cubes) -> np.ndarray:
+    # --------------------------------------------------------------
+    # Generate forecast period
+    # --------------------------------------------------------------
+
+    time_coord = cubes[0].coord("time")
+    frt_coord = cubes[0].coord("forecast_reference_time")
+
+    frt_date = frt_coord.units.num2date(frt_coord.points[0])
+
+    fp_hours = []
+    for dt in time_coord.units.num2date(time_coord.points):
+        fp_hours.append((dt - frt_date).total_seconds() / 3600)
+
+    fp_hours = np.asarray(fp_hours)
+
+    return fp_hours
+
+
+def _make_aggregated_obs_cube(
+    cubes: CubeList,
+    station_lookup: StationLookup,
+    common_stations: list[str],
+    fp_hours: np.ndarray,
+) -> Cube:
+    time_coord = cubes[0].coord("time")
+
+    # --------------------------------------------------------------
+    # Stack data
+    # --------------------------------------------------------------
+
+    data = np.stack(station_lookup.subset_data, axis=0)
+
+    # shape:
+    #
+    # (forecast_reference_time,
+    #  forecast_period,
+    #  station)
+
+    # --------------------------------------------------------------
+    # Output coordinates
+    # --------------------------------------------------------------
+
+    frt_out = DimCoord(
+        station_lookup.frt_points,
+        standard_name="forecast_reference_time",
+        units=cubes[0].coord("forecast_reference_time").units,
+    )
+
+    fp_out = DimCoord(
+        fp_hours,
+        standard_name="forecast_period",
+        units="hours",
+    )
+
+    station_out = DimCoord(
+        np.arange(len(common_stations)),
+        long_name="station",
+    )
+
+    cube_out = Cube(
+        data,
+        standard_name=cubes[0].standard_name,
+        long_name=cubes[0].long_name,
+        var_name=cubes[0].var_name,
+        units=cubes[0].units,
+        attributes=cubes[0].attributes.copy(),
+        dim_coords_and_dims=[
+            (frt_out, 0),
+            (fp_out, 1),
+            (station_out, 2),
+        ],
+    )
+
+    # --------------------------------------------------------------
+    # Preserve station metadata coordinates
+    # --------------------------------------------------------------
+
+    ref_cube = cubes[0]
+
+    ref_names = ref_cube.coord("Station_Name").points
+
+    ref_lookup = {name: idx for idx, name in enumerate(ref_names)}
+
+    common_idx = [ref_lookup[name] for name in common_stations]
+
+    # skip coords as we have awkward station and station_0 arbitrary monotonic arrays.
+    for coord in ref_cube.aux_coords:
+        dims = ref_cube.coord_dims(coord)
+
+        # only coords attached solely to station axis
+        if dims != (1,):
+            continue
+
+        values = coord.points[common_idx]
+
+        # verify same in every cube
+        for cube in cubes[1:]:
+            cube_names = cube.coord("Station_Name").points
+
+            cube_lookup = {name: idx for idx, name in enumerate(cube_names)}
+
+            idx = [cube_lookup[name] for name in common_stations]
+
+            other_values = cube.coord(coord.name()).points[idx]
+
+            if not np.array_equal(
+                values,
+                other_values,
+            ):
+                raise ValueError(f"Station metadata differs for coord '{coord.name()}'")
+
+        aux = AuxCoord(
+            values,
+            standard_name=coord.standard_name,
+            long_name=coord.long_name,
+            var_name=coord.var_name,
+            units=coord.units,
+            attributes=coord.attributes.copy(),
+        )
+
+        cube_out.add_aux_coord(aux, (2,))
+
+    # --------------------------------------------------------------
+    # Add valid-time auxiliary coord
+    # --------------------------------------------------------------
+
+    time_2d = np.vstack(station_lookup.time_points)
+
+    cube_out.add_aux_coord(
+        AuxCoord(
+            time_2d,
+            standard_name="time",
+            units=time_coord.units,
+        ),
+        (0, 1),
+    )
+
+    return cube_out
diff --git a/src/CSET/operators/imageprocessing.py b/src/CSET/operators/imageprocessing.py
index 6d767a231..06e9f550c 100644
--- a/src/CSET/operators/imageprocessing.py
+++ b/src/CSET/operators/imageprocessing.py
@@ -23,7 +23,7 @@
 
 from CSET._common import is_increasing
 from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
-from CSET.operators.misc import _extract_common_time_points
+from CSET.operators.time_utils import _extract_common_time_points
 from CSET.operators.regrid import regrid_onto_cube
 
 logger = logging.getLogger(__name__)
diff --git a/src/CSET/operators/misc.py b/src/CSET/operators/misc.py
index aae2fcd58..ffd41e230 100644
--- a/src/CSET/operators/misc.py
+++ b/src/CSET/operators/misc.py
@@ -27,6 +27,7 @@
 from CSET._common import is_increasing, iter_maybe
 from CSET.operators._utils import fully_equalise_attributes, get_cube_yxcoordname
 from CSET.operators.regrid import regrid_onto_cube
+from CSET.operators.time_utils import _extract_common_time_points
 
 logger = logging.getLogger(__name__)
 
@@ -443,8 +444,10 @@ def difference(cubes: CubeList):
         other.data = np.flip(other.data, other.coord(other_lat_name).cube_dims(other))
 
     # Extract just common time points.
+
     base, other = _extract_common_time_points(base, other)
 
+
     # Equalise attributes so we can merge.
     fully_equalise_attributes([base, other])
     logger.debug("Base: %s\nOther: %s", base, other)
@@ -467,55 +470,6 @@ def difference(cubes: CubeList):
     difference.data = other.data - base.data
     return difference
 
-
-def _extract_common_time_points(base: Cube, other: Cube) -> tuple[Cube, Cube]:
-    """Extract common time points from cubes to allow comparison."""
-    # Get the name of the first non-scalar time coordinate.
-    time_coord = next(
-        (
-            coord.name()
-            for coord in filter(
-                lambda coord: coord.shape > (1,) and coord.name() in ["time", "hour"],
-                base.coords(),
-            )
-        ),
-        None,
-    )
-    if not time_coord:
-        logger.debug("No time coord, skipping equalisation.")
-        return (base, other)
-    base_time_coord = base.coord(time_coord)
-    other_time_coord = other.coord(time_coord)
-    logger.debug("Base: %s\nOther: %s", base_time_coord, other_time_coord)
-    if time_coord == "hour":
-        # We directly compare points when comparing coordinates with
-        # non-absolute units, such as hour. We can't just check the units are
-        # equal as iris automatically converts to datetime objects in the
-        # comparison for certain coordinate names.
-        base_times = base_time_coord.points
-        other_times = other_time_coord.points
-        shared_times = set.intersection(set(base_times), set(other_times))
-    else:
-        # Units don't match, so converting to datetimes for comparison.
-        base_times = base_time_coord.units.num2date(base_time_coord.points)
-        other_times = other_time_coord.units.num2date(other_time_coord.points)
-        shared_times = set.intersection(set(base_times), set(other_times))
-    logger.debug("Shared times: %s", shared_times)
-    time_constraint = iris.Constraint(
-        coord_values={
-            time_coord: lambda cell, shared_times=shared_times: (
-                cell.point in shared_times
-            )
-        }
-    )
-    # Extract points matching the shared times.
-    base = base.extract(time_constraint)
-    other = other.extract(time_constraint)
-    if base is None or other is None:
-        raise ValueError("No common time points found!")
-    return (base, other)
-
-
 def convert_units(cubes: iris.cube.Cube | iris.cube.CubeList, units: str):
     """Convert the units of a cube.
 
@@ -638,6 +592,7 @@ def extract_common_points(cubes: iris.cube.CubeList, coordinate: str):
         CubeList containing the two cubes sliced to common points
         for the given coordinate.
     """
+
     # Check type of input
     if type(cubes) is not iris.cube.CubeList:
         raise TypeError(f"Not a CubeList, got type {type(cubes)}")
diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py
index 6138512ea..d7da4bdf9 100644
--- a/src/CSET/operators/plot.py
+++ b/src/CSET/operators/plot.py
@@ -64,7 +64,7 @@
     validate_cubes_coords,
 )
 from CSET.operators.collapse import collapse
-from CSET.operators.misc import _extract_common_time_points
+from CSET.operators.time_utils import _extract_common_time_points
 from CSET.operators.regrid import regrid_onto_cube
 
 logger = logging.getLogger(__name__)
diff --git a/src/CSET/operators/regrid.py b/src/CSET/operators/regrid.py
index 9edb68741..e63cf87e5 100644
--- a/src/CSET/operators/regrid.py
+++ b/src/CSET/operators/regrid.py
@@ -20,11 +20,13 @@
 import iris
 import iris.coord_systems
 import iris.cube
+from iris.cube import Cube, CubeList
 import numpy as np
 from iris.analysis.cartography import rotate_pole
 
 from CSET._common import iter_maybe
 from CSET.operators._utils import get_cube_yxcoordname
+from CSET.operators.time_utils import _extract_common_time_points, _extract_common_time_points_multiple_cubes
 
 logger = logging.getLogger(__name__)
 
@@ -407,7 +409,9 @@ def transform_lat_long_points(lon, lat, cube):
 
 
 def interpolate_to_point_cube(
-    fld: iris.cube.Cube | iris.cube.CubeList, point_cube: iris.cube.Cube, **kwargs
+    fld: iris.cube.Cube | iris.cube.CubeList,
+    point_cube: iris.cube.Cube | iris.cube.CubeList,
+    **kwargs,
 ) -> iris.cube.Cube | iris.cube.CubeList:
     """Interpolate a 2D field in cube or CubeList to a set of points.
 
@@ -432,34 +436,17 @@ def interpolate_to_point_cube(
     # Empty CubeList To store regridded cubes.
     regridded_cubes = iris.cube.CubeList()
 
+    # Generate array of point cube lat and lon points.
+    point_lat_name, point_lon_name = get_cube_yxcoordname(point_cube)
+    point_lats = point_cube.coord(point_lat_name).points
+    point_lons = point_cube.coord(point_lon_name).points
+
+    #extract common time points between fld cubes
+    fld = _extract_common_time_points_multiple_cubes(fld)
     # Iterate over all cubes and regrid.
     for cube in iter_maybe(fld):
-        # Ensure matching times in fld cube and point_cube
-        base_time_coord = point_cube.coord("time")
-        other_time_coord = cube.coord("time")
-        base_times = base_time_coord.units.num2date(base_time_coord.points)
-        other_times = other_time_coord.units.num2date(other_time_coord.points)
-        shared_times = set.intersection(set(base_times), set(other_times))
-        logger.debug("Shared times: %s", shared_times)
-        time_constraint = iris.Constraint(
-            coord_values={
-                "time": lambda cell, shared_times=shared_times: (
-                    cell.point in shared_times
-                )
-            }
-        )
-
-        # Extract points matching the shared times.
-        cube = cube.extract(time_constraint)
-        point_cube = point_cube.extract(time_constraint)
-        if cube is None or point_cube is None:
-            raise ValueError("No common time points found!")
-
-        # Generate array of point cube lat and lon points.
-        point_lat_name, point_lon_name = get_cube_yxcoordname(point_cube)
-        point_lats = point_cube.coord(point_lat_name).points
-        point_lons = point_cube.coord(point_lon_name).points
-
+        base, cube = _extract_common_time_points(point_cube, cube)
+        # Get forecast field cube spatial names.
         y_coord, x_coord = get_cube_yxcoordname(cube)
 
         # Rotate point_cube coords if required to match model coord rotation.
@@ -488,53 +475,119 @@ def interpolate_to_point_cube(
             ]
 
         # Interpolate fld cube to required sample points
+
         fld_point_cube = cube.interpolate(
             sample_points, iris.analysis.Linear(extrapolation_mode="mask")
         )
 
-        # Retain only diagonal elements of 2D interpolated cube to vector points
-        od_index = point_cube.coord_dims("station")[0]
+        diag_data = np.diagonal(
+            fld_point_cube.data,
+            axis1=-2,
+            axis2=-1,
+        )
+
         fv_cube = iris.cube.Cube(
-            fld_point_cube.data.diagonal(axis1=od_index, axis2=od_index + 1),
+            diag_data,
             standard_name=cube.standard_name,
             long_name=cube.long_name,
             units=cube.units,
         )
-        # Copy all non-lat/lon coordinates and cube attributes
-        if "time" in [coord.name() for coord in fld_point_cube.coords(dim_coords=True)]:
+
+        #
+        # Add all non-horizontal dimension coordinates
+        # from the forecast cube.
+        #
+        out_dim = 0
+        for coord in cube.coords(dim_coords=True):
+            if coord.name() in (
+                "latitude",
+                "longitude",
+                "grid_latitude",
+                "grid_longitude",
+            ):
+                continue
+
             fv_cube.add_dim_coord(
-                point_cube.coord("time"), point_cube.coord_dims("time")[0]
+                coord.copy(),
+                out_dim,
             )
+
+            out_dim += 1
+
+        #
+        # Station dimension is always last.
+        #
+        station_dim = out_dim
+
+        fv_cube.add_dim_coord(
+            point_cube.coord("station").copy(),
+            station_dim,
+        )
+
+        #
+        # Copy forecast auxiliary coordinates.
+        #
         for coord in cube.coords():
-            if coord.name() not in [
-                "time",
+            if coord.name() in (
                 "latitude",
                 "longitude",
                 "grid_latitude",
                 "grid_longitude",
-            ] and coord.name() not in [coord.name() for coord in fv_cube.coords()]:
-                fv_cube.add_aux_coord(coord.copy(), cube.coord_dims(coord))
+            ):
+                continue
+
+            if coord.name() in [c.name() for c in fv_cube.coords()]:
+                continue
+
+            dims = cube.coord_dims(coord)
+
+            if dims:
+                fv_cube.add_aux_coord(
+                    coord.copy(),
+                    dims,
+                )
+            else:
+                fv_cube.add_aux_coord(
+                    coord.copy(),
+                )
+
+        #
+        # Copy observation auxiliary coordinates.
+        #
         for coord in point_cube.coords():
-            if coord.name() not in [
-                "time",
+            if coord.name() in (
+                "station",
                 "forecast_period",
                 "forecast_reference_time",
-                "realization",
-                "station",
-            ] and coord.name() not in [coord.name() for coord in fv_cube.coords()]:
-                fv_cube.add_aux_coord(coord.copy(), point_cube.coord_dims(coord))
-        fv_cube.add_dim_coord(point_cube.coord("station"), od_index)
+            ):
+                continue
+
+            if coord.name() in [c.name() for c in fv_cube.coords()]:
+                continue
+
+            dims = point_cube.coord_dims(coord)
+
+            if dims:
+                fv_cube.add_aux_coord(
+                    coord.copy(),
+                    dims,
+                )
+            else:
+                fv_cube.add_aux_coord(
+                    coord.copy(),
+                )
+
         fv_cube.attributes = cube.attributes.copy()
         fv_cube.cell_methods = cube.cell_methods
         fv_cube.units = cube.units
+
         regridded_cubes.append(fv_cube)
 
-    # Preserve returning a cube if only a cube has been supplied to regrid.
+    # Preserve returning a cube if only a cube supplied.
     if len(regridded_cubes) == 1:
         return regridded_cubes[0]
-    else:
-        return regridded_cubes
 
+    return regridded_cubes
 
 def vertical_interpolation(
     cubes: iris.cube.Cube | iris.cube.CubeList,
@@ -576,3 +629,4 @@ def vertical_interpolation(
         return interpolated_cubes[0]
     else:
         return interpolated_cubes
+
diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py
index 922972ae0..41d83bfa4 100644
--- a/src/CSET/operators/scoreswrappers.py
+++ b/src/CSET/operators/scoreswrappers.py
@@ -35,7 +35,7 @@
     generate_realization_constraint,
     generate_remove_single_ensemble_member_constraint,
 )
-from CSET.operators.misc import _extract_common_time_points
+from CSET.operators.time_utils import _extract_common_time_points
 from CSET.operators.read import _realization_callback
 from CSET.operators.regrid import regrid_onto_cube
 
@@ -45,7 +45,7 @@
 def scores_rmse(
     cubes: CubeList,
     preserved_coordinates: list[str] | str | None = None,
-) -> CubeList:
+) -> CubeList | Cube:
     r"""Calculate the Root Mean Square Error (RMSE) using scores.
 
     Acts as a wrapper around the RMSE calculation from ``scores`` ([scoresa]_, [scoresb]_).
@@ -71,12 +71,9 @@ def scores_rmse(
         A cubelist containing the RMSE between the base and other cube.
     """
     scores_cubelist = CubeList()
-
     base, others = _split_base_and_other(cubes)
 
     for other in others:
-        base, other = _process_cubes_for_verification(base, other)
-
         scores_cube = _make_scores_cube(base, other, "rmse", preserved_coordinates)
 
         scores_cube.rename(f"RMSE_of_{base.name()}")
@@ -114,8 +111,6 @@ def scores_mae(
     base, others = _split_base_and_other(cubes)
 
     for other in others:
-        base, other = _process_cubes_for_verification(base, other)
-
         scores_cube = _make_scores_cube(base, other, "mae", preserved_coordinates)
 
         scores_cube.rename(f"MAE_of_{base.name()}")
@@ -155,8 +150,6 @@ def scores_additive_bias(
     base, others = _split_base_and_other(cubes)
 
     for other in others:
-        base, other = _process_cubes_for_verification(base, other)
-
         scores_cube = _make_scores_cube(
             base, other, "additive_bias", preserved_coordinates
         )
@@ -195,8 +188,6 @@ def scores_correlation_pearsonr(
     base, others = _split_base_and_other(cubes)
 
     for other in others:
-        base, other = _process_cubes_for_verification(base, other)
-
         scores_cube = _make_scores_cube(
             base, other, "pearson_correlation", preserved_coordinates
         )
@@ -670,12 +661,16 @@ def _make_scores_cube(
 
 
     """
+    if not base.long_name or not "observed" in base.long_name:
+        base, other = _process_cubes_for_verification(base, other)
+
     other_xr = xr.DataArray.from_iris(other)
     base_xr = xr.DataArray.from_iris(base)
     preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates)
 
     # Scores operates on xarray data arrays, so we transform the iris cube into an array,
     # apply scores, and then transform it back.
+
     if metric == "rmse":
         scores_cube = xr.DataArray.to_iris(
             scores.continuous.rmse(other_xr, base_xr, preserve_dims=preserve_dims)
@@ -805,6 +800,7 @@ def _process_cubes_for_verification(base: Cube, other: Cube) -> tuple[Cube, Cube
     # on variable type. Linear regridding can in general be appropriate for smooth
     # variables. Care should be taken with interpretation of differences
     # given this dependency on regridding.
+
     if (
         base.coord(base_lat_name).shape != other.coord(other_lat_name).shape
         or base.coord(base_lon_name).shape != other.coord(other_lon_name).shape
@@ -931,7 +927,7 @@ def _attach_scaler_time_coord_maybe(scores_cube: Cube, base: Cube) -> None:
 
     """
     try:
-        if not scores_cube.coords("time"):
+        if not scores_cube.coords("time") and not scores_cube.coords("forecast_period"):
             base_time = base.coord("time")
             time_vals = (
                 base_time.bounds.flatten()
@@ -953,10 +949,15 @@ def _attach_scaler_time_coord_maybe(scores_cube: Cube, base: Cube) -> None:
                     attributes=base_time.attributes.copy(),
                 )
             )
+
     except iris.exceptions.CoordinateNotFoundError:
         pass
 
 
+def _get_obs_cube(cubes: CubeList):
+    return [cb for cb in cubes if "observed" in (cb.long_name or "")]
+
+
 def _split_base_and_other(cubes: CubeList):
     r"""Split the cube into base and other cubes.
 
@@ -976,7 +977,7 @@ def _split_base_and_other(cubes: CubeList):
         A tuple containing a base cube, and other cube/cubelist.
 
     """
-    obs_cube = [cb for cb in cubes if "observed" in (cb.long_name or "")]
+    obs_cube = _get_obs_cube(cubes)
     if obs_cube:
         if len(obs_cube) > 1:
             raise ValueError(
diff --git a/src/CSET/operators/time_utils.py b/src/CSET/operators/time_utils.py
new file mode 100644
index 000000000..9d1063a3b
--- /dev/null
+++ b/src/CSET/operators/time_utils.py
@@ -0,0 +1,213 @@
+
+import logging
+
+import iris
+import iris.analysis.calculus
+import numpy as np
+from iris.cube import Cube, CubeList
+
+logger = logging.getLogger(__name__)
+
+
+def _extract_common_time_points_multiple_cubes(cubes:  CubeList | Cube) -> CubeList | Cube:
+    """Equalise time points across all cubes."""
+
+    if isinstance(cubes,Cube) or len(cubes) < 2:
+        return cubes
+    if cubes[0].coords("forecast_reference_time"):
+        return _extract_common_time_points_multiple_cubes_with_frt(cubes)
+    else:
+        return _extract_common_time_points_multiple_cubes_no_frt(cubes)
+
+def _extract_common_time_points(base: Cube, other: Cube) -> tuple[Cube, Cube]:
+    """Extract common time points from cubes to allow comparison."""
+    # Get the name of the first non-scalar time coordinate.
+
+    if base.coords("forecast_reference_time") and np.size(base.coords("forecast_reference_time"))>1:
+        return _extract_common_time_points_with_frt(base, other)
+    else:
+        return _extract_common_time_points_no_frt(base, other)
+
+
+def _extract_common_time_points_multiple_cubes_with_frt(cubes: CubeList) -> CubeList:
+
+
+    # Check all cubes have identical FRTs.
+    reference_frts = cubes[0].coord("forecast_reference_time").points
+
+    for cube in cubes[1:]:
+        cube_frts = cube.coord("forecast_reference_time").points
+
+        if not np.array_equal(reference_frts, cube_frts):
+            raise ValueError(
+                "Cubes do not share the same "
+                "forecast_reference_time values."
+            )
+
+    # Start from the first cube's forecast periods.
+    shared_periods = set(
+        cubes[0].coord("forecast_period").points
+    )
+
+    # Find intersection across all cubes.
+    for cube in cubes[1:]:
+        shared_periods &= set(
+            cube.coord("forecast_period").points
+        )
+
+    if not shared_periods:
+        raise ValueError("No common forecast periods found.")
+
+    logger.debug(
+        "Common forecast periods: %s",
+        sorted(shared_periods),
+    )
+
+    constraint = iris.Constraint(
+        forecast_period=lambda cell, shared_periods=shared_periods: (
+                cell.point in shared_periods
+        )
+    )
+
+    output = CubeList()
+
+    for cube in cubes:
+        extracted = cube.extract(constraint)
+
+        if extracted is None:
+            raise ValueError(
+                f"No common forecast periods remain for {cube.name()}"
+            )
+
+        output.append(extracted)
+
+    return output
+
+def _extract_common_time_points_multiple_cubes_no_frt(cubes: CubeList) -> CubeList:
+    base = cubes[0]
+    time_coord = next(
+        (
+            coord.name()
+            for coord in base.coords()
+            if coord.shape > (1,) and coord.name() in ("time", "hour")
+        ),
+        None,
+    )
+    if not time_coord:
+        logger.debug("No time coord, skipping equalisation.")
+        return cubes
+
+    base_time_coord = base.coord(time_coord)
+    shared_times = set(_get_time_points(base_time_coord))
+    logger.debug("Shared times: %s", shared_times)
+
+    for other in cubes[1:]:
+        other_time_coord = other.coord(time_coord)
+        logger.debug("Base: %s\nOther: %s", base_time_coord, other_time_coord)
+        shared_times &= set(_get_time_points(other_time_coord))
+        logger.debug("Shared times: %s", shared_times)
+
+    if not shared_times:
+        raise ValueError("No common time points found!")
+
+    time_constraint = iris.Constraint(
+        coord_values={
+            time_coord: lambda cell, shared_times=shared_times: cell.point in shared_times
+        }
+    )
+    return cubes.extract(time_constraint)
+
+def _extract_common_time_points_no_frt(base: Cube, other: Cube) -> tuple[Cube, Cube]:
+    time_coord = next(
+        (
+            coord.name()
+            for coord in filter(
+            lambda coord: coord.shape > (1,) and coord.name() in ["time", "hour"],
+            base.coords(),
+        )
+        ),
+        None,
+    )
+    if not time_coord:
+        logger.debug("No time coord, skipping equalisation.")
+        return base, other
+    base_time_coord = base.coord(time_coord)
+    other_time_coord = other.coord(time_coord)
+    logger.debug("Base: %s\nOther: %s", base_time_coord, other_time_coord)
+
+    base_times = _get_time_points(base_time_coord)
+    other_times = _get_time_points(other_time_coord)
+
+    shared_times = set.intersection(set(base_times), set(other_times))
+
+    logger.debug("Shared times: %s", shared_times)
+    time_constraint = iris.Constraint(
+        coord_values={
+            time_coord: lambda cell, shared_times=shared_times: (
+                    cell.point in shared_times
+            )
+        }
+    )
+
+    # Extract points matching the shared times.
+    base = base.extract(time_constraint)
+    other = other.extract(time_constraint)
+    if base is None or other is None:
+        raise ValueError("No common time points found!")
+    return base, other
+
+def _get_time_points(coord):
+    """Comparable time points for a coord. 'hour' uses raw points since
+    it has non-absolute units; iris auto-converts other time coords to
+    datetimes for comparison, so match that here."""
+    if coord.name() == "hour":
+        return coord.points
+    return coord.units.num2date(coord.points)
+
+
+def _extract_common_time_points_with_frt(base: Cube,other: Cube) -> tuple[Cube, Cube]:
+    """Equalise time points across all cubes."""
+
+    # Check all cubes have identical FRTs.
+    reference_frts = base.coord("forecast_reference_time").points
+    cube_frts = other.coord("forecast_reference_time").points
+
+    if not np.array_equal(reference_frts, cube_frts):
+        raise ValueError(
+            "Cubes do not share the same "
+            "forecast_reference_time values."
+        )
+
+    constraint = _make_shared_period_constraint(base,other)
+
+    base = base.extract(constraint)
+    other = other.extract(constraint)
+
+    return base,other
+
+
+def _make_shared_period_constraint(base, other) -> iris.Constraint:
+    # Start from the first cube's forecast periods.
+    shared_periods = set(
+        base.coord("forecast_period").points
+    )
+
+    # Find intersection across all cubes.
+
+    shared_periods &= set(
+        other.coord("forecast_period").points
+    )
+
+    if not shared_periods:
+        raise ValueError("No common forecast periods found.")
+
+    logger.debug(
+        "Common forecast periods: %s",
+        sorted(shared_periods),
+    )
+
+    return iris.Constraint(
+        forecast_period=lambda cell, shared_periods=shared_periods: (
+                cell.point in shared_periods
+        )
+    )
diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml
index 0c8c0634e..41bc4d1fe 100644
--- a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml
+++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE.yaml
@@ -2,7 +2,7 @@ category: Surface Model verus Observation Time Series
 title: "$VARNAME scores mae time series for model vs observation points $SUBAREA_NAME"
 description: Extracts and plot a time series of $VARNAME for all times based on RMSE between given models and observation points.
 
- The RMSE is calculated based on that used in the
+  The MAE is calculated based on that used in the
   package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mae).
 
   This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the MAE is
diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE_aggregation.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE_aggregation.yaml
new file mode 100644
index 000000000..0cde9ed2a
--- /dev/null
+++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_MAE_aggregation.yaml
@@ -0,0 +1,66 @@
+category: Surface Model versus Observation Time Series Aggregation
+title: "$VARNAME\nscores MAE time series for model vs observation points $SUBAREA_NAME aggregation by $AGGREGATION_METHOD"
+description: |
+  The MAE is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mae).
+
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the MAE is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+
+  A larger MAE implies a greater error than a smaller MAE. An
+  MAE of zero indicates the two fields most likely match. The MAE is calculated
+  on the grid point and thus a spatial view of the MAE provides useful
+  information about where the differences are, or if placement errors
+  are domininating the score (usually indicated by dipoles). MAE is a fair
+  measure of error when forecasting the median and is not sensitive to outliers.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - operator: misc.combine_cubes_into_cubelist
+    first:
+      operator: aggregate.combine_obs_across_forecasts
+      cubes:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: observed_$VARNAME
+    second:
+      operator: regrid.interpolate_to_point_cube
+      fld:
+        operator: aggregate.ensure_aggregatable_across_cases
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: $VARNAME
+
+      point_cube:
+        operator: aggregate.combine_obs_across_forecasts
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: observed_$VARNAME
+
+  - operator: scoreswrappers.scores_mae
+    preserved_coordinates: $AGGREGATION_METHOD
+
+  - operator: plot.plot_line_series
+    series_coordinate: $AGGREGATION_METHOD
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_RMSE_aggregation.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_RMSE_aggregation.yaml
new file mode 100644
index 000000000..1f88bacfa
--- /dev/null
+++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_RMSE_aggregation.yaml
@@ -0,0 +1,71 @@
+category: Surface Model versus Observation Time Series Aggregation
+title: "$VARNAME\nscores rmse time series for model vs observation points $SUBAREA_NAME aggregation by $AGGREGATION_METHOD"
+description: |
+  Extracts and plot a time series of $VARNAME for all times based on RMSE between given models and observation points.
+
+  The RMSE is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.rmse).
+
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the RMSE is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+
+  A larger RMSE implies a greater error than a smaller RMSE. An
+  RMSE of zero indicates the two fields match. The RMSE is calculated
+  on the grid point and thus a spatial view of the RMSE provides useful
+  information about where the differences are, or if placement errors
+  are domininating the score (usually indicated by dipoles).
+
+  For RMSE values computed across an entire casse study, the computation
+  is done within the scores RMSE module, rather than trying to mean across
+  each timestep. This preserves the nonlinear nature of the RMSE metric.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - operator: misc.combine_cubes_into_cubelist
+    first:
+      operator: aggregate.combine_obs_across_forecasts
+      cubes:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: observed_$VARNAME
+    second:
+      operator: regrid.interpolate_to_point_cube
+      fld:
+        operator: aggregate.ensure_aggregatable_across_cases
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: $VARNAME
+
+      point_cube:
+        operator: aggregate.combine_obs_across_forecasts
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: observed_$VARNAME
+
+  - operator: scoreswrappers.scores_rmse
+    preserved_coordinates: $AGGREGATION_METHOD
+
+  - operator: plot.plot_line_series
+    series_coordinate: $AGGREGATION_METHOD
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias_aggregation.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias_aggregation.yaml
new file mode 100644
index 000000000..020a953d4
--- /dev/null
+++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_additive_bias_aggregation.yaml
@@ -0,0 +1,64 @@
+category: Surface Model versus Observation Time Series Aggregation
+title: "$VARNAME\nscores additive bias time series for model vs observation points $SUBAREA_NAME aggregation by $AGGREGATION_METHOD"
+description: |
+  Extracts and plots the Mean Error (also known as the Additive Bias) in $VARNAME
+  computed over the domain for each timestep. The ME is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.mean_error).
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the ME is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+  A larger ME implies a greater error than a smaller ME. Although a near zero ME could be indicative
+  of a good match between two fields, this could also arise if there are large compensating positive and
+  negative  errors. The ME is calculated on the grid point and thus a spatial view of the ME provides useful
+  information about whether a field has a bias for over or under-forecasting.
+
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - operator: misc.combine_cubes_into_cubelist
+    first:
+      operator: aggregate.combine_obs_across_forecasts
+      cubes:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: observed_$VARNAME
+    second:
+      operator: regrid.interpolate_to_point_cube
+      fld:
+        operator: aggregate.ensure_aggregatable_across_cases
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: $VARNAME
+
+      point_cube:
+        operator: aggregate.combine_obs_across_forecasts
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: observed_$VARNAME
+
+  - operator: scoreswrappers.scores_additive_bias
+    preserved_coordinates: $AGGREGATION_METHOD
+
+  - operator: plot.plot_line_series
+    series_coordinate: $AGGREGATION_METHOD
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
diff --git a/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr_aggregation.yaml b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr_aggregation.yaml
new file mode 100644
index 000000000..87e61dbc5
--- /dev/null
+++ b/src/CSET/recipes/verification/timeseries_surface_scores_model_vs_obs_correlation_pearsonr_aggregation.yaml
@@ -0,0 +1,62 @@
+category: Surface Model versus Observation Time Series Aggregation
+title: "$VARNAME\nscores pearson correlation time series for model vs observation points $SUBAREA_NAME aggregation by $AGGREGATION_METHOD"
+description: |
+  Extracts and plots the Pearson's Correlation coefficient in $VARNAME
+  computed over the domain for each timestep. The PC is calculated based on that used in the
+  package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.correlation.pearsonr).
+  This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the PC is
+  collapsed over all other coordinates in the cube and calculated for every timestep.
+
+  The PC coefficient provides information on the linear relationship between two fields.
+  Perfectly correlated fields would yield a PC coefficient of 1, whereas perfectly anti-correlated
+  fields would yield a PC coefficient of -1. A value of 0 would indicate that the fields are non-correlated.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+    constraint:
+      operator: constraints.generate_var_constraint
+      varname: ['observed_$VARNAME', '$VARNAME']
+
+  - operator: filters.filter_multiple_cubes
+    constraint:
+        operator: constraints.generate_cell_methods_constraint
+        cell_methods: []
+
+  - operator: misc.combine_cubes_into_cubelist
+    first:
+      operator: aggregate.combine_obs_across_forecasts
+      cubes:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: observed_$VARNAME
+    second:
+      operator: regrid.interpolate_to_point_cube
+      fld:
+        operator: aggregate.ensure_aggregatable_across_cases
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: $VARNAME
+
+      point_cube:
+        operator: aggregate.combine_obs_across_forecasts
+        cubes:
+          operator: filters.filter_multiple_cubes
+          constraint:
+            operator: constraints.generate_var_constraint
+            varname: observed_$VARNAME
+
+  - operator: scoreswrappers.scores_correlation_pearsonr
+    preserved_coordinates: $AGGREGATION_METHOD
+
+  - operator: plot.plot_line_series
+    series_coordinate: $AGGREGATION_METHOD
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
diff --git a/tests/conftest.py b/tests/conftest.py
index bbd604b56..06618217b 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -25,6 +25,7 @@
 from uuid import uuid4
 
 import cf_units
+import cftime
 import iris
 import iris.coord_systems
 import iris.coords
@@ -1355,3 +1356,178 @@ def _make_cube(data, long_name, model_name=None):
         return cube
 
     return _make_cube
+
+
+def _make_test_cube_stations(
+    shape: tuple[int, int],
+    seed: int,
+    long_name: str,
+    standard_name: str | None = None,
+    model_name: str | None = None,
+):
+    rng = np.random.default_rng(seed)
+
+    fp = DimCoord(
+        np.arange(shape[0]),
+        standard_name="forecast_period",
+        units="hours",
+    )
+    station = DimCoord(
+        np.arange(shape[1]),
+        long_name="station",
+        units="no_unit",
+    )
+
+    station_name_coord = iris.coords.AuxCoord(
+        points=np.array([f"st{i}" for i in range(shape[1])]),
+        long_name="Station_Name",
+        units="unknown",
+    )
+
+    latitude_coord = iris.coords.AuxCoord(
+        points=50 + (0.5 * np.arange(shape[1])),
+        standard_name="latitude",
+        units="degrees",
+    )
+
+    longitude_coord = iris.coords.AuxCoord(
+        points=-5 + (0.2 * np.arange(shape[1])),
+        standard_name="longitude",
+        units="degrees",
+    )
+
+    data = rng.normal(loc=280, scale=5, size=shape)
+    return Cube(
+        data,
+        long_name=long_name,
+        standard_name=standard_name,
+        units="K",
+        dim_coords_and_dims=[(fp, 0), (station, 1)],
+        attributes={"model_name": model_name},
+        aux_coords_and_dims=[
+            (station_name_coord, 1),
+            (latitude_coord, 1),
+            (longitude_coord, 1),
+        ],
+    )
+
+
+@pytest.fixture
+def dummy_cubelist_obs_3_common_stations():
+    """CubeList of [obs_cube1, obs_cube2] with time (per-row) and scalar forecast reference time coords."""
+    time_units = cf_units.Unit("hours since 1970-01-01", calendar="360_day")
+
+    def add_time_coords(cube, time_start, frt_dt, hour_step=1):
+        """Add a per-row time coord (data_dims=0) and scalar frt coord to a cube, in place."""
+        n_times = cube.shape[0]
+        time_datetimes = [
+            time_start + datetime.timedelta(hours=hour_step * i) for i in range(n_times)
+        ]
+        cube.add_aux_coord(
+            AuxCoord(
+                points=time_units.date2num(time_datetimes),
+                standard_name="time",
+                units=time_units,
+            ),
+            data_dims=0,
+        )
+        cube.add_aux_coord(
+            AuxCoord(
+                points=time_units.date2num(frt_dt),
+                standard_name="forecast_reference_time",
+                units=time_units,
+            )
+        )
+        return cube
+
+    obs_cube1 = add_time_coords(
+        _make_test_cube_stations(
+            shape=(10, 3), seed=1, long_name="observed_temperature_at_screen_level"
+        ),
+        time_start=cftime.datetime(2024, 1, 1, 6, 0, calendar="360_day"),
+        frt_dt=cftime.datetime(2024, 1, 1, 0, 0, calendar="360_day"),
+    )
+    obs_cube2 = add_time_coords(
+        _make_test_cube_stations(
+            shape=(10, 5), seed=1, long_name="observed_temperature_at_screen_level"
+        ),
+        time_start=cftime.datetime(2024, 1, 2, 6, 0, calendar="360_day"),
+        frt_dt=cftime.datetime(2024, 1, 2, 0, 0, calendar="360_day"),
+    )
+
+    return CubeList([obs_cube1, obs_cube2])
+
+
+def _make_test_cube_multi_forecasts(
+    shape: tuple[int, int, int],
+    seed: int,
+    long_name: str,
+    standard_name: str | None = None,
+    model_name: str | None = None,
+):
+    rng = np.random.default_rng(seed)
+
+    frt = DimCoord(
+        np.arange(shape[0]),
+        standard_name="forecast_reference_time",
+        units="hours since 1970-01-01",
+    )
+    fp = DimCoord(
+        np.arange(shape[1]),
+        standard_name="forecast_period",
+        units="hours",
+    )
+    station = DimCoord(
+        np.arange(shape[2]),
+        long_name="station",
+        units="no_unit",
+    )
+
+    station_name_coord = iris.coords.AuxCoord(
+        points=np.array([f"st{i}" for i in range(shape[2])]),
+        long_name="Station_Name",
+        units="unknown",
+    )
+
+    latitude_coord = iris.coords.AuxCoord(
+        points=rng.uniform(50.0, 55.0, size=shape[2]),
+        standard_name="latitude",
+        units="degrees",
+    )
+
+    longitude_coord = iris.coords.AuxCoord(
+        points=rng.uniform(-5.0, -2.0, size=shape[2]),
+        standard_name="longitude",
+        units="degrees",
+    )
+
+    data = rng.normal(loc=280, scale=5, size=shape)
+
+    return Cube(
+        data,
+        long_name=long_name,
+        standard_name=standard_name,
+        units="K",
+        dim_coords_and_dims=[(frt, 0), (fp, 1), (station, 2)],
+        attributes={"model_name": model_name},
+        aux_coords_and_dims=[
+            (station_name_coord, 2),
+            (latitude_coord, 2),
+            (longitude_coord, 2),
+        ],
+    )
+
+
+@pytest.fixture
+def dummy_cubelist_model_obs_multiple_forecasts():
+    """CubeList of [obs_cube, model_cube] with forecast reference time."""
+    obs_cube = _make_test_cube_multi_forecasts(
+        shape=(2, 10, 12), seed=1, long_name="observed_temperature_at_screen_level"
+    )
+    model_cube = _make_test_cube_multi_forecasts(
+        shape=(2, 10, 12),
+        seed=2,
+        long_name="temperature_at_screen_level",
+        model_name="model_a",
+    )
+    return CubeList([obs_cube, model_cube])
diff --git a/tests/operators/test_aggregate.py b/tests/operators/test_aggregate.py
index ce3e97f7f..eb0008dfc 100644
--- a/tests/operators/test_aggregate.py
+++ b/tests/operators/test_aggregate.py
@@ -18,6 +18,7 @@
 import iris.cube
 import numpy as np
 import pytest
+from iris.cube import CubeList
 
 from CSET.operators import aggregate
 
@@ -159,3 +160,76 @@ def test_rolling_window_time_aggregation_cubelist(long_forecast):
         assert cube_a.shape == cube_b.shape
         assert cube_a.shape[0] == cube_c.shape[0] - 23
         assert np.allclose(cube_a.data, cube_b.data, rtol=1e-6, atol=1e-2)
+
+
+def test_get_common_stations(dummy_cubelist_obs_3_common_stations):
+    """Test _get_common_stations basic functionality."""
+    common_stations = aggregate._get_common_stations(
+        dummy_cubelist_obs_3_common_stations
+    )
+    assert len(common_stations) == 3
+    assert common_stations == ["st0", "st1", "st2"]
+
+
+def test_get_common_stations_one_cube_fail(dummy_cubelist_obs_3_common_stations):
+    """Test _get_common_stations basic functionality."""
+    with pytest.raises(
+        ValueError, match="Need at least two cubes to find common stations, but got 1"
+    ):
+        aggregate._get_common_stations(
+            CubeList([dummy_cubelist_obs_3_common_stations[0]])
+        )
+
+
+def test_build_station_lookup(dummy_cubelist_obs_3_common_stations):
+    """Test _build_station_lookup basic functionality."""
+    common_stations = aggregate._get_common_stations(
+        dummy_cubelist_obs_3_common_stations
+    )
+    station_lookup = aggregate._build_station_lookup(
+        dummy_cubelist_obs_3_common_stations, common_stations
+    )
+
+    assert np.shape(station_lookup.subset_data) == (2, 10, 3)
+
+    expected_frt_points = [
+        cb.coord("forecast_reference_time").points.item()
+        for cb in dummy_cubelist_obs_3_common_stations
+    ]
+    assert station_lookup.frt_points == expected_frt_points
+
+    expected_time_points = [
+        cb.coord("time").points for cb in dummy_cubelist_obs_3_common_stations
+    ]
+    np.testing.assert_array_equal(station_lookup.time_points, expected_time_points)
+
+
+def test_generate_forecast_period(dummy_cubelist_obs_3_common_stations):
+    """Test _generate_forecast_period basic functionality."""
+    forecast_period = aggregate._generate_forecast_period(
+        dummy_cubelist_obs_3_common_stations
+    )
+
+    expected = np.array([6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0])
+    np.testing.assert_allclose(forecast_period, expected)
+
+
+def test_make_aggregated_obs_cube(dummy_cubelist_obs_3_common_stations):
+    """Test _make_aggregated_obs_cube basic functionality."""
+    common_stations = aggregate._get_common_stations(
+        dummy_cubelist_obs_3_common_stations
+    )
+    station_lookup = aggregate._build_station_lookup(
+        dummy_cubelist_obs_3_common_stations, common_stations
+    )
+    forecast_period = aggregate._generate_forecast_period(
+        dummy_cubelist_obs_3_common_stations
+    )
+
+    agg_obs_cube = aggregate._make_aggregated_obs_cube(
+        dummy_cubelist_obs_3_common_stations,
+        station_lookup,
+        common_stations,
+        forecast_period,
+    )
+    assert agg_obs_cube.shape == (2, 10, 3)
diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py
index dc2cc1873..2d476666e 100644
--- a/tests/operators/test_scoreswrappers.py
+++ b/tests/operators/test_scoreswrappers.py
@@ -975,3 +975,300 @@ def test_pfd_gt_complete_miss(make_cube_categorical_testing):
 
     assert len(result) == 1
     assert np.allclose(result[0].data, 1.0, atol=1e-2, rtol=1e-6)
+
+
+def test_rmse_multiple_forecasts_preserve_forecast_reference_time(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Testing RMSE aggregated by frt."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_rmse = []
+    for i in range(obs.coord("forecast_reference_time").shape[0]):
+        calculate_rmse.append(
+            np.sqrt(np.mean((data_obs[i, :, :] - data_model[i, :, :]) ** 2))
+        )
+
+    rmse_scores = scoreswrappers.scores_rmse(input_cubelist, "forecast_reference_time")
+    assert np.allclose(rmse_scores.data, calculate_rmse, atol=1e-2, rtol=1e-6)
+
+
+def test_rmse_multiple_forecasts_preserve_forecast_period(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test RMSE aggregated by forecast period."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_rmse = []
+    for i in range(obs.coord("forecast_period").shape[0]):
+        calculate_rmse.append(
+            np.sqrt(np.mean((data_obs[:, i, :] - data_model[:, i, :]) ** 2))
+        )
+
+    rmse_scores = scoreswrappers.scores_rmse(input_cubelist, "forecast_period")
+    assert np.allclose(rmse_scores.data, calculate_rmse, atol=1e-2, rtol=1e-6)
+
+
+def test_rmse_multiple_forecasts_preserve_none(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test RMSE preserving no coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_rmse = np.sqrt(np.mean((data_obs[:, :, :] - data_model[:, :, :]) ** 2))
+
+    rmse_scores = scoreswrappers.scores_rmse(input_cubelist)
+    assert np.allclose(rmse_scores.data, calculate_rmse, atol=1e-2, rtol=1e-6)
+
+
+def test_rmse_multiple_forecasts_preserve_lat_lon(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test RMSE preserving lat/lon coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_rmse = []
+    for i in range(obs.coord("station").shape[0]):
+        calculate_rmse.append(
+            np.sqrt(np.mean((data_obs[:, :, i] - data_model[:, :, i]) ** 2))
+        )
+
+    rmse_scores = scoreswrappers.scores_rmse(
+        input_cubelist, preserved_coordinates=["latitude", "longitude"]
+    )
+
+    assert np.allclose(rmse_scores.data, calculate_rmse, atol=1e-2, rtol=1e-6)
+
+
+def test_mae_multiple_forecasts_preserve_forecast_reference_time(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Testing MAE aggregated by frt."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_mae = []
+    for i in range(obs.coord("forecast_reference_time").shape[0]):
+        calculate_mae.append(np.mean(np.abs(data_obs[i, :, :] - data_model[i, :, :])))
+
+    mae_scores = scoreswrappers.scores_mae(input_cubelist, "forecast_reference_time")
+    assert np.allclose(mae_scores.data, calculate_mae, atol=1e-2, rtol=1e-6)
+
+
+def test_mae_multiple_forecasts_preserve_forecast_period(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test MAE aggregated by forecast period."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_mae = []
+    for i in range(obs.coord("forecast_period").shape[0]):
+        calculate_mae.append(np.mean(np.abs(data_obs[:, i, :] - data_model[:, i, :])))
+
+    mae_scores = scoreswrappers.scores_mae(input_cubelist, "forecast_period")
+    assert np.allclose(mae_scores.data, calculate_mae, atol=1e-2, rtol=1e-6)
+
+
+def test_mae_multiple_forecasts_preserve_none(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test MAE preserving no coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_mae = np.mean(np.abs(data_obs[:, :, :] - data_model[:, :, :]))
+
+    mae_scores = scoreswrappers.scores_mae(input_cubelist)
+    assert np.allclose(mae_scores.data, calculate_mae, atol=1e-2, rtol=1e-6)
+
+
+def test_mae_multiple_forecasts_preserve_lat_lon(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test mae preserving lat/lon coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_mae = []
+    for i in range(obs.coord("station").shape[0]):
+        calculate_mae.append(np.mean(np.abs(data_obs[:, :, i] - data_model[:, :, i])))
+
+    mae_scores = scoreswrappers.scores_mae(
+        input_cubelist, preserved_coordinates=["latitude", "longitude"]
+    )
+    assert np.allclose(mae_scores.data, calculate_mae, atol=1e-2, rtol=1e-6)
+
+
+def test_additive_bias_multiple_forecasts_preserve_forecast_reference_time(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Testing additive bias aggregated by frt."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_bias = []
+    for i in range(obs.coord("forecast_reference_time").shape[0]):
+        calculate_bias.append(np.mean(data_model[i, :, :] - data_obs[i, :, :]))
+
+    bias_scores = scoreswrappers.scores_additive_bias(
+        input_cubelist, "forecast_reference_time"
+    )
+    assert np.allclose(bias_scores.data, calculate_bias, atol=1e-2, rtol=1e-6)
+
+
+def test_additive_bias_multiple_forecasts_preserve_forecast_period(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test additive bias aggregated by forecast period."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_bias = []
+    for i in range(obs.coord("forecast_period").shape[0]):
+        calculate_bias.append(np.mean(data_model[:, i, :] - data_obs[:, i, :]))
+
+    bias_scores = scoreswrappers.scores_additive_bias(input_cubelist, "forecast_period")
+    assert np.allclose(bias_scores.data, calculate_bias, atol=1e-2, rtol=1e-6)
+
+
+def test_additive_bias_multiple_forecasts_preserve_none(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test additive bias preserving no coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_bias = np.mean(data_model[:, :, :] - data_obs[:, :, :])
+
+    bias_scores = scoreswrappers.scores_additive_bias(input_cubelist)
+    assert np.allclose(bias_scores.data, calculate_bias, atol=1e-2, rtol=1e-6)
+
+
+def test_additive_bias_multiple_forecasts_preserve_lat_lon(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test additive bias preserving lat/lon coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_bias = []
+    for i in range(obs.coord("station").shape[0]):
+        calculate_bias.append(np.mean(data_model[:, :, i] - data_obs[:, :, i]))
+
+    bias_scores = scoreswrappers.scores_additive_bias(
+        input_cubelist, preserved_coordinates=["latitude", "longitude"]
+    )
+    assert np.allclose(bias_scores.data, calculate_bias, atol=1e-2, rtol=1e-6)
+
+
+def test_correlation_pearsonr_multiple_forecasts_preserve_forecast_reference_time(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Testing Pearson correlation aggregated by frt."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_corr = []
+    for i in range(obs.coord("forecast_reference_time").shape[0]):
+        calculate_corr.append(
+            np.corrcoef(data_obs[i, :, :].flatten(), data_model[i, :, :].flatten())[
+                0, 1
+            ]
+        )
+
+    corr_scores = scoreswrappers.scores_correlation_pearsonr(
+        input_cubelist, "forecast_reference_time"
+    )
+    assert np.allclose(corr_scores.data, calculate_corr, atol=1e-2, rtol=1e-6)
+
+
+def test_correlation_pearsonr_multiple_forecasts_preserve_forecast_period(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test Pearson correlation aggregated by forecast period."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_corr = []
+    for i in range(obs.coord("forecast_period").shape[0]):
+        calculate_corr.append(
+            np.corrcoef(data_obs[:, i, :].flatten(), data_model[:, i, :].flatten())[
+                0, 1
+            ]
+        )
+
+    corr_scores = scoreswrappers.scores_correlation_pearsonr(
+        input_cubelist, "forecast_period"
+    )
+    assert np.allclose(corr_scores.data, calculate_corr, atol=1e-2, rtol=1e-6)
+
+
+def test_correlation_pearsonr_multiple_forecasts_preserve_none(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test Pearson correlation preserving no coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_corr = np.corrcoef(data_obs.flatten(), data_model.flatten())[0, 1]
+
+    corr_scores = scoreswrappers.scores_correlation_pearsonr(input_cubelist)
+    assert np.allclose(corr_scores.data, calculate_corr, atol=1e-2, rtol=1e-6)
+
+
+def test_correlation_pearsonr_multiple_forecasts_preserve_lat_lon(
+    dummy_cubelist_model_obs_multiple_forecasts,
+):
+    """Test Pearson correlation preserving lat/lon coordinates."""
+    input_cubelist = dummy_cubelist_model_obs_multiple_forecasts
+    obs, model = input_cubelist
+    data_obs = obs.data
+    data_model = model.data
+
+    calculate_corr = []
+    for i in range(obs.coord("station").shape[0]):
+        calculate_corr.append(
+            np.corrcoef(data_obs[:, :, i].flatten(), data_model[:, :, i].flatten())[
+                0, 1
+            ]
+        )
+
+    corr_scores = scoreswrappers.scores_correlation_pearsonr(
+        input_cubelist, preserved_coordinates=["latitude", "longitude"]
+    )
+    assert np.allclose(corr_scores.data, calculate_corr, atol=1e-2, rtol=1e-6)
diff --git a/tests/operators/test_time_utils.py b/tests/operators/test_time_utils.py
new file mode 100644
index 000000000..e69de29bb