Parameter identification of solute transport with spatial fractional advection-dispersion equation via Tikhonov regularization

Accepted Manuscript Title: Parameter identification of solute transport with spatial fractional advection-dispersion equation via Tikhonov regularization Author: Hui Wei Wenjing Dang Song Wei PII: DOI: Reference:

S0030-4026(16)31219-0 http://dx.doi.org/doi:10.1016/j.ijleo.2016.10.052 IJLEO 58320

To appear in: Received date: Accepted date:

14-8-2016 18-10-2016

Please cite this article as: Hui Wei, Wenjing Dang, Song Wei, Parameter identification of solute transport with spatial fractional advection-dispersion equation via Tikhonov regularization, Optik - International Journal for Light and Electron Optics http://dx.doi.org/10.1016/j.ijleo.2016.10.052 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Parameter identification of solute transport with spatial fractional advection-dispersion equation via Tikhonov regularization

Hui WEI1,*, Wenjing Dang1,Song WEI2

1 College of Science, Anhui University of Science and Technology, Huainan, Anhui, China, 232001 2 Institute of Soft Matter Mechanics, Department of Engineering Mechanics, Hohai University, Nanjing, China, 210098

Abstract This paper employs Tikhonov regularization method to estimate the parameters of solute transport in soil with spatial fractional derivative advection-dispersion equation. The estimated parameters include the fractional derivative order, the dispersion coefficient and the average pore-water velocity. Several examples are included to test the presented method, and the method presented is efficient, feasible and easy to develop for the problem whose analytical solutions are difficult to be obtained. Keywords: Parameter identification; Fractional derivative advection-dispersion equation; Tikhonov regularization

1 Introduction Many studies have indicated that the solute transport processes can be descried by the conventional advection-dispersion equation (ADE) model based on Fick’s law, especially to simulate the contaminant transport in homogenous media. However, *

E-mail: [email protected] 1

most natural porous media, .i.e., natural soils or aquifers usually are heterogeneous. Due to the heterogeneity, solute transport processes in soil may be anomalous diffusion. It no longer follows Fick’s second law. The spatial distribution of the contaminant concentration is no longer in Gaussian form. The measured concentrations are usually much higher than those estimated by ADE at the early stage of the breakthrough curves, this phenomena was so called as anomalous or super-Fickian or non-Fickian transport[1]. The classical ADE fails to model the anomalous character of the solute transport in heterogeneous soil and other medium[2]. To more accurately describe the non-local property of anomalous diffusion (super-diffusion) in soil, fractional derivative has become a promising approach in recent years [2-7]. Non-local diffusion processes can be governed by a generalized space-fractional diffusion equation. It is obtained from the standard linear diffusion equation by replacing the second-order space derivative with a fractional derivative operator [8]. Based on Lévy motion theory, Benson et al. [9-11] described the spatial and temporal distribution of contaminant concentration by a fractional advection-dispersion equation (FADE). The one-dimensional FADE for describing non-reactive contaminant transport is given by C C 1  C 1  C  v  1    D   1    D , t x 2 x 2 ( x)

(1)

where C is the solute concentration, v is the average pore-water velocity, x is the spatial coordinate, t is the time, D is the dispersion coefficient with dimension of

[ L T 1 ] ,  ( 1    1 ) is the skewness, i.e. the relative weight of solute particle forward versus backward transition probability,  ( 1    2 ) is the order of fractional derivative, the FADE is reduced to ADE when  is equal to 2. The definitions of the fractional derivative operator can be [12]  m  x  C 1 C ( , t )  d ,  m    x (m   )  x   x    m 1

(2)

 C (1)m   m   C ( , t )  d .   ( x) (m   )  x m  x   x  m 1

(3)

2

The model has been used to simulate the non-Fickian process for conservative solute by Pachepsky[13] and Huang[14]. The more popular used FADE model with the symmetrical diffusion   0

can be written as follow[14]

C C  C  v D  . t x x

(4)

