Formal adjoint of the covariant derivative

Ad 1: Yes, there is. The formula is $$\nabla^*(X^\flat \otimes u) = - \nabla_X u -\mathrm{div}(X) \cdot u,$$ as can easily seen by local computation. Here, $X$ is a vector field and $X^\flat$ is the dual one form w.r.t. the metric.

Note that unless you have a scalar product on the bundle $E$, too, the dual operator will be an operator $\Gamma^\infty(T^*M \otimes E^*) \longrightarrow \Gamma^\infty(E^*)$. The above formula holds in both cases (depending on which case you are in, $u$ will be a section of $E$ or $E^*$, and the connection on the right will be either your given connection on $E$ or the corresponding dual connection on $E^*$ (this connection is always well-defined by the formula $(\nabla u)(s) = d(u(s)) - u(\nabla(s))$).

A similar formula holds even in the case that you don't have a metric on $M$, in which case the dual operator will send sections of $TM \otimes E^* \otimes |\Lambda|$ to sections of $E^* \otimes |\Lambda|$, where $|\Lambda|$ is the density bundle.

Ad 2: Let me first say that in my opinion, one uses the word "formal" only to indicate that one doesn't bother about any functional analytic meaning of the word adjoint whatsoever, and one does not automatically want to say that there has to be any setting where there is an "actual" adjoint.

However, to give a more constructive answer: Yes, this operator is always an adjoint operator in the functional analytic sense.

First, this is true for elliptic differential operators on manifolds. On compact manifolds, they have a unique closed extension, whose adjoint operator is given on smooth functions by the formal adjoint.

In general, consider for example the space $\mathscr{D}=\mathscr{D}(M, E)$ of test sections (which is not a Fréchet space unless $M$ is compact!). It can be embedded into its dual $\mathscr{D^\prime}(M, E^*)$. Now any differential operator $A$ is a continuous operator on $\mathscr{D}$, and its adjoint $A^\prime$ is an operator on $\mathscr{D}^\prime$. However, on the elements $\mathscr{D}\subset \mathscr{D}^\prime$ $A^\prime$ acts just as the formal adjoint of $A$. In this sense, the formal adjoint is also an "actual" adjoint in some sense.

However, this is quite formal and not at all deep. It follows at once from the definitions of distributions etc.

Being a bit late for the party, here is nevertheless a small answer. In fact, there is a rather explicit way to compute adjoints of every differential operator (any order) between (smooth, compactly supported) sections of vector bundles. Of course, this is (in general) not the Hilbert space adjoint, here you need a bit more analysis ;)

The main idea is to use a symbol calculus based on a covariant derivative. I explain here the scalar version (trivial vector bundles) but the whole thing can be done in full generality as well. First, you choose a covariant derivative, say torsion-free and you choose a density on your manifold in order to have an integration measure. As you probably know, the symmetrized covariant derivative allows you to establish a $C^\infty(M)$-linear bijection between symbols, i.e. smooth functions on the cotangent bundle being polynomial in the fiber directions, and differential operators. Note that this is a real bijection, not just taking into account the leading symbol. Of course, the symbol depends on the chosen covariant derivative.

In a second step you compute once and for all the adjoint of a differential operator $D$ with symbol $f$ by zillions of integrations by parts. The funny thing is that there is a fairly simple way how the symbol of the adjoint looks like. You need two ingredients for that:

First, the covariant derivative allows you to define a horizontal lift which in turn determines a maximally indefinite pseudo Riemannian metric on the cotangent bundle (horizontal spaces are in bijection to tangent spaces at the base point, vertical spaces are in bijection to the cotangent space, thus there is a natural pairing). This metric has a Laplace operator (better: d'Alembert operator) $\Delta$ with which you can act on the symbol $f$. In the flat situation this is just \begin{equation} \Delta_{\mathrm{flat}} = \frac{\partial^2}{\partial q^i \partial p_i} \end{equation} for a Darboux chart on $T^*M$ induced by a chart on $M$. In general, there are a couple of Christoffel symbols needed to make this globally defined ;)

Second, the density $\mu$ of your integration might not be covariantly constant. In any case, it defines a one-form by \begin{equation} \alpha(X) = \frac{\nabla_X \mu}{\mu}, \end{equation} which is now be used to cook up a new differential operator on $T^*M$. You can lift $\alpha$ vertically to a vector field $F(\alpha)$ on $T^*M$, completely canonical.

Having these two ingredients, the adjoint of $D^*$ has the following symbol \begin{equation} f^* = \exp\left(\frac{1}{2i}(\Delta + F(\alpha))\right) \overline{f}. \end{equation} The prefactor in the exponential depends a bit on your conventions concerning the assignment of symbol to operator. With this formula it is typically really just a computation to get adjoints of all kind of operators. In many cases, you can chose your density to be covariantly constant, so $\alpha = 0$.

You can find all this in much detail in my book on Poisson geometry, based on some old papers on quantization of cotangent bundles in the late 90s (together with Bordemann and Neumaier).

In the case of interesting bundles, the formula is essentially the same: you only have to choose covariant derivatives for the two vector bundles in question and modify the Laplace operator accordingly. Then you can proceed in the same way. This generalization is in a paper of mine with Bordemann, Neumaier and Pflaum.