Let $u$ harmonic in $\Omega $ and $B(x,r)\subset \subset \Omega $. Then, $\partial _iu$ also harmonic, and thus by mean value property, $$\partial _iu=\frac{1}{|B(x,r)|}\int_{B(x,r)}\partial _iu\mathrm d y.$$ Using divergence theorem on $\partial _iu \boldsymbol e_i=div(u\boldsymbol e_i)$ (where $\boldsymbol e_i=(0,...,0,1,0,...,0)$ where the 1 is at the $i-$th position), we get
$$\partial _iu=\frac{1}{|B(x,r)|}\int_{\partial B(x,r)}u\nu_i,$$ where $\nu_i$ is the $i-$th component of the unit normal vector of $\partial B(x,r)$.
Then,
$$|\partial _iu|\leq \frac{1}{|B(x,r)|}\int_{\partial B(x,r)}|u\nu_i|.$$
Now, I don't understand how we get to $$|\partial _i u|\leq \frac{C}{r^{d}}r^{d-1}\|u\|_{L^1(B(x,4r)},$$ for a constant $C$.
The key thing you're missing is the preliminary estimate $$ |u(y)| \le \frac{1}{\omega_n s^n} \int_{B(y,s)} |u(z)| dz, $$ which follows directly from the mean-value property. Apply this for each $y \in \partial B(x,r)$ with $s=r$ to get rid of the dependence on $y$ in the ball: $$ |u(y)| \le \frac{1}{\omega_n s^n} \int_{B(y,r)} |u(z)| dz \le \frac{1}{\omega_n r^n} \int_{B(x,2r)} |u(z)| dz. $$ Thus $$ \sup_{y \in \partial B(x,r)} |u(y)| \le \frac{1}{\omega_n r^n} \int_{B(x,2r)} |u(z)| dz. $$
Now we combine with your last estimate: $$ |\partial_i u(x)| \le \frac{1}{\omega_n r^n} \int_{\partial B(x,r)} |u(y)| d\sigma(y) \le \frac{1}{\omega_n r^n} \int_{\partial B(x,r)} \sup_{ \partial B(x,r)} |u| d\sigma(y) = \frac{\alpha_n r^{n-1}}{\omega_n r^n} \sup_{ \partial B(x,r)} |u| \\ \le \frac{\alpha_n}{\omega_n r }\frac{1}{\omega_n r^n} \int_{B(x,2r)} |u(z)| dz = \frac{C_n}{r^{n+1}} \int_{B(x,2r)} |u(z)| dz. $$
This is not the inequality that you have stated, but it is the correct form of the derivative estimates for harmonic functions. I suspect you have a typo in your statement. In general one can prove $$ |D^k u(x)| \le \frac{C_{n,k}}{r^{n+k}} \int_{B(x,2r)} |u|, $$ i.e. every time we take a derivative we increase the power of $r$ in the denominator on the right side.