We apply symmetry group methods to find the group of transformations of the dependent and independent variables that leave the thermocline equations unchanged, These transformations lead to an optimal subset of sixteen forms of similarity solution, Each form obeys an equation with one fewer dependent variable than the original thermocline equations. Previously obtained similarity solutions, which are based solely upon scaling symmetries, are special cases of just three of these forms. Two of the sixteen forms lead to linear, two-dimensional, advection-diffusion equations for the temperature, Bernoulli functional or potential vorticity. Simple exact solutions contain "internal boundary layers" that resemble the thermocline in subtropical gyres.