To gain insights into generalization error trends for a different type of dataset and ML model, we utilize Mat-Fold on our dataset of 58,332 work functions, ϕ, of surfaces (generated from 3,716 bulk crystals that have a zero band gap) calculated by density functional theory (DFT). On average, each unique bulk crystal structure has ∼ 15 derived surfaces. The dataset contains work functions of elemental (261), binary (14,623), and ternary (43,448) crystalline surfaces. The ML model trained on this dataset was based on a random forest (RF) model and a physics-motivated custom featurization of the topmost three atomic surface layers considering their electron affinities, atomic radii, ionization energy, Mendeleev number as well as structural information in the form of area packing fraction and interlayer spacing (details explained in our previous work ). The final RF model trained with 15 features has a test-MAE of 0.09 eV utilizing a random 90/10 split and 5-fold CV for hyperparameter optimization. This MAE is about 4-5 times better than the best benchmarking model and more than six times better than the random baseline. The model enabled the discovery of surfaces with extreme work functions for thermionic energy conversion and high-brightness photocathodes. Studying this dataset with MatFold is especially interesting as it significantly differs from the defect dataset in size, classes of materials, and model architecture.