# Modeling Heavy Metal Sorption Kinetics Using Fractional Calculus

Authors: V. C. Friesen; D. P. Leitoles; G. Gonçalves; E. K. Lenzi; M. K. Lenzi
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Heavy metals are commonly regarded as environmentally aggressive and hazardous to human health. Among the different metals, lead plays an important economic role due to its large use in the automotive industry, being an essential component of batteries. Different approaches have been reported in the literature aimed at lead removal, and among them a very successful one considers the use of water hyacinths for sorption-based operation. The modeling of the metal sorption kinetics is a fundamental step towards in-depth studies and proper separation equipment design and optimization. Fractional calculus represents a novel approach and a growing research field for process modeling, which is based on the successful use of derivatives of arbitrary order. This paper reports the modeling of the kinetics of lead sorption by water hyacinths (Eichhornia crassipes) using a fractional calculus. A general procedure on error analysis is also employed to prove the actual fractional nature of the proposed model by the use of parametric variance analysis, which was carried out using two different approaches (with the complete Hessian matrix and with a simplified Hessian matrix). The joint parameter confidence regions were generated, allowing to successfully show the fractional nature of the model and the sorption process.


---

## Body

## 1. Introduction

Heavy metals, particularly lead, are widely present in industrial applications such as in paints [1] and in batteries [2]. The presence of lead represents an important environmental issue due to its high toxicity [3, 4]. Therefore, several techniques have been successfully proposed to remove lead and other heavy metals from wastewater streams allowing a safe discharge. The most common techniques are based on chemical precipitation [5], ion exchange [6], electrochemical methods [2], and sorption/biosorption [7]. Over the past decade, several applications have reported the use of water hyacinths [7, 8] aimed at lead removal, providing a self-sustainable approach.An essential tool for in-depth process studies, design, optimization, and control of processes concerning lead removal is given by mathematical models. They play a very important role as actual experimental runs cannot necessarily be carried out because accurate predictions and deviations can be obtained from simulation of different operating conditions [9]. This turns out to be a very efficient and economic approach for process studies, if the model has been properly derived and tuned by adequate parameter estimation and error analysis tasks [10].Modeling tasks in heavy metal sorption are usually splitted into two parts: (i) equilibrium/isotherm modeling and (ii) kinetics modeling. Equilibrium/isotherm modeling usually concerns the use of classical models based on Langmuir [11], Freundlich [12], and Redlich and Peterson [13], among others [14]. On the other hand, kinetics modeling usually employs either first and pseudo-first order models [15] or second and pseudo-second order models [16].Fractional calculus represents an innovative approach and a very fast growing research field for process modeling, based on derivatives of arbitrary order [17]. The literature reports a broad range of applications ranging from molecular diffusion processes [18] to industrial systems engineering [19]. Regarding heavy metals sorption modeling with the aid of fractional calculus, Gomes et al. [10] reported a novel isotherm to describe the lead sorption behavior and also provide an in-depth methodology to analyze the actual fractional feature of the estimated order of the fractional derivative. Concerning fractional kinetics modeling, some milestone reported works provide the fundamental knowledge of the theoretical fractional kinetics [20, 21], while others report the application of the fractional modeling approach to different heavy metal sorption systems [22–25].This work reports the modeling of the kinetic behavior of experimental lead sorption data usingEichhornia crassipes as the sorption medium. In order to describe the experimental behavior, a fractional kinetics model is considered and a comparison to integer order kinetic model is also performed. Differently from previously reported manuscripts, a parametric analysis [10] is carried out and joint parameter confidence regions are generated, being successfully able to show that the fractional nature of the model is statistically valid.


## 2. Experimental Data

