  This work proposes a thermophysical model for realistic surface layers on
airless small bodies (RSTPM), for the use of interpreting their multi-epoch
thermal lightcurves (e.g WISE/NEOWISE). RSTPM considers real orbital cycle,
rotation cycle, rough surface, temperature dependent thermal parameters, as
well as contributions of sunlight reflection to observations, hence being able
to produce precise temperature distribution and thermal emission of airless
small bodies regarding the variations in orbital time scales. Details of the
physics, mathematics and numerical algorithms of RSTPM are presented. When used
to interpret multi-epoch thermal lightcurves by WISE/NEOWISE, RSTPM can give
constraints on the spin orientation and surface physical properties, like mean
thermal inertia or mean size of dust grains, roughness fraction, albedo and so
on via radiometric procedure. As an application example, we apply this model to
the main-belt object (24) Themis, the largest object of the Themis family,
which is believed to be the source region of many main-belt comets. We find
multi-epoch (2010, 2014-2018) observations of Themis by WISE/NEOWISE, yielding
18 thermal lightcurves. By fitting these data with RSTPM, best-fit spin
orientation of Themis is derived to be ($\lambda=137^\circ$, $\beta=59^\circ$)
in ecliptic coordinates, the mean radius of dust grains on the surface is
estimated to be $\tilde{b}=140^{+500}_{-114}(6\sim640)~\mu$m, indicating the
surface thermal inertia to vary from $\sim3\rm~Jm^{-2}s^{-0.5}K^{-1}$ to
$\sim60\rm~Jm^{-2}s^{-0.5}K^{-1}$ due to seasonal temperature variation.
Further analysis found that thermal light curves of Themis show a weak
rotation-phase dependent feature, indicative of heterogeneous thermal
properties or imperfections of lightcurve inversion shape model.
