Motivation for covariant derivative axioms in the context of General Relativity In General Relativity the idea of a covariant derivative on a manifold is quite important and it is usually defined by a set of axioms:

Let $M$ be a smooth manifold. A covariant derivative $\nabla$ on $M$ is a map $\nabla $ which takes a vector field $X$ and an $(r,s)$-tensors $T$ for all $r,s$ producing the $(r,s)$ tensor $\nabla_X T$ and satisfying:
  
  
*
  
*$\nabla _X f = Xf,$ when $f\in \mathcal{C}^\infty(M)$,
  
*$\nabla_X (T+S)=\nabla_X T + \nabla_X S$
  
*$\nabla_X T(\omega, Y)=(\nabla_X T)(\omega,Y)+T(\nabla_X \omega, Y)+T(\omega, \nabla_X Y)$ and similarly for all $(r,s)$-tensor
  
*$\nabla_{f X+g Y}T = f\nabla_X T+g\nabla_Y T$
  

I know that in a more general contexts, this covariant derivative can be recovered from a connection on a principal bundle.
That's not what I am talking here. What I am talking here is exactly this definition of covariant derivative, usually the one that is used in General Relativity.
The problem is: the covariant derivative is highly important in the context of General Relativity, however, the definition with these axioms is overly abstract.
Is there any way in which we can motivate this definition in the context of General Relativity? Or more generaly, is there any way to motivate this definition of a covariant derivative from the point of view of Physics?
Again, I could just accept the axioms and move forward, but since this is overly abstract, and I'm dealing with Physics and not math, I'd like to gain a little motivation and insight if possible.
 A: Sure. Let me start with the story up 'til this point (someday I'll write it down in some central spot), so that we can have examples dotted all the way through.
