diff --git a/src/Calculators/Rate.jl b/src/Calculators/Rate.jl
index c873c58a0..65d165a30 100644
--- a/src/Calculators/Rate.jl
+++ b/src/Calculators/Rate.jl
@@ -7,35 +7,38 @@ export AbstractRate
 abstract type AbstractFalloffRate <: AbstractRate end
 export AbstractFalloffRate
 
-@with_kw struct Arrhenius{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct Arrhenius{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
-@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-arr.Ea/(R*T))
+@inline (arr::Arrhenius)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
+@inline (arr::Arrhenius)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
 export Arrhenius
 
-@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct StickingCoefficient{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
-@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-arr.Ea/(R*T)),1.0)
+@inline (arr::StickingCoefficient)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages),1.0)
+@inline (arr::StickingCoefficient)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath min(arr.A*T^arr.n*exp(-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages),1.0)
 export StickingCoefficient
 
-@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,B} <: AbstractRate
+@with_kw struct Arrheniusq{N<:Real,K<:Real,Q<:Real,P<:AbstractRateUncertainty,B,H<:AbstractRateCoverageDependence} <: AbstractRate
         A::N
         n::K
         Ea::Q
         q::B = 0.0
         unc::P = EmptyRateUncertainty()
+        covdep::H = EmptyRateCoverageDependence()
 end
-@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
-@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-arr.Ea-arr.q*F*phi)/(R*T))
+@inline (arr::Arrheniusq)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))-arr.q*F*phi)/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
+@inline (arr::Arrheniusq)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath arr.A*T^arr.n*exp((-(arr.Ea+getcovdepactivationbarriercorrection(arr.covdep,T,coverages))-arr.q*F*phi)/(R*T))*getcovdepfactorcorrection(arr.covdep,T,coverages)
 export Arrheniusq
 
 @with_kw struct PdepArrhenius{T<:Real,Q<:AbstractRateUncertainty,Z<:AbstractRate} <: AbstractRate
@@ -45,7 +48,7 @@ export Arrheniusq
 end
 PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} = PdepArrhenius(sort(Ps),arrs)
 
-@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrhenius)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,V<:Real,S<:Real}
     inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
 
     if inds[2] == -1
@@ -60,15 +63,16 @@ PdepArrhenius(Ps::Array{Q,1},arrs::Array{Z,1}) where {Q<:Real,Z<:AbstractRate} =
 end
 export PdepArrhenius
 
-@with_kw struct MultiArrhenius{Q<:AbstractRateUncertainty} <: AbstractRate
+@with_kw struct MultiArrhenius{Q<:AbstractRateUncertainty,H<:AbstractRateCoverageDependence} <: AbstractRate
     arrs::Array{Arrhenius,1}
     unc::Q = EmptyRateUncertainty()
+    covdep::H = EmptyRateCoverageDependence()
 end
 
-@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (marr::MultiArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for arr in marr.arrs
-        @fastmath out += arr(T)
+        @fastmath out += arr(T;coverages=coverages)
     end
     return out
 end
@@ -79,7 +83,7 @@ export MultiArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (parr::MultiPdepArrhenius)(;T::Q,P::R=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     out = 0.0
     for pdar in parr.parrs
         @fastmath out += pdar(T=T,P=P)
@@ -95,7 +99,7 @@ export MultiPdepArrhenius
     unc::Q = EmptyRateUncertainty()
 end
 
-(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
+(tbarr::ThirdBody)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real} = C*(tbarr.arr(T))
 export ThirdBody
 
 @with_kw struct Lindemann{N<:Integer,K<:AbstractFloat,Q<:AbstractRateUncertainty} <: AbstractFalloffRate
@@ -106,7 +110,7 @@ export ThirdBody
     unc::Q = EmptyRateUncertainty()
 end
 
-@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (lnd::Lindemann)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = lnd.arrlow(T=T)
     kinf = lnd.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
@@ -126,7 +130,7 @@ export Lindemann
     unc::R = EmptyRateUncertainty()
 end
 
-@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troe)(;T::Q,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = tr.arrlow(T=T)
     kinf = tr.arrhigh(T=T)
     @fastmath Pr = k0*C/kinf
@@ -195,7 +199,7 @@ export getredtemp
 end
 export getredpress
 
-@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshev)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,coverages=nothing) where {N<:Real,B<:Real,Q<:Real}
     k = 0.0
     Tred = getredtemp(ch,T)
     Pred = getredpress(ch,P)
@@ -210,7 +214,12 @@ end
 export Chebyshev
 
 @inline function getkineticstype(kin::B) where {B<:AbstractRate}
-    return extracttypename(typeof(kin).name)
+    cdep = [string(x) for x in typeof(kin).parameters if occursin("CoverageDependence",string(x))]
+    if length(cdep) == 0 || cdep[1] == "EmptyRateCoverageDependence"
+        return extracttypename(typeof(kin).name)
+    else
+        return extracttypename(typeof(kin).name) * "_" * cdep[1]
+    end
 end
 
 @inline function getkineticstype(kin::PdepArrhenius)
diff --git a/src/Calculators/RateCoverageDependence.jl b/src/Calculators/RateCoverageDependence.jl
new file mode 100644
index 000000000..f88f446b9
--- /dev/null
+++ b/src/Calculators/RateCoverageDependence.jl
@@ -0,0 +1,72 @@
+using Parameters
+
+abstract type AbstractRateCoverageDependence end
+export AbstractRateCoverageDependence
+
+struct EmptyRateCoverageDependence <: AbstractRateCoverageDependence end
+export EmptyRateCoverageDependence
+
+getcovdepactivationbarriercorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 0.0
+getcovdepfactorcorrection(covdep::EmptyRateCoverageDependence,T,coverages) = 1.0
+
+"""
+Polynomial rate coefficient coverage dependence
+correction form is: 10^(sum(a_i*theta_k^i,i)) * theta_k^m * exp(sum(E_i*theta_k^i,i)/(RT))
+dictionaries map the name of the species (k) to polynomials with ordering a_1 + a_2 * theta + a_3 * theta^2 etc.
+Epolys corresponds to E_i (array of polynomial coefficients), avals corresponds to a_i (array of polynomial coefficients), 
+ms corresponds to m (float)
+"""
+@with_kw struct PolynomialRateCoverageDependence <: AbstractRateCoverageDependence
+    avals::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    Epolys::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    ms::Dict{String,Vector{Float64}} = Dict{String,Float64}()
+    indavals::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+    indEpolys::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+    indms::Dict{Int64,Vector{Float64}} = Dict{Int64,Float64}()
+end
+@inline function PolynomialRateCoverageDependence(Epolys,ms=Dict{String,Float64}(),avals=Dict{String,Vector{Float64}}())
+    return PolynomialRateCoverageDependence(;Epolys=Epolys,
+            ms=ms,avals=avals,indavals=Dict{Int64,Vector{Float64}}(),indEpolys=Dict{Int64,Vector{Float64}}(),
+            indms=Dict{Int64,Float64}())
+end
+@inline function getcovdepactivationbarriercorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    if coverages === nothing
+        return 0.0
+    elseif length(coverages) == 0
+        return zero(typeof(coverages).parameters[1])
+    end
+    E = 0.0
+    
+    for (ind,Epoly) in covdep.indEpolys
+        E += evalpoly(coverages[ind],Epoly)
+    end
+    return E
+end
+export getcovdepactivationbarriercorrection
+@inline function getcovdepfactorcorrection(covdep::PolynomialRateCoverageDependence,T,coverages)
+    if coverages === nothing
+        return 1.0
+    elseif length(coverages) == 0
+        return one(typeof(coverages).parameters[1])
+    end
+    av = 0.0
+    for (ind,a) in covdep.indavals
+       av += coverages[ind]*a
+    end
+
+    if av != 0
+        A = 10.0^av
+    else 
+        A = 1.0
+    end
+
+    for (ind,m) in covdep.indms 
+        A *= coverages[ind]^m
+    end
+
+    return A
+end
+export PolynomialRateCoverageDependence
+
+export getcovdepactivationbarriercorrection
+export getcovdepfactorcorrection
\ No newline at end of file
diff --git a/src/Calculators/Ratevec.jl b/src/Calculators/Ratevec.jl
index 7b5045743..ec65c8af9 100644
--- a/src/Calculators/Ratevec.jl
+++ b/src/Calculators/Ratevec.jl
@@ -10,7 +10,7 @@ export AbstractRatevec
         n::K
         Ea::Q
 end
