We solve numerically the equations describing the transfer of heat through the lithosphere of Io by a mixture of conduction and volcanic advection as proposed by O’Reilly and Davies (O’Reilly, T.C., Davies, G.F. [1981]. Geophys. Res. Lett. 8, 313–316), removing the requirement that average material properties must be used. As expected, the dominance of advective heat transfer by volcanic eruptions means that Io’s geothermal gradient well away from volcanic centres is very small, of order 1 K km−1. This result is independent of any reasonable assumptions about the radiogenic heating rate in the lithosphere. The lithosphere temperature does not increase greatly above the surface temperature until the base of the lithosphere is approached, except in limited areas around shallow magma bodies. As a consequence, solid volatile sulphur compounds mobilized by volcanic processes and re-deposited on the surface of Io commonly remain solid until they reach great depths as they are progressively buried by ongoing activity. For current estimates of the volcanic heat transfer rate, melting of SO2 does not begin until a depth of ∼20 km and sulphur remains solid to a depth of ∼26 km in a 30 km thick lithosphere. Rising magmas can incorporate fluids from these deep sulphur compound aquifers, and we quantify the major influence that this can have on the bulk density of the magma and hence the resulting possible intrusion and eruption styles.