Simulation and systematic analysis of the surfaces of amorphous materials is a challenge for computational chemistry. For example, silica has found widespread industrial use as an adsorbent and catalyst support but available models for use with periodic DFT are limited in variety and representativeness of realistic materials. Herein we present a generic approach for the systematic construction of ensembles of amorphous materials surface models with varied roughness and termination characteristics. The power of the approach is shown with silica as the representative example. By combining MD simulations and Fourier-series-based randomization, bulk amorphous silica was modeled and cleaved to produce surfaces with systematically varied roughness and surface saturation. An automated saturation procedure resulted in surface models with silanol densities typical of high-temperature activation protocols in the range 0.35–2.00 OH nm−2, in excellent agreement with the experimental data on surface chemistry of dehydroxylated silica materials.