diff --git a/opm/models/discretization/common/tpsalinearizer.hpp b/opm/models/discretization/common/tpsalinearizer.hpp
index 91f7e039b09..eb534ab9276 100644
--- a/opm/models/discretization/common/tpsalinearizer.hpp
+++ b/opm/models/discretization/common/tpsalinearizer.hpp
@@ -258,6 +258,83 @@ class TpsaLinearizer
         }
     }
 
+    /*!
+    * \brief Compute the stress (or traction) for mechanics output
+    */
+    void updateStressInfo()
+    {
+        OPM_TIMEBLOCK(updateStressInfoTPSA);
+
+        if (stressInfo_.empty()) {
+            return;
+        }
+
+        const auto& geoMechModel = geoMechModel_();
+        auto& problem = problem_();
+        const unsigned int numCells = fullDomain_.cells.size();
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        // Interior traction of every cell, one per face to its neighbors
+        for (unsigned ii = 0; ii < numCells; ++ii) {
+            const unsigned globI = fullDomain_.cells[ii];
+            const MaterialState& materialStateIn =
+                geoMechModel.materialState(globI, /*timeIdx=*/0);
+
+            short loc = 0;
+            for (const auto& nbInfo : neighborInfo_[globI]) {
+                const unsigned globJ = nbInfo.neighbor;
+                const MaterialState& materialStateEx =
+                    geoMechModel.materialState(globJ, /*timeIdx=*/0);
+
+                // Compute local face term
+                ADVectorBlock adres(0.0);
+                LocalResidual::computeFaceTerm(adres,
+                                               materialStateIn,
+                                               materialStateEx,
+                                               problem,
+                                               globI,
+                                               globJ);
+                adres *= nbInfo.faceArea;
+
+                // Insert interior traction vector
+                // OBS: Assume traction vector is the three first entries in residual!
+                stressInfo_[globI][loc].faceNormal = problem.cellFaceNormal(globI, globJ);
+                stressInfo_[globI][loc].faceArea = nbInfo.faceArea;
+                for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
+                    stressInfo_[globI][loc].traction[tractionIdx] = adres[tractionIdx].value();
+                }
+                ++loc;
+            }
+        }
+
+        // Boundary traction
+        for (const auto& bdyInfo : boundaryInfo_) {
+            const unsigned globI = bdyInfo.cell;
+            const MaterialState& materialStateIn = geoMechModel.materialState(globI, /*timeIdx=*/0);
+
+            // Compute local boundary condition
+            ADVectorBlock adres(0.0);
+            LocalResidual::computeBoundaryTerm(adres,
+                                               materialStateIn,
+                                               bdyInfo.bcdata,
+                                               problem,
+                                               globI);
+            adres *= bdyInfo.bcdata.faceArea;
+
+            // Insert boundary traction vector
+            // OBS: Assume traction vector is the three first entries in residual!
+            const short loc = neighborInfo_[globI].size() + bdyInfo.bfIndex;
+            stressInfo_[globI][loc].faceNormal =
+                problem.cellFaceNormalBoundary(globI, bdyInfo.bfIndex);
+            stressInfo_[globI][loc].faceArea = bdyInfo.bcdata.faceArea;
+            for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
+                stressInfo_[globI][loc].traction[tractionIdx] = adres[tractionIdx].value();
+            }
+        }
+    }
+
     // ///
     // Public get and set functions
     // ///
@@ -496,7 +573,6 @@ class TpsaLinearizer
                 OPM_TIMEBLOCK_LOCAL(faceCalculationForEachCellTPSA, Subsystem::Assembly);
 
                 // Loop over neighboring cells
-                short loc = 0;
                 for (auto& nbInfo : nbInfos) {
                     OPM_TIMEBLOCK_LOCAL(calculationForEachFaceTPSA, Subsystem::Assembly);
 
@@ -536,16 +612,6 @@ class TpsaLinearizer
                     // SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
                     bMat *= -1.0;
                     *nbInfo.matBlockAddress += bMat;
-
-                    // Insert interior traction vector
-                    // OBS: Assume traction vector is the three first entries in residual!
-                    const auto& faceNormal = problem.cellFaceNormal(globI, globJ);
-                    stressInfo_[globI][loc].faceNormal = faceNormal;
-                    stressInfo_[globI][loc].faceArea = nbInfo.faceArea;
-                    for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
-                        stressInfo_[globI][loc].traction[tractionIdx] = res[tractionIdx];
-                    }
-                    ++loc;
                 }
             }
 
