The Krafla area in north Iceland hosts a high-temperature geothermal system within a volcanic caldera. Temperature measurements from boreholes drilled for power generation reveal enigmatic contrasts throughout the drilled area. While wells in the western part of the production field indicate a 0.5-1 km thick near-isothermal (~210 °C) liquid-dominated reservoir underlain by a deeper boiling reservoir, wells in the east indicate boiling conditions extending from the surface to the maximum depth of drilled wells (~2 km). Understanding these systematic temperature contrasts in terms of the subsurface permeability structure and overall dynamics of fluid flow has remained challenging. Here, we present a new numerical model of the natural, pre-exploitation state of the Krafla system, incorporating a new geologic/conceptual model and a version of TOUGH2 extending to supercritical conditions. The model shows how the characteristic temperature distribution results from structural partitioning of the system by a rift-parallel eruptive fissure and an aquitard at the transition between deeper basement intrusions and high-permeability extrusive volcanic rocks. As model calibration is performed using a Bayesian framework, the posterior results reveal significant uncertainty in the inferred permeability values for the different rock types, often exceeding two orders of magnitude. While the model shows how zones of single-phase vapor develop above the deep intrusive heat source, more data from deep wells is needed to better constrain the extent and temperature of the deep vapor zones. However, the model suggests the presence of a significant untapped resource at Krafla.