I think you can still make it work, you might just need to shape your data differently. Judging by the shape of your matrices, it looks like you have 88 conditions, which differ by the 2D X coordinate, and 400 observations for each condition. What are these observations?
If these observations are something continuous like a time series or spectra, then you just need to add that as an extra dimension to your model. Let’s call that dimension t. Your X matrix becomes shape (35200, 3), where X[:,0] is X_0, X[:,1] is X_1, and X[:,2] is t, and your Y matrix becomes shape (35200,). Note you’ll also need to include this additional dimension in your covariance structure.
If these observations are realizations of some random variable, and the GP represents some latent function, just reshape your X to shape (35200,2) and your Y to shape (35200,). You end up repeating the same coordinate pair over and over again for repeated observations, but that’s fine. If your random variable likelihood is Gaussian, you can stick with a marginal GP like you have there. If your RV is non-Guassian, you have to use a latent GP, but otherwise the structure is the same and there are several examples in the docs about how to do that.