diff --git a/.github/workflows/benchmark.yml b/.github/workflows/benchmark.yml
index df4a4fc9da..9ebc7a4c45 100644
--- a/.github/workflows/benchmark.yml
+++ b/.github/workflows/benchmark.yml
@@ -15,7 +15,7 @@ jobs:
     env:
       IRIS_TEST_DATA_LOC_PATH: benchmarks
       IRIS_TEST_DATA_PATH: benchmarks/iris-test-data
-      IRIS_TEST_DATA_VERSION: "2.13"
+      IRIS_TEST_DATA_VERSION: "2.14"
       # Lets us manually bump the cache to rebuild
       ENV_CACHE_BUILD: "0"
       TEST_DATA_CACHE_BUILD: "2"
diff --git a/.github/workflows/ci-docs-tests.yml b/.github/workflows/ci-docs-tests.yml
index 9e200c124e..c2e3fcfcf8 100644
--- a/.github/workflows/ci-docs-tests.yml
+++ b/.github/workflows/ci-docs-tests.yml
@@ -39,7 +39,7 @@ jobs:
         session: ["doctest", "gallery", "linkcheck"]
 
     env:
-      IRIS_TEST_DATA_VERSION: "2.13"
+      IRIS_TEST_DATA_VERSION: "2.14"
       ENV_NAME: "ci-docs-tests"
 
     steps:
diff --git a/.github/workflows/ci-tests.yml b/.github/workflows/ci-tests.yml
index 31cccfb0bb..18cb1f5e21 100644
--- a/.github/workflows/ci-tests.yml
+++ b/.github/workflows/ci-tests.yml
@@ -39,7 +39,7 @@ jobs:
         session: ["tests"]
 
     env:
-      IRIS_TEST_DATA_VERSION: "2.13"
+      IRIS_TEST_DATA_VERSION: "2.14"
       ENV_NAME: "ci-tests"
 
     steps:
diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
index a16ce598cc..0e4d51ffec 100644
--- a/docs/src/whatsnew/latest.rst
+++ b/docs/src/whatsnew/latest.rst
@@ -72,6 +72,10 @@ This document explains the changes made to Iris for this release
    :func:`numpy.percentile` keywords through the :obj:`~iris.analysis.PERCENTILE`
    aggregator. (:pull:`4791`)
 
+#. `@wjbenfold`_ and `@bjlittle`_ (reviewer) implemented
+   :func:`iris.plot.fill_between` and :func:`iris.quickplot.fill_between`.
+   (:issue:`3493`, :pull:`4647`)
+
 
 🐛 Bugs Fixed
 =============
diff --git a/lib/iris/plot.py b/lib/iris/plot.py
index 74e5d5788c..47c63dc173 100644
--- a/lib/iris/plot.py
+++ b/lib/iris/plot.py
@@ -645,7 +645,30 @@ def _u_object_from_v_object(v_object):
 
 
 def _get_plot_objects(args):
