Bentonite is a fine-grained geologic material consisting mainly of montmorillonite clay. It presents a low permeability, a high swelling pressure, and a strong capacity to retain radionuclides that make it an important component in current efforts to design engineered barrier systems for the isolation of radioactive waste. In these barriers, the thermal gradient generated by radioactive decay is expected to lead to coupled thermal-hydrologic-mechanical-chemical (THMC) processes that may impact barrier performance. However, constitutive relations characterizing the THMC coupled properties of bentonite in variable temperature, aqueous chemistry, and dry density conditions remain incompletely understood. Here, we use high-performance molecular dynamics (MD) simulations to gain insight into the THMC constitutive relations of compacted montmorillonite clay. Specifically, we report large-scale MD simulations of water-saturated clay assemblages containing 27 montmorillonite particles performed using the codes GROMACS and LAMMPS (Fig. 1). Simulations were carried out using the replica-exchange MD (REMD) technique, with 96 replicas of the system with a wide range of temperatures up to 100 °C. In addition, simulated systems were progressively dehydrated to examine a range of dry densities. Results were analyzed to determine a series of properties including hydraulic conductivity, water and ion self-diffusivity, heat capacity, thermal expansion, and swelling pressure as a function of temperature, dry density, and the type of exchangeable cations (Na, K, Ca). Finally, simulation predictions were validated and refined by benchmarking against experimental results and previous MD simulation predictions. This research provides new insight into the coupled THMC properties of clay barrier systems and advances efforts to predict the performance of engineered clay barriers over a long timescale.