Mesh orientation effect for the temperature propagation

Describe the bug
The temperature propagation is mesh-oriented for a hexahedral inclined mesh (tstep 1):

There are 3 meshes used: hexahedral rectangular, hexahedral inclined, and tetrahedral.
Neither of them work for the purpose:

• hexahedral inclined - has the mesh-orientation effect
• hexahedral rectangular - is not aligned to materials boundary, and it was made just to check the mesh orientation effect
• tetrahedral - doesn't converge

Results for the tetrahedral mesh (only timesteps before the simulation divergence are available):

Results for the hexahedral rectangular mesh:

To Reproduce
Steps to reproduce the behavior:

1. Run the models (the link is provided below)
2. Load the resulting vtk files into Paraview and observe the issue

Expected behavior

1. The temperature propagation front is perpendicular to Z-axis

Platform (please complete the following information):
Machine: RHEL 8.10
Compiler: gcc 13.1.0
GEOS Version: 1.1.0
GEOS Version sha: develop https://github.com/GEOS-DEV/GEOS/commit/1e617be8614817d92f0a7a159994cbed1661ff98 , tpls: 9687680c61fa3378f7b2ce7eea8327736809fb63

Models and logs
You can find them in this repository: https://gitlab.com/i-l-s-h-a-t/geos_models/-/tree/main/THM/triaxial/thermal_issue/mesh_effect?ref_type=heads

Additional context
The idea of the model is to simulate a triaxial experiment for the heterogenous rock.
Model uses 2 materials (without contact mechanics), currently they have the same properties. It uses singlePhase singleComponent physics with TPFA discretization.

Next boundary conditions are defined: no load from the top and side boundaries, lower boundary has a roller condition except the central part where it is fixed. Cooling is imposed at the top side by setting the temperature lower than the initial.

I thought that HybridFVM might help and found there are number of features not implemented for the case:

1. There is no iterative linear solver for the SinglePhaseHybridFVM:
   PoroelasticitySolver: MGR strategy is not implemented for thermal SinglePhasePoromechanics/SinglePhaseHybridFVM
2. After switching to the direct linear solver: The thermal option is not supported by SinglePhaseHybridFVM

Checking the code for CompositionalMultiphaseHybridFVM showed there is no such an exception so maybe adding a dummy component and phase is a solution, but I have not checked yet. Anyway, it would be good to resolve the issue for the SinglePhase physics.

This issue could be related to this https://github.com/GEOS-DEV/GEOS/issues/3739.

Reproduce the models from the linked geos_models repository under THM/triaxial/thermal_issue/mesh_effect, then inspect the SinglePhase physics path using TPFA discretization. Compare the resulting VTK files in ParaView for the inclined, rectangular, and tetrahedral meshes. Done means the temperature propagation front is perpendicular to the Z-axis without the reported mesh-orientation effect.