The work function is the key surface property that determines the energy required to extract an electron from the surface of a material. This property is crucial for thermionic energy conversion, band alignment in heterostructures, and electron emission devices. This work presents a high‐throughput workflow using density functional theory (DFT) to calculate the work function and cleavage energy of 33,631 slabs (58,332 work functions) that are created from 3,716 bulk materials. The number of calculated surface properties surpasses the previously largest database by a factor of ≈27. Several surfaces with an ultra‐low (<2 eV) and ultra‐high (>7 eV) work function are identified. Specifically, the (100)‐Ba‐O surface of BaMoO3 and the (001)‐F surface of Ag2F have record‐low (1.25 eV) and record‐high (9.06 eV) steady‐state work functions. Based on this database a physics‐based approach to featurize surfaces is utilized to predict the work function. The random forest model achieves a test mean absolute error (MAE) of 0.09 eV, comparable to the accuracy of DFT. This surrogate model enables rapid predictions of the work function (≈ 105 faster than DFT) across a vast chemical space and facilitates the discovery of material surfaces with extreme work functions for energy conversion and electronic device applications.