Full disclosure, I posed this question on MSE about a week ago, but realized it may be a better fit here after it sat dormant. For completeness, I'll include the full body of the question here:
Let $\boldsymbol{u}:\mathbb R^n\to \mathbb R^n$ be a $C^2$ vector field, representing the velocity of a (steady) fluid flow. If we let $\Phi_t(\boldsymbol x)$ be the flow map for the field $\boldsymbol u$, i.e. $$\partial_t\Phi_t(\boldsymbol x) = \boldsymbol u(\Phi_t(\boldsymbol x)),\quad \Phi_0(\boldsymbol x) = \boldsymbol x$$ for every $\boldsymbol x$, then $(\boldsymbol x,t)\mapsto \Phi_t(\boldsymbol x)$ is $C^2$ and $$\frac{d}{dt} D\Phi_t(\boldsymbol x) = Du(\Phi_t(\boldsymbol x))D\Phi_t(\boldsymbol x).$$ So, by Jacobi's formula, we find that $\det D\Phi_t$ is $C^1$ with $$\frac{d}{dt}\det D\Phi_t = (\det D\Phi_t(\boldsymbol x))(\nabla\cdot \boldsymbol u)(\Phi_t(\boldsymbol x)).$$
Let $\Omega\subseteq \mathbb R^n$ open bounded connected with $C^1$ boundary, so that $\Omega$ is a nice packet of fluid. Set $\Omega_t = \Phi_t[\Omega]$ for $t \in \mathbb R$. Then, for $t > 0$, a change of variables (the details of which I'll spare the reader) yields the Reynolds transport theorem $$\frac{d}{dt}\int_{\Omega_t}\rho(\boldsymbol x, t)\,d\boldsymbol x = \int_{\Omega_t} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\boldsymbol u)\,d\boldsymbol x$$ for any $\rho: \mathbb R^n\times [0,\infty)\to \mathbb R$ which is (say) jointly $C^1$ in $(\boldsymbol x, t)$ and $C^1$ in $t$ uniformly in $\boldsymbol x$. In particular, plugging in $\rho = 1$, we find $$\frac{d}{dt} \mathcal L^n(\Omega_t) = \int_{\Omega_t}\nabla\cdot\boldsymbol u\,d\boldsymbol x.$$ This little argument proves the perhaps intuitive claim that "the flow of $\boldsymbol u$ preserves volume globally iff $\boldsymbol u$ preserves volume locally, i.e. $\nabla\cdot\boldsymbol u = 0$, $\boldsymbol u$ is incompressible."
A similar argument, made for appropriate vector fields $\boldsymbol F$, yields $$\frac{d}{dt}\int_{\partial \Omega_t} \boldsymbol F(\boldsymbol x, t)\cdot \boldsymbol \nu(\boldsymbol x, t)\,d\mathcal H^{n-1}(\boldsymbol x) = \int_{\partial \Omega_t}\left(\frac{\partial \boldsymbol F}{\partial t}+(\nabla\cdot\boldsymbol F)\boldsymbol u\right)\cdot \boldsymbol \nu\,d\mathcal H^{n-1}(\boldsymbol x)$$ where $\boldsymbol \nu(\boldsymbol x, t)$ is a continuous outward normal to $\partial \Omega_t$ at a point $\boldsymbol x$.
Here's my main question. Try as I might, I cannot figure out how to extend this style of reasoning to understand the quantity $$\frac{d}{dt}\int_{\partial \Omega_t}\rho(\boldsymbol x, t)\,d\mathcal H^{n-1}(\boldsymbol x).$$ In particular, I'm interested in which flows $\boldsymbol u$ preserve surface area. This seems much more restrictive than volume preservation – an informal argument using infinitesimal spheres suggests that $\boldsymbol u$ should probably be incompressible: I find it hard to believe that the volume of (the forward image of) a small ball around a point with nonzero divergence can change while the surface area stays constant without causing some spikes in the boundary (which can't exist because $\partial\Omega_t$ is the image of a $C^1$ manifold under a $C^1$ diffeomorphism.)
I think some argument using the isoperimetric inequality $\mathcal H^{n-1}(\partial \Omega_t) \geq c\mathcal L^n(\Omega_t)^{1-1/n}$ gets us somewhere in this direction. This shows that if the flow of $\boldsymbol u$ preserves surface area, then the volume of $\Omega_t$ is bounded for all $t > 0$ by a (dimensional) constant multiple of the surface area of $\partial \Omega_0$. Maybe then you show that the surface can't collapse in on itself without creating spikes? Unfortunately isoperimetry only goes one way...
However, in 2D any flow of the form $(f(y), 0)$ is incompressible and can stretch surface areas wildly, so surface-area-preservation ought to be much stronger.
A concrete example of a 2D incompressible flow which doesn't preserve surface area is $\boldsymbol u(x,y) = (y,0).$ The flow map induced by $\boldsymbol u$ is $\Phi_t(x,y) = (x+ty,y),$ which stretches a square (perimeter $4$) to a parallelogram with perimeter $2+2\sqrt{1+t^2}$ at time $t$ if I did my computation correctly.
In conversation with some friends, we suspect that this condition implies $\boldsymbol u(\boldsymbol x) = A\boldsymbol x+\boldsymbol b$ for $\boldsymbol b \in \mathbb R^n$ and $A\in \mathbb R^{n\times n}$ antisymmetric, i.e. $A^T = -A$. Indeed, these flows are surface-area-preserving since $\exp(At)$ is orthogonal when $A$ is antisymmetric; they're also incompressible for the same reason. However, we don't have a rigorous proof for general $n$ that these are all of them. Any ideas? An expression for $\frac{d}{dt}\int_{\partial\Omega_t} \rho\,d\mathcal H^{n-1}$ in the same style as our other two derivations would be fantastic, though I'm perfectly happy with resolving the original question in the title, classifying all surface-area-preserving flows as precisely these guys.
I've proven that the flow of a $C^1$ vector field $\boldsymbol u$ is an isometry if and only if the field $\boldsymbol u$ has the form we suspect. So we suspect the following is true, but aren't sure how to prove it: if $\boldsymbol F:\mathbb R^n\to \mathbb R^n$ is a $C^1$ diffeomorphism which preserves the $(n-1)$-dimensional Hausdorff measure of embedded $C^1$ manifolds (with boundary), then $\boldsymbol F$ is an isometry. Any ideas towards proving or disproving this would also be incredibly helpful.