While I am preparing new paper, I found that I need to use the parabolic analogue result of the following proposition.
Proposition. Let $\Omega$ be a bounded domain in $\mathbb{R}^d$, $d\geq 2$ and $1<p<\infty$. If $f\in W^{-1}_{p}(\Omega)$, then there exists $\mathbf{F} \in L_p(\Omega)$ such that $f=\mathrm{div}\,\mathbf{F}$ in $\Omega$, i.e., $$\left<f,\varphi \right>=-\int_\Omega \mathbf{F}\cdot\nabla \varphi$$ for all $\varphi \in C_0^\infty(\Omega)$.

Proof of this fact requires Riesz representation theorem on $W^{-1}_p(\Omega)$ and Calderon-Zygmund estimates for Poisson equations, see Kim-Tsai [Lemma 3.9,KT20] A similar strategy also works for negative parabolic Sobolev spaces: recall that for $1<s,q<\infty$, we define