In general it’s pretty straightforward to add most new link functions or other functions - just write out the analytic expression as you would in R. The problem with the (multivariate) CDF of a bivariate normal is that it doesn’t have an analytic expression.
mvtnorm::pmvnorm uses a Monte Carlo simulation to approximate the value of the CDF.
That means it’s not differentiable, so can’t be used with HMC. It might be possible to use it with some other MCMC algorithm that doesn’t use gradient information. But in any case, to use it in greta would require a Tensorflow (or C, linked to Tensorflow) version of the Genz-Bretz algorithm they use to do the efficient Monte Carlo simulation. That’s doable, but very hard work, and you still wouldn’t be able to use HMC, so probably not worth doing in greta.
An alternative might be too recast your model to explicitly consider (rather than marginalise) the latent bivariate Gaussian. Not sure if that’s feasible though.