$\Delta u(x) \ge 0\iff u(x_0) \le \frac{1}{\omega_N}\int_{S_1(0)}u(x_0+ry)d\sigma(y)$ 
Let $\Omega\subset \mathbb{R}^N$ open and $u\in C^2(\Omega)$. Show
  that the proprierties below are equivalent: 
a) $\Delta u(x) \ge 0$ for
  all $x\in \Omega$ 
b) For all $x_0\in\Omega$ and all $r>0$ such that
  $\overline{B_r(x_0)}\subset \Omega$ it's true that $$u(x_0) \le
 \frac{1}{\omega_N}\int_{S_1(0)}u(x_0+ry)d\sigma(y)$$

Suggestions:

Suggestion: for $1\to 2$, study the derivative of $$g(t) =
 \frac{1}{\omega_N}\int_{S_1(0)}u(x_0+ty)d\sigma(y), 0\le t\le r$$ 
For $2\to 1$, write the taylor expansion of rder $2$ of $u$ around
  $x_0$ and use:
Let $\Omega$ be a bounded open of $\mathbb{R}^N$ and consider a sequence >$\{u_j\}$ of harmonic functions in $\Omega$, each one of them continuous >in $\overline{\Omega}$. Suppose that $$\max_{y\in \partial\Omega}|u_j(y)-u_k(y)|\le \frac{1}{j}+\frac{1}{k}, \ k,k=1,2,\cdots$$ Then $\{u_j\}$ converges uniformly in $\overline{\Omega}$ for a function $u\in C(\overline{\Omega})$, which is still harmonic in $\Omega$

So to solve $a)\implies b)$ I must study
$$g'(t) = \frac{1}{\omega_N}\int_{S_1(0)}\sum_i\partial u_i(x_0+ty)y_i d\sigma(y)$$ 
can I do that? Take the derivative inside the integral. Well, I can see that the sum of the first derivatives appear inside the integral. Maybe if I take the derivative again I arrive at the laplacian inside the integral but there's nowhere to go from there.
To solve $b)\implies a)$ I did:
$$u(x_0+ry) = u(x_0) + r\nabla u(x_0)^Ty + \frac{1}{2}r^2y^t\nabla^2 u(x_0+try)y$$
for some $t\in (0,1)$
I think I must consider two versions of this expansion and thnk about the max of the difference of their absolute sum. I know that taylor expansions have bounds for the error terms, but I don't think the difference of two expansions would have these $\frac{1}{j}+\frac{1}{k}$ bounds.
UPDATE:
$\rightarrow$
Let's try the hint gave below:
$$g'(t)=\frac{1}{\omega_N}\int_{\partial \mathbb{S}_1(0)}\nabla u(x_0+ty)\cdot yd\sigma(y) = \frac{1}{\omega_N}\int_{B_1(0)}\Delta u(x_0+ty)\ d\sigma(y)\ge 0$$
So $g'$ is crescent. $g(0) =\frac{1}{\omega_N}\int_{S_1(0)}u(x_0)\ d\sigma (y) = u(x_0)$. So for $r>0$ we have $g'(0)\le g'(r)$ which implies the result.
$\leftarrow$
Now, for the converse, by hypothesis and by the taylor expansion we have:
$$u(x_0) \le \frac{1}{\omega_N}\int_{S_1(0)} u(x_0+ry) = \\ \frac{1}{\omega_N}\int_{S_1(0)} \left(u(x_0) + r\nabla u(x_0)^Ty + \frac{1}{2}r^2y^t\nabla^2 u(x_0+try)y\right)\ d\sigma(y)$$
for some $t\in (0,1)$
I must somehow find $u_j$ and $u_k$ both harmonic in order to use that result and find another harmonic function
 A: Since it looks like you've worked out the $(a) \implies (b)$ part, I'll focus on $(b) \implies (a)$.  Suppose $u \in C^2(\Omega)$ satisfies condition (b).  Then, by Taylor's theorem,
$$u(x_0+ry) = u(x_0) + r\nabla u(x_0)^Ty + \frac{1}{2}r^2 y^T \nabla^2u(x_0)y + o(r^2).$$
Now, integrating gives
$$u(x_0) \le \frac{1}{\Omega_n} \int_{S^1} u(x_0) + r\nabla u(x_0)^Ty + \frac{1}{2}r^2 y^T \nabla^2u(x_0)y + o(r^2)\;d\sigma(y)$$
I claim that
$$\frac{1}{\Omega_n} \int_{S^1} u(x_0) + r\nabla u(x_0)^Ty + \frac{1}{2}r^2 y^T \nabla^2u(x_0)y + o(r^2)\;d\sigma(y) = u(x_0) + C\frac{N}{2}r^2 \Delta u(x_0) + o(r^2)$$
to see why, notice that by symmetry,
$$\int_{S^1} \nabla u(x_0)^Ty d\sigma(y) = \int_{S^1} \nabla u(x_0)^T(-y) d\sigma(y)$$
since the sphere is symmetric under inversions.  Similarly, by considering reflections of the coordinate axes, we can see that
$$\int_{S^1} \partial_{x_ix_j} u(x_0) y_iy_j d\sigma(y) = 0 \text{ for } i \neq j$$
so
$$\int_{S^1} y^T \nabla^2u(x_0) y \;d\sigma(y) = \Delta u(x_0) \int_{S^1} y_i^2\;d\sigma(y) = C\Delta u(x_0)$$
Thus,
$$u(x_0) \leq u(x_0) + Cr^2\Delta u(x_0) + o(r^2).$$
Now, suppose for contradiction that $\Delta u(x_0) < 0$.  Then, for $r$ sufficiently small, the $o(r^2)$ term is less than $-Cr^2 \Delta u(x_0)$, contradicting condition (b).  It follows that $\Delta u(x_0) \ge 0$ for every $x_0$, which is precisely condition (a).
