Higher integrability for Sobolev functions Let $u \in W^{1,2}(\mathbb{R}^2)$ be a given function satisfying $$\frac{1}{|B_r|}\int_{B_r(x)} |\nabla u|^2 dy \leq \frac{1}{r^{\delta}}$$ for all $r \leq 1$, $x \in \mathbb{R}^2$ and some fixed $\delta \in (0,1)$.
Hueristically, this looks like that the gradient cannot go to infinity too quickly. In this sense, does it rigorously imply that $\nabla u \in L^{2+\epsilon}_{loc}$ for some $\epsilon >0$?
 A: Since the OP asked for a discussion of features, I provide one by way of an explanation of Christian Remling's counterexample:
Holder's inequality states that
$$ | \int fg | \leq \| f\|_p \|g\|_q $$
if $p^{-1} + q^{-1} = 1$, and the dual characterization of $L^p$ states that
$$ f \in L^p \iff \exists M \forall g\in L^q, \int fg \leq M \|g\|_q $$
Your question can be interpreted as "does the dual characterization of $L^p$ work if, instead of all $g$ in $L^q$, we only look at those $g$s which are characteristic functions of balls (here we use that for $g$ the characteristic function of $B(x,r)$ in $\mathbb{R}^n$, its $L^q$ norm is $\approx r^{n/q}$).
While one may expect that the set of characteristic functions of balls form a nice "basis" for $L^q$, the problem however is that we don't see some convexity properties: knowing that the bounds hold for $g_1$ and $g_2$ only gives
$$ \int f (g_1 + g_2) \leq M (\|g_1\|_q  + \|g_2 \|_q) $$
the right hand side of which is not $\leq \|g_1 + g_2\|_q$. In Christian Remling's answer his $f$ is built out of functions
$$ \chi_{B(0, 2^{-k} + 2^{-k^2 - k})} - \chi_{B(0,2^{-k})} $$
We have for this function, in $\mathbb{R}^n$
$$ \|\chi_{B(0, 2^{-k} + 2^{-k^2 - k})} - \chi_{B(0,2^{-k})}\|_q \approx 2^{-nk/q - k^2/q} \ll 2^{-nk/q} \approx \|\chi_{B(0, 2^{-k} + 2^{-k^2 - k})}\|_q + \|\chi_{B(0,2^{-k})}\|_q $$
These functions are, in some sense, "orthogonal" to the characteristic functions of balls, and so their largeness are not detectable by integration against balls.

Edit: since there are some questions about the actual construction of the counterexample, here's an outline of how one can construct one.

*

*We will work in $\mathbb{R}$ for simplicity. You can do the same on $\mathbb{R}^2$ using balls instead of intervals.

*We will study $f = |\nabla u|^2$ to simplify notation. Throughout $f$ will be non-negative. The requirements are now
$$ \int_{I} f \leq r^{1-\delta}  \tag{req}$$
for all $r\in (0,1]$, and $I$ an interval of width $r$.

