A wide range of different LJ parameters for interactions between hematite and various fluids can be found in the literature. Since none of these have been verified for the nanoscale thermal transport, we systematically vary the LJ potential parameters for the hematite atoms, ε F e/O and σ F e/O , drawing from different sources, including; ClayFF, INTERFACE-FF, Berro et al. and Savio et al. These parameters were combined, using the geometric mean mixing rules, with the LJ parameters for the fluid, ε C/H and σ C/H , which were taken from the LOPLS-AA force field. In total, we compared eight different LJ potentials, with varying interaction strength. The Fe and O LJ and partial charge parameters used for each of the solid-liquid potentials are shown in Table . Previous comparisons have suggested that, compared to ClayFF, the stronger LJ interactions proposed by Savio et al. more accurately reproduced force-distance curves obtained from DFT for hydrocarbons adsorbed on hematite surfaces. To ensure the stability of the surfaces using the solid Morse potential, the partial charges remained constant for all of the interfacial force fields at q F e = +1.800 and q O = -1.200. This also meant that any differences in the ITR were due to changes to the LJ parameters. Note that the original INTERFACE-FF (q F e = +1.740, q O = -1.160), ClayFF (q F e = +1.575, q O = -1.050), and Berro et al. (q F e = +0.771, q O = -0.514) force fields used smaller partial charges. Tests were also performed for ClayFF using the original LJ parameters for the solid-solid interactions and partial charges. For the relatively nonpolar PAO molecules considered here, we found that the ITR was identical to that obtained with the modified partial charges within the uncertainty of the NEMD simulations. These comparisons are shown in the Supplementary Information (Figure ). The NEMD simulations to calculate the ITR were performed in LAMMPS 56 using a velocity-Verlet integrator and a timestep of 0.5 fs. This timestep was sufficiently short to ensure energy conservation during the NEMD simulations. Five independent repeats were performed for each of the eight solid-liquid potentials, with each trajectory being given its own randomly generated initial velocity seed. Each system was energy minimised, followed by a 1 ns equilibration in the NPT ensemble at 300 K with an anisotropic pressure coupling of 1 atm in the z-axis, using a Nosé-Hoover thermostat and barostat with temperature and pressure damping coefficients of 0.1 ps and 1.0 ps respectively. Finally, the systems underwent a 0.25 ns equilibration in the NVT ensemble at 300K using a Nosé-Hoover thermostat. To calculate the ITR, a heat flux was applied along the z-axis, as described in the previous sections. The hot and cold thermostats regions had a thickness of 10 Å and atoms falling within these regions were thermostatted to 325 K and 275 K respectively, using a Langevin thermostat with a temperature damping coefficient of 1.0 ps. Previous studies have highlighted the importance of appropriate thermostatting when calculating the ITR from NEMD simulations. In particular, for stochastic thermostats such as Langevin, if the coupling of the thermostat and the system is too strong (short temperature damping coefficient), then spurious temperature jumps can be created between the solid and fluid re-gions, leading to overestimated ITR. Therefore, we ensure that our choice of temperature damping coefficient was sufficiently large such that it did not influence our ITC results.