This study addresses a key scientific issue that remains unresolved in the past three decades for recovering geothermal energy from flooded mines. This issue is that no scientific explanation is available for the layering phenomenon of both temperature and salinity in large bodies of subterranean water (e.g., mine water), i.e., thermohaline stratification, which is commonly observed in mine water. Such a layering phenomenon, however, is very significant to the geothermal application by determining the temperature distribution and consequently the energy reserve and efficiency. For this purpose, multiphysics simulation with unique non-isothermal and non-isosolutal hydrodynamics is adopted to predict the formation and evolution of thermohaline stratifications in large bodies of mine water. The multiphysics simulation, for the first time, succeeded in reproducing the formation and evolution of thermohaline stratifications with a theory assuming lateral double-diffusive intrusions to mine water. The simulation results revealed that the evolution of thermohaline stratifications involves the layer-merging event, in which several small layers gradually merge to form layers with a larger thickness. The results also indicated that the buoyancy ratio is a key parameter for producing clear thermohaline stratifications in large bodies of mine water and its critical value was suggested to be 1.0. To successfully reproduce thermohaline stratifications, the required condition was concluded to be the lateral heat flux with a difference between the lateral heat fluxes, while the lateral salinity flux was not required. It is the first time, to the best of our knowledge, that the layering phenomenon in large-scale subterranean water bodies has been successfully reproduced and explained scientifically. This study will provide a solid scientific basis for the efficient and sustainable use of large bodies of subterranean water in flooded mines for geothermal energy recovery.