diff --git a/src/Applications/cavity.jl b/src/Applications/cavity.jl
index 38f36aa..1dec97e 100644
--- a/src/Applications/cavity.jl
+++ b/src/Applications/cavity.jl
@@ -38,15 +38,16 @@ function _cavity(;
   distribute=nothing,
   np=1,
   nc=(4,4,4),
-  ν=1.0,
-  ρ=1.0,
-  σ=1.0,
-  ζ=0.0,
-  B=VectorValue(0.0, 0.0, 10.0),
-  f=VectorValue(0.0, 0.0, 0.0),
-  L=1.0,
-  u0=1.0,
-  B0=norm(B),
+  ν = 1.0,
+  ρ = 1.0,
+  σ = 1.0,
+  ζ = 0.0,   # Augmented Lagrangian weight  
+  μ = 0,     # Stabilization weight
+  B = VectorValue(0.0, 0.0, 10.0),
+  f = VectorValue(0.0, 0.0, 0.0),
+  L = 1.0,
+  u0 = 1.0,
+  B0 = norm(B),
   order = 2,
   order_j = order,
   formulation = :mhd,
@@ -140,6 +141,10 @@ function _cavity(;
     )
   end
 
+  if μ > 0
+    params[:bcs][:stabilization] = Dict(:μ=>μ)
+  end
+
   if !uses_petsc(params[:solver])
     xh,fullparams,info = main(params;output=info)
   else
diff --git a/src/Applications/expansion.jl b/src/Applications/expansion.jl
index 5903894..c0d2759 100644
--- a/src/Applications/expansion.jl
+++ b/src/Applications/expansion.jl
@@ -48,16 +48,16 @@ function _expansion(;
   verbose = true,
   solver  = :julia,
   formulation = :mhd,
-  Z = 4.0,                  # Expansion Ratio, it has to be consistent with the mesh
-  b = 0.2,                  # Outlet channel aspect ratio, it has to be consistent with the mesh
+  Z = 4.0,    # Expansion Ratio, it has to be consistent with the mesh
+  b = 0.2,    # Outlet channel aspect ratio, it has to be consistent with the mesh
   N  = 1.0,
   Ha = 1.0,
   cw = 0.028,
   τ  = 100,
-  ζ  = 0.0,
+  ζ  = 0.0,   # Augmented Lagrangian weight  
+  μ  = 0,     # Stabilization weight
   order = 2,
   order_j = order,
-  μ = 0,
   inlet = :parabolic,
   initial_value = :zero,
   solid_coupling = :none,
diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl
index 0aa5677..bee6374 100644
--- a/src/Applications/hunt.jl
+++ b/src/Applications/hunt.jl
@@ -47,7 +47,8 @@ function _hunt(;
   σ=1.0,
   B=(0.0,10.0,0.0),
   f=(0.0,0.0,1.0),
-  ζ=0.0,
+  ζ  = 0.0, # Augmented Lagrangian weight  
+  μ  = 0,   # Stabilization weight
   L=1.0,
   u0=1.0,
   B0=norm(VectorValue(B)),
@@ -162,6 +163,10 @@ function _hunt(;
     :j => Dict(:tags=>"insulating"),
   )
 
+  if μ > 0
+    params[:bcs][:stabilization] = Dict(:μ=>μ)
+  end
+
   toc!(t,"pre_process")
 
   # Solve it