Perhaps the most iconic feature of melting Arctic sea ice is the distinctive ponds that form on its surface. The geometrical patterns describing how melt water is distributed over the surface largely determine the solar reflectance and transmittance of the sea ice cover, which are key parameters in climate modeling and upper ocean ecology. In order to help develop a predictive theoretical approach to studying melting sea ice, and the resulting patterns of light and dark regions on its surface in particular, we look to the statistical mechanics of phase transitions and introduce a two-dimensional random field Ising model which accounts for only the most basic physics in the system. The ponds are identified as metastable states in the model, where the binary spin variable corresponds to the presence of melt water or ice on the sea ice surface. With the lattice spacing determined by snow topography data as the only measured parameter input into the model, energy minimization drives the system toward realistic pond configurations from an initially random state. The model captures the essential mechanism of pattern formation of Arctic melt ponds, with predictions that agree very closely with observed scaling of pond sizes and transition in pond fractal dimension.