The key fact is the condition that $\gamma$ is a geodesic is strictly stronger than the condition that the image of $\gamma$ is totally geodesic.
This is clear from the equations - saying $\gamma$ is a geodesic means $\nabla_{\dot{\gamma}} \dot{\gamma} = 0$ while the image of $\gamma$ being totally geodesic means $$\nabla_{\dot{\gamma}} \dot{\gamma} - \frac{g(\nabla_{\dot{\gamma}}\dot{\gamma}, \dot{\gamma} )}{g(\dot{\gamma}, \dot{\gamma})}\dot{\gamma} = 0.$$
If the image of $\gamma$ is totally geodesic it does follow that by reparamaterizing $\gamma$, it does become a geodesic.
Now, suppose $\Gamma = (\gamma_1(t), \gamma_2(t))$. If $\gamma_1$ and $\gamma_2$ are both geodesics (meaning $\nabla_{\dot{\gamma}_i} \dot{\gamma}_i = 0$ for $i = 1,2$), then by your above equations, $\Gamma$ is a geodesic.
On the other hand, if $\gamma_1$ is a geodesic and $\gamma_2$ merely has totally geodesic image, then $\Gamma = (\gamma_1(t), \gamma_2(t))$ need not be totally geodesic. Your equations show exactly why you shouldn't expect it to be, but here's a concrete counterexample.
Consider $M_1 = M_2 = \mathbb{R}$. Let $\gamma_1(t) = t$ and let $\gamma_2 = t^3$. Then $\gamma_1$ is a geodesic and the image of $\gamma_2$ is totally geodesic because there is no normal direction at all. (If this is too trivial, take $M_1 = M_2 = \mathbb{R}^2$, $\gamma_1(t) = (t,0)$ and $\gamma_2(t) = (t^3, 0)$ - the rest of the argument will work in either case). Then in $M_1\times M_2$, the image of $\Gamma(t) = (\gamma_1(t), \gamma_2(t))$ is the graph of a cubic polynomial, so is NOT a straight line, so is not a geodesic at all. Nor is the subset totally geodesic.
I recently answer my question by finding a formula I didn't know.
Let $\nabla$ be a torsion-free connexion and $X$ a vector field. Then $\mathcal{L}_X\nabla$ is a tensor and
\begin{align}
\mathcal{L}_X\nabla &= -i_X\circ R^{\nabla} + \nabla^2X
\end{align}
where $R^{\nabla}(U,V) = \nabla_{[U,V]} - [\nabla_U,\nabla_V]$ is the curvature tensor of $\nabla$, and $\nabla_{U,V}^2X = \nabla_U\nabla_VX - \nabla_{\nabla_UV}X$. Applying this to $\nabla_{\partial_j}S$ we get
\begin{align}
\mathcal{L}_{\partial_r}\left(\nabla_{\partial_j}S\right) &= \mathcal{L}_{\partial_r}(\nabla)(\partial_j,S) + \nabla_{[\partial_r,\partial_j]}S + \nabla_{\partial_j}(\mathcal{L}_{\partial_r}S)
\end{align}
and using the above formula and the Riccati equation for $S$ leads to a linear differential equation.
Best Answer
The following is an important fact which basically means that Levi-Civita connections behave well under isometric embeddings. (It can actually be stated for isometric immersions as well).
Fact (or exercise, if you like): Let $\widetilde{M}$ be a Riemannian manifold, and let $\widetilde{\nabla}$ denote the Levi-Civita connection of $\widetilde{M}$. Let $M$ be a submanifold of $\widetilde{M}$, equipped with the induced Riemannian metric, and let $X,Y$ be vector fields on $M$. Let $\widetilde{X},\widetilde{Y}$ be vector fields on $\widetilde{M}$ which extend $X,Y$, respectively. Then:
i) The orthogonal projection of $\widetilde{\nabla}_{\widetilde{X}}\widetilde{Y}$ on $TM$ is independent of the extensions $\widetilde{X},\widetilde{Y}$.
ii) Let $\nabla_XY$ denote the orthogonal projection of i). Then $\nabla$ is the Levi-Civita connection of $M$.
Now, to your question. The Levi-Civita connection of $\mathbb{R}^3$ (with the standard metric) is just the usual derivative. As follows from the above fact, what described in the post is covariant derivation with respect to the Levi-Civita connection on the surface $S$.
Historical remark: Gauss and his gang knew how to take covariant derivatives on surfaces many years before Levi-Civita was born, and before the term "connection" was ever used in its current geometrical meaning. However, the approach described above is the one which is most commonly used today. It helps one to put all the pieces together.