Data-carrying reference signals are a type of reference signal (RS) constructed on the Grassmann manifold, which allows for simultaneous data transmission and channel estimation to achieve boosted spectral efficiency at high signal-to-noise ratios (SNRs). However, they do not improve spectral efficiency at low to middle SNRs compared with conventional RSs. To address this problem, we propose a numerical optimization-based Grassmann constellation design on the Grassmann manifold that accounts for both data transmission and channel estimation. In our numerical optimization, we derive an upper bound on the normalized mean squared error (NMSE) of estimated channel matrices and a lower bound on the noncoherent average mutual information (AMI), and these bounds are optimized simultaneously by using a Bayesian optimization technique. The proposed objective function outperforms conventional design metrics in obtaining Pareto-optimal constellations for NMSE and AMI. The constellation obtained by our method achieves an NMSE comparable to conventional non-data-carrying RSs while enabling data transmission, resulting in superior AMI performance and improved spectral efficiency even at middle SNRs.