Nowadays, parameter determination has attracted more and more attention with the development of the fractional derivative model. Variety of software have been developed to estimate the parameter in ADE model, i.e., CXTFIT[15] . But to the best of our knowledge, there is very limited literatures on parameter estimation can be found for FADE model. Huang[14] developed a program, named as FADEMain based on FORTRAN, to estimate the parameters of Eq. (4), include the fractional order  , the dispersive coefficient D and the average pore-water velocity v. FADEMain is based on the nonlinear least square fitting algorithm[16]. Huang [17] employed the analytical solution of Eq. (4) to get the Hesse matrix in their method. However, analytical solutions of fractional derivative models are usually difficult to obtain or too complex to use. Therefore, the method presented by Huang is difficult to be extended into solving other problems. Hereby, we develop a numerical method to estimate the unknown parameters in Eq.(4), based on the Tikhonov regularization method. The new method does not need to involve the analytical solution of direct problem.

2 Methodology

2.1 Parameter Identification Problem We consider unabsorbed solute transport in a soil column with FADE model in a finite domain

C C  C  v  D  , 0  x  L,1    2 , t x x

(5)

subject to the following initial and boundary conditions

C (0, t )  0 , 3

(6)

C (0, t )  C0  .  C  0  x xL 

(7)

If all the parameters of (5)-(7) are given, we can solve the concentration distribution with time and space, it is so called direct problem. The numerical solution of this problem can be given by employing numerical methods, such as finite difference scheme [18-21], finite element method [22]. However, some parameters cannot be determined in prior or measured directly. Thus, we need to estimate the parameters via mathematical algorithms by adding additional conditions, which is so-called inverse problem, namely, parameter identification. From the experience of parameter identification of ADE, the breakthrough curves of solute should be considered as additional conditions. The observed concentration data at a particular observation point can be denoted as C obs (ti ), i  1, 2,

N,

(8)

Denote the unknown parameters as a vector p   v, D,   . Then the parameter T

identification is converted into a nonlinear optimization problem min F (p)  Cp  Cobs ,

(9)

p

where





Cobs  C obs (t1 ), C obs (t2 ),

Cp  C p (p, t1 ), C p (p, t2 ),

, C obs (t N )

, C p (p, t N )



T



T

represents

the

observed

data,

represents the computed results with the

estimated parameter vector p   v, D,   , . represents a norm, when we choose T

2-norm, the problem is reduced to a nonlinear least square problem 2

min F (p)  Cp  Cobs . p

(10)

However, this problem is generally ill-posed. To overcome the ill-posedness, we can employ the regularization approach to get a best fitted solution. By using the Tikhonov regularization, we can convert the problem (10) into the following form

4

min F (p)  Cp  Cobs    p  , 2

p

(11)

where  is the regularization parameter,   p  represents a stable functional. 2.2 Tikhonov regularization In this paper, we consider the zero-th Tikhonov regularization. Therefore, the problem (11) can be rewritten as 2

2

min F (p)  Cp  Cobs   p  p* , p

(12)

where p* is the exact solution. To minimize the problem (12), the first-order derivative of relation (12) should satisfy the following condition

Cp Cp  Cobs     p  p*   0 .  p 

(13)

Let p n1 be the (n+1)-th iteration parameter vector, expanding Cp (pn1 , ti ) in Taylor series near the point p n

Cp (p n 1 , ti )  Cp (p n , ti ) 

Cp p

p

n 1

 pn  ,

(14)

n

Substituting the relation (14) into (13) and replacing p* by p n , we have

Cp  Cp  p n  p Denote G 

Cp p

p n

n 1

  p n   Cp  Cobs     p n 1  p n   0 , 

(15)

and dp  pn1  pn , rewriting the relation (15) n

GT Gdp  GT Cp  Cobs    dp  0 ,

(16)

For the derivatives of the fitted concentration with respect to the parameter p j ( j  1, 2,

, M ) are evaluated according to

5

C p (ti ) C p ( p1 , p2 ,  p j

, p j  ,

, pM , ti )  C p ( p1 , p2 ,

, pj,

, pM , t i )



(17)

The current setting for  is 0.01 for all parameters, which we found to be appropriate for most cases[15]. Hence, G is a N  M matrix and Gij 

C p (ti ) . We p j

