Numerical Computation of Kinetic Parameters for
Reductive Dechlorination of Chlorinated Solvents
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This problem was presented in the paper by Antonia Balazova, Marian Slodicka and Roger Van Keer, Parameter determination for reductive dechlorination of chlorinated solvents, accepted for publication in Transport in Porous Media, 2006.

The pathway of degradation of PCE can be expressed by the following scheme :

PCE -> TCE -> DCE ->  VC ->  ETH

where PCE is tetrachlorethylene and TCE is trichlorethylene, 
DCE is dichloroethene, VC is vinyl chloride and finally ETH is ethylene.

Governing equations are given by the Michaelis-Menten’s model :

C1'[t]=-((k1 C1[t])/(K1 + C1[t]))
C2'[t]=\(\(k1\ C1[t]\)\/\(K1 + C1[t]\) - \(k2\ C2[t]\)\/\(K2 + C2[t]\)\)
C3'[t]=\(\(k2\ C2[t]\)\/\(K2 + C2[t]\) - \(k3\ C3[t]\)\/\(K3 + C3[t]\)\)
C4'[t]=\(\(k3\ C3[t]\)\/\(K3 + C3[t]\) - \(k4\ C4[t]\)\/\(K4 + C4[t]\)\)
C5'[t]=\(\(k4\ C4[t]\)\/\(K4 + C4[t]\) - \(k5\ C5[t]\)\/\(K5 + C5[t]\)\)


where C1, C2, C3, C4 and C5 are the molar concentrations of parent (PCE),
intermediate compounds (TCE,DCE and VC) and the final product (ETH), respectively.

Off[General::"spell1"]

Off[NDSolve::"ndsz"]

Off[InterpolatingFunction::"dmval"]

Off[NDSolve::"ndnum"]

Off[FindMinimum::"lstol"]

Solution of the governing equations with K1=2; K2=5; K3=4; K4=6; K5=10; k1=1; k2=1; k3=2; k4=1 and k5=0.5.

K1 = 2; K2 = 5; K3 = 4; K4 = 6; K5 = 10;

k1 = 1; k2 = 1; k3 = 2; k4 = 1; k5 = 0.5;

sol = 
      NDSolve[{D[
              C1[t], {t, 
                1}] ==  -k1 C1[
                  t]/(K1 + C1[t]), 
D[
              C2[t], {t, 1}] ==  -k2 C2[t]/(K2 + C2[t]) + 
              k1 C1[t]/(K1 + C1[t]), 
D[
              C3[t], {t, 1}] ==  -k3 C3[t]/(K3 + C3[t]) + 
              k2 C2[t]/(K2 + C2[t]), 
D[
              C4[t], {t, 1}] ==  -k4 C4[t]/(K4 + C4[t]) + 
              k3 C3[t]/(K3 + C3[t]), 
D[
              C5[t], {t, 1}] ==  -k5 C5[t]/(K5 + C5[t]) + 
              k4 C4[t]/(K4 + C4[t]), C1[0] ==  30, 
          C2[0] ==  0, C3[0] ==  0, C4[0] ==  0, 
          C5[0] ==  0}, 
{C1[t], C2[t], C3[t], 
          C4[t], C5[t]}, {t, 0, 100}]

