  Let $\Omega \subset \mathbb{R}^3$ be a Lipschitz domain, and consider a
harmonic map $v: \Omega \rightarrow \mathbb{S}^2$ with boundary data
$v|\partial\Omega = \varphi$ which minimises the Dirichlet energy. For $p\geq
2$, we show that any energy minimiser $u$ whose boundary map $\psi$ has a small
$W^{1,p}$-distance to $\varphi$ is close to $v$ in H\"{o}lder norm modulo
bi-Lipschitz homeomorphisms, provided that $v$ is the unique minimiser
attaining the boundary data. The index $p=2$ is sharp: the above stability
result fails for $p<2$ due to the constructions by Almgren--Lieb \cite{al} and
Mazowiecka--Strzelecki \cite{ms}.
