diff --git a/parcel/parcel.jl b/parcel/parcel.jl
index ddb1c8981e..8122ddd2b7 100644
--- a/parcel/parcel.jl
+++ b/parcel/parcel.jl
@@ -76,17 +76,27 @@ function parcel_model(dY, Y, p, t)
     dN_act_dt_depo = FT(0)
     dqi_dt_new_depo = FT(0)
     if "MohlerAF_Deposition" in ice_nucleation_modes
-        AF = CMI_het.dust_activated_number_fraction(
-            aerosol,
-            ip.deposition,
-            S_i,
-            T,
-        )
+        if S_i < ip.deposition.Sᵢ_max
+            println("Supersaturation exceeds the allowed value. No dust will be activated.")
+            AF = FT(0)
+        else
+            AF = CMI_het.dust_activated_number_fraction(
+                aerosol,
+                ip.deposition,
+                S_i,
+                T,
+            )
+        end
         dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt
         dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
     end
     if "MohlerRate_Deposition" in ice_nucleation_modes
-        dN_act_dt_depo = CMI_het.MohlerDepositionRate(aerosol, ip.deposition, S_i, T, (dY[1] * ξ), N_aer)
+        if S_i < ip.deposition.Sᵢ_max
+            println("Supersaturation exceeds the allowed value. No dust will be activated.")
+            dN_act_dt_depo = FT(0)
+        else
+            dN_act_dt_depo = CMI_het.MohlerDepositionRate(aerosol, ip.deposition, S_i, T, (dY[1] * ξ), N_aer)
+        end
         dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
     end