Is metric compatibility an on-shell statement? In General Relativity (GR), it is standard to endow the spacetime manifold with a metric compatible connection, i.e. $\nabla_{\alpha}g_{\mu\nu}=0$. My question is: is metric compatibility of the connection an on-shell statement, i.e. does it only apply to metrics that satisfy the Einstein field equations (EFEs)? 
The reason I ask, is that when on derives the EFEs from the Einstein-Hilbert action, one considers an arbitrary variation of the metric around a metric that extremizes the action. Indeed, let $g_{\mu\nu}$ extremize the Einstein-Hilbert action (i.e. it is a putative solution to the EFEs), then consider a perturbation about this metric $\delta g_{\mu\nu} = \bar{g}_{\mu\nu} - g_{\mu\nu}$, where $\bar{g}_{\mu\nu}$ is a perturbed metric. Now, if the “on-shell” metric $g_{\mu\nu}$ satisfies $\nabla_{\alpha}g_{\mu\nu}=0$, it surely can’t be that $\nabla_{\alpha}\bar{g}_{\mu\nu}=0$, since this would imply that $\nabla_{\alpha}\delta g_{\mu\nu}=0$, and the equations of motion would look vastly different to the EFEs. So it appears that there are cases in which the variational operator $\delta$, and the covariant derivative $\nabla_{\mu}$ do not commute.
If my understanding is correct, then if we consider a variation of the metric $g_{\mu\nu}$, that is any solution to the EFEs, such that $\delta g_{\mu\nu} = \bar{g}_{\mu\nu} - g_{\mu\nu}$, however, this time $\bar{g}_{\mu\nu}$ is diffeomorphically related $g_{\mu\nu}$, i.e. $\bar{g}_{\mu\nu}=(\phi^{\ast}g)_{\mu\nu}$, then it should be that $\nabla_{\alpha}g_{\mu\nu}=0=\nabla_{\alpha}\bar{g}_{\mu\nu}$. 
In this sense, does the variational operator $\delta$ only commute with the covariant derivative $\nabla_{\mu}$ if the perturbed metric $\bar{g}_{\mu\nu} = g_{\mu\nu} + \delta g_{\mu\nu}$ is related to $g_{\mu\nu}$ by a diffeomorphism, such that they both satisfy the EFEs?
Edit
I think I can show that the variational operator $\delta$, and the covariant derivative $\nabla_{\mu}$ do not commute in general. Not sure if it’s a circular argument though. 
Consider a general tensor $A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n}$. Acting on it with the covariant derivative gives: $$\nabla_{\alpha}A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n} = \partial_{\alpha}A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n} + A^{\lambda,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n}\Gamma^{\mu_1}_{\;\alpha\lambda} +\cdots\\ \qquad\qquad\qquad + A^{\mu_1,\ldots,\alpha}_{\;\;\;\nu_1,\ldots,\nu_n}\Gamma^{\mu_n}_{\;\alpha\lambda} - A^{\mu_1,\ldots,\mu_n}_{\;\;\;\lambda,\ldots,\nu_n}\Gamma^{\lambda}_{\;\alpha\nu_1} - \cdots\\  - A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\lambda}\Gamma^{\lambda}_{\;\alpha\nu_n}$$ Then, acting from the left on both sides with the variational operator $\delta$ gives $$\delta\nabla_{\alpha}A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n} = \nabla_{\alpha}\delta A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n} + A^{\lambda,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\nu_n}\delta\Gamma^{\mu_1}_{\;\alpha\lambda} +\cdots\\ \qquad\qquad\qquad +  A^{\mu_1,\ldots,\alpha}_{\;\;\;\nu_1,\ldots,\nu_n}\delta\Gamma^{\mu_n}_{\;\alpha\lambda} -  A^{\mu_1,\ldots,\mu_n}_{\;\;\;\lambda,\ldots,\nu_n}\delta\Gamma^{\lambda}_{\;\alpha\nu_1} - \\  -  A^{\mu_1,\ldots,\mu_n}_{\;\;\;\nu_1,\ldots,\lambda}\delta\Gamma^{\lambda}_{\;\alpha\nu_n}$$ where we have used that the variational operator commutes with the partial derivative operator. We see them that $\delta$ and $\nabla_{\mu}$ do not commute in general.
 A: No, the (mathematical) property of non-metricity ($\nabla_\alpha g_{\mu \nu} \neq 0$) is a fundamental assumption you impose before you even write down the action (gravity + matter) whose physics you wish to describe. You can't follow the logic backwards, i.e., you can't use Einstein's equations to justify/prove non-metricity. You can't prove a fundamental mathematical assumption (which is completely oblivious to physics) using physical equations of motion of gravity.
