Recent advances in additive manufacturing enable opportunities for realizing composite materials and novel parts with unique structures and properties [1]. Here, graphene, a highly conductive two-dimensional material with an in-plane conductivity in the range of [3000 − 5000] W/mK, embedded in a copper matrix, is studied for its heat transfer properties. A large variability in the thermal properties of the novel graphene-copper composites has been experimentally measured, ranging from 200 to 800 W/mK [1,2]. The distribution of the graphene, interface dynamics, and defects have been hypothesized to be responsible for this range, all influenced by the production process. In order to understand and eventually optimize the thermal properties of these graphene-copper composites, fundamental interface dynamics are investigated using Non-Equilibrium Molecular Dynamics (NEMD) simulations, performed in LAMMPS [3] following state-of-the-art approaches [4]. Using these simulations, we predict the thermal interface resistance and effective conductivity due to phonon dynamics of single-layer graphene embedded in various geometries. To test the robustness of this method and investigate its sensitivity, a systematic study of the response surface of NEMD simulations in the case of a high-dimensional parameter space for thermal interface resistance is conducted, comparing various Uncertainty Quantification (UQ) strategies [5]. In this contribution, we present estimates of the thermal properties of various graphene-copper configurations obtained via NEMD simulations of graphene-copper interfaces. Additionally, we will demonstrate its robustness and parameter-sensitivity using high-dimensional UQ strategies. Using UQ, we quantify the expectation that force field parameters can induce significant uncertainty in the predicted effective conductivity, the effective anisotropy, and the interplay of length scales and the thermostat resolved in the phonon mode spectrum. We aim to obtain reliable microscale heat properties of graphene-copper composites from NEMD simulations and to identify parameters critical to the model outcome, related to key system dynamics.