The experimental data used for the modeling purposes were reported by dos Santos and Lenzi [8], where further details regarding the experimental procedures can be obtained. Table 1 presents a summary of the experimental datasets used. In both sets, the concentration of the water hyacinth was the same and equal to 20 g of hyacinth per liter of solution. This amount of hyacinth was submitted to two different initial concentrations of lead. In experimental run RUN01, this concentration was equal to 10 mg of lead (Pb)/L, while in experimental run RUN02, the initial concentration was equal to 15 mg of lead (Pb)/L. Prior to the parameter estimation tasks, the data were normalized by dividing all values by the initial condition (concentration). Finally it is worth mentioning that the concentration values present in Table 1 describe the dynamic behavior of the decrease of the concentration of lead in the solution in contact with the hyacinth.Table 1
Experimental data for the lead sorption.
RUN01 –y0 = 10 mgPb/L
y
(
t
)
/
y
0
1
0.479
0.292
0.214
0.142
0.075
0.049
0.049
t – [day]
0
0.125
0.25
0.375
0.5
1
2
3
RUN02 –y0 = 15 mgPb/L
y
(
t
)
/
y
0
1
0.519
0.389
0.217
0.161
0.088
0.041
0.041
t – [day]
0
0.125
0.25
0.375
0.5
1
2
3


## 3. Theoretical Framework

As already mentioned, the goal of this work is the description of the dynamic behavior of lead sorption. In order to accomplish this task, the fractional kinetic model given by (1) a long with a correspondent initial condition is considered. This model resembles the classical one of first order used to describe processes with exponential relaxation; however, the differential operator is extended to a noninteger order, α, with the aid of the fractional calculus generalization in order to cover a broad of situations which are not typically characterized by exponential processes. The choice of this model occurred because of the experimental data behavior the drop of concentration of lead in water resembles an exponential one, but with different asymptotic behavior:(1)dαytdtα=-k·yt,yt=0=y0.Although there is a fractional order derivate, the model equation given by (1) is linear; therefore, the analytical solution can be obtained by Laplace transform technique [9]. By considering the Caputo approach [26] for the fractional derivate, (2) can be obtained after the Laplace transform inversion. This solution resembles the Mittag-Leffler function [27] with two parameters [28]:(2)yt=y0·∑j=0∞-kj·tα·j+1-1Γα·j+1.For the sake of comparison, the classical first order model, given by (3) has also been taken into account. It is important to stress that (3) can be obtained by setting α equal to 1 in (2). This happens due to the generalization capability provided by the Mittag-Leffler based function present in the solution of (1):(3)dytdt=-k·yt,yt=0=y0⟶yt=y0·e-k·t.Depending on the value of the orderα of the derivative, the Mittag-Leffler function kernel present in (2) turns into a different mathematical function, as mentioned; if α=1, an exponential function shows up. Therefore, one must pay attention that, as mentioned above, α is not an ordinary parameter, as it not only provides a better the data fitting but it also turns the sum of (2) into the best mathematical function to describe the experimental behavior.The statistical validation methodology steps used were reported by Gomes et al. [10], where further details can be obtained. The parameter estimation task involved the determination of the parameters α and k of the fractional model, (2), and the parameter k of the integer order model, (3). This task considered an objective function given by (4). A genetic algorithm was used to avoid a local minimum value of the objective function:(4)FOBJ=∑p=1NR∑i=1NEyp,iEt-yp,iMt2.The experimental variance matrix,VyE, is given by (5), where the hypothesis of constant variance, δyE2, was considered valid. Besides, experimental covariance was not taken into account:(5)VyENE·NR×NE·NR=δyE2·10⋯001⋯0⋮⋮⋱⋮00⋯1;δyE2=∑p=1NR∑i=1NEyp,iEt-yp,iMt2NR·NE-NP.In this work, the parametric variance analysis considered two different approaches, with the complete Hessian matrixHθ and with a simplified Hessian matrix [B], given by (6) and (7), respectively. The elements of the matrix Hθ are given by (8), while the elements of [B] consider only the first order derivates of (8). In the following expressions, θ refers to the parameter set {α;k} to be estimated. Further details can be found in the literature [10]: (6)VθNP×NP=Hθ-1NP×NP·GyNP×NE·NR·VyENE·NR×NE·NRVyENE·NR×NE·NR·GyTNE·NR×NP·Hθ-1NP×NP,(7)VθNP×NP=BTNP×NE·NR·VyE-1NE·NR×NE·NRVyE-1NE·NR×NE·NR·BNE·NR×NP-1,(8)hij=∂2FOBJ∂θi∂θJ=2·∂yM∂θiT·VyE-1·∂yM∂θj-2·∂2yM∂θi∂θjT·Vy-1·yE-yM.Finally, joint parameter confidence regions were generated. As reported by Gomes et al. [10], (9) refers to joint confidence regions obtained by the linearization of the objective function, which are typically of ellipsoidal shape. On the other hand, (10) presents the formula used to obtain joint confidence region considering the nonlinear features of the objective function. This region is not necessarily ellipsoidal and it is usually much larger than the one obtained by (9). In these expressions, θ is the set of parameters, θ∗ is the set of estimated parameters, FNP,NE-NP95% is the value of the Fischer distribution with 95% of confidence and NP and (NE-NP) are degrees of freedom, NE is the number of experimental points of a given experimental run, and NR is the number of experimental runs:(9)θ-θ∗T1×NP·Vθ-1NP×NP·θ-θ∗NP×1≤FOBJθ∗·NPNE-NP·FNP,NE-NP95%,(10)FOBJθ≤FOBJθ∗·1+NPNE-NP·FNP,NE-NP95%.The modeling work was carried out considering two different approaches. In the first one, named APPROACH01, the experimental set RUN01 was used for parameter estimation of each model, that is, (2) and (3), and the resulting equations were used to validate the behavior exhibited by the experimental data set RUN 02. In the second one, APPROACH02, both experimental sets were used to estimate a global set of parameters for each model. It must be mentioned that, for each of these approaches, the parametric variance was calculated using the two different approaches mentioned in the last paragraph, that is, (6) and (7).


