Local smoothness of harmonic heat flow

Question

Assume that $ (M,g) $ is a smooth closed manifold and $ \mathbb{S}^{L-1} $ is a unit sphere with dimension $ L-1 $ in $ \mathbb{R}^{L} $. Consider the equation of harmonic heat flow $$ \partial_tu-\Delta_g u=|\nabla u|^2u, $$ where $ u=(u^1,u^2,...,u^{L}) $. For the initial problem $$ \left\{\begin{matrix} \partial_tu-\Delta_g u=|\nabla u|^2u&\text{in}&M\times(0,T),\\ u(\cdot,0)=u_0&\text{on}&M. \end{matrix}\right.\quad(*) $$ A well-known result is that if $ u_0\in C^{\infty}(M) $ and $ T $ is sufficiently small, then the above initial problem has a unique smooth solution. Such result can be find in the book "The analysis of harmonic maps and their heat flows" by F. Lin and C. Wang, Page 111. The method in that book is by using the contraction mapping theorem. I want to go another way. Firstly, standard arguments of Galerkin methods imply that $ (*) $ has weak solution for any $ T>0 $, i.e. $ \partial_tu\in L^2(0,T;L^2(M)) $, $ \nabla u\in L^{\infty}(0,T;H^1(M)) $, $ $ for any $ \phi\in C_c^{\infty}(M\times(0,T)) $, we have $$ \int_0^T\int_{M}\partial_tu\phi \mathrm{d}\operatorname{Vol}_g\mathrm{d}t+\int_{0}^T\int_M\langle\nabla u,\nabla\phi\rangle_g\mathrm{d}\operatorname{Vol}_g\mathrm{d}t=\int_0^T\int_{M}|\nabla u|^2u\phi \mathrm{d}\operatorname{Vol}_g\mathrm{d}t. $$ I want to enhance the regularity of this weak solution for sufficiently small $ T $ but I do not know how to start. Can you give me some hints or references?