diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index 9f23be3b..ac1e073a 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -119,6 +119,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             {
                 phis[i] = phi;
                 phis[i] += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf();
+
+                if (relaxDrift < 1.0 && relaxDrift > 0.0)
+                {
+                    phis[i] *= relaxDrift;
+                }
             }
 
             if (ions.contains(Y[i].name()))
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index d48e25cc..1bfe5605 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -12,6 +12,11 @@ IOdictionary physicalProperties
     )
 );
 
+scalar relaxDrift
+(
+    physicalProperties.lookupOrDefault("relaxDrift", 1.0)
+);
+
 dimensionedScalar epsilon0
 (
     physicalProperties.lookup("epsilon0")
diff --git a/applications/solvers/combustion/plasmaReactingFoam/neEqn.H b/applications/solvers/combustion/plasmaReactingFoam/neEqn.H
index 5c551986..dd1978f5 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/neEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/neEqn.H
@@ -33,6 +33,10 @@
 
 
     Udrift = -(mue/ng)*E;
+    if (relaxDrift < 1.0 && relaxDrift > 0.0)
+    {
+        Udrift *= relaxDrift;
+    }
     Uthermal = -((De/ng/Te)*fvc::grad(Te));
     ve = (linearInterpolate(Udrift+Uthermal) & mesh.Sf()) + q;