-function Arrheniusvec(arrs::T) where {T<:AbstractArray}
+function Arrheniusvec(arrs)
     A = zeros(length(arrs)) 
     n = zeros(length(arrs)) 
     Ea = zeros(length(arrs))
@@ -21,8 +21,8 @@ function Arrheniusvec(arrs::T) where {T<:AbstractArray}
     end
     return Arrheniusvec(A=A,n=n,Ea=Ea)
 end
-@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
-@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
+@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
 export Arrheniusvec
 
 @with_kw struct Chebyshevvec{T<:AbstractArray,Q<:Real,S<:Real,V<:Real,B<:Real} <: AbstractRate
@@ -86,7 +86,7 @@ export getredtemp
 end
 export getredpress
 
-@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
+@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0,coverages=nothing) where {N<:Real,B<:Real,Q<:Real}
     k = zeros(N,size(ch.coefs)[1])
     Tred = getredtemp(ch,T)
     Pred = getredpress(ch,P)
@@ -170,7 +170,7 @@ function Troevec(troes::T) where {T<:AbstractArray}
 end
 export Troevec
     
-@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
+@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0,coverages=nothing) where {Q<:Real,R<:Real,S<:Real}
     k0 = tr.arrlow(T=T)
     kinf = tr.arrhigh(T=T)
     @fastmath Pr = k0.*C./kinf
@@ -201,7 +201,7 @@ function PdepArrheniusvec(pdeparrs::T) where {T<:AbstractArray}
 end
 export PdepArrheniusvec
 
-@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
+@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0,coverages=nothing) where {Q<:Real,V<:Real,S<:Real}
     inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
     if inds[2] == -1
         return @inbounds parr.arrvecs[inds[1]](T=T)
diff --git a/src/Calculators/Thermo.jl b/src/Calculators/Thermo.jl
index 534873a56..936c69446 100644
--- a/src/Calculators/Thermo.jl
+++ b/src/Calculators/Thermo.jl
@@ -3,18 +3,19 @@ using Parameters
 abstract type AbstractThermo end
 export AbstractThermo
 
-@inline function getGibbs(th::P,T::N) where {N<:Number,P<:AbstractThermo}
-    return @fastmath getEnthalpy(th,T)-T*getEntropy(th,T)
+@inline function getGibbs(th::P,T::N;coverages=nothing) where {N<:Number,P}
+    return @fastmath (getEnthalpy(th,T)+getcovdepenthalpycorrection(th.covdep,T,coverages))-T*(getEntropy(th,T)+getcovdepentropycorrection(th.covdep,T,coverages))
 end
 
-@with_kw struct NASApolynomial <: AbstractThermo
+@with_kw struct NASApolynomial{H<:AbstractThermoCoverageDependence} <: AbstractThermo
     coefs::Array{Float64,1}
     Tmin::Float64
     Tmax::Float64
+    covdep::H = EmptyThermoCoverageDependence()
 end
 export NASApolynomial
 
-@inline function calcHSCpdless(poly::NASApolynomial,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomial,T::X;coverages=nothing) where {X<:Real}
     if length(poly.coefs) != 7
         Tpoly0 = T
         Tpoly1 = T*T
@@ -58,43 +59,51 @@ export NASApolynomial
         hdivRT =ct2+0.5*ct3+0.33333333333*ct4+0.25*ct5+0.2*ct6+poly.coefs[6]*Tpoly4
         sdivR = Tpoly6*ct2 + ct3 + 0.5*ct4 + 0.33333333333*ct5 + 0.25*ct6 + poly.coefs[7]
     end
+    hdivRT += getcovdepenthalpycorrection(th.covdep,T,coverages)/(R*T)
+    sdivR += getcovdepentropycorrection(th.covdep,T,coverages)/R
+    cpdivR += getcovdepheatcapacitycorrection(th.covdep,T,coverages)/R
     return (cpdivR,hdivRT,sdivR)::Tuple{Float64,Float64,Float64}
 end
 export calcHSCpdless
 
-@inline function getHeatCapacity(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getHeatCapacity(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:7])/T^2*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:7])/T^2*R + getcovdepheatcapacitycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:5])*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:5])*R + getcovdepheatcapacitycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@inline function getEntropy(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getEntropy(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath ((-poly.coefs[1]/(2*T)-poly.coefs[2])/T+poly.coefs[3]*log(T)+T*evalpoly(T,poly.coefs[4:end-2]./(1:4))+poly.coefs[end])*R
+        return @views @inbounds @fastmath ((-poly.coefs[1]/(2*T)-poly.coefs[2])/T+poly.coefs[3]*log(T)+T*evalpoly(T,poly.coefs[4:end-2]./(1:4))+poly.coefs[end])*R + getcovdepentropycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath (poly.coefs[1]*log(T)+T*evalpoly(T,poly.coefs[2:end-2]./(1:4))+poly.coefs[end])*R
+        return @views @inbounds @fastmath (poly.coefs[1]*log(T)+T*evalpoly(T,poly.coefs[2:end-2]./(1:4))+poly.coefs[end])*R + getcovdepentropycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@inline function getEnthalpy(poly::NASApolynomial,T::N) where {N<:Number}
+@inline function getEnthalpy(poly::NASApolynomial,T::N,coverages=nothing) where {N<:Number}
     if length(poly.coefs) == 9
-        return @views @inbounds @fastmath ((-poly.coefs[1]/T+poly.coefs[2]*log(T))/T+evalpoly(T,poly.coefs[3:end-2]./(1:5)))*R*T+poly.coefs[end-1]*R
+        return @views @inbounds @fastmath ((-poly.coefs[1]/T+poly.coefs[2]*log(T))/T+evalpoly(T,poly.coefs[3:end-2]./(1:5)))*R*T+poly.coefs[end-1]*R+th.covdep(T,coverages) + getcovdepenthalpycorrection(poly.covdep,T,coverages)
     elseif length(poly.coefs) == 7
-        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:end-2]./(1:5))*R*T+poly.coefs[end-1]*R
+        return @views @inbounds @fastmath evalpoly(T,poly.coefs[1:end-2]./(1:5))*R*T+poly.coefs[end-1]*R + getcovdepenthalpycorrection(poly.covdep,T,coverages)
     else
         throw(error("NASA polynomial has a number of coefficients not equal to 9 or 7"))
     end
 end
 
