diff --git a/libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh b/libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh
index f584e5be03..8449d2f2b4 100644
--- a/libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh
+++ b/libsrc/pylith/fekernels/IsotropicLinearPoroelasticity.hh
@@ -1850,11 +1850,16 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
 
         // Rheology Auxiliaries
         const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
 
-        g0[0] -= biotCoefficient * trace_strain_t;
+        // Rheological Auxiliaries
+
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
     } // g0p_implicit
 
     // ----------------------------------------------------------------------
@@ -1897,15 +1902,20 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
 
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
 
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        // Rheological Auxiliaries
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
     } // g0p_source
 
     // ----------------------------------------------------------------------
@@ -1948,15 +1958,19 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
 
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
 
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
     } // g0p_source_body
 
     // ----------------------------------------------------------------------
@@ -1999,15 +2013,19 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
 
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
 
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
     } // g0p_source_grav
 
     // ----------------------------------------------------------------------
@@ -2050,15 +2068,19 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
 
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
 
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
     } // g0p_source_grav_body
 
     // -----------------------------------------------------------------------------
@@ -3891,11 +3913,16 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
 
         // Rheology Auxiliaries
         const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
+
+        // Rheological Auxiliaries
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
 
-        g0[0] -= biotCoefficient * trace_strain_t;
     } // g0p_implicit
 
     // ----------------------------------------------------------------------
@@ -3938,13 +3965,18 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
+
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
     } // g0p_source
 
     // ----------------------------------------------------------------------
@@ -3987,13 +4019,19 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
+
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
+
     } // g0p_source_body
 
     // ----------------------------------------------------------------------
@@ -4036,13 +4074,18 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
+
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
     } // g0p_source_grav
 
     // ----------------------------------------------------------------------
@@ -4085,13 +4128,19 @@ public:
 
         // Solution Variables
         const PylithScalar trace_strain_t = poroelasticContext.trace_strain_t;
+        const PylithScalar pressure_t = poroelasticContext.pressure_t;
+
+        // Rheology Auxiliaries
+        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
+        const PylithScalar biotModulus = rheologyContext.biotModulus;
+
         // Poroelastic Auxiliaries
         const PylithScalar source = poroelasticContext.sourceDensity;
-        // Rheologic Auxiliaries
-        const PylithScalar biotCoefficient = rheologyContext.biotCoefficient;
 
-        g0[0] += source;
-        g0[0] -= biotCoefficient * trace_strain_t;
+        g0[0] -= source;
+        g0[0] += pressure_t / biotModulus;
+        g0[0] += biotCoefficient * trace_strain_t;
+
     } // g0p_source_grav_body
 
     // -----------------------------------------------------------------------------
diff --git a/libsrc/pylith/materials/Poroelasticity.cc b/libsrc/pylith/materials/Poroelasticity.cc
index d704589636..836164b502 100644
--- a/libsrc/pylith/materials/Poroelasticity.cc
+++ b/libsrc/pylith/materials/Poroelasticity.cc
@@ -457,20 +457,20 @@ pylith::materials::Poroelasticity::_setKernelsResidual(pylith::feassemble::Integ
         // Displacement
         const PetscPointFunc f0u = pylith::fekernels::Poroelasticity::f0u_explicit;
         const PetscPointFunc f1u = NULL;
-        const PetscPointFunc g0u = NULL; // pylith::fekernels::Poroelasticity::g0u;
+        const PetscPointFunc g0u = NULL;
         const PetscPointFunc g1u = NULL;
 
         // Pressure
-        const PetscPointFunc f0p = _rheology->getKernelf0p_explicit(coordsys);
-        const PetscPointFunc f1p = NULL;
-        const PetscPointFunc g0p = _rheology->getKernelg0p(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
-        const PetscPointFunc g1p = _rheology->getKernelg1p_explicit(coordsys, _gravityField); // Darcy velocity
+        const PetscPointFunc f0p = _rheology->getKernelg0p(coordsys, _useBodyForce, _gravityField, _useSourceDensity);
+        const PetscPointFunc f1p = _rheology->getKernelg1p_explicit(coordsys, _gravityField); // Darcy velocity;
+        const PetscPointFunc g0p = NULL;
+        const PetscPointFunc g1p = NULL;
 
         // Velocity
         const PetscPointFunc f0v = pylith::fekernels::Poroelasticity::f0v_explicit;
-        const PetscPointFunc f1v = NULL;
-        const PetscPointFunc g0v = r0;
-        const PetscPointFunc g1v = _rheology->getKernelg1v_explicit(coordsys);
+        const PetscPointFunc f1v = _rheology->getKernelg1v_explicit(coordsys);
+        const PetscPointFunc g0v = r0; // NEED TO MOVE TO LHS
+        const PetscPointFunc g1v = NULL;
 
         kernels.resize(6);
         kernels[0] = ResidualKernels("displacement", pylith::feassemble::Integrator::LHS, f0u, f1u);