The efficient and accurate simulation of coupled multi-physics phenomena is a fundamental challenge in various scientific and engineering disciplines. Such problems, which involve the interplay of multiple physical processes, often result in large and potentially ill-conditioned linear systems. The computational demand associated with solving these systems can be prohibitively high. As a consequence, researchers and practitioners have sought innovative solution strategies to address these computational challenges. In this work, we inspect the efficient solution of the four-field formulation of the non-linear thermo-poroelastic problem. The study includes preliminary results and theoretical findings on different strategies for addressing the solution of large - possibly ill-conditioned - linear system. For the spatial discretization, we design a high-order symmetric weighted interior penalty scheme that supports general polytopal grids and is robust with respect to strong heterogeneities in the model coefficients. When designing the solution strategy, a particular focus is devoted to the treatment of the non-linear convective transport term in the energy conservation equation proposing suitable stabilization techniques that make the scheme robust for advection-dominated regimes and trying to exploit the solution strategy for reducing the complexity of the problem. A broad set of numerical simulations is presented in order to validate the theoretical analysis, to inspect numerically the robustness properties, and to test the capability of the proposed method in a practical scenario inspired by a geothermal problem.