*Given an interval $I$ with with $|I| \leq 1$. Let $J$ be the interval with the same center as $I$, but width $|J| = \frac1{3^{1/(1-\delta)}} |I|^{1/(1-\delta)} \leq\frac13 |I|$. Observe that as $I$ shrinks in length, the interval $J$ shrinks faster.  Denote by
$$ f_I(x) = \begin{cases} \frac{1}{|J|^\delta} & x\in J \\ 0 & x\not\in J \end{cases} $$
One can check that $f_I$ satisfies the requirement (req), as given any interval $I'$ we have $$\int_{I'} f_I \leq |J|^{-\delta} \min(|I'|,|J|) \leq |I'|^{1-\delta}$$

*Now take $\mathscr{I}$ a collection of disjoint intervals, then $f = \sum_{I\in \mathscr{I}} f_I$ is a function that satisfies (req): Given an arbitrary interval $I'$, if it only intersects one of the $J$s from our construction, then by the previous part we have that the inequality holds. When $I'$ spans multiple $J$s, let $\mathscr{I}'$ be the set of the corresponding $I$s, we have
$$ \int_{I'} f \leq \sum_{I\in \mathscr{I}'} |J|^{1-\delta} = \sum_{I\in \mathscr{I}'} \frac13 |I| $$
But the length of $I'$ must be at least $\frac13 \sum_{I\in \mathscr{I}'} |I|$ as the $I$ in $\mathscr{I}'$ are disjoint intervals, and the $J$s are in their centers.

*The $L^1$ norm of $f_I$ is exactly $\frac13 |I|$. For $p > 1$, the $L^p$ norm of $f_I$ is
$$ \| f_I\|_{p}^p = \left(\frac13|I|\right)^{\frac{1-\delta p }{1-\delta}} $$
Note that the exponent is $< 1$.

*Finally, take $w_n$ a sequence of numbers, each $\leq \frac13$, such that $\sum w_n$ converges but $\sum (w_n)^{\alpha}$ diverges for every $\alpha < 1$ (so something like $w_n = ( n \ln(n)^2)^{-1}$). Choose a disjoint family of intervals $I_n$ such that $|I_n| = w_n$, and setting $f = \sum f_{I_n}$. This gives an example of an $L^1$ function satisfying condition (req) that is not in any $L^p$ for $p > 1$.

A: No. Consider a function $f\in L^1(\mathbb R)$, $f\ge 0$, with
$f(x) = 2^{n^2}$ on $2^{-n}<x<2^{-n}+2^{-n^2-n}$ and essentially $f=0$ otherwise.
Then $\int_{-r}^r f(x)\, dx \simeq \sum_{n\gtrsim (-\log r)} 2^{-n}\lesssim r^{\alpha}$, but $f\notin L^p$ for $p>1$.
We can also make $f$ continuous at all $x\not=0$, so $(1/2r)\int_{x-r}^{x+r} f(t)\, dt$ stays bounded at these points (it converges to $f(x)$).
A radial function of this type gives a counterexample in any dimension.
A: $\newcommand\R{\mathbb R}\newcommand\ep\epsilon\newcommand{\de}{\delta} $For $(s,t)\in\R^2$, let
\begin{equation}
    u(s,t):=\sum c_k g\Big(\frac{R-r_k}{h_k}\Big),
\end{equation}
where $g(z):=\max(0,1-|z|)$ for real $z$,
$R:=\sqrt{s^2+t^2}$,
\begin{equation}
    c_k:=h_k^{1-\de/2},\quad h_k:=k^{-3/(1-\de)}, \quad r_k:=k^2/\ln^2k, 
\end{equation}
and $\sum:=\sum_{k\ge k_0}$, where in turn $k_0$ is an integer large enough
so that $k\ge2$ and for all $k\ge k_0$ we have $0<r_k-h_k<r_k+h_k<r_{k+1}-h_{k+1}-2$.
Then
\begin{equation}
    |\nabla u(s,t)|^2=\sum \frac{c_k^2}{h_k^2}\, 1(|R-r_k|<h_k)
    =\sum h_k^{-\de}\, 1(|R-r_k|<h_k)
\end{equation}
almost everywhere (a.e.).
So, for each $x=(s,t)\in\R^2$ there is some integer $k\ge k_0$ such that for all $r\in(0,1]$
\begin{equation}
\frac1{|B_r|}\int_{B_r(x)} |\nabla u|^2 
\ll \frac1{r^2} h_k^{-\de}\,\min(h_k,r)r\le \frac1{r^\de};  
\end{equation}
here and in what follows, $A\ll B$ means $A\le CB$ for some universal real constant $C>0$. So, the condition displayed in the OP holds (up to a universal positive real constant factor, which can obviously be removed by rescaling $u$).
Also,
\begin{equation}
    \int_{\R^2} |\nabla u|^2\ll \sum  \frac{c_k^2}{h_k^2}\,r_k h_k
    =\sum h_k^{-\de}\,r_k h_k=\sum\frac1{k\ln^2k}<\infty
\end{equation}
and
\begin{equation}
    \int_{\R^2} |u|^2\ll \sum c_k^2\,r_k h_k\le\sum  \frac{c_k^2}{h_k^2}\,r_k h_k <\infty,
\end{equation}
so that $u\in W^{1,2}(\R^2)$.
Thus, all the conditions on $u$ hold. However, for any real $\ep>0$ there some real $\eta>0$ such that
\begin{equation}
    \int_{\R^2} |\nabla u|^{2+\ep}\asymp \sum  \frac{c_k^{2+\ep}}{h_k^{2+\ep}}\,r_k h_k
    =\sum h_k^{-(1+\ep/2)\de}\,r_k h_k
    =\sum r_k h_k^{(1-\de)(1-\eta)}=
    \sum\frac1{k^{1-3\eta}\ln^2k}=\infty,
\end{equation}
so that $\nabla u\notin L^{2+\ep}$. (This answers the original version of the question, before the replacement of $L^{2+\ep}$ by $L^{2+\ep}_{loc}$.)
