  Given an increasing number of gamma-ray bursts accompanied by potential
kilonovae there is a growing importance to advance modelling of kilonova
afterglows. In this work, we investigate how the presence of two electron
populations that follow a Maxwellian (thermal) and a power-law (non-thermal)
distributions affect kilonova afterglow light curves. We employ semi-analytic
afterglow model, $\texttt{PyBlastAfterglow}$. We consider kilonova ejecta
profiles from ab-initio numerical relativity binary neutron star merger
simulations, targeted to GW170817. We do not perform model selection. We find
that the emission from thermal electrons dominates at early times. If the
interstellar medium density is high (${\simeq}0.1\,\ccm$) it adds an early time
peak to the light curve. As ejecta decelerates the spectral and temporal
indexes change in a characteristic way that, if observed, can be used to
reconstruct the ejecta velocity distribution. For the low interstellar medium
density, inferred for GRB 170817A, the emission from the non-thermal electron
population generally dominates. We also assess how kilonova afterglow light
curves change if the interstellar medium has been partially removed and
pre-accelerated by laterally expanding gamma-ray burst ejecta. We find that the
main effect is the emission suppression at early time ${\lesssim}10^{3}\,$days,
and at its maximum it reaches ${\sim}40\%$ when the fast tail of the kilonova
ejecta moves subsonically through the wake of laterally spreading gamma-ray
burst ejecta. The subsequent rebrightening, when these ejecta break through and
shocks form, is very mild (${\lesssim}10\%$), and may not be observable.
