We used the bond order-based ReaxFF potential for all of the NEMD simulations. Similar systems consisting of phosphoruscontaining additives confined between iron surfaces have previously been studied with NEMD simulations using first principles methods. However, bond order potentials are several orders of magnitude cheaper than first principles methods, allowing much larger systems to be simulated under experimentallyrelevant conditions. ReaxFF was originally developed by van Duin et al. to study hydrocarbon reactivity and has been parameterised to model a wide range of chemical systems and processes. Of particular relevance to this study, ReaxFF-based NEMD simulations are now routinely used to study tribochemical reactions. The C/H/O/Fe/P ReaxFF parameters used in this study were developed by Khajeh et al. . The point charges on the atoms vary dynamically during the simulations and were calculated using the charge equilibration (QEq) method. The parameters for Fe/O/H were originally developed by Aryanpour et. al. , Fe/C by Zou et. al. , P/O/C/H by Verlackt et. al. , and C/H/O by Chenoweth et. al. . The combined C/H/O/Fe/P ReaxFF parameters have been successfully applied to study the thermal and mechanochemical decomposition of phosphate esters with different alkyl and aryl substituents on several ferrous surfaces. The ReaxFF parameters have also been validated against first principles calculations for the adsorption energy and dissociation energy for trialkyl phosphates on α-Fe(110). The ReaxFF parameters for Fe/O/H have been shown to accurately reproduce the experimental lattice constants for α-Fe (within 1 %), Fe 3 C (within 0.5 %), and Fe 3 O 4 (within 3 %). The ReaxFF parameters have also been shown to perform favourably compared to other many-body force fields in describing the mechanical properties of α-Fe. At 300 K, the experimental elastic modulus is 166 GPa for α-Fe, 174 GPa for Fe 3 C, and 186 GPa for Fe 3 O 4 (100). With the ReaxFF parameters used here, the bulk elastic modulus at 300 K is 136 GPa for α-Fe, 148 GPa for Fe 3 C, and 129 GPa for Fe 3 O 4 . Thus, at 300 K, the experi-mental bulk elastic modulus is underpredicted by 18 % for α-Fe, 15 % for Fe 3 C, and 31 % for Fe 3 O 4 . Since previous density functional theory (DFT) calculations have highlighted the importance of contact stiffness to mechanochemical processes, we also calculated the elastic properties of thin single slabs of the three materials used in the NEMD simulations. For this calculation, the lowermost layer of atoms in each slab of material was fixed in the z-direction, while a constant normal force was applied the topmost layer of atoms to compress the system to 2 GPa. The temperature was varied between 300-900 K using the same thermostat as applied in the NEMD simulations. The corresponding elastic constant (C 33 ) was calculated from the strain fluctuations at constant normal stress. The elastic constant results for the thin slabs are shown in the Electronic Supplementary Information (ESI) (Table ). For the α-Fe(110) slabs, the elastic constant C 33 decreases from 206 GPa at 300 K to 196 GPa at 900 K. The elastic constant C 33 also decreases with increasing temperature for the Fe 3 C(010) (from 200 GPa to 181 GPa) and Fe 3 O 4 (001) (from 220 GPa to 196 GPa) slabs. The higher stiffness of thin solid slabs compared to bulk solids was also observed in previous MD simulations and was attributed to geometrical confinement. The decreasing stiffness with increasing temperature is also consistent with previous experiments. 3 Results and Discussion