  We prove boundary higher integrability for the (spatial) gradient of
\emph{very weak} solutions of quasilinear parabolic equations of the form $$u_t
- \text{div}\,\mathcal{A}(x,t, \nabla u)=0 \quad \text{on} \ \Omega \times
\mathbb{R},$$ where the non-linear structure $\text{div}\,\mathcal{A}(x,
t,\nabla u)$ is modelled after the $p$-Laplace operator. To this end, we prove
that the gradients satisfy a reverse H\"older inequality near the boundary. In
order to do this, we construct a suitable test function which is Lipschitz
continuous and preserves the boundary values. \emph{These results are new even
for linear parabolic equations on domains with smooth boundary and make no
assumptions on the smoothness of $\mathcal{A}(x,t,\nabla u)$}. These results
are also applicable for systems as well as higher order parabolic equations.
