Changes in soil pore volume and shape in response to internal and external mechanical stresses alter key soil hydrologic and transport properties. The extent of these changes is dependent on details of pore shape and size evolution. We present a model for quantifying rates of deformation and shape evolution of idealized spheroidal pores as functions of macroscopic stresses and soil rheological properties. Previous solutions for shrinkage of spherical pores embedded in a viscoplastic matrix under isotropic stress were extended to spheroidal pore shapes and biaxial stresses using Eshelby's classical theory. Bulk soil behavior was obtained from upscaling of detailed single pore deformation. Results show that pore closure rates increase with decreasing initial aspect ratio (i.e., oblate pores close faster than spherical pores), and with higher deviatoric stress. Incomplete pore closure is attributed to soil hardening due to pore shape accommodation under biaxial stresses. The model provides a means for approximating pore deformation as input to predictive models for soil hydraulic properties.