diff --git a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
index bc3247f59dc78bde2ee78711947cdba56b9c8677..e41f4ca124d3522c6f83bcf7391c1b73fddba99b 100644
--- a/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
@@ -19,6 +19,7 @@
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/CreateNewtonRaphsonSolverParameters.h"
 #include "ParameterLib/Utils.h"
+#include "ProcessLib/Common/HydroMechanics/CreateInitialStress.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 #include "RichardsMechanicsProcess.h"
@@ -166,12 +167,9 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
     DBUG("Media properties verified.");
 
     // Initial stress conditions
-    auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
-        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__initial_stress}
-        config, "initial_stress", parameters,
-        // Symmetric tensor size, 4 or 6, not a Kelvin vector.
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
-        &mesh);
+    auto const initial_stress =
+        ProcessLib::createInitialStress<DisplacementDim>(config, parameters,
+                                                         mesh);
 
     std::optional<MicroPorosityParameters> micro_porosity_parameters;
     if (auto const micro_porosity_config =
diff --git a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
index fd8d870dafed3d7d40290fd759ed0e53806d1c93..a5d6b6b2d17fe499302f00a480568d6c6ca08401 100644
--- a/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
+++ b/ProcessLib/RichardsMechanics/LocalAssemblerInterface.h
@@ -76,13 +76,13 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
                 element_.getID());
         }
 
-        if (name == "sigma" && process_data_.initial_stress != nullptr)
+        if (name == "sigma" && process_data_.initial_stress.value != 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);
+                process_data_.initial_stress.value->name);
         }
 
         if (name.starts_with("material_state_variable_"))
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 84be48989be244f3402b24a0fc6a967262893bd8..bf2798ca7b07b9eee4f346f6646d6eb4ccf0a64e 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -219,6 +219,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     auto const& solid_phase = medium->phase("Solid");
 
+    auto const& identity2 = MathLib::KelvinVector::Invariants<
+        MathLib::KelvinVector::kelvin_vector_dimensions(
+            DisplacementDim)>::identity2;
+
     unsigned const n_integration_points =
         this->integration_method_.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -257,6 +261,38 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         S_L_prev = medium->property(MPL::PropertyType::saturation)
                        .template value<double>(variables, x_position, t, dt);
 
+        if (this->process_data_.initial_stress.isTotalStress())
+        {
+            auto const alpha_b =
+                medium->property(MPL::PropertyType::biot_coefficient)
+                    .template value<double>(variables, x_position, t, dt);
+
+            variables.liquid_saturation = S_L_prev;
+            double const chi_S_L =
+                medium->property(MPL::PropertyType::bishops_effective_stress)
+                    .template value<double>(variables, x_position, t, dt);
+
+            // Initial stresses are total stress, which were assigned to
+            // sigma_eff in
+            // RichardsMechanicsLocalAssembler::initializeConcrete().
+            auto& sigma_eff =
+                std::get<ProcessLib::ThermoRichardsMechanics::
+                             ConstitutiveStress_StrainTemperature::
+                                 EffectiveStressData<DisplacementDim>>(
+                    this->current_states_[ip]);
+
+            auto& sigma_eff_prev = std::get<
+                PrevState<ProcessLib::ThermoRichardsMechanics::
+                              ConstitutiveStress_StrainTemperature::
+                                  EffectiveStressData<DisplacementDim>>>(
+                this->prev_states_[ip]);
+
+            // Reset sigma_eff to effective stress
+            sigma_eff.sigma_eff.noalias() +=
+                chi_S_L * alpha_b * (-p_cap_ip) * identity2;
+            sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
+        }
+
         if (medium->hasProperty(MPL::PropertyType::saturation_micro))
         {
             MPL::VariableArray vars;
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index 4cec7a14c8e6836cc0867b6d09d585cd295c3056..559c902bfc98cf3ea7bdfb4af22018d155b8b7ca 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -131,17 +131,20 @@ public:
                                                                 ip_data.N_u))};
 
             /// Set initial stress from parameter.
-            if (this->process_data_.initial_stress != nullptr)
+            if (this->process_data_.initial_stress.value)
             {
                 std::get<ProcessLib::ThermoRichardsMechanics::
                              ConstitutiveStress_StrainTemperature::
                                  EffectiveStressData<DisplacementDim>>(SD)
                     .sigma_eff =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
-                        DisplacementDim>((*this->process_data_.initial_stress)(
-                        std::numeric_limits<
-                            double>::quiet_NaN() /* time independent */,
-                        x_position));
+                        DisplacementDim>(
+                        // The data in process_data_.initial_stress.value can
+                        // be total stress or effective stress.
+                        (*this->process_data_.initial_stress.value)(
+                            std::numeric_limits<
+                                double>::quiet_NaN() /* time independent */,
+                            x_position));
             }
 
             double const t = 0;  // TODO (naumov) pass t from top
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 5ad61ae913f91b159120cdff8239cc24c70601a2..d237ed6c01cc8fb630097a2e831752fb628232e0 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -17,6 +17,7 @@
 #include "ComputeMicroPorosity.h"
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "ParameterLib/Parameter.h"
+#include "ProcessLib/Common/HydroMechanics/InitialStress.h"
 
 namespace MaterialLib
 {
@@ -42,9 +43,7 @@ struct RichardsMechanicsProcessData
                       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;
+    InitialStress const initial_stress;
 
     /// Specific body forces applied to solid and fluid.
     /// It is usually used to apply gravitational forces.