  We characterize the model spaces $K_\Theta$ in which functions with smooth
boundary extensions are dense. It turns out that such approximation is possible
if and only if the singular measure associated to the singular inner factor of
$\Theta$ is essentially concentrated on the so-called Beurling-Carleson sets.
However, in absence of this assumption, approximation is not possible even by
functions with lower order of regularity, such as the H\"older classes and some
analytic Sobolev spaces. In contrast to earlier results on density of functions
with continuous extensions to the boundary in $K_\Theta$ and related spaces,
our method is not based on duality arguments, and is largely constructive.
