We present an atomistic model that accounts for a range of anomalous thermodynamic properties of the fcc $\ensuremath{\delta}$ phase of Pu-Ga alloys in terms of the Invar mechanism. Two modified embedded atom method potentials are employed to represent competing electronic states in $\ensuremath{\delta}$-Pu, each of which has an individual configuration dependence as well as distinct interactions with gallium. Using classical Monte Carlo simulations, we compute the temperature dependence of various thermodynamic properties for different dilute gallium concentrations. The model reproduces the observed effects of excessive volume reduction along with a rapid shift in thermal expansion from negative to positive values with increasing gallium concentration. It also predicts progressive stiffening upon dilute-gallium alloying, while the calculated thermal softening is nearly independent of the gallium concentration in agreement with resonant ultrasound spectroscopy measurements in the literature. Analysis of the local structure predicted by the model indicates that the distribution of the gallium atoms is not completely random in the $\ensuremath{\delta}$ phase due to the presence of short-range order associated with the Invar mechanism. This effect is consistent with the nanoscale heterogeneity in local gallium concentration which is observed in recent extended x-ray absorption fine structure spectroscopy experiments. Implications of the Invar effect for phase stability and physical interpretations of the two states are also discussed.