  We prove boundary higher integrability for the (spatial) gradient of
\emph{very weak} solutions of quasilinear parabolic equations of the form $$
\left\{
  \begin{array}{ll} u_t - div \mathcal{A}(x,t,\nabla u) = 0 &\quad \text{on} \
\Omega \times (-T,T), \\ u = 0 &\quad \text{on} \ \partial \Omega \times
(-T,T),
  \end{array} \right. $$ where the non-linear structure $\mathcal{A}(x,
t,\nabla u)$ is modelled after the variable exponent $p(x,t)$-Laplace operator
given by $|\nabla u|^{p(x,t)-2} \nabla u$. To this end, we prove that the
gradients satisfy a reverse H\"older inequality near the boundary by
constructing a suitable test function which is Lipschitz continuous and
preserves the boundary values. In the interior case, such a result was proved
in \cite{bogelein2014very} provided $p(x,t) \geq \mathfrak{p}^- \geq 2$ holds
and was then extended to the singular case $\frac{2n}{n+2}< \mathfrak{p}^-\leq
p(x,t)\leq \mathfrak{p}^+ \leq 2$ in \cite{li2017very}. This restriction was
necessary because the intrinsic scalings for quasilinear parabolic problems are
different in the case $\mathfrak{p}^+ \leq 2$ and $\mathfrak{p}^-\geq 2$.
  In this paper, we develop a new unified intrinsic scaling, using which, we
are able to extend the results of \cite{bogelein2014very,li2017very} to the
full range $\frac{2n}{n+2} < \mathfrak{p}^- \leq p(x,t)\leq
\mathfrak{p}^+<\infty$ and also obtain analogous results upto the boundary.
\emph{The main novelty of this paper is that our methods are able to handle
both the singular case and degenerate case simultaneously.} To simplify the
exposition, we will only prove the higher integrability result near the
boundary, provided the domain $\Omega$ satisfies a uniform measure density
condition. Our techniques are also applicable to higher order equations as well
as systems.
