diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index 5ed001d5..081bd630 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -21,7 +21,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     forAll(ions, k) // ion-neutral pair
     {
         const word nIon(ions[k]);
-        const word nNeu(neutrals[k]);
 
         const volScalarField& Di = composition.D(nIon);
         const scalar z(composition.z(composition.species()[nIon]));
@@ -29,13 +28,9 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         // P_Reflex list for the ion
         const scalarList &rK = reflexes[k];
 
-        surfaceScalarField::GeometricBoundaryField &bfIonFlux
-            = ionFluxBFs[k];
-        surfaceScalarField::GeometricBoundaryField &bfNeuFlux
-            = neutralFluxBFs[k];
-
+        surfaceScalarField::GeometricBoundaryField&
+            bfIonFlux = ionFluxBFs[k];
         bfIonFlux = phi.boundaryField();
-        bfNeuFlux = phi.boundaryField();
 
         // Adding drift flux to boundary patches
         forAll (bfIonFlux, pidx)
@@ -50,40 +45,55 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         }
 
         const scalar WIon(composition.W(composition.species()[nIon]));
-        const scalar WNeu(composition.W(composition.species()[nNeu]));
-
         const scalar MIon(WIon / NA.value() / 1000.0);
-        const scalar MNeu(WNeu / NA.value() / 1000.0);
 
         const volScalarField& Yion = composition.Y(nIon);
-        const volScalarField& Yneu = composition.Y(nNeu);
 
-        forAll(wallPatcheIDs, pidx) // loop over wall patches
+        const hashedWordList &targets = targetList[k];
+
+        forAll (targets, tidx)
         {
-            label patchID = wallPatcheIDs[pidx];
+            const word &nNeu(targets[tidx]);
 
-            // Probability of ion reflex
-            const scalar pReflex = max(min(rK[pidx],1.0),0.0);
+            surfaceScalarField::GeometricBoundaryField &
+                bfNeuFlux = neutralFluxBFs[neutrals[nNeu]];
+            bfNeuFlux = phi.boundaryField();
 
-            scalarField &wallFluxIon = bfIonFlux[patchID];
-            scalarField &wallFluxNeu = bfNeuFlux[patchID];
+            const scalar WNeu(composition.W(composition.species()[nNeu]));
+            const scalar MNeu(WNeu / NA.value() / 1000.0);
 
-            const scalarField &wallMSf  =  msf.boundaryField()[patchID];
-            const scalarField &wallT    =    T.boundaryField()[patchID];
-            const scalarField &wallYion = Yion.boundaryField()[patchID];
-            const scalarField &wallYneu = Yneu.boundaryField()[patchID];
+            const volScalarField& Yneu = composition.Y(nNeu);
 
-            scalarField vt(sqrt(8.0*kB.value()/pi/MIon*wallT) / 4.0);
+            forAll(wallPatcheIDs, pidx) // loop over wall patches
+            {
+                const label patchID = wallPatcheIDs[pidx];
 
-            // remove negative wallFlux value (flux from wall)
-            wallFluxIon = max(wallFluxIon, 0.0);
+                // Probability of ion reflex
+                const scalar pReflex = max(min(rK[pidx],1.0),0.0);
 
-            // add flux by thermal velocity
-            wallFluxIon += vt * wallMSf;
-            wallFluxIon *= (1.0 - pReflex);
+                scalarField &wallFluxIon = bfIonFlux[patchID];
 
-            // add flux by ion neutralization
-            wallFluxNeu -= wallFluxIon * wallYion / wallYneu / (WIon / WNeu);
+                const scalarField &wallMSf  =  msf.boundaryField()[patchID];
+                const scalarField &wallT    =    T.boundaryField()[patchID];
+                const scalarField &wallYion = Yion.boundaryField()[patchID];
+
+                const scalarField vt(sqrt(8.0*kB.value()/pi/MIon*wallT) / 4.0);
+
+                // remove negative wallFlux value (flux from wall)
+                wallFluxIon = max(wallFluxIon, 0.0);
+
+                // add flux by thermal velocity
+                wallFluxIon += vt * wallMSf;
+                wallFluxIon *= (1.0 - pReflex);
+
+                // add flux by ion neutralization
+                scalarField &wallFluxNeu = bfNeuFlux[patchID];
+
+                const scalarField &wallYneu = Yneu.boundaryField()[patchID];
+
+                wallFluxNeu -= wallFluxIon * wallYion / wallYneu
+                                              / (WIon / WNeu);
+            }
         }
     }
 
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index e7ec79a8..439e1c34 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -425,25 +425,73 @@ dictionary wallIonFluxes
     physicalProperties.subDict("wallIonFluxes")
 );
 
+const label nMaxTargets = 5;
 const hashedWordList ions(wordList(wallIonFluxes.lookup("ions")));
-wordList neutrals_ (ions.size(), "");
+wordList neutrals_ (ions.size()*nMaxTargets);
+PtrList<hashedWordList> targetList (ions.size());
 PtrList<scalarList> reflexes (ions.size());
 PtrList<surfaceScalarField::GeometricBoundaryField> ionFluxBFs (ions.size());
-PtrList<surfaceScalarField::GeometricBoundaryField> neutralFluxBFs (ions.size());
+
+label nTargets = 0;
 
 forAll (ions, iidx)
 {
     const dictionary &wallIonFlux = wallIonFluxes.subDict(ions[iidx]);
 
-    neutrals_[iidx] = word(wallIonFlux.lookup("neutral"));
-    reflexes.set(iidx,
-            new scalarList(wallIonFlux.lookup("wallReflexes")));
+    const Switch noTarget (wallIonFlux.lookupOrDefault("noTarget", false));
+
+    const wordList targets (wallIonFlux.lookupOrDefault("neutrals", wordList()));
+
+    if (noTarget)
+    {
+        targetList.set(iidx, new hashedWordList());
+    }
+    else if (targets.empty())
+    {
+        word neutralName(wallIonFlux.lookup("neutral"));
+        neutrals_[nTargets] = neutralName;
+        nTargets += 1;
+
+        targetList.set(iidx, new hashedWordList(wordList(1, neutralName)));
+    }
+    else
+    {
+        wordList::subList sub (neutrals_, targets.size(), nTargets);
+        forAll (targets, tidx)
+        {
+            sub[tidx] = targets[tidx];
+        }
+        nTargets += targets.size();;
+
+        targetList.set(iidx, new hashedWordList(targets));
+    }
+
+    reflexes.set(iidx, new scalarList(wallIonFlux.lookup("wallReflexes")));
+
     ionFluxBFs.set(iidx,
-            new surfaceScalarField::GeometricBoundaryField
-            (phi.boundaryField()));
-    neutralFluxBFs.set(iidx,
-            new surfaceScalarField::GeometricBoundaryField
-            (phi.boundaryField()));
+            new surfaceScalarField::GeometricBoundaryField(phi.boundaryField()));
 }
 
-const hashedWordList neutrals(neutrals_);
+hashedWordList neutrals;
+
+label nNeutrals = 0;
+
+for (label nidx = 0; nidx < nTargets; nidx++)
+{
+    const word n(neutrals_[nidx]);
+
+    if (!neutrals.contains(n))
+    {
+        neutrals.append(n);
+        nNeutrals += 1;
+    }
+
+}
+
+PtrList<surfaceScalarField::GeometricBoundaryField> neutralFluxBFs (neutrals.size());
+
+forAll (neutrals, iidx)
+{
+    neutralFluxBFs.set(iidx,
+            new surfaceScalarField::GeometricBoundaryField (phi.boundaryField()));
+}