A nonlinear, low-order physics-based model for the dynamics of forced convection wall heat transfer in pulsating flow is formulated, based on the proper orthogonal decomposition technique. In a multivariate approach, proper orthogonal decomposition modes are constructed from computational fluid dynamics data for laminar flow and heat transfer over a flat plate in pulsating flow, spanning a range of pulsation frequencies and amplitudes. Then, the conservation equations for mass, momentum, and energy are projected onto the proper orthogonal decomposition modes, such that a system of ordinary differential equations for the modal amplitudes is obtained. The forcing at the inlet is written explicitly in the ordinary differential equations of the low-order model. The contribution of the nonvanishing pressure term resulting from the incompressible Navier-Stokes equation is included with a calibration method. The accuracy and stability of the low-order model are evaluated by comparison with computational fluid dynamics data. Possible applications of this heat source model to the computation of a describing function or the prediction of limit cycle amplitudes of thermoacoustic instabilities are discussed.