Reynolds transport theorem: link with the Lie derivative?

Instead of $\rho$, think about the three form $\rho ~\mathrm{d}x$. (Like I tell my calculus students, you should always write the integration element in an integral.) Here $\mathrm{d}x$ is the standard volume form.

By definition of the divergence of a vector field you have that $$ \mathrm{div}(\rho \vec{v}) ~\mathrm{d}x = \mathrm{d}( \iota_{\rho\vec{v}} ~\mathrm{d}x ) = \mathrm{d} (\iota_{\vec{v}} \rho ~\mathrm{d}x ) = L_{\vec{v}} \left(\rho~\mathrm{d}x\right) $$ which is what you expect.

There is indeed a connection. Recall the definition of the Lie derivative for a $k$-form $\omega$ along a vector field $X$. If $\chi$ is the one-parameter family of diffeomorphisms defined (at least locally) by $X$, then $$L_X\omega = \left.\frac{d}{d\tau}\right|_0(\chi_{\tau}^{*}\omega)$$ Integrating over some $k$-submanifold $\Omega$, we get $$ \int_\Omega L_X\omega = \int_{\Omega}\left.\frac{d}{d\tau}\right|_0(\chi_{\tau}^{*}\omega) = \left.\frac{d}{d\tau}\right|_0\int_{\Omega}\chi_{\tau}^{*}\omega = \left.\frac{d}{d\tau}\right|_0\int_{\chi_\tau[\Omega]}\omega $$ Note that we first used the Leibniz integral rule in one variable, and then the fact that diffeomorphisms preserve the integral. I claim Reynold's transport theorem is a special case of the formula $$\left.\frac{d}{d\tau}\right|_0\int_{\chi_\tau[\Omega]}\omega = \int_\Omega L_X\omega$$ we just obtained.

Say your space is an $n$-dimensional manifold $M$, and consider an interval $I\subseteq \mathbb R$ (say with $0\in I$) which for us will represent time. In the context of Reynold's transport theorem you have a submanifold $\Omega$ of $M$ (say of dimension $k$) which "varies (smoothly) over time". We can interpret this as having a fixed manifold $\Omega$ and a family of embeddings $i_t:\Omega\to M$. Using these embeddings, we can trace the "worldsheet" of $\Omega$ inside the "spacetime" $I\times M$. Define another family of embeddings $j_t:\Omega\to I\times M$ given by $$j_t(x)=(t,i_t(x)),$$ and define $\Omega_t:=j_t[\Omega]\subseteq I\times M$. Now we want to get the "spacetime velocity" of $\Omega$, which will be a vector field defined on $\Theta:=\bigcup_{t\in I}\Omega_t$. To do this, we observe that the embeddings $j_t$ are smooth, so they define a family of smooth curves $\gamma_x:I\to I\times M$ given for every $x\in\Omega$ by $\gamma_x(t)=(t,i_t(x))$. The vectors tangent to this family of curves give us a vector field $Y$ defined by $$Y_{(t,i_t(x))}=\gamma'_x(t)$$ for every $(t,i_t(x))\in\Theta=\bigcup_{t\in I}\Omega_t$. This is the spacetime velocity we wanted. Using the extension lemma for vector fields on submanifolds, we get a vector field $X$ on $I\times M$ which agrees with $M$ on $\Theta$ (we actually need to show that $\Theta$ is a submanifold, but I think this follows from the smoothness of all the maps).

Similarly, if you had a $k$-form $\bar{\omega}$ on $M$ which depends on time $t\in I$, then it can be regarded as a $k$-form $\omega$ on $M\times I$ by defining $$\omega_{(t,m)}(v_{(t,m)}) := \bar{\omega(t)_m}(\text{proj}_{TM}(v_{(t,m)}))$$ for every tangent vector $v_{(t,m)}\in T_{(t,m)}(I\times M)$.

Now applying the formula we first derived to the submanifold $\Omega_0$, we have $$\left.\frac{d}{dt}\right|_0\int_{\chi_t[\Omega_0]}\omega = \int_{\Omega_0} L_X\omega$$ where $\chi$ is the one-parameter family of diffeomorphisms induced by $X$ (we might need to show that $\chi$ is actually defined for all $t\in I$, but I think this is true at least in a neighbourhood of $\Omega_0$). Since $X$ extends $Y$, I think one can show that $\chi_t[\Omega_0]=\Omega_t$. Also, by Cartan's magic formula, we have $L_X\omega = d(\iota_X\omega) + \iota_X(d\omega)$. Hence

$$\left.\frac{d}{dt}\right|_0\int_{\Omega_t}\omega = \int_{\Omega_0} d(\iota_X\omega) + \int_{\Omega_0}\iota_X(d\omega)$$

When integrating the right hand side we can take charts consisting of a chart on $I$ times a chart on $M$, so $\omega$ has no $dt$ component because we defined it using $\bar\omega$. This way, letting $V=\text{proj}_{TM}Y$ we have $\iota_Y\omega=\iota_V\omega$ and even $\iota_X\omega=\iota_Y\omega=\iota_V\omega$ on $\Omega_0$. On the other hand, since $\Omega_0$ is contained in $\{0\}\times M\subseteq I\times M$ we have $dt|_{\Omega_0}=0$. These observations give the first two equalities in $$\int_{\Omega_0} d(\iota_X\omega)=\int_{\Omega_0} d(\iota_V\omega)=\int_{\Omega_0} d_M(\iota_V\omega)=\int_{\partial\Omega_0}\iota_V\omega$$ and the third one comes from Stokes' theorem.

Similarly, over $\Omega_0$ we have $X=Y=\frac{\partial}{\partial t}+V$, and since $\omega=\omega_Idx^I$ (where $I$ is a multiindex where $t$ does not appear) we have $$\begin{aligned} d\omega &=\partial_t\omega_Idt\wedge dx^I+\sum_i \partial_i\omega_Idx^i\wedge dx^I \\ &=\partial_t\omega_Idt\wedge dx^I+d_M\omega \end{aligned}$$ and the second integral on the right hand side becomes $$\begin{aligned} \int_{\Omega_0}\iota_X(d\omega) &= \int_{\Omega_0}\iota_{\frac{\partial}{\partial t}+V}(\partial_t\omega_Idt\wedge dx^I+d_M\omega) \\ &= \int_{\Omega_0}\left[\iota_{\frac{\partial}{\partial t}}(\partial_t\omega_Idt\wedge dx^I)+\iota_{\frac{\partial}{\partial t}}(d_M\omega)+\iota_V(\partial_t\omega_Idt\wedge dx^I)+\iota_V(d_M\omega)\right] \\ &= \int_{\Omega_0}\left[\partial_t\omega_Idx^I+\iota_V(d_M\omega)\right] \\ &= \int_{\Omega_0}\dot\omega+\int_{\Omega_0}\iota_V(d_M\omega) \end{aligned}$$ (note the second and third terms vanish because they pair "spacelike" basis forms with "timelike" basis vectors and viceversa).

Putting al this together, we have $$\left.\frac{d}{dt}\right|_0\int_{\Omega_t}\omega = \int_{\partial\Omega_0}\iota_V\omega+\int_{\Omega_0}\dot\omega+\int_{\Omega_0}\iota_V(d_M\omega)$$ as desired.