Coupled tensor approximation has recently emerged as a promising approach for the fusion of hyperspectral and multispectral images, reconciling state of the art performance with strong theoretical guarantees. However, tensor-based approaches previously proposed assume that the different observed images are acquired under exactly the same conditions. A recent work proposed to accommodate inter-image spectral variability in the image fusion problem using a matrix factorization-based formulation, but did not account for spatially-localized variations. Moreover, it lacks theoretical guarantees and has a high associated computational complexity. In this paper, we consider the image fusion problem while accounting for both spatially and spectrally localized changes in an additive model. We first study how the general identifiability of the model is impacted by the presence of such changes. Then, assuming that the high-resolution image and the variation factors admit a Tucker decomposition, two new algorithms are proposed -- one purely algebraic, and another based on an optimization procedure. Theoretical guarantees for the exact recovery of the high-resolution image are provided for both algorithms. Experimental results show that the proposed method outperforms state-of-the-art methods in the presence of spectral and spatial variations between the images, at a smaller computational cost.