In addition to other multi-query problems, global optimization of design variables in multi-physics engineering problems requires the resolution of coupled systems of equations for many sets of design variables. Such coupled systems are often resolved using the partitioned approach, whereby each individual system of equations is solved separately using its numerical scheme, and the coupling equations are enforced through fixed point iteration. This optimization procedure can be made more efficient by using, for example, Gaussian Process based surrogate models. Gaussian Processes are attractive for this purpose, since they are flexible, have a strong theoretical foundation and the variance of the Gaussian Processes can be used to estimate the uncertainty of the surrogate solution. This estimate can be used to improve the accuracy of the surrogate model in areas of the parameter space where it is inaccurate. Previous work has explored the potential of applying this procedure to each distinct system of equations separately, and calculating the distribution of the solutions of the resulting coupled system through fixed point iteration of the sample paths of the Gaussian Processes. This algorithm was successfully used in Bayesian and gradient-based optimization. Conversely, the Gaussian Process modelling framework has recently been extended to handle various types of constraints. Novel strategies have been developed to ensure that Gaussian Process sample paths be, for example, monotonic, or convex. Constraining Gaussian Process sample paths to satisfy a given PDE has also been a popular line of inquiry. The aim of the present work is to extend these results concerning structured Gaussian Processes to the multiphysics case. First of all, ways for increasing the regularity of the Gaussian Process sample paths are investigated, with the objective of ensuring that the sample paths of the Gaussian Processes satisfy certain conditions that the solutions of the PDEs are known to satisfy. Certain strategies from the physics-constrained Gaussian Process literature are extended to the multiphysics case. Finally, we present some theoretical results concerning our work, and discuss and compare our work with others in the literature.