@@ -627,17 +693,6 @@ class TpsaLinearizer
             // Insert contribution to (globI, globI) sub-block
             // SparseAdapter syntax: jacobian_->addToBlock(globI, globI, bMat);
             *diagMatAddress_[globI] += bMat;
-
-            // Insert boundary traction vector
-            // OBS: Assume traction vector is the three first entries in residual!
-            const auto& nbInfos = neighborInfo_[globI];
-            short loc = nbInfos.size() + bdyInfo.bfIndex;
-            const auto& bndyNormal = problem.cellFaceNormalBoundary(globI, bdyInfo.bfIndex);
-            stressInfo_[globI][loc].faceNormal = bndyNormal;
-            stressInfo_[globI][loc].faceArea = bdyInfo.bcdata.faceArea;
-            for (unsigned tractionIdx = 0; tractionIdx < 3; ++tractionIdx) {
-                stressInfo_[globI][loc].traction[tractionIdx] = res[tractionIdx];
-            }
         }
     }
 
diff --git a/opm/models/tpsa/tpsanewtonmethod.hpp b/opm/models/tpsa/tpsanewtonmethod.hpp
index 387838a5c21..01c4ef41495 100644
--- a/opm/models/tpsa/tpsanewtonmethod.hpp
+++ b/opm/models/tpsa/tpsanewtonmethod.hpp
@@ -26,23 +26,18 @@
 #define TPSA_NEWTON_METHOD_HPP
 
 #include <opm/common/Exceptions.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/models/tpsa/tpsabaseproperties.hpp>
 #include <opm/models/tpsa/tpsanewtonmethodparams.hpp>
 #include <opm/models/utils/timer.hpp>
 #include <opm/models/utils/timerguard.hh>
 
-#include <opm/simulators/linalg/PropertyTree.hpp>
-
-#include <algorithm>
 #include <cmath>
 #include <cstddef>
-#include <exception>
-#include <iostream>
 #include <iomanip>
 #include <sstream>
 #include <string>
