I would like to numerically solve a hyperbolic PDE of the form

$\frac{\partial\theta_t}{\partial t}(x,y)+\frac{\partial\left[\theta_t \gamma_t^x\right]}{\partial x}(x,y)+\frac{\partial\left[\theta_t \gamma_t^y\right]}{\partial y}(x,y)=0,$

which is very similar to the 2D advection equation, except that the partial derivatives w.r.t $x$ and $y$ are products of my dependent variable $\theta_t(x,y)$ with another function $\gamma_t^x(x,y)$ or $\gamma_t^y(x,y)$.

Since I can calculate the $\gamma_t^x(x,y)$ and $\gamma_t^y(x,y)$ functions and their derivatives exactly, I don't need to approximate these and can just use the product rule to expand the terms with these in. My PDE is therefore of the form

$\frac{\partial\theta_t}{\partial t}(x,y)+\gamma_t^x(x,y)\frac{\partial\theta_t}{\partial x}(x,y)+\gamma_t^y(x,y)\frac{\partial \theta_t}{\partial y}(x,y) +\theta_t(x,y)\frac{\partial\gamma_t^x}{\partial x}(x,y)+\theta_t(x,y)\frac{\partial \gamma_t^x}{\partial y}(x,y) =0,$

I need to use a numeric scheme to approximate the derivatives with $\theta_t(x,y)$ in them. Numeric schemes which work for solving a 2D advection equation (where the gammas are absent and the partial derivatives are simply multiplied by a constant) do not seem to work for me. I have tried several (upwind, downwind, adaptive, etc.), but can't find one which yields conservation of mass for $\theta_{t}$.

Any advice would be greatly appreciated.