Stream temperature is a critical ecosystem indicator and plays an important role in river ice formation and breakup. Typically, stream temperature is simulated either using statistical models or discrete Eulerian energy balance models. These physically-based models often rely upon a grid-based finite difference or finite volume approach for simulating the advection-dominated in-stream energy balance. Such methods are conditionally stable and thus introduce constraints upon time and space step (e.g., Courant and Peclet constraints). They are also not conceptually consistent with commonly applied convolution-based hydrologic routing approaches, which don’t require reach discretization. Here, a novel semi-analytical technique for simulating advection, latent and sensible heat transfer, heat generation from friction, and hyporheic exchange with groundwater is introduced. It assumes a discrete convolution (i.e., transfer function or unit hydrograph) method is used for routing flows through the reach. The energy balance is applied to discrete parcels of water that travel along the reach exchanging energy with their surroundings; the parcel-based energy balance is solved exactly. The method is unconditionally stable, runs at the native time step of the hydrological model, and requires no spatial discretization. It can also easily simulate edge cases that are exceedingly difficult using discrete methods, such as a near-infinite convective exchange coefficient. The only approximation errors are associated with the model time step used to represent the inflow time series, the linearization of the Stefan-Boltzmann equation for longwave radiation flux, and the full mixing assumed at nodes of the stream network. The stream temperature model, as implemented within the Raven hydrological modelling framework, is demonstrated at several test catchments in the North American Rocky Mountains.