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$.