diff --git a/lib/mpl_toolkits/mplot3d/axes3d.py b/lib/mpl_toolkits/mplot3d/axes3d.py
index 55b204022fb9..bc327d7c490b 100644
--- a/lib/mpl_toolkits/mplot3d/axes3d.py
+++ b/lib/mpl_toolkits/mplot3d/axes3d.py
@@ -2138,8 +2138,10 @@ def fill_between(self, x1, y1, z1, x2, y2, z2, *,
         self.auto_scale_xyz([x1, x2], [y1, y2], [z1, z2], had_data)
         return polyc
 
-    def plot_surface(self, X, Y, Z, *, norm=None, vmin=None,
-                     vmax=None, lightsource=None, axlim_clip=False, **kwargs):
+    def plot_surface(
+    self, X, Y, Z, *, norm=None, vmin=None, vmax=None, lightsource=None,
+    axlim_clip=False, interpolation=None, **kwargs
+):
         """
         Create a surface plot.
 
@@ -2249,6 +2251,28 @@ def plot_surface(self, X, Y, Z, *, norm=None, vmin=None,
 
         fcolors = kwargs.pop('facecolors', None)
 
+        if fcolors is not None and interpolation is not None:
+            try:
+                from scipy.ndimage import zoom
+            except ImportError:
+                raise ImportError("scipy is required for interpolation of facecolors")
+            fc = np.asarray(fcolors)
+            zoom_factors = [
+                X.shape[0] / fc.shape[0],
+                X.shape[1] / fc.shape[1]
+            ]
+            order = (
+                1 if interpolation == 'bilinear'
+                else 3 if interpolation == 'bicubic'
+                else 0
+            )
+            fc_interp = zoom(fc, zoom_factors, order=order)
+            if fc_interp.shape[2] not in (3, 4):
+                raise ValueError(
+                    "Interpolated facecolors must have 3 or 4 channels (RGB or RGBA)"
+                )
+            fcolors = fc_interp.astype(float)
+
         cmap = kwargs.get('cmap', None)
         shade = kwargs.pop('shade', cmap is None)
         if shade is None:
@@ -2308,7 +2332,7 @@ def plot_surface(self, X, Y, Z, *, norm=None, vmin=None,
 
         if fcolors is not None:
             polyc = art3d.Poly3DCollection(
-                polys, edgecolors=colset, facecolors=colset, shade=shade,
+                polys, facecolors=fcolors, shade=shade,
                 lightsource=lightsource, axlim_clip=axlim_clip, **kwargs)
         elif cmap:
             polyc = art3d.Poly3DCollection(polys, axlim_clip=axlim_clip, **kwargs)
@@ -2503,12 +2527,15 @@ def plot_trisurf(self, *args, color=None, norm=None, vmin=None, vmax=None,
         shade : bool, default: True
             Whether to shade the facecolors.  Shading is always disabled when
             *cmap* is specified.
+
         lightsource : `~matplotlib.colors.LightSource`, optional
             The lightsource to use when *shade* is True.
+
         axlim_clip : bool, default: False
             Whether to hide patches with a vertex outside the axes view limits.
 
             .. versionadded:: 3.10
+
         **kwargs
             All other keyword arguments are passed on to
             :class:`~mpl_toolkits.mplot3d.art3d.Poly3DCollection`
diff --git a/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py b/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
index cd45c8e33a6f..8f43d1949851 100644
--- a/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
+++ b/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
@@ -652,6 +652,25 @@ def test_surface3d():
     fig.colorbar(surf, shrink=0.5, aspect=5)
 
 
+@mpl3d_image_comparison(['surface3d_facecolors_interp.png'], style='mpl20')
+def test_surface3d_facecolors_interpolation():
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+    X = np.linspace(-2, 2, 5)
+    Y = np.linspace(-2, 2, 5)
+    X, Y = np.meshgrid(X, Y)
+    Z = np.sin(X**2 + Y**2)
+    # Facecolors: low-res RGB grid
+    C = np.zeros((4, 4, 3))
+    C[..., 0] = np.linspace(0, 1, 4)[:, None]  # Red gradient
+    C[..., 1] = np.linspace(1, 0, 4)[None, :]  # Green gradient
+    C[..., 2] = 0.5  # Blue constant
+
+    # Test bicubic interpolation
+    ax.plot_surface(X, Y, Z, facecolors=C, interpolation='bicubic', shade=False)
+    ax.set_title("plot_surface with facecolors and bicubic interpolation")
+
+
 @image_comparison(['surface3d_label_offset_tick_position.png'], style='mpl20')
 def test_surface3d_label_offset_tick_position():
     plt.rcParams['axes3d.automargin'] = True  # Remove when image is regenerated