  We provide a construction of multiscale systems on a bounded domain $\Omega
\subset \mathbb{R}^2$ coined boundary shearlet systems, which satisfy several
properties advantageous for applications to imaging science and the numerical
analysis of partial differential equations. More precisely, we construct
boundary shearlet systems that form frames for the Sobolev spaces
$H^s(\Omega),s\in \mathbb{N} \cup \{0\},$ with controllable frame bounds and
admit optimally sparse approximations for functions, which are smooth apart
from a curve-like discontinuity. We show that the constructed systems allow
incorporating boundary conditions. Furthermore, for $s \geq 0$ and $f\in
H^s(\Omega)$ we prove that weighted $\ell^2$ norms of the $L^2-$analysis
coefficients of $f$ are equivalent to its $H^s(\Omega)$ norm. This yields in
particular, that the reweighted systems are frames also for $H^{-s}(\Omega)$.
Moreover, we demonstrate numerically, that the associated $L^2-$synthesis
operator is also stable as a map to $H^s(\Omega)$ which, in combination with
the previous result, strongly indicates that these systems constitute so-called
Gelfand frames for $(H^s(\Omega), L^2(\Omega), H^{-s}(\Omega))$.
