Time derivative of flux We have a time and even "position" invariant vector field and a surface. If the surface is moving with constant velocity, is the flux through the surface should constant in time?
Also, is there an easy to follow proof for the formula
$$\frac{d}{dt} \iint_{S(t)}\overline{V}(\overline{r},t) \cdot d\overline{A}  $$
 A: Here is a proof, though I'm not sure it is easy to follow.
In below derivation, we will use the adjective "nice" for anything sufficiently regular
to make the argument works. If I'm not mistaken, continuous differentiable up to second order is sufficient.
Let


*

*$\vec{v}(\vec{x},t)$ be a "nice" velocity field.

*$\vec{B}(\vec{x},t)$ be a "nice" vector field which we want to compute the flux.

*$S$ be a compact "nice" surface with piecewise "nice" boundary in $\mathbb{R}^3$.  

*$S(t)$ be a family of surfaces generated by $S$ using the velocity field $\vec{v}$.
More precisely, for each $\vec{x} \in S$, consider the initial value problem:
$$\frac{d}{dt} \gamma_{\vec{x}}(t) = \vec{v}(t,\gamma_{\vec{x}}(t)),\quad t \in (-\epsilon, \epsilon )
\quad\text{ with initial condition }\quad \gamma_{\vec{x}}(0) = \vec{x}.
$$
Define a function by 
$$\Gamma : S \times (-\epsilon, \epsilon ) \quad\mapsto\quad \Gamma(\vec{x},t) = \gamma_{\vec{x}}(t)\in\mathbb{R}^3$$
Standard theory of ODE tell us $\Gamma$ is a "nice" function (as least for small enough $\epsilon$).
$S(t)$ is simply the image $\Gamma(S \times \{t\})$.
As a geometric object, I'm not 100% sure $S(t)$ is "nice".
For the sake to get a result, let's assume they are.
The number we want to compute is following time derivative 
$$\mathscr{F}_{t} = \frac{d}{dt}\int_{S(t)} \vec{B}\cdot d\vec{A}$$
where $d\vec{A}$ is the area element associated with $S(t)$.
To carry out the computation, let us assume $S$ is small enough to fit into a single coordinate chart 
$$\varphi : \mathbb{R}^2 \supset \mathcal{O} \ni (r,s) \quad\mapsto\quad \varphi(r,s) \in S \subset \mathbb{R^3}$$
We will further assume the domain $\mathcal{O}$ has piecewise "nice" boundary${}^{\color{blue}{[1]}}$.
Using $\Gamma$, we can extend it to a function
$$\vec{X} :  \mathcal{O} \times (-\epsilon, \epsilon ) \ni (r, s, t) \quad\mapsto\quad
\vec{X}(r,s,t) = \Gamma(\varphi(r,s),t)$$
For any function (or expression) $\psi$ that depends on $(r,s,t)$, we will adopt
the shorthand 
$$\psi_r = \frac{\partial\psi}{\partial r},\quad
\psi_s = \frac{\partial\psi}{\partial s},\quad\text{ and }\quad
\psi_t = \frac{\partial\psi}{\partial t}$$
In terms of $\vec{X}$, the area element $d\vec{A}$ is simply $\vec{X}_r \times \vec{X}_s dr ds$ and
$$\begin{align}
\mathcal{F}_t 
&= \frac{d}{dt} \int_{\mathcal{O}} \vec{B}\cdot (\vec{X}_r \times \vec{X}_s) dr ds\\
&= \int_{\mathcal{O}}\left[
\frac{d\vec{B}}{dt} \cdot (\vec{X}_r \times \vec{X}_s)
+ \vec{B} \cdot \left( (\vec{X}_r)_t \times \vec{X}_s +
\vec{X}_r \times (\vec{X}_s)_t \right)
\right] dr ds\\
&= \int_{\mathcal{O}}\left[
\left(\vec{B}_t + \color{firebrick}{(\vec{v}\cdot\vec{\nabla})\vec{B}}\right) \cdot 
\color{firebrick}{(\vec{X}_r \times \vec{X}_s)}
+ \vec{B} \cdot\left( \vec{v}_r \times \vec{X}_s + \vec{X}_r \times \vec{v}_s\right)
\right] dr ds
\end{align}
$$
Notice the second term in the integrand
$$
 \vec{B} \cdot\left( \vec{v}_r \times \vec{X}_s + \vec{X}_r \times \vec{v}_s\right)
= \vec{B} \cdot
\left( ( \vec{v} \times \vec{X}_s)_r - ( \vec{v} \times \vec{X}_r)_s \right)$$
can be rewritten as
$$
( \vec{B}\cdot( \vec{v} \times \vec{X}_s) )_r -
( \vec{B}\cdot( \vec{v} \times \vec{X}_r) )_s
+ \color{red}{((\vec{X}_r\cdot\vec{\nabla})B) \cdot (\vec{X}_s \times \vec{v})}
+ \color{red}{((\vec{X}_s\cdot\vec{\nabla})B) \cdot (\vec{v}\times\vec{X}_r)}
$$
Now for any five vectors $\vec{a},\vec{b},\vec{c}, \vec{p}, \vec{q} \in \mathbb{R}^3$, we have an identity${}^{\color{blue}{[2]}}$:
$$
  (\vec{p}\cdot\vec{q})(\vec{a}\cdot(\vec{b}\times\vec{c}))
