Temperature doesn't propagate properly with the Global stabilization enabled

Describe the bug
There is no convergence for a THM model with a tetrahedral mesh.
Enabling the Global stabilization helps the simulation to converge. However, the temperature propagates only to the first layer of cells:

I tried the Local stabilization just in case, but it is not supported:

geos/geos-develop/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp:155
***** Error cause: m_stabilizationType == stabilization::StabilizationType::Local
***** Rank 0: PoroelasticitySolver/stabilizationType (triaxial_base.xml, l.11): Local stabilization has been temporarily disabled

Changing the mesh to hexahedral makes the simulation converge without the stabilization. The temperature propagates. However with a mesh orientation effect, described in https://github.com/GEOS-DEV/GEOS/issues/3933:

To Reproduce
Steps to reproduce the behavior:

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

Expected behavior

1. The temperature propagation for both hexahedral/tetrahedral meshes when Global Stabilization is enabled
2. Converging simulation without the stabilization
3. Local Stabilization support with a proper temperature propagation, which Global Stabilization lacks

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

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 am wondering whether there are some ways to resolve the convergence issues with or without stabilization.

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

Start by running the triaxial models from the linked geos_models repository and inspecting the VTK output and logs in ParaView. Read geos/geos-develop/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp around line 155, then compare tetrahedral and hexahedral results. Done means the THM simulation converges and temperature propagates through both mesh types with the documented stabilization settings.