diff --git a/src/ConstitutiveRelationships.jl b/src/ConstitutiveRelationships.jl
index c984f0663..bd849e359 100644
--- a/src/ConstitutiveRelationships.jl
+++ b/src/ConstitutiveRelationships.jl
@@ -70,7 +70,7 @@ export param_info,
 for myType in (
         :LinearViscous, :LinearMeltViscosity, :ViscosityPartialMelt_Costa_etal_2009, :GiordanoMeltViscosity,
         :DiffusionCreep, :DislocationCreep, :ConstantElasticity, :DruckerPrager, :ArrheniusType,
-        :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous,
+        :GrainBoundarySliding, :PeierlsCreep, :NonLinearPeierlsCreep, :PowerlawViscous, :HerschelBulkley,
     )
     @eval begin
         @inline compute_εII(a::$(myType), TauII, args) = compute_εII(a, TauII; args...)
diff --git a/src/CreepLaw/CreepLaw.jl b/src/CreepLaw/CreepLaw.jl
index 644aada13..850d96857 100644
--- a/src/CreepLaw/CreepLaw.jl
+++ b/src/CreepLaw/CreepLaw.jl
@@ -16,7 +16,7 @@
 abstract type AbstractCreepLaw{T} <: AbstractConstitutiveLaw{T} end
 
 export isvolumetric,
-    LinearViscous, PowerlawViscous, CorrectionFactor, AbstractCreepLaw, ArrheniusType
+    LinearViscous, PowerlawViscous, CorrectionFactor, AbstractCreepLaw, ArrheniusType, HerschelBulkley
 
 isvolumetric(a::AbstractCreepLaw) = false
 
@@ -75,6 +75,7 @@ include("PeierlsCreep.jl")
 include("NonLinearPeierlsCreep.jl")
 include("CustomRheology.jl")
 include("MeltViscosity.jl")
