I have been in meetings lately where colleagues are interested in jointly modelling the spatial variability of over one risk factor for cancer across a geographical map. The dataset has approximately 2000 spatial observations. CARBayes (an R package) supports "multivariate CAR " models. CAR models are sort of GPs, but the kernel is only non-zero if the two observations are spatial neighbours. This allows one to attempt either integrated nested Laplace approximations (INLA) or MCMC methods. In practice, I have found using any MCMC-style algorithm other than HMC means waiting forever to get any reasonable form of convergence.
Also have talked with people using gstats (another R package) for co-kriging, which is essentially the same as multi-output GPs (I think?). They were looking at modelling multiple mineral measurements at the same mining sites. Not too sure how they found it.
