diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index f64ec0152854ab908941d047af478fa5c711ea2e..6cfa2b4c248bffcb34084ef6fda7fa1e730ec158 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -243,7 +243,7 @@ void RichardsMechanicsLocalAssembler<
              displacement_index, displacement_index)
             .noalias() += B.transpose() * C * B * w;
 
-        double const rho = rho_SR * (1 - porosity) + porosity * rho_LR;
+        double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
         rhs.template segment<displacement_size>(displacement_index).noalias() +=
             N_u_op.transpose() * rho * b * w;
 
@@ -452,7 +452,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
                 displacement_index, displacement_index)
             .noalias() += B.transpose() * C * B * w;
 
-        double const rho = rho_SR * (1 - porosity) + porosity * rho_LR;
+        double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
         local_rhs.template segment<displacement_size>(displacement_index)
             .noalias() -=
             (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;