We study the temperature dependence of energy diffusion in two chaotic gapped quantum spin chains, a tilted-field Ising model and an XZ model, using an open-system approach. We introduce an energy imbalance by coupling the chain to thermal baths at its boundary and study the nonequilibrium steady states of the resulting Lindblad dynamics using a matrix product operator ansatz for the density matrix. We define an effective local temperature profile by comparing local reduced density matrices in the steady state with those of a uniform thermal state. We then measure the energy current for a variety of driving temperatures and extract the temperature dependence of the energy diffusion constant. For the Ising model, we are able to study temperatures well below the energy gap and find a regime of dilute excitations where rare three-body collisions control energy diffusion. A kinetic model correctly predicts the observed exponential increase of the energy diffusion constant at low temperatures. For the XZ model, we are only able to access intermediate to high temperatures relative to the energy gap, and we show that the data are well described by an expansion around the infinite temperature limit. We also discuss the limitations of the particular driving scheme and suggest that lower temperatures can be accessed using larger baths.