## 4. Results and Discussions

### 4.1. APPROACH01 Results

Table2 presents a summary of the results of the parameter estimation. FRAC01 refers to the parametric analysis of the fractional model, (2), using the complete hessian matrix (6) to evaluate the variance of the estimated parameters. FRAC02 refers to the parametric analysis of the fractional model, (2), using the simplified hessian matrix (7) to evaluate the variance of the estimated parameters. INT01 refers to the parametric analysis of the integer order model, (3), using the complete hessian matrix (6) to evaluate the variance of the estimated parameters. INT02 refers to the parametric analysis of the integer order model, (3), using the simplified hessian matrix (7) to evaluate the variance of the estimated parameters. From the results reported in Table 1, it is possible to observe that the choice of calculation of parametric variance has a little influence on the parameter deviation. This probably happened because as the data fit was good, the term yE-yM, corresponding to the difference between experimental results and model predictions, has a small value.Table 2
Parameters and confidence interval for APPROACH01.
Model
Parameter
Parameter standard deviationδθ
Parameter confidence intervalθi±z95%·δθi
r
F
OBJ
RUN01
FRAC01
k: 0.728
0.239
0.73 ± 0.47
0.99
3.64·10−3
α: 0.203
0.027
0.20 ± 0.05
FRAC02
k: 0.728
0.221
0.73 ± 0.43
0.99
3.64·10−3
α: 0.203
0.025
0.20 ± 0.05
INT01
k: 4.775
0.498
4.8 ± 1.0
0.98
1.91·10−2
INT02
k: 4.775
0.438
4.8 ± 0.8
0.98
1.91·10−2The objective function (FOBJ) has a low value for fractional model when compared to the integer order model; however, for the fractional model, the value is one order of magnitude lower. This probably happened because the fractional model could better describe the memory effects [10, 21] present in the sorption process. Also, one can note that the parameter α, when estimated, also shapes the sum to a specific mathematical function, providing a better fit than the integer order model. It is also important to stress that due to its intrinsic features, the fractional order model also considers memory effects, which happens due to the convolution integral of the fractional derivate approach considered. This type of effect is constantly present in irregular, disordered, and fractal media diffusion and mass transfer [29]. Taking the memory effects into account significantly improves the modeling results, such as in systems like the irregular shape of the hyacinth leaves present in the interface with the lead aqueous solution.It is worth mentioning that the correlation coefficient (r) was very good in both cases, presenting values equal to 0.99. Moreover, the estimated values of parameters k and α are greater than their respective parametric deviation for all models; consequently these parameters can be regarded as statistically significant.Figures1 and 2 show a comparison between experimental data and model predictions plotted against time for model FRAC01 and INT01, where parameter variance was obtained using the complete Hessian matrix. Again, the fractional model described the experimental data more adequately than the integer order model, as expected after analyzing Table 1.Figure 1
Experimental and predicted values against time Experiment  01, fractional order model, APPROACH01.Figure 2
Experimental and predicted values against time Experiment  01, integer order model, APPROACH01.In these figures, one can see the confidence interval of the model predictions. It can be observed that for both models, the experimental data lie within these regions, but in the fractional model analysis the obtained regions are smaller showing that the model predictions tend to be more accurate.Figures3 and 4 present the prediction/validation of experimental set RUN02 using the model parameters estimated using only experimental RUN01. It can be clearly seen that the model predictions of the fractional model are better than the predictions of the integer order model, as the experimental points are closer to the model prediction and the confidence region is smaller. Finally, it is worth mentioning that this result shows the extrapolation capability of operating conditions as the amount of lead in experiment RUN02 is 50% higher than in experiment RUN01.Figure 3
Experimental and predicted values against time Experiment  02, fractional order model, APPROACH01.Figure 4
Experimental and predicted values against time Experiment  02, integer order model, APPROACH01.Figure5 shows the joint parameter confidence (ellipsoid) region, considering a confidence level of 95% and the linearization of the estimation problem as mentioned before (see (9)). Figure 6 shows the joint parameter confidence region obtained considering the nonlinear feature of the estimation problem; see (10). In both cases, it can be seen that parameter α does not reach a value equal to 1, proving that the fractional order estimated is statistically not integer, allowing one to conclude that the model is in fact fractional.Figure 5
Confidence region (ellipsoid) for simplified and complete Hessian matrix, APPROACH01.Figure 6
Confidence region (nonlinear), APPROACH01.


