  We use a Bayesian-calibrated multistage viscous hydrodynamic model to explore
deuteron yield, mean transverse momentum and flow observables in LHC Pb-Pb
collisions. We explore theoretical uncertainty in the production of deuterons,
including (i) the contribution of thermal deuterons, (ii) models for the
subsequent formation of deuterons (hadronic transport vs coalescence) and (iii)
the overall sensitivity of the results to the hydrodynamic model -- in
particular to bulk viscosity, which is often neglected in studies of deuteron
production. Using physical parameters set by a comparison to only light hadron
observables, we find good agreement with measurements of the mean transverse
momentum $\langle p_T \rangle$ and elliptic flow $v_2$ of deuterons; however,
tension is observed with experimental data for the deuteron multiplicity in
central collisions. The results are found to be sensitive to each of the
mentioned theoretical uncertainties, with a particular sensitivity to bulk
viscosity, indicating that the latter is an important ingredient for an
accurate treatment of deuteron production.
