Consider the system of polynomial equations $$ \begin{cases} a(X^2+Y^2)^2-2b(X^2+Y^2)+4(bX-cY)X \equiv 0 &\pmod{p^k} \\ d(X^2+Y^2)^2+2c(X^2+Y^2)+4(bX-cY)Y \equiv 0 &\pmod{p^k} \end{cases} $$ in two variables $X,Y$, with $k\in\mathbb{Z}_{\geq1}$ and parameters $a,b,c,d\in \mathbb{Z}/p^k\mathbb{Z}$. How does one count or obtain an upper bound for the number of solutions $X,Y\in(\mathbb{Z}/p^{k}\mathbb{Z})^\times$?
So far I have thought of the following. Fixing a value for $X$, if the parameters $a,b,c,d$ are "nice enough", there can be at most 4 solutions for a value of $Y$ per equation due to Hensel's lemma. But it is not clear to me how one controls the number of admissible values for $X$. Is there a specific technique that one can use to get an idea of how the number of solutions behaves, perhaps in terms of the $p$-adic valuations of the parameters $a,b,c,d$?