The box model, originally introduced to account for the nonresonant hole burning (NHB) dielectric experiments in supercooled liquids, is compared to the measurements of the third harmonics P3 of the polarisation, reported recently in glycerol, close to the glass transition temperature Tg [C. Crauste-Thibierge, C. Brun, F. Ladieu, D. L'Hôte, G. Biroli, and J.-P. Bouchaud, Phys. Rev. Lett. 104, 165703 (2010)]. In this model, each box is a distinct dynamical relaxing entity (hereafter called dynamical heterogeneity (DH)) which follows a Debye dynamics with its own relaxation time τdh. When it is submitted to a strong electric field, the model posits that a temperature increase δTdh, depending on τdh, arises due to the dissipation of the electrical power. Each DH has thus its own temperature increase, on top of the temperature increase of the phonon bath δTph. Contrary to the 'fast' hole burning experiments where δTph is usually neglected, the P3 measurements are, from a thermal point of view, fully in a stationary regime, which means that δTph can no longer be neglected a priori. This is why the version of the box model that we study here takes δTph into account, which implies that the δTdh of the DHs are all coupled together. The value of P3, including both the 'intrinsic' contribution of each DH as well as the 'spurious' one coming from δTph, is computed within this box model and compared to the P3 measurements for glycerol, in the same range of frequencies and temperatures T. Qualitatively, we find that this version of the box model shares with experiments some nontrivial features, e.g., the existence of a peak at finite frequency in the modulus of P3 as well as its order of magnitude. Quantitatively, however, some experimental features are not accounted for by this model. We show that these differences between the model and the experiments do not come from δTph but from the 'intrinsic' contribution of the DHs. Finally, we show that the interferences between the 3ω response of the various DHs are the most important issue leading to the discrepancies between the box model prediction and the experiments. We argue that this could explain why the box model is quite successful to account for some kinds of nonlinear experiments (such as NHB) performed close to Tg, even if it does not completely account for all of them (such as the P3 measurements). This conclusion is supported by an analytical argument which helps understanding how a 'space-free' model as the box model is able to account for some of the experimental nonlinear features. [ABSTRACT FROM AUTHOR]