can obtain the best fitted solution of the original problem by iterating the relation(16). However, there still exists an important problem to be considered, how to choose the regularization parameter. In this study, we choose L-curve approach to determine regularization parameter. The detail of this method can be seen in ref. [23]. Strategy: (1) Given an initial guess parameter vector p (0) , and the iteration termination criterion  , set k  0 ; (2) Using the implicit finite difference to compute the concentration

C (p( k ) , ti ) i  1, 2,

, N , and compute G ( k )

by the relation (17), if

Cp  Cobs   , then p ( k ) should be the finial regularization solution of the

original problem, otherwise go to the step (3); (3) By using the L-curve approach, select an appropriate regularization parameter

 ( k ) , compute the value of dp( k ) by using the relation (16), and set k=k+1; (4) Let p( k 1)  p( k )  dp( k ) , then go to the step (2). Note that we employ the implicit finite difference method proposed by Meerschaert [19] to solve the direct problem in this study. All the programs are run in Matlab 2011b environment, Windows 7, 32 bits. The convergence and feasibility of the Tikhonov regularization have been studied well. Therefore, we no longer discuss the details in this study. Readers can refer to the literatures [24-28].

3 Application Example In this study, we consider a laboratory experiment conducted through 1250cm long, horizontally placed column packed with heterogeneity sandy soil [29]. NaCl was used as the tracer, the concentrations of Cl were measured with electrical conductivity at 6

100cm intervals in the column. We consider the second experiment, i.e., a tracer injection (transport) experiment in the heterogeneous sandy soil column by replacing inflowing tap water with a NaCl solution of concentration C0 = 6 g/L at the coarse end. The details of this experiment can be referred to [29]. In term of the analysis by Gao [30] and Pachepsky [13], we can model the solute transport by the fractional advection-dispersion equation mentioned above. To reflect the goodness-of –fit, we employ the coefficient of determination r 2 and the root mean square error (RMSE) N

r2  1

 C i 1 N

p

(ti )  C obs (ti ) 

 C

p

i 1

RMSE 

(ti )  C obs 

2

2

(18)

2 1 N C p (ti )  C obs (ti )    N i 1 

(19)

where C obs is the average value of the observations C obs (ti ), i  1, 2,

N , N is the

number of observed concentration data at a particular observation point. Compared results of Gao [30] and the method presented in this study is shown in the Table 1 and Figures 1-4. From the inverse results, we can see that the method presented in this study work as well as the FADEMain software.

The determination

coefficients is greater than 0.9. Because of the iteration error when applying the implicit difference, all values of RMSE are lightly greater than the FADEMain’s. However, without getting the analytical solution, it is more convenient to complement than FADEMain and easy to extend to other complex models.

4 Discussion The Tikhonov regularization method is sensitive to initial guess value. Thus, we need an initial guess value near the exact solution. If the initial guess value is too far away from the exact value, the iteration may be divergent and the best fitted value cannot be found. Another troubling issue is to choose appropriate regularization parameter. Variety 7

of principles has been proposed by the prior researchers on selecting regularization parameter, and the key step is choosing the best. Based on this study, we believe that the L-curve principle suitable to solve the considered problem in this study. However, further research is necessary to propose appropriate methods for choosing regularization parameter.

Acknowledgments The authors are grateful of Dr. HongGuang Sun for his valuable help in the preparation of this paper. The work described in this paper was supported by NSFC Grant No. 61472003.