-@with_kw struct NASA{T<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct NASA{T<:AbstractThermoUncertainty,H<:AbstractThermoCoverageDependence} <: AbstractThermo
     polys::Array{NASApolynomial,1}
     unc::T = EmptyThermoUncertainty()
+    covdep::H = EmptyThermoCoverageDependence()
+end
+function NASA(polys,unc=EmptyThermoUncertainty(),covdep=EmptyThermoCoverageDependence())
+    ps = [polys.covdep isa EmptyThermoCoverageDependence ? NASApolynomial(p.coefs,p.Tmin,p.Tmax,covdep) : p for p in polys]
+    return NASA(ps,unc,covdep)
 end
 export NASA
 
@@ -111,13 +120,13 @@ export NASA
 end
 export selectPoly
 
-@inline getHeatCapacity(nasa::NASA,T::N) where {N<:Number} = getHeatCapacity(selectPoly(nasa,T),T)
+@inline getHeatCapacity(nasa::NASA,T::N) where {N<:Number} = getHeatCapacity(selectPoly(nasa,T),T) #we do not add covdep corrections in NASA because it is also added in the NASAPolynomial, so it would double count
 @inline getEntropy(nasa::NASA,T::N) where {N<:Number} = getEntropy(selectPoly(nasa,T),T)
 @inline getEnthalpy(nasa::NASA,T::N) where {N<:Number} = getEnthalpy(selectPoly(nasa,T),T)
 @inline getGibbs(nasa::NASA,T::N) where {N<:Number} = getGibbs(selectPoly(nasa,T),T)
 @inline calcHSCpdless(nasa::NASA,T::N) where {N<:Real}= calcHSCpdless(selectPoly(nasa,T),T)
 
-@with_kw struct Wilhoit{N,Q,T,P,U,R<:Number,M<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct Wilhoit{N,Q,T,P,U,R<:Number,M<:AbstractThermoUncertainty,H<:AbstractThermoCoverageDependence} <: AbstractThermo
     Cp0::N
     Cpinf::T
     coefs::Array{Q,1}
@@ -125,15 +134,16 @@ export selectPoly
     S0::U
     B::R
     unc::M = EmptyThermoUncertainty()
+    covdep::H = EmptyThermoCoverageDependence()
 end
 export Wilhoit
 
-@inline function getHeatCapacity(w::Wilhoit,T::N) where {N<:Number}
+@inline function getHeatCapacity(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
-    return @fastmath w.Cp0 + (w.Cpinf-w.Cp0)*y^2*(1+(y-1)*evalpoly(y,w.coefs))
+    return @fastmath w.Cp0 + (w.Cpinf-w.Cp0)*y^2*(1+(y-1)*evalpoly(y,w.coefs)) + getcovdepheatcapacitycorrection(th.covdep,T,coverages)
 end
 
-@inline function getEnthalpy(w::Wilhoit,T::N) where {N<:Number}
+@inline function getEnthalpy(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
     return @views @fastmath @inbounds w.H0 + w.Cp0 * T - (w.Cpinf - w.Cp0) * T * (
             y * y * ((3 * w.coefs[1] + sum(w.coefs[2:end])) / 6. +
@@ -141,12 +151,12 @@ end
                      (5 * w.coefs[3] + w.coefs[4]) * y^2 / 20. +
                      w.coefs[4] * y^3 / 5.) +
             (2 + sum(w.coefs)) * (y / 2. - 1 + (1.0 / y - 1.) * log(w.B + T))
-        )
+        ) + getcovdepenthalpycorrection(th.covdep,T,coverages)
 end
 
-@inline function getEntropy(w::Wilhoit,T::N) where {N<:Number}
+@inline function getEntropy(w::Wilhoit,T::N,coverages=nothing) where {N<:Number}
     @fastmath y = T/(T+w.B)
-    return @fastmath w.S0 + w.Cpinf*log(T)-(w.Cpinf-w.Cp0)*(log(y)+y*(1+y*evalpoly(y,w.coefs./(2:5))))
+    return @fastmath w.S0 + w.Cpinf*log(T)-(w.Cpinf-w.Cp0)*(log(y)+y*(1+y*evalpoly(y,w.coefs./(2:5)))) + getcovdepentropycorrection(th.covdep,T,coverages)
 end
 
 @with_kw struct ConstantG{B<:Number,J} <: AbstractThermo
@@ -155,7 +165,7 @@ end
 end
 export ConstantG
 
-@inline function getGibbs(cg::ConstantG,T::B) where {B<:Number}
+@inline function getGibbs(cg::ConstantG,T::B,coverages=nothing) where {B<:Number}
     return cg.G
 end
 
diff --git a/src/Calculators/ThermoCoverageDependence.jl b/src/Calculators/ThermoCoverageDependence.jl
new file mode 100644
index 000000000..5d5c10c1c
--- /dev/null
+++ b/src/Calculators/ThermoCoverageDependence.jl
@@ -0,0 +1,37 @@
+using Parameters
+
+abstract type AbstractThermoCoverageDependence end
+export AbstractThermoCoverageDependence
+
+struct EmptyThermoCoverageDependence <: AbstractThermoCoverageDependence end
+@inline getcovdepenthalpycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+@inline getcovdepentropycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+@inline getcovdepheatcapacitycorrection(covdep::EmptyThermoCoverageDependence,T,coverages) = zero(T)
+export EmptyThermoCoverageDependence
+
+"""
+Polynomial energy coverage dependence
+correction is applied to enthalpy
+poly is a dictionary mapping species names to polynomial coverage dependence corrections to their energy/enthalpy
+polynomial ordering is a_1 + a_2 * theta + a_3 * theta^2 etc.
+"""
+@with_kw struct PolynomialThermoEnergyCoverageDependence <: AbstractThermoCoverageDependence
+    polys::Dict{String,Vector{Float64}} = Dict{String,Vector{Float64}}()
+    indpolys::Dict{Int64,Vector{Float64}} = Dict{Int64,Vector{Float64}}()
+end
+@inline function PolynomialThermoEnergyCoverageDependence(polys)
+    return PolynomialThermoEnergyCoverageDependence(polys,Dict{Int64,Vector{Float64}}())
+end
+@inline function getcovdepenthalpycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages)
+    if length(coverages) == 0
+        return zero(typeof(coverages).parameters[1])
+    end
+    E = zero(coverages[1])
+    for (ind,poly) in covdep.indpolys
+        E += evalpoly(coverages[ind],poly)
+    end
+    return E
+end
+@inline getcovdepentropycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages) = zero(coverages[1])
+@inline getcovdepheatcapacitycorrection(covdep::PolynomialThermoEnergyCoverageDependence,T,coverages) = zero(coverages[1])
+export PolynomialThermoEnergyCoverageDependence
\ No newline at end of file
diff --git a/src/Calculators/Thermovec.jl b/src/Calculators/Thermovec.jl
index 1c99db8bf..bc2a951b6 100644
--- a/src/Calculators/Thermovec.jl
+++ b/src/Calculators/Thermovec.jl
@@ -5,19 +5,22 @@ export AbstractThermovec
 
 include("Thermo.jl")
 
-@with_kw struct NASApolynomialvec <: AbstractThermo
+@with_kw struct NASApolynomialvec{Q} <: AbstractThermo
     coefs::Array{Float64,2}
     Tmin::Float64
     Tmax::Float64
+    covdeps::Q = nothing
 end
 
-@with_kw struct NASAvec{T<:AbstractThermoUncertainty} <: AbstractThermo
+@with_kw struct NASAvec{T<:AbstractThermoUncertainty,Q} <: AbstractThermo
     polys::Array{NASApolynomialvec,1}
     unc::T = EmptyThermoUncertainty()
+    covdeps::Q = nothing
 end
 function NASAvec(nasas::B) where {B<:Array}
     Tmin = nasas[1].polys[1].Tmin
     Tmax = nasas[1].polys[end].Tmax
+    covdeps = [nasa.covdep for nasa in nasas]
     Ts = Array{Float64,1}()
     for nasa in nasas
         for (i,poly) in enumerate(nasa.polys)
@@ -53,12 +56,10 @@ function NASAvec(nasas::B) where {B<:Array}
                 polyvec[:,j] = nasapoly.coefs
             end
         end
-        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1])
+        nasapvec = NASApolynomialvec(polyvec,Ts[i],Ts[i+1],covdeps)
         push!(polyvecs,nasapvec)
     end
