Gaussian process analysis of processes with multiple outputs is limited by the fact that far fewer good classes of covariance functions exist compared with the scalar (single-output) case. To address this difficulty, we turn to covariance function models that take a form consistent in some sense with physical laws that govern the underlying simulated process. Models that incorporate such information are suitable when performing uncertainty quantification or inferences on multidimensional processes with partially known relationships among different variables, also known as co-kriging. One example is in atmospheric dynamics where pressure and wind speed are driven by geostrophic assumptions (wind ∝ ∂/∂x pressure). In this study we develop both analytical and numerical auto-covariance and cross-covariance models that are consistent with physical constraints or can incorporate automatically sensible assumptions about the process that generated the data. We also determine high-order closures, which are required for nonlinear dependencies among the observables. We use these models to study Gaussian process regression for processes with multiple outputs and latent processes (i.e., processes that are not directly observed and predicted but interrelate the output quantities). Our results demonstrate the effectiveness of the approach on both synthetic and real data sets.