The ablation hollows or "suncups" that form on the surface of snowfields in summer are a wonderful example of pattern formation in nature. Suncups reduce the albedo of the snow and set a characteristic length for interaction of wind with the snowpack. They also contain information about the properties of the snow and its ablation rate, which could be extracted if we had a more quantitative understanding of how suncups form. A mathematical model is proposed that explains the shape, size, and dynamical behavior of suncups in terms of the interaction of solar radiation with the snowpack. Using a perturbation method, we derive a nonlinear partial differential equation for the time-dependent shape of the snow surface from an approximate physical model for the interaction of solar radiation with snow. The resulting equation, which is similar to the Kuramoto-Sivashinsky equation in fluid mechanics, has solutions with characteristic length and amplitude. We find expressions for the characteristic size of suncups in terms of the spectrally averaged diffusion length of solar radiation in snow. The model correctly describes the shape of suncups, with their spatially ordered patterns of parabolic valleys and V-shaped ridges. It is also in remarkably good agreement with the observed length scales and growth rates. Depending on the relative values of the coefficients of the nonlinear terms in the differential equation, the suncup patterns can be either stationary in time or chaotic.