I tested 8 possible kernel_size from [1, 1, 1] to [3, 3, 3] and found that only [3, 1, 1] has the issue.
So, as a workaround, I considered to use Conv3d with [1, 1, 3] kernel and transpose the input/output tensor before/after the convolution inside the forward method as below.
class Conv3d(nn.Module):
    def __init__(self, in_channels, out_channels, kernel_size, stride=1, padding=0, dilation=1, groups=1, bias=True,
              padding_mode='zeros'):
        super(Conv3d, self).__init__()
        print('Use Conv3d Module')
        self.permute = False
        ks, sr, pad, dl = kernel_size, stride, padding, dilation
        if isinstance(kernel_size, (list, tuple)):
            if kernel_size == (3, 1, 1) or kernel_size == [3, 1, 1]:
                self.permute = True
                ks = ks[::-1]
                if isinstance(sr, (list, tuple)): sr = sr[::-1]
                if isinstance(pad, (list, tuple)): pad = pad[::-1]
                if isinstance(dl, (list, tuple)):  dl = dl[::-1]
        self.conv = nn.Conv3d(in_channels, out_channels, ks, stride=sr, padding=pad, dilation=dl,
                  groups=groups, bias=bias, padding_mode=padding_mode)

    def forward(self, x):
        if self.permute:
            return self.conv(x.permute(0,1,4,3,2)).permute(0,1,4,3,2)
        else:
            return self.conv(x)

It seems the MSE is reduced as normal in the previous test code.
@ptrblck, do you think it can be a possible workaround? Will the permutation increase the computation time a lot?
If the operation is functionally same as the original Conv3d, I’m considering testing it in my model, though the computation time will increase somewhat.