### 4.2. APPROACH02 Results

Table3 presents a summary of the results of the parameter estimation. FRAC01 refers to calculation using complete hessian matrix (6) of the fractional model and FRAC02 refers to the simplified hessian matrix (7). The same observation is valid for the integer order model (INT01 and INT02). From this estimation task result, it is possible to observe that the choice of calculation of parametric variance has practically no influence on the parameter deviation for this case. This probably happened because, as both sets were simultaneously used, the larger number of experimental points led to a lower experimental variance (5).Table 3
Parameters and confidence interval for APPROACH02.
Model
Parameter
Parameter standard deviationδθ
Parameter confidence intervalθi±z95%·δθi
r
F
OBJ
RUN01 + RUN02
FRAC01
k: 0.650
0.231
0.65 ± 0.45
0.99
1.78·10−2
α: 0.201
0.029
0.20 ± 0.05
FRAC02
k: 0.650
0.226
0.65 ± 0.44
0.99
1.78·10−2
α: 0.201
0.028
0.20 ± 0.05
INT01
k: 4.412
0.256
4.4 ± 1.0
0.99
4.00·10−2
INT02
k: 4.412
0.201
4.4 ± 0.8
0.99
4.00·10−2The objective function (FOBJ) has a lower value for fractional model when compared to the integer order model. It is worth mentioning that the correlation coefficient (r) was very good in both cases, presenting values over 0.99. Besides, the parameters k, e, and α were greater than their respective parametric deviation in all models; consequently the parameters are statistically significant.Figure7 shows a comparison of experimental data and model predictions plotted against the independent variable for fractional order model. Figure 8 presents the same plot but with considering the integer order model. The deviation was smaller than in previous approach. This is probably because both datasets were used in the estimation problem. However, in both cases it can be observed that the integer order model prediction was slightly worse and has more deviations.Figure 7
Experimental and predicted values against time Experiments  01 and 02, fractional order model, APPROACH02.Figure 8
Experimental and predicted values against time Experiments  01 and 02, INT model, APPROACH02.Figure9 presents the confidence region. As can be seen, independently of the Hessian matrix considered, that is, complete or simplified, the two regions are very close. This happens because the term containing the second order derivate in (8) is closer to zero due to the better fit in APPROACH02 as both datasets were simultaneously used. Finally, Figure 10 shows the joint confidence region considering the nonlinear features of the problem. Again, this region is larger than the one obtained for the linear problem and it also shows that the order of the derivative is smaller than 1, indicating a fractional process.Figure 9
Confidence region with simplified and complete Hessian matrix, APPROACH02.Figure 10
Confidence region (nonlinear), APPROACH02.


