diff --git a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf
index 3032d9d7e..1c8c5f184 100644
--- a/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf
+++ b/src/CSET/cset_workflow/meta/diagnostics/rose-meta.conf
@@ -118,6 +118,25 @@ type=python_boolean,python_boolean,python_boolean,python_boolean
 compulsory=true
 sort-key=0surface4
 
+[template variables=TIMESERIES_SURFACE_FIELD_LAND_MASK]
+ns=Diagnostics/Fields
+description=Create time series plots for specific surface fields filtered
+            with a land mask. Calculated as domain (or sub-area) mean.
+            Requires: land_binary_mask.
+type=python_boolean
+compulsory=true
+sort-key=0surface4
+
+[template variables=TIMESERIES_SURFACE_FIELD_SEA_MASK]
+ns=Diagnostics/Fields
+description=Create time series plots for specific surface fields filtered
+            with a sea mask. Calculated as domain (or sub-area) mean.
+            Requires: land_binary_mask.
+type=python_boolean
+compulsory=true
+sort-key=0surface4
+
+
 [template variables=HISTOGRAM_SURFACE_FIELD]
 ns=Diagnostics/Fields
 description=Create histogram plots of specified surface fields.
diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example
index ed16ab88a..c734660cb 100644
--- a/src/CSET/cset_workflow/rose-suite.conf.example
+++ b/src/CSET/cset_workflow/rose-suite.conf.example
@@ -175,6 +175,8 @@ TIMESERIES_PLEVEL_FIELD=False
 TIMESERIES_PLEVEL_FIELD_AGGREGATION=False,False,False,False
 TIMESERIES_SURFACE_FIELD=False
 TIMESERIES_SURFACE_FIELD_AGGREGATION=False,False,False,False
+TIMESERIES_SURFACE_FIELD_LAND_MASK=False
+TIMESERIES_SURFACE_FIELD_SEA_MASK=False
 !!USE_WMO_STATION_NUMBERS=False
 !!VERTICAL_COORDINATE_A=[]
 !!VERTICAL_COORDINATE_B=[]
diff --git a/src/CSET/loaders/timeseries.py b/src/CSET/loaders/timeseries.py
index e2ace9781..c36a21379 100644
--- a/src/CSET/loaders/timeseries.py
+++ b/src/CSET/loaders/timeseries.py
@@ -41,6 +41,40 @@ def load(conf: Config):
                 aggregation=False,
             )
 
