The common method for purification of macromolecular bioproducts is preparative packed-bed chromatography using polymer-based, compressible, viscoelastic resins. Because of a downstream processing bottleneck, the chromatography equipment is often operated at its hydrodynamic limit. In this case, the resins may exhibit a complex behavior which results in compression-relaxation hystereses. Up to now, no modeling approach of transient flow through a chromatography packing has been made considering the viscoelasticity of the resins. The aim of the present work was to develop a novel model and compare model calculations with experimental data of two agarose-based resins. Fluid flow and bed permeability were modeled by Darcy's law and the Kozeny-Carman equation, respectively. Fluid flow was coupled to solid matrix stress via an axial force balance and a continuity equation of a deformable packing. Viscoelasticity was considered according to a Kelvin-Voigt material. The coupled equations were solved with a finite difference scheme using a deformable mesh. The model boundary conditions were preset transient pressure drop functions which resemble simulated load/elution/equilibration cycles. Calculations using a homogeneous model (assuming constant variables along the column height) gave a fair agreement with experimental data with regard to predicted flow rate, bed height, and compression-relaxation hysteresis for symmetric as well as asymmetric pressure drop functions. Calculations using an inhomogeneous model gave profiles of the bed porosity as a function of the bed height. In addition, the influence of medium wall support and intraparticle porosity was illustrated. The inhomogeneous model provides insights that so far are not easily experimentally accessible.