Surface phase diagrams (SPDs) are essential for understanding the dependence of surface chemistry on reaction condition. For multi-component systems such as metal alloys, the derivation of such diagrams often relies on separate first-principles global optimization tasks under different reaction conditions. Here we show that this can be significantly accelerated by leveraging the fact that all tasks essentially share a unified configurational search space, and only a single expensive electronic structure calculation is required to evaluate the stabilities of a surface structure under all considered reaction conditions. As a general solution, we propose a Bayesian evolutionary multitasking (BEM) framework combining Bayesian statistics with evolutionary multitasking, which allows efficient mapping of SPDs even for very complex surface systems. As proofs of concept, we showcase the performance of our methods in deriving the alloy SPDs for two heterogeneous catalytic systems: the electrochemical oxygen reduction reaction (ORR) and the gas phase steam methane reforming (SMR) reaction.