To explore the geothermal system in the Gonghe basin located between Kunlun and Qinlin Orogen in northeast Tibetan Plateau, 310 Magnetotelluric (MT) sites were deployed in the basin. The MT sites have an average spacing of 1 similar to 2 km and almost cover the whole basin with an area of similar to 100 km by 100 km. We first analyzed the MT data using phase tensor analysis to qualitatively estimate the dimensionality of the subsurface resistivity structure and the resistivity anomalies. We then determined a 3D resistivity model from the surface to the depth of similar to 50 km through a 3D MT inversion of MT data at the 272 sites with relatively higher quality. Overall, the inverted resistivity model fits the MT data very well. Our 3D resistivity model reveals significantly conductive anomalies at depths above similar to 3 km, which are interpreted as widespread hot dry rocks or hydrothermal reservoirs. Furthermore, the 3D resistivity model shows pronounced low resistivity zones in the middle crust (similar to 15-35 km). These anomalies trend a NW-SE direction, consistent with the Kunlun Fault strike. Between the near-surface conductive zone and the middle crust conductive zone, lie several vertical narrow conductive zones, which are generally associated with hidden faults identified in the basin. Synthetic tests are used to validate the resolution of large low resistivity anomalies in the middle crust and the vertically elongated low resistivity anomalies in the model. These faults may serve as the channels for transporting heat from the heat source in the middle crust to the shallow zone around the depth of 3 km. Our 3D MT imaging results characterize for the first time the 3D distribution of the geothermal system in the Gonghe Basin.