Abstract This paper investigates a few mathematical aspects of computer simulation of salt gradient solar pond’s thermal behavior. The basic equation governing heat flow in the non-convective zone of solar pond is solved by finite difference approach using the Crank–Nicholsen method. Stability and convergence of the method, specifically for the case of solar pond, is examined over a wide range of depth difference (Δx) and time difference (Δt). It is observed that the mesh ratio parameter ( r = K Δ t / ρ C p Δ x 2 ) which is used to define the stability and convergence of the method does not have an absolute value, rather its value varies with Δx. While using an actual set of Δx and Δt, the stability must be tested with reference to the set being used. Few other mathematical aspects pertaining to the actual application of the method are also investigated. Also, the effect of fineness of ambient input data on long term performance of the pond is investigated. It is observed that the diurnal variation of ambient input data yields the same accuracy as the hourly variation. Different approaches of calculating the heat losses from upper convective zone are compared for long term performance of the pond. A simple method is suggested to calculate the radiation flux at a depth which results due to multiple reflections between bottom and surface of the pond. The method saves computational time when used for simulation and is also suitable for hand calculations.