diff --git a/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_damping.md b/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_damping.md
new file mode 100644
index 0000000000000000000000000000000000000000..e38215482d0860aca7b18434c2c5e6ec064ef0a1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/nonlinear_solvers/nonlinear_solver/t_damping.md
@@ -0,0 +1,5 @@
+A positive damping factor.
+
+The default value 1.0 gives a non-damped Newton method. Common values are in the
+range 0.5 to 0.7 for somewhat conservative method and seldom become smaller than
+0.2 for very conservative approach.
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 96c899ee32087815a2d4c4801c96c48805a9296d..ae2e52b30f0926274c9e1590b2fba45227925f74 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -232,7 +232,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
             auto& x_new =
                 NumLib::GlobalVectorProvider::provider.getVector(
                     x, _x_new_id);
-            LinAlg::axpy(x_new, -_alpha, minus_delta_x);
+            LinAlg::axpy(x_new, -_damping, minus_delta_x);
 
             if (postIterationCallback)
                 postIterationCallback(iteration, x_new);
@@ -321,10 +321,20 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
     }
     if (type == "Newton")
     {
+        //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__damping}
+        auto const damping = config.getConfigParameter<double>("damping", 1.0);
+        if (damping <= 0)
+        {
+            OGS_FATAL(
+                "The damping factor for the Newon method must be positive, got "
+                "%g.",
+                damping);
+        }
         auto const tag = NonlinearSolverTag::Newton;
         using ConcreteNLS = NonlinearSolver<tag>;
         return std::make_pair(
-            std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
+            std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
+            tag);
     }
     OGS_FATAL("Unsupported nonlinear solver type");
 }
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index c725de023aece2293cd7fc3f3cceff96a5057ea4..ba84bb5f3258b58a96eabd01811d06d4614fd6dc 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -86,12 +86,12 @@ public:
      * \param linear_solver the linear solver used by this nonlinear solver.
      * \param maxiter the maximum number of iterations used to solve the
      *                equation.
+     * \param damping \copydoc _damping
      */
-    explicit NonlinearSolver(
-        GlobalLinearSolver& linear_solver,
-        const unsigned maxiter)
-        : _linear_solver(linear_solver),
-          _maxiter(maxiter)
+    explicit NonlinearSolver(GlobalLinearSolver& linear_solver,
+                             const unsigned maxiter,
+                             double const damping = 1.0)
+        : _linear_solver(linear_solver), _maxiter(maxiter), _damping(damping)
     {
     }
 
@@ -117,8 +117,11 @@ private:
     ConvergenceCriterion* _convergence_criterion = nullptr;
     const unsigned _maxiter;  //!< maximum number of iterations
 
-    double const _alpha =
-        1;  //!< Damping factor. \todo Add constructor parameter.
+    //! A positive damping factor. The default value 1.0 gives a non-damped
+    //! Newton method. Common values are in the range 0.5 to 0.7 for somewhat
+    //! conservative method and seldom become smaller than 0.2 for very
+    //! conservative approach.
+    double const _damping;
 
     std::size_t _res_id = 0u;            //!< ID of the residual vector.
     std::size_t _J_id = 0u;              //!< ID of the Jacobian matrix.