A reliable description of surfaces structures in a reactive environment is crucial to understand materials' functions. We present a first-principles theory of replica-exchange grand-canonical-ensemble molecular dynamics and apply it to evaluate phase equilibria of surfaces in a reactive gas-phase environment. We identify the different surface phases and locate phase boundaries including triple and critical points. The approach is demonstrated by addressing open questions for the Si(100) surface in contact with a hydrogen atmosphere. In the range from 300 to 1000 K, we find 25 distinct thermodynamically stable surface phases, for which we also provide microscopic descriptions. Most of the identified phases, including few order-disorder phase transitions, have not yet been observed experimentally. Furthermore, we show that the dynamic Si-Si bonds forming and breaking is the driving force behind the phase transition between 3×1 and 2×1 adsorption patterns.