In this paper, we study the fundamental solution $\varGamma(t,x;\tau,\xi)$ of the parabolic operator $L_t=\partial_t-\Delta+b(t,x)\cdot\nabla$, where for every $t$, $b(t,\cdot)$ is a divergence-free vector field; and we consider the case that $b$ belongs to the Lebesgue space $L^{\ell}(0,T;L^q(\mathbb{R}^n))$. The regularity of weak solutions to the parabolic equation $L_tu=0$ depends critically on the value of the parabolic exponent $\gamma=2/\ell+n/q$. Without the divergence-free condition on $b$, the regularity of weak solutions has been established when $\gamma\leq1$, and the heat kernel estimate has been obtained as well, except for the case that $\ell=\infty$, $q=n$. The regularity of weak solutions was deemed not true for the critical case $L^{\infty}(0,T;L^n(\mathbb{R}^n))$ for a general $b$, while it is true for the divergence-free case, and a written proof can be deduced from the results in [Yu. A. Semenov, \emph{Regularity theorems for parabolic equations}, J. Funct. Anal. \textbf{231} (2006), no. 2, 375--417]. One of the results obtained in the present paper establishes the Aronson-type estimate for critical and supercritical cases and for vector fields $b$ which are divergence free. We will prove the best possible lower and upper bounds for the fundamental solution one can derive under the current approach. The significance of the divergence-free condition has entered the study of parabolic equations rather recently, mainly because of the discovery of the compensated compactness. The interest for the study of such parabolic equations comes from its connections with Leray's weak solutions of the Navier-Stokes equations and the Taylor diffusion associated with a vector field where the heat operator $L_t$ appears naturally.