## 4.1. APPROACH01 Results

Table2 presents a summary of the results of the parameter estimation. FRAC01 refers to the parametric analysis of the fractional model, (2), using the complete hessian matrix (6) to evaluate the variance of the estimated parameters. FRAC02 refers to the parametric analysis of the fractional model, (2), using the simplified hessian matrix (7) to evaluate the variance of the estimated parameters. INT01 refers to the parametric analysis of the integer order model, (3), using the complete hessian matrix (6) to evaluate the variance of the estimated parameters. INT02 refers to the parametric analysis of the integer order model, (3), using the simplified hessian matrix (7) to evaluate the variance of the estimated parameters. From the results reported in Table 1, it is possible to observe that the choice of calculation of parametric variance has a little influence on the parameter deviation. This probably happened because as the data fit was good, the term yE-yM, corresponding to the difference between experimental results and model predictions, has a small value.Table 2
Parameters and confidence interval for APPROACH01.
Model
Parameter
Parameter standard deviationδθ
Parameter confidence intervalθi±z95%·δθi
r
F
OBJ
RUN01
FRAC01
k: 0.728
0.239
0.73 ± 0.47
0.99
3.64·10−3
α: 0.203
0.027
0.20 ± 0.05
FRAC02
k: 0.728
0.221
0.73 ± 0.43
0.99
3.64·10−3
α: 0.203
0.025
0.20 ± 0.05
INT01
k: 4.775
0.498
4.8 ± 1.0
0.98
1.91·10−2
INT02
k: 4.775
0.438
4.8 ± 0.8
0.98
1.91·10−2The objective function (FOBJ) has a low value for fractional model when compared to the integer order model; however, for the fractional model, the value is one order of magnitude lower. This probably happened because the fractional model could better describe the memory effects [10, 21] present in the sorption process. Also, one can note that the parameter α, when estimated, also shapes the sum to a specific mathematical function, providing a better fit than the integer order model. It is also important to stress that due to its intrinsic features, the fractional order model also considers memory effects, which happens due to the convolution integral of the fractional derivate approach considered. This type of effect is constantly present in irregular, disordered, and fractal media diffusion and mass transfer [29]. Taking the memory effects into account significantly improves the modeling results, such as in systems like the irregular shape of the hyacinth leaves present in the interface with the lead aqueous solution.It is worth mentioning that the correlation coefficient (r) was very good in both cases, presenting values equal to 0.99. Moreover, the estimated values of parameters k and α are greater than their respective parametric deviation for all models; consequently these parameters can be regarded as statistically significant.Figures1 and 2 show a comparison between experimental data and model predictions plotted against time for model FRAC01 and INT01, where parameter variance was obtained using the complete Hessian matrix. Again, the fractional model described the experimental data more adequately than the integer order model, as expected after analyzing Table 1.Figure 1
Experimental and predicted values against time Experiment  01, fractional order model, APPROACH01.Figure 2
Experimental and predicted values against time Experiment  01, integer order model, APPROACH01.In these figures, one can see the confidence interval of the model predictions. It can be observed that for both models, the experimental data lie within these regions, but in the fractional model analysis the obtained regions are smaller showing that the model predictions tend to be more accurate.Figures3 and 4 present the prediction/validation of experimental set RUN02 using the model parameters estimated using only experimental RUN01. It can be clearly seen that the model predictions of the fractional model are better than the predictions of the integer order model, as the experimental points are closer to the model prediction and the confidence region is smaller. Finally, it is worth mentioning that this result shows the extrapolation capability of operating conditions as the amount of lead in experiment RUN02 is 50% higher than in experiment RUN01.Figure 3
Experimental and predicted values against time Experiment  02, fractional order model, APPROACH01.Figure 4
Experimental and predicted values against time Experiment  02, integer order model, APPROACH01.Figure5 shows the joint parameter confidence (ellipsoid) region, considering a confidence level of 95% and the linearization of the estimation problem as mentioned before (see (9)). Figure 6 shows the joint parameter confidence region obtained considering the nonlinear feature of the estimation problem; see (10). In both cases, it can be seen that parameter α does not reach a value equal to 1, proving that the fractional order estimated is statistically not integer, allowing one to conclude that the model is in fact fractional.Figure 5
Confidence region (ellipsoid) for simplified and complete Hessian matrix, APPROACH01.Figure 6
Confidence region (nonlinear), APPROACH01.


