Remark that the distribution of Gaussian processes $\mathcal{X}\to\mathbb{R}$ is the extension of multivariate Gaussian for possibly infinite $\mathcal{X}$. Thus, you can use the KL divergence between the GP probability distributions by integrating over $\mathbb{R}^\mathcal{X}$ :

$$D_{KL}(P|Q)=\int_{\mathbb{R}^\mathcal{X}} \log \frac{dP}{dQ} dP\,.$$

You can use MC methods to approximate numerically this quantity over a discretized $\mathcal{X}$ by repeatedly sampling processes according to their GP distribution. I don't know if the convergence speed is sufficiently good...

Remark that the distribution of Gaussian processes $\mathcal{X}\to\mathbb{R}$ is the extension of multivariate Gaussian for possibly infinite $\mathcal{X}$. Thus, you can use the KL divergence between the GP probability distributions by integrating over $\mathbb{R}^\mathcal{X}$ :

$$D_{KL}(P|Q)=\int_{\mathbb{R}^\mathcal{X}} \log \frac{dP}{dQ} dP\,.$$

You can use MC methods to approximate numerically this quantity over a discretized $\mathcal{X}$ by repeatedly sampling processes according to their GP distribution. I don't know if the convergence speed is sufficiently good...

Remark that if $\mathcal{X}$ is finite with $|\mathcal{X}|=n$, then you fall back to the usual KL divergence for multivariate Normal distributions: $$D_{KL}\big(\mathcal{GP}(\mu_1,K_1), \mathcal{GP}(\mu_2,K_2)\big) = \frac 1 2 \Big(tr(K_2^{-1}K_1) + (\mu_2\!-\!\mu_1)^\top K_2^{-1}(\mu_2\!-\!\mu_1)-n+\log\frac{|K_2|}{|K_1|}\Big)$$