Symmetric teleparallel gravity  is a formulation where the effect of gravity is described through non-metricity, while having vanishing curvature ('teleparallel') and vanishing torsion ('symmetric').
And yes, $\delta$ and $\nabla_\alpha$ do not commute (see your own comment in Sounak's answer). In fact, $[\delta,\nabla_\alpha] X^{\cdots}_{\ \ \ \ \cdots}$ is quite complicated, and depends on the tensorial structure of $X^{\cdots}_{\ \ \ \ \cdots}$. (Dots indicate Lorentz indices).
A: Normally, while deriving field equations from a Lagrangian, $L_m(\phi_A,\nabla_a\phi_A)$ with an action, $A=\int d^4x\sqrt{-g}L_m$ for matter fields, $\phi_A(x)$, a one-parameter family of field configurations, $\phi_A(\lambda,x), \lambda\in[0,1]$ is chosen, such that $\lambda=0$  corresponds to the classical field configurations and $\phi_A(\lambda,x)=\phi_A(0,x)$ on the boundary. Then upon differentiating the action w.r.t. $\lambda$ at zero, we get the field equations. Basically, in the following analysis I have replaced $\delta$ by $\frac{d}{d\lambda}$. Now, while following this procedure one gets a term like this, $$\frac{\partial L}{\partial\nabla_a\phi_A}\frac{d\nabla_a\phi_A}{d\lambda}$$, by the usual chain rule of functional derivatives. Since for standard matter field theories, the covariant derivative does not depend on the matter fields, one usually commutes the $\frac{d}{d\lambda}$ past the $\nabla$ and ending up with a surface term (which is zero since $\phi$ is a constant function of $\lambda$ on the boundary) and the usual Euler-Lagrange equations. However, this procedure does not work for gravity, mainly because now the covariant derivatives and the metric determinant all depend on $g$. Also, the $\frac{d}{d\lambda}$ and the $\nabla$ do not commute. For example, $$\nabla_c\frac{dg_{ab}}{d\lambda}=\lim_{\lambda\rightarrow0}\frac{1}{\lambda}(\nabla_c\bar g_{ab}(\lambda)-\nabla_cg_{ab})\neq0$$, while $$\frac{d}{d\lambda}\nabla_cg_{ab}=\lim_{\lambda\rightarrow0}\frac{1}{\lambda}(\bar\nabla_c\bar g_{ab}(\lambda)-\nabla_cg_{ab})=0$$. However, since $$A_G=\alpha\int d^4x (-g)^{1/2}g^{ab}R_{ab}+GHY$$, we can still write, $$\frac{dA_G}{d\lambda}=\alpha\int d^4x(\frac{dg^{ab}(-g)^{1/2}}{d\lambda}R_{ab}+(-g)^{1/2}g^{ab}\frac{dR_{ab}}{d\lambda})+...$$ by the chain rule. Since the first term does not involve any commuting of derivatives, we get the Einstein tensor times $\frac{dg^{ab}}{d\lambda}$. The second term however, requires some care. Let the connection associated with $g$ be $\nabla$, while the one associated with $\bar g$ be $\bar\nabla$. Since $\nabla_ag_{bc}=0$, we can always write,$$C^c_{ab}=\frac{1}{2}\bar g^{cd}(\lambda)(\nabla_a\bar g_{bd}+\nabla_b\bar g_{ad}-\nabla_d\bar g_{ab})$$, where $$\bar\nabla_ag_{bc}=\nabla_ag_{bc}+C^d_{ab}g_{dc}+C^d_{ac}g_{bd}$$ (I'm following Wald's notation so please consult his book if something is unclear, or leave a comment). Also, $$\frac{dC^a_{bc}}{d\lambda}=\lim_{\lambda\rightarrow0}\frac{C^a_{bc}(\lambda)}{\lambda}$$ since $C^a_{bc}(0)=0$. The relation $$\bar R_{abc}^{\quad d}=R_{abc}^{\quad d}-2\nabla_{[a}C^d_{b]c}+2C^e_{c[a}C^d_{b]e}$$, follows directly from the definition of the Riemann tensor in terms of the commutator of covariant derivates (just substitute for $\bar\nabla$ in terms of $\nabla$ and $C$). The last term, being higher order in $\lambda$ cancels out (keep in mind $C(0)=0$). Remember that you really can't contract the two sides of the above equation as you need to to multiply the L.H.S. with $\bar g(\lambda)$ while the R.H.S. needs to be multiplied with $g$. But you surely can sum over a and d, resulting in,
$$\bar R_{ab}=R_{ab}-2\nabla_{[a}C^b_{b]c}$$. Bringing the unperturbed Ricci tensor to the L.H.S., dividing by lambda and taking the limit gives,$$\dot R_{ab}=-2\nabla_{[a}\dot C^b_{b]c}$$. The expression for $\dot C$ which has already been given above, can be re-expressed as,$$\dot C^c_{ab}=\frac{1}{2}g^{cd}(\nabla_a\dot g_{bd}+\nabla_b\dot g_{ad}-\nabla_d\dot g_{ab})$$. Where dot denotes differentiation w.r.t. lambda and I have assumed $\nabla g=0$. I am allowed to take the one upon lambda inside the covariant derivative because it does not have any spatial variations, only variations over field configurations. The derivative of the Ricci tensor is, $$\dot R_{ac}=-\frac{1}{2}g^{bd}\nabla_a\nabla_c\dot g_{bd}--\frac{1}{2}g^{bd}\nabla_b\nabla_d\dot g_{ac}+g^{bd}\nabla_b\nabla_{(c}\dot g_{a)d}$$ (I copy pasted all of this from Wald's book) and after contracting it with the unperturbed metric (which can go inside the unperturbed covariant derivatives) we get a pure surface term (after playing around with the indices for a while obviously). It's variations are cancelled out by the variations of the GHY term at the boundary and also because $\dot g$ is zero there. Hence, since the metric derivatives w.r.t. $\lambda$ are arbitrary, we get the EFE  without assuming the commutativity of $\delta$ and $\nabla$. 
I hope this answers your question!
