diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 26d0b2c431c6355e2d2ceb01b486ab3978ef2c57..2dca2852594c06020d50c4390a5820eb60c77710 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -496,23 +496,30 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
                 Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
             solver_ = std::make_unique<
                 details::EigenDirectLinearSolver<SolverType>>();
+            can_solve_rectangular_ = false;
             return;
         }
         case EigenOption::SolverType::BiCGSTAB:
         case EigenOption::SolverType::BiCGSTABL:
         case EigenOption::SolverType::CG:
-        case EigenOption::SolverType::LeastSquareCG:
         case EigenOption::SolverType::GMRES:
         case EigenOption::SolverType::IDRS:
         case EigenOption::SolverType::IDRSTABL:
             solver_ = details::createIterativeSolver(option_.solver_type,
                                                      option_.precon_type);
+            can_solve_rectangular_ = false;
+            return;
+        case EigenOption::SolverType::LeastSquareCG:
+            solver_ = details::createIterativeSolver(option_.solver_type,
+                                                     option_.precon_type);
+            can_solve_rectangular_ = true;
             return;
         case EigenOption::SolverType::PardisoLU:
         {
 #ifdef USE_MKL
             using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
             solver_.reset(new details::EigenDirectLinearSolver<SolverType>);
+            can_solve_rectangular_ = false;
             return;
 #else
             OGS_FATAL(
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
index 1955f1fd8495e3609af9556119b066bc3e725887..2f58fc83738f5137132abb4778f6b7cb7e615b19 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.h
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
@@ -71,9 +71,13 @@ public:
                MathLib::LinearSolverBehaviour const linear_solver_behaviour =
                    MathLib::LinearSolverBehaviour::RECOMPUTE);
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return can_solve_rectangular_; }
+
 protected:
     EigenOption option_;
     std::unique_ptr<EigenLinearSolverBase> solver_;
+    bool can_solve_rectangular_ = false;
     void setRestart();
     void setL();
     void setS();
diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
index 23763277d0a9241931e0537cb5330f1eb52f5654..3543d65533131c26542f07f35a48713e88a9368e 100644
--- a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
+++ b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
@@ -51,6 +51,9 @@ public:
 
     bool solve(EigenMatrix& A, EigenVector& b, EigenVector& x);
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return false; }
+
 private:
     bool solve(LisMatrix& A, LisVector& b, LisVector& x);
     std::string lis_options_;
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
index 7310848f3ea60f118df33340b5dca8ea1ed13060..1fc14a8bd588b019e9d3dbf442787d5759dd2154 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -40,6 +40,10 @@ PETScLinearSolver::PETScLinearSolver(std::string const& prefix,
 
     KSPSetInitialGuessNonzero(solver_, PETSC_TRUE);
     KSPSetFromOptions(solver_);  // set run-time options
+
+    KSPType ksp_type;
+    KSPGetType(solver_, &ksp_type);
+    can_solve_rectangular_ = canSolverHandleRectangular(ksp_type);
 }
 
 bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector& b, PETScVector& x)
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.h b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
index 0dfcbc359560f5945d8996c3b1ac0ff1a04ac58c..59c1ad0d8d2ba9d16a41aa6a339b7ee9d7512287 100644
--- a/MathLib/LinAlg/PETSc/PETScLinearSolver.h
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
@@ -18,6 +18,8 @@
 
 #include <petscksp.h>
 
+#include <algorithm>
+#include <array>
 #include <string>
 
 #include "PETScMatrix.h"
@@ -62,10 +64,26 @@ public:
     /// Get elapsed wall clock time.
     double getElapsedTime() const { return elapsed_ctime_; }
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return can_solve_rectangular_; }
+
 private:
+    bool canSolverHandleRectangular(std::string_view ksp_type)
+    {
+        // List of KSP types that can handle rectangular matrices
+        constexpr auto rectangular_solvers = std::to_array<std::string_view>(
+            {KSPGMRES, KSPFGMRES, KSPBCGS, KSPBCGSL, KSPTFQMR, KSPCGS});
+
+        return std::ranges::any_of(rectangular_solvers,
+                                   [&](const auto& solver)
+                                   { return solver == ksp_type; });
+    }
+
     KSP solver_;  ///< Solver type.
     PC pc_;       ///< Preconditioner type.
 
+    bool can_solve_rectangular_ = false;
+
     double elapsed_ctime_ = 0.0;  ///< Clock time
 };
 
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index eb533380fd16d64f63b39752765ace7c16bba02c..516ffa207d8dced39eeda68dc48909cd6a7ffb3a 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -170,8 +170,18 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         sys.getRhs(*x_prev[process_id], rhs);
 
         // Normalize the linear equation system, if required
-        if (sys.requiresNormalization())
+        if (sys.requiresNormalization() &&
+            !_linear_solver.canSolveRectangular())
+        {
             sys.getAandRhsNormalized(A, rhs);
+            WARN(
+                "The equation system is rectangular, but the current linear "
+                "solver only supports square systems. "
+                "The system will be normalized, which lead to a squared "
+                "condition number and potential numerical issues. "
+                "It is recommended to use a solver that supports rectangular "
+                "equation systems for better numerical stability.");
+        }
 
         INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());