Let $\Sigma$ be an $n$-dimensional smooth closed manifold ($n\ge 3$) with a non-continuous metric $g\in W^{2,\frac{n}2}\cap L^{\infty}(\Sigma)$. Let $g'$ be a fixed smooth metric on $\Sigma$, there exists $\Lambda>1$ such that $\Lambda^{-1} g'_p(X,X)\le g_p(X,X)\le \Lambda g'_p(X,X)$ for a.e. $p\in \Sigma$ and any $X\in T_p(\Sigma)$. For a differential form $\alpha\in L^2(\Sigma,\Omega^k(\Sigma))$, assume we have $d^{*}\alpha=0$ and $d^*d\alpha=0$, where $d^*$ is the codifferential with respect to $g$. If $\alpha$ is a function, then it's a solution of the divergence-form elliptic PDE $\Delta_g \alpha=0$. A duality argument as in divergence-form regularity implies that $\alpha\in W^{1,p}$ for any $p<\infty$. And we can further prove that $\alpha\in W^{3,q}$ for any $q<\frac n2$.
If $\alpha$ is a $k$-form for $k\ge 1$, I wonder if we can get the same improvement of regularity. We have $\Delta_g \alpha=0$, which in local coordinates is a non-divergence-form elliptic system of equations. But it seems we can only get $\alpha\in W^{1,\frac{2n}{2+n}}$ and $\alpha\in W^{2,p}$ for any $p\ge 1$ satisfying $\frac 1p\ge \frac 12+\frac 2n$. Maybe this is really the best we can do, otherwise there should be other ways to use the condition $d^*\alpha=0$.