  Variables contained within the global oceans can detect and reveal the
effects of the warming climate as the oceans absorb huge amounts of solar
energy. Hence, information regarding the joint spatial distribution of ocean
variables is critical for climate monitoring. In this paper, we investigate the
spatial correlation structure between ocean temperature and salinity using data
harvested from the Argo program and construct a model to capture their
bivariate spatial dependence from the surface to the ocean's interior. We
develop a flexible class of multivariate nonstationary covariance models
defined in 3-dimensional (3D) space (longitude x latitude x depth) that allows
for the variances and correlation to change along the vertical pressure
dimension. These models are able to describe the joint spatial distribution of
the two variables while incorporating the underlying vertical structure of the
ocean. We demonstrate that proposed cross-covariance models describe the
complex vertical cross-covariance structure well, while existing
cross-covariance models including bivariate Mat\'{e}rn models poorly fit
empirical cross-covariance structure. Furthermore, the results show that using
one more variable significantly enhances the prediction of the other variable
and that the estimated spatial dependence structures are consistent with the
ocean stratification.
