3D Convolution with sum of multiple 2D Convolution operation

Question

Looking at the code below, why res3d == res2d?

In other words, why I can simulate a 3D convolution operation as the sum of multiple 2D convolution operations?

import numpy as np
from scipy import signal 

X,Y,Z = 10,10,10  # can be anything
k_x,k_y,k_z = 3,3,3  # can be anything

img3d = np.random.rand(X,Y,Z)

f = np.random.randint(10,size=(1,k_x))
g = np.random.randint(10,size=(k_y,1))
h = np.random.randint(10,size=(k_z,1,1))
flt = g*f*h

########################################
# Conv3D 
res3d = signal.convolve(img3d, flt, mode='same')

########################################
# Multiple Conv 2D to replicate Conv 3D
(l_i, d_i, h_i) = img3d.shape
(l_k, d_k, h_k) = flt.shape

res2d = np.zeros(img3d.shape)

for i in range(l_k):
  for j in range(l_i):
    o = j - (i - l_k//2)
    if o < 0 or o >= l_i: continue
    res2d[o,...] += signal.convolve2d(img3d[j,...],flt[l_k-i-1,...],mode='same')

PS. I've already asked the question here but with no answers...