Basics
Points and scalar fields
So you start with a set of objects $\mathcal M,$ and we're not actually going to peek at the structure of the objects themselves (except, perhaps, for equality) so we just call them "points" to indicate that we don't care about their internal structure. Kind of like in category theory, we're going to treat them as black boxes and describe their structure by adding a set of functions: in this case, the scalar fields $\mathcal S \subseteq (\mathcal M \to \mathbb R),$ which we want to be "smooth." To get this smoothness, we reinterpret functions $\mathbb R^k \to \mathbb R$ as functions $(\mathcal M \to \mathbb R)^k\to(\mathcal M \to \mathbb R)$ by applying them "pointwise". Let me denote this dual nature with square brackets (scalar field side) and parentheses (function side), formally $$f[s_1, s_2, \dots s_k] = p\mapsto f\big(s_1(p), s_2(p), \dots s_k(p)\big),$$ where $\mapsto$ constructs a function from a symbol ($p$) and a symbolic expression it maps to ($f(\dots)$). I have not seen a good name for this dual-interpretation of these smooth functions: so, when a smooth function in $C^\infty(\mathbb R^k, \mathbb R)$ is interpreted this way, I like to call them $k$-functors because there's a sorta-cool category diagram out there. 
So we've got a set of points with another set $\mathcal S$ of smooth scalar fields defined over it, and $\mathcal S$ is closed under $k$-functors for all $k$. These actually do a ton of work up-front; $\operatorname{plus}(a, b) = a + b$ and $\operatorname{times}(a, b) = a \cdot b$ are both 2-functors and so under our axiom that $\mathcal S$ is closed under $k$-functors (they map $\mathcal S^k \to \mathcal S$) these are both allowed pointwise operations on scalar fields. Even better: define that a subset of $\mathcal M$ is closed if it is a kernel for a field in $\mathcal S$ or open if its complement is closed, and you have a natural topology: pointwise multiplication gives a union operation, pointwise addition of squares gives an intersection, and you can allow infinite intersections and finite unions without any problem. Using bump functions, you can even prove that all of the scalar fields in $\mathcal S$ are continuous maps to $\mathbb R$ on this topology. As an example of this theory-point: we can now demand that the space is connected, which in topology means "the whole space is not a union of two disjoint open sets." Working back through the definitions, we first rewrite the claim to its complement; if $A \cup B=\mathcal M$ with $A,B$ disjoint $\mathcal M - A$ being open means $B$ is closed. So it's equivalently not the union of two disjoint closed sets. So axiomatically we're saying that if $s_1 \cdot s_2 = 0$ is the zero-field (which must exist because it's a $0$-functor!) then there is some point $p$ such that $s_1(p) = s_2(p) = 0.$ And that's a nice property because these scalar fields do lack this crucial property that we're so used to, "$ab = 0$ implies either $a=0$ or $b=0.$" Scalar fields can each be zero on non-overlapping subsets to multiply together to form zero. But as long as the space is connected, at least we recover something similar.
So, for actual examples of scalar fields, on the surface of the sphere the points are in fact $\{(x, y, z) : x^2 + y^2 + z^2=1\},$ but we refuse to peek inside directly. Instead we start with scalar fields $x,y,z$ that happen to extract these components, and close over smooth functions to get the full set of scalar fields. On the flip side, say, $\theta$ (the azimuthal angle -- the polar angle I think is OK) is not a valid smooth scalar field, because it has this nasty discontinuity which takes us away from the obvious topology that we'd like to use. You can also see that "locally" this will look like $\mathbb R^2$ and will have similar open sets. We could do a similar treatment with the torus etc.
The overlapping coordinate sets.
Then we have one of our most important axioms: the statement that around any point $p$ there is an open set containing $p$ and $D$ scalar fields which can (a) be used to distinguish points in that open set, and (b) can be used to expand scalar fields, so that every scalar field on that open set can be expanded as a $D$-functor of the coordinate fields. So again, on the sphere, we can use the fields $x, y$ as our coordinates in the North or South hemispheres (which are open sets if we don't include the equator: use a bump function on $z$ to see this). Similarly we have overlapping hemispheres with respect to $y$ and $z$ which do not include their respective "equators". However even if some point is on two of these equators, we can see that it's not on the third: so each point has an open set, and two "coordinate" fields, and on that subset all of the scalar fields can be written as functions $f(x, y)$ or what have you. This means $D=2$ and the sphere is 2-dimensional. Easy peasy.
Vector fields, tensor fields.
Now we introduce the vector fields which are a set $\mathcal V \subset (\mathcal S \to \mathcal S)$ obeying the Leibniz law. Say $f_{(m)}$ is the partial derivative of $f$ (which is some function in $C^\infty(\mathbb R^k, \mathbb R)$, mind) with respect to its $m^\text{th}$ argument. This Leibniz law says that for any $k$-functor $f,$ $$V f[s_1, s_2, s_3, \dots s_k] = \sum_{m=1}^k f_{(m)}[s_1, \dots s_k] ~\cdot~ V s_m.$$ 
If this seems to come out of nowhere, keep in mind that it actually is very logically related to everything we said above. Not only does the closure axiom create these operations $\sum$ and $\cdot$, but the coordinate axiom means that now on some open set every scalar $s$ is secretly a $D$-functor $s[c_1, c_2, \dots c_D].$ Define on this subset the scalar fields $v_i = V c_i$, now you have straightforwardly that $V s = \sum_{i=1}^D v_i \cdot \partial_i s.$ So that's why these Leibniz linear maps are "vector fields"; locally they are directional derivatives of scalar fields. But, they are defined geometrically: they aren't defined by these components $v_i,$ they just happen to be representable that way locally. It's not hard to see that $U + V$ is well-defined or that $s V$ is well-defined, but there is no obvious $U \cdot V$ from the above definition. However there is a Lie bracket: $[U, V] = U \circ V - V\circ U$ must be Leibniz if both $U$ and $V$ are. (Furthermore this is not a "vector space" in the normal mathematical sense; it is a "module." This is just because the scalar fields are not a "field" in the normal mathematical sense: just like you can't divide by 0, you can't divide by a scalar field which is 0 in some places, so an axiom [the existence of multiplicative inverses for every nonzero element] fails.)
Once we have vector fields, we have covector fields (call this $\bar {\mathcal V}$), the linear maps $\mathcal V \to \mathcal S$. and then we can introduce the $[a, b]$ tensor fields as the multilinear maps from $(\mathcal V^m, \bar {\mathcal V}^n) \to \mathcal S.$ Call this $\mathcal V_m^n$ for natural numbers $m, n.$ Now there is a geometric version of Einstein notation, where we just create a lot of copies of this tensor space $\mathcal V_m^n$ and annotate it with a new letter $\mathcal T$ plus $n$ distinct upper symbols and $m$ mutually distinct lower symbols. We also annotate any residents of one of these spaces with the corresponding symbols, and we may need to specify those symbols to be in a tensor-dependent order (i.e. not all tensors are symmetric). Outer products are defined the obvious way, e.g. a map from $\mathcal T^a \times \mathcal T^b \to \mathcal T^{ab}$. As I recall we need an additional axiom that says that every tensor in, say, $\mathcal T^{abc}_{de}$ can be written as a sum of outer products of terms in $\mathcal T^a \times \mathcal T^b \times \mathcal T^c \times T_d \times T_e,$ but this is (if memory serves) apparently a consequence of paracompactness or the existence of the metric or something. The point is that each tensor is officially "any multilinear map from vectors and covectors to scalars", but is secretly a finite sum of outer products of vectors and covectors.
Anyway, the reason that this last axiom is important is that it lets you do index contractions: expand in terms of the finite sum, then you can apply one of the terms of $\bar {\mathcal V}$ to the corresponding term of $\mathcal V$ to get a scalar field. And as you can expect, we can symbolize this by repeating an index among the top and bottom vectors, to say "these are being joined." So $v^{abc}_{bd}$ lives in $\mathcal T^{ac}_d$ and has a purely geometric interpretation, there is no "implicit summation" of "components."
At this point we also have an automatic gradient operation on scalar fields; $\nabla_\bullet s = V \mapsto V s$ maps any scalar field to a covector field. We also introduce the metric tensor, a special $[0, 2]$ and $[2, 0]$ tensor which contract to the identity $[1,1]$ tensor and show a special bijection between covector fields and vector fields. 
The connection
OK, so once we have this whole story, the obvious question is whether there is a meaningful generalization of $\nabla_a$ to vectors, as it's uniquely defined for scalars. And the answer is, "Well, it's not quite so unique, but yeah, in many cases that exists." 
But we basically just start with the axioms. For example, we start from $\nabla_a v^b$ being meaningful, and then we want to generalize to $\nabla_a (k v^b)$ with the Leibniz rule, and we find that it should be $v^b \nabla_a k + k \nabla_a v^b.$ Similarly we want $\nabla_a (u^b + v^b) = \nabla_a u^b + \nabla_a v^b$ as a straightforward linearity requirement. Our definition for its action on covectors is also really straightforward; recall that the contraction $u_b u^b$ is a scalar, and we expect $\nabla_a(u_b v^b) = u_b \nabla_a v^b + v^b\nabla_a u_b.$ Since the first and second terms are already well-defined, we merely define the action of $\nabla_a$ on a covector as the difference of those two terms, and we get this equation for free. So we assume that some generalization of this form exists.
Your equations all concern this operator $\nabla_a.$ The connection is easily seen when you remember that $v^a (\nabla_a s)$ is defined as $V s$ by the geometric definition of the covector $\nabla_a.$ Your expression $\nabla_V$ is therefore equivalent to $v^a\nabla_a,$  and we're generalizing $\nabla_a$ to operate on vectors so it makes sense that then $v^a \nabla_a$ gets generalized as well. Your first axiom is just "the ungeneralized form still needs to do what the scalar gradient does, no messing with that please." Your second axiom is "this is a linear operator" and your third axiom is "this is a Leibniz operator", and your fourth is just a direct consequence of the fact that the $v^a$ premultiplier and contraction operation are also linear on $v^a,$ or in other words $\nabla_a$ maps $\mathcal T^\bullet \to \mathcal T^\bullet_a.$
An intuitive understanding for the degeneracy
The basic reason that this is not unique in general, is not too hard to understand either. Parallel transport of a scalar makes sense; if you're going in the direction of the gradient it's increasing, in the opposite direction it's decreasing, and it's just a number at the end of the day, so you can believe that you always get to the same number no matter how you walk. But parallel transport of a vector is harder. Let's say that I am in Kansas City in the US and I face North and stretch out my right arm as a vector pointing East. I now walk to the North pole, I'm pointing out South (of course I am, all directions are South from the North Pole), roughly towards Madrid. But suppose that I first side-step East, I should run more or less into Washington, D.C. : now if I walk North to the pole I will be pointing instead at Rome. The path you take matters, and you can roughly predict that it involves 3 tensor indices; there's something there about "you're taking as input a vector field, and a direction (which is also a vector field) and giving as output a new vector field" that seems to relate 3 different vector fields, 2 as input and 1 as output. In other words it looks something like a $[1, 2]$-tensor field.
Let's do this out formally with the geometry. Suppose you have two different connections $\nabla$ and $\nabla'$. Form the difference operator between them, $$\Delta_a = \nabla'_a - \nabla_a.$$ Recall that these both map scalar fields to the same value -- there was no ambiguity about that scalar gradient field! So $\Delta_a s = 0.$ But that means something very powerful, because $\Delta_a$ is Leibniz: it means that $\Delta_a (s~v^b) = s~\Delta_a v^b$. So it's a linear mapping of vector fields to tensor fields. In particular, this means that $u^a \Delta_a$ maps a vector field to another vector field linearly. Add in a covector $w_b$ and you get $u^a w_b \Delta_a v^b$ being a linear mapping from two vector fields $u, v$ and one covector field $w$ to a scalar: and that was precisely our definition of a $[1, 2]$ tensor field. So, in fact, there exists a tensor $D$ such that $u^a w_b \Delta_a v^b = D_{ac}^b u^a v^c w_b.$ Since this holds for all $u, w$ we can remove those and equivalently say, $$\nabla'_a v^b = \nabla_a v^b + D_{ac}^b v^c.$$ And this argument that $D$ needs to exist can also be run backwards, "assume we add this tensor term to $\nabla$, then we get another connection." So this is both necessary and sufficient.
Then of course, we exploit this freedom to get a case where $\nabla_a \nabla_b = \nabla_b \nabla_a$ and $\nabla_a g_{bc} = 0$ where $g$ is the metric tensor, and that is the Levi-Civita connection. But this is already a very long answer. I will give you a hint: define $\Delta_{ab} = \nabla_a\nabla_b-\nabla_b\nabla_a$ and use a slightly more interesting version of the above argument to argue that this is actually a derivation on scalars and therefore takes the form $T^c_{ab} \nabla_c$, this $T$ is the torsion tensor. What does changing our connection by adding $D$ do to it? 
A: A non technical or satisfying (even for me), but simple approach is to firstly consider physical laws in flat spacetimes. There, the directional derivative of a tensor (with components $T^{\alpha\beta}$ along a vector field (with components $x^\mu$) has components $x^\mu\partial_\mu T^{\alpha \beta}$. In flat spacetimes, this term transforms as a tensor under the coordinate transformations that leave the (Minkowski) metric tensor invariant, which are the Lorentz transformations. They are linear, i.e. ${\Lambda^\mu}_\nu\equiv\partial x^\mu/\partial x^\nu$ are constants. Then \begin{equation}\partial_\mu({\Lambda^\alpha}_\nu\,{\Lambda^\beta}_\rho T^{\nu\rho})={\Lambda^\alpha}_\nu\,{\Lambda^\beta}_\rho\partial_\mu T^{\nu\rho}.\tag{1}\end{equation}
If you are now interested in writing similar terms but that preserve the tensor character under general transformations that leave a general metric tensor invariant you need a new object (let's call it $\tilde\partial_\mu$) such that:
(i) recovers the ordinary derivative in a locally inertial frame, since we want general relativity to be true in non gravity situations, and
(ii) satisfies in any reference frame the same properties that satisfies in locally inertial frames. This would imply linearity and Leibniz rule, the usual things for derivatives.
Clearly $\tilde\partial_\mu$ can't be equal to $\partial_\mu$ in any reference frame, since in general coordinate transformations the components ${\Lambda^\mu}_\nu$ are coordinate dependent, and (1) is no longer true. You can write $\tilde\partial_\mu=\partial_\mu+D_\mu$, where $D_\mu$ is metric tensor-dependent, and $D_\mu=0$ in flat spacetimes. Now, I guess that imposing compatibility of this derivative with a metric tensor would give that $D_\mu$ are related to the usual Christoffer symbols for the Levi-Civita connection, and you'd be able to do the reverse process and get your coordinate free and metric tensor independent definition for the directional covariant derivative that you gave in your question.
In summary, I'd say that you just want something that behaves as a derivative but that its action on a tensor is also a tensor, since that doesn't happen with the ordinary derivative. 
A: Note: When I originally wrote this post, I misread you, and basically missed the whole thing about the motivation being "physical". I however spent a very long time writing this answer up and I am not gonna delete it. Hopefully, this will be useful to you but if not to you then to somebody else finding this question. With that said, I added a section at the end that gives "physical" motivation for the covariant derivative. This section is indicated by the boldface starting sentence.
The motivation is that when you move to a manifold instead of a vector space, you lose the ability to differentiate tensor fields.
If $T$ is some tensor field with components $T^{\mu_1...\mu_r}_{\nu_1...\nu_s}$, then the derivative $\partial_\sigma T^{\mu_1...\mu_r}_{\nu_1...\nu_s}$ does not transform as a tensor. The reasons why this is so is usually discussed in the literature.
If we want to skip any axiomatic definition of a differential operator, we still have some options. One is to realize that the reason "usual differentiation" fails is because a vector located at $x$ is an element of $T_x M$ and a vector located at $y$ is an element of $T_yM$, these are separate vector spaces, comparision is impossible.
We then introduce the notion of parallel transport. If $\gamma:\mathbb{R}\rightarrow M$ is a smooth curve, then let $P_\gamma(t_1,t_0):T_{\gamma(t_0)}M\rightarrow T_{\gamma(t_1)}M$ be a parallel transport map called a parallel propagator associated with the curve. It moves a vector located at $T_{\gamma(t_0)}M$ to $T_{\gamma(t_1)}M$.
Some axiomatics are needed here:


*

*We want parallel transport to be a linear transformation.

*We want parallel transport to be invertible.

*We want $P_\gamma(t_0,t_0)=\text{Id}$ and $P_\gamma(t_1,t_0)^{-1}=P_\gamma(t_0,t_1)$.

*We want $P_\gamma(t_1,t')P_\gamma(t',t_0)=P_\gamma(t_1,t_0)$.

*We want $P_\gamma$ to depend smoothly on both $t_1$ and $t_0$, and we want "$P$" to depend smoothly on $\gamma$, the latter being quite hard to actually describe mathematically.


Once we have this, we can define the following: If $V$ is a vector field along $\gamma$ (rigorously speaking, it is a "section" of the form $V:\mathbb{R}\rightarrow TM$ such that $\pi\circ V=\gamma$), then we define the covariant derivative of $V$ along $\gamma$ at $t_0$ as $$ \left.\frac{d^\nabla V}{dt}\right|_{t_0}=\lim_{t\rightarrow t_0}\frac{P_\gamma(t_0+t,t_0)^{-1}V(t_0+t)-V(t_0)}{t-t_0}. $$
To evaluate this map explicitly, we need to do some modifications.
We let $(U,\psi)$ be a local chart in the neighborhood of $\gamma(t_0)=x$, and we denote coordinates as $x^\mu$. Since $P_\gamma(t_1,t_0)$ is a linear transform between finite dimensional spaces, it is representable as a matrix, provided bases are chosen in both vector spaces. The local chart gives us a chosen basis, so we have for $v=v^\mu\partial_\mu|_{\gamma(t_0)}$, $P_\gamma(t_1,t_0)v=P_\gamma(t_1,t_0)^\mu_{\ \nu}v^\nu\ \partial_\mu|_{\gamma(t_1)}$. To ensure that $P$ maps invariant vectors to invariant vectors, we need the upper index on $P$'s matrix representation to transform as a vector at $\gamma(t_1)$ and the lower index to transform as a vector at $\gamma(t_0)$, so $P_\gamma(t_1,t_0)$ is essentially a two-point tensor.
The actual modifications happen now. Instead of considering a single curve $\gamma$, consider a vector field $X$ and its flow $\phi^X$, where $\phi^X(x_0,t)$ is the instruction to move along the integral curve that starts at $x_0$ for the time period $t$.
Let $P_X(x_0,t)$ denote $P_\gamma(t,0)$, where $\gamma$ is the integral curve that starts at $x_0$. What we actually have here is the following dependencies: $P$ is actually a composite function in the way $P_X=P\circ\phi^X$, so we have $P_X(x_0,t)=P(\phi^X(x_0,t))$. If $P_X(x_0,t)^\mu_{\ \nu}$ is a matrix representation, we have $$ \left.\frac{d}{dt}\right|_{t=0}P_X(x_0,t)^\mu_{\ \nu}=\left.\frac{\partial P^\mu_{\ \nu}}{\partial (\phi^X)^\sigma}\right|_{\phi^X=\phi^X(x_0,0)}\left.\frac{d(\phi^X)^\sigma}{dt}\right|_{x=x_0,t=0}.  $$
This is confusing because pretty much all notation for derivatives is terrible in some ways, but the flow $\phi^X$ is always the identity for $t=0$, so we actually have $\phi^X(x_0,0)=x_0$, so the first derivative could be written as $\partial P/\partial x_0^\sigma$, which is absolutely terrible, because $P$ is not something that actually depends directly on positions but for the sake of readability I'll write it that way. We have, then $$ \frac{\partial P^\mu_{\ \nu}}{\partial x^\sigma}X^\sigma(x_0), $$ since the time derivative of the flow is the vector field itself.
All this is needed to finally be able to have, for a $V$ that is no extended to be defined on a suitable open region, instead of just along a curve, $$ \left.\frac{d^\nabla V}{dt}\right|_{t=0,x=x_0}=\nabla_X V|_{x=x_0}=\lim_{t\rightarrow 0}\frac{P_X(x_0,t)^{-1}V(\phi^X(x_0,t))-V(x_0)}{t}=\left.\frac{d}{dt}\right|_{t=0}[P_X(x_0,t)^{-1}V(\phi^X(x_0,t))]. $$
We want to express this in terms of local coordinates. Before we do that we note that if $A(t)$ is a $t$-dependent matrix that is invertible for all $t$s, and $A(0)=I$, then we have $$ \frac{d}{dt}(A^{-1})|_{t=0}=-\frac{d}{dt}A|_{t=0}, $$ you can verify that yourself by differentiating the $I=A(t)A^{-1}(t)$ expression at zero.
Also, we a priori name $\frac{\partial P^\mu_{\ \nu}}{\partial x^\sigma}$ as $-\Gamma^\mu_{\sigma\nu}$.
Local coordinate expressions follow as $$ \nabla_XV|_{x=x_0}=\left.\frac{d}{dt}[\left.P_X(x_0,t)^{-1}\right.^\mu_{\ \nu}V^\nu(\phi^X(x_0,t))]\right|_{t=0}\ \partial_\mu|_{x_0}= \\ =\left(\frac{\partial \left.P^{-1}\right.^\mu_{\ \nu}}{\partial x^\sigma}X^\sigma(x_0)V^\nu(x_0)+\delta^\mu_\nu\frac{\partial V^\nu}{\partial x^\sigma}X^\sigma\right)\partial_\mu|_{x_0}= \\ = \left(\Gamma^\mu_{\sigma\nu}X^\sigma V^\nu+\partial_\sigma V^\mu X^\sigma\right)\partial_\mu, $$ where in the last line all expressions are to be evaluated at $x_0$ and in the middle line the Kronecker delta appeared because $P_X^{-1}$ at $t=0$ is just the identity.
From this expression, we can read off all properties of the covariant derivative, for example that it is tensorial in $X$ and that it still makes sense if $V$ is only defined along a curve.

Remarks: As you can see, this approach if far more laborous than defining an algebraic differential operator. And my statement that $P_X=P\circ\phi^X$ is actually somewhat iffy. It is believable but I honestly do not know how to make this derivation without this "iffy" statement or even do this the coordinate-free way. The parallel propagator's actual functional dependencies are extremely non-trivial.
But this approach has the advantage that we start off from an easy-to motivate concept of parallel translating vectors along curves, and the familiar covariant derivative fell out nicely at the end.
If you are curious about motivating the Levi-Civita covariant derivative, we can add to the list of requirements of the parallel transport that parallel transport preserves the lengths and angles of vectors. When you define covariant derivatives of tensors of arbitrary rank, then this requirement naturally implies that the metric tensor is parallel-transported along all curves. Torsionlessness cannot be motivated this easily though.
This motivation wasn't based on any sort of physics, though, instead I tried to make the covariant derivative intuitive by starting from the fact that we can parallel transport vectors in euclidean space, but you cannot in manifolds in general. So we, knowing what properties does good ol' parallel transport have, we put it in by hand.
If you want some really physical motivation, the best we can have is to follow Weinberg, and base GR on the equivalence principle instead of Riemannian geometry. The two are actually equivalent because equivalence principle $\Longleftrightarrow$ Riemannian normal coordinates $\Longleftrightarrow$ Riemannian geometry, and the implications are all two-way.
According to the equivalence principle, at around any $x$ spacetime event it is possible to set up coordinates, for which at $x$ and in its first-order infinitesimal neighborhood, the laws of special relativity apply.
Let $\xi^0,...,\xi^3$ be these special coordinates, and let $x^0,...,x^3$ be completely general coordinates. In addition, let primed indices refer to the special coordinate system and unprimed indices refer to the general coordinate system.
If $V^\mu$ is some vector field, then the expression $\partial_\mu V^\nu$ is valid in special relativity, and it only contains first derivatives, so let us interpret this expression to be made in the special coordinate system at the point $x$, and let us write it as $\partial_{\mu'}V^{\nu'}$. By the equivalence principle, this expression is valid.
Let us introduce the notation $\partial_{\mu'}V^{\nu'}=\nabla_{\mu'}V^{\nu'}$ for the primed indices, and let $\nabla_\mu V^\nu$ mean the tensor - transformed form of this expression in the general coordinate system, so $$ \nabla_\mu V^\nu=\frac{\partial \xi^{\mu'}}{\partial x^\mu}\frac{\partial x^\nu}{\partial \xi^{\nu'}}\nabla_{\mu'}V^{\nu'}.$$
We would like to relate the expression $\nabla_\mu V^\nu$ to the partial derivatives of $V^\nu$ in the general coordinate system.
Note that $$ \partial_\mu V^\nu=\frac{\partial \xi^{\mu'}}{\partial x^\mu}\partial_{\mu'}\left(\frac{\partial x^\nu}{\partial \xi^{\nu'}}V^{\nu'}\right)=\frac{\partial \xi^{\mu'}}{\partial x^\mu}\frac{\partial^2 x^\nu}{\partial \xi^{\mu'}\partial \xi^{\nu'}}V^{\nu'}+\frac{\partial \xi^{\mu'}}{\partial x^\mu}\frac{\partial x^\nu}{\partial \xi^{\nu'}}\partial_{\mu'}V^{\nu'}, $$ and here the second term on the RHS is essentially $\nabla_\mu V^\nu$, so we substract the first term on the RHS from the expression with the substitution $V^{\nu'}=\frac{\partial\xi^{\nu'}}{\partial x^\sigma}V^\sigma$.
What we get is $$ \nabla_\mu V^\nu=\partial_\mu V^\nu -\frac{\partial^2x^\nu}{\partial\xi^{\mu'}\partial\xi^{\nu'}}\frac{\partial\xi^{\mu'}}{\partial x^\mu}\frac{\partial\xi^{\nu'}}{\partial x^\sigma}V^\sigma=\partial_\mu V^\nu+\Gamma^\nu_{\mu\sigma}V^\sigma, $$ where we named $$ \Gamma^\nu_{\mu\sigma}=-\frac{\partial^2x^\nu}{\partial\xi^{\mu'}\partial\xi^{\nu'}}\frac{\partial\xi^{\mu'}}{\partial x^\mu}\frac{\partial\xi^{\nu'}}{\partial x^\sigma}. $$
Notes:


*

*All expressions are evaluated at the chosen point $x$, since these special coordinates are only "special relativistic" at that one point.

*This reasoning is more "physical", because the equivalence principle is essentially the main physical postulate behind GR.

*This approach has the advantage that the covariant derivative is immediately torsionless and metric-compatible, however it has the disadvantage that there is not closed-form expression for the Christoffel symbols that only reference the general coordinate system. This can be remedied by using the metric compatibility condition to derive the usual expression for $\Gamma$.
