For the emulator, I really don’t need a Bayesian calibration, but just one perfect parameter combination.
Would it be possible to use the GP implementations in pymc, but used fixed hyperparameter values (like lengthscale), taken from your PyTorch model? That would help computationally.
I’m not familiar with glacier mass balance models, but if you did change your mind and want to pursue getting it into Aesara, there are lots of folks here who could help! I know people have embedded differential equation solvers into PyMC models.