I have the following coupled equations for the conservation of mass and momentum of a compressible fluid
\begin{equation} \rho_t + (\rho u)_z = 0, \end{equation}
$$ (\rho u)_t + (\rho u^2)_z + \tfrac{1}{k}\rho_z = \left(\frac{dp}{dz}\right)_f + \rho g $$ over the domain $[0, 2L]\times[0, \infty)$ for density and velocity $\rho = \rho(z, t)$, $u = u(z, t)$ respectively. Note that $$ g = g(z) = \begin{cases} 9.81 & \text{if } 0 \le z \le L \\ -9.81 & \text{if } L < z \le 2L \end{cases} $$ and that $\left(\frac{dp}{dz}\right)_f$ is a frictional pressure term that is non-linear in $\rho$ and $u$.
The weak variational problem is: $\textit{find $u\in V_1, \rho\in V_2$ such that}$ \begin{equation} \int_0^{2L} v\rho_t - v_z\rho u dz = 0 \end{equation} $$ \int_0^{2L}q(\rho u)_t - q_z\rho u^2 - \tfrac{1}{k}q_z\rho - q\left(\frac{dp}{dz}\right)_f - q\rho gdz = 0 $$ $\textit{for all test functions $v, q$ in spaces $V_1$, $V_2$ respectively.}$
I am then discretising these equations in time using an implicit backward Euler discretisation for implementation in the finite element solver software Firedrake (this is a standard approach when using Firedrake).
The problem I'm having is that the Firedrake solver does not converge with the inclusion of $\left(\frac{dp}{dz}\right)_f$, and gives the wrong solution when $\left(\frac{dp}{dz}\right)_f$ is excluded.
So in order to narrow down the causes of this divergence, my question is: are there specific spaces for $V_1$ and $V_2$ that should be used for this kind of variational problem?