Recent experimental and theoretical work has debated whether N-heterocyclic carbenes (NHCs) are natively present in imidazolium-based ionic liquids (ILs) such as 1-ethyl-3-methylimidazolium acetate ([EMIM+][OAc-]) at room temperature. Because NHCs are powerful catalysts, determining their presence within imidazolium-based ILs is important, but experimental characterization is difficult due to the transient nature of the carbene species. Because the carbene formation reaction involves acid-base neutralization of two ions, ion solvation will largely dominate the reaction free energy and thus must be considered in any quantum chemical investigation of the reaction. To computationally study the NHC formation reaction, we develop physics-based, neural network reactive force fields to enable free energy calculations for the reaction in bulk [EMIM+][OAc-]. Our force field explicitly captures the formation of NHC and acetic acid by deprotonation of a EMIM+ molecule by acetate and in addition describes the dimerization of acetic acid and acetate. Using umbrella sampling, we compute reaction free energy profiles within the bulk IL and at the liquid/vapor interface to understand the influence of the environment on ion solvation and reaction free energies. Compared to reaction of the EMIM+/OAc- dimer in the gas phase, the bulk environment destabilizes formation of the NHC as expected due to the large ion solvation energies. Our simulations reveal a preference for the product acetic acid to share its proton with an acetate in solution and at the interface. We predict NHC content in bulk [EMIM+][OAc-] to be on the order of parts-per-million (ppm) levels, with order-of-magnitude enhancement of NHC concentration at the liquid/vapor interface. The interfacial enhancement of NHC content is due to both poorer solvation of the ionic reactants and solvophobic stabilization of the neutral NHC molecule at the liquid/vapor interface.