The thermal behavior of the (010), (110) and (111) copper surfaces is studied by molecular dynamics simulation. We have used a many-body potential based on the tight-binding model in order to describe the Cu-Cu interaction. The calculations we have performed correspond to simulations in the temperature range between 600 and 1800 K. The observed order in the stability follows the same order as in the packing density, i. e., (110), (010) and (111). The (110) disorder results from anharmonic effects and by vacancy-adatom formation. On the other end, the (111) surface is very stable, and remains so up to temperatures of the order of the bulk melting point. The melting proceeds by a layer-by-layer mechanism.