We provide a construction of multiscale systems on a bounded domain Ω ⊂ 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 (Ω), s ∈ N ∪ {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 to incorporate boundary conditions. Furthermore, for s ≥ 0 and f ∈ H s (Ω) we prove that weighted ℓ 2 norms of the L 2 −analysis coefficients of f are equivalent to its H s (Ω) norm. This yields in particular, that the reweighted systems are frames also for H −s (Ω). Moreover, we demonstrate numerically, that the associated L 2 −synthesis operator is also stable as a map to H s (Ω) which, in combination with the previous result, strongly indicates that these systems constitute so-called Gelfand frames for (H s (Ω), L 2 (Ω), H −s (Ω)).