+include("HerschelBulkley.jl")
 
 # Linear viscous rheology ------------------------------------------------
 """
diff --git a/src/CreepLaw/HerschelBulkley.jl b/src/CreepLaw/HerschelBulkley.jl
new file mode 100644
index 000000000..501604611
--- /dev/null
+++ b/src/CreepLaw/HerschelBulkley.jl
@@ -0,0 +1,209 @@
+export HerschelBulkley,
+    compute_εII!,
+    compute_εII,
+    compute_τII!,
+    compute_τII,
+    compute_hb_viscosity_εII,
+    compute_hb_viscosity_τII
+
+struct HerschelBulkley{T, U1, U2, U3} <: AbstractCreepLaw{T}
+    n::T  # shear thinning exponent
+    η0::GeoUnit{T, U1} # "rigid" viscosity
+    τ0::GeoUnit{T, U2} # critical stress
+    ηr::GeoUnit{T, U1} # reference viscosity at the critical strain rate, which is given by 0.5*τ0/η0 and the critical temperature
+    Q::GeoUnit{T, U3} # temperature dependence of ηr, activation energy divided by R, unit is K
+    Tr::GeoUnit{T, U3} # reference temperature
+    function HerschelBulkley(;
+            n = 3.0,
+            η0 = 1.0e24Pa * s,
+            τ0 = 100.0e6Pa,
+            ηr = 1.0e20Pa * s,
+            Q = 0.0K,
+            Tr = 1273K,
+        )
+        # Convert to GeoUnits
+        η0U = convert(GeoUnit, η0)
+        τ0U = convert(GeoUnit, τ0)
+        ηrU = convert(GeoUnit, ηr)
+        QU = convert(GeoUnit, Q)
+        Tr = convert(GeoUnit, Tr)
+        # Extract struct types
+        T = typeof(η0U).types[1]
+        U1 = typeof(η0U).types[2]
+        U2 = typeof(τ0U).types[2]
+        U3 = typeof(Tr).types[2]
+        # Create struct
+        return new{T, U1, U2, U3}(
+            n, η0U, τ0U, ηrU, QU, Tr
+        )
+    end
+
+    function HerschelBulkley(n, η0, τ0, ηr, Q, Tr)
+        return HerschelBulkley(; n = n, η0 = η0, τ0 = τ0, ηr = ηr, Q = Q, Tr = Tr)
+    end
+end
+
+function compute_εII(a::HerschelBulkley, TauII; T = one(precision(a)), kwargs...)
+    η = compute_hb_viscosity_τII(a, TauII; T = T)
+    EpsII = 0.5 * TauII / η
+    return EpsII
+end
+
+function compute_εII(a::HerschelBulkley, TauII::Quantity; T = 1K, kwargs...)
+    η = compute_hb_viscosity_τII(a, TauII; T = T)
+    EpsII = 0.5 * TauII / η
+    return EpsII
+end
+
+"""
+    compute_εII!(EpsII::AbstractArray{T, N}, a::HerschelBulkley, TauII::AbstractArray{T, N}; T = ones(size(TauII)), kwargs...)
+
+In-place function for the second invariant of the strain rate for Herschel-Bulkley rheology.
+"""
+function compute_εII!(
+        EpsII::AbstractArray{_T, N},
+        a::HerschelBulkley,
+        TauII::AbstractArray{_T, N};
+        T = ones(size(TauII))::AbstractArray{_T, N},
+        kwargs...,
+    ) where {_T, N}
+    @inbounds for i in eachindex(EpsII)
+        EpsII[i] = compute_εII(a, TauII[i]; T = T[i])
+    end
+
+    return nothing
+end
+
+"""
+    compute_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...)
+
+"""
+function compute_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...)
+    η = compute_hb_viscosity_εII(a, EpsII; T = T)
+    TauII = 2 * EpsII * η
+    return TauII
+end
+
+function compute_τII(a::HerschelBulkley, EpsII::Quantity; T = 1K, kwargs...)
+    η = compute_hb_viscosity_εII(a, EpsII; T = T)
+    TauII = 2 * EpsII * η
+    return TauII
+end
+
+"""
+    compute_τII!(TauII::AbstractArray{T, N}, a::HerschelBulkley, EpsII::AbstractArray{T, N}; T = ones(size(EpsII)), kwargs...)
+
+In-place function for the second invariant of the stress for Herschel-Bulkley rheology.
+"""
+function compute_τII!(
+        TauII::AbstractArray{_T, N},
+        a::HerschelBulkley,
+        EpsII::AbstractArray{_T, N};
+        T = ones(size(EpsII)),
+        kwargs...,
+    ) where {_T, N}
+    @inbounds for i in eachindex(TauII)
+        TauII[i] = compute_τII(a, EpsII[i]; T = T[i])
+    end
+
+    return nothing
+end
+
+
+"""
+    compute_hb_viscosity_εII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...)
+
+function to compute the viscosity if EpsII is given
+"""
+@inline function compute_hb_viscosity_εII(v::HerschelBulkley, εII; T = 1.0, kwargs...)
+    η0, τ0, ηr, Q, Tr = if εII isa Quantity
+        @unpack_units η0, τ0, ηr, Q, Tr = v
+        η0, τ0, ηr, Q, Tr
+    else
+        @unpack_val η0, τ0, ηr, Q, Tr = v
+        η0, τ0, ηr, Q, Tr
+    end
+    (; n) = v
+
+    ηT = ηr * exp(Q * (1 / T - 1 / Tr)) # temperature dependence
+    εr = 0.5 * τ0 / η0 # strain rate at which the Bingham yield stress is reached, this is defined as the reference strain rate
+    η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1))
+    return η
+end
+
+
+"""
+compute_hb_viscosity_τII(a::HerschelBulkley, EpsII; T = one(precision(a)), kwargs...)
+
+function to compute the viscosity if TauII is given
+"""
+
+@inline function compute_hb_viscosity_τII(v::HerschelBulkley, τII; T = one(precision(v)), kwargs...)
+
+    η0, τ0, ηr, Q, Tr = if τII isa Quantity
+        @unpack_units η0, τ0, ηr, Q, Tr = v
+        η0, τ0, ηr, Q, Tr
+    else
+        @unpack_val η0, τ0, ηr, Q, Tr = v
+        T = ustrip(T)
+        η0, τ0, ηr, Q, Tr
+    end
+    (; n) = v
+
+    ηT = ηr * exp(Q * (1 / T - 1 / Tr))
+    εr = 0.5 * τ0 / η0
+
+    # initial guess
+    η = if τII < τ0
+        η0
+    elseif τII == τ0
+        (1 - exp(-one(η0))) * (η0 + ηT)
+    else
+        @pow (ηT / (1 - τ0 / τII))^n * (τII / (2 * εr))^(1 - n)
+    end
+
+    εII = 0.5 * τII / η
+    εII_unit = τII isa Quantity ? unit(εII) : one(εII)
+
+    # strip ALL quantities to plain floats before Newton iteration
+    # so that ForwardDiff never sees Quantity{Dual} types
+    τII_s = ustrip(τII)
+    η0_s = ustrip(η0)
+    τ0_s = ustrip(τ0)
+    ηT_s = ustrip(ηT)
+    εr_s = ustrip(εr)
+    εII_s = ustrip(εII)
+
+    tol = 1.0e-10
+    it_max = 100
+
+    for _ in 1:it_max
+        f, dfdε = value_and_partial(
+            ε -> fres_hb(ε, τII_s, η0_s, τ0_s, ηT_s, n, εr_s),  # all plain floats
+            εII_s
+        )
+        Δε = f / dfdε
+        εII_s -= Δε
+        abs(Δε) / (abs(εII_s) + eps(typeof(εII_s))) < tol && break
+    end
+
+    εII = εII_s * εII_unit
+    η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1))
+
+    return η
+end
+
+# non-dimensional residual
+function fres_hb(εII::Number, τII::Number, η0::Number, τ0::Number, ηT::Number, n::Number, εr::Number)
+    η = @pow (1.0 - exp(-2.0 * η0 * εII / τ0)) * (0.5 * τ0 / εII + ηT * (εII / εr)^(one(n) / n - 1))
+    return 2.0 * η * εII - τII
+end
+
+
+# print info
+function show(io::IO, g::HerschelBulkley)
+    return print(
+        io,
+        "Hershel Bulkley viscosity: η0=$(Value(g.η0)), τ0=$(Value(g.τ0)), ηr=$(Value(g.ηr)), n=$(g.n), Q=$(Value(g.Q)), Tr=$(Value(g.Tr))",
+    )
+end
diff --git a/src/GeoParams.jl b/src/GeoParams.jl
index 598db8002..e46d8c0be 100644
--- a/src/GeoParams.jl
+++ b/src/GeoParams.jl
@@ -174,6 +174,7 @@ module GeoParams
         GiordanoMeltViscosity,
         PowerlawViscous,
         ArrheniusType,
+        HerschelBulkley,
         CustomRheology,
         DislocationCreep,
         SetDislocationCreep,
diff --git a/src/Viscosity/Viscosity.jl b/src/Viscosity/Viscosity.jl
index f7aa12209..5fb72067c 100644
--- a/src/Viscosity/Viscosity.jl
+++ b/src/Viscosity/Viscosity.jl
@@ -18,6 +18,9 @@ end
 
 @inline compute_viscosity_εII(v::LinearViscous, εII, args) = v.η.val
 @inline compute_viscosity_εII(v::ConstantElasticity, εII, args) = v.G * args.dt
+@inline compute_viscosity_εII(v::HerschelBulkley, εII, args) = compute_hb_viscosity_εII(v, εII; args...)
+
+
 
 # compute effective "creep" viscosity from deviatoric stress tensor
 """
