diff --git a/src/relations.jl b/src/relations.jl
index 5745b16c..64b19a43 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -222,6 +222,24 @@ vapor_specific_humidity(q::PhasePartition) = max(0, q.tot - q.liq - q.ice)
 vapor_specific_humidity(param_set::APS, ts::ThermodynamicState) =
     vapor_specific_humidity(PhasePartition(param_set, ts))
 
+"""
+    ∂q_vap_sat_∂T(param_set, ts)
+
+The change in saturation vapor specific humidity with temperature given by the Clausius Clapeyron relation
+ - `param_set` an `AbstractParameterSet`, see the [`Thermodynamics`](@ref) for more details
+ - `ts` ThermodynamicState
+"""
+function ∂q_vap_sat_∂T(param_set::APS, ts::ThermodynamicState)
+    R_v::FT = ICP.R_v(param_set)
+    L_v = latent_heat_vapor(param_set, ts)
+    L_s = latent_heat_sublim(param_set, ts)
+    λ = liquid_fraction(param_set, ts)
+    L = λ * L_v + (1 - λ) * L_s
+    T = air_temperature(param_set, ts)
+    q_vap_sat = vapor_specific_humidity(param_set, ts)
+    return L * q_vap_sat / (R_v * T^2)
+end
+
 """
     cp_m(param_set, q::PhasePartition)
 
diff --git a/test/relations.jl b/test/relations.jl
index 02bc4665..19b14a43 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -516,6 +516,16 @@ end
         @test all(saturated.(param_set, ts[RH_sat_mask]))
         @test !any(saturated.(param_set, ts[RH_unsat_mask]))
 
+        # test Clausius Clapeyron relation
+        ts = PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot)
+        L_v = TD.latent_heat_vapor.(param_set, ts)
+        L_s = TD.latent_heat_sublim.(param_set, ts)
+        λ = TD.liquid_fraction.(param_set, ts)
+        L = λ .* L_v .+ (1 .- λ) .* L_s
+        T = TD.air_temperature.(param_set, ts)
+        q_vap_sat = TD.vapor_specific_humidity(param_set, ts)
+        @test all(TD.∂q_vap_sat_∂T(param_set, ts) .≈ L .^ 2 .* q_vap_sat ./ (R_v .* T.^2))
+
         # PhaseEquil (freezing)
         _T_freeze = FT(TP.T_freeze(param_set))
         e_int_upper =