We compare the effect of using four different formulations for the surface thermal boundary condition on a primitive equation eddy-resolving model. The first formulation is the conventional restoring boundary condition. This relaxes the surface temperature of the model to a specified ''restoring temperature'' on a timescale of 30 to 60 days. The second formulation calculates the surface heat flux interactively by coupling the ocean model to a simple atmospheric model with an effective restoring time of several hundred days. The third formulation (Rahmstorf and Willebrand, 1995) (RW0) is a simplified energy balance model without atmospheric heat transport. The fourth formulation (Rahmstorf and Willebrand, 1995) (RW1) is a linearized energy balance model with atmospheric heat transport parameterized as a diffusion term. The biggest impact is on the vertical structure of the temperature variance. Under the restoring condition the maxima in this variance always occur beneath the surface. Under the other three boundary conditions, maxima are found at the surface and/or subsurface levels, depending on geographical location and in closer agreement with observations. There is also an increase in the magnitude and eastward extension of both the eddy and mean kinetic energy at midlatitudes and in the subpolar gyre region with the use of less constraining surface boundary conditions. We suggest that the use of a conventional restoring surface boundary condition acts to suppress mesoscale variability in eddy-resolving models. The northward heat transport is also increased by using the RW0 and RW1 formulations. The main reason for the enhancement of eddy variability and northward heat transport using the RW0 and RW1 surface boundary conditions is the release of the mean state. ¿ American Geophysical Union 1995