diff --git a/lib/matplotlib/image.py b/lib/matplotlib/image.py
index c1846f92608c..8a55ea868850 100644
--- a/lib/matplotlib/image.py
+++ b/lib/matplotlib/image.py
@@ -512,8 +512,10 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
                     if A.shape[2] == 3:  # image has no alpha channel
                         A = np.dstack([A, np.ones(A.shape[:2])])
                 elif np.ndim(alpha) > 0:  # Array alpha
-                    # user-specified array alpha overrides the existing alpha channel
-                    A = np.dstack([A[..., :3], alpha])
+                    if A.shape[2] == 3:  # image has no alpha channel
+                        A = np.dstack([A, alpha])
+                    else:  # blend with existing alpha channel
+                        A = np.dstack([A[..., :3], A[..., 3:] * alpha[..., np.newaxis]])
                 else:  # Scalar alpha
                     if A.shape[2] == 3:  # broadcast scalar alpha
                         A = np.dstack([A, np.full(A.shape[:2], alpha, np.float32)])
diff --git a/lib/matplotlib/tests/test_image.py b/lib/matplotlib/tests/test_image.py
index 330a2fab503d..fd2e9a8bd985 100644
--- a/lib/matplotlib/tests/test_image.py
+++ b/lib/matplotlib/tests/test_image.py
@@ -1809,19 +1809,21 @@ def test_interpolation_stage_rgba_respects_alpha_param(fig_test, fig_ref, intp_s
     axs_ref[0][2].imshow(im_rgba, interpolation_stage=intp_stage)
 
     # When the image already has an alpha channel, multiply it by the
-    # scalar alpha param, or replace it by the array alpha param
+    # scalar alpha param, or blend it with the array alpha param
     axs_tst[1][0].imshow(im_rgba)
     axs_ref[1][0].imshow(im_rgb, alpha=array_alpha)
     axs_tst[1][1].imshow(im_rgba, interpolation_stage=intp_stage, alpha=scalar_alpha)
     axs_ref[1][1].imshow(
-        np.concatenate(  # combine rgb channels with scaled array alpha
-            (im_rgb, scalar_alpha * array_alpha.reshape((ny, nx, 1))), axis=-1
-        ), interpolation_stage=intp_stage
-    )
+        np.concatenate(
+            (im_rgb,
+             (scalar_alpha * array_alpha).reshape((ny, nx, 1))),
+            axis=-1),
+        interpolation_stage=intp_stage)
     new_array_alpha = np.random.rand(ny, nx)
     axs_tst[1][2].imshow(im_rgba, interpolation_stage=intp_stage, alpha=new_array_alpha)
     axs_ref[1][2].imshow(
-        np.concatenate(  # combine rgb channels with new array alpha
-            (im_rgb, new_array_alpha.reshape((ny, nx, 1))), axis=-1
-        ), interpolation_stage=intp_stage
-    )
+        np.concatenate(
+            (im_rgb,
+             (array_alpha * new_array_alpha).reshape((ny, nx, 1))),
+            axis=-1),
+        interpolation_stage=intp_stage)