diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index a283e400a3cfb627616d3827ac2d8dcf8c4a743a..0f0126813e34a50d8e702d57cf6dd9cc619697d5 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -910,6 +910,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         _integration_method.getNumberOfPoints();
 
     double saturation_avg = 0;
+    double porosity_avg = 0;
 
     using KV = MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
     KV sigma_avg = KV::Zero();
@@ -936,12 +937,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         S_L = medium->property(MPL::PropertyType::saturation)
                   .template value<double>(variables, x_position, t, dt);
         saturation_avg += S_L;
+        porosity_avg +=
+            _ip_data[ip].porosity;  // Note, this is not updated, because needs
+                                    // xdot and dt to be passed.
         sigma_avg += _ip_data[ip].sigma_eff;
     }
     saturation_avg /= n_integration_points;
+    porosity_avg /= n_integration_points;
     sigma_avg /= n_integration_points;
 
     (*_process_data.element_saturation)[_element.getID()] = saturation_avg;
+    (*_process_data.element_porosity)[_element.getID()] = porosity_avg;
 
     Eigen::Map<KV>(&(*_process_data.element_stresses)[_element.getID() *
                                                       KV::RowsAtCompileTime]) =
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index bfbf25f5d340a79a71f77d7a0b5f6359942b5699..c05497d0809f20d9153040cd3fd4a2e68909fdce 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -206,6 +206,10 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
         const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
         MeshLib::MeshItemType::Cell, 1);
 
+    _process_data.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
+        MeshLib::MeshItemType::Cell, 1);
+
     _process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
         const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
         MeshLib::MeshItemType::Cell,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
index 97dfcfb101e0d6503fb654863cb32e719284d27d..93511f2e8fccc5d77b84e297a244331ef0aabeb5 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
@@ -58,6 +58,7 @@ struct RichardsMechanicsProcessData
     bool const apply_mass_lumping;
 
     MeshLib::PropertyVector<double>* element_saturation = nullptr;
+    MeshLib::PropertyVector<double>* element_porosity = nullptr;
     MeshLib::PropertyVector<double>* element_stresses = nullptr;
     MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;