Does the covariant derivative commute with the variational derivative? Say we have a (compactly supported normal) variation of a surface $\Sigma$ immersed in $\mathbb{R}^3$,  $\mathbf{r}: \Sigma \times \mathbb{R} \to \mathbb{R}^3$ with $\frac{\partial{\mathbf{r}}}{\partial t} = u \mathbf{N}$ where $\mathbf{N}$ is the unit normal vector field of the surface and $u: \Sigma \times \mathbb{R} \to \mathbb{R}$ is some smooth function.
Let $f_{ij}$ denote the twice-covariant derivative of a smooth function $f: \Sigma \times \mathbb{R} \to \mathbb{R}$, with respect to $i$ followed by $j$.  In computing stuff like $\frac{\partial}{\partial t}\,(f_{ij})$, does it make sense to write $$\frac{\partial}{\partial t}\,(f_{ij}) =^? (f_t)_{ij}$$ where $f_t$ is the $t$-partial derivative of $f$?
My hunch is that it's more complicated than this, but I don't have  a good reason for my skepticism.  Any help is much appreciated.  Thanks in advance.
 A: I am not completely sure about the context of your question and whether you mean by $f_ij$ the $ij$ coordinate of $\nabla \operatorname{d}f$ or whether you firstly apply $\nabla_i$ on $f$ followed by application of $\nabla_j$ (here $\nabla$ refers to the covariant derivative). This answer is focused on the second case which is not tensorial (unlike the first case!). Hopefuly my answer won't be missleading and will provide some help. 
$\Sigma$ is some embedded surface in $\mathbb{R}^3$ and thus has dimension $2$ so we need (and have, at least locally) $2$ coordinate functions that describes $\Sigma$. Since the domain of $f$ is $\Sigma \times \mathbb{R}$, we are working with $3$ coordinates and hence it makes sense to consider derivatives in different coordinate directions $i,j$ and $t$. Important point: the covariant derivative of a function is simply the directinal derivative (classical partial derivative). That is, according to the interpretation of your notation discussed above, $f_{ij} = \frac{\partial^2 f}{\partial x^i \partial x^j} $. From the Schwartz theorem (see this wiki article) we know that the partial derivatives of a function commute. Putting these two facts together we have 
$$\frac{\partial}{\partial t} f_{ij} = f_{ijt} = f_{itj} = f_{tij} = (f_t)_{ij} \ .$$
A: Since the metric on $\Sigma$ is changing in time, it is indeed more complicated than this. I think it's easiest here to do a coordinate calculation; so we'll start with the formula $$f_{;ij} = \partial_i \partial_j f - \Gamma_{ij}^k \partial_k f$$for the second covariant derivative. Differentiating this in time and commuting partial derivatives we find 
\begin{align}
\partial_t f_{ij} &= \partial_i \partial_j \partial_tf-\Gamma^k_{ij}\partial_k\partial_tf-(\partial_t\Gamma^k_{ij})\partial_kf\\
&=(\partial_tf)_{ij} - (\partial_t \Gamma^k_{ij})\partial_kf; \tag{1}
\end{align}
so your proposed formula needs a correction term due to the variation of the induced connection. To get a formula for the variation of $\Gamma$, we will first  need to know the variation of the induced metric, which we can find with the product rule: we have
\begin{align}
\def\r{\mathbf{r}}\def\v{\mathbf{v}}\def\N{\mathbf{N}}
\partial_t g_{ij} &= \partial_t\left(\partial_i \r \cdot \partial_j \r\right)
\\ &=\partial_i \v\cdot\partial_j\r+\partial_i\r\cdot\partial_j\v
\end{align}
where $\v = \partial_t \r = u \N$ is the velocity. Expanding this using the facts $\N\cdot \partial_j r= 0 $ and $\partial_i \N \cdot \partial_j r = A_{ij}$ (here $A$ is the second fundamental form) we find
$$h_{ij} := \partial_t g_{ij} = 2u A_{ij}.$$
Plugging this in to the formula $$\partial_t \Gamma^k_{ij} = \frac 12 g^{kl}\left(-\nabla_l h_{ij} + \nabla_i h_{jl} + \nabla_j h_{li} \right) \tag{2}$$ for the variation of the connection yields
$$\partial_t \Gamma^k_{ij} = g^{kl}\left( -A_{ij}u_l +A_{jl}u_i + A_{li} u_j-u\nabla_l A_{ij} + u\nabla_i A_{jl} + u\nabla_j A_{li}\right).$$ (You can prove (2) easily using geodesic normal coordinates - if you get stuck, it's Prop 3.6 here.)
Sticking this back into (1) gives you the formula you want. You can simplify it a little using the fact that $\nabla A$ is totally symmetric (from the Codazzi equation), which lets you cancel the fourth and fifth term.
