Harmonic average of Diffusion Tensors in Finite Volume Method

I want to implement a Bilinear Finite Volume discretisation of the anisotropic diffusion problem:

$$\frac{du}{dt} = \nabla \cdot (\textbf{D} \nabla u)$$

Both my degrees of freedom as well as the Diffusion Tensors are given on the vertices of my quadrilateral grid. I can, therefore, use bilinear test/trialfunctions on each cell to calculate the full gradient.

I am wondering and am deeply puzzled as to how I am to approximate the Diffusion tensor on an internal sub-face of my control volume.

In the 1D case (TPFA), the diffusion coefficient is modeled to be constant over one cell and the degree of freedom is thought to be the same everywhere within the cell. The flux over an interface is made consistent by using the harmonic average of the two given diffusivities:

$$f_{A->B}= \frac{2}{\frac{1}{D_A}+ \frac{1}{D_B}} ~\nabla_x u$$

Usually, the gradient is then approximated by the simple finite difference of the values modeled to be in the cell centers. sketch to clarify: So far so good. Now I am not in a setting where the Grid is aligned with the principal directions of the diffusion tensors in involved (not "K-orthogonal"). The Two-point flux approximation would, therefore, fail, because it can not capture the off-diagonal entries of the diffusion tensors involved, because only the x-derivative can be calculated.

In my 2D setting, I do reconstruct the full gradient. I want to calculate four sub-fluxes within my cell as indicated by the following sketch: Here the black dots indicate where my degrees of freedom lie. The dashed line shows the area where I model the diffusion tensor to be constant (dual grid), the friendly green dot marks the location where I want to evaluate the bottom-most of the four sub-fluxes to properly assemble. Following the procedure from the 1D case, I could use the harmonic average of the diffusion tensors of the involved dual cells A,B. Let's call it

option 1:

$$D_{eff} = \big(\frac{\textbf{A}^{-1 } + \textbf{B}^{-1} }{2}\big)^{-1}$$

Here the "-1" in the exponents indicates the matrix inverse. What puzzles me so much is that by calculating the matrix inverse, the diagonal entries become coupled with the off-diagonal ones. To my understanding, these are different physical processes and I feel awkward to have, say, The component $D_{xx}$ depend on $D_{xy}$ etc.

The second option I see is to do a harmonic average component-wise, to have:

option 2

$$D_{eff} = \begin{bmatrix} \frac{2}{\frac{1}{D^{A}_{xx}} + \frac{1}{D^{B}_{xx}}} & \frac{2}{\frac{1}{D^{A}_{xy}} + \frac{1}{D^{B}_{xy}}} \\ \frac{2}{\frac{1}{D^{A}_{yx}} + \frac{1}{D^{B}_{yx}}} & \frac{2}{\frac{1}{D^{A}_{yy}} + \frac{1}{D^{B}_{yy}}} \end{bmatrix}$$

If i do it this way, the diagonal component $D_{xx}$ is only dependend on the corresponding entries in $D_A$ and $D_B$.

What is the correct way to do average the Diffusion Tensors at the sub-face marked with the green dot?

UPDATE: Investigating option two I found that as the offdiagonal entries of the Diffusion Matrices may be negative, I face the situation where I take the harmonic average of a negative and a positive number which, according to wikipdeia, is not defined properly. Help appreciated.