-    
-    return NASAvec(polys=polyvecs) 
-        
+    return NASAvec(polys=polyvecs,covdeps=covdeps)
 end
 @inline function selectPoly(nasa::NASAvec,T::N) where {N<:Real}
     """
@@ -73,7 +74,7 @@ end
 end
 export selectPoly
 
-@inline function calcHSCpdless(poly::NASApolynomialvec,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomialvec{Array{1,EmptyThermoCoverageDependence}},T::X;coverages=nothing) where {X<:Real}
     if size(poly.coefs)[1] != 7
         Tpoly0 = T
         Tpoly1 = T*T
@@ -118,9 +119,63 @@ export selectPoly
     return (cpdivR,hdivRT,sdivR)
 end
 
-@inline function calcHSCpdless(nasavec::NASAvec,T::X) where {X<:Real}
+@inline function calcHSCpdless(poly::NASApolynomialvec{Q},T::X;coverages=nothing) where {X<:Real,Q}
+    hdivRT = similar(coverages)
+    sdivR = similar(coverages)
+    cpdivR = similar(coverages)
+    if size(poly.coefs)[1] != 7
+        Tpoly0 = T
+        Tpoly1 = T*T
+        Tpoly2 = Tpoly1*T
+        Tpoly3 = Tpoly2*T
+        Tpoly4 = 1.0/T
+        Tpoly5 = Tpoly4*Tpoly4
+        Tpoly6 = log(T)
+
+        ct0 = poly.coefs[1,:].*Tpoly5
+        ct1 = poly.coefs[2,:].*Tpoly4
+        ct2 = poly.coefs[3,:]
+        ct3 = poly.coefs[4,:].*Tpoly0
+        ct4 = poly.coefs[5,:].*Tpoly1
+        ct5 = poly.coefs[6,:].*Tpoly2
+        ct6 = poly.coefs[7,:].*Tpoly3
+
+        cpdivR .= ct0 .+ ct1 .+ ct2 .+ ct3 .+ ct4 .+ ct5 .+ ct6
+        hdivRT .= -ct0.+Tpoly6.*ct1.+ct2.+0.5.*ct3.+0.33333333333.*ct4.+0.25.*ct5+0.2.*ct6+poly.coefs[8,:].*Tpoly4
+        sdivR .= -0.5.*ct0 .- ct1 .+ Tpoly6.*ct2 .+ ct3 .+ 0.5.*ct4 .+ 0.33333333333.*ct5 .+ 0.25*ct6 .+ poly.coefs[9,:]
+    else
+        Tpoly0 = T
+        Tpoly1 = T*T
+        Tpoly2 = Tpoly1*T
+        Tpoly3 = Tpoly2*T
+        Tpoly4 = 1.0/T
+        #Tpoly5 = Tpoly4*Tpoly4
+        Tpoly6 = log(T)
+
+        #ct0 = 0.0
+        #ct1 = 0.0
+        ct2 = poly.coefs[1,:]
+        ct3 = poly.coefs[2,:].*Tpoly0
+        ct4 = poly.coefs[3,:].*Tpoly1
+        ct5 = poly.coefs[4,:].*Tpoly2
+        ct6 = poly.coefs[5,:].*Tpoly3
+
+        cpdivR .= ct2 .+ ct3 .+ ct4 .+ ct5 .+ ct6
+        hdivRT .= ct2.+0.5.*ct3.+0.33333333333.*ct4.+0.25.*ct5.+0.2.*ct6.+poly.coefs[6,:].*Tpoly4
+        sdivR .= Tpoly6.*ct2 .+ ct3 .+ 0.5.*ct4 .+ 0.33333333333.*ct5 .+ 0.25.*ct6 .+ poly.coefs[7,:]
+    end
+    for i in 1:length(hdivRT)
+        hdivRT .+= getcovdepenthalpycorrection(poly.covdeps[i],T,coverages)/(R*T)
+        sdivR .+= getcovdepentropycorrection(poly.covdeps[i],T,coverages)/R 
+        cpdivR .+= getcovdepheatcapacitycorrection(poly.covdeps[i],T,coverages)/R
+    end
+    return (cpdivR,hdivRT,sdivR)
+end
+
+
+@inline function calcHSCpdless(nasavec::NASAvec,T::X;coverages=nothing) where {X<:Real}
     poly = selectPoly(nasavec,T)
-    return calcHSCpdless(poly,T)
+    return calcHSCpdless(poly,T;coverages=coverages)
 end
 
 @with_kw struct ConstantGvec{B,J} <: AbstractThermo
diff --git a/src/Domain.jl b/src/Domain.jl
index 82d2e05ea..4d897b1c2 100644
--- a/src/Domain.jl
+++ b/src/Domain.jl
@@ -656,6 +656,7 @@ mutable struct ConstantTAPhiDomain{N<:AbstractPhase,S<:Integer,W<:Real,W2<:Real,
     t::MArray{Tuple{1},W2,1,1}
     p::Array{W,1}
     thermovariabledict::Dict{String,Int64}
+    iscovdep::Bool
 end
 function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspecies::Array{X3,1}=Array{String,1}(),
     sparse::Bool=false, sensitivity::Bool=false, stationary::Bool=false) where {E<:Real,E2<:AbstractPhase,W<:Real,X,X2,X3}
@@ -683,6 +684,7 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
     @assert T != 0.0
     ns = y0
     N = sum(ns)
+    coverages = ns ./ phase.sitedensity
 
     if length(constantspecies) > 0
         spcnames = getfield.(phase.species, :name)
@@ -691,15 +693,22 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
         constspcinds = Array{Int64,1}()
     end
     efficiencyinds = [rxn.index for rxn in phase.reactions if typeof(rxn.kinetics) <: AbstractFalloffRate && length(rxn.kinetics.efficiencies) > 0]
-    Gs = calcgibbs(phase, T)
+    Gs = calcgibbs(phase, T; coverages=coverages)
     if :solvent in fieldnames(typeof(phase)) && typeof(phase.solvent) != EmptySolvent
         mu = phase.solvent.mu(T)
     else
         mu = 0.0
     end
     C = 0.0 #this currently shouldn't matter here, on a surface you shouldn't have pdep
-    kfs, krevs = getkfkrevs(phase, T, 0.0, C, N, ns, Gs, [], A, phi)
-    p = vcat(Gs, kfs)
+    iscovdep = any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in phase.species]) || any([!(rxn.kinetics.covdep isa EmptyRateCoverageDependence) for rxn in phase.reactions])
+    kfs, krevs = getkfkrevs(phase, T, 0.0, C, N, ns, Gs, [], A, phi; coverages=coverages)
+    if !iscovdep
+        p = vcat(Gs, kfs)
+        alternativepformat = false 
+    else 
+        p = vcat(zeros(length(phase.species)), ones(length(phase.reactions)))
+        alternativepformat = true 
+    end 
     if sparse
         jacobian = spzeros(typeof(T), length(phase.species), length(phase.species))
     else
@@ -707,7 +716,7 @@ function ConstantTAPhiDomain(; phase::E2, initialconds::Dict{X,X2}, constantspec
     end
     rxnarray = getreactionindices(phase)
     return ConstantTAPhiDomain(phase, [1, length(phase.species)], [1, length(phase.species) + length(phase.reactions)], constspcinds,
-        T, A, phi, kfs, krevs, efficiencyinds, Gs, rxnarray, mu, Array{Float64,1}(), jacobian, sensitivity, false, MVector(false), MVector(0.0), p, Dict{String,Int64}()), y0, p
+        T, A, phi, kfs, krevs, efficiencyinds, Gs, rxnarray, mu, Array{Float64,1}(), jacobian, sensitivity, alternativepformat, MVector(false), MVector(0.0), p, Dict{String,Int64}(), iscovdep), y0, p
 end
 export ConstantTAPhiDomain
 
@@ -922,7 +931,7 @@ export ConstantTLiqFilmDomain
     for ind in d.efficiencyinds #efficiency related rates may have changed
         d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0)
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:Array{Float64,1},Q<:Float64} #uses parameter input
@@ -943,13 +952,13 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     elseif nothermochg
         d.kfs = kfps
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     else #need to handle thermo changes
         d.kfs .= kfps
         if !d.alternativepformat
@@ -961,7 +970,7 @@ end
         for ind in d.efficiencyinds #efficiency related rates may have changed
             d.kfs[ind], d.krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; f=kfps[ind])
         end
-        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+        return ns, cs, d.T, d.P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
     end
 end
 
@@ -982,7 +991,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J,Q} #Autodiff p
@@ -1002,7 +1011,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Autodiff p
@@ -1024,7 +1033,7 @@ end
     for ind in d.efficiencyinds #efficiency related rates may have changed
         kfs[ind], krevs[ind] = getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; f=kfs[ind])
     end
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J,Q} #Tracker/reversediff
@@ -1041,7 +1050,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},W<:IdealGas,Y<:Integer,J<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q} #Tracker/reversediff
@@ -1058,7 +1067,7 @@ end
         kfs = [ind in d.efficiencyinds ? getkfkrev(d.phase.reactions[ind], d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0)[1] * d.p[length(d.phase.species)+ind] * p[d.parameterindexes[1]-1+length(d.phase.species)+ind] : p[d.parameterindexes[1]-1+length(d.phase.species)+ind] for ind in 1:length(d.phase.reactions)]
     end
     krevs = getkfkrevs(d.phase, d.T, d.P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfs)[2]
-    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, d.P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1083,7 +1092,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, d.V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1109,7 +1118,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1135,7 +1144,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, d.V, 0.0)
-    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, d.V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1159,7 +1168,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1185,9 +1194,9 @@ end
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
     if p != SciMLBase.NullParameters()
-        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+        return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
     else
-        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+        return ns, cs, T, d.P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
     end
 end
 
@@ -1213,7 +1222,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, d.P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, d.P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1238,7 +1247,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1264,7 +1273,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1290,7 +1299,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Us, Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1315,7 +1324,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1341,7 +1350,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 @inline function calcthermo(d::ParametrizedPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
     P = d.P(t)
@@ -1366,7 +1375,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Hs, Array{Float64,1}(), Gs, diffs, Cvave, cpdivR, 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1387,7 +1396,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, T, P, V, C, N, mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1409,7 +1418,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTConstantVDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1431,7 +1440,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, d.phi)
-    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi
+    return @views @fastmath ns, cs, T, P, V, C, N, mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:SciMLBase.NullParameters,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1456,7 +1465,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2<:Array{Float64,1},W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1482,7 +1491,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ParametrizedTPDomain{W,Y}, y::J, t::Q, p::W2=SciMLBase.NullParameters()) where {W2,W<:IdealGas,Y<:Integer,J<:AbstractArray,Q}
@@ -1508,7 +1517,7 @@ end
         diffs = Array{Float64,1}()
     end
     kfs, krevs = getkfkrevs(d.phase, T, P, C, N, ns, Gs, diffs, V, 0.0)
-    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0
+    return @views @fastmath ns, cs, T, P, V, C, N, 0.0, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, diffs, 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q}
@@ -1517,7 +1526,7 @@ end
     cs = ns ./ d.V
     C = N / d.V
     P = 1.0e8
-    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1530,31 +1539,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, d.V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
         end
     end
 
@@ -1574,7 +1583,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTVDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1591,7 +1600,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, d.V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}()
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q}
@@ -1600,7 +1609,19 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    if d.iscovdep
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
+    else 
+        kfs = d.kfs 
+        krevs = d.krevs
+        Gs = d.Gs
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+    end
+    
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealSurface,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1609,26 +1630,35 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        if nothermochg
-            return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    if !d.iscovdep
+        if !d.alternativepformat
+            @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            if nothermochg
+                return ns, cs, d.T, P, d.A, C, N, d.mu, p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            else
+                d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+                d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+                d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            end
         else
-            d.kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-            d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-            d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            if nothermochg
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            else
+                d.kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+                d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+                d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
+                return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+            end
         end
     else
-        @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        if nothermochg
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)], d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, 0.0, Array{Float64,1}(), d.phi
-        else
-            d.kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-            d.Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-            d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, d.V, d.phi; kfs=d.kfs)[2]
-            return ns, cs, d.T, P, d.A, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
-        end
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R *d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
 end
 
@@ -1638,15 +1668,24 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = convert(typeof(y), p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
-    else
-        Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+    if !d.iscovdep
+        if !d.alternativepformat
+            Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = convert(typeof(y), p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+        else
+            Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
+        end
+        krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2])
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+    else 
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * d.T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
-    krevs = convert(typeof(y), getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2])
-    return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
 end
 
 @inline function calcthermo(d::ConstantTAPhiDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealSurface,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1655,15 +1694,24 @@ end
     cs = ns ./ d.A
     C = N / d.A
     P = 0.0
-    if !d.alternativepformat
-        Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
-    else
-        Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
-        kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+    if !d.iscovdep
+        if !d.alternativepformat
+            Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+        else
+            Gs = d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
+            kfs = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
+        end
+        krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2]
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, similar(cs,0,0)
+    else 
+        coverages = cs ./ d.phase.sitedensity
+        cpdivR, hdivRT, sdivR = calcHSCpdless(d.phase.vecthermo, d.T; coverages=coverages)
+        @fastmath @views hdivRT .+= p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)] ./ (R * T)
+        @fastmath Gs = (hdivRT .- sdivR) * (R * d.T)
+        kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, Array{Float64,1}(), d.A, d.phi; coverages=coverages)
+        return ns, cs, d.T, P, d.A, C, N, d.mu, kfs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], krevs .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(kfs)], Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, coverages
     end
-    krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, d.A, d.phi; kfs=kfs)[2]
-    return ns, cs, d.T, P, d.A, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:SciMLBase.NullParameters,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q}
@@ -1674,7 +1722,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
     cs = ns ./ V
     C = N / V
     P = 1.0e8
-    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:FragmentBasedIdealFilm,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1689,31 +1737,31 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}()
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, 0.0; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
         end
     end
 
@@ -1735,7 +1783,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::Array{W2,1}, t::
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
 end
 
 function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:FragmentBasedIdealFilm,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1754,7 +1802,7 @@ function calcthermo(d::FragmentBasedConstantTrhoDomain{W,Y}, y::J, t::Q, p::Q2=S
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, 0.0; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), 0.0, Array{Float64,1}
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2<:Array{Float64,1},W<:IdealDiluteSolution,Y<:Integer,J<:Array{Float64,1},Q}
@@ -1768,31 +1816,31 @@ end
         @views nothermochg = d.Gs == p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         else
             d.kfsnondiff = p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs = p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         end
     else
         @views nothermochg = d.Gs == d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
         @views nokfchg = d.kfsnondiff == d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
         if nothermochg && nokfchg
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         elseif nothermochg
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         else
             d.kfsnondiff .= d.p[length(d.phase.species)+1:length(d.phase.species)+length(d.phase.reactions)] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
             d.Gs .= d.p[1:length(d.phase.species)] .+ p[d.parameterindexes[1]-1+1:d.parameterindexes[1]-1+length(d.phase.species)]
             d.kfs, d.krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, d.Gs, d.diffusivity, V, d.phi; kfs=d.kfsnondiff)
-            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+            return ns, cs, d.T, P, V, C, N, d.mu, d.kfs, d.krevs, Array{Float64,1}(), Array{Float64,1}(), d.Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
         end
     end
 
@@ -1813,7 +1861,7 @@ end
         kfsnondiff = convert(typeof(y), d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)])
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
 end
 
 @inline function calcthermo(d::ConstantTLiqFilmDomain{W,Y}, y::J, t::Q, p::Q2=SciMLBase.NullParameters()) where {Q2,W<:IdealDiluteSolution,Y<:Integer,J<:AbstractArray,Q} #autodiff p
@@ -1831,7 +1879,7 @@ end
         kfsnondiff = d.p[length(d.phase.species)+1:end] .* p[d.parameterindexes[1]-1+length(d.phase.species)+1:d.parameterindexes[1]-1+length(d.phase.species)+length(d.phase.reactions)]
     end
     kfs, krevs = getkfkrevs(d.phase, d.T, P, C, N, ns, Gs, d.diffusivity, V, d.phi; kfs=kfsnondiff)
-    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi
+    return ns, cs, d.T, P, V, C, N, d.mu, kfs, krevs, Array{Float64,1}(), Array{Float64,1}(), Gs, Array{Float64,1}(), 0.0, Array{Float64,1}(), d.phi, Array{Float64,1}
 end
 
 export calcthermo
diff --git a/src/Interface.jl b/src/Interface.jl
index 466022abc..42f0da98e 100644
--- a/src/Interface.jl
+++ b/src/Interface.jl
@@ -35,6 +35,7 @@ struct ReactiveInternalInterface{T,B,C,C2,N,Q<:AbstractReaction,X} <: AbstractRe
     forwardability::Array{Bool,1}
 end
 function ReactiveInternalInterface(domain1, domain2, reactions, A)
+    reactions = upgradekinetics(reactions, domain1, domain2)
     vectuple, vecinds, otherrxns, otherrxninds, posinds = getveckinetics(reactions)
     rxns = vcat(reactions[vecinds], reactions[otherrxninds])
     rxns = [ElementaryReaction(index=i, reactants=rxn.reactants, reactantinds=rxn.reactantinds, products=rxn.products,
@@ -49,20 +50,28 @@ function ReactiveInternalInterface(domain1, domain2, reactions, A)
 end
 export ReactiveInternalInterface
 
-function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, cstot::Array{Q,1}) where {Q}
-    kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
+function getkfskrevs(ri::ReactiveInternalInterface, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot::Array{Q,1}) where {Q}
+    if !(coverages1 === nothing) && coverages2 === nothing
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1; coverages=coverages1)
+    elseif !(coverages2 === nothing) && coverages1 === nothing
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1; coverages=coverages2)
+    elseif coverages1 === nothing && coverages2 === nothing 
+        kfs = getkfs(ri, T1, 0.0, 0.0, Array{Q,1}(), ri.A, phi1)
+    else 
+        throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+    end
     Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), T1, phi1)
     krevs = kfs ./ Kc
     return kfs, krevs
 end
 
-function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
-    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
+function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p::W) where {W<:SciMLBase.NullParameters}
+    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
     addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
 end
 
-function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p)
-    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
+function evaluate(ri::ReactiveInternalInterface, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p)
+    kfs, krevs = getkfskrevs(ri, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
     addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs .* p[ri.parameterindexes[1]:ri.parameterindexes[2]], krevs .* p[ri.parameterindexes[1]:ri.parameterindexes[2]], ri.A)
 end
 export evaluate
@@ -85,12 +94,15 @@ struct ReactiveInternalInterfaceConstantTPhi{J,N,B,B2,B3,C,C2,Q<:AbstractReactio
     p::Array{Float64,1}
     reversibility::Array{Bool,1}
     forwardability::Array{Bool,1}
+    iscovdep::Bool
+    phi::Float64
 end
 function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A, phi=0.0)
     @assert domain1.T == domain2.T
     reactions = upgradekinetics(reactions, domain1, domain2)
     rxnarray = getinterfacereactioninds(domain1, domain2, reactions)
-    kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi)
+    iscovdep = any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in domain1.phase.species]) || any([!(spc.thermo.covdep isa EmptyThermoCoverageDependence) for spc in domain2.phase.species]) || any([!(rxn.kinetics.covdep isa EmptyRateCoverageDependence) for rxn in reactions])
+    kfs = getkf.(reactions, nothing, T, 0.0, 0.0, Ref([]), A, phi; coverages=nothing)
     Kc = getKc.(reactions, domain1.phase, domain2.phase, Ref(domain1.Gs), Ref(domain2.Gs), T, phi)
     krevs = kfs ./ Kc
     M, Nrp1, Nrp2 = getstoichmatrix(domain1, domain2, reactions)
@@ -102,33 +114,87 @@ function ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions, T, A
     if isa(kfs, Vector{Any})
         kfs = convert(Vector{Float64}, kfs)
     end
-    return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
-        rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
-        [0, 1], kfs[1:end], reversibility, forwardability), kfs[1:end]
+    if !iscovdep
+        return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
+            rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
+            [0, 1], kfs[1:end], reversibility, forwardability, iscovdep, phi), kfs[1:end]
+    else 
+        return ReactiveInternalInterfaceConstantTPhi(domain1, domain2, reactions,
+            rxnarray, M, Nrp1, Nrp2, kfs, krevs, T, A, [1, length(reactions)],
+            [0, 1], kfs[1:end], reversibility, forwardability, iscovdep, phi), ones(length(kfs))
+    end
 end
 export ReactiveInternalInterfaceConstantTPhi
 
-function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi, T1, T2, phi1, phi2, Gs1, Gs2, cstot)
-    return ri.kfs, ri.krevs
-end
-
-function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p::W) where {W<:SciMLBase.NullParameters}
-    addreactionratecontributions!(dydt, ri.rxnarray, cstot, ri.kfs, ri.krevs, ri.A)
+function getkfskrevs(ri::ReactiveInternalInterfaceConstantTPhi, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot)
+    if ri.iscovdep
+        if length(coverages1) != 0 && length(coverages2) == 0
+            coverages = coverages1
+        elseif length(coverages2) != 0 && length(coverages1) == 0
+            coverages = coverages2
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, ri.phi)
+        krevs = kfs ./ Kc
+        return kfs,krevs
+    else
+        return ri.kfs, ri.krevs
+    end
 end
 
-function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, cstot, p)
-    if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
-        kfs = ri.kfs
+function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p::W) where {W<:SciMLBase.NullParameters}
+    if ri.iscovdep
+        if length(coverages1) != 0 && length(coverages2) == 0
+            coverages = coverages1
+        elseif length(coverages2) != 0 && length(coverages1) == 0
+            coverages = coverages2
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages)
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, phi1)
+        krevs = kfs ./ Kc
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
     else
-        kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, ri.kfs, ri.krevs, ri.A)
     end
-    if length(Gs1) == 0 || length(Gs2) == 0 || (all(Gs1 .== ri.domain1.Gs) && all(Gs2 .== ri.domain2.Gs))
-        krevs = ri.krevs
-    else
-        Kc = getKcs(ri, T1, Gs1, Gs2)
+end
+
+function evaluate(ri::ReactiveInternalInterfaceConstantTPhi, dydt, domains, T1, T2, phi1, phi2, Gs1, Gs2, coverages1, coverages2, cstot, p)
+    if ri.iscovdep
+        if length(coverages1) != 0 && length(coverages2) == 0
+            coverages = coverages1
+        elseif length(coverages2) != 0 && length(coverages1) == 0
+            coverages = coverages2
+        elseif length(coverages1) == 0 && length(coverages2) == 0 
+            coverages = coverages1
+        else
+            throw(DomainError("Coverage dependence of an interface between two surfaces is indeterminate (which coverage should be used?)"))
+        end
+        kfs = getkf.(ri.reactions, nothing, ri.T, 0.0, 0.0, Ref([]), ri.A, ri.phi; coverages=coverages) .* p[ri.parameterindexes[1]:ri.parameterindexes[2]] 
+        Kc = getKc.(ri.reactions, ri.domain1.phase, ri.domain2.phase, Ref(Gs1), Ref(Gs2), ri.T, ri.phi)
         krevs = kfs ./ Kc
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
+    else
+        if p[ri.parameterindexes[1]:ri.parameterindexes[2]] == ri.kfs
+            kfs = ri.kfs
+        else
+            kfs = p[ri.parameterindexes[1]:ri.parameterindexes[2]]
+        end
+        if length(Gs1) == 0 || length(Gs2) == 0 || (all(Gs1 .== ri.domain1.Gs) && all(Gs2 .== ri.domain2.Gs))
+            krevs = ri.krevs
+        else
+            Kc = getKcs(ri, T1, Gs1, Gs2)
+            krevs = kfs ./ Kc
+        end
+        addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
     end
-    addreactionratecontributions!(dydt, ri.rxnarray, cstot, kfs, krevs, ri.A)
 end
 export evaluate
 
@@ -204,6 +270,31 @@ function upgradekinetics(rxns, domain1, domain2)
     end
     newrxns = Array{ElementaryReaction,1}(undef, length(rxns))
     for (i, rxn) in enumerate(rxns)
+        if hasproperty(rxn.kinetics,:covdep) && !(rxn.kinetics.covdep isa EmptyRateCoverageDependence)
+            if domain1.phase isa IdealSurface && !(domain2.phase isa IdealSurface)
+                names = getfield.(domain1.phase.species,:name)
+            elseif domain2.phase isa IdealSurface && !(domain1.phase isa IdealSurface)
+                names = getfield.(domain2.phase.species,:name)
+            else 
+                throw(DomainError(rxn.kinetics.covdep,"Unsure which surface domain should be used for coverage dependence"))
+            end
+            if rxn.kinetics.covdep isa PolynomialRateCoverageDependence
+                for (name,v) in rxn.kinetics.covdep.Epolys
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indEpolys[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.avals
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indavals[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.ms
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indms[ind] = v
+                end 
+            else
+                throw(TypeError(rxn.kinetics.covdep,"Kinetic Coverage Dependence Type Not Understood"))
+            end
+        end
         if isa(rxn.kinetics, StickingCoefficient)
             spc = [spc for spc in rxn.reactants if !(spc in surfdomain.phase.species)]
             @assert length(spc) == 1
@@ -228,7 +319,7 @@ function stickingcoefficient2arrhenius(sc, sitedensity, N, mw; Tmin=300.0, Tmax=
     @assert fit.converged
     p = fit.param
     p[1] = abs(p[1])
-    return Arrhenius(; A=p[1], n=p[2], Ea=p[3])
+    return Arrhenius(; A=p[1], n=p[2], Ea=p[3], covdep=sc.covdep)
 end
 
 struct Inlet{Q<:Real,S,V<:AbstractArray,U<:Real,X<:Real,FF<:Function} <: AbstractBoundaryInterface
diff --git a/src/Parse.jl b/src/Parse.jl
index 3698081f3..5251e06bb 100644
--- a/src/Parse.jl
+++ b/src/Parse.jl
@@ -4,7 +4,9 @@ using PyCall
 using StaticArrays
 
 module Calc
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
@@ -15,6 +17,7 @@ module Calc
     include("Calculators/kLAkH.jl")
 end
 module Spc
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
@@ -22,7 +25,9 @@ module Spc
     include("Species.jl")
 end
 module Rxn
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
diff --git a/src/Phase.jl b/src/Phase.jl
index 784ed5c80..5165f7456 100644
--- a/src/Phase.jl
+++ b/src/Phase.jl
@@ -317,9 +317,23 @@ export iterate
 Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
 export broadcastable
 
-function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
+function getreactionindices(spcs,rxns)
     arr = zeros(Int64,(8,length(rxns)))
     names = [spc.name for spc in spcs]
+    for spc in spcs
+        if hasproperty(spc.thermo,:covdep)
+            if spc.thermo.covdep isa PolynomialThermoEnergyCoverageDependence
+                for (name,v) in spc.thermo.covdep.polys
+                    ind = findfirst(isequal(name),names)
+                    spc.thermo.covdep.indpolys[ind] = v
+                end 
+            elseif spc.thermo.covdep isa EmptyThermoCoverageDependence
+
+            else 
+                throw(TypeError(spc.thermo.covdep,"Thermo Coverage Dependence Type Not Understood"))
+            end
+        end
+    end
     for (i,rxn) in enumerate(rxns)
         inds = [findfirst(isequal(spc),spcs) for spc in rxn.reactants]
         for (j,spc) in enumerate(rxn.reactants)
@@ -343,6 +357,28 @@ function getreactionindices(spcs,rxns) where {Q<:AbstractPhase}
                 end
             end
         end
+        if hasproperty(rxn.kinetics,:covdep) 
+            if rxn.kinetics.covdep isa PolynomialRateCoverageDependence
+                for (name,v) in rxn.kinetics.covdep.Epolys
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indEpolys[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.avals
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indavals[ind] = v
+                end 
+                for (name,v) in rxn.kinetics.covdep.ms
+                    ind = findfirst(isequal(name),names)
+                    rxn.kinetics.covdep.indms[ind] = v
+                end 
+            elseif rxn.kinetics.covdep isa EmptyRateCoverageDependence
+
+            else 
+                throw(TypeError(rxn.kinetics.covdep,"Kinetic Coverage Dependence Type Not Understood"))
+            end
+
+                
+        end
     end
     return arr
 end
diff --git a/src/PhaseState.jl b/src/PhaseState.jl
index ff2683e8f..f0922d0d0 100644
--- a/src/PhaseState.jl
+++ b/src/PhaseState.jl
@@ -5,23 +5,23 @@ using Tracker
 using ReverseDiff
 using RecursiveArrayTools
 
-@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
-    return getGibbs.(getfield.(ph.species,:thermo),T)
+@inline function calcgibbs(ph::U,T::W;coverages=nothing) where {U<:IdealPhase,W<:Real}
+    return getGibbs.(getfield.(ph.species,:thermo),T,coverages=coverages)
 end
 export calcgibbs
 
-@inline function calcenthalpyinternalgibbs(ph::U,T::W,P::Z,V::Q) where {U<:IdealPhase,W,Z,Q<:Real}
-    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
+@inline function calcenthalpyinternalgibbs(ph::U,T::W,P::Z,V::Q;coverages=nothing) where {U<:IdealPhase,W,Z,Q<:Real}
+    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T,coverages=coverages)
     Us = Hs .- P*V
-    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
+    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T,coverages=coverages)
     return Hs,Us,Gs
 end
 
 
-@inline function calcenthalpyinternalgibbs(ph::Union{IdealGas,IdealSurface},T::W,P::Z,V::Q) where {W,Z,Q<:Real}
-    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
+@inline function calcenthalpyinternalgibbs(ph::Union{IdealGas,IdealSurface},T::W,P::Z,V::Q;coverages=nothing) where {W,Z,Q<:Real}
+    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T,coverages)
     Us = Hs .- R*T
-    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
+    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T,coverages)
     return Hs,Us,Gs
 end
 
@@ -50,27 +50,27 @@ end
 
 export makespcsvector
 
-@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi)
+@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi;coverages=nothing)
     if isdefined(rxn.kinetics,:efficiencies) && length(rxn.kinetics.efficiencies) > 0
         @views @inbounds @fastmath C += sum([ns[i]*val for (i,val) in rxn.kinetics.efficiencies])/V
     end
-    return rxn.kinetics(T=T,P=P,C=C,phi=phi)
+    return rxn.kinetics(T=T,P=P,C=C,phi=phi,coverages=coverages)
 end
 export getkf
 
-@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
+@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi;coverages=nothing) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
     kfs = similar(ns,length(ph.reactions))
     i = 1
     oldind = 1
     ind = 0
     while i <= length(ph.veckineticsinds) #vectorized kinetics
         @inbounds ind = ph.veckineticsinds[i]
-        @inbounds kfs[oldind:ind] = ph.veckinetics[i](;T=T,P=P,C=C)
+        @inbounds kfs[oldind:ind] = ph.veckinetics[i](;T=T,P=P,C=C,coverages=coverages)
         oldind = ind+1
         i += 1
     end
     @simd for i in ind+1:length(ph.reactions)
-        @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi)
+        @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi,coverages=coverages)
     end
     return kfs
 end
@@ -181,12 +181,12 @@ export getKcs
 Calculates the forward and reverse rate coefficients for a given reaction, phase and state
 Maintains diffusion limitations if the phase has diffusionlimited=true
 """