References [1] Lévy, M. and B. Berkowitz, Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. Journal of Contaminant Hydrology, 2003. 4(3-4): p. 205-226. [2] Fomin, S., C. Chulsky, and T. Hashida, Mathematical modeling of anomalous diffusion in porous media. Fractional differential calculus, 2011. 1(1): p. 1-28. [3] Klafter, J., A. Blumen, and M.F. Shlesinger, Stochastic pathway to anomalous diffusion. Physical Review A, 1987. 35(7): p. 3081-3085. [4] Berkowitz, B. and H. Scher, Anomalous transport in random fracture networks. Physical Review Letters, 1997. 79(20): p. 4038-4041. [5] Cortis, A. and B. Berkowitz, Anomalous transport in "classical" soil and sand columns. Soil Science Society of America Journal, 2004. 68: p. 1539-1548. [6] Sun H.G, W. Chen, et al. Finite difference schemes for variable-order time fractional diffusion equation. International Journal of Bifurcation and Chaos (2012), 22 (4): 1250085. [7] Metzler, R. and J. Klafter, The random walk's guide to anomalous diffusion: a fractional dynamics approach. Physics Reports, 2000. 339(1): p. 1-77. [8] Gorenflo, R. and F. Mainardi, Random walk models for space-fractional diffusion processes. Fractional Calculus & Applied Analysis, 1998. 1(2): p. 167-191. [9] Benson, D.A., The fractional advection-dispersion equation:development and application. 1998, University of Nevada Reno. [10] Benson, D.A., et al., Fractional dispersion, Lévy motion, and the MADE tracer tests. Transport in Porous Media, 2001. 42(1): p. 211-240. [11] Benson, D.A., S.W. Wheatcraft, and M.M. Meerschaert, Application of a fractional advection-dispersion equation. Water Resources Research, 2000. 36(6): p. 1403-1412. [12] Samko, S.G., A.A. Kilbas, and O.I. Marchev, Fractional integrals and derivatives: theory and applications. 1993, New York: Gordon and Breach Science. [13] Pachepsky, Y., B. David, and R. Walter, Simulatingscale-dependent contaminant transport in soils with the fractional advective-dispersive equation. Soil Science Society of America Journal, 2000. 64: p. 8

1234-1243. [14] Huang, G., Q. Huang, and H. Zhan, et al.,, Modeling contaminant transport in homogenuous porous media with fractional advection-dispersion equation. Science in China, Ser. D Earth Sciencies, 2005. 48(Suppl.II): p. 295-302. [15] Toride, N., F. Leij, and M.T. Van Genuchten, The CXTFIT Code for Estimating Transport Parameters from Laboratory Or Filed Tracer Experiments. 1995: US Salinity Laboratory Riverside, CA. [16] Press, W.H., et al., Numerical Recipes in FORTRAN 77: Volume 1, Volume 1 of Fortran Numerical Recipes: The Art of Scientific Computing. Vol. 1. 1992: Cambridge university press. [17] Guanhua, H., Modeling adsorbing solute transport with fractional advection dispersion equation. 2003, China Agricultural University: Beijing. [18] Lin, Y. and C. Xu, Finite difference/spectral approximationa for the time-fractional diffusion equation. Journal of Computational Physics, 2007. 225(2): p. 1533-1553. [19] Meerschaert, M.M. and C. Tadjeran, Finite difference approximations

for fractional

advection-dispersion flow equations. Journal of Computation and Applied Mathematics, 2004. 172: p. 65-77. [20] Su, L., W. Wang, and Z. Yang, Finite difference approximations for the fractional advection–diffusion equation. Physics Letters A, 2009. 373(48): p. 4405-4408. [21] Yuste, S.B., Finite difference approximations for the fractional advection-diffusion equation. Journal of Computational Physics, 2006. 216: p. 264-274. [22] Deng, W., Finite element method for the space and time fractional Fokker-Planck equation. SIAM Journal on Numerical Analysis, 2008. 47(1): p. 204-226. [23] Hansen, P., REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems. Numerical Algorithms, 1994. 6: p. 1-35. [24] Tikhonov, A.N., Solution of incorrectly formulated problems and the regularization method. Soviet Math. Dokl., 1963. 4: p. 1035-1038. [25] Groetsch, C.W., The theory of Tikhonov regularization for Fredholm equations of the first kind. Vol. 105. 1984: Pitman Boston. [26] Golub, G.H., P.C. Hansen, and D.P. O'Leary, Tikhonov regularization and total least squares. SIAM Journal on Matrix Analysis and Applications, 1999. 21(1): p. 185-194. [27] Gfrerer, H., An a posteriori parameter choice for ordinary and iterated Tikhonov regularization of ill-posed problems leading to optimal convergence rates. Mathematics of computation, 1987. 49(180): p. 507-522, S5. [28] Natterer, F., Error bounds for Tikhonov regularization in Hilbert scales. Applicable Analysis, 1984. 18(1-2): p. 29-37. [29] Huang, K., N. Toride, and M.T. van Genuchten, Experimental investigation of solute transport in large, homogeneous and hetergeneous, saturated soil columns. Transport in Porous Media, 1995. 18: p. 283-302. [30] Guangyao, G., et al., Comparison of alternative models for simulating solute transport in heterogeneous soil at large scale. SHUILI XIAOBAO, 2010. 41(2): p. 164-172.

