In BOLD fMRI, stimulus related phase changes have been repeatedly observed in humans. However, virtually all fMRI processing utilizes the magnitude information only, while ignoring the phase. This results in an unnecessary loss of physiological information and signal-to-noise efficiency. A widely held view is that the BOLD phase change is zero for a voxel containing randomly orientated blood vessels and that phase changes are only due to the presence of large vessels. Based on a previously developed theoretical model, we show through simulations and experimental human BOLD fMRI data that a non-zero phase change can be present in a region with randomly oriented vessels. Using simulations of the model, we first demonstrate that a spatially distributed susceptibility results in a non-zero phase distribution. Next, experimental data in a finger-tapping experiment show consistent bipolar phase distribution across multiple subjects. This model is then used to show that in theory a bipolar phase distribution can also be produced by the model. Finally, we show that the model can produce a bipolar phase pattern consistent with that observed in the experimental data. Understanding of the mechanisms behind the experimentally observed phase changes in BOLD fMRI would be an important step forward and will enable biophysical model based methods for integrating the phase and magnitude information in BOLD fMRI experiments.