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