diff --git a/src/actors/current/current_actor.jl b/src/actors/current/current_actor.jl
index 45ac1839d..e3eb7b9c1 100644
--- a/src/actors/current/current_actor.jl
+++ b/src/actors/current/current_actor.jl
@@ -5,7 +5,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorCurrent{T<:Real} <: ParametersAc
     _parent::WeakRef = WeakRef(nothing)
     _name::Symbol = :not_set
     _time::Float64 = NaN
-    model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :none], "-", "Current actor to run"; default=:SteadyStateCurrent)
+    model::Switch{Symbol} = Switch{Symbol}([:SteadyStateCurrent, :QED, :none], "-", "Current actor to run"; default=:QED)
     allow_floating_plasma_current::Entry{Bool} = Entry{Bool}("-", "Zero loop voltage if non-inductive fraction exceeds 100% of the target Ip"; default=true)
     #== data flow parameters ==#
     ip_from::Switch{Symbol} = switch_get_from(:ip)
diff --git a/src/actors/current/qed_actor.jl b/src/actors/current/qed_actor.jl
index 3b04d9887..90fe50d1c 100644
--- a/src/actors/current/qed_actor.jl
+++ b/src/actors/current/qed_actor.jl
@@ -14,6 +14,7 @@ Base.@kwdef mutable struct FUSEparameters__ActorQED{T<:Real} <: ParametersActor{
     #== data flow parameters ==#
     ip_from::Switch{Symbol} = switch_get_from(:ip)
     vloop_from::Switch{Symbol} = switch_get_from(:vloop)
+    qmin_desired::Entry{Float64} = Entry{Float64}("-", "Desired minimum magnitude of q-profile"; default=1.0)
 end
 
 mutable struct ActorQED{D,P} <: SingleAbstractActor{D,P}
@@ -120,8 +121,20 @@ function _step(actor::ActorQED)
             Ip = nothing
             Vedge = IMAS.get_from(dd, Val{:vloop}, par.vloop_from)
         end
+        B0 = eqt.global_quantities.vacuum_toroidal_field.b0
 
         actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]); Vedge, Ip)
+        j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0
+
+        # using Jardin's model for stationary sawteeth to raise q>1
+        # requires running QED twice (the first time to evaluate q, the second to relax using updated q profile)
+        if par.qmin_desired > 0  
+            qval = 1.0 ./ abs.(actor.QO.ι.(cp1d.grid.rho_tor_norm))
+            i_qdes = findlast(qval .< par.qmin_desired)
+
+            actor.QO = QED.steady_state(actor.QO, η_imas(dd.core_profiles.profiles_1d[]; i_qdes); Vedge, Ip)
+            j_total = QED.JB(actor.QO; ρ=cp1d.grid.rho_tor_norm) ./ B0
+        end
     end
 
     return actor
@@ -191,8 +204,10 @@ function qed_init_from_imas(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_p
     return QED.initialize(rho_tor, B0, gm1, f, dvolume_drho_tor, q, j_tor, gm9; ρ_j_non_inductive, ρ_grid)
 end
 
-function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true)
-    rho = cp1d.grid.rho_tor_norm
+function η_imas(cp1d::IMAS.core_profiles__profiles_1d; use_log::Bool=true, i_qdes::Int=0)
     η = 1.0 ./ cp1d.conductivity_parallel
-    return QED.η_FE(rho, η; use_log)
+    if i_qdes != 0
+        η[1:i_qdes] .= η[i_qdes]
+    end
+    return QED.η_FE(cp1d.grid.rho_tor_norm, η; use_log)
 end