-#include <unistd.h>
 
 
 namespace Opm {
@@ -52,6 +47,10 @@ namespace Opm {
 *
 * Generates the Jacobian matrix, J(u^n) and residual vector, R(u^n), with a solution vector, u^n, at iteration n.
 * Subsequently the linear system J(u^n)\delta u^n = -R(u^n) is solved to get u^{n+1} = u^n + \Delta u^n.
+*
+* By default only a single iteration is done to solve the linear elasticity system. Note that,
+* in this case, only linear solver convergence is checked, and Newton error is only used to check if
+* the linearization step produced a low enough residual to jump out before a linear solve.
 */
 template <class TypeTag>
 class TpsaNewtonMethod
@@ -82,8 +81,10 @@ class TpsaNewtonMethod
         , linearSolver_(simulator)
         , error_(1e100)
         , lastError_(1e100)
+        , initialError_(1e100)
         , numIterations_(0)
         , numLinearizations_(0)
+        , numTotLinearIterations_(0)
     {
         // Read runtime/default Newton parameters
         params_.read();
@@ -154,7 +155,6 @@ class TpsaNewtonMethod
                 linearSolver_.getResidual(residual);
                 solveTimer_.stop();
 
-                // The preSolve_() method usually computes the errors, but it can do something else in addition.
                 // TODO: should its costs be counted to the linearization or to the update?
                 updateTimer_.start();
                 preSolve_(currentSolution, residual);
@@ -174,12 +174,13 @@ class TpsaNewtonMethod
                 solveTimer_.start();
                 solutionUpdate = 0.0;
                 const bool conv = linearSolver_.solve(solutionUpdate);
+                numTotLinearIterations_ += linearSolver_.iterations();
                 solveTimer_.stop();
 
                 if (!conv) {
                     solveTimer_.stop();
                     if (verbosity_() > 0) {
-                        std::cout << "TPSA: Linear solver did not converge!" << std::endl;
+                        OpmLog::warning("TPSA: Linear solver did not converge!");
                     }
 
                     prePostProcessTimer_.start();
@@ -203,8 +204,8 @@ class TpsaNewtonMethod
         catch (const Dune::Exception& e)
         {
             if (verbosity_() > 0) {
-                std::cout << "TPSA: Newton method caught exception: \""
-                          << e.what() << "\"\n" << std::flush;
+                OpmLog::error("TPSA: Newton method caught exception: \"" +
+                              std::string(e.what()) + "\"");
             }
 
             prePostProcessTimer_.start();
@@ -216,8 +217,8 @@ class TpsaNewtonMethod
         catch (const NumericalProblem& e)
         {
             if (verbosity_() > 0) {
-                std::cout << "TPSA: Newton method caught exception: \""
-                          << e.what() << "\"\n" << std::flush;
+                OpmLog::error("TPSA: Newton method caught exception: \"" +
+                              std::string(e.what()) + "\"");
             }
 
             prePostProcessTimer_.start();
@@ -229,34 +230,30 @@ class TpsaNewtonMethod
 
         // print the timing summary of the time step
         if (verbosity_() > 0) {
-            Scalar elapsedTot =
-                linearizeTimer_.realTimeElapsed() +
-                solveTimer_.realTimeElapsed() +
-                updateTimer_.realTimeElapsed();
-            const auto default_precision{std::cout.precision()};
-            std::cout << std::setprecision(2)
-                      << "TPSA: "
-                      << "Newton iter = " << numIterations() << " (error="
-                      << error_ << ") | "
-                      << "linearization = "
-                      << linearizeTimer_.realTimeElapsed() << "s ("
-                      << 100 * linearizeTimer_.realTimeElapsed() / elapsedTot << "%) | "
-                      << "solve = "
-                      << solveTimer_.realTimeElapsed() << "s ("
-                      << 100 * solveTimer_.realTimeElapsed() / elapsedTot << "%) | "
-                      << "update = "
-                      << updateTimer_.realTimeElapsed() << "s ("
-                      << 100 * updateTimer_.realTimeElapsed() / elapsedTot << "%)"
-                      << "\n" << std::flush;
-            std::cout << std::setprecision(default_precision); // restore default output width
+            std::ostringstream oss;
+            oss << std::setprecision(2)
+                << "TPSA: "
+                << "Newton iter = " << numIterations();
+            if (!singleIteration_()) {
+                oss << " (error=" << error_ << ")";
+            }
+            oss << " | "
+                << "Linearizations = "
+                << numLinearizations() << " ("
+                << linearizeTimer_.realTimeElapsed() << "s) | "
+                << "Linear iter = "
+                << numTotLinearIterations() << " ("
+                << solveTimer_.realTimeElapsed() << "s)";
+            OpmLog::info(oss.str());
         }
 
-        // if we're not converged, tell the implementation that we've failed
-        if (!converged()) {
+        // if we're not converged, tell the implementation that we've failed; ignored for
+        // max. iteration = 1.
+        if (!singleIteration_() && !converged()) {
             prePostProcessTimer_.start();
             failed_();
             if (verbosity_() > 0) {
-                std::cout << "TPSA: Newton iterations did not converge!" << std::endl;
+                OpmLog::warning("TPSA: Newton iterations did not converge!");
             }
             prePostProcessTimer_.stop();
             return false;
@@ -268,10 +265,29 @@ class TpsaNewtonMethod
     * \brief Returns true if the error of the solution is below the tolerance.
     *
     * \returns Bool indicating if convergence has been achived
+    *
+    * \note Decides whether another iteration is worth doing, but not whether apply() succeeded when
+    *       only a single iteration is done, see singleIteration_()
     */
     bool converged() const
     { return error_ <= tolerance(); }
 
+    /*!
+    * \brief Returns the error of the first linearization of the last apply().
+    *
+    * \returns Weighted maximum norm of the residual the Newton method started from
+    */
+    Scalar initialError() const
+    { return initialError_; }
+
+    /*!
+    * \brief Returns true if the state handed to the last apply() already solved the system.
+    *
+    * \returns Bool indicating if the first linearization was below the tolerance
+    */
+    bool initiallyConverged() const
+    { return initialError_ <= tolerance(); }
+
     /*!
     * \brief Returns a reference to the object describing the current physical problem.
     *
@@ -336,6 +352,14 @@ class TpsaNewtonMethod
     int numLinearizations() const
     { return numLinearizations_; }
 
+    /*!
+    * \brief Returns the number of linear solver iterations done since the Newton method was invoked.
+    *
+    * \returns Number of linear iterations, summed over the Newton iterations
+    */
+    int numTotLinearIterations() const
+    { return numTotLinearIterations_; }
+
     /*!
     * \brief Return the current tolerance at which the Newton method considers itself to be converged.
     *
@@ -395,6 +419,17 @@ class TpsaNewtonMethod
     int verbosity_() const
     { return simulator_.gridView().comm().rank() == 0 ? params_.verbosity_ : 0; }
 
+    /*!
+    * \brief Whether the Newton method is limited to a single iteration.
+    *
+    * \returns Bool indicating if at most one linearization and linear solve is done
+    *
+    * If the elasticity system is linear, a single solve is exact and more than one Newton error is
+    * unnecessary.
+    */
+    bool singleIteration_() const
+    { return params_.maxIterations_ <= 1; }
+
     /*!
     * \brief Called before the Newton method is applied to an non-linear system of equations.
     */
@@ -402,7 +437,9 @@ class TpsaNewtonMethod
     {
         numIterations_ = 0;
         numLinearizations_ = 0;
+        numTotLinearIterations_ = 0;
         error_ = 1e100;
+        initialError_ = 1e100;
     }
 
     /*!
@@ -449,6 +486,11 @@ class TpsaNewtonMethod
         // Take the other processes into account
         error_ = simulator_.gridView().comm().max(error_);
 
+        // Remember what the state handed to the Newton method was worth, see initiallyConverged()
+        if (numLinearizations_ == 1) {
+            initialError_ = error_;
+        }
+
         // Make sure that the error never grows beyond the maximum allowed one
         if (error_ > newtonMaxError) {
             throw NumericalProblem("TPSA: Newton error " + std::to_string(double(error_)) +
@@ -519,9 +561,8 @@ class TpsaNewtonMethod
 
         // Output error info
         if (verbosity_() > 1) {
-            std::cout << "TPSA: End Newton iteration " << numIterations_ << ""
-                      << " with error = " << error_
-                      << std::endl;
+            OpmLog::info("TPSA: End Newton iteration " + std::to_string(numIterations_) +
+                         " with error = " + std::to_string(error_));
         }
     }
 
@@ -532,6 +573,11 @@ class TpsaNewtonMethod
     */
     bool proceed_() const
     {
+        // Exactly one linearization and one linear solve, whatever the error did
+        if (singleIteration_()) {
+            return numIterations() < 1;
+        }
+
         if (numIterations() < params_.minIterations_) {
             return true;
         }
@@ -555,7 +601,7 @@ class TpsaNewtonMethod
     * \brief Called if the Newton method broke down.
     */
     void failed_()
-    { numIterations_ = params_.targetIterations_ * 2; }
+    { }
 
     Simulator& simulator_;
     LinearSolverBackend linearSolver_;
@@ -567,10 +613,12 @@ class TpsaNewtonMethod
 
     Scalar error_;
     Scalar lastError_;
+    Scalar initialError_;
     TpsaNewtonMethodParams<Scalar> params_;
 
     int numIterations_;
     int numLinearizations_;
+    int numTotLinearIterations_;
 };  // class TpsaNewtonMethod
 
 }  // namespace Opm
diff --git a/opm/models/tpsa/tpsanewtonmethodparams.cpp b/opm/models/tpsa/tpsanewtonmethodparams.cpp
index 789d796d605..ba5c1f57419 100644
--- a/opm/models/tpsa/tpsanewtonmethodparams.cpp
+++ b/opm/models/tpsa/tpsanewtonmethodparams.cpp
@@ -40,8 +40,6 @@ void TpsaNewtonMethodParams<Scalar>::registerParameters()
 {
     Parameters::Register<Parameters::TpsaNewtonVerbosity>
         ("Verbosity level of TPSA Newton solver: 0 (=none), 1 (=basic), 2 (=all)");
-    Parameters::Register<Parameters::TpsaNewtonTargetIterations>
-        ("The 'optimum' number of TPSA Newton iterations");
     Parameters::Register<Parameters::TpsaNewtonMaxIterations>
         ("The maximum number of TPSA Newton iterations");
     Parameters::Register<Parameters::TpsaNewtonMinIterations>
@@ -59,7 +57,6 @@ template<class Scalar>
 void TpsaNewtonMethodParams<Scalar>::read()
 {
     verbosity_ = Parameters::Get<Parameters::TpsaNewtonVerbosity>();
-    targetIterations_ = Parameters::Get<Parameters::TpsaNewtonTargetIterations>();
     minIterations_ = Parameters::Get<Parameters::TpsaNewtonMinIterations>();
     maxIterations_ = Parameters::Get<Parameters::TpsaNewtonMaxIterations>();
     tolerance_ = Parameters::Get<Parameters::TpsaNewtonTolerance<Scalar>>();
diff --git a/opm/models/tpsa/tpsanewtonmethodparams.hpp b/opm/models/tpsa/tpsanewtonmethodparams.hpp
index 9885b86bd05..35a01a48d7a 100644
--- a/opm/models/tpsa/tpsanewtonmethodparams.hpp
+++ b/opm/models/tpsa/tpsanewtonmethodparams.hpp
@@ -32,15 +32,12 @@ namespace Opm::Parameters {
 template<class Scalar>
 struct TpsaNewtonMaxError { static constexpr Scalar value = 1e100; };
 
-// Number of maximum iterations for the Newton method
-struct TpsaNewtonMaxIterations { static constexpr int value = 20; };
+// Number of maximum iterations for the Newton method.
+struct TpsaNewtonMaxIterations { static constexpr int value = 1; };
 
 // Number of minimum iterations for the Newton method
 struct TpsaNewtonMinIterations { static constexpr int value = 1; };
 
-// Target number of iterations
-struct TpsaNewtonTargetIterations { static constexpr int value = 10; };
-
 // Convergence tolerance
 template<class Scalar>
 struct TpsaNewtonTolerance { static constexpr Scalar value = 1e-3; };
@@ -63,7 +60,6 @@ struct TpsaNewtonMethodParams
 
     int verbosity_;
     bool writeConvergence_;
-    int targetIterations_;
     int minIterations_;
     int maxIterations_;
     Scalar tolerance_;
diff --git a/opm/simulators/flow/FlowProblemTPSA.hpp b/opm/simulators/flow/FlowProblemTPSA.hpp
index 16a577b893b..217e2f7b7c3 100644
--- a/opm/simulators/flow/FlowProblemTPSA.hpp
+++ b/opm/simulators/flow/FlowProblemTPSA.hpp
@@ -228,12 +228,13 @@ class FlowProblemTPSA : public FlowProblemBlackoil<TypeTag>
         Scalar avgSmodulus = 0.0;
         const auto& gridView = this->gridView();
         ElementContext elemCtx(this->simulator());
+        unsigned numDof = 0;
         for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
             elemCtx.updatePrimaryStencil(elem);
             int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
             avgSmodulus += this->shearModulus(elemIdx);
+            ++numDof;
         }
-        std::size_t numDof = this->model().numGridDof();
         const auto& comm = this->simulator().vanguard().grid().comm();
         avgSmodulus = comm.sum(avgSmodulus);
         Scalar numTotalDof = comm.sum(numDof);
@@ -265,6 +266,19 @@ class FlowProblemTPSA : public FlowProblemBlackoil<TypeTag>
         }
     }
 
+    /*!
+     * \brief Called by simulator at the end of each timestep
+     */
+    void endTimeStep() override
+    {
+        // Update info for mechanics output
+        // OBS: Must be done before ParentClass::endTimeStep!
+        geoMechModel().linearizer().updateStressInfo();
+
+        // Call parent class endTimeStep()
+        ParentType::endTimeStep();
+    }
+
     /*!
     * \brief Organize mechanics boundary conditions
     *
diff --git a/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp
index c4d33969828..38ee7c541bd 100644
--- a/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp
+++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoirTPSA.hpp
@@ -145,10 +145,10 @@ class NonlinearSystemBlackOilReservoirTPSA : public NonlinearSystemBlackOilReser
             ++seqIter_;
 
             // Fixed-stress convergence check:
-            // If the initial residual error, hence check for no. linearizations == 1, was small enough, we have
-            // convergence in the fixed-stress iterations
+            // If the initial residual error was small enough, we have convergence in the
+            // fixed-stress iterations
             if (tpsaConv
-                && this->simulator_.problem().geoMechModel().newtonMethod().numLinearizations() == 1
+                && this->simulator_.problem().geoMechModel().newtonMethod().initiallyConverged()
                 && seqIter_ >= minSeqIter) {
                 // Info
                 std::string msg = fmt::format("TPSA: Fixed-stress scheme converged in {} iterations", seqIter_);
diff --git a/opm/simulators/linalg/ISTLSolverTPSA.hpp b/opm/simulators/linalg/ISTLSolverTPSA.hpp
index 0333f98c9a6..c614e71f06d 100644
--- a/opm/simulators/linalg/ISTLSolverTPSA.hpp
+++ b/opm/simulators/linalg/ISTLSolverTPSA.hpp
@@ -49,7 +49,9 @@
 
 #include <algorithm>
 #include <any>
+#include <array>
 #include <cctype>
+#include <cmath>
 #include <cstddef>
 #include <functional>
 #include <memory>
@@ -57,6 +59,8 @@
 #include <string>
 #include <vector>
 
+#include <fmt/format.h>
+
 namespace Opm {
 
 /*!
@@ -105,6 +109,10 @@ class ISTLSolverTPSA :
           , matrix_(nullptr)
           , rhs_(nullptr)
     {
+        // Init. scaling factors
+        rowFactor_.fill(Scalar(1.0));
+        colFactor_.fill(Scalar(1.0));
+
         // Init parameters
         parameters_.init();
 
@@ -181,6 +189,10 @@ class ISTLSolverTPSA :
         if (isParallel() && type != "paroverilu0") {
             matrix_->makeOverlapRowsInvalid(overlapRows_);
         }
+
+        // Scale linear system (the right-hand side is scaled in solve())
+        computeScalingFactors_();
+        matrix_->scaleFields(rowFactor_, colFactor_);
     }
 
     /*!
@@ -219,6 +231,9 @@ class ISTLSolverTPSA :
 
         x = 0.0;
 
+        // Scale right-hand side before solve
+        rhs_->scaleFields(rowFactor_);
+
         // Solve linear system
         Dune::InverseOperatorResult result;
         assert(solver_);
@@ -227,6 +242,9 @@ class ISTLSolverTPSA :
         // Store no. linear iterations
         iterations_ = result.iterations;
 
+        // x solves the scaled system, so recover the update of the original one
+        x.scaleFields(colFactor_);
+
         // Return result for convergence check (boolean)
         return checkConvergence(result);
     }
@@ -330,6 +348,76 @@ class ISTLSolverTPSA :
 #endif
     }
 
+    /*!
+     * \brief Compute scaling factors for matrix and right-hand-side vector
+     */
+    void computeScalingFactors_()
+    {
+        rowFactor_.fill(Scalar(1.0));
+        colFactor_.fill(Scalar(1.0));
+
+        const auto& mode = parameters_.scale_linear_system_;
+        if (mode == "none") {
+            return;
+        }
+
+        if (mode != "eqweight" && mode != "user") {
+            OPM_THROW(std::runtime_error,
+                      fmt::format("TPSA: \"{}\" is not a valid setting for "
+                                  "--tpsa-scale-linear-system. Use none, eqweight or user.",
+                                  mode));
+        }
+
+        // Helper array for eq. indices to matrix/vetor field indices
+        static constexpr std::array<unsigned, Linear::numTpsaFields> fieldEq {
+            0,
+            1,
+            2,
+            3 * Linear::numDispDofs,
+            3 * Linear::numDispDofs + Linear::numRotDofs
+        };
+
+        // eqweight and user applied in the same manner, see
+        // FlowProblemTPSA::computeAndSetEqWeights_()
+        if (mode == "eqweight") {
+            const auto& model = simulator_.problem().geoMechModel();
+            for (std::size_t fieldIdx = 0; fieldIdx < Linear::numTpsaFields; ++fieldIdx) {
+                rowFactor_[fieldIdx] = model.eqWeight(/*dofIdx=*/0, fieldEq[fieldIdx]);
+                colFactor_[fieldIdx] = rowFactor_[fieldIdx];
+            }
+        }
+        else {
+            const Scalar factor = userFactor_();
+            for (std::size_t fieldIdx = 0; fieldIdx < Linear::numTpsaFields; ++fieldIdx) {
+                rowFactor_[fieldIdx] = (fieldEq[fieldIdx] < 3 * Linear::numDispDofs)
+                    ? Scalar(1.0) / factor
+                    : factor;
+                colFactor_[fieldIdx] = rowFactor_[fieldIdx];
+            }
+        }
+    }
+
+    /*!
+     * \brief The factor of the "user" field scaling with some checks
+     *
+     * \returns Factor
+     */
+    Scalar userFactor_() const
+    {
+        const Scalar factor = parameters_.scale_linear_system_factor_;
+
+        // The solution of the scaled system is recovered by multiplying by these
+        // factors, so anything that cannot be inverted is a mistake, not a scaling.
+        if (!std::isfinite(factor) || factor <= Scalar(0.0)) {
+            OPM_THROW(std::runtime_error,
+                      fmt::format("TPSA: --tpsa-scale-linear-system-factor must be finite "
+                                  "and positive, but {} was given",
+                                  factor));
+        }
+
+        return factor;
+    }
+
     /*!
     * \brief Check for linear solver convergence
     *
@@ -433,6 +521,10 @@ class ISTLSolverTPSA :
 
     SolverType* solver_ = nullptr;
     PrecondType* precond_ = nullptr;
+
+    // Row and column factors for scaling
+    std::array<Scalar, Linear::numTpsaFields> rowFactor_ {};
+    std::array<Scalar, Linear::numTpsaFields> colFactor_ {};
 };
 
 }  // namespace Opm
diff --git a/opm/simulators/linalg/TPSALinearSolverParameters.cpp b/opm/simulators/linalg/TPSALinearSolverParameters.cpp
index 03c4d08171e..efae27e931e 100644
--- a/opm/simulators/linalg/TPSALinearSolverParameters.cpp
+++ b/opm/simulators/linalg/TPSALinearSolverParameters.cpp
@@ -48,6 +48,8 @@ void TpsaLinearSolverParameters::init()
     ignoreConvergenceFailure_ = Parameters::Get<Parameters::TpsaLinearSolverIgnoreConvergenceFailure>();
     linsolver_ = Parameters::Get<Parameters::TpsaLinearSolver>();
     linear_solver_print_json_definition_ = Parameters::Get<Parameters::TpsaLinearSolverPrintJsonDefinition>();
+    scale_linear_system_ = Parameters::Get<Parameters::TpsaScaleLinearSystem>();
+    scale_linear_system_factor_ = Parameters::Get<Parameters::TpsaScaleLinearSystemFactor>();
 
     // Hardcode use of CPU linear solvers (?)
     linear_solver_accelerator_ = Parameters::LinearSolverAcceleratorType::CPU;
@@ -83,6 +85,12 @@ void TpsaLinearSolverParameters::registerParameters()
     Parameters::Register<Parameters::TpsaLinearSolverPrintJsonDefinition>
         ("Print JSON formatted configuration of the TPSA linear solver. Can be used to make configuration JSON file "
         "for --tpsa-linear-solver");
+    Parameters::Register<Parameters::TpsaScaleLinearSystem>
+        ("Scaling of the TPSA linear system. Valid options are: none, eqweight, "
+         "or user. The user option uses the value in --tpsa-scale-linear-system-factor.");
+    Parameters::Register<Parameters::TpsaScaleLinearSystemFactor>
+        ("Factor of the 'user' field scaling of the TPSA linear system;"
+         "factor = 1.0 means system is unscaled.");
 }
 
 /*!
@@ -92,8 +100,8 @@ void TpsaLinearSolverParameters::registerParameters()
 */
 void TpsaLinearSolverParameters::reset()
 {
-    linear_solver_reduction_ = 1e-3;
-    relaxed_linear_solver_reduction_ = 1e-3;
+    linear_solver_reduction_ = 1e-5;
+    relaxed_linear_solver_reduction_ = 1e-5;
     linear_solver_maxiter_ = 200;
     linear_solver_restart_ = 40;
     linear_solver_verbosity_ = 0;
@@ -101,8 +109,10 @@ void TpsaLinearSolverParameters::reset()
     ilu_fillin_level_ = 0;
     newton_use_gmres_ = false;
     ignoreConvergenceFailure_ = false;
-    linsolver_ = "hypre";
+    linsolver_ = "ilu0";
     linear_solver_print_json_definition_ = false;
+    scale_linear_system_ = "eqweight";
+    scale_linear_system_factor_ = 1.0;
 }
 
 }  // namespace Opm
diff --git a/opm/simulators/linalg/TPSALinearSolverParameters.hpp b/opm/simulators/linalg/TPSALinearSolverParameters.hpp
index 76802c0e34d..43bc33f9b14 100644
--- a/opm/simulators/linalg/TPSALinearSolverParameters.hpp
+++ b/opm/simulators/linalg/TPSALinearSolverParameters.hpp
@@ -33,8 +33,8 @@
 // Default runtime parameters
 namespace Opm::Parameters {
 
-struct TpsaLinearSolverReduction { static constexpr double value = 1e-3; };
-struct TpsaRelaxedLinearSolverReduction { static constexpr double value = 1e-3; };
+struct TpsaLinearSolverReduction { static constexpr double value = 1e-5; };
+struct TpsaRelaxedLinearSolverReduction { static constexpr double value = 1e-5; };
 struct TpsaLinearSolverMaxIter { static constexpr int value = 200; };
 struct TpsaLinearSolverRestart { static constexpr int value = 40; };
 struct TpsaLinearSolverVerbosity { static constexpr int value = 0; };
@@ -42,8 +42,10 @@ struct TpsaIluRelaxation { static constexpr double value = 0.9; };
 struct TpsaIluFillinLevel { static constexpr int value = 0; };
 struct TpsaUseGmres { static constexpr bool value = false; };
 struct TpsaLinearSolverIgnoreConvergenceFailure { static constexpr bool value = false; };
-struct TpsaLinearSolver { static constexpr auto value = "hypre"; };
+struct TpsaLinearSolver { static constexpr auto value = "ilu0"; };
 struct TpsaLinearSolverPrintJsonDefinition { static constexpr auto value = false; };
+struct TpsaScaleLinearSystem { static constexpr auto value = "eqweight"; };
+struct TpsaScaleLinearSystemFactor { static constexpr double value = 1.0; };
 
 }  // namespace Opm::Parameters
 
@@ -57,6 +59,9 @@ struct TpsaLinearSolverParameters : public FlowLinearSolverParameters
     void init();
     static void registerParameters();
     void reset();
+
+    std::string scale_linear_system_;
+    double scale_linear_system_factor_;
 };
 
 }  // namespace Opm
diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp
index 29047733943..97f421d33ce 100644
--- a/opm/simulators/linalg/setupPropertyTree.cpp
+++ b/opm/simulators/linalg/setupPropertyTree.cpp
@@ -477,17 +477,16 @@ setupTpsa(std::string conf, const FlowLinearSolverParameters& p)
         OPM_THROW(std::invalid_argument,
                   fmt::format(
                       "No valid settings found for --tpsa-linear-solver={}! Valid preset "
-                      "options are: default, hypre, ilu0, dilu or amg. Alternatively, give "
+                      "options are: hypre, ilu0, dilu or amg. Alternatively, give "
                       "the name of a .json file holding a full solver configuration.", conf)
             );
     }
 
 #if ! HAVE_HYPRE
     if (conf == "hypre") {
-        OpmLog::warning(
-            "--tpsa-linear-solver=hypre requires a build with Hypre support (USE_HYPRE=ON). "
-            "Switching to ilu0!");
-        conf = "ilu0"s;
+        OPM_THROW(std::invalid_argument,
+                  "--tpsa-linear-solver=hypre requires a build with Hypre support "
+                  "(USE_HYPRE=ON).");
     }
 #endif