Conventional ultrasound (US) imaging relies on delay-and-sum (DAS) beamforming which retrieves a radio- frequency (RF) image, a blurred estimate of the tissue reflectivity function (TRF). Despite the non-stationarity of the blur induced by propagation effects, most state-of-the-art US restoration approaches exploit shift-invariant models and are inaccurate in realistic situations. Recent techniques approximate the shift- variant blur using sectional methods resulting in improved accuracy. But such methods assume shift-invariance of the blur in the lateral dimension which is not valid in many US imaging configurations. In this work, we propose a physical model of the non-stationary blur, which accounts for the diffraction effects related to the propagation. We show that its evaluation results in the sequential application of a forward and an adjoint propagation operators under some specific assumptions that we define. Taking into account this sequential structure, we exploit efficient formulations of the operators in the discrete domain and provide an evaluation strategy which exhibits linear complexity with respect to the grid size. We also show that the proposed model can be interpreted in terms of common simplification strategies used to model non-stationary blur. Through simulations and in vivo experimental data, we demonstrate that using the proposed model in the context of maximum-a-posteriori image restoration results in higher image quality than using state-of-the-art shift-invariant models. The supporting code is available on github: https://github.com/LTS5/us-non-stationary-deconv.