I think the issue you’re having is that you didn’t use the contraction $u_bu^b=-1$ along the way. Anyway, below is the full answer.
In general, if you’re given a unit timelike vector (i.e a vector $u$ such that $g(u,u)=-1$), then the orthogonal projection of a vector $\xi$ onto the subspace spanned by $u$ is $P_{\parallel}(\xi)=-g(\xi,u)u$. A similar statement holds for covectors.
In index notation, if $u_au^b=-1$, then the projection of $\xi^a$ onto $u^b$ is $-(\xi_au^a)u^b$, where index raising and lowering is done with $g_{ab}$. The covector version is of course that the projection of $\xi_a$ onto $u_b$ is $-(\xi_au^a)u_b$.
Ok, so now coming to the stress energy tensor, let’s write it as $T_{ab}=(\rho+P)u_au_b+Pg_{ab}$, and let’s define $\xi$ to be the divergence:
\begin{align}
\xi_b&=\nabla^a(T_{ab})\\
&=\nabla^a(\rho+P)\,u_au_b+(\rho+P)(\nabla^au_a)u_b+(\rho+P)u_a(\nabla^au_b)+(\nabla^aP)g_{ab}+\underbrace{P\nabla^ag_{ab}}_{=0}
\end{align}
where I have used metric compatibility to say the last term vanishes. Now, as mentioned above, the projection onto $u$ is $-(\xi_au^a)u_b$, but since $\xi_b=0$ (the stress tensor is divergence free) this (co)vector equation is equivalent to the scalar equation $\xi_bu^b=0$. We can easily evaluate $\xi_bu^b$ using the above formula and the normalization condition $u_bu^b=-1$:
\begin{align}
0&=\nabla^a(\rho+P)\,u_a\cdot (-1)+(\rho+P)(\nabla^au_a)\cdot (-1)+(\rho+P)u_a(\nabla^au_b)u^b+\underbrace{(\nabla^aP)g_{ab}u^b}_{=(\nabla^aP)u_a}\\
&=-(\nabla^a\rho)u_a-(\rho+P)(\nabla^au_a)+(\rho+P)u_a(\nabla^au_b)u^b,
\end{align}
where I have cancelled the term $(\nabla^aP)u_a$ which appears twice with opposite sign. Lastly, this third term vanishes by differentiating the normalization condition $u_bu^b=-1$ (in index-free notation, $g(u,u)=-1$ so applying $\nabla_u$ to both sides and using metric compatibility, and symmetry of $g$, we get $2g(\nabla_uu,u)=0$). Thus, we end up with the equation
\begin{align}
0&=-\xi_bu^b=(\nabla^a\rho)u_a+(\rho+P)(\nabla^au_a),
\end{align}
which is your equation (3).
To get the perpendicular component, you simply look at $\xi_b+(\xi_au^a)u_b$, which simplifies to
\begin{align}
0&=\xi_b+(\xi_au^a)u_b\\
&=\xi_b-\left[(\nabla^a\rho)u_a+(\rho+P)(\nabla^au_a)\right]u_b\\
&=(\nabla^aP)u_au_b+(\rho+P)u_a(\nabla^au_b)+(\nabla^aP)g_{ab}\\
&=(\rho+P)u_a(\nabla^au_b)+(g_{ab}+u_au_b)\nabla^aP,
\end{align}
which gives you equation (4).
