When dealing with complex engineering systems, obtaining many field solutions in multi-query scenarios such as design optimization, uncertainty quantification or coupled simulations can become prohibitively expensive. A possible remedy for this problem are so-called surrogate models, which provide a less accurate but computationally cheaper approximation of the field solution. One common type of surrogate model are Gaussian processes (GPs), which often possess advantages in terms of training cost and quality of approximation compared to standard neural network approaches. In previous work, we developed an approach for constructing GP surrogates of parameter dependent stiff ODE solutions using reparameterizations. The aim of this thesis is to extend that approach to the 2D case, illustrating its flexibility and enabling a larger range of applications.