It seems that the following assertion is widely accepted:

For $k\in\mathbb N$, $p\geq 2$, $\Omega \subset \mathbb R^n$ bounded with $\partial\Omega\in C^{k+2}$ and $f\in W^{k,p}(\Omega)$, the weak solution $u\in H^1_0(\Omega)$ of the problem
$$
\begin{equation}
\left\{
\begin{aligned}
-\Delta u =f \text{ in } \Omega\\
u=0 \text{ on } \partial\Omega
\end{aligned}
\right.
\end{equation}
$$
satisfies $u\in W^{k,p}(\Omega)$ and $\|u\|_{W^{k+2,p}}\leq C_{\Omega,k,p}\|f\|_{W^{k,p}}$ for some $C_{\Omega,k,p}>0$.

The above is proved in Evans using difference quotients for $p=2$. For $k=0$ it appears to be true due to an interpolation argument (Theorem 7.1 of Giaquinta's and Martinazzi's book on regularity theory). For Hölder continuous domains one can use the classical Schauder theory. But is there a reference for the complete result?