Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2011 (2011), Article ID 679531, 14 pages
Research Article

A Nonlinear Solute Transport Model and Data Reconstruction with Parameter Determination in an Undisturbed Soil-Column Experiment

1Institute of Applied Mathematics, Shandong University of Technology, Shandong, Zibo 255049, China
2Institute of Mining Technology, Inner Mongolia University of Technology, Hohhot, Inner Mongolia 010051, China
3Center of Analysis and Testing, Shandong University of Technology, Shandong, Zibo 255049, China

Received 22 June 2011; Accepted 19 October 2011

Academic Editor: J. Rodellar

Copyright © 2011 Gongsheng Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A real undisturbed soil-column infiltrating experiment in Zibo, Shandong, China, is investigated, and a nonlinear transport model for a solute ion penetrating through the column is put forward by using nonlinear Freundlich's adsorption isotherm. Since Freundlich's exponent and adsorption coefficient and source/sink terms in the model cannot be measured directly, an inverse problem of determining these parameters is encountered based on additional breakthrough data. Furthermore, an optimal perturbation regularization algorithm is introduced to determine the unknown parameters simultaneously. Numerical simulations are carried out and then the inversion algorithm is applied to solve the real inverse problem and reconstruct the measured data successfully. The computational results show that the nonlinear advection-dispersion equation discussed in this paper can be utilized by hydrogeologists to research solute transport behaviors with nonlinear adsorption in porous medium.

1. Introduction

Soil and groundwater pollution has become a serious threat to sustainable development throughout the world. It is important to characterize physical/chemical reactions quantitatively in the solute transport processes in the soil and groundwater. To understand the behaviors of the soil in the presence of infiltrating contaminants, soil-column experiments are often performed in laboratory or in field. There are disturbed and undisturbed soil-column experiments. For an undisturbed soil-column experiment, the structure of the soil layers and any contaminants within the column is preserved when the soil is transferred to the experimental apparatus. Generally speaking, for an undisturbed soil-column experiment, some complicated physical/chemical reactions could happen in the liquid phase and the solid phase when solute ions are being transported through the column with the inflow.

