diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
index 02505c7d..b8d7af64 100644
--- a/doc/source/changelog.rst
+++ b/doc/source/changelog.rst
@@ -22,6 +22,8 @@ v1.6.2
   
   * Fixed bug with interactive widget contour plot and update of colorbar.
 
+* ``MatplotlibBackend`` can now combine 3d plots with contour plots.
+
 * Fixed bug with addition of interactive plots.
 
 
diff --git a/spb/backends/matplotlib.py b/spb/backends/matplotlib.py
index 4c90837a..5f05db54 100644
--- a/spb/backends/matplotlib.py
+++ b/spb/backends/matplotlib.py
@@ -249,11 +249,13 @@ def _create_figure(self):
 
             is_3D = [s.is_3D for s in self.series]
             if any(is_3D) and (not all(is_3D)):
-                raise ValueError(
-                    "The matplotlib backend can not mix 2D and 3D.")
+                # allow sum of 3D plots with contour plots
+                if not all(s.is_3D or s.is_contour for s in self.series):
+                    raise ValueError(
+                        "MatplotlibBackend can not mix 2D and 3D.")
 
             kwargs = dict(aspect=aspect)
-            if all(is_3D):
+            if any(is_3D):
                 kwargs["projection"] = "3d"
             elif any(s.is_2Dline and s.is_polar for s in self.series):
                 kwargs["projection"] = "polar"
diff --git a/tests/test_backends.py b/tests/test_backends.py
index 9d0abf8e..5963a16b 100644
--- a/tests/test_backends.py
+++ b/tests/test_backends.py
@@ -4295,3 +4295,33 @@ def test_plot3d_list_interactive():
 
     p = _p2(KB) + _p1(KB)
     p.backend._update_interactive({t: 1})
+
+
+def test_matplotlib_contour_and_3d():
+    # verify that it's possible to combine contour and 3d plots, but that
+    # combining a 2d line plot with contour and 3d plot raises an error.
+
+    x, y = symbols("x, y")
+    expr = cos(x * y) * exp(-0.05 * (x**2 + y**2))
+    ranges = (x, -5, 5), (y, -5, 5)
+
+    p1 = plot3d(expr, *ranges,
+        show=False, legend=True, zlim=(-2, 1), n=4)
+    p2 = plot_contour(
+        expr, *ranges,
+        {"zdir": "z", "offset": -2, "levels": 5},
+        show=False, is_filled=False, legend=True, n=4)
+    p3 = plot(cos(x), (x, 0, 2*pi), adaptive=False, n=5, show=False)
+
+    p = p1 + p2
+    p.process_series()
+    p = p2 + p1
+    p.process_series()
+    p = p2 + p3
+    p.process_series()
+    p = p1 + p3
+    raises(ValueError, lambda : p.process_series())
+    p = p1 + p2 + p3
+    raises(ValueError, lambda : p.process_series())
+    p = p2 + p1 + p3
+    raises(ValueError, lambda : p.process_series())