Existing ultrasound deconvolution approaches unrealistically assume, primarily for computational reasons, that the convolution model relies on a spatially invariant kernel and circulant boundary conditions. We discard both restrictions and introduce an image formation model applicable to ultrasound imaging and deconvolution based on an axially varying kernel, that accounts for arbitrary boundary conditions. Our model has the same computational complexity as the one employing spatially invariant convolution and has negligible memory requirements. To accommodate state-of-the-art deconvolution approaches when applied to a variety of inverse problem formulations, we also provide an equally efficient adjoint expression of our model. Simulation results confirm the tractability of our model for the deconvolution of large images. Moreover, in terms of accuracy metrics, the quality of reconstruction using our model is superior to that obtained using spatially invariant convolution.