-@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
+@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;coverages=nothing,kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
     if signbit(kf) 
         if signbit(f)
-            kf = getkf(rxn,ph,T,P,C,ns,V,phi)
+            kf = getkf(rxn,ph,T,P,C,ns,V,phi;coverages=coverages)
         else
-            kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
+            kf = getkf(rxn,ph,T,P,C,ns,V,phi;coverages=coverages)*f
         end
     end
     Kc = getKc(rxn,ph,T,Gs,phi)
@@ -223,9 +223,9 @@ Maintains diffusion limitations if the phase has diffusionlimited=true
 end
 export getkfkrev
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else 
@@ -241,14 +241,14 @@ export getkfkrev
         len = length(phase.reactions)
         krev = zeros(typeof(N),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
     else
         len = length(phase.reactions)
         kfs = zeros(typeof(N),len)
         krev = zeros(typeof(N),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     kfs .*= phase.forwardability
@@ -256,9 +256,9 @@ export getkfkrev
     return kfs,krev
 end
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else 
@@ -275,7 +275,7 @@ end
         krev = similar(kfs)
         kfsdiff = similar(kfs)
         @simd for i = 1:len
-           @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
         return kfsdiff, krev
     else
@@ -283,15 +283,15 @@ end
         kfs = zeros(typeof(Gs[1]),len)
         krev = zeros(typeof(Gs[1]),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     return kfs,krev
 end
 
-@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
+@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;coverages=nothing,kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
     if !phase.diffusionlimited && kfs === nothing
-        kfs = getkfs(phase,T,P,C,ns,V,phi)
+        kfs = getkfs(phase,T,P,C,ns,V,phi;coverages=coverages)
         if phi == 0.0
             krev = @fastmath kfs./getKcs(phase,T,Gs)
         else 
@@ -307,14 +307,14 @@ end
         len = length(phase.reactions)
         krev = similar(kfs)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages,kf=kfs[i])
         end
     else
         len = length(phase.reactions)
         kfs = zeros(typeof(Gs[1]),len)
         krev = zeros(typeof(Gs[1]),len)
         @simd for i = 1:len
-           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
+           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;coverages=coverages)
         end
     end
     kfs .*= phase.forwardability
diff --git a/src/ReactionMechanismSimulator.jl b/src/ReactionMechanismSimulator.jl
index 6e30c01dd..40e520141 100644
--- a/src/ReactionMechanismSimulator.jl
+++ b/src/ReactionMechanismSimulator.jl
@@ -52,7 +52,9 @@ module ReactionMechanismSimulator
     end
     include("Constants.jl")
     include("Tools.jl")
+    include("Calculators/RateCoverageDependence.jl")
     include("Calculators/RateUncertainty.jl")
+    include("Calculators/ThermoCoverageDependence.jl")
     include("Calculators/ThermoUncertainty.jl")
     include("Calculators/Thermo.jl")
     include("Calculators/Diffusion.jl")
diff --git a/src/Reactor.jl b/src/Reactor.jl
index 549e9ead4..53b79aaf8 100644
--- a/src/Reactor.jl
+++ b/src/Reactor.jl
@@ -647,7 +647,7 @@ end
     cstot .= 0.0
     dydt .= 0.0
     domain = domains[1]
-    ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi = calcthermo(domain, y, t, p)
+    ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi, coverages = calcthermo(domain, y, t, p)
     vns = Array{Any,1}(undef, length(domains))
     vns[1] = ns
     vcs = Array{Any,1}(undef, length(domains))
@@ -683,11 +683,13 @@ end
     vcpdivR[1] = cpdivR
     vphi = Array{Any,1}(undef, length(domains))
     vphi[1] = phi
+    vcoverages = Array{Any,1}(undef, length(domains))
+    vcoverages[1] = coverages
     addreactionratecontributions!(dydt, domain.rxnarray, cstot, kfs, krevs)
     @views dydt[domain.indexes[1]:domain.indexes[2]] .*= V
     for (i, domain) in enumerate(@views domains[2:end])
         k = i + 1
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k] = calcthermo(domain, y, t, p)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(domain, y, t, p)
         cstot[domain.indexes[1]:domain.indexes[2]] .= vcs[k]
         addreactionratecontributions!(dydt, domain.rxnarray, cstot, vkfs[k], vkrevs[k])
         @views dydt[domain.indexes[1]:domain.indexes[2]] .*= vV[k]
