We have applied our ab initio simulation approach for the photoemission process at solid surfaces to calculate two-photon photoemission spectra from the p(2 x 2)-reconstructed Si(001) surface. In this approach, the ground-state electronic structure of the surface is obtained within density functional theory. The subsequent time-dependent simulation is carried through at frozen effective potential, while an optical potential is applied to account for inelastic scattering in the excited state. We have derived normal emission spectra for s-and p-polarized light with photon energies in the range (h) over bar omega = 3.85-4.75 eV. The dependence of the theoretical spectra on photon energy and polarization is analyzed and compared to experimental spectra from the literature. To unravel the role of the unoccupied states between Fermi energy and the vacuum level which are acting as intermediate states in the excitation process, we investigate the expression for the two-photon photocurrent from perturbation theory. The scattering states, which serve as the final states of photoemission, are obtained from a time-dependent simulation of a LEED-type experiment. The evaluation of the dipole matrix elements allows us to identify the relevant bulk band transitions and to address the influence of surface states. DOI: 10.1103/PhysRevB.86.235134