The thermodynamic nonequilibrium (TNE) effects in a coalescence process of two initially static bubbles under thermal conditions are investigated by a discrete Boltzmann model. The spatial distributions of the typical nonequilibrium quantity, i.e., nonorganized momentum fluxes (NOMFs), during evolutions are investigated in detail. The density-weighted statistical method is used to highlight the relationship between the TNE effects and the morphological and kinetics characteristics of bubble coalescence. The results show that the xx component and yy component of NOMFs are antisymmetrical; the xy component changes from an antisymmetric internal and external double quadrupole structure to an outer octupole structure during the coalescence process. Moreover, the evolution of the averaged xx component of NOMFs provides two characteristic instants, which divide the nonequilibrium process into three stages. The first instant, when the averaged xx component of the NOMFs reaches its first local minimum, corresponds to the moment when the mean coalescence speed gets the maximum, and at this time the ratio of minor and major axes is about 1/2. The second instant, when the averaged xx component of the NOMFs gets its second local maximum, corresponds to the moment when the ratio of minor and major axes becomes 1 for the first time. It is interesting to find that the three quantities, TNE intensity, acceleration of coalescence, and the slope of boundary length, show a high degree of correlation and attain their maxima simultaneously. The surface tension and the heat conduction accelerate the process of bubble coalescence, while the viscosity delays it. Both the surface tension and the viscosity enhance the global nonequilibrium intensity, whereas the heat conduction restrains it. These TNE features and findings present some insights into the kinetics of bubble coalescence.