Abstract-Due to the increasing complexity of metamaterial geometric structures, direct optimisation of these designs using conventional approaches, such as Gradient-based and evolutionary algorithms, are often impractical and limited. This is in part due to the inherently high computational cost associated with running multiple expensive high-fidelity full-wave simulations, commonly required to optimise the constitutive parameters of a single metamaterial particle. In order to alleviate this issue, we propose an efficient optimisation approach which exploits the Co-Kriging methodology, such that we can successfully couple varying levels of discretisation and solver accuracy obtained from a 3d full-wave numerical solver suite. In contrast to other optimisation strategies, we investigate the improvement in efficiency of optimisation through the use of the LOLA-Voronoi, in conjunction with Expected Improvement and the embedding of a trustregion framework within our optimisation algorithm, to accelerate the convergence of Co-Kriging. Finally, the effectiveness of the outlined algorithm will be demonstrated by a quantitative evaluation of the performance of optimised planar 2D negative index of refraction structures.