We impose a surface forcing on the 2D, Boussinesq, thermohaline equations in a rectangular domain, in the form of equatorially symmetric cosine distributions of salinity flux and temperature. This system may be seen as an idealization of the ocean thermohaline circulation on the global scale over intervals of centuries or millenia. Multiple steady states are found numerically. They reflect the competition between the opposite signs of the temperature and salinity-driven equatorially symmetric circulations. There are also pole-to-pole, equatorially asymmetric circulations. In the control space of the temperature and salinity-flux forcing amplitudes, these equilibria form two cusp catastrophes, and transitions between stable equilibria occur through several distinct bifurcations. These catastrophes can be reproduced in simple box models connecting stirred reservoirs through capillary pipes. This steady-state analysis may provide a framework for a better understanding of climatic transitions between different stable regimes of the ocean-atmosphere system.