@@ -33,6 +36,7 @@ end
 
 @inline compute_viscosity_τII(v::LinearViscous, τII, args) = v.η.val
 @inline compute_viscosity_τII(v::ConstantElasticity, τII, args) = v.G * args.dt
+@inline compute_viscosity_τII(v::HerschelBulkley, εII, args) = compute_hb_viscosity_τII(v, εII; args...)
 
 for fn in (:compute_viscosity_εII, :compute_viscosity_τII)
     @eval begin
diff --git a/test/test_CreepLaw.jl b/test/test_CreepLaw.jl
index 1ed916467..d0f64831d 100644
--- a/test/test_CreepLaw.jl
+++ b/test/test_CreepLaw.jl
@@ -22,6 +22,9 @@ using GeoParams
 
     x = GiordanoMeltViscosity()
     @test isbits(x)
+
+    x = HerschelBulkley()
+    @test isbits(x)
     # This tests the MaterialParameters structure
     CharUnits_GEO = GEO_units(; viscosity = 1.0e19, length = 1000km)
 
@@ -343,6 +346,147 @@ using GeoParams
     end
     T_plots = 10000 ./ ustrip.(T)
 
+
+    # Herschel-Bulkley viscosity (Regularized Bingham rheology)
+    x1_D = HerschelBulkley()
+    @test isDimensional(x1_D) == true
+    @test x1_D.n == 3
+    @test x1_D.η0.val == 1.0e24
+    @test x1_D.τ0.val == 100.0e6
+    @test x1_D.ηr.val == 1.0e20
+    @test x1_D.Q.val == 0
+    @test x1_D.Tr.val == 1273
+
+    @test x1_D.η0.unit == Pas
+    @test x1_D.τ0.unit == Pa
+    @test x1_D.ηr.unit == Pas
+    @test x1_D.Q.unit == K
+    @test x1_D.Tr.unit == K
+
+    x1_ND = nondimensionalize(x1_D, CharUnits_GEO)
+    @test isDimensional(x1_ND) == false
+    @test x1_ND.η0 * 1.0 == 1.0e5
+    @test x1_ND.τ0 * 1.0 ≈ 10.0
+    @test x1_ND.ηr * 1.0 == 1.0e1
+    @test x1_ND.Q * 1.0 == 0.0
+    @test x1_ND.Tr * 1.0 ≈ 1.0 rtol = 1.0e-3
+    @test dimensionalize(x1_ND, CharUnits_GEO) == x1_D
+
+    # below yield case - dimensional
+    εt = 0.5e-20 * s^-1
+    τ1273 = compute_τII(x1_D, εt, (; T = 1273.0K))
+
+    @test ustrip(τ1273) ≈ 2 * x1_D.η0.val * ustrip(εt) rtol = 1.0e-4
+    ε1273 = compute_εII(x1_D, τ1273, (; T = 1273.0K))
+    @test εt ≈ ε1273 rtol = 1.0e-6 # see if we can recover the strain rate
+
+    # below yield case - nondimensional
+    εr = 0.5 * x1_ND.τ0 / x1_ND.η0 # yield strain rate
+    εt = 1.0e-4 * εr # well below the yield
+
+    τ1 = compute_τII(x1_ND, εt, (; T = 1.0))
+    @test τ1 ≈ 2 * x1_ND.η0.val * ustrip(εt) rtol = 1.0e-4
+    ε1 = compute_εII(x1_ND, τ1, (; T = 1.0))
+    @test εt ≈ ε1 rtol = 1.0e-6 # see if we can revocer the strain rate
+
+    εt = [1.0e-4 * εr; 1.0e-4 * εr]
+    τ1 = [1.0e3, 1.0e3] # initialize τ1
+    compute_τII!(τ1, x1_ND, εt, (; T = ones(size(εt))))
+    @test τ1 ≈ [2 * x1_ND.η0.val * ustrip(εt[1]), 2 * x1_ND.η0.val * ustrip(εt[2])] rtol = 1.0e-4
+
+    # above yield case, dimensional
+    x1_D = HerschelBulkley(
+        n = 3.0,
+        η0 = 1.0e24Pa * s,
+        τ0 = 100.0e6Pa,
+        ηr = 1.0e20Pa * s,
+        Q = 0.0K,
+        Tr = 1273K,
+    )
+
+    εt = 0.5e-9 * s^-1
+    εr = 0.5 * x1_D.τ0.val / x1_D.η0.val # yield strain rate
+
+    η_approx = 0.5 * x1_D.τ0.val / ustrip(εt) + x1_D.ηr.val * (ustrip(εt) / εr)^(one(x1_D.n) / x1_D.n - 1)
+    τ1273 = compute_τII(x1_D, εt; T = 1273.0K)
+    ε1273 = compute_εII(x1_D, τ1273; T = 1273.0K)
+
+    @test ustrip(τ1273) ≈ 2 * η_approx * ustrip(εt) rtol = 1.0e-6
+    @test ε1273 ≈ εt rtol = 1.0e-6
+    #
+    #
+    x1_ND = nondimensionalize(x1_D, CharUnits_GEO)
+    εt_ND = nondimensionalize(εt, CharUnits_GEO)
+    #
+    τ1_ND = compute_τII(x1_ND, εt_ND; T = 1.0)
+    ε1_ND = compute_εII(x1_ND, τ1_ND; T = 1.0)
+
+    @test ε1_ND ≈ εt_ND rtol = 1.0e-6
+
+
+    # temperature dependence test
+    x1_Q = HerschelBulkley(
+        n = 3.0,
+        η0 = 1.0e24Pa * s,
+        τ0 = 100.0e6Pa,
+        ηr = 1.0e20Pa * s,
+        Q = 63772.0K,   # E/R for olivine, E ≈ 530 kJ/mol
+        Tr = 1273K,
+    )
+
+    # compute ηT at both temperatures analytically
+    Q_val = 63772.0
+    Tr_val = 1273.0
+    T1_val = 1273.0   # reference temperature → ηT = ηr
+    T2_val = 1473.0   # higher temperature   → ηT < ηr
+
+    ηT1 = x1_Q.ηr.val * exp(Q_val * (1 / T1_val - 1 / Tr_val))  # = 1e20
+    ηT2 = x1_Q.ηr.val * exp(Q_val * (1 / T2_val - 1 / Tr_val))  # << 1e20
+
+    εr_Q = 0.5 * x1_Q.τ0.val / x1_Q.η0.val  # reference strain rate
+
+    # --- above yield, dimensional ---
+    εt_Q = 0.5e-9 * s^-1  # well above yield
+
+    # at T = 1273K (reference): ηT = ηr, same as Q=0 case
+    τ_T1 = compute_τII(x1_Q, εt_Q; T = T1_val * K)
+    ε_T1 = compute_εII(x1_Q, τ_T1; T = T1_val * K)
+    η_approx_T1 = 0.5 * x1_Q.τ0.val / ustrip(εt_Q) + ηT1 * (ustrip(εt_Q) / εr_Q)^(one(x1_Q.n) / x1_Q.n - 1)
+    @test ustrip(τ_T1) ≈ 2 * η_approx_T1 * ustrip(εt_Q) rtol = 1.0e-6
+    @test ε_T1 ≈ εt_Q rtol = 1.0e-6
+
+    # at T = 1473K: ηT < ηr → lower viscosity → lower stress at same strain rate
+    τ_T2 = compute_τII(x1_Q, εt_Q; T = T2_val * K)
+    ε_T2 = compute_εII(x1_Q, τ_T2; T = T2_val * K)
+    η_approx_T2 = 0.5 * x1_Q.τ0.val / ustrip(εt_Q) + ηT2 * (ustrip(εt_Q) / εr_Q)^(one(x1_Q.n) / x1_Q.n - 1)
+    @test ustrip(τ_T2) ≈ 2 * η_approx_T2 * ustrip(εt_Q) rtol = 1.0e-6
+    @test ε_T2 ≈ εt_Q rtol = 1.0e-6
+
+    # monotonicity: higher T → lower τII at same εII (thermal softening)
+    @test ustrip(τ_T2) < ustrip(τ_T1)
+
+    # monotonicity: higher T → higher εII at same τII (thermal softening)
+    τ_fixed = compute_τII(x1_Q, εt_Q; T = T1_val * K)
+    ε_at_T1 = compute_εII(x1_Q, τ_fixed; T = T1_val * K)
+    ε_at_T2 = compute_εII(x1_Q, τ_fixed; T = T2_val * K)
+    @test ustrip(ε_at_T2) > ustrip(ε_at_T1)
+
+    # nondimensional temperature dependence
+    x1_Q_ND = nondimensionalize(x1_Q, CharUnits_GEO)
+    εt_ND = nondimensionalize(εt_Q, CharUnits_GEO)
+    T1_ND = nondimensionalize(T1_val * K, CharUnits_GEO)
+    T2_ND = nondimensionalize(T2_val * K, CharUnits_GEO)
+
+    τ_T1_ND = compute_τII(x1_Q_ND, εt_ND; T = T1_ND)
+    τ_T2_ND = compute_τII(x1_Q_ND, εt_ND; T = T2_ND)
+
+    # ND results should be consistent with dimensional results after rescaling
+    @test ustrip(nondimensionalize(τ_T1, CharUnits_GEO)) ≈ τ_T1_ND rtol = 1.0e-6
+    @test ustrip(nondimensionalize(τ_T2, CharUnits_GEO)) ≈ τ_T2_ND rtol = 1.0e-6
+
+    # thermal softening should hold in ND as well
+    @test τ_T2_ND < τ_T1_ND
+
     # fig = Figure(size = (800, 600))
     # ax = Axis(fig[1, 1], xlabel = "T [C]", ylabel = "log10(η [Pas])")
     # lines!(ax, ustrip.(T) .- 273.15, log10.(ustrip.(η_MORB)), label = "MORB, εII=1e-5/s")