The well-known Bochner formula and refined Kato inequality (for harmonic function) tells us that for a harmonic function $u$ in $B_{2}\subset\mathbb{R}^{n}$,
$$
\frac{1}{n}|\nabla^{2}u|^{2}\leqslant |\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2}= |\nabla u|\Delta|\nabla u|\leqslant |\nabla^{2}u|^{2}.
$$
It seems that there is a certain connection between $\Delta|\nabla u|$ and $|\nabla^{2}u|^{2}$. My question is how to obtain an analogous $L^{p}$ estimate for $\Delta|\nabla u|$. Let me fully state my question below.

For any $\delta>0$, does there exist $\epsilon>0$ such that if $u$ is a harmonic function in $B_{2}\subset\mathbb{R}^{n}$ with 
$$
\begin{split}
&|\nabla u|\leqslant 1, \\
&\int_{B_{2}}||\nabla u|-1| < \epsilon\\
&\int_{B_{2}}|\nabla^{2}u|^{2p} < \epsilon
\end{split}
$$ for some integer $p\geqslant1$, then 
$$
\int_{B_{1}}(\Delta|\nabla u|)^{p}<\delta\; ?
$$
Note that we only need the estimate in a smaller ball, so one may use a cutoff function. And it is easy when $p=1$, since we can multiply a good cutoff function with bounded laplacian and then integrate by parts. But for integer $p\geqslant2$, I can’t see how to do it.