Climate simulations are typically performed at horizontal spatial scales of hundreds of kilometers and larger, whereas snow cover can vary at spatial scales as small as a few meters. Because of the significant unresolved heterogeneity at subgrid spatial scales, the uniform gridcell snow cover implicit in the vertical one-dimensional CLM4 model structure tends to overestimate surface albedo. To better simulate surface albedo when snow is present, CLM4 employs a snow-covered fraction (SCF), parameterization. The snow depth–SCF relationship inferred from NOAA’s Snow Data Assimilation System (SNODAS) and NASA's MODIS SCF data exhibits significant differences from the snow depth–SCF relationship used in CLM4. Based on the SNODAS-MODIS data, we develop an analytical snow depth–SCF parameterization that reproduces the general features of the observed relationship, and is straightforward to implement in a land surface model. This study shows that partitioning the surface energy fluxes based on snow cover can have a significant impact on interactions between the atmosphere and the land surface, depending on when and how often conditions of incomplete snow cover occur. In the fall, increased direct exposure of the soil column to the atmosphere leads to greater emitted longwave radiation and turbulent heat fluxes, resulting in cooler soil temperatures. In the spring, a similar effect occurs, but is offset by higher absorbed solar radiation, causing the soils to warm more rapidly. Thus, the manner in which fractional snow-covered area and surface energy fluxes are calculated affects the model’s predictions of subsurface temperatures, which in turn lead to changes in the amount and vertical distribution of soil moisture. The resulting feedbacks can alter the timing of snowmelt, the partitioning of runoff and infiltration, and vegetation dynamics. By better capturing these thermal-hydrological interactions, it is expected that the model will simulate a more realistic surface and subsurface climate, and therefore positively impact the simulation of carbon and nutrient cycling. Subsequent research will incorporate these model improvements in CLM along with parameterizations for excess ground ice and thermokarst development, prognostic wetland distributions, methane emissions, and vertically resolved soil biogeochemistry. The resulting model will then be used to study and project the evolution of permafrost under future climate change as well as the potential fate of the extensive perennially frozen soil carbon stores found in high-latitude regions.