Definition
Suppose $X,Y$ are locally convex topological vector spaces and $U\subseteq X$ is open. We say a functional \begin{align} J:\ U &\longrightarrow Y \\ u &\longmapsto J[u] \end{align} is Gâteaux differentiable at $u$ if forall $h\in X$, there exist the limit \begin{align} \lim_{\varepsilon\to0} \dfrac{J[u+\varepsilon h]-J[u]}{\varepsilon} =: \delta J[u;h]. \end{align}
In the Riemannian geometry setting
Let $M$ be a Riemannian manifold, $Met(M)$ denotes the space of all smooth metric tensor fields on $M$ and $S^0_2(M)$ denotes the space of all smooth symmetric $(0,2)-$tensor field on $M$. Then the Levi-Civita Theorem shows that for each $g\in Met(M),\ h\in S^0_2(M)$ and for $\varepsilon\in\mathbb R$ small enough so that $g+\varepsilon h\in Met(M)$, there exists unique Levi-Civita connections $\nabla[g+\varepsilon h]$ and $\nabla[g]$.
Proposition (4.13 page 94, John M. Lee - Introduction to Riemannian Manifolds) asserts that for all $\varepsilon$, the difference $\nabla[g+\varepsilon h]-\nabla[g] $ forms a $(1,2)-$tensor field. Hence, for $\varepsilon\ne0$, $\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon}$ is also a $(1,2)-$tensor field.
Then, I interest that if the Gâteaux derrivative $\delta\nabla[g;h] = \displaystyle\lim_{\varepsilon\to0}\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon} $ is also a $(1,2)-$tensor field ? In this context, the spaces $U\subseteq X$ and $Y$ in the precede definition are $Met(M)\subseteq S^0_2(M)$ and $ C^\infty(M) $.
So I found in Lee's Introduction to smooth manifolds the Lemma 12.24 page 318 and Introduction to Riemmanian manifolds the Lemma B.6 page 398 states that the map $\delta\nabla[g;h]:\ \mathfrak X^*(M)\times\mathfrak X(M)\times\mathfrak X(M)\longrightarrow C^\infty(M)$ is a $(1,2)-$tensor field if and only if it is multilinear over $C^\infty(M)$.
I deduce that since $\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon}$ is multilinear over $C^\infty(M)$, and \begin{align} \left(\lim_{\varepsilon\to0}\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon}\right)(\omega,X,Y) &= \lim_{\varepsilon\to0}\left(\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon}(\omega,X,Y)\right) \end{align} for all $\omega\in\mathfrak X^*(M), X, Y\in\mathfrak X(M) $, it follows that $\displaystyle \lim_{\varepsilon\to0}\dfrac{\nabla[g+\varepsilon h]-\nabla[g]}{\varepsilon}$ is multilinear over $C^\infty(M)$, which is justified by the above Lemmas to be a $(1,2)-$tensor field on $M$.
Is my above argument correct ? And is this enough for one to be able to apply the covariant derivative for the Gâteaux derivative $\delta\nabla[g;h]$ ? I would appreciate any help. Thanks.