diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md b/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_non_equilibrium_initial_residuum.md
similarity index 59%
rename from Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md
rename to Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_non_equilibrium_initial_residuum.md
index 96a72b7a0b46e99d4abc8229254c406f4363ca1b..b5e74bf01fb56a9973ad96c8a38cbe25978d849e 100644
--- a/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_initial_out_of_balance_forces.md
+++ b/Documentation/ProjectFile/prj/time_loop/processes/process/t_compensate_non_equilibrium_initial_residuum.md
@@ -1 +1 @@
-\copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_initial_out_of_balance_forces
+\copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_non_equilibrium_initial_residuum
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 59f406cf3b003bec73effc92a1660f9dd98e6870..0e97ab08bfb99c85c7bad85bed5eb1335471c7c0 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -27,11 +27,14 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
     _equation_system->assemble(x, process_id);
 }
 
-void NonlinearSolver<NonlinearSolverTag::Picard>::calculateOutOfBalanceForces(
-    std::vector<GlobalVector*> const& x, int const process_id)
+void NonlinearSolver<NonlinearSolverTag::Picard>::
+    calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
+                                           int const process_id)
 {
-    if (!_compensate_initial_out_of_balance_forces)
+    if (!_compensate_non_equilibrium_initial_residuum)
+    {
         return;
+    }
 
     auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
     auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
@@ -39,10 +42,10 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::calculateOutOfBalanceForces(
     _equation_system->getA(A);
     _equation_system->getRhs(rhs);
 
-    // F_oob = A * x - rhs
-    _f_oob = &NumLib::GlobalVectorProvider::provider.getVector();
-    MathLib::LinAlg::matMult(A, *x[process_id], *_f_oob);
-    MathLib::LinAlg::axpy(*_f_oob, -1.0, rhs);  // res -= rhs
+    // r_neq = A * x - rhs
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
+    MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs);  // res -= rhs
 }
 
 NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
@@ -91,10 +94,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         sys.getRhs(rhs);
         INFO("[time] Assembly took %g s.", time_assembly.elapsed());
 
-
-        // Apply out-of-balance forces if set.
-        if (_f_oob != nullptr)
-            LinAlg::axpy(rhs, -1, *_f_oob);
+        // Subract non-equilibrium initial residuum if set
+        if (_r_neq != nullptr)
+        {
+            LinAlg::axpy(rhs, -1, *_r_neq);
+        }
 
         timer_dirichlet.start();
         sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
@@ -214,15 +218,18 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
     //      equation every time and could not forget it.
 }
 
-void NonlinearSolver<NonlinearSolverTag::Newton>::calculateOutOfBalanceForces(
-    std::vector<GlobalVector*> const& x, int const process_id)
+void NonlinearSolver<NonlinearSolverTag::Newton>::
+    calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
+                                           int const process_id)
 {
-    if (!_compensate_initial_out_of_balance_forces)
+    if (!_compensate_non_equilibrium_initial_residuum)
+    {
         return;
+    }
 
     _equation_system->assemble(x, process_id);
-    _f_oob = &NumLib::GlobalVectorProvider::provider.getVector();
-    _equation_system->getResidual(*x[process_id], *_f_oob);
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    _equation_system->getResidual(*x[process_id], *_r_neq);
 }
 
 NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
@@ -273,9 +280,9 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
         sys.getJacobian(J);
         INFO("[time] Assembly took %g s.", time_assembly.elapsed());
 
-        // Apply out-of-balance forces if set.
-        if (_f_oob != nullptr)
-            LinAlg::axpy(res, -1, *_f_oob);
+        // Subract non-equilibrium initial residuum if set
+        if (_r_neq != nullptr)
+            LinAlg::axpy(res, -1, *_r_neq);
 
         minus_delta_x.setZero();
 
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index dcb1b5c754393ca2e466e2de843fcc57b5635f4f..83303f469afdbacdcac64193a6c99c83bdd1d985 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -49,7 +49,7 @@ public:
     virtual void assemble(std::vector<GlobalVector*> const& x,
                           int const process_id) const = 0;
 
