Using 1D and 2D cooling codes we study thermal emission from neutron stars with steady state internal heaters of various intensities and geometries (blobs or spherical layers) located at different depths in the crust. The generated heat tends to propagate radially, from the heater down to the stellar core and up to the surface; it is also emitted by neutrinos. In local regions near the heater the results are well described with the 1D code. The heater's region projects onto the stellar surface forming a hot spot. There are two heat propagation regimes. In the first, conduction outflow regime (realized at heat rates \(H_0 \lesssim 10^{20}\) erg cm\(^{-3}\) s\(^{-1}\) or temperatures \(T_\mathrm{h} \lesssim 10^9\) K in the heater) the thermal surface emission of the star depends on the heater's power and neutrino emission in the stellar core. In the second, neutrino outflow regime (\(H_0 \gtrsim 10^{20}\) erg cm\(^{-3}\) s\(^{-1}\) or \(T_\mathrm{h} \gtrsim 10^9\) K) the surface thermal emission becomes independent of heater's power and the physics of the core. The largest (a few per cent) fraction of heat power is carried to the surface if the heater is in the outer crust and the heat regime is intermediate. The results can be used for modeling young cooling neutron stars (prior to the end of internal thermal relaxation), neutron stars in X-ray transients, magnetars and high-\(B\) pulsars, as well as merging neutron stars.