## 4.2. APPROACH02 Results

Table3 presents a summary of the results of the parameter estimation. FRAC01 refers to calculation using complete hessian matrix (6) of the fractional model and FRAC02 refers to the simplified hessian matrix (7). The same observation is valid for the integer order model (INT01 and INT02). From this estimation task result, it is possible to observe that the choice of calculation of parametric variance has practically no influence on the parameter deviation for this case. This probably happened because, as both sets were simultaneously used, the larger number of experimental points led to a lower experimental variance (5).Table 3
Parameters and confidence interval for APPROACH02.
Model
Parameter
Parameter standard deviationδθ
Parameter confidence intervalθi±z95%·δθi
r
F
OBJ
RUN01 + RUN02
FRAC01
k: 0.650
0.231
0.65 ± 0.45
0.99
1.78·10−2
α: 0.201
0.029
0.20 ± 0.05
FRAC02
k: 0.650
0.226
0.65 ± 0.44
0.99
1.78·10−2
α: 0.201
0.028
0.20 ± 0.05
INT01
k: 4.412
0.256
4.4 ± 1.0
0.99
4.00·10−2
INT02
k: 4.412
0.201
4.4 ± 0.8
0.99
4.00·10−2The objective function (FOBJ) has a lower value for fractional model when compared to the integer order model. It is worth mentioning that the correlation coefficient (r) was very good in both cases, presenting values over 0.99. Besides, the parameters k, e, and α were greater than their respective parametric deviation in all models; consequently the parameters are statistically significant.Figure7 shows a comparison of experimental data and model predictions plotted against the independent variable for fractional order model. Figure 8 presents the same plot but with considering the integer order model. The deviation was smaller than in previous approach. This is probably because both datasets were used in the estimation problem. However, in both cases it can be observed that the integer order model prediction was slightly worse and has more deviations.Figure 7
Experimental and predicted values against time Experiments  01 and 02, fractional order model, APPROACH02.Figure 8
Experimental and predicted values against time Experiments  01 and 02, INT model, APPROACH02.Figure9 presents the confidence region. As can be seen, independently of the Hessian matrix considered, that is, complete or simplified, the two regions are very close. This happens because the term containing the second order derivate in (8) is closer to zero due to the better fit in APPROACH02 as both datasets were simultaneously used. Finally, Figure 10 shows the joint confidence region considering the nonlinear features of the problem. Again, this region is larger than the one obtained for the linear problem and it also shows that the order of the derivative is smaller than 1, indicating a fractional process.Figure 9
Confidence region with simplified and complete Hessian matrix, APPROACH02.Figure 10
Confidence region (nonlinear), APPROACH02.


## 5. Conclusion

The main purpose of this work was to present a model to describe lead sorption kinetics using water hyacinths taking into account fractional calculus issues and to carry out an error analysis to validate fractional nature of the model. The parameter estimation tasks were made by two different approaches, one of them considering one experimental dataset and the other one for validation, and the second approach with all experimental datasets.Parametric variance matrix was calculated in two different ways, but for this problem the results were roughly the same. From parametric variance, the confidence regions could be obtained considering a linearization of the problem or the nonlinear problem itself.It was concluded that the fractional order model presented better results than the integer order model. Besides, the noninteger order model was in fact fractional, as statistically proven by the error analysis and parameter confidence region. The model in APPROACH01 was satisfactorily used for extrapolation of experimental conditions, condition not reached by integer order model. In this way, the developed model can be applied with success in the heavy metal sorption modeling.



---

*Source: 549562-2015-03-26.xml*