Maximum depth is crucial for many lake processes and biota, but attempts to explain its variation have achieved little predictive power. In this paper, we describe the probability distribution of maximum depths based on recent developments in the theory of fractal Brownian motions. The theoretical distribution is right-tailed and adequately captures variations in maximum depth in a dataset of 8164 lakes (maximum depths 0.1–135 m) from the northeastern United States. Maximum depth increases with surface area, but with substantial random variation—the 95% prediction interval spans more than an order of magnitude for lakes with any specific surface area. Our results explain the observed variability in lake maximum depths, capture the link between topographic characteristics and lake bathymetry, and provide a means to upscale maximum depth-dependent processes, which we illustrate by upscaling the diffusive flux of methane from northern lakes to the atmosphere.