@essicolo For your priors on lengthscale and variance, I might suggest using zero-avoiding priors (Gamma or InvGamma) as recommended here: Robust Gaussian Processes in Stan. That should avoid some of the issues you saw with HalfCauchy. You can also take a principled approach to defining them by scaling them to the observed quantiles of your input and output data.
As for the “meaning” of the W matrix, I’d suggest watching this: Multi Output Gaussian Processes, Mauricio Alvarez - YouTube and/or reading this: http://gpss.cc/gpss17/slides/multipleOutputGPs.pdf. As far as I understand it (and I’m no expert), coregionalization uses a single latent GP and considers the different outputs to be linear combinations of different samples from this GP. The rank of W is essentially the number of distinct samples that are combined in this way, so (I think) higher rank → more samples → more idiosyncratic behavior between the outputs.
Also, I’d recommend using a different color palette for the uncertainty vs the concentration.
It’s a great notebook!