The Ghil-Sellers model, a diffusive one-dimensional energy balance model of Earth's climate, features -for a considerable range of the parameter descriptive of the intensity of the incoming radiation -two stable climate states (a warm one, taken to represent the present-day climate at the appropriate solar strength, and a cold one representative of snowball conditions), where the bistability results from the celebrated ice-albedo feedback. The unstable solution is obtained and characterized in this paper. We find such unstable states by applying for the first time in a geophysical context the so-called edge tracking method that has been used for studying multiple coexisting states in shear flows. This method has a great potential for studying the global instabilities in multistable systems, and for providing crucial information on the possibility of transitions when forcing is present. We examine robustness, efficiency, and accuracy properties of the edge tracking algorithm. We find that the procedure is the most efficient when taking a single bisection per cycle. Due to the strong diffusivity of the system trajectories of transient dynamics, initialized between the stable states with respect to the mean temperature, are confined to the heteroclininc trajectory, one which connects the fixed unstable and stable states, after relatively short transient times. This constraint dictates a functional relationship between observables. We characterize such a relationship between the global average temperature and a descriptor of nonequilibrium thermodynamics, the large scale temperature gradient between low and high latitudes. We find that a maximum of the temperature gradient is realized at the same value of the average temperature, about 270 K, largely independent of the strength of incoming solar radiation. Due to this maximum, a transient increase and nonmonotonic evolution of the temperature gradient is possible and not untypical. We also examine the structural properties of the system defined by bifurcation diagrams describing the equilibria depending on a system parameter of interest, here the solar strength. We construct new bifurcation diagrams in terms of quantities relevant for describing the thermodynamic disequilibrium, such as the temperature gradient and the material entropy production due to heat transport. We compare our results for the EBM to results for the intermediate complexity GCM PlaSim and find an interesting qualitative agreement.