[ADDITION] Added a variant of the Herschel Bulkley rheology

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 ------------------------------------------------
 """

src/CreepLaw/HerschelBulkley.jl

@@ -0,0 +1,38 @@
+
+
+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 = 1e24Pa*s,
+        τ0 = 100e6Pa,
+        ηr = 1e20Pa*s,
+        Q  = 0K,
+        Tr = 1600K,
+        )
+        # 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

src/GeoParams.jl

@@ -174,6 +174,7 @@ module GeoParams
         GiordanoMeltViscosity,
         PowerlawViscous,
         ArrheniusType,
+        HerschelBulkley,
         CustomRheology,
         DislocationCreep,
         SetDislocationCreep,

src/Viscosity/Viscosity.jl

@@ -18,6 +18,25 @@ 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_viscosity_εII(v, εII; args...)
+
+@inline function compute_viscosity_εII(v::HerschelBulkley, εII; T = 0.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(B * (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
+    η = @fastpow (1.0 - exp(-2.0*η0*εII/τ0)) * (0.5*τ0/εII  + ηT*(εII./εr)^(1/n - 1))
+    return  η
+end
+
+
 
 # compute effective "creep" viscosity from deviatoric stress tensor
 """