Fit Gaussian Process Model when output is not a 1d

In Julia, GP(x,y,mean,kernel,logNoise), the output y must be a vector. But I'm doing a Gaussian Process Simulation and my y is a 100x100 matrix. How should I fitting them into Gaussian Process?

 x=rand(100)
 dis = [abs(i-j) for i in x, j in x]

 exp(-dis)
 σ2= 1
 g = 1
 l = Matrix(I,100,100)
 μ = zeros(100)
 Σ = (σ2*exp(-dis/g))+0.1l
 d = MvNormal(μ,Σ)
 Y = rand(d,100)
 (xtrain, xtest), (ytrain, ytest) = MLJBase.partition((x,Y),0.75,multi=true)