= (\vec{p}\cdot\vec{a})(\vec{q}\cdot(\vec{b}\times\vec{c}))
+ (\vec{p}\cdot\vec{b})(\vec{q}\cdot(\vec{c}\times\vec{a}))
+ (\vec{p}\cdot\vec{c})(\vec{q}\cdot(\vec{a}\times\vec{b}))
$$ 
$\vec{\nabla} \otimes \vec{B}$ is a rank 2 tensor, we can decompose it as a sum of outer product of vectors
$$\vec{\nabla} \otimes \vec{B} = \vec{p}_1 \otimes \vec{q}_1 + \vec{p}_2 \otimes \vec{q}_2 + \cdots$$
If we substitute $\vec{a},\vec{b},\vec{c}$ in above identity by $\vec{v}$, $\vec{X}_r$ and $\vec{X}_s$ respectively, it is not hard to deduce:
$$\color{firebrick}{((\vec{v}\cdot\vec{\nabla})\vec{B}) \cdot (\vec{X}_r \times \vec{X}_s)}
+ \color{red}{((\vec{X}_r\cdot\vec{\nabla})B) \cdot (\vec{X}_s \times \vec{v})
+ ((\vec{X}_s\cdot\vec{\nabla})B) \cdot (\vec{v}\times\vec{X}_r)}
= \color{green}{( \vec{\nabla}\cdot \vec{B}) (\vec{v}\cdot(\vec{X}_r \times \vec{X}_s)}.$$ 
Substitute this back into the integrand for $\mathcal{F}_t$, we get
$$
\mathcal{F}_t 
= \int_{\mathcal{O}}
\left[
\left((\vec{B}_t + \color{green}{(\vec{\nabla}\cdot \vec{B}) \vec{v}} )\cdot 
\color{green}{( \vec{X}_r \times \vec{X}_s )}\right)_{\color{blue}{\verb/I/}}
+ 
\left(
( \vec{B}\cdot( \vec{v} \times \vec{X}_s) )_r -
( \vec{B}\cdot( \vec{v} \times \vec{X}_r) )_s
\right)_{\color{blue}{\verb/II/}}
\right] dr ds\\
$$
The integrand split into two pieces. It is obvious what Piece $\color{blue}{\verb/I/}$ is. 
For piece $\color{blue}{\verb/II/}$, we 
can transform it first to a line integral along $\partial\mathcal{O}$ using the classical Green's theorem in $\mathbb{R}^2$. We then re-express it as a line integral along $\partial S(t)$ in $\mathbb{R}^3$. Finally, we convert it back to a surface integral over $S(t)$ using Kevin-Stokes theorem:
$$\begin{align}
\text{Piece }\color{blue}{\verb/I/} &= \int_{S(t)}\left( \vec{B}_t + (\vec{\nabla}\cdot\vec{B})\vec{v}\right)\cdot d\vec{A}\\
\\
\text{Piece }\color{blue}{\verb/II/} &= 
\int_{\partial S(t)} ( \vec{B}\times\vec{v} ) \cdot d\vec{X}
=\int_{S(t)} \vec{\nabla} \times (\vec{B} \times \vec{v}) \cdot d\vec{A}
\end{align}
$$
Combine this, we get
$$\mathcal{F}_t = \int_{S(t)} \left( \vec{B}_t + (\vec{\nabla}\cdot\vec{B})\vec{v} - 
\vec{\nabla} \times (\vec{v} \times \vec{B})\right)\cdot d\vec{A}\tag{*1}$$
When $S$ is too large to fit in a single coordinate chart, we can split $S$ into finite many
pieces that fit. We then apply $(*1)$ to the individual piece and sum their contributions. 
At the end, $(*1)$ continue to work...
Notes


*

*$\color{blue}{[1]}$ $\partial\mathcal{O}$ need to be nice enough to apply classical Green's theorem.

*$\color{blue}{[2]}$ When $\vec{a},\vec{b},\vec{c}$ are linear independent, 
they form a basis of $\mathbb{R}^3$. The correspond dual basis consists of
$\frac{\vec{b} \times \vec{c}}{\Delta}$, $\frac{\vec{c} \times \vec{a}}{\Delta}$, and $\frac{\vec{a} \times \vec{b}}{\Delta}$ where $\Delta = \vec{a}\cdot(\vec{b}\times\vec{c})$.
The identity in main text is nothing special but a formula of the dot product for this particular pair of basis/dual basis.

A: ${\color{#66f}{\large\tt\mbox{This is from $\tt\color{#c00000}{G\&R}$ Table,}\ 7^{\rm\underline{a}}\ \mbox{ed.}}}$

