We considered the following parabolic system:

$$\{\begin{array}{ll}& u_t=d_1\mathrm{\Delta }u-a(x)\cdot \mathrm{\nabla }u+f(u,v),x\in \mathrm{\Omega },t>0,\\ & v_t=d_2\mathrm{\Delta }v-b(x)\cdot \mathrm{\nabla }v+g(u,v),x\in \mathrm{\Omega },t>0,\\ & u(x,0)=u_0(x),v(x,0)=v_0(x),x\in \mathrm{\Omega },\end{array}$$

subject to Dirichlet (or Neumann) boundary conditions. Here $\mathrm{\Omega }\subset {\mathbb{R}}^N(N\ge 1)$ is a bounded smooth domain. In addition to some results on blowup and global existence of the solution, we found some more interesting results as follows: (1) There exists double thresholds for blowup and global existence of the solution. Under certain conditions, if $f(u,v)=f_1(v)g_1(u)$ and $g(u,v)=f_2(v)g_2(u)$, then the first watershed is

$${\int }_{c_1}^{+\mathrm{\infty }}\frac{du}{g_1(u)}=+\mathrm{\infty }\mathrm{a}\mathrm{n}\mathrm{d}{\int }_{c_2}^{+\mathrm{\infty }}\frac{dv}{f_2(v)}=+\mathrm{\infty },$$

and the second watershed is

$${\int }_{{\tilde{c}}_1}^{+\mathrm{\infty }}\frac{dU}{\tilde{f}({\tilde{F}}^{-1}(K\tilde{G}(U)))}=+\mathrm{\infty }\mathrm{a}\mathrm{n}\mathrm{d}{\int }_{{\tilde{c}}_2}^{+\mathrm{\infty }}\frac{dV}{\tilde{g}({\tilde{G}}^{-1}(\frac{1}{ϵ}\tilde{F}(V)))}=+\mathrm{\infty }.$$

Here $\tilde{f},\tilde{g},\tilde{F}$ and $\tilde{G}$ will be defined in Section 2.2. (2) If there exist nonnegative smooth functions $h(u)$, $l(v)$ and $H(s)$ such that

$$f(u,v)h^{\prime }(u)l(v)+g(u,v)h(u)l^{\prime }(v)=H[h(u)l(v)]\ge 0,$$

then the watershed for blowup in finite time and global existence of the solution is

$${\int }_0^{+\mathrm{\infty }}\frac{ds}{H(s)}=+\mathrm{\infty }.$$