Numerical, time domain, models of broadband acoustic propagation using k-space methods will be described and applied to the problem of image reconstruction in photoacoustic imaging. k-space methods are a subset of time-stepping pseudospectral methods, which use FFTs to calculate field gradients, with an additional correction to make the solutions exact in the homogeneous case. Furthermore, they are closely related to a number of methods in the mathematical, engineering and physical sciences literature, including non-standard finite difference schemes, exponential integrators, wavenumber domain reverse time migration, and wave propagator models developed to solve the time-dependent Schroedinger equation. These connections will be discussed and used to illuminate the advantages of the k-space approach for large scale modelling of broadband acoustic waves. To illustrate the method, a fluid-based k-space model will be applied to the inverse acoustic initial value problem encountered in tomographic photoacoustic imaging. It is well known that this can be solved, in the linear case, using a time-reversed numerical model. Extensions from this basic case to absorbing and non-linear media will be explored.