Divergence in Riemannian Geometry (General Relativity) I'm taking a course in General Relativity and I'm having some problems with the notation. 
I know that Einstein's tensor verifies $\nabla_aG^{ab}=0$. In physics textbooks this consequence of Bianchi identity is phrased as "the tensor has 0 divergence". I don't understand this because for me that identity means:
If you take Einstein's tensor $G$ and take the covariant derivative $\nabla_a G$ then $(\nabla_a G)^{ab}=G^{ab}_{\quad;a}=0$. I cannot see the divergence in there.
I studied that the divergence and all of those classical differential operators could be understood through the external derivative which makes sense.. 
So my question: Is $\nabla_a u^{a}$ called divergence because if you take the summation conventions it looks like a divergence? I come from a mathematical background and I think index notation is powerful but sometimes I feel like I'm missing the point.
 A: The divergence of a vector field $V = V^a \partial _a$ on a (pseudo-)Riemannian manifold is given by $\operatorname{div} V = V^a{}_{;a}$. In words, this is obtained by taking the trace of the total covariant derivative. In the special case of $\mathbb R^n$ with a flat Euclidean or pseudo-Euclidean metric, this yields the usual calculus formula for the divergence. 
By extension, it is common to define the divergence of an arbitrary tensor field as the trace of its total covariant derivative on (usually) the last two indices. So if $G$ is the (contravariant) Einstein tensor, then its divergence would be the vector field $\operatorname{div} G = G^{ba}{}_{;a} \partial_b$.  Because $G$ is symmetric, this is also equal to $G^{ab}{}_{;a}\partial _b$.
You can also apply this to the covariant Einstein tensor with components $G_{ab}$; its divergence is the $1$-form $G_{ba;}{}^{a}dx^b$.
A: Since you said you have a mathematical background and know about the unification in terms of the exterior derivative, I thought I should add some information about the connection to the exterior derivative here.
As it is often the case, the generalisation of concepts can be performed in various ways, because the aspects they rely on fall together in special cases. For the divergence of vector fields $V$ we have two coinciding definitions, namely as the trace of $\nabla V$, which is locally $\nabla_a V^a$ and motivated by the divergence on $\mathbb{R}^n$, and via the canoncial identification with 1-forms, namely $\operatorname{div} V=*d*V^\flat$. Here $V^\flat=g(V,\cdot)$ and $*$ is the Hodge dual.
Differential forms are only a subset of all tensor fields, so it is natural to use the first definition to define the divergence of a $(p,q)$-tensor $T$ by taking a trace of $\nabla T$. How this is precisely done depends on the author. In Riemannian geometry you sometimes encounter the definition
$$
\operatorname{div}T(Y_1,\dots,Y_{q-1}) = \operatorname{tr}(X\mapsto (\nabla_XT(\cdot,Y_1,\dots,Y_{q-1}))^\sharp)=\sum_i\nabla_{E_i}T(E_i,Y_1,\dots,Y_{q-1}).
$$
Here the dependance on forms is implicit. The last equality holds locally after choosing a local frame $\{E_i\}$. With this you can show that $\operatorname{div}(fg)=df$ for any function $f$ and that $2\operatorname{div} \operatorname{Ric}=dS$, so that
$$
\operatorname{div}\left(\operatorname{Ric}-\frac{1}{2}Sg\right)=0,
$$
which is the mentioned relation for the Einstein tensor (here $\operatorname{Ric}$ is the Ricci tensor and $S$ the Ricci scalar).
This also sometimes coincides with the GR convention, noting that the Einstein, energy-momentum and field strength tensors are naturally covariant tensors and you first have to raise indices to compute the divergence.
For $p$-forms, this tensor definition again coincides with the definition using the exterior derivative, so that we find, e.g. for the field strength tensor $F$:
$$
\operatorname{div} F = *d*F,
$$
which plays a role in the equations of motion of Einstein-Maxwell theory.
