If $\Delta u=0$ in $\Omega\subset\mathbb{R}^n (n\geq2)$, then for $p>\frac{n-2}{n-1}$, $|Du|^p$ is subharmonic.
Proof:
For $|Du|(x_0)\neq 0$, we have
\begin{align}
\Delta |Du|^p&=p(p-2)|Du|^{p-4}D_kuD_{kj}uD_iuD_{ij}u+p|Du|^{p-2}|D^2u|^2\nonumber\\
&=p|Du|^{p-2}\Big((p-2)\sum_j(\frac{D_iu}{|Du|}D_{ij}u)^2+|D^2u|^2\Big).
\end{align}
Note that
\begin{align}
\sum_j(\frac{D_iu}{|Du|}D_{ij}u)^2=\Big(\frac{Du}{|Du|}\Big)^T(D^2u)^2\frac{Du}{|Du|}
\end{align}
is a quadratic form, and $\Delta u=0$, then we can assume that $D^2u$ is diagonal and $D^2u(x_0)=diag(\lambda_1,\lambda_2,...,\lambda_n)$. It follows that
\begin{align}
\sum_i\lambda_i=0.
\end{align}
Without loss of generality, we assume $|\lambda_n|=\max_i|\lambda_i|$, then
\begin{align}
|D^2u|^2=\sum_{i=1}^n\lambda_i^2\geq\lambda_n^2+ \frac{1}{n-1}(\sum_{i=1}^{n-1}\lambda_i)^2=\frac{n}{n-1}\lambda_n^2
\end{align}
We will show that if $p>\frac{n-2}{n-1}$( In fact, if $n>2$, we may choose $p\geq \frac{n-2}{n-1}$; if $n=2$, $p=\frac{n-2}{n-1}=0$ is a trivial case), then $\Delta|Du|^p\geq 0$. In fact, let $\xi=\frac{Du}{|Du|}$, then
\begin{align}
\Delta |Du|^p&\geq p|Du|^{p-2}\left(\sum_i\lambda_i^2-\frac{n}{n-1}\sum_i\lambda_i^2\xi_i^2\right)\nonumber\\
&\geq p|Du|^{p-2}\left(\frac{n}{n-1}\lambda_n^2-\frac{n}{n-1}\sum_i\lambda_i^2\xi_i^2\right)\nonumber\\
&\geq \frac{np}{n-1}|Du|^{p-2}\left(\lambda_n^2-\sum_i\lambda_n^2\xi_i^2\right)=0.
\end{align}
This completes the proof in the case $|Du(x_0)|\neq 0$. In general case, for any $0<\epsilon\leq 1$, follows the argument as above, we can prove that $\Delta (|Du|^2+\epsilon)^{\frac{p}{2}}\geq0$ if $p>\frac{n-2}{n-1}$. Then $(|Du|^2+\epsilon)^{\frac{p}{2}}$ satisfies the mean value inequality. Note that $(|Du|^2+\epsilon)^{\frac{p}{2}}$ is locally uniformly bounded and pointwise converge to $\rightarrow |Du|^p$. By Lebesgue CCT, then $|Du|^p$ is a continuous weakly subharmonic function (in viscosity sense or generalized sense) .
If $p=1$, we have that if $\nabla u(x_0)\neq 0$, then we get the Kato's inequality (this inequality can be generalized to Riemannian manifold, and applied in proof of gradient estimate by S.T.Yau.)
\begin{align}
|D^2u|^2\geq |\nabla|\nabla u||^2.
\end{align}
It follows that $|\nabla u|$ is subharmonic (which has been used in Alt and Caffarelli's paper).
