diff --git a/Documentation/ProjectFile/prj/processes/process/SMALL_DEFORMATION/t_initial_stress.md b/Documentation/ProjectFile/prj/processes/process/SMALL_DEFORMATION/t_initial_stress.md
new file mode 100644
index 0000000000000000000000000000000000000000..9832658f886be42d2bcca2c5134338db7402ffea
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/SMALL_DEFORMATION/t_initial_stress.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::SmallDeformation::SmallDeformationProcessData::initial_stress
diff --git a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
index 2c16a7b9f2ae26e634d6fb75a38eb9198b1687c4..24d6e9a4083629372d60fc3ffb7662a0e3df4b2f 100644
--- a/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
+++ b/ProcessLib/SmallDeformation/CreateSmallDeformationProcess.cpp
@@ -12,6 +12,7 @@
 #include <cassert>
 
 #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MathLib/KelvinVector.h"
 #include "ParameterLib/Utils.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
@@ -102,9 +103,30 @@ std::unique_ptr<Process> createSmallDeformationProcess(
         config.getConfigParameter<double>(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
+    // Initial stress conditions
+    ParameterLib::Parameter<double> const* initial_stress = nullptr;
+    auto const initial_stress_parameter_name =
+        //! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__initial_stress}
+        config.getConfigParameterOptional<std::string>("initial_stress");
+    if (initial_stress_parameter_name)
+    {
+        initial_stress = ParameterLib::findParameterOptional<double>(
+            *initial_stress_parameter_name, parameters,
+            // Symmetric tensor size, 4 or 6, not a Kelvin vector.
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value,
+            &mesh);
+    }
+    if (initial_stress != nullptr)
+    {
+        DBUG("Use '%s' as initial stress parameter.",
+             initial_stress->name.c_str());
+    }
+
     SmallDeformationProcessData<DisplacementDim> process_data{
-        materialIDs(mesh), std::move(solid_constitutive_relations),
-        solid_density, specific_body_force, reference_temperature};
+        materialIDs(mesh),   std::move(solid_constitutive_relations),
+        initial_stress,      solid_density,
+        specific_body_force, reference_temperature};
 
     SecondaryVariableCollection secondary_variables;
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 7d195f6d899cab7a26abfe2b87ad72b2b9aedc4e..7f0a9110aae5c9f50c682a7220e69110195d908e 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -176,6 +176,14 @@ public:
 
         if (name == "sigma_ip")
         {
+            if (_process_data.initial_stress != nullptr)
+            {
+                OGS_FATAL(
+                    "Setting initial conditions for stress from integration "
+                    "point data and from a parameter '%s' is not possible "
+                    "simultaneously.",
+                    _process_data.initial_stress->name.c_str());
+            }
             return setSigma(values);
         }
 
@@ -186,10 +194,28 @@ public:
     {
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
-
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            _ip_data[ip].pushBackState();
+            auto& ip_data = _ip_data[ip];
+
+            /// Set initial stress from parameter.
+            if (_process_data.initial_stress != nullptr)
+            {
+                ParameterLib::SpatialPosition const x_position{
+                    boost::none, _element.getID(), ip,
+                    MathLib::Point3d(interpolateCoordinates<ShapeFunction,
+                                                            ShapeMatricesType>(
+                        _element, ip_data.N))};
+
+                ip_data.sigma =
+                    MathLib::KelvinVector::symmetricTensorToKelvinVector<
+                        DisplacementDim>((*_process_data.initial_stress)(
+                        std::numeric_limits<
+                            double>::quiet_NaN() /* time independent */,
+                        x_position));
+            }
+
+            ip_data.pushBackState();
         }
     }
 
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcessData.h b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
index d48be7dabf7b7d2f0f2627c3a87cccf02e9d0b47..25864000a6eb554f6b23e5c5d041978e0b0d3a45 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcessData.h
@@ -37,6 +37,11 @@ struct SmallDeformationProcessData
         int,
         std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
         solid_materials;
+
+    /// Optional, initial stress field. A symmetric tensor, short vector
+    /// representation of length 4 or 6, ParameterLib::Parameter<double>.
+    ParameterLib::Parameter<double> const* const initial_stress;
+
     /// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
     ParameterLib::Parameter<double> const& solid_density;
 
diff --git a/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj b/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
index d1c2b783f2c1c2a37bf08af422e49796712a7c72..614368396cbd9402ea42c9e35e35b319c8c25915 100644
--- a/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
+++ b/Tests/Data/Mechanics/InitialStates/non_equilibrium_initial_state.prj
@@ -26,15 +26,14 @@
                 <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
                 <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
             </secondary_variables>
-            <nonequilibrium_state_variables>
-                <stress>nonequilibrium_stress</stress>
-            </nonequilibrium_state_variables>
+            <initial_stress>nonequilibrium_stress</initial_stress>
         </process>
     </processes>
     <time_loop>
         <processes>
             <process ref="SD">
                 <nonlinear_solver>basic_newton</nonlinear_solver>
+                <compensate_initial_out_of_balance_forces>true</compensate_initial_out_of_balance_forces>
                 <convergence_criterion>
                     <type>DeltaX</type>
                     <norm_type>NORM2</norm_type>