We present GAUSSN, a Bayesian semiparametric Gaussian Process (GP) model for time-delay estimation with resolved systems of gravitationally lensed supernovae (glSNe). GAUSSN models the underlying light curve non-parametrically using a GP. Without assuming a template light curve for each SN type, GAUSSN fits for the time delays of all images using data in any number of wavelength filters simultaneously. We also introduce a novel time-varying magnification model to capture the effects of microlensing alongside time-delay estimation. In this analysis, we model the time-varying relative magnification as a sigmoid function, as well as a constant for comparison to existing time-delay estimation approaches. We demonstrate that GAUSSN provides robust time-delay estimates for simulations of glSNe from the Nancy Grace Roman Space Telescope and the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (Rubin-LSST). We find that up to 43.6 per cent of time-delay estimates from Roman and 52.9 per cent from Rubin-LSST have fractional errors of less than 5  per cent. We then apply GAUSSN to SN Refsdal and find the time delay for the fifth image is consistent with the original analysis, regardless of microlensing treatment. Therefore, GAUSSN maintains the level of precision and accuracy achieved by existing time-delay extraction methods with fewer assumptions about the underlying shape of the light curve than template-based approaches, while incorporating microlensing into the statistical error budget. GAUSSN is scalable for time-delay cosmography analyses given current projections of glSNe discovery rates from Rubin-LSST and Roman.