  Aims. We aim to investigate whether a multilayer ice model can be as
successful as a bulk ice model in reproducing the observed abundances of
various deuterated gas-phase species toward starless cores. Methods. We
calculate abundances for various deuterated species as functions of time
adopting fixed physical conditions. We also estimate abundance gradients by
adopting a modified Bonnor-Ebert sphere as a core model. In the multilayer ice
scenario, we consider desorption from one or several monolayers on the surface.
Results. We find that the multilayer model predicts abundances of $\rm DCO^+$
and $\rm N_2D^+$ that are about an order of magnitude lower than observed,
caused by the trapping of CO and $\rm N_2$ into the grain mantle. As a result
of the mantle trapping, deuteration efficiency in the gas phase increases and
we find stronger deuterium fractionation in ammonia than what has been
observed. Another distinguishing feature of the multilayer model is that $\rm
D_3^+$ becomes the main deuterated ion at high density. The bulk ice model is
generally easily reconciled with observations. Conclusions. Our results
underline that more theoretical and experimental work is needed to understand
the composition and morphology of interstellar ices, and the desorption
processes that can act on them. With the current constraints, the bulk ice
model appears to be better in reproducing observations than the multilayer ice
model. According to our results, the $\rm H_2D^+$ to $\rm N_2D^+$ abundance
ratio is higher than 100 in the multilayer model, while only a few $\times$ 10
in the bulk model, and so observations of this ratio could provide information
on the ice morphology in starless cores. Observations of the abundance of $\rm
D_3^+$ compared to $\rm H_2D^+$ and $\rm D_2H^+$ would provide additional
constraints for the models.
