Earth system models account for seasonal snow cover, but many do not accommodate the deeper snowpack on ice sheets (aka firn) that slowly transforms to ice under accumulating snowfall. To accommodate and resolve firn depths of up to 60 m in the Energy Exascale Earth System Model's land surface model (ELM), we add 11 layers to its snowpack and evaluate three dry snow compaction equations in multi-century simulations. After comparing results from ELM simulations (forced with atmospheric reanalysis) with empirical data, we find that implementing into ELM a two-stage firn densification model produces more accurate dry firn densities at intermediate depths of 20–60 m. Compared to modeling firn using the equations in the (12 layer) Community Land Model (version 5), switching to the two-stage firn densification model (with 16 layers) significantly decreases root-mean-square errors in upper 60 m dry firn densities by an average of 41 kg m−3 (31%). Simulations with three different firn density parameterizations show that the two-stage firn densification model should be used for applications that prioritize accurate upper 60 m firn air content (FAC) in regions where the mean annual surface temperature is greater than roughly −31°C. Because snow metamorphism, firn density, and FAC are major components in modeling ice sheet surface albedo, melt water retention, and climatic mass balance, these developments advance broader efforts to simulate the response of land ice to atmospheric forcing in Earth system models.