Implicit methods for variable coefficients based on equations of state

There are many ways to do this. Your "predictor-corrector" approach is one way. A better conceptual framework is to do a time discretization first and then look at what kind of problem you have. For example, if you wanted to do an implicit Euler scheme, you'd have to solve the following problem in each time step: $$ frac{rho^n-rho^{n-1}}{Delta t} = -nablacdot (rho^n u) + nabla cdot(D(rho^n, T) nabla rho^n) + rho_s, $$ ignoring the other variables for the moment. That's a nonlinear problem to be solved in $rho^n$, and all techniques for solving nonlinear problems are on the table now. For example, you can do a Newton iteration that tries to find a sequence $rho^{n,k}$ ($k$ is the Newton iteration index in the $n$th time step) that you need to start with some $rho^{n,0}$. A good first guess would be $rho^{n,0}=rho^{n-1,ast}$, the Newton final iterate in the $(n-1)$st time step. An even better Newton iterate would be $rho^{n,0}=rho^{n-1,ast}+Delta t frac{rho^{n-1,ast}-rho^{n-2,ast}}{k} approx rho^{n-1,ast}+Delta t frac{partial rho(t)}{partial t}$. Another option would be to start with $rho^{n,0}=rho^{n,ee}$ where the latter term is the explicit Euler solution of the time step.

Many schemes can be seen as variations of any of these ideas. Maybe they do exactly one Newton iteration with any of these starting guesses, or they do one fixed point ("Picard") iteration with any of these starting guesses. But ultimately, if you want to solve the problem implicitly, it helps to first write it down as a nonlinear problem for $rho^n$.