9

(a)

(b)

0

10

0.9

1 0.8

0.9

-1

10 0.7

0.8 -2

0.7

10

0.6

0.5

C/C0

C/C0

0.6

0.4

-3

10

0.5 0.4

-4

0.3

10

0.3

0.2

0.2 -5

0.1 0

10

Measured FADEMain Tikhonov Regularization 0

200

400

600

0.1 0

-6

10

800

t [min]

0

0

0.2 200

Measured FADEMain Tikhonov Regularization 0.4 0.6 0.8 1 400 600 800

t [min]

Figure 1 Comparison of NaCl concentration fitted by FADEMain and TR and the measured data at x=400cm for the tracer injection experiment in the heterogeneous soil column in which water flowed from the coarse-textured end to the fine-textured end: (a) linear axes and (b) semi-log axes.

10

(a)

(b)

0

1

10

0.9 -1

0.8

10

0.7 -2

10

C/C0

C/C0

0.6 0.5

-3

0.4

10

0.3 -4

0.2 0.1 0

10

Measured FADEMain Tikhonov Regularization

Measured FADEMain Tikhonov Regularization

-5

0

500

1000

1500

10

2000

0

500

1000

1500

2000

t [min]

t [min]

Figure 2 Comparison of NaCl concentration fitted by FADEMain and TR and the measured data at x=600cm for the tracer injection experiment in the heterogeneous soil column in which water flowed from the coarse-textured end to the fine-textured end: (a) linear axes and (b) semi-log axes.

11

(a)

(b)

0

1

10

0.9 -1

10

0.8 0.7

-2

10

C/C0

C/C0

0.6 0.5

-3

10

0.4 -4

10

0.3 0.2 0.1 0

-5

0

500

1000

1500

2000

Measured FADEMain Tikhonov Regualrization

10

Measured FADEMain Tikhonov Regualrization

-6

2500

t[min]

10

0

500

1000

1500

2000

2500

t[min]

Figure 3 Comparison of NaCl concentration fitted by FADEMain and TR and the measured data at x=800cm for the tracer injection experiment in the heterogeneous soil column in which water flowed from the coarse-textured end to the fine-textured end: (a) linear axes and (b) semi-log axes.

12

(a)

(b)

0

10

1 0.9

-1

10 0.8 0.7

-2

10

C/C0

C/C0

0.6 0.5

-3

10

0.4 -4

10 0.3 0.2

-5

0.1 0

0

500

1000

1500

2000

Measured FADEMain Tikhonov Regularization

10

Measured FADEMain Tikhonov Regularization

-6

10

2500

0

500

1000

1500

2000

2500

t [min]

t [min]

Figure 4 Comparison of NaCl concentration fitted by FADEMain and TR and the measured data at x=1000cm for the tracer injection experiment in the heterogeneous soil column in which water flowed from the coarse-textured end to the fine-textured end: (a) linear axes and (b) semi-log axes.

13

Table 1 The results of FADEMain and the method presented in this study Distance (cm) 400

600

800

1000

D



r2

RMSE

FADEMain 1.000

7.451

1.798

0.949

0.0763

TR

1.198

8.061

1.838

0.938

0.0837

FADEMain 1.140

11.01

1.660

0.986

0.0458

TR

10.20

1.776

0.968

0.0616

FADEMain 1.050

29.590

1.750

0.985

0.0430

TR

30.088

1.899

0.979

0.0522

FADEMain 0.900

22.98

1.780

0.992

0.0431

TR

19.01

1.924

0.982

0.0480

Method

v

1.128

0.976

0.813

TR-Tikhonov regularization

14