-    if len(args) > 1 and isinstance(
+    if len(args) > 2 and isinstance(
+        args[2], (iris.cube.Cube, iris.coords.Coord)
+    ):
+        # three arguments
+        u_object, v_object1, v_object2 = args[:3]
+        u1, v1 = _uv_from_u_object_v_object(u_object, v_object1)
+        _, v2 = _uv_from_u_object_v_object(u_object, v_object2)
+        args = args[3:]
+        if u1.size != v1.size or u1.size != v2.size:
+            msg = "The x and y-axis objects are not all compatible. They should have equal sizes but got ({}: {}), ({}: {}) and ({}: {})"
+            raise ValueError(
+                msg.format(
+                    u_object.name(),
+                    u1.size,
+                    v_object1.name(),
+                    v1.size,
+                    v_object2.name(),
+                    v2.size,
+                )
+            )
+        u = u1
+        v = (v1, v2)
+        v_object = (v_object1, v_object2)
+    elif len(args) > 1 and isinstance(
         args[1], (iris.cube.Cube, iris.coords.Coord)
     ):
         # two arguments
@@ -823,6 +846,52 @@ def _draw_1d_from_points(draw_method_name, arg_func, *args, **kwargs):
     return result
 
 
+def _draw_two_1d_from_points(draw_method_name, arg_func, *args, **kwargs):
+    """
+    This function is equivalend to _draw_two_1d_from_points but expects two
+    y-axis variables rather than one (such as is required for .fill_between). It
+    can't be used where the y-axis variables are string coordinates. The y-axis
+    variable provided first has precedence where the two differ on whether the
+    axis should be inverted or whether a map should be drawn.
+    """
+    # NB. In the interests of clarity we use "u" to refer to the horizontal
+    # axes on the matplotlib plot and "v" for the vertical axes.
+
+    # retrieve the objects that are plotted on the horizontal and vertical
+    # axes (cubes or coordinates) and their respective values, along with the
+    # argument tuple with these objects removed
+    u_object, v_objects, u, vs, args = _get_plot_objects(args)
+
+    v_object1, _ = v_objects
+    v1, v2 = vs
+
+    # if both u_object and v_object are coordinates then check if a map
+    # should be drawn
+    if (
+        isinstance(u_object, iris.coords.Coord)
+        and isinstance(v_object1, iris.coords.Coord)
+        and _can_draw_map([v_object1, u_object])
+    ):
+        # Replace non-cartopy subplot/axes with a cartopy alternative and set
+        # the transform keyword.
+        kwargs = _ensure_cartopy_axes_and_determine_kwargs(
+            u_object, v_object1, kwargs
+        )
+
+    axes = kwargs.pop("axes", None)
+    draw_method = getattr(axes if axes else plt, draw_method_name)
+    if arg_func is not None:
+        args, kwargs = arg_func(u, v1, v2, *args, **kwargs)
+        result = draw_method(*args, **kwargs)
+    else:
+        result = draw_method(u, v1, v2, *args, **kwargs)
+
+    # Invert y-axis if necessary.
+    _invert_yaxis(v_object1, axes)
+
+    return result
+
+
 def _replace_axes_with_cartopy_axes(cartopy_proj):
     """
     Replace non-cartopy subplot/axes with a cartopy alternative
@@ -1599,6 +1668,45 @@ def scatter(x, y, *args, **kwargs):
     return _draw_1d_from_points("scatter", _plot_args, *args, **kwargs)
 
 
+def fill_between(x, y1, y2, *args, **kwargs):
+    """
+    Plots y1 and y2 against x, and fills the space between them.
+
+    Args:
+
+    * x: :class:`~iris.cube.Cube` or :class:`~iris.coords.Coord`
+        A cube or a coordinate to plot on the x-axis.
+
+    * y1: :class:`~iris.cube.Cube` or :class:`~iris.coords.Coord`
+        First cube or a coordinate to plot on the y-axis.
+
+    * y2: :class:`~iris.cube.Cube` or :class:`~iris.coords.Coord`
+        Second cube or a coordinate to plot on the y-axis.
+
+    Kwargs:
+
+    * axes: :class:`matplotlib.axes.Axes`
+        The axes to use for drawing.  Defaults to the current axes if none
+        provided.
+
+    See :func:`matplotlib.pyplot.fill_between` for details of additional valid
+    keyword arguments.
+
+    """
+    # here we are more specific about argument types than generic 1d plotting
+    if not isinstance(x, (iris.cube.Cube, iris.coords.Coord)):
+        raise TypeError("x must be a cube or a coordinate.")
+    if not isinstance(y1, (iris.cube.Cube, iris.coords.Coord)):
+        raise TypeError("y1 must be a cube or a coordinate.")
+    if not isinstance(y1, (iris.cube.Cube, iris.coords.Coord)):
+        raise TypeError("y2 must be a cube or a coordinate.")
+    args = (x, y1, y2) + args
+    _plot_args = None
+    return _draw_two_1d_from_points(
+        "fill_between", _plot_args, *args, **kwargs
+    )
+
+
 # Provide convenience show method from pyplot
 show = plt.show
 
diff --git a/lib/iris/quickplot.py b/lib/iris/quickplot.py
index 2c4a94b1d0..14f9e5d2d5 100644
--- a/lib/iris/quickplot.py
+++ b/lib/iris/quickplot.py
@@ -311,5 +311,19 @@ def scatter(x, y, *args, **kwargs):
     return result
 
 
+def fill_between(x, y1, y2, *args, **kwargs):
+    """
+    Draws a labelled fill_between plot based on the given cubes or coordinates.
+
+    See :func:`iris.plot.fill_between` for details of valid arguments and
+    keyword arguments.
+
+    """
+    axes = kwargs.get("axes")
+    result = iplt.fill_between(x, y1, y2, *args, **kwargs)
+    _label_1d_plot(x, y1, axes=axes)
+    return result
+
+
 # Provide a convenience show method from pyplot.
 show = plt.show
diff --git a/lib/iris/tests/results/imagerepo.json b/lib/iris/tests/results/imagerepo.json
index 5ae8046c5b..28d6f0bb03 100644
--- a/lib/iris/tests/results/imagerepo.json
+++ b/lib/iris/tests/results/imagerepo.json
@@ -68,6 +68,10 @@
     "iris.tests.test_mapping.TestLowLevel.test_simple.0": "faa0e558855f9de7857a1ab16a85a51d36a1e55a854e58a5c13837096e8fe17a",
     "iris.tests.test_mapping.TestMappingSubRegion.test_simple.0": "b9913d90c66eca6ec66ec2f3689195aecf5b2f00392cb3496495e21da4db6c92",
     "iris.tests.test_mapping.TestUnmappable.test_simple.0": "fa81b54a817eca37817ec701857e3e64943e7bb41b806f996e817e006ee1b19b",
+    "iris.tests.test_plot.Test1dFillBetween.test_coord_coord.0": "f31432798cebcd87723835b4a5c5c2dbcf139c6c8cf4730bf3c36d801e380378",
+    "iris.tests.test_plot.Test1dFillBetween.test_coord_cube.0": "ea17352b92f0cbd42d6c8d25e59d36dc3a538d2bb2e42d26c6d2c2c8e4a1ce99",
+    "iris.tests.test_plot.Test1dFillBetween.test_cube_coord.0": "aff8e44af2019b3d3d03e0d1865e272cc1643de292db4b98c53c7ce5b0c37b2c",
+    "iris.tests.test_plot.Test1dFillBetween.test_cube_cube.0": "ea1761f695a09c0b70cc938d334b4e4f4c3671f2cd8b7996973c2c68e1c39e26",
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_coord.0": "8bfec2577e01a5a5ed013b4ac4521c94817d4e6d91ff63369c6d61991e3278cc",
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_coord_coord.0": "8fff941e7e01e1c2f801c878a41e5b0d85cf36e1837e2d9992c62f21769e6a4d",
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_coord_coord_map.0": "bbe0c214cd979dc3b05e4b68db0771b48698961b7962d2446e8ca5bb36716c6e",
@@ -75,6 +79,10 @@
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_cube.0": "8fffc1dc7e019c70f001b70ee4386de1814e7938837b6a7f84d07c9f15b02f21",
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_cube_coord.0": "8fffc1dc7e019c70f001b70ee4386de1814e7938837b6a7f84d07c9f15b02f21",
     "iris.tests.test_plot.Test1dPlotMultiArgs.test_cube_cube.0": "8ff8c0567a01b296e4019d2ff10b464bd4da6391943678e5879f7e3903e63f1c",
+    "iris.tests.test_plot.Test1dQuickplotFillBetween.test_coord_coord.0": "f314b2798ce3cd87723835a4a5c5c2dbcf139c6c8cf4730bd3c36d801c3c6378",
+    "iris.tests.test_plot.Test1dQuickplotFillBetween.test_coord_cube.0": "ea17352bd2f0cbd4256c8da5e59c36dc1a538d2b92e41d26ced2c2c8eca1ce99",
+    "iris.tests.test_plot.Test1dQuickplotFillBetween.test_cube_coord.0": "a3ffe44af6009b3d2907c8f1f6588f2cc96619e290fb4b88cd2c3ce590e3770c",
+    "iris.tests.test_plot.Test1dQuickplotFillBetween.test_cube_cube.0": "ea17e1f695a09c0b60cc938d334b4e4f4c3671f2cd8b7996973c2c69e1c31e26",
     "iris.tests.test_plot.Test1dQuickplotPlotMultiArgs.test_coord.0": "83fec2777e002427e801bb4ae65a1c94813dcec999db4bbc9ccd79991f3238cc",
     "iris.tests.test_plot.Test1dQuickplotPlotMultiArgs.test_coord_coord.0": "83ff9d9f7e01e1c2b001c8f8f63e1b1d81cf36e1837e258982ce6f215c9a626c",
     "iris.tests.test_plot.Test1dQuickplotPlotMultiArgs.test_coord_coord_map.0": "bbe0c214cd979dc3b05e4b68db0771b48698961b7962d2446e8ca5bb36716c6e",
diff --git a/lib/iris/tests/test_image_json.py b/lib/iris/tests/test_image_json.py
index 68f70753bf..b5213156f8 100644
--- a/lib/iris/tests/test_image_json.py
+++ b/lib/iris/tests/test_image_json.py
@@ -30,17 +30,17 @@ def test_json(self):
         missing_from_json = test_data_name_set - repo_name_set
         if missing_from_json:
             amsg = (
-                "Missing images: Image names are referenced in "
-                "imagerepo.json, that are not present in the iris-test-data "
-                "repo"
+                "Missing images: Images are present in the iris-test-data "
+                "repo, that are not referenced in imagerepo.json"
             )
             # Always fails when we get here: report the problem.
             self.assertEqual(missing_from_json, set(), msg=amsg)
         missing_from_test_data = repo_name_set - test_data_name_set
         if missing_from_test_data:
             amsg = (
-                "Missing images: Images are present in the iris-test-data "
-                "repo, that are not referenced in imagerepo.json"
+                "Missing images: Image names are referenced in "
+                "imagerepo.json, that are not present in the iris-test-data "
+                "repo"
             )
             # Always fails when we get here: report the problem.
             self.assertEqual(missing_from_test_data, set(), msg=amsg)
diff --git a/lib/iris/tests/test_plot.py b/lib/iris/tests/test_plot.py
index 458616a6fb..0c47bd6d3a 100644
--- a/lib/iris/tests/test_plot.py
+++ b/lib/iris/tests/test_plot.py
@@ -16,6 +16,7 @@
 import numpy as np
 
 import iris
+import iris.analysis
 import iris.coords as coords
 import iris.tests.stock
 
@@ -341,6 +342,108 @@ def test_circular_changes(self):
         self.check_graphic()
 
 
+class Test1dFillBetween(tests.GraphicsTest):
+    def setUp(self):
+        super().setUp()
+        self.cube = iris.load_cube(
+            tests.get_data_path(
+                ("NetCDF", "testing", "small_theta_colpex.nc")
+            ),
+            "air_potential_temperature",
+        )[0, 0]
+        self.draw_method = iplt.fill_between
+
+    def test_coord_coord(self):
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube.coord("surface_altitude")[:, 0]
+        y2 = self.cube.coord("surface_altitude")[:, 1]
+        self.draw_method(x, y1, y2)
+        self.check_graphic()
+
+    def test_coord_cube(self):
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)
+        self.draw_method(x, y1, y2)
+        self.check_graphic()
+
+    def test_cube_coord(self):
+        x = self.cube.collapsed("grid_longitude", iris.analysis.MEAN)
+        y1 = self.cube.coord("surface_altitude")[:, 0]
+        y2 = y1 + 10
+        self.draw_method(x, y1, y2)
+        self.check_graphic()
+
+    def test_cube_cube(self):
+        x = self.cube.collapsed("grid_longitude", iris.analysis.MEAN)
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)
+        self.draw_method(x, y1, y2)
+        self.check_graphic()
+
+    def test_incompatible_objects_x_odd(self):
+        # cubes/coordinates of different sizes cannot be plotted
+        x = self.cube.coord("grid_latitude")[:-1]
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)
+        with self.assertRaises(ValueError):
+            self.draw_method(x, y1, y2)
+
+    def test_incompatible_objects_y1_odd(self):
+        # cubes/coordinates of different sizes cannot be plotted
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)[:-1]
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)
+        with self.assertRaises(ValueError):
+            self.draw_method(x, y1, y2)
+
+    def test_incompatible_objects_y2_odd(self):
+        # cubes/coordinates of different sizes cannot be plotted
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)[:-1]
+        with self.assertRaises(ValueError):
+            self.draw_method(x, y1, y2)
+
+    def test_incompatible_objects_all_odd(self):
+        # cubes/coordinates of different sizes cannot be plotted
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)[:-1]
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)[:-2]
+        with self.assertRaises(ValueError):
+            self.draw_method(x, y1, y2)
+
+    def test_multidimensional(self):
+        # multidimensional cubes/coordinates are not allowed
+        x = self.cube.coord("grid_latitude")
+        y1 = self.cube
+        y2 = self.cube
+        with self.assertRaises(ValueError):
+            self.draw_method(x, y1, y2)
+
+    def test_not_cube_or_coord(self):
+        # inputs must be cubes or coordinates
+        x = np.arange(self.cube.shape[0])
+        y1 = self.cube.collapsed("grid_longitude", iris.analysis.MIN)
+        y2 = self.cube.collapsed("grid_longitude", iris.analysis.MAX)
+        with self.assertRaises(TypeError):
+            self.draw_method(x, y1, y2)
+
+
+@tests.skip_data
+@tests.skip_plot
+class Test1dQuickplotFillBetween(Test1dFillBetween):
+    def setUp(self):
+        tests.GraphicsTest.setUp(self)
+        self.cube = iris.load_cube(
+            tests.get_data_path(
+                ("NetCDF", "testing", "small_theta_colpex.nc")
+            ),
+            "air_potential_temperature",
+        )[0, 0]
+        self.draw_method = qplt.fill_between
+
+
 @tests.skip_data
 @tests.skip_plot
 class TestAttributePositive(tests.GraphicsTest):