As we know, there are a lot of researches on solute transport behaviors in soil-column infiltrating experiments since the 1980s. Typical work could belong to the researching group of Nielsen and Van Genuchten. Nielsen et al. [1] put forward general equations of advection-dispersion type to describe solute transport behaviors, and Van Genuchten and Wagenet [2] constructed solute transport models of two sites/two regions in the soils. With development of computational tool and technique, numerical methods and software packages based on convection-dispersion equations are widely utilized on researches of soil-column infiltrating experiments (see, e.g., [3,4,5,6,7,8,9,10). Recently, the authors have even considered undisturbed soil-column experiments from two different aspects respectively. One aspect is to consider single-solute transport and identify source/sink parameters based on linear or nonlinear adsorption [11, 12]; the other is to deal with multicomponent solute transport and determine multiparameters based on hydrochemical analysis with advection dispersion mechanism [13].

It is obvious that research difficulties for soil-column infiltrating experiments lie in the construction of a suitable solute transport model and determination of model parameters. For an undisturbed soil-column experiment, the situation is more complicated due to difficulty of describing physical/chemical reactions for real problems. However, quite a few models employed linear adsorption isotherm leading to ordinary advection-dispersion equation with linear source/sink terms. As for other adsorption principles, such as Freundlich’s and Langmuir’s isotherms, there seem to be few applications on researches of real solutes transportation problems in mathematics. The reason maybe comes from that there are some parameters unknown in the adsorbing process, and these parameters can not be measured directly by ordinary experiments or they cannot be obtained except for spending much more cost.

On the other hand, there is a possible approach to get nonlinear adsorption parameters and source coefficients with less cost that is to apply inverse problem methods and effective inversion algorithms (see, e.g., [1420]). Actually, employing Freundlich’s adsorption principle, a nonlinear transport model can be obtained which is a nonlinear parabolic type equation in mathematics. However, Freundlich’s adsorption coefficient and exponent source/sink terms cannot be measured directly by the experiment; then an inverse problem of determining these model parameters is encountered based on measured breakthrough data. This paper will deal with a solute transport problem in an undisturbed soil-column infiltrating experiment based on nonlinear Fredunlich’s adsorption isotherm. A similar work has been presented in [12], but there is a problem left that is to determine all unknown transport parameters by using suitable inversion algorithms. For example, linear adsorption coefficient was regarded as the nonlinear adsorption coefficient in paper [12] which seems to be unscientific to some extent, and Fredunlich’s exponent was utilized empirically by testification.

In this paper, we will not only determine source/sink parameters but also determine the nonlinear adsorption coefficient and Fredunlich’s exponent simultaneously by the optimal perturbation regularization algorithm. Numerical inversions are carried out with which reasonable explanations to the experiment are obtained, and the measured data are reconstructed successfully. The inversion algorithm itself seems to be standard, but the paper will present detailed analysis on its concrete realization, for example, numerical convergence is testified, and several factors having important impacts on the algorithm are discussed.

The paper is arranged as follows. In Section 2, a nonlinear transport equation based on Fredunlich’s adsorption isotherm is introduced by which a mathematical model describing solute transport behaviors in an undisturbed soil-column experiment is put forward. In Section 3, an inverse problem of determining Freundlich’s adsorption coefficients and source/sink terms simultaneously is considered with the measured breakthrough data, and numerical simulations are carried out by applying an optimal perturbation regularization algorithm. In Section 4, the optimal inversion algorithm is implemented to determine the model parameters for the real soil-column experiment and then the measured data are reconstructed. Finally several concluding remarks are given.

2. A Mathematical Model and an Undisturbed Soil-Column Experiment

This paper is limited to utilize deterministic models to describe solute transport behaviors in soil-column infiltrating experiments. In general, solute transport in porous media often satisfies an advection-dispersion-reaction equation, but it is really complicated for undisturbed soil-column experiments due to unknown adsorbing process and uncontrollable physical/chemical reactions occurring in the column. However, in mathematics, nonlinear adsorption can be expressed by Fredunlich’s or Langmuir’s principles, and nonlinear reaction processes could be described with some nonlinear source/sink terms.

Denote by 𝑐 (ML−3) a solute concentration in liquid phase and by 𝑠𝑒 (MM−1) as the solute concentration in the adsorbed phase; then by mass conservation, there is (see, e.g., [3])𝜕𝜌𝜕𝑡𝑏𝑠𝑒𝜕+𝜃𝑐=𝐷𝜕𝑥𝜕(𝜃𝑐)𝜕𝑥𝑣𝜕(𝜃𝑐)𝜕𝑥𝜇𝑙𝜃𝑐𝜌𝑏𝜇𝑠,𝑒𝑠𝑒+𝛾𝑙𝜃𝑐+𝜌𝑏𝛾𝑠,𝑒𝑠𝑒,(2.1) where 𝜌𝑏 (ML−3) is the bulk density, 𝜃 (L3L−3) is the volumetric water content, 𝐷 (L2T−1) is the dispersion coefficient, 𝑣 (LT−1) is the average pore-water velocity, 𝜇𝑙 (T−1) is the first-order decay coefficient for the liquid phase, 𝜇𝑠,𝑒 (T−1) is the first-order decay coefficient for the adsorbed phase, 𝛾𝑙 (T−1) is the first-order production coefficient for the liquid phase, and 𝛾𝑠,𝑒 (T−1) is the first-order production coefficient for the adsorbed phase.

Suppose that nonlinear Fredunlich’s adsorption occurs between the solid and liquid phases in the column, that is, there is (see, e.g., [21])𝑠𝑒=𝐾𝐹𝑐𝑛,(2.2) where 𝐾𝐹 is Fredunlich’s adsorption coefficient and 𝑛>0 is called Freundlich’s exponent. By substituting expression (2.2) into (2.1), we can get𝑅(𝑐)𝜕𝑐𝜕𝜕𝑡=𝐷2𝑐𝜕𝑥2𝑣𝜕𝑐𝜕𝑥+𝜇1𝑐+𝜇2𝜌𝑏𝐾𝐹𝜃𝑐𝑛,(2.3) where𝜌𝑅(𝑐)=1+𝑏𝐾𝐹𝜃𝑛𝑐𝑛1(2.4) is called nonlinear retardation factor, and𝜇1=𝛾𝑙𝜇𝑙,𝜇(2.5)2=𝛾𝑠,𝑒𝜇𝑠,𝑒,(2.6) denote source/sink coefficients, respectively.

Equation (2.3) is a new mathematical model for describing 1D solute transport process with nonlinear adsorption behaviors in porous medium. If 𝑛=1, (2.3) is reduced to a linear advection-dispersion-reaction equation, in which case the linear adsorbing coefficient can be testified by lab experiment. However, if employing Fredunlich’s adsorption law, there are no effective ways to get the nonlinear adsorption coefficient 𝐾𝐹 and adsorption exponent 𝑛 by the experiment. In what follows, (2.3) will be applied to describe a solute transport process in an undisturbed soil-column experiment with help of an optimal perturbation regularization algorithm to determine the unknowns.

Let us first investigate an undisturbed soil-column experiment. (The experiment was supplied by The Inspecting Station of Geology and Environment in Zibo, Shandong.) The experiment was carried out in a laboratory in Zibo, Shandong, China, by taking an undisturbed soil-column nearby a coal mine region and infiltrating with the coal mine water. As we know, acid mine pollutants, for example, SO42, Cl, and Ca2+, Mg2+, are rich in coal-mine water. The aim of doing this experiment was mainly to reveal transport behaviors of the sulfate ions through the soil column.

This paper will take Ca2+ as an example to study its transport behaviors when it penetrates through the column. As for the experimental parameters, what we can obtain directly by experience and other experiments are listed as follows: length of the column, 𝑙 (cm), average pore-water velocity in the column, 𝑣 (cm/s), dispersion coefficient, 𝐷 (cm2/s), bulk density, 𝜌𝑏 (g/cm3), volumetric water content, 𝜃 (no dimension), and total experimental time infiltrating with the coal mine water, 𝑇1 (h). The values of the above-mentioned parameters are given in Table 1.

Table 1: Basic parameters in the soil-column experiment.

Before doing the experiment, the initial concentration of Ca2+ in the inflow was measured and denoted by 𝑐0=338.28 [mg/L]. On the other hand, throughout the experiment, the fluid that reaches the bottom of the column was collected and analyzed by which the so-called breakthrough data were obtained. In what follows, the measured breakthrough data for Ca2+ at the outflow are listed in Table 2.

Table 2: The measured breakthrough data (𝑡: time (h); 𝑐: concentration (mg/L)).

By Table 2, we find that the solute concentration in the first outflow at 𝑡=0.5 (h) is 𝑐(𝑙,0.5)=794.79 (mg/L), which is double that of 𝑐0=338.28 (mg/L) in the inflow, and the breakthrough data go down rapidly at the initial stage from 𝑡=0.5 (h) to 𝑡=4.1 (h) and then decrease gradually with the time going on. Maybe there was a rapid dissolution of ion species in the solid phase into the liquid phase at the initial stage. After the transient dissolution stage, that is, after 𝑡=0.5 (h), nonlinear adsorption reactions may play an important role in the solute transportation, and the solute concentration in the out-flow has a decreasing trend. So, we will utilize (2.3) based on nonlinear Fredunlich’s adsorption isotherm as the dominating model, and suppose that the primary source/sink coefficient 𝜇2 given by (2.6) is a time-dependent function, that is, 𝜇2=𝜇2(𝑡).

Now, let us give initial boundary value conditions for (2.3). Since the breakthrough data are available from the time of 𝑡=0.5 (h), we will take 𝑡=0.5 (h) as the initial moment. However, at the moment of 𝑡=0.5 (h), the solute concentration distribution in the liquid phase can not be measured directly throughout the column. What we know about the initial condition is that it may be monotonously increasing through the column, and 𝑐(0,0.5)=𝑐0 (mg/L), and 𝑐(𝑙,0.5)=794.79 (mg/L). According to the above information and by interpolation, we can get𝑐(𝑥,0.5)=𝑐0+794.79𝑐0𝑙𝑚𝑥𝑚,0𝑥𝑙,(2.7) where 𝑚>0 is called initial distribution index referring to characteristics of the solute distribution at initial moment of the first outflow.

As for boundary conditions, we will utilize ordinary conditions usually used for 1D soil-column experiment given as𝑐(0,𝑡)=𝑐0,𝑐𝑥(𝑙,𝑡)=0,0.5𝑡𝑇1.(2.8)

Obviously, if the retardation factor 𝑅(𝑐), Freundlich’s exponent 𝑛, the adsorption coefficient 𝐾𝐹, the source/sink coefficients 𝜇1 and 𝜇2(𝑡) in equation (2.3), and the hydraulic parameters 𝐷,𝑣 and the initial index 𝑚 are all known, the problem (2.3) with (2.7) and (2.8) is just an ordinary initial boundary value problem of partial differential equation of parabolic type which is called a forward problem. The forward problem can be solved numerically by using Matlab software. However, as stated in the above, the adsorption coefficient 𝐾𝐹 and exponent 𝑛 and the source/sink coefficients 𝜇1 and 𝜇2(𝑡) including the initial index 𝑚 cannot be obtained directly by the experiment. One thing we can do in mathematics is to determine them by some inversion algorithms with help of the measured breakthrough data which often leads to an inverse problem of parameter determination.

3. The Inverse Problem and the Inversion Algorithm

3.1. The Inverse Problem

For convenience of solving numerically, we will transform the model to a dimensionless form.

Denote 𝐶=𝑐/𝑐0, 𝑍=𝑥/𝑙, and 𝑇=𝑣𝑡/𝑙; substituting them into (2.3), we have𝑅(𝐶)𝜕𝐶=1𝜕𝑇𝑃𝜕2𝐶𝜕𝑍2𝜕𝐶+𝑙𝜕𝑍𝑣𝜇1𝐶+𝑙𝜌𝑏𝐾𝐹𝑐0𝑛1𝑣𝜃𝜇2(𝑇)𝐶𝑛,0<𝑍<1,𝑇0<𝑇<𝑇,(3.1) where 𝑅(𝐶)=1+(𝜌𝑏𝑐0𝑛1/𝜃)𝐾𝐹𝐶𝑛1, 𝜇2(𝑇)=𝜇2(𝑙/𝑣𝑇), and 𝑃=𝑙𝑣/𝐷, 𝑇0=0.5𝑣/𝑙, 𝑇=𝑇1𝑣/𝑙. Meantime, the initial condition (2.7) is transformed to𝐶𝑍,𝑇0=1+794.79𝑐0𝑍1𝑚,0𝑍1,(3.2) and the boundary conditions are transformed to𝐶(0,𝑇)=1;𝜕𝐶𝜕𝑍(1,𝑇)=0,𝑇0𝑇𝑇.(3.3)

The additional condition we will utilize is the measured breakthrough data listed in Table 2. Also by dimensionlessness to the real data, we have𝐶1,𝑇𝑘=𝐶𝑘,𝑘=1,2,,𝐾,(3.4) here 𝐾=16 according to Table 2.

As a result, an inverse problem of simultaneously determining the nonlinear exponent 𝑛, the adsorption coefficient 𝐾𝐹, the initial index 𝑚, and the source/sink coefficients 𝜇1 and 𝜇2(𝑇) is formulated by (3.1) with initial boundary conditions (3.2)-(3.3), and the overposed condition (3.4). In what follows, an optimal perturbation regularization algorithm (see, e.g., [11, 15]) will be introduced to solve the above inverse problem numerically, and several numerical simulations will be presented to verify the inversion algorithm validity.

3.2. The Optimal Perturbation Regularization Algorithm

For determining the above five kinds of parameters, 𝐾𝐹 and 𝑛, 𝜇1 and 𝜇2(𝑇), and 𝑚, we need a suitable approximate space to simulate 𝜇2(𝑇) due to its dependence on time variable. With similar method as used in [11], we will take lower-order polynomials as basis functions. Suppose that 𝜇2(𝑇) can be approximated by a quadratic polynomial given as follows:𝜇2(𝑇)=𝑎0+𝑎1𝑇+𝑎2𝑇2.(3.5) Then the model parameters we want to determine can be denoted by a vector𝐫=𝑛,𝐾𝐹,𝜇1,𝑎0,𝑎1,𝑎2,𝑚𝑇,(3.6) where 𝑛 is Freundlich’s exponent, 𝐾𝐹 is the nonlinear adsorption coefficient, 𝜇1 is the source/sink coefficient given by (2.5), and 𝑎0, 𝑎1, and 𝑎2 are referred to as in (3.5), and 𝑚 is the initial index coming out in (3.2). Obviously, the vector 𝐫 belongs to Euclid space 𝐑7, and in concrete computations we will set a bounded ball 𝑆𝐸={𝐫𝐑7𝐫𝐸} as an admissible set for the unknown parameters, where 𝐸 is a positive constant.

For any prescribed 𝐫𝑗 (𝑗=0,1,), set𝐫𝑗+1=𝐫𝑗+𝜖𝑗,𝑗=0,1,.(3.7) Thus, in order to get 𝐫𝑗+1 from the given 𝐫𝑗, we only need to compute an optimal perturbation 𝜖𝑗=(𝜖1𝑗,𝜖2𝑗,,𝜖7𝑗)𝑇𝑆𝐸. Denote by 𝐶(𝑍,𝑇;𝐫𝑗) the unique solution of the forward problem (3.1)–(3.3) for any given vector 𝐫𝑗𝑆𝐸. Taking Taylor’s expansion for 𝐶(𝑍,𝑇;𝐫𝑗+𝜖𝑗) at 𝐫𝑗 and ignoring higher-order terms, we can get𝐶𝑍,𝑇;𝐫𝑗+𝜖𝑗𝐶𝑍,𝑇;𝐫𝑗+𝑇𝐶𝑍,𝑇;𝐫𝑗𝜖𝑗.(3.8) So, taking 𝑍=1 and 𝑇=𝑇𝑘 (𝑘=1,2,,𝐾) in expression (3.8) and in view of optimality and stability, to determine a perturbation 𝜖𝑗 can be reduced to minimize the following multivariable error function of 𝜖𝑗𝑆𝐸:𝜖err𝑗=𝐶(1,𝑇𝑘;𝐫𝑗)+𝑇𝐶(1,𝑇𝑘;𝐫𝑗)𝜖𝑗𝐂22𝜖+𝛼𝑗22,(3.9) where the norm 2 denotes ordinary Euclid norm, 𝛼>0 is regularization parameter, and 𝐶𝐶=(1,𝐶2𝐶,,𝐾)𝑇 is the additional data vector given by (3.4). Note that𝑇𝐶1,𝑇𝑘;𝐫𝑗𝜖𝑗=7𝑖=1𝜕𝐶1,𝑇𝑘;𝐫𝑗𝜕𝑟𝑖𝑗𝜖𝑖𝑗7𝑖=1𝐶1,𝑇𝑘;𝐫𝑗+𝜏𝑖𝐞𝑖𝐶1,𝑇𝑘;𝐫𝑗𝜏𝑖𝜖𝑖𝑗,(3.10) where 𝐞𝑖=(0,,1,,0)𝑇 (𝑖=1,2,,7) are basis functions of 𝐑7 and 𝜏𝑖 (𝑖=1,2,,7) is numerical differential step. Furthermore, for given 𝐫𝑗, denote𝐔𝑗=𝐶1,𝑇1;𝐫𝑗,𝐶1,𝑇2;𝐫𝑗,,𝐶1,𝑇𝐾;𝐫𝑗𝑇,(3.11) and for 𝑘=1,2,,𝐾, 𝑖=1,2,,7, denote𝑔𝑗𝑘𝑖=𝐶1,𝑇𝑘;𝐫𝑗+𝜏𝑖𝐞𝑖𝐶1,𝑇𝑘;𝐫𝑗𝜏𝑖,𝐆𝑗=𝑔𝑗𝑘𝑖𝐾×7.(3.12) Then the error function (3.9) can be transformed to the following form:𝜖err𝑗=𝐆𝑗𝜖𝑗(𝐂𝐔𝑗)22𝜖+𝛼𝑗22.(3.13) By the general Tikhonov regularization theory (see, e.g., [22]), we know that the above minimization problem has one unique solution which can be expressed as𝜖𝛼𝑗=𝛼𝐈+𝐆𝑇𝑗𝐆𝑗1𝐆𝑇𝑗𝐂𝐔𝑗(3.14) for 𝑗=0,1,, where the regularization parameter 𝛼>0 should be chosen suitably.

Therefore, an optimal coefficient vector can be obtained approximately by iteration scheme (3.7) as long as a perturbation satisfyies a given convergent precision. This is the principal idea of optimal perturbation regularization algorithms. The detailed steps to implement the above algorithm are given as follows.

Step 1. Given initial iteration vector 𝐫𝑗, numerical differentiation step vector 𝜏=(𝜏1,𝜏2,,𝜏7)7, and convergent precision eps, and the additional measured data vector 𝐂.

Step 2. Solve the forward problem (3.1)–(3.3) with Matlab to get 𝐶(1,𝑇𝑘;𝐫𝑗) and 𝐶(1,𝑇𝑘;𝐫𝑗+𝜏𝑖𝑒𝑖), and then obtain the vector 𝐔𝑗 and the matrix 𝐆𝑗=(𝑔𝑗𝑘𝑖) by formula (3.12),

Step 3. Choose suitable regularization parameter 𝛼>0, and get an optimal perturbation vector 𝜖𝛼𝑗 by using formula (3.14).

Step 4. If there is 𝜖𝛼𝑗eps, then the inversion algorithm can be terminated, and 𝐫𝑗+1=𝐫𝑗+𝜖𝛼𝑗 is taken as the coefficient solution, what we just want to determine; otherwise, go to Step 2 by replacing 𝐫𝑗 with 𝐫𝑗+1.

4. Numerical Simulations

In order to verify numerical convergence of the inversion algorithm, several simulations will be presented by setting 𝐫true=(2,0.1,0.5,0.1,0.5,0.01,2)𝑇 as a true parameter vector in this section. We will reconstruct the true solution by applying the above optimal perturbation regularization algorithm. The additional data here are attained by solving the forward problem (3.1)–(3.3) with the above true parameters vector. It is noticeable that if choosing regularization parameter 𝛼=0, that is, without using regularization, the inversion computations always fail. In other words, regularization strategy must be utilized such that the inversion algorithm can be performed, and the regularization parameter is selected by numerical testification in this paper. In addition, all of the computations are performed in a PC of Dell Dimension 9200.

4.1. Solution Errors with Number of Iteration Times

Let us first investigate numerical convergence of the algorithm when the number of iteration times goes to infinity for given regularization parameters. Set the initial iteration as 𝐫0=𝟎 and the numerical differential step vector as 𝜏=(1𝑒2,1𝑒2,1𝑒3,1𝑒5,1𝑒7,1𝑒9,1𝑒11)𝑇. Table 3 and Figure 1 can show some convergence of the inversion algorithm. In Table 3, 𝑗 denotes the number of iteration times, 𝛼 is regularization parameter, and the solutions error is an average absolute error expressed byErrsolu=𝐫true𝐫inv2=7𝑖=117||𝑟true𝑖𝑟inv𝑖||21/2,(4.1) where 𝐫inv denotes an inversion reconstruction solution corresponding to the true solution 𝐫true.

Table 3: Solution errors with iteration numbers and regularization parameters.
Figure 1: Solution errors with iteration numbers for different regularization parameters.

By Table 3 and Figure 1, we can see that the inversion algorithm is of numerical convergence for suitable regularization parameters. The solutions errors become smaller with larger number of iteration times for given regularization parameters. As for similar iteration times, the solution errors are also smaller with smaller regularization parameters. For example, in the case of 𝑗=1000, the solution error is Errsolu=1.04143𝑒2 for 𝛼=1𝑒3, but it becomes Errsolu=1.73988𝑒6 for 𝛼=1𝑒4. In addition, we also find that solutions errors become very small and have little changes after thousands of iteration times for given regularization parameters.

4.2. Solution Errors with Regularization Parameters

In this subsection, we will investigate changes of solution errors with regularization parameters for given convergent precision. Also, set initial iteration and differential step vector as used in the last subsection and take convergent precision as eps=1𝑒8. Figure 2 plots the solution errors with regularization parameters, where the straight line represents a linear fitting of all of the solution errors, whose equation is expressed byErrsolu[](𝛼)=15.4158𝛼+0.0056,𝛼2𝑒5,3𝑒3.(4.2)

Figure 2: Solution errors with regularization parameters and fitting line.

By the computations, we find that numerical inversions should be performed for regularization parameters taking values in the interval of 𝛼[2𝑒5,3𝑒3]. Even though the solution errors become small as regularization parameters become small, regularization parameters cannot be chosen too small. On the other hand, since data noises are ignored in performing the inversion algorithm, the solution errors should be in a linear relation with regularization parameters by general regularization theory. Actually, by the general Tikhonov regularization theory, there should be𝑟inv𝑟true2𝐸𝛼,(4.3) where 𝑟inv represents the inversion solution, 𝑟true represents the true solution, and 𝐸 is a positive constant. Fortunately, by the above computations, we can find that the inversion results basically coincide with the theoretical analysis of the Tikhonov regularization.

4.3. Solution Errors with Convergent Precision

As stated in Section 3.2, the inversion algorithm can be terminated if there is 𝜖𝛼𝑗2eps, where eps denotes convergent precision. We will show that solution errors should become small and go to zero as convergent precision goes to zero. By choosing regularization parameter as 𝛼=2𝑒5 and initial iteration and differential steps as before, the solution errors are listed in Table 4 and plotted in Figure 3, respectively, where eps also denotes convergent precision, Errsolu denotes solutions error defined by (4.1), abscissa represents logarithmic of eps, and longitudinal coordinates denote solutions errors in Figure 3.

Table 4: Solution errors with convergent precision.
Figure 3: Solution errors with convergent precision.

Finally, by choosing regularization parameter as 𝛼=2𝑒5, convergent precision as eps=1𝑒16, and other parameters as before, the true parameter vector can be reconstructed by the inversion algorithm given as 𝐫inv=(2.00000,0.100000,0.500000,0.999998𝑒1,0.500000,0.100000𝑒1,2.00000),(4.4) which is very close to the true solution and the solutions error is Errsolu=7.145𝑒8.

5. Inversion Results for the Inverse Problem

Now we will apply the inversion algorithm to solve the real inverse problem (3.1)–(3.3). By the above numerical simulations, we know that regularization parameter and convergent precision both play important roles in the algorithm realization. However, for real problems, regularization parameter and convergent precision often have to be chosen suitably larger than artificial simulations due to noises of real data. Based on the above numerical simulations, we will perform the inversion algorithm on the real inverse problem (3.1)–(3.3) utilizing larger parameters.

Also, setting initial iteration and numerical differentiation step vector as before and choosing regularization parameter as 𝛼=0.08 and convergent precision as eps=5𝑒6, an optimal solution of the inverse problem (3.1)–(3.3) was worked out by the inversion algorithm by 1351-time iterations which cost 1422.9 seconds of CPU time. The inversion coefficient solution is given as follows:𝐫inv=(0.8649,1.753𝑒2,7.103,19.77,4.863𝑒2,6.736𝑒4,1.475)𝑇.(5.1) That is, the five coefficients in the dimensionless model (3.1)–(3.3) are given as𝑛=𝑟inv1=0.8649,𝐾𝐹=𝑟inv2𝜇=0.01753,(5.2)1=𝑟inv3=7.103,𝑚=𝑟inv7=1.475,(5.3)𝜇2(𝑇)=19.77+0.04863𝑇+0.0006736𝑇2,(5.4) respectively. Furthermore, in order to see the inversion results visibly, we substitute the above inversion parameters into the forward problem (3.1)–(3.3) and get the reconstruction data of the solute concentration at 𝑍=1, which is plotted in Figure 4 compared with the actually measured breakthrough data.

Figure 4: Reconstruction data and the measured breakthrough data.

By Figure 4, we can see that the computational reconstruction data coincide with the measured breakthrough data very well. On the other hand, to quantify the goodness-of-fit, the absolute error in Euclid norm and the corresponding relative error are worked out and expressed asErrinv=𝐶(1,𝑇;𝐫inv𝐂)2=𝐾𝑘=11𝐾||𝐶1,𝑇𝑘;𝐫inv𝐶𝑘||21/2=0.008696,Errrel=Errinv𝐂2=0.001507,(5.5) respectively, where 𝐶(1,𝑇𝑘;𝐫inv) (𝑘=1,2,,𝐾) denote reconstruction data by solving the forward problem with the inversion solution 𝐫inv, 𝐶𝑘 (𝑘=1,2,,𝐾) also denote the measured breakthrough data, and 𝐶𝐂=(1,𝐶2𝐶,,𝐾)𝑇.

6. Concluding Remarks

Remark 6.1. This paper deals with numerical inversion and simulation to a real inverse problem arising from a soil-column infiltrating experiment. By employing Freundlich’s adsorption principle, a nonlinear transport model is put forward, and the measured breakthrough data are reconstructed successfully by applying an optimal perturbation regularization algorithm to determine the unknown parameters. The inversion results show that the solute transport process in the column can be described by the proposed nonlinear equation (2.3) with suitable initial boundary value conditions.

Remark 6.2. The optimal perturbation regularization algorithm is efficient to the inverse problem of determining model parameters numerically. By performing the algorithm, we find that there are several factors having important impacts on the algorithm realization, which are regularization parameter, convergent precision or number of iteration times, numerical differential steps, initial iterations, computational errors of the forward problem, and so forth. However, for the inverse problem studied here, numerical differential steps and initial iterations both have little impact on the algorithm, but regularization parameter and convergent precision both have important impact on the algorithm realization. An approximate solution could be obtained by the inversion algorithm when taking suitable regularization parameters and making number of iteration times go to infinity. Furthermore, for given convergent precision or iteration number, the solutions errors become small with approximate linearity as regularization parameter tends to zero due to free noises of data, which just coincides with general theory of the Tikhonov regularization.

Remark 6.3. Consider hydrogeological meanings of the inversion parameters. By the measured breakthrough data, we know that the solute concentration reaches its peak at the time of 𝑡=0.5 (h), and after that the soil adsorbing capability plays a dominating role in the solute transport process. The adsorption principle basically agrees with the nonlinear Fredunlich’s isotherm of 𝑠𝑒=0.01753𝑐0.8649. As for the source/sink coefficients, by the inversion results (5.4) and (5.3) and noting that 𝜇2(𝑇)=𝜇2(𝑡) and 𝑇=𝑣𝑡/𝑙, we can get the source coefficient in real-time dimension given as 𝜇2(𝑡)=19.77+0.01460𝑡+0.00006095𝑡2,(6.1) and the sink parameter is 𝜇1=7.102.(6.2) Since 𝜇2(𝑡)=𝛾𝑠,𝑒𝜇𝑠,𝑒0 for 𝑡>0.5 (h) and 𝜇1=𝛾𝑙𝜇𝑙<0, we can deduce that there are strong reactions in the adsorbed phase where the ion production rate is much higher than its decay rate, and it is on the contrary for the liquid phase where the ion production rate is lower than the decay rate.

Finally, we will briefly discuss for the uniqueness of the source function 𝜇2(𝑡). By expression (6.1) in dimension form, we find that the coefficient of the third term of 𝜇2(𝑡) is 0.00006095, which is very small as compared with coefficients of the first two terms. Actually, by numerical computations we find that if 𝜇2(𝑡) takes higher-order polynomials, its expansion coefficients of higher-order terms also go to zero showing a numerical uniqueness of it.


This project is supported by the National Natural Science Foundation of China (nos. 11071148, 10926194, and 10471080).


  1. D. R. Nielsen, M. T. van Genuchten, and J. W. Biggar, “Water flow and solute transport process in the unsaturated zone,” Water Resources Research, vol. 22, no. 9, pp. 89S–108S, 1986. View at Google Scholar · View at Scopus
  2. M. T. Van Genuchten and R. J. Wagenet, “Two-site/two-region models for pesticide transport and degradation: theoretical development and analytical solutions,” Soil Science Society of America Journal, vol. 53, no. 5, pp. 1303–1310, 1989. View at Google Scholar
  3. N. Toride, F. J. Leij, and M. T. Van Genuchten, “The CXTFIT code for estimating transport parameters from laboratory or field tracer experiments,” Version 2.0, Research Report 137, U.S. Department of Agriculture, 1995. View at Google Scholar
  4. Z. Torsten, “Capability of convection-dispersion transport models to predict transient water and solute movement in undisturbed soil columns,” Journal of Contaminant Hydrology, vol. 30, no. 1-2, pp. 101–128, 1998. View at Publisher · View at Google Scholar · View at Scopus
  5. L. Pang and M. E. Close, “Non-equilibrium transport of Cd in alluvial gravels,” Journal of Contaminant Hydrology, vol. 36, no. 1-2, pp. 185–206, 1999. View at Publisher · View at Google Scholar · View at Scopus
  6. M. Inoue, J. Šimůnek, S. Shiozawa, and J. W. Hopmans, “Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments,” Advances in Water Resources, vol. 23, no. 7, pp. 677–688, 2000. View at Publisher · View at Google Scholar
  7. S. K. Kamra, B. Lennartz, M. T. Van Genuchten, and P. Widmoser, “Evaluating non-equilibrium solute transport in small soil columns,” Journal of Contaminant Hydrology, vol. 48, no. 3-4, pp. 189–212, 2001. View at Publisher · View at Google Scholar · View at Scopus
  8. K. Cui, B. Y. Li, X. S. Li, and G. W. Yang, “Model parameter inversion for cadmium ion transport through unsaturated soils,” Advances in Water Science, vol. 15, no. 6, pp. 700–705, 2004. View at Google Scholar · View at Scopus
  9. J. M. Köhne, B. P. Mohanty, and J. Šimůnek, “Inverse dual-permeability modeling of preferential water flow in a soil column and implications for field-scale solute transport,” Vadose Zone Journal, vol. 5, no. 1, pp. 59–76, 2006. View at Publisher · View at Google Scholar
  10. D. A. Barry, “Effect of nonuniform boundary conditions on steady flow in saturated homogeneous cylindrical soil columns,” Advances in Water Resources, vol. 32, no. 4, pp. 522–531, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. G. S. Li, J. Cheng, D. Yao, H. L. Liu, and J. J. Liu, “One-dimensional equilibrium model and source parameter determination for soil-column experiment,” Applied Mathematics and Computation, vol. 190, no. 2, pp. 1365–1374, 2007. View at Publisher · View at Google Scholar · View at Scopus
  12. G. S. Li, Y. J. Tan, D. Yao, X. Q. Wang, and H. L. Liu, “A non-linear mathematical model for an undisturbed soil-column experiment and source parameter identification,” Inverse Problems in Science and Engineering, vol. 16, no. 7, pp. 885–901, 2008. View at Publisher · View at Google Scholar
  13. G. S. Li, D. Yao, Y. Z. Wang, and H. Y. Jiang, “Numerical inversion of multi-parameters in multi-components reactive solutes transportation in an undisturbed soil-column experiment,” Computer Modeling in Engineering & Sciences, vol. 51, no. 1, pp. 53–72, 2009. View at Google Scholar
  14. N. Z. Sun, Inverse Problem in Groundwater Modelling, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1994.
  15. C. W. Su, Numerical Methods and Applications of Inverse Problems in PDE, Northwestern Polytechnical University Press, Xi'an, China, 1995.
  16. J. Atmadja and A. C. Bagtzoglou, “State of the art report on mathematical methods for groundwater pollution source identification,” Environmental Forensics, vol. 2, no. 3, pp. 205–214, 2001. View at Publisher · View at Google Scholar · View at Scopus
  17. P. S. Mahar and B. Datta, “Optimal identification of ground-water pollution sources and parameter estimation,” Journal of Water Resources Planning and Management, vol. 127, no. 1, pp. 20–29, 2001. View at Publisher · View at Google Scholar · View at Scopus
  18. J. Zueco and F. Alhama, “Simultaneous inverse determination of temperature-dependent thermophysical properties in fluids using the network simulation method,” International Journal of Heat and Mass Transfer, vol. 50, no. 15-16, pp. 3234–3243, 2007. View at Publisher · View at Google Scholar · View at Scopus
  19. A. Khlaifi, A. Ionescu, and Y. Candau, “Pollution source identification using a coupled diffusion model with a genetic algorithm,” Mathematics and Computers in Simulation, vol. 79, no. 12, pp. 3500–3510, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  20. F. Yang and C.-L. Fu, “The method of simplified Tikhonov regularization for dealing with the inverse time-dependent heat source problem,” Computers & Mathematics with Applications, vol. 60, no. 5, pp. 1228–1236, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. C. M. Zheng and G. D. Bennett, Applied Contaminant Transport Modeling, John Wiley & Sons, New York, NY, USA, 2nd edition, 2002.
  22. A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, vol. 120, Springer, New York, NY, USA, 1996.