  The prevalence of spatially referenced multivariate data has impelled
researchers to develop a procedure for the joint modeling of multiple spatial
processes. This ordinarily involves modeling marginal and cross-process
dependence for any arbitrary pair of locations using a multivariate spatial
covariance function. However, building a flexible multivariate spatial
covariance function that is nonnegative definite is challenging. Here, we
propose a semiparametric approach for multivariate spatial covariance function
estimation with approximate Mat\'ern marginals and highly flexible
cross-covariance functions via their spectral representations. The flexibility
in our cross-covariance function arises due to B-spline based specification of
the underlying coherence functions, which in turn allows us to capture
non-trivial cross-spectral features. We then develop a likelihood-based
estimation procedure and perform multiple simulation studies to demonstrate the
performance of our method, especially on the coherence function estimation.
Finally, we analyze particulate matter concentrations ($\text{PM}_{2.5}$) and
wind speed data over the North-Eastern region of the United States, where we
illustrate that our proposed method outperforms the commonly used full
bivariate Mat\'ern model and the linear model of coregionalization for spatial
prediction.