-    virtual void calculateOutOfBalanceForces(
+    virtual void calculateNonEquilibriumInitialResiduum(
         std::vector<GlobalVector*> const& x, int const process_id) = 0;
 
     /*! Assemble and solve the equation system.
@@ -116,8 +116,8 @@ public:
     void assemble(std::vector<GlobalVector*> const& x,
                   int const process_id) const override;
 
-    void calculateOutOfBalanceForces(std::vector<GlobalVector*> const& x,
-                                     int const process_id) override;
+    void calculateNonEquilibriumInitialResiduum(
+        std::vector<GlobalVector*> const& x, int const process_id) override;
 
     NonlinearSolverStatus solve(
         std::vector<GlobalVector*>& x,
@@ -125,9 +125,9 @@ public:
             postIterationCallback,
         int const process_id) override;
 
-    void compensateInitialOutOfBalanceForces(bool const value)
+    void compensateNonEquilibriumInitialResiduum(bool const value)
     {
-        _compensate_initial_out_of_balance_forces = value;
+        _compensate_non_equilibrium_initial_residuum = value;
     }
 
 private:
@@ -144,18 +144,19 @@ private:
     //! conservative approach.
     double const _damping;
 
-    GlobalVector* _f_oob = nullptr;      //!< Out-of-balance forces vector.
+    GlobalVector* _r_neq = nullptr;      //!< non-equilibrium initial residuum.
     std::size_t _res_id = 0u;            //!< ID of the residual vector.
     std::size_t _J_id = 0u;              //!< ID of the Jacobian matrix.
     std::size_t _minus_delta_x_id = 0u;  //!< ID of the \f$ -\Delta x\f$ vector.
     std::size_t _x_new_id =
         0u;  //!< ID of the vector storing \f$ x - (-\Delta x) \f$.
 
-    /// Enables computation of the out-of-balance forces \f$ F_{\rm oob} \f$
-    /// before the first time step. The forces are zero if the external forces
-    /// are in equilibrium with the initial state/initial conditions. During the
-    /// simulation the new residuum reads \f$ \tilde r = r - F_{\rm oob} \f$.
-    bool _compensate_initial_out_of_balance_forces = false;
+    /// Enables computation of the non-equilibrium initial residuum \f$ r_{\rm
+    /// neq} \f$ before the first time step. The forces are zero if the external
+    /// forces are in equilibrium with the initial state/initial conditions.
+    /// During the simulation the new residuum reads \f$ \tilde r = r - r_{\rm
+    /// neq} \f$.
+    bool _compensate_non_equilibrium_initial_residuum = false;
 };
 
 /*! Find a solution to a nonlinear equation using the Picard fixpoint iteration
@@ -193,8 +194,8 @@ public:
     void assemble(std::vector<GlobalVector*> const& x,
                   int const process_id) const override;
 
-    void calculateOutOfBalanceForces(std::vector<GlobalVector*> const& x,
-                                     int const process_id) override;
+    void calculateNonEquilibriumInitialResiduum(
+        std::vector<GlobalVector*> const& x, int const process_id) override;
 
     NonlinearSolverStatus solve(
         std::vector<GlobalVector*>& x,
@@ -202,9 +203,9 @@ public:
             postIterationCallback,
         int const process_id) override;
 
-    void compensateInitialOutOfBalanceForces(bool const value)
+    void compensateNonEquilibriumInitialResiduum(bool const value)
     {
-        _compensate_initial_out_of_balance_forces = value;
+        _compensate_non_equilibrium_initial_residuum = value;
     }
 
 private:
@@ -215,15 +216,15 @@ private:
     ConvergenceCriterion* _convergence_criterion = nullptr;
     const int _maxiter;  //!< maximum number of iterations
 
-    GlobalVector* _f_oob = nullptr;  //!< Out-of-balance forces vector.
+    GlobalVector* _r_neq = nullptr;  //!< non-equilibrium initial residuum.
     std::size_t _A_id = 0u;      //!< ID of the \f$ A \f$ matrix.
     std::size_t _rhs_id = 0u;    //!< ID of the right-hand side vector.
     std::size_t _x_new_id = 0u;  //!< ID of the vector storing the solution of
                                  //! the linearized equation.
 
     /// \copydoc
-    /// NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_initial_out_of_balance_forces
-    bool _compensate_initial_out_of_balance_forces = false;
+    /// NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_non_equilibrium_initial_residuum
+    bool _compensate_non_equilibrium_initial_residuum = false;
 };
 
 /*! Creates a new nonlinear solver from the given configuration.
diff --git a/ProcessLib/CreateProcessData.cpp b/ProcessLib/CreateProcessData.cpp
index 0c7955c5198b11fc056610e50ec5829ea5e30a52..6d959c9086db6a69e7bf168c8f80a7795c119e74 100644
--- a/ProcessLib/CreateProcessData.cpp
+++ b/ProcessLib/CreateProcessData.cpp
@@ -24,7 +24,7 @@ static std::unique_ptr<ProcessData> makeProcessData(
     Process& process,
     std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
     std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit,
-    bool const compensate_initial_out_of_balance_forces)
+    bool const compensate_non_equilibrium_initial_residuum)
 {
     using Tag = NumLib::NonlinearSolverTag;
 
@@ -32,8 +32,8 @@ static std::unique_ptr<ProcessData> makeProcessData(
             dynamic_cast<NumLib::NonlinearSolver<Tag::Picard>*>(
                 &nonlinear_solver))
     {
-        nonlinear_solver_picard->compensateInitialOutOfBalanceForces(
-            compensate_initial_out_of_balance_forces);
+        nonlinear_solver_picard->compensateNonEquilibriumInitialResiduum(
+            compensate_non_equilibrium_initial_residuum);
         return std::make_unique<ProcessData>(
             std::move(timestepper), *nonlinear_solver_picard,
             std::move(conv_crit), std::move(time_disc), process_id, process);
@@ -42,8 +42,8 @@ static std::unique_ptr<ProcessData> makeProcessData(
             dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
                 &nonlinear_solver))
     {
-        nonlinear_solver_newton->compensateInitialOutOfBalanceForces(
-            compensate_initial_out_of_balance_forces);
+        nonlinear_solver_newton->compensateNonEquilibriumInitialResiduum(
+            compensate_non_equilibrium_initial_residuum);
         return std::make_unique<ProcessData>(
             std::move(timestepper), *nonlinear_solver_newton,
             std::move(conv_crit), std::move(time_disc), process_id, process);
@@ -92,10 +92,10 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
             //! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
             pcs_config.getConfigSubtree("convergence_criterion"));
 
-        auto const compensate_initial_out_of_balance_forces =
-            //! \ogs_file_param{prj__time_loop__processes__process__compensate_initial_out_of_balance_forces}
+        auto const compensate_non_equilibrium_initial_residuum =
+            //! \ogs_file_param{prj__time_loop__processes__process__compensate_non_equilibrium_initial_residuum}
             pcs_config.getConfigParameter<bool>(
-                "compensate_initial_out_of_balance_forces", false);
+                "compensate_non_equilibrium_initial_residuum", false);
 
         //! \ogs_file_param{prj__time_loop__processes__process__output}
         auto output = pcs_config.getConfigSubtreeOptional("output");
@@ -114,7 +114,7 @@ std::vector<std::unique_ptr<ProcessData>> createPerProcessData(
         per_process_data.emplace_back(
             makeProcessData(std::move(timestepper), nl_slv, process_id, pcs,
                             std::move(time_disc), std::move(conv_crit),
-                            compensate_initial_out_of_balance_forces));
+                            compensate_non_equilibrium_initial_residuum));
         ++process_id;
     }
 
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index af939682793e5554f421ef6ac39f376220bf5b2f..b2dbfedd2ac35531ab4ec89f00423c5055ce6f76 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -167,11 +167,12 @@ std::vector<GlobalVector*> setInitialConditions(
     return process_solutions;
 }
 
-void calculateOutOfBalanceForces(
+void calculateNonEquilibriumInitialResiduum(
     std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
-    std::vector<GlobalVector*> process_solutions)
+    std::vector<GlobalVector*>
+        process_solutions)
 {
-    INFO("Calculate out-of-balance forces.");
+    INFO("Calculate non-equilibrium initial residuum.");
     for (auto& process_data : per_process_data)
     {
         auto& ode_sys = *process_data->tdisc_ode_sys;
@@ -187,8 +188,8 @@ void calculateOutOfBalanceForces(
         double const t = 0;
         double const dt = 1;
         time_disc.nextTimestep(t, dt);
-        nonlinear_solver.calculateOutOfBalanceForces(process_solutions,
-                                                     process_data->process_id);
+        nonlinear_solver.calculateNonEquilibriumInitialResiduum(
+            process_solutions, process_data->process_id);
     }
 }
 
@@ -470,7 +471,7 @@ bool TimeLoop::loop()
     bool const is_staggered_coupling =
         !isMonolithicProcess(*_per_process_data[0]);
 
-    bool out_of_balance_forces_computed = false;
+    bool non_equilibrium_initial_residuum_computed = false;
     double t = _start_time;
     std::size_t accepted_steps = 0;
     std::size_t rejected_steps = 0;
@@ -498,10 +499,11 @@ bool TimeLoop::loop()
                 t, process_data->process_id);
         }
 
-        if (!out_of_balance_forces_computed)
+        if (!non_equilibrium_initial_residuum_computed)
         {
-            calculateOutOfBalanceForces(_per_process_data, _process_solutions);
-            out_of_balance_forces_computed = true;
+            calculateNonEquilibriumInitialResiduum(_per_process_data,
+                                                   _process_solutions);
+            non_equilibrium_initial_residuum_computed = true;
         }
 
         if (is_staggered_coupling)
diff --git a/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj b/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
index 614368396cbd9402ea42c9e35e35b319c8c25915..a7a0c9682a182713b304ebfbe021fa1b414336a6 100644
--- a/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
+++ b/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
@@ -33,7 +33,7 @@
         <processes>
             <process ref="SD">
                 <nonlinear_solver>basic_newton</nonlinear_solver>
-                <compensate_initial_out_of_balance_forces>true</compensate_initial_out_of_balance_forces>
+                <compensate_non_equilibrium_initial_residuum>true</compensate_non_equilibrium_initial_residuum>
                 <convergence_criterion>
                     <type>DeltaX</type>
                     <norm_type>NORM2</norm_type>
diff --git a/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj b/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
index 79bdbc0c4ae9858d18af4678fa8fcb22dae18189..aec9f3636171f33bced373c0e0696df510c31959 100644
--- a/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
+++ b/Tests/Data/ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj
@@ -37,7 +37,7 @@
         <processes>
             <process ref="TM">
                 <nonlinear_solver>basic_newton</nonlinear_solver>
-                <compensate_initial_out_of_balance_forces>true</compensate_initial_out_of_balance_forces>
+                <compensate_non_equilibrium_initial_residuum>true</compensate_non_equilibrium_initial_residuum>
                 <convergence_criterion>
                     <type>DeltaX</type>
                     <norm_type>NORM2</norm_type>