@@ -696,7 +698,7 @@ end
         if isa(inter, FragmentBasedReactiveFilmGrowthInterfaceConstantT)
             evaluate(inter, dydt, vV[inter.domaininds[1]], cstot)
         elseif isa(inter, AbstractReactiveInternalInterface)
-            evaluate(inter, dydt, domains, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot, p)
+            evaluate(inter, dydt, domains, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot, p)
         elseif isa(inter, VaporLiquidMassTransferInternalInterfaceConstantT)
             evaluate(inter, dydt, vV[inter.domaininds[1]], vV[inter.domaininds[2]], vT[inter.domaininds[1]], vT[inter.domaininds[2]], vN[inter.domaininds[1]], vN[inter.domaininds[2]], vP[inter.domaininds[1]], vP[inter.domaininds[2]], vCvave[inter.domaininds[1]], vCvave[inter.domaininds[2]], vns[inter.domaininds[1]], vns[inter.domaininds[2]], vUs[inter.domaininds[1]], vUs[inter.domaininds[2]], cstot, p)
         end
diff --git a/src/Simulation.jl b/src/Simulation.jl
index f80f614d4..653ad5911 100644
--- a/src/Simulation.jl
+++ b/src/Simulation.jl
@@ -53,7 +53,7 @@ function Simulation(sol::Q, domain::W, reducedmodelmappings::ReducedModelMapping
         end
         @inbounds @views yunlumped[end-length(domain.thermovariabledict)+1:end] .= y[end-length(domain.thermovariabledict)+1:end]
 
-        ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi = calcthermo(domain, yunlumped, t, p)
+        ns, cs, T, P, V, C, N, mu, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi, coverages = calcthermo(domain, yunlumped, t, p)
 
         reducedmodelmappings.qssc!(qssc, cs, kfs, krevs)
         @inbounds yunlumped[reducedmodelmappings.qssindexes] .= qssc .* V
@@ -272,10 +272,12 @@ function rops(ssys::SystemSimulation, t)
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     ropmat = spzeros(Nrxns, Nspcs)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropmat, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start)
         start += length(vkfs[k])
