Whilst the diffusivity of hydrogen (H) in forsterite (Mg2SiO4) has been extensively studied, there still remain some puzzling observations. Firstly, experiments measuring ostensibly the same process have provided different results in terms of diffusion coefficients. Secondly, despite H diffusion in forsterite generally being associated with diffusion of M-vacancies charge compensated by 2H(+), a plethora of H-bearing point defects have been observed, including those associated with Si vacancies, trivalent cations and tetravalent cations in the form of so-called 'clinohumite-type' point defects. This has been tentatively associated with some form of inter-site reaction, such as one in which a M-site vacancy associated with 2H(+) reacts with tetrahedrally coordinated Ti4+. Equivalent reactions can be constructed to form all point defects mentioned above. Here, we present a series of numerical models in which these processes are simulated. In the models, the mobility of H is described using a diffusion coefficient (D*) for the hydrogenated M-site vacancies and an equilibrium constant (K) for the relevant inter-site reaction(s). Reevaluation of published data shows that the extracted D* and K values are consistent between some different datasets, even in situations where the phenomenological (chemical) diffusion coefficient (D) over tilde extracted simply using solutions of Fick's second law, did not agree. The 'true' mobility (D*) of the M-site vacancy associated with 2H(+) must be between 1 and 2 orders of magnitude greater than previously determined (D) over tilde in order to form measurable profiles of all point defects observed by vibrational spectroscopy. Density functional theory calculations of the K of each of the inter-site reactions implemented in our model show good agreement (within an order of magnitude) with those determined experimentally for the reactions forming 'clinohumite-type' defects, but considerable disagreement (similar to 3 orders of magnitude) for the defects involving trivalent cations, potentially due to assumptions related to binding of different components within individual defects. Overall, the first-order implication is that H diffusion profiles that we observe in natural and experimental samples are unlikely to be formed by simple diffusion alone. These models provide a new methodological framework for further understanding of complex ionic diffusion mechanisms in olivine and the other nominally anhydrous minerals.