The magnetotelluric (MT) method is important for exploration of geothermal systems. The information on the Earth's resistivity obtained with MT methods has been valuable in imaging the hydrothermal alteration of such systems. Given its ability to recover complex resistivity models for the Earth, three-dimensional (3D) MT inversion has become a common practice in geothermal exploration. However, 3D inversion is a time-consuming a nd complicated procedure that relies on computer algorithms to search for a model that can explain the measured data to a sufficient level. Furthermore, many elements of inversion require input from the practitioner, which can easily bias the results. Consequently, final 3D MT results depend on various factors, including the inversion code, the model mesh used to represent the Earth, data quality and processing, and constraints imposed during the inversion procedure. In this paper, to explore how this variability in 3D MT modeling impacts the final model, we invert MT data sets from the Krafla and Hengill geothermal areas in Iceland, using three different inversion codes. In each case, the modelers had the freedom to select a subset of the data and implement the inversion for the respective code in an optimized way. We compare the results from all the inversion codes, as well as consider the setup and assumptions made during the inversion process, all of which helps enhance the robustness and quality of the results. The comparison is done in multiple ways, using visual comparison of the recovered resistivity models, as well as comparing the structural similarities of the models by employing a structural correlation metric based on cross-gradients and other types of metrics for structural correlation. This approach highlights structures that are common in all three models, and implies that these structures are independent of the inversion code and necessary to fit the data. All modeling results from both Krafla and Hengill are consistent to first order, recovering a conductive layer on top of a resistive core typical of high temperature geothermal systems. For Hengill, the models show strong structural agreement, with all inversions recovering a moderately layered resistivity model but adding detail to previous work done in the area. Major differences are found in areas with coarse data coverage and hence questionable model resolution. Where the recovered structures in different models coincide, our confidence that these structures are well-constrained by the data is elevated, in spite of the different setup and assumptions in the codes these structures are required; so they can be interpreted in terms of geology with more certainty. Results from Krafla are not as consistent as results for Hengill, related in part to the Krafla data being nosier than the Hengill data. The models from Krafla have coinciding larger structures, but small-scale structures there are less coherent. One of the consistent structures in all the models is a conductive zone reaching from a depth of 5 km to shallower depths in the northern part of the area. (C) 2015 Elsevier Ltd. All rights reserved.