@@ -286,7 +288,7 @@ function rops(ssys::SystemSimulation, t)
             rops!(ropmat, nothing, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], start)
             start += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
             start += length(kfs)
         end
@@ -316,7 +318,9 @@ function massrops(ssys::SystemSimulation, t)
     vGs = Array{Any,1}(undef, length(domains))
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     if any([domain isa FragmentBasedConstantTrhoDomain for domain in domains])
         ropmat = spzeros(Nrxns, Nspcs + 1)
     else
@@ -326,7 +330,7 @@ function massrops(ssys::SystemSimulation, t)
     ropmat = spzeros(Nrxns, Nspcs)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropmat, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start)
         start += length(vkfs[k])
@@ -337,7 +341,7 @@ function massrops(ssys::SystemSimulation, t)
             rops!(ropmat, massropvec, inter.rxnarray, inter.fragmentbasedrxnarray, cstot, kfs, krevs, vV[inter.domaininds[1]], inter.Mws, inter.domainfilm.indexes[1]:inter.domainfilm.indexes[2], start)
             start += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropmat, inter.rxnarray, cstot, kfs, krevs, inter.A, start)
             start += length(kfs)
         end
@@ -392,17 +396,19 @@ function rops(ssys::SystemSimulation, name, t)
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     ropvec = spzeros(Nrxns)
     start = 0
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rops!(ropvec, sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k], vV[k], start, ind)
         start += length(vkfs[k])
     end
     for inter in ssys.interfaces
         if hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rops!(ropvec, inter.rxnarray, cstot, kfs, krevs, inter.A, start, ind)
             start += length(kfs)
         end
