diff --git a/src/rheology/GeoParams.jl b/src/rheology/GeoParams.jl
index ca65c7401..e8cfad774 100644
--- a/src/rheology/GeoParams.jl
+++ b/src/rheology/GeoParams.jl
@@ -7,11 +7,11 @@ function get_bulk_modulus(args::Vararg{Any, N}) where {N}
 end
 
 function get_shear_modulus(args::Vararg{Any, N}) where {N}
-    Kb = GeoParams.get_G(args...)
-    if isnan(Kb) || iszero(Kb)
+    G = GeoParams.get_G(args...)
+    if isnan(G) || iszero(G)
         return Inf
     end
-    return Kb
+    return G
 end
 
 get_thermal_expansion(args::Vararg{Any, N}) where {N} = get_α(args...)
diff --git a/src/stokes/PressureKernels.jl b/src/stokes/PressureKernels.jl
index fb1556f15..23410927a 100644
--- a/src/stokes/PressureKernels.jl
+++ b/src/stokes/PressureKernels.jl
@@ -7,9 +7,9 @@
 end
 
 ## Compressible
-@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, dt, r, θ_dτ)
+@parallel_indices (I...) function compute_P!(P, P0, RP, ∇V, Q, η, K, G, dt, r, θ_dτ)
     @inbounds RP[I...], P[I...] = _compute_P!(
-        P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], dt, r, θ_dτ
+        P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K[I...], G[I...], dt, r, θ_dτ
     )
     return nothing
 end
@@ -20,7 +20,8 @@ end
         P, P0, RP, ∇V, Q, η, rheology::NTuple{N, MaterialParams}, phase, dt, r, θ_dτ, args
     ) where {N}
     K = get_bulk_modulus(rheology, phase[I...])
-    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ)
+    G = get_shear_modulus(rheology, phase[I...])
+    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ)
     return nothing
 end
 
@@ -95,7 +96,8 @@ end
         ::Nothing,
     ) where {N, C <: JustRelax.CellArray}
     K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...]))
-    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ)
+    G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...]))
+    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ)
     return nothing
 end
 
@@ -115,7 +117,8 @@ end
         melt_fraction,
     ) where {N, C <: JustRelax.CellArray}
     K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...]))
-    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, dt, r, θ_dτ)
+    G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...]))
+    @inbounds RP[I...], P[I...] = _compute_P!(P[I...], P0[I...], ∇V[I...], Q[I...], η[I...], K, G, dt, r, θ_dτ)
     return nothing
 end
 
@@ -136,9 +139,10 @@ end
     ) where {N, C <: JustRelax.CellArray}
     @inbounds phase_ratio_I = phase_ratio[I...]
     @inbounds K = fn_ratio(get_bulk_modulus, rheology, phase_ratio_I)
+    @inbounds G = fn_ratio(get_shear_modulus, rheology, phase_ratio_I)
     @inbounds α = fn_ratio(get_thermal_expansion, rheology, phase_ratio_I)
     @inbounds RP[I...], P[I...] = _compute_P!(
-        P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ
+        P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, G, dt, r, θ_dτ
     )
     return nothing
 end
@@ -159,9 +163,10 @@ end
         melt_fraction,
     ) where {N, C <: JustRelax.CellArray}
     @inbounds K = fn_ratio(get_bulk_modulus, rheology, @cell(phase_ratio[I...]))
+    @inbounds G = fn_ratio(get_shear_modulus, rheology, @cell(phase_ratio[I...]))
     @inbounds α = fn_ratio(get_thermal_expansion, rheology, @cell(phase_ratio[I...]), (; ϕ = melt_fraction[I...]))
     @inbounds RP[I...], P[I...] = _compute_P!(
-        P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, dt, r, θ_dτ
+        P[I...], P0[I...], ∇V[I...], Q[I...], ΔTc[I...], α, η[I...], K, G, dt, r, θ_dτ
     )
     return nothing
 end
@@ -174,21 +179,23 @@ function _compute_P!(P, ∇V, Q, η, dt, r, θ_dτ)
     return RP, P
 end
 
-function _compute_P!(P, P0, ∇V, Q, η, K, dt, r, θ_dτ)
+function _compute_P!(P, P0, ∇V, Q, η, K, G, dt, r, θ_dτ)
     _Kdt = inv(K * dt)
+    _Gdt = inv(G * dt)
     _dt = inv(dt)
     RP = muladd(-(P - P0), _Kdt, (-∇V + (Q * _dt)))
-    ψ = inv(inv(r / θ_dτ * η) + _Kdt)
+    ψ = inv(inv(η) + _Gdt) * r / θ_dτ
     P = ((muladd(P0, _Kdt, (-∇V + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ)
 
     return RP, P
 end
 
-function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, dt, r, θ_dτ)
+function _compute_P!(P, P0, ∇V, Q, ΔTc, α, η, K, G, dt, r, θ_dτ)
     _Kdt = inv(K * dt)
+    _Gdt = inv(G * dt)
     _dt = inv(dt)
     RP = muladd(-(P - P0), _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt)))
-    ψ = inv(inv(r / θ_dτ * η) + _Kdt)
+    ψ = inv(inv(η) + _Gdt) * r / θ_dτ
     P = ((muladd(P0, _Kdt, (-∇V + (α * (ΔTc * _dt)) + (Q * _dt)))) * ψ + P) / (1 + _Kdt * ψ)
 
     return RP, P
diff --git a/src/stokes/Stokes2D.jl b/src/stokes/Stokes2D.jl
index 3af252bd8..81213852e 100644
--- a/src/stokes/Stokes2D.jl
+++ b/src/stokes/Stokes2D.jl
@@ -203,7 +203,7 @@ function _solve!(
             @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), _di)
 
             @parallel compute_P!(
-                stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, dt, r, θ_dτ
+                stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, ητ, K, G, dt, r, θ_dτ
             )
 
             @parallel (@idx ni .+ 1) compute_strain_rate!(
@@ -328,6 +328,7 @@ function _solve!(
     # end
 
     Kb = get_Kb(rheology)
+    G = get_G(rheology)
 
     # errors
     err_it1 = 1.0
@@ -368,7 +369,7 @@ function _solve!(
 
             @parallel (@idx ni) compute_∇V!(stokes.∇V, @velocity(stokes), _di)
             @parallel compute_P!(
-                stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, η, Kb, dt, r, θ_dτ
+                stokes.P, stokes.P0, stokes.R.RP, stokes.∇V, stokes.Q, η, Kb, G, dt, r, θ_dτ
             )
 
             update_ρg!(ρg[2], rheology, args)
diff --git a/src/stokes/Stokes3D.jl b/src/stokes/Stokes3D.jl
index 661e0301b..a3fdee728 100644
--- a/src/stokes/Stokes3D.jl
+++ b/src/stokes/Stokes3D.jl
@@ -81,6 +81,7 @@ function _solve!(
                 stokes.Q,
                 η,
                 K,
+                G,
                 dt,
                 pt_stokes.r,
                 pt_stokes.θ_dτ,