Is the Covariant Derivative a phenomenological attempt? I am trying to self study QFT and i am very confused about the covariant derivative.
When we require our theory to be invariant under local gauge transformations we kind of "guess" that we need to change our differential operator so that the term  and find e.g
$$
\displaystyle D_{\mu }:=\partial _{\mu }-iqA_{\mu } \quad.
$$
We then see that our theory describes the desired/observed behavior of nature.
But from mathematics I know if I have a vector-field living on a somehow curved manifold, and I want to know how a vector changes along a small geodesics of the manifold that I "subtract the curvature" of the manifold, since this change belongs to the manifold itself, and not to the vector-field.
So I have two understandings of two seemingly similar things and I do not manage to merge them in my physical understanding of the theory.
 A: Suppose that $M$ is a manifold you are working on and for each $x\in M$ you have a vector space $E_x$ attached to that point. We demand that all of these vector spaces are of the same dimension $k$, and that in a suitable sense, they "vary smoothly" as you vary the point $x$ they are attached to.
Suppose now that you have a map $\psi$, which assigns to each $x\in M$ and element $\psi_x\in E_x$. The mathematical buzzwords here are "vector bundle" and "section".
In general you can cover the space $M$ with open neighborhoods $U_\alpha$ ($\alpha$ is some element of an indexing set), and in each of these $U_\alpha$ neighborhoods you have a list of $k$ pointwise linearly independent elements $u_{(a),\alpha}$ such that if $U_\alpha$ and $U_\beta$ are overlapping domains, on the overlap you have $$ u_{(b),\beta}=\varphi_{\beta\alpha}{}^a_{\ b}u_{(a),\alpha}, $$ where $\varphi_{\beta\alpha}$ is an $k\times k$ invertible matrix that is also a function of the points of $M$ in the overlap $U_\alpha\cap U_\beta$.
Then one can take components of sections by $\psi=\psi^a_\alpha u_{(a),\alpha}$ (no sum on the neighborhood indices $\alpha$), then on the overlap, the components will change as $$ \psi^a_\beta=\varphi_{\alpha\beta}{}^a_{\ b}\psi^b_{\alpha}. $$
Note that $\varphi_{\alpha\beta}=\varphi_{\beta\alpha}^{-1}$ and in gauge theories it is pretty common that the $\varphi_{\alpha\beta}$ matrices do not take values in the entirety of $\mathrm{GL}(k,\mathbb R)$, but only in a subgroup $G$ of it.
Just like it is the case for tangent vectors, since the matrices $\varphi_{\alpha\beta}$ depend on the points of the manifold $M$, it is impossible to compare the values of $\psi_x$ and $\psi_y$ for $x\neq y$, therefore a rule is needed to parallel transport the field values. This rule can be given infinitesimally by demanding that it has the form of an infinitesimal basis transformation, e.g. when you move from $x$ to $x+dx$ the components $\psi^a$ change as $$ \delta\psi^a=-A_\mu{}^a_{\ b}\psi^bdx^\mu, $$ and the covector-of-matrices $A_\mu{}^a_{\ b}$ actually take value in the Lie algebra of $G$ because we want this to be an infinitesimal basis transformation.
Then the covariant derivative of $\psi$ is given in components as $$ D\psi^a=\psi^a(x+dx)-(\psi^a(x)+\delta\psi^a(x))=(\partial_\mu\psi^a(x)+A_\mu{}^a_{b}(x)\psi^b(x))dx^\mu=D_\mu\psi^adx^\mu, $$ where $$ D_\mu\psi^a=\partial_\mu\psi^a+A_\mu{}^a_{\ b}\psi^b. $$
So it is actually the same differential geometric structure that underlies the covariant derivative in gauge theory and in general relativity. The primary difference is that the general relativity covariant derivative describes the parallel transport of so-called "natural geometric objects" such as vectors and tensors, while the parallel transport in gauge theory describes the parallel transport of so-called "gauge natural objects", which are objects that cannot be derived from the manifold structure alone but have to come from additional structures defined above the manifold by hand, such as a principal fibre bundle.