@@ -800,9 +806,11 @@ function rates(ssys::Q, t::X) where {Q<:SystemSimulation,X<:Real}
     vdiffs = Array{Any,1}(undef, length(domains))
     vCvave = Array{Any,1}(undef, length(domains))
     vphi = Array{Any,1}(undef, length(domains))
+    vcpdivR = Array{Any,1}(undef, length(domains))
+    vcoverages = Array{Any,1}(undef, length(domains))
     index = 1
     for (k, sim) in enumerate(ssys.sims)
-        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vphi[k] = calcthermo(sim.domain, ssys.sol(t), t)
+        vns[k], vcs[k], vT[k], vP[k], vV[k], vC[k], vN[k], vmu[k], vkfs[k], vkrevs[k], vHs[k], vUs[k], vGs[k], vdiffs[k], vCvave[k], vcpdivR[k], vphi[k], vcoverages[k] = calcthermo(sim.domain, ssys.sol(t), t)
         cstot[sim.domain.indexes[1]:sim.domain.indexes[2]] = vcs[k]
         rts[index:index+length(vkfs[k])-1] .= getrates(sim.domain.rxnarray, cstot, vkfs[k], vkrevs[k]) .* getdomainsize(sim, t)
         index += length(vkfs[k])
@@ -813,7 +821,7 @@ function rates(ssys::Q, t::X) where {Q<:SystemSimulation,X<:Real}
             @views rts[index:index+length(kfs)-1] = getrates(inter.rxnarray, cstot, kfs, krevs) .* vV[inter.domaininds[1]]
             index += length(kfs)
         elseif hasproperty(inter, :reactions)
-            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], cstot)
+            kfs, krevs = getkfskrevs(inter, vT[inter.domaininds[1]], vT[inter.domaininds[2]], vphi[inter.domaininds[1]], vphi[inter.domaininds[2]], vGs[inter.domaininds[1]], vGs[inter.domaininds[2]], vcoverages[inter.domaininds[1]], vcoverages[inter.domaininds[2]], cstot)
             rts[index:index+length(kfs)-1] = getrates(inter.rxnarray, cstot, kfs, krevs) .* inter.A
             index += length(kfs)
         end
diff --git a/src/rmstest.jl b/src/rmstest.jl
index 1ffa8bcce..f1a939c4f 100644
--- a/src/rmstest.jl
+++ b/src/rmstest.jl
@@ -19,7 +19,9 @@ copy!(pydot,pyimport_conda("pydot","pydot","rmg"))
 
 include("Constants.jl")
 include("Tools.jl")
+include("Calculators/RateCoverageDependence.jl")
 include("Calculators/RateUncertainty.jl")
+include("Calculators/ThermoCoverageDependence.jl")
 include("Calculators/ThermoUncertainty.jl")
 include("Calculators/Thermo.jl")
 include("Calculators/Diffusion.jl")