+    # Surface time series: land mask.
+    if conf.TIMESERIES_SURFACE_FIELD_LAND_MASK:
+        for field in conf.SURFACE_FIELDS:
+            yield RawRecipe(
+                recipe="land_mask_for_surface_domain_mean_time_series.yaml",
+                variables={
+                    "VARNAME": field,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=False,
+            )
+
+    # Surface time series: sea mask.
+    if conf.TIMESERIES_SURFACE_FIELD_SEA_MASK:
+        for field in conf.SURFACE_FIELDS:
+            yield RawRecipe(
+                recipe="sea_mask_for_surface_domain_mean_time_series.yaml",
+                variables={
+                    "VARNAME": field,
+                    "MODEL_NAME": [model["name"] for model in models],
+                    "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None,
+                    "SUBAREA_EXTENT": conf.SUBAREA_EXTENT
+                    if conf.SELECT_SUBAREA
+                    else None,
+                },
+                model_ids=[model["id"] for model in models],
+                aggregation=False,
+            )
+
     # Pressure level fields.
     if conf.TIMESERIES_PLEVEL_FIELD:
         for field, plevel in itertools.product(
diff --git a/src/CSET/operators/filters.py b/src/CSET/operators/filters.py
index e7cf461e9..c2baeb01f 100644
--- a/src/CSET/operators/filters.py
+++ b/src/CSET/operators/filters.py
@@ -25,22 +25,22 @@
 
 
 def apply_mask(
-    original_field: iris.cube.Cube,
-    mask: iris.cube.Cube,
-) -> iris.cube.Cube:
+    original_field: iris.cube.Cube | iris.cube.CubeList,
+    mask: iris.cube.Cube | iris.cube.CubeList,
+) -> iris.cube.Cube | iris.cube.CubeList:
     """Apply a mask to given data as a masked array.
 
     Parameters
     ----------
-    original_field: iris.cube.Cube
-        The field to be masked.
-    mask: iris.cube.Cube
-        The mask being applied to the original field.
+    original_field: iris.cube.Cube | iris.cube.CubeList
+        The field(s) to be masked.
+    mask: iris.cube.Cube | iris.cube.CubeList
+        The mask(s) being applied to the original field(s).
 
     Returns
     -------
-    masked_field: iris.cube.Cube
-        A cube of the masked field.
+    masked_field: iris.cube.Cube | iris.cube.CubeList
+        A cube or cubelist of the masked field(s).
 
     Notes
     -----
@@ -54,17 +54,24 @@ def apply_mask(
     --------
     >>> land_points_only = apply_mask(temperature, land_mask)
     """
-    # Ensure mask is only 1s or NaNs.
-    mask.data[mask.data == 0] = np.nan
-    mask.data[~np.isnan(mask.data)] = 1
-    logging.info(
-        "Mask set to 1 or 0s, if addition of multiple masks results"
-        "in values > 1 these are set to 1."
-    )
-    masked_field = original_field.copy()
-    masked_field.data *= mask.data
-    masked_field.attributes["mask"] = f"mask_of_{original_field.name()}"
-    return masked_field
+    masked_fields = iris.cube.CubeList([])
+    for M, F in zip(iter_maybe(mask), iter_maybe(original_field), strict=True):
+        # Ensure mask data are floats and only 1s or NaNs.
+        M.data = np.float64(M.data)
+        M.data[M.data == 0.0] = np.nan
+        M.data[~np.isnan(M.data)] = 1.0
+        logging.info(
+            "Mask set to 1 or 0s, if addition of multiple masks results"
+            "in values > 1 these are set to 1."
+        )
+        masked_field = F.copy()
+        masked_field.data *= M.data
+        masked_field.attributes["mask"] = f"mask_of_{F.name()}"
+        masked_fields.append(masked_field)
+    if len(masked_fields) == 1:
+        return masked_fields[0]
+    else:
+        return masked_fields
 
 
 def filter_cubes(
@@ -201,17 +208,17 @@ def generate_mask(
         mask.data[:] = 0.0
         match condition:
             case "eq":
-                mask.data[cube.data == value] = 1
+                mask.data[cube.data == value] = 1.0
             case "ne":
-                mask.data[cube.data != value] = 1
+                mask.data[cube.data != value] = 1.0
             case "gt":
-                mask.data[cube.data > value] = 1
+                mask.data[cube.data > value] = 1.0
             case "ge":
-                mask.data[cube.data >= value] = 1
+                mask.data[cube.data >= value] = 1.0
             case "lt":
-                mask.data[cube.data < value] = 1
+                mask.data[cube.data < value] = 1.0
             case "le":
-                mask.data[cube.data <= value] = 1
+                mask.data[cube.data <= value] = 1.0
             case _:
                 raise ValueError("""Unexpected value for condition. Expected eq, ne,
                                   gt, ge, lt, le. Got {condition}.""")
diff --git a/src/CSET/recipes/surface_fields/land_mask_for_surface_domain_mean_time_series.yaml b/src/CSET/recipes/surface_fields/land_mask_for_surface_domain_mean_time_series.yaml
new file mode 100644
index 000000000..51fb8d1c3
--- /dev/null
+++ b/src/CSET/recipes/surface_fields/land_mask_for_surface_domain_mean_time_series.yaml
@@ -0,0 +1,47 @@
+category: Surface Time Series
+title: Land mask domain mean $VARNAME time series
+description: |
+  Creates a land mask and applies it to $VARNAME. It then creates and plots
+  a time series of the domain horizontal mean.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: filters.apply_mask
+    original_field:
+      operator: filters.filter_multiple_cubes
+      constraint:
+        operator: constraints.combine_constraints
+        varname_constraint:
+          operator: constraints.generate_var_constraint
+          varname: $VARNAME
+        cell_methods_constraint:
+          operator: constraints.generate_cell_methods_constraint
+          cell_methods: []
+          varname: $VARNAME
+        pressure_level_constraint:
+          operator: constraints.generate_level_constraint
+          coordinate: "pressure"
+          levels: []
+    mask:
+      operator: filters.generate_mask
+      mask_field:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: "land_binary_mask"
+      condition: 'eq'
+      value: 1
+
+  - operator: collapse.collapse
+    coordinate: [grid_latitude, grid_longitude]
+    method: MEAN
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
+
+  - operator: plot.plot_line_series
diff --git a/src/CSET/recipes/surface_fields/sea_mask_for_surface_domain_mean_time_series.yaml b/src/CSET/recipes/surface_fields/sea_mask_for_surface_domain_mean_time_series.yaml
new file mode 100644
index 000000000..bfc42f11e
--- /dev/null
+++ b/src/CSET/recipes/surface_fields/sea_mask_for_surface_domain_mean_time_series.yaml
@@ -0,0 +1,47 @@
+category: Surface Time Series
+title: Sea mask domain mean $VARNAME time series
+description: |
+  Creates a sea mask and applies it to $VARNAME. It then creates and plots
+  a time series of the domain horizontal mean.
+
+steps:
+  - operator: read.read_cubes
+    file_paths: $INPUT_PATHS
+    model_names: $MODEL_NAME
+    subarea_type: $SUBAREA_TYPE
+    subarea_extent: $SUBAREA_EXTENT
+
+  - operator: filters.apply_mask
+    original_field:
+      operator: filters.filter_multiple_cubes
+      constraint:
+        operator: constraints.combine_constraints
+        varname_constraint:
+          operator: constraints.generate_var_constraint
+          varname: $VARNAME
+        cell_methods_constraint:
+          operator: constraints.generate_cell_methods_constraint
+          cell_methods: []
+          varname: $VARNAME
+        pressure_level_constraint:
+          operator: constraints.generate_level_constraint
+          coordinate: "pressure"
+          levels: []
+    mask:
+      operator: filters.generate_mask
+      mask_field:
+        operator: filters.filter_multiple_cubes
+        constraint:
+          operator: constraints.generate_var_constraint
+          varname: "land_binary_mask"
+      condition: 'eq'
+      value: 0
+
+  - operator: collapse.collapse
+    coordinate: [grid_latitude, grid_longitude]
+    method: MEAN
+
+  - operator: write.write_cube_to_nc
+    overwrite: True
+
+  - operator: plot.plot_line_series
diff --git a/tests/operators/test_filters.py b/tests/operators/test_filters.py
index 7e5e10580..674cf69d5 100644
--- a/tests/operators/test_filters.py
+++ b/tests/operators/test_filters.py
@@ -185,7 +185,7 @@ def test_generate_mask_equal_to(cube):
     """Generates a mask with values equal to a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data == 276] = 1
+    mask.data[cube.data == 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "eq", 276).data,
         mask.data,
@@ -198,7 +198,7 @@ def test_generate_mask_not_equal_to(cube):
     """Generates a mask with values not equal to a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data != 276] = 1
+    mask.data[cube.data != 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "ne", 276).data,
         mask.data,
@@ -211,7 +211,7 @@ def test_generate_mask_greater_than(cube):
     """Generates a mask with values greater than a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data > 276] = 1
+    mask.data[cube.data > 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "gt", 276).data,
         mask.data,
@@ -224,7 +224,7 @@ def test_generate_mask_greater_equal_to(cube):
     """Generates a mask with values greater than or equal to a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data >= 276] = 1
+    mask.data[cube.data >= 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "ge", 276).data,
         mask.data,
@@ -237,7 +237,7 @@ def test_generate_mask_less_than(cube):
     """Generates a mask with values less than a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data < 276] = 1
+    mask.data[cube.data < 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "lt", 276).data,
         mask.data,
@@ -250,7 +250,7 @@ def test_generate_mask_less_equal_to(cube):
     """Generates a mask with values less than or equal to a specified value."""
     mask = cube.copy()
     mask.data = np.zeros(mask.data.shape)
-    mask.data[cube.data <= 276] = 1
+    mask.data[cube.data <= 276] = 1.0
     assert np.allclose(
         filters.generate_mask(cube, "le", 276).data,
         mask.data,
@@ -267,7 +267,7 @@ def test_generate_mask_cube_list(cubes):
     for cube in cubes:
         mask = cube.copy()
         mask.data[:] = 0.0
-        mask.data[cube.data <= 276] = 1
+        mask.data[cube.data <= 276] = 1.0
         masks_calc.append(mask)
     for cube, mask in zip(masks, masks_calc, strict=True):
         assert np.allclose(cube.data, mask.data, rtol=1e-06, atol=1e-02)
@@ -276,8 +276,8 @@ def test_generate_mask_cube_list(cubes):
 def test_apply_mask(cube):
     """Apply a mask to a cube."""
     mask = filters.generate_mask(cube, "eq", 276)
-    mask.data[mask.data == 0] = np.nan
-    mask.data[~np.isnan(mask.data)] = 1
+    mask.data[mask.data == 0.0] = np.nan
+    mask.data[~np.isnan(mask.data)] = 1.0
     test_data = cube.copy()
     test_data.data *= mask.data
     assert np.allclose(
@@ -289,6 +289,21 @@ def test_apply_mask(cube):
     )
 
 
+def test_apply_mask_cubelist(cube):
+    """Apply a mask to a cube list."""
+    mask = filters.generate_mask(cube, "eq", 276)
+    mask.data[mask.data == 0.0] = np.nan
+    mask.data[~np.isnan(mask.data)] = 1.0
+    test_data = cube.copy()
+    test_data.data *= mask.data
+    expected_list = iris.cube.CubeList([test_data, test_data])
+    mask_list = iris.cube.CubeList([mask, mask])
+    input_list = iris.cube.CubeList([cube, cube])
+    actual_cubelist = filters.apply_mask(input_list, mask_list)
+    for cube, mask in zip(expected_list, actual_cubelist, strict=True):
+        assert np.allclose(cube.data, mask.data, rtol=1e-06, atol=1e-02, equal_nan=True)
+
+
 def test_generate_single_ensemble_member_constraint_reduced_member(ensemble_cube):
     """Remove a single ensemble member from a cube."""
     remove_member_constraint = (