The geostatistical analysis of multivariate spatial data for inference as well as spatial predictions (co-kriging) ordinarily relies on modeling of the marginal and cross-covariance functions. While the former quantifies the spatial dependence within variables, the latter quantifies the spatial dependence across distinct variables. The marginal covariance functions are always symmetric; however, the cross-covariance functions often exhibit asymmetries in the real data. Here, we propose a novel approach to introduce flexible asymmetries in the cross-covariances of a particular class of multivariate covariance functions. The proposed approach involves modeling the phase spectrum to allow for asymmetric cross-covariances. We show the capability of our proposed model to recover the right cross-dependence structure and improve spatial predictions against traditionally used models through simulation studies. Additionally, we illustrate our approach on a real trivariate dataset of particulate matter concentration (PM2.5), wind speed and relative humidity. The real data example shows that our approach outperforms the traditionally used models, in terms of model fit and spatial predictions.