Let $n\ge2$, $s\in(0,1)$, and $\Omega\subset\mathbb{R}^n$ be a bounded Lipschitz domain. In this paper, we investigate the global (higher-order) Sobolev regularity of weak solutions to the fractional Dirichlet problem

\begin{equation*}\begin{cases}(-\Delta)^su=f \ \ & \text{in}\ \ \Omega,\\u=0 \ \ & \text{in}\ \ \mathbb{R}^n\setminus\Omega.\end{cases}\end{equation*}

Precisely, we prove that there exists a positive constant $\varepsilon\in(0,s]$ depending on n, s, and the Lipschitz constant of Ω such that, for any $t\in[\varepsilon,\min\{1+\varepsilon,2s\})$, when $f\in L^q(\Omega)$ with some $q\in(\frac{n}{2s-t},\infty]$, the weak solution u satisfies

\begin{equation*}\|u\|_{W^{t,p}(\mathbb{R}^n)}\le C\|f\|_{L^q(\Omega)}\end{equation*}

for all $p\in[1,\frac{1}{t-\varepsilon})$. In particular, when Ω is a bounded C1 domain or a bounded Lipschitz domain satisfying the uniform exterior ball condition, the aforementioned global regularity estimates hold with $\varepsilon=s$ and they are sharp in this case. Moreover, if Ω is a bounded $C^{1,\kappa}$ domain with $\kappa\in(0,s)$ or a bounded Lipschitz domain satisfying the uniform exterior ball condition, we further show the global BMO-Sobolev regularity estimate

\begin{equation*}\left\|(-\Delta)^{\frac{s}{2}}u\right\|_{\mathrm{BMO}(\mathbb{R}^n)}+\left\|\nabla^{s}u\right\|_{\mathrm{BMO}(\mathbb{R}^n)}\le C\|f\|_{L^q(\Omega)}\end{equation*}

for some $q\in(\frac{n}{s},\infty]$, which is sharp in the sense that the BMO norm can not be improved to the $L^\infty$ norm.