Abstract. A thorough understanding of ice thermodynamics is essential for an accurate description of glaciers, ice sheets and ice shelves. Yet there exists a significant gap in our theoretical knowledge of the time-dependent behaviour of ice temperatures due to the inevitable compromise between mathematical tractability and the accurate description of physical phenomena. In order to bridge this shortfall, we have analytically solved the 1D time-dependent advective–diffusive heat problem including additional terms due to strain heating and depth-integrated horizontal advection. Newton's law of cooling is applied as a Robin-type top boundary condition to consider potential non-equilibrium temperature states across the ice–air interface. The solution is expressed in terms of confluent hypergeometric functions following a separation of variables approach. Non-dimensionalization reduces the parameter space to four numbers that fully determine the shape of the solution at equilibrium: surface insulation, effective geothermal heat flow, the Péclet number and the Brinkman number. The initial temperature distribution exponentially converges to the stationary solution. Transient decay timescales are only dependent on the Péclet number and the surface insulation, so higher advection rates and lower insulating values imply shorter equilibration timescales, respectively. In contrast, equilibrium temperature profiles are mostly independent of the surface insulation parameter. We have extended our study to a broader range of vertical velocities by using a general power-law dependence on depth, unlike prior studies limited to linear and quadratic velocity profiles. Lastly, we present a suite of benchmark experiments to test numerical solvers. Four experiments of gradually increasing complexity capture the main physical processes for heat propagation. Analytical solutions are then compared to their numerical counterparts upon discretization over unevenly spaced coordinate systems. We find that a symmetric scheme for the advective term and a three-point asymmetric scheme for the basal boundary condition best match our analytical solutions. A further convergence study shows that n≥15 vertical points are sufficient to accurately reproduce the temperature profile. The solutions presented herein are general and fully applicable to any problem with an equivalent set of boundary conditions and any given initial temperature distribution.