Let $\rho$ be a density on some space $M$. For all practical purposes some subspace of $\mathbb R^n$. Let $v$ be a smooth vector field with flow $\Phi_t$. The transport equation, as far as I understand, governs the time evolution of $\rho_t$, which is defined as the pushforward density of $\rho$ along the flow $\Phi_t$: \begin{equation} \rho_t:=\rho\circ\Phi_{t}^{-1}\; \qquad (=\rho\circ\Phi_{-t}) \end{equation} So in particular $\rho_0=\rho$.
I have seen two forms of the transport equation: $0=\partial_t \rho_t+\nabla \cdot(\rho_t v)$ and $something=\partial_t\rho_t+v\cdot \nabla \rho_t$. I only seem to be able to derive the latter with $something=0$ and the former only by the heuristic argument that the mass change of density in a volume is equal to the total drainage out of that volume (using Gauss' law). I'd be helpful for guidance or hints. Good literature is rare, it seems.
So here is my derivation of the second which relies also on a result that I have not seen the proof of. Let $f$ be a smooth function of compact support on $M$. \begin{align} \int_Mdx \;f\;\partial_t\big|_s\rho_t =&\partial_t\big|_s\int_M dx \; f(x)\;\rho(\Phi_{-t}(x))\\ =&\partial_t\big|_s\int_{\Phi_t(M)=M} dy \;\left|\frac{\partial \Phi_t}{\partial y}\right| f(\Phi_t(y)) \rho(y)\\ =& \int_{M} dy \;\rho(y)\partial_t\big|_s\left|\frac{\partial \Phi_t}{\partial y}\right| f(\Phi_s(y)) + \int_{M} dy \; \rho(y)\left|\frac{\partial \Phi_s}{\partial y}\right| \partial_t\big|_sf(\Phi_t(y))\\ =& \int_{M} dy \; \rho(y)\partial_t\big|_s\left|\frac{\partial \Phi_t}{\partial y}\right| f(\Phi_s(y)) + \int_{M} dy \; \rho(y)\left|\frac{\partial \Phi_s}{\partial y}\right| \big(v\cdot\nabla f\big)(\Phi_s(y)) \\ =& \int_{M} dy \; \rho(y)\left(\nabla\cdot v\right)(\Phi_s(y))\left|\frac{\partial \Phi_s}{\partial y}\right| f(\Phi_s(y)) \\ &+\int_{M} dy \; \rho(y)\left|\frac{\partial \Phi_s}{\partial y}\right| \big(v\cdot\nabla f\big)(\Phi_s(y))\\ =& \int_{\Phi_{-s}(M)=M} dx \; \rho(\Phi_{-s}(x))\big(\nabla\cdot v\big)(x)f(x) \\ &+\int_{\Phi_{-s}(M)=M} dx \; \rho(\Phi_{-s}(x)) \big(v\cdot\nabla f\big)(x)\\ =& \int_M dx\; \nabla\big(f\;v\rho_s\big)-\int_M dx \; f \;v\cdot \nabla \rho_s \end{align} where we used that $\Phi_t$ is a diffeomorphism (i.e. $\Phi_t(M)=M$). Now we just need to apply Stokes' theorem (i.e. Gauss' law in this case) and that $f$ has compact support (unbounded $M$) or $v\big|_{\partial M}\cdot d\vec S\leq 0$ or something like that to kill the first term.
Question 1: Is this derivation correct? (I guess not)
Question 2: How do i get the weak equation $0=\partial_t\big|_s \rho_t+\nabla \cdot(v\rho_s)$
Question 3: Where do I find a proof of $\partial_t\big|_s\left|\frac{\partial \Phi_t}{\partial y}\right|=\left(\nabla\cdot v\right)(\Phi_s(y))\left|\frac{\partial \Phi_s}{\partial y}\right|$ (which I found for dim=3 here ) and also stated in another similar post (Trouble with the derivation of the Reynolds Transport Theorem).
Question 4: Any hints at good literature?
Let $\mu$ be a measure with density $\rho$, so $d\mu(x)=\rho(x)dx$, and let $v$ and $\Phi_t$ as above. Then the density of $\mu_t:=(\Phi_t)_*\mu$ can be inferred from the the definition of the push forward measure: \begin{equation} \int f \,d\,(\Phi_t)_*\mu =\int f\circ \Phi_t\; d\mu = \int dx \, f(\Phi_t(x)) \, \rho(x)=\int dy\, \left|\frac{\partial \Phi_{-t}}{\partial y}\right|\!(y)\,\rho(\Phi_{-t}(y)) f(y) \end{equation} where $f$ is a test function as above. The density of $\mu_t$ is therefore $\rho_t(x)=\left|\frac{\partial \Phi_{-t}}{\partial x}\right|\!(x)\,\rho\big(\Phi_{-t}(x)\big)$. We can check whether this density satisfies the continuity equation. We don't even need Jacobi's formula (stated in the question).
Proposition: $\rho_t$ satisfies the weak continuity equation: $0=\int dx \Big[\partial_t|_s\rho_t+\nabla(v\rho_s)\Big]\,f$.
Proof: \begin{align} \int dx \;v\cdot \left(\nabla\rho_s\right) f=& -\int dx \rho_s\nabla\cdot(vf)+\int dx\;\nabla\cdot\left(v\rho_s\,f\right)\\ =&-\int dx \rho_s\nabla\cdot(vf) \end{align} by Gauss' law. Then we change coordinates $y=\Phi_{-s}(x)$: \begin{align} \int dx \;v\cdot \left(\nabla\rho_s\right) f=&-\int dy \;\rho(y) \big(\nabla \cdot vf\big)\circ \Phi_s(y)\\ =&-\int dy \;\rho(y) \big(\nabla \cdot v\big)(\Phi_s(y))f(\Phi_s(y)))-\int dy \;\rho(y) \big(v\cdot\nabla f\big)\circ \Phi_s(y)\\ =& -\int dy \;\rho(y) \big(\nabla \cdot v\big)(\Phi_s(y))f(\Phi_s(y))) -\int dy \;\rho(y) \frac{\partial}{\partial t}\bigg|_{s}\,f(\Phi_t(y))\\ =& -\int dx \; \rho_s(x)\big(\nabla \cdot v\big)(x)\; f(x) -\frac{\partial}{\partial t}\bigg|_{s}\,\int dy \;\rho(y) f(\Phi_t(y))\\ =& -\int dx \; \rho_s\big(\nabla \cdot v\big) \; f -\frac{\partial}{\partial t}\bigg|_{s}\,\int dx \;\rho_t \;f\;, \end{align} where we changed back to the original coordinate system $x=\Phi_{s}(y)$ in the last and penultimate line. Pulling the time derivative inside the integral (imposing the relevant conditions needed to do that) and rearranging yields the result: \begin{align} 0&=-\int dx \;v\cdot \left(\nabla\rho_s\right) f-\int dx \; \rho_s\big(\nabla \cdot v\big) \; f -\int dx \;\frac{\partial}{\partial t}\bigg|_{s}\,\rho_t \;f\;,\\ &=-\int dx \Big[\partial_t|_s\rho_t+\nabla(v\rho_s)\Big]\,f \end{align}