plt1 = 
      Plot[Evaluate[C1[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[1/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt2 = 
      Plot[Evaluate[C2[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[2/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt3 = 
      Plot[Evaluate[C3[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[3/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt4 = 
      Plot[Evaluate[C4[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[4/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt5 = 
      Plot[Evaluate[C5[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[5/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

Show[plt1, plt2, plt3, plt4, plt5, 
      DisplayFunction ->  $DisplayFunction]

Determination of kinetic parameters K1, K2, K3, K4, K5, k1, k2, k3, k4 and k5 using data without noise and data with 10% noise. Data is obtained from the above solution. 10% Noise represents experimental errors.

tbl1 = 
      Table[First[Evaluate[C1[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl2 = 
      Table[First[Evaluate[C2[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl3 = 
      Table[First[Evaluate[C3[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl4 = 
      Table[First[Evaluate[C4[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl5 = 
      Table[First[Evaluate[C5[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl1n = 
      Table[First[Evaluate[C1[t] /. sol] /. t -> 10 i] + 
          0.1 (1 - 2 Random[]) 
            First[Evaluate[C1[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl2n = 
      Table[First[Evaluate[C2[t] /. sol] /. t -> 10 i] + 
          0.1 (1 - 2 Random[]) 
            First[Evaluate[C2[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl3n = 
      Table[First[Evaluate[C3[t] /. sol] /. t -> 10 i] + 
          0.1 (1 - 2 Random[]) 
            First[Evaluate[C3[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl4n = 
      Table[First[Evaluate[C4[t] /. sol] /. t -> 10 i] + 
          0.1 (1 - 2 Random[]) 
            First[Evaluate[C4[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

tbl5n = 
      Table[First[Evaluate[C5[t] /. sol] /. t -> 10 i] + 
          0.1 (1 - 2 Random[]) 
            First[Evaluate[C5[t] /. sol] /. t -> 10 i], {i, 0, 
          10}]

K1 =. ; K2 =. ; K3 =. ; K4 =. ; K5 =. ; k1 =. ; k2 =. ; k3 =. ; k4 =. ; 
    k5 =. ;

Sol[K1_, K2_, K3_, K4_, K5_, k1_, k2_, k3_, k4_, k5_] := 
      NDSolve[{D[
              C1[t], {t, 
                1}] ==  -k1 C1[
                  t]/(K1 + C1[t]), 
D[
              C2[t], {t, 1}] ==  -k2 C2[t]/(K2 + C2[t]) + 
              k1 C1[t]/(K1 + C1[t]), 
D[
              C3[t], {t, 1}] ==  -k3 C3[t]/(K3 + C3[t]) + 
              k2 C2[t]/(K2 + C2[t]), 
D[
              C4[t], {t, 1}] ==  -k4 C4[t]/(K4 + C4[t]) + 
              k3 C3[t]/(K3 + C3[t]), 
D[
              C5[t], {t, 1}] ==  -k5 C5[t]/(K5 + C5[t]) + 
              k4 C4[t]/(K4 + C4[t]), C1[0] ==  30, 
          C2[0] ==  0, C3[0] ==  0, C4[0] ==  0, 
          C5[0] ==  0}, 
{C1[t], C2[t], C3[t], 
          C4[t], C5[t]}, {t, 0, 100}]

Determination of kinetic parameters K1, K2, K3, K4, K5, k1, k2, k3, k4 and k5 using data without noise. We find K1=2; K2=5; K3=4; K4=6; K5=10; k1=1; k2=1; k3=2; k4=1 and k5=0.5. Computational time is 15 min using a Pentium III 800 MHZ.

Obj[K1_?NumericQ, K2_?NumericQ, K3_?NumericQ, K4_?NumericQ, 
        K5_?NumericQ, k1_?NumericQ, k2_?NumericQ, k3_?NumericQ, k4_?NumericQ, 
        k5_?NumericQ] := 
      Sum[(tbl1[[i + 1]] - 
                First[Evaluate[
                      C1[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl2[[i + 1]] - 
                First[Evaluate[
                      C2[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl3[[i + 1]] - 
                First[Evaluate[
                      C3[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl4[[i + 1]] - 
                First[Evaluate[
                      C4[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl5[[i + 1]] - 
                First[Evaluate[
                      C5[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}]

FindMinimum[
      Obj[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5], {K1, 2.2, 1.8, 2.5}, {K2, 
        5.2, 4.8, 5.5}, {K3, 4.2, 3.8, 4.5}, {K4, 6.2, 5.8, 6.5}, {K5, 10.2, 
        9.8, 10.5}, {k1, 1.2, 0.8, 1.5}, {k2, 1.2, 0.8, 1.5}, {k3, 2.2, 1.8, 
        2.5}, {k4, 1.2, 0.8, 1.5}, {k5, 0.7, 0.3, 0.8}]

Determination of kinetic parameters K1, K2, K3, K4, K5, k1, k2, k3, k4 and k5 using data with 10% noise. We find K1->1.93686, K2->3.47416, K3->5.24849, K4->9.59076,
K5->10.5118, k1->1.02308, k2->0.86059, k3->2.42092, k4->1.3047 and k5->0.495326. Computational time is 15 min using a Pentium III 800 MHZ. Values of parameters will be different for different runs because random data used to add 10% noise will be different from run to run. User may need to change the value of the initial guess of FindMinimum to get correct solution.

Obj[K1_?NumericQ, K2_?NumericQ, K3_?NumericQ, K4_?NumericQ, 
        K5_?NumericQ, k1_?NumericQ, k2_?NumericQ, k3_?NumericQ, k4_?NumericQ, 
        k5_?NumericQ] := 
      Sum[(tbl1n[[i + 1]] - 
                First[Evaluate[
                      C1[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl2n[[i + 1]] - 
                First[Evaluate[
                      C2[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl3n[[i + 1]] - 
                First[Evaluate[
                      C3[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl4n[[i + 1]] - 
                First[Evaluate[
                      C4[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}] + 
        Sum[(tbl5n[[i + 1]] - 
                First[Evaluate[
                      C5[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> 10 i])^2, {i, 0, 10}]

SOL = 
      FindMinimum[
          Obj[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5], {K1, 2.2}, {K2, 
            5.2}, {K3, 4.2}, {K4, 6.2}, {K5, 10.2}, {k1, 1.2}, {k2, 
            1.2}, {k3, 2.2}, {k4, 1.2}, {k5, 0.7}, 
          MaxIterations ->  1000] // Timing

sol = 
      Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5] /. 
        SOL[[2, 2]]

plt1 = 
      Plot[Evaluate[C1[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[1/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt2 = 
      Plot[Evaluate[C2[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[2/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt3 = 
      Plot[Evaluate[C3[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[3/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt4 = 
      Plot[Evaluate[C4[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[4/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt5 = 
      Plot[Evaluate[C5[t] /. sol], {t, 0, 100}, PlotStyle -> Hue[5/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

Show[plt1, plt2, plt3, plt4, plt5, 
      DisplayFunction ->  $DisplayFunction]

Tbl1n = Transpose[{Table[i 10, {i, 0, 10}], tbl1n}]

Tbl2n = Transpose[{Table[i 10, {i, 0, 10}], tbl2n}]

Tbl3n = Transpose[{Table[i 10, {i, 0, 10}], tbl3n}]

Tbl4n = Transpose[{Table[i 10, {i, 0, 10}], tbl4n}]

Tbl5n = Transpose[{Table[i 10, {i, 0, 10}], tbl5n}]

plt1n = 
      ListPlot[Tbl1n, PlotStyle -> {Hue[1/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt2n = 
      ListPlot[Tbl2n, PlotStyle -> {Hue[2/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt3n = 
      ListPlot[Tbl3n, PlotStyle -> {Hue[3/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt4n = 
      ListPlot[Tbl4n, PlotStyle -> {Hue[4/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt5n = 
      ListPlot[Tbl5n, PlotStyle -> {Hue[5/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

Obtained  kinetics parameters give concentration profiles that are in agreement with the noisy data used for their determination.

Show[plt1, plt1n, plt2, plt2n, plt3, plt3n, plt4, plt4n, plt5, plt5n, 
      DisplayFunction ->  $DisplayFunction]

Determination of kinetic parameters K1, K2, K3, K4, K5, k1, k2, k3, k4 and k5 using experimental data from W.M. Wu, J. Nye, R.F. Hickey, M.K. Jain, and J.G. Zeikus. Dechlorination of PCE and TCE to ethylene using an anaerobic microbial consortium. In R.E. Hinchee, A. Leeson, and L. Semprini, editors, Bioremediation of chlorinated solvents, pages 45-52. Batelle Pr. Columbus, 1995.

time={0, 2.5, 8, 13, 30, 40, 54, 77, 90}*50/90.

pce = {60, 50, 30, 21, 0, 0, 0, 0, 0}/40.

tce = {0, 2, 11, 20, 26, 16, 8, 2, 0}/40.

dce = {0, 1, 2, 4, 8, 8, 4, 2, 0}/40.

vc = {0, 0, 1, 2, 12, 15, 15, 4, 0}/40.

eth = {0, 0, 1, 3, 11, 18, 33, 55, 59}/40.

Sol[K1_, K2_, K3_, K4_, K5_, k1_, k2_, k3_, k4_, k5_] := 
      NDSolve[{D[
              C1[t], {t, 
                1}] ==  -k1 C1[
                  t]/(K1 + C1[t]), 
D[
              C2[t], {t, 1}] ==  -k2 C2[t]/(K2 + C2[t]) + 
              k1 C1[t]/(K1 + C1[t]), 
D[
              C3[t], {t, 1}] ==  -k3 C3[t]/(K3 + C3[t]) + 
              k2 C2[t]/(K2 + C2[t]), 
D[
              C4[t], {t, 1}] ==  -k4 C4[t]/(K4 + C4[t]) + 
              k3 C3[t]/(K3 + C3[t]), 
D[
              C5[t], {t, 1}] ==  -k5 C5[t]/(K5 + C5[t]) + 
              k4 C4[t]/(K4 + C4[t]), C1[0] ==  1.5, 
          C2[0] ==  0, C3[0] ==  0, C4[0] ==  0, 
          C5[0] ==  0}, 
{C1[t], C2[t], C3[t], 
          C4[t], C5[t]}, {t, 0, 50}]

Obj[K1_?NumericQ, K2_?NumericQ, K3_?NumericQ, K4_?NumericQ, 
        K5_?NumericQ, k1_?NumericQ, k2_?NumericQ, k3_?NumericQ, k4_?NumericQ, 
        k5_?NumericQ] := 
      Sum[(pce[[i + 1]] - 
                First[Evaluate[
                      C1[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> time[[i + 1]]])^2, {i, 0, 8}] + 
        Sum[(tce[[i + 1]] - 
                First[Evaluate[
                      C2[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> time[[i + 1]]])^2, {i, 0, 8}] + 
        Sum[(dce[[i + 1]] - 
                First[Evaluate[
                      C3[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> time[[i + 1]]])^2, {i, 0, 8}] + 
        Sum[(vc[[i + 1]] - 
                First[Evaluate[
                      C4[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> time[[i + 1]]])^2, {i, 0, 8}] + 
        Sum[(eth[[i + 1]] - 
                First[Evaluate[
                      C5[t] /. Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5]] /. 
                    t -> time[[i + 1]]])^2, {i, 0, 8}]

SOLEXP = 
      FindMinimum[
          Obj[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5], {K1, 0.02}, {K2, 
            0.15}, {K3, 0.04}, {K4, 0.2}, {K5, 0.2}, {k1, 0.1}, {k2, 
            0.05}, {k3, 0.05}, {k4, 0.05}, {k5, 10^-4}, 
          MaxIterations ->  1000] // Timing

sol = 
      Sol[K1, K2, K3, K4, K5, k1, k2, k3, k4, k5] /. 
        SOLEXP[[2, 2]]

plt1 = 
      Plot[Evaluate[C1[t] /. sol], {t, 0, 50}, PlotStyle -> Hue[1/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt2 = 
      Plot[Evaluate[C2[t] /. sol], {t, 0, 50}, PlotStyle -> Hue[2/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt3 = 
      Plot[Evaluate[C3[t] /. sol], {t, 0, 50}, PlotStyle -> Hue[3/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt4 = 
      Plot[Evaluate[C4[t] /. sol], {t, 0, 50}, PlotStyle -> Hue[4/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

plt5 = 
      Plot[Evaluate[C5[t] /. sol], {t, 0, 50}, PlotStyle -> Hue[5/6], 
        PlotRange -> All, DisplayFunction ->  Identity]

Show[plt1, plt2, plt3, plt4, plt5, 
      DisplayFunction ->  $DisplayFunction]

Tbl1n = Transpose[{Table[time[[i]], {i, 1, 9}], pce}]

Tbl2n = Transpose[{Table[time[[i]], {i, 1, 9}], tce}]

Tbl3n = Transpose[{Table[time[[i]], {i, 1, 9}], dce}]

Tbl4n =  Transpose[{Table[time[[i]], {i, 1, 9}], vc}]

Tbl5n = Transpose[{Table[time[[i]], {i, 1, 9}], eth}]

plt1n = 
      ListPlot[Tbl1n, PlotStyle -> {Hue[1/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt2n = 
      ListPlot[Tbl2n, PlotStyle -> {Hue[2/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt3n = 
      ListPlot[Tbl3n, PlotStyle -> {Hue[3/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt4n = 
      ListPlot[Tbl4n, PlotStyle -> {Hue[4/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

plt5n = 
      ListPlot[Tbl5n, PlotStyle -> {Hue[5/6], PointSize[0.025]}, 
        DisplayFunction ->  Identity]

Obtained  kinetics parameters give concentration profiles that are in agreement with the experimental data used for their determination.

Show[plt1, plt1n, plt2, plt2n, plt3, plt3n, plt4, plt4n, plt5, plt5n, 
      DisplayFunction ->  $DisplayFunction]