  We present a novel 2D flux density model for observed HI emission lines
combined with a Bayesian stacking technique to measure the baryonic
Tully-Fisher relation below the nominal detection threshold. We simulate a
galaxy catalogue, which includes HI lines described either with Gaussian or
busy function profiles, and HI data cubes with a range of noise and survey
areas similar to the MeerKAT International Giga-Hertz Tiered Extragalactic
Exploration (MIGHTEE) survey. With prior knowledge of redshifts, stellar masses
and inclinations of spiral galaxies, we find that our model can reconstruct the
input baryonic Tully-Fisher parameters (slope and zero point) most accurately
in a relatively broad redshift range from the local Universe to $z = 0.3$ for
all the considered levels of noise and survey areas, and up to $z = 0.55$ for a
nominal noise of $90\,\mu$Jy/channel over 5 deg$^{2}$. Our model can also
determine the $M_{\rm HI} - M_{\star}$ relation for spiral galaxies beyond the
local Universe, and account for the detailed shape of the HI emission line,
which is crucial for understanding the dynamics of spiral galaxies. Thus, we
have developed a Bayesian stacking technique for measuring the baryonic
Tully-Fisher relation for galaxies at low stellar and/or HI masses and/or those
at high redshift, where the direct detection of HI requires prohibitive
exposure times.
