  We use SHARK, a semi-analytic galaxy formation model, to investigate the
physical processes involved in dictating the shape, scatter and evolution of
the HI-halo mass relation at $0\leq z \leq 2$. We compare SHARK with HI
clustering and spectral stacking of the HI-halo mass relation derived from
observations finding excellent agreement with the former and a deficiency of HI
in SHARK at $M_{\rm vir}\approx 10^{12-13} M_{\odot}$ in the latter, but
otherwise great agreement below and above that mass threshold. In SHARK, we
find that the HI mass increases with the halo mass up to a critical mass of
$\approx 10^{11.8} M_{\odot}$; between $\sim 10^{11.8}-10^{13}M_{\odot}$, the
scatter in the relation increases by 0.7 dex and the HI mass decreases with the
halo mass on average; at $M_{\rm vir} \geq 10^{13} M_{\odot}$, the HI content
continues to increase with halo mass. We find that the critical halo mass of
$\approx 10^{12} M_{\odot}$ is largely set by feedback from Active Galactic
Nuclei (AGN), and the exact shape and scatter of the HI-halo mass relation
around that mass is extremely sensitive to how AGN feedback is modelled, with
other physical processes playing a less significant role. We determine the main
secondary parameters responsible for the scatter of the HI-halo mass relation,
namely the halo spin parameter at $M_{\rm vir}\leq 10^{11.8} M_{\odot}$, and
the fractional contribution from substructure to the total halo mass for
$M_{\rm vir}\geq 10^{13} M_{\odot}$. The scatter at $10^{11.8}<M_{\rm
vir}<10^{13} M_{\odot}$ is best described by the black-hole-to-stellar mass
ratio of the central galaxy, reflecting the AGN feedback relevance. We present
a numerical model to populate dark matter-only simulations with HI at $0\leq z
\leq 2$ based solely on halo parameters that are measurable in such
simulations.
