The question is: given a Riemannian manifold with a non-degenerate metric g and an inner product $u_{\beta}u^{\beta}=1$, is $\nabla_{\mu} (u_{\alpha}u_{\beta})=0$ without demanding the trivial solution $\nabla_{\mu}u_{\alpha}=0$? There are two ways to support this.
1) Direct calculation. $\nabla_{\mu} (u_{\alpha}u_{\beta})=\partial_{\mu}(u_{\alpha}u_{\beta})-\Gamma^{\lambda}_{\mu\alpha}u_{\lambda}u_{\beta}-\Gamma^{\lambda}_{\mu\beta}u_{\lambda}u_{\alpha}$=$\partial_{\mu}(u_{\alpha}u_{\beta})-\frac{1}{2}u^{\rho}u_{\beta}(\partial_{\mu}g_{\varrho\alpha}+\partial_{\alpha}g_{\varrho\mu}-\partial_{\varrho}g_{\mu\alpha})-\frac{1}{2}u^{\rho}u_{\alpha}(\partial_{\mu}g_{\varrho\beta}+\partial_{\beta}g_{\varrho\mu}-\partial_{\varrho}g_{\mu\beta})$.
If we now evaluate this in an orthonormal basis $(e_{\alpha})$ at a point $p\in M$ with $ e_{0}=u. $ Then, $ u^{0}u_{0}=1 $, $ u^{i}u_{i}=0 $, $g_{\alpha\beta}=\delta_{\alpha\beta}$ and it follows that $\nabla_{\mu} (u_{\alpha}u_{\beta})=0$. Since it is a tensor, if it vanishes in one frame, it vanishes in all frames.
2)$\nabla_{\mu}(u_{\beta}u^{\beta})=0$ which means $ g^{\alpha\beta}\nabla_{\mu}(u_{\alpha}u_{\beta})=0 $. Since g is non-degenerate, $ \mid \nabla_{\mu}(u_{\alpha}u_{\beta})\mid=0 $ and $\nabla_{\mu}(u_{\alpha}u_{\beta})=0 $ is a solution that does not require the trivial solution $\nabla_{\mu}u_{\alpha}=0 $. The trivial solution is sufficient but not necessary. In fact, $ u_{\beta}\nabla_{\mu}u_{\alpha}+u_{\alpha}\nabla_{\mu}u_{\beta}=0 $ means $ u^{\beta}\nabla_{\mu}u_{\beta}=0 $ so in general, the vector is orthogonal to the covariant derivative of its covector.
Am I missing something? Please advise in detail if so.
There is no reason why this should be true, just think about what this maens for the flat connection on an open subset $U\subset\mathbb R^n$: In this case $u$ is just a function $u:U\to\mathbb R^n$ and the condition you impose is that $\sum_i(u_i(x))^2=1$ for all $x\in U$, where $u_i$ are the components of $u$. So this just means that $u$ has values in the unit sphere. In contrast, the equation $\nabla_\mu(u_\alpha u_\beta)=0$ would say that for each choice of $i$ and $j$ all partial derivatives of the $x\mapsto u_i(x)u_j(x)$ vanish, so all these funcitons have to be constant. So already the function $x\mapsto x/|x|$ on $\mathbb R^2\setminus\{0\}$ gives a counter example.
Both your agruments are flawed for different reason. In 1) you have to compute in local coordinates, so you cannot require $g_{\alpha\beta}=\delta_{\alpha\beta}$. And you not only need this at a point, since you have to take partial derivatives. You could compute in an orthonormal frame, but there the description of the covariant derivative takes a completely different form.
In 2) there is no interpretation of the contraction of two indices of a tensor with three indices as a norm. To get a norm, you would need two copies of the tensor and theree copies of the metric, so this would be $g^{\mu\nu}g^{\alpha\gamma}g^{\beta\rho}(\nabla_\mu (u_\alpha u_\beta))( \nabla_\nu (u_\gamma u_\rho))$.