ABSTRACT In this article, the Polynomial Chaos–Kriging (PCK) metamodel is proposed to surrogate the physical model for model-assisted probability of detection (MAPoD) analysis in eddy current non-destructive testing (ECNDT). Two popular non-intrusive metamodels, the Polynomial Chaos Expansion (PCE) and Kriging, are combined properly to replace the time-consuming boundary element method-based ECNDT forward solver for accelerating the uncertainty propagation within MAPoD analysis in NDT. In PCK, the global behaviour of the model is estimated by the trend function of PCE, and the remaining local changes of it are captured by Kriging. Both the accuracy and efficiency of the proposed PCK metamodel are demonstrated on 3D ECNDT MAPoD cases involving surface flaws on the conducting thick plate using the detecting coil with rectangular cross section. Numerical results show improved performance of PCK metamodel over PCE and Kriging metamodels, respectively, while maintaining at least two digits of accuracy of the PoD metrics.