This article presents a new variational data assimilation (VDA) approach for the formal treatment of bias in both model outputs and observations. This approach relies on the Wasserstein metric, stemming from the theory of optimal mass transport, to penalize the distance between the probability histograms of the analysis state and an a priori reference dataset, which is likely to be more uncertain but less biased than both model and observations. Unlike previous bias-aware VDA approaches, the new Wasserstein metric VDA (WM-VDA) treats systematic biases of unknown magnitude and sign dynamically in both model and observations, through assimilation of the reference data in the probability domain, and can recover the probability histogram of the analysis state fully. The performance of WM-VDA is compared with the classic three-dimensional VDA (3D-Var) scheme for first-order linear dynamics and the chaotic Lorenz attractor. Under positive systematic biases in both model and observations, we consistently demonstrate a significant reduction in the forecast bias and unbiased root-mean-squared error.