I'm trying to solve 2D unsteady compressible Navier-Stokes equations with finite-difference or finite-volume method. Here is the system, it's pretty standard:

$$
\frac{\partial U}{\partial t} + 
\frac{\partial F}{\partial x} +
\frac{\partial G}{\partial y} = 0;
$$

$$
U = \left( \begin{array} {c} \rho \\ \rho u \\ \rho v \\ e \end{array} \right), \;\;
F = \left( \begin{array} {c} \rho u \\ \rho u^2 + p -\tau\_{xx} \\ \rho uv -\tau\_{xy} \\ (e + p)u -u\tau\_{xx} -v\tau\_{xy} - k \frac{\partial T}{\partial x}\end{array} \right), \;\;
G = \left( \begin{array} {c} \rho v \\ \rho uv -\tau\_{xy} \\ \rho v^2 + p -\tau\_{yy} \\ (e + p)v -u\tau\_{xy} -v\tau\_{yy} - k \frac{\partial T}{\partial y}\end{array} \right)
$$
where $\tau$ is viscous stress tensor:
$$
\tau\_{xx} = \frac{4}{3}\mu \frac{\partial u}{\partial x} - \frac{2}{3}\mu \frac{\partial v}{\partial y} \\
\tau\_{xy} = \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \\
\tau\_{yy} = \frac{4}{3}\mu \frac{\partial v}{\partial y} - \frac{2}{3}\mu \frac{\partial u}{\partial x}
$$

Suppose that we discretize non-viscous flux terms (i.e. those parts of F and G vectors that do not contain derivatives) independently with some finite difference or finite-volume scheme, maybe with WENO reconstruction. Then viscous fluxes may be discretized and added as source terms.

With some of derivatives, it's pretty straightforward, for example:

$$
\frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial x} \right) = \mu \frac{\partial^2 u}{\partial x^2} \approx \mu \frac{u\_{i+1} - 2u\_{i} + u\_{i+1}}{\Delta x^2}
$$
(second order approximation for uniform $\mu$).

But there are more complex derivatives in the viscous fluxes, such as:
$$
\frac{\partial}{\partial x} \left( \mu \frac{\partial v}{\partial y} \right), \;
\frac{\partial}{\partial x} \left( \mu u \frac{\partial u}{\partial x} \right), \;
\frac{\partial}{\partial x} \left( \mu u \frac{\partial v}{\partial y} \right), \;
\frac{\partial}{\partial x} \left( \mu v \frac{\partial u}{\partial y} \right), \;
\frac{\partial}{\partial x} \left( \mu v \frac{\partial v}{\partial x} \right)
$$
and so on.
And here are the questions:

1. What is the most simple way to discretize all these derivatives under finite-difference approach? I assume it would yield 2nd-order approximation.
2. How to use high-order central differences for them under FD approach?
3. What is the simplest way under finite-volume approach? I'm aware that we should use Gauss theorem in some form; there is some info in [Blazek's book](https://www.sciencedirect.com/book/9780080999951/computational-fluid-dynamics-principles-and-applications), but it's not detailed enough to be easily understandable for me. Are there more detailed books/papers on this?
4. How to apply high-order schemes in FV approach? It seems that method given in Blazek's book gives only 2nd-order approximation.
5. How to discretize these derivatives for non-uniform viscosity $\mu$ and heat conductivity $k$?