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

Authors: Gongsheng Li; De Yao; Yongzai Wang; Xianzheng Jia
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

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.


---

## Body

## 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., [14–20]). 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 Section2, 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 byc (ML−3) a solute concentration in liquid phase and by se (MM−1) as the solute concentration in the adsorbed phase; then by mass conservation, there is (see, e.g., [3])(2.1)∂∂t(ρbse+θc)=D∂∂x∂(θc)∂x-v∂(θc)∂x-μlθc-ρbμs,ese+γlθc+ρbγs,ese,
where ρb (ML−3) is the bulk density, θ (L3L−3) is the volumetric water content, D (L2T−1) is the dispersion coefficient, v (LT−1) is the average pore-water velocity, μl (T−1) is the first-order decay coefficient for the liquid phase, μs,e (T−1) is the first-order decay coefficient for the adsorbed phase, γl (T−1) is the first-order production coefficient for the liquid phase, and γs,e (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)se=KFcn,
where KF is Fredunlich’s adsorption coefficient and n>0 is called Freundlich’s exponent. By substituting expression (2.2) into (2.1), we can get(2.3)R(c)∂c∂t=D∂2c∂x2-v∂c∂x+μ1c+μ2ρbKFθcn,
where(2.4)R(c)=1+ρbKFθncn-1
is called nonlinear retardation factor, and(2.5)μ1=γl-μl,(2.6)μ2=γs,e-μs,e,
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 n=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 KF and adsorption exponent n 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 takeCa2+ 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, l (cm), average pore-water velocity in the column, v (cm/s), dispersion coefficient, D (cm2/s), bulk density, ρb (g/cm3), volumetric water content, θ (no dimension), and total experimental time infiltrating with the coal mine water, T1 (h). The values of the above-mentioned parameters are given in Table 1.Table 1
Basic parameters in the soil-column experiment.
lvDρbθT1453.76e-33.76e-31.480.20120Before doing the experiment, the initial concentration ofCa2+ in the inflow was measured and denoted by c0=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 (t: time (h); c: concentration (mg/L)).
t0.51.12.14.17.1101725c794.79628.65544.49510.42498.60490.60480.56464.53t364758698295108120c442.48430.46410.42396.39388.38382.36378.36376.35By Table2, we find that the solute concentration in the first outflow at t=0.5 (h) is c(l,0.5)=794.79 (mg/L), which is double that of c0=338.28 (mg/L) in the inflow, and the breakthrough data go down rapidly at the initial stage from t=0.5 (h) to t=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 t=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(t).Now, let us give initial boundary value conditions for (2.3). Since the breakthrough data are available from the time of t=0.5 (h), we will take t=0.5 (h) as the initial moment. However, at the moment of t=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 c(0,0.5)=c0 (mg/L), and c(l,0.5)=794.79 (mg/L). According to the above information and by interpolation, we can get(2.7)c(x,0.5)=c0+794.79-c0lmxm,0≤x≤l,
where m>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(2.8)c(0,t)=c0,cx(l,t)=0,0.5≤t≤T1.Obviously, if the retardation factorR(c), Freundlich’s exponent n, the adsorption coefficient KF, the source/sink coefficients μ1 and μ2(t) in equation (2.3), and the hydraulic parameters D,v and the initial index m 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 KF and exponent n and the source/sink coefficients μ1 and μ2(t) including the initial index m 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.DenoteC=c/c0, Z=x/l, and T=vt/l; substituting them into (2.3), we have(3.1)R(C)∂C∂T=1P∂2C∂Z2-∂C∂Z+lvμ1C+lρbKFc0n-1vθμ¯2(T)Cn,0<Z<1,T0<T<T¯,
where R(C)=1+(ρbc0n-1/θ)KFCn-1, μ¯2(T)=μ2(l/vT), and P=lv/D, T0=0.5v/l, T¯=T1v/l. Meantime, the initial condition (2.7) is transformed to(3.2)C(Z,T0)=1+(794.79c0-1)Zm,0≤Z≤1,
and the boundary conditions are transformed to(3.3)C(0,T)=1;∂C∂Z(1,T)=0,T0≤T≤T¯.The additional condition we will utilize is the measured breakthrough data listed in Table2. Also by dimensionlessness to the real data, we have(3.4)C(1,Tk)=Ĉk,k=1,2,…,K,
here K=16 according to Table 2.As a result, an inverse problem of simultaneously determining the nonlinear exponentn, the adsorption coefficient KF, the initial index m, and the source/sink coefficients μ1 and μ¯2(T) 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,KF and n, μ1 and μ¯2(T), and m, we need a suitable approximate space to simulate μ¯2(T) 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(T) can be approximated by a quadratic polynomial given as follows:(3.5)μ¯2(T)=a0+a1T+a2T2.
Then the model parameters we want to determine can be denoted by a vector(3.6)r=(n,KF,μ1,a0,a1,a2,m)T,
where n is Freundlich’s exponent, KF is the nonlinear adsorption coefficient, μ1 is the source/sink coefficient given by (2.5), and a0, a1, and a2 are referred to as in (3.5), and m is the initial index coming out in (3.2). Obviously, the vector r belongs to Euclid space R7, and in concrete computations we will set a bounded ball SE={r∈R7:∥r∥≤E} as an admissible set for the unknown parameters, where E is a positive constant.For any prescribedrj (j=0,1,…), set(3.7)rj+1=rj+ϵj,j=0,1,….
Thus, in order to get rj+1 from the given rj, we only need to compute an optimal perturbation ϵj=(ϵj1,ϵj2,…,ϵj7)T∈SE. Denote by C(Z,T;rj) the unique solution of the forward problem (3.1)–(3.3) for any given vector rj∈SE. Taking Taylor’s expansion for C(Z,T;rj+ϵj) at rj and ignoring higher-order terms, we can get(3.8)C(Z,T;rj+ϵj)≈C(Z,T;rj)+∇TC(Z,T;rj)⋅ϵj.
So, taking Z=1 and T=Tk (k=1,2,…,K) in expression (3.8) and in view of optimality and stability, to determine a perturbation ϵj can be reduced to minimize the following multivariable error function of ϵj∈SE:(3.9)err⁡(ϵj)=‖C(1,Tk;rj)+∇TC(1,Tk;rj)⋅ϵj-Ĉ‖22+α‖ϵj‖22,
where the norm ∥·∥2 denotes ordinary Euclid norm, α>0 is regularization parameter, and C=(Ĉ1,Ĉ2,⋯,ĈK)T is the additional data vector given by (3.4). Note that(3.10)∇TC(1,Tk;rj)⋅ϵj=∑i=17∂C(1,Tk;rj)∂rjiϵji≈∑i=17C(1,Tk;rj+τiei)-C(1,Tk;rj)τiϵji,
where ei=(0,…,1,…,0)T (i=1,2,…,7) are basis functions of R7 and τi (i=1,2,…,7) is numerical differential step. Furthermore, for given rj, denote(3.11)Uj=(C(1,T1;rj),C(1,T2;rj),…,C(1,TK;rj))T,
and for k=1,2,…,K, i=1,2,…,7, denote(3.12)gkij=C(1,Tk;rj+τiei)-C(1,Tk;rj)τi,Gj=(gkij)K×7.
Then the error function (3.9) can be transformed to the following form:(3.13)err⁡(ϵj)=‖Gjϵj-(Ĉ-Uj)‖22+α‖ϵj‖22.
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(3.14)ϵjα=(αI+GjTGj)-1GjT(Ĉ-Uj)
for j=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 vectorrj, numerical differentiation step vector τ=(τ1,τ2,…,τ7)7, and convergent precision eps, and the additional measured data vector Ĉ.Step 2.
Solve the forward problem (3.1)–(3.3) with Matlab to get C(1,Tk;rj) and C(1,Tk;rj+τiei), and then obtain the vector Uj and the matrix Gj=(gkij) by formula (3.12),Step 3.
Choose suitable regularization parameterα>0, and get an optimal perturbation vector ϵjα by using formula (3.14).Step 4.
If there is∥ϵjα∥≤eps, then the inversion algorithm can be terminated, and rj+1=rj+ϵjα is taken as the coefficient solution, what we just want to determine; otherwise, go to Step 2 by replacing rj with rj+1.


## 3.1. The Inverse Problem

For convenience of solving numerically, we will transform the model to a dimensionless form.DenoteC=c/c0, Z=x/l, and T=vt/l; substituting them into (2.3), we have(3.1)R(C)∂C∂T=1P∂2C∂Z2-∂C∂Z+lvμ1C+lρbKFc0n-1vθμ¯2(T)Cn,0<Z<1,T0<T<T¯,
where R(C)=1+(ρbc0n-1/θ)KFCn-1, μ¯2(T)=μ2(l/vT), and P=lv/D, T0=0.5v/l, T¯=T1v/l. Meantime, the initial condition (2.7) is transformed to(3.2)C(Z,T0)=1+(794.79c0-1)Zm,0≤Z≤1,
and the boundary conditions are transformed to(3.3)C(0,T)=1;∂C∂Z(1,T)=0,T0≤T≤T¯.The additional condition we will utilize is the measured breakthrough data listed in Table2. Also by dimensionlessness to the real data, we have(3.4)C(1,Tk)=Ĉk,k=1,2,…,K,
here K=16 according to Table 2.As a result, an inverse problem of simultaneously determining the nonlinear exponentn, the adsorption coefficient KF, the initial index m, and the source/sink coefficients μ1 and μ¯2(T) 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,KF and n, μ1 and μ¯2(T), and m, we need a suitable approximate space to simulate μ¯2(T) 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(T) can be approximated by a quadratic polynomial given as follows:(3.5)μ¯2(T)=a0+a1T+a2T2.
Then the model parameters we want to determine can be denoted by a vector(3.6)r=(n,KF,μ1,a0,a1,a2,m)T,
where n is Freundlich’s exponent, KF is the nonlinear adsorption coefficient, μ1 is the source/sink coefficient given by (2.5), and a0, a1, and a2 are referred to as in (3.5), and m is the initial index coming out in (3.2). Obviously, the vector r belongs to Euclid space R7, and in concrete computations we will set a bounded ball SE={r∈R7:∥r∥≤E} as an admissible set for the unknown parameters, where E is a positive constant.For any prescribedrj (j=0,1,…), set(3.7)rj+1=rj+ϵj,j=0,1,….
Thus, in order to get rj+1 from the given rj, we only need to compute an optimal perturbation ϵj=(ϵj1,ϵj2,…,ϵj7)T∈SE. Denote by C(Z,T;rj) the unique solution of the forward problem (3.1)–(3.3) for any given vector rj∈SE. Taking Taylor’s expansion for C(Z,T;rj+ϵj) at rj and ignoring higher-order terms, we can get(3.8)C(Z,T;rj+ϵj)≈C(Z,T;rj)+∇TC(Z,T;rj)⋅ϵj.
So, taking Z=1 and T=Tk (k=1,2,…,K) in expression (3.8) and in view of optimality and stability, to determine a perturbation ϵj can be reduced to minimize the following multivariable error function of ϵj∈SE:(3.9)err⁡(ϵj)=‖C(1,Tk;rj)+∇TC(1,Tk;rj)⋅ϵj-Ĉ‖22+α‖ϵj‖22,
where the norm ∥·∥2 denotes ordinary Euclid norm, α>0 is regularization parameter, and C=(Ĉ1,Ĉ2,⋯,ĈK)T is the additional data vector given by (3.4). Note that(3.10)∇TC(1,Tk;rj)⋅ϵj=∑i=17∂C(1,Tk;rj)∂rjiϵji≈∑i=17C(1,Tk;rj+τiei)-C(1,Tk;rj)τiϵji,
where ei=(0,…,1,…,0)T (i=1,2,…,7) are basis functions of R7 and τi (i=1,2,…,7) is numerical differential step. Furthermore, for given rj, denote(3.11)Uj=(C(1,T1;rj),C(1,T2;rj),…,C(1,TK;rj))T,
and for k=1,2,…,K, i=1,2,…,7, denote(3.12)gkij=C(1,Tk;rj+τiei)-C(1,Tk;rj)τi,Gj=(gkij)K×7.
Then the error function (3.9) can be transformed to the following form:(3.13)err⁡(ϵj)=‖Gjϵj-(Ĉ-Uj)‖22+α‖ϵj‖22.
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(3.14)ϵjα=(αI+GjTGj)-1GjT(Ĉ-Uj)
for j=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 vectorrj, numerical differentiation step vector τ=(τ1,τ2,…,τ7)7, and convergent precision eps, and the additional measured data vector Ĉ.Step 2.
Solve the forward problem (3.1)–(3.3) with Matlab to get C(1,Tk;rj) and C(1,Tk;rj+τiei), and then obtain the vector Uj and the matrix Gj=(gkij) by formula (3.12),Step 3.
Choose suitable regularization parameterα>0, and get an optimal perturbation vector ϵjα by using formula (3.14).Step 4.
If there is∥ϵjα∥≤eps, then the inversion algorithm can be terminated, and rj+1=rj+ϵjα is taken as the coefficient solution, what we just want to determine; otherwise, go to Step 2 by replacing rj with rj+1.


## 4. Numerical Simulations

In order to verify numerical convergence of the inversion algorithm, several simulations will be presented by settingrtrue=(2,0.1,0.5,0.1,-0.5,0.01,2)T 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 asr0=0 and the numerical differential step vector as τ=(1e-2,1e-2,1e-3,1e-5,1e-7,1e-9,1e-11)T. Table 3 and Figure 1 can show some convergence of the inversion algorithm. In Table 3, j denotes the number of iteration times, α is regularization parameter, and the solutions error is an average absolute error expressed by(4.1)Errsolu=‖rtrue-rinv‖2=(∑i=1717|ritrue-riinv|2)1/2,
where rinv denotes an inversion reconstruction solution corresponding to the true solution rtrue.Table 3
Solution errors with iteration numbers and regularization parameters.
jα=1e-3α=5e-4α=1e-4α=5e-5100.2603770.3073080.1722560.356620200.1467210.1211664.31409e-20.286733408.28802e-24.87719e-22.77992e-20.1978731003.71804e-21.04950e-21.50932e-26.41182e-22002.43731e-25.33821e-35.50282e-38.79467e-34001.91484e-23.53564e-37.33347e-41.60040e-46001.56244e-22.35879e-39.78162e-52.90766e-68001.27555e-21.57379e-31.30485e-55.39138e-810001.04143e-21.05009e-31.73988e-63.26469e-912008.50343e-37.00688e-42.31128e-73.11220e-9Figure 1
Solution errors with iteration numbers for different regularization parameters.By Table3 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 j=1000, the solution error is Errsolu=1.04143e-2 for α=1e-3, but it becomes Errsolu=1.73988e-6 for α=1e-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 aseps=1e-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 by(4.2)Errsolu(α)=15.4158α+0.0056,α∈[2e-5,3e-3].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α∈[2e-5,3e-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(4.3)‖rinv-rtrue‖2≤Eα,
where rinv represents the inversion solution, rtrue represents the true solution, and E 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 Section3.2, the inversion algorithm can be terminated if there is ∥ϵjα∥2≤eps, 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 α=2e-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.
eps1e-41e-61e-81e-101e-121e-141e-16Errsolu7.383e-27.206e-37.427e-47.294e-57.164e-67.382e-77.145e-8Figure 3
Solution errors with convergent precision.Finally, by choosing regularization parameter asα=2e-5, convergent precision as eps=1e-16, and other parameters as before, the true parameter vector can be reconstructed by the inversion algorithm given as (4.4)rinv=(2.00000,0.100000,0.500000,0.999998e-1,-0.500000,0.100000e-1,2.00000),
which is very close to the true solution and the solutions error is Errsolu=7.145e-8.


## 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 asr0=0 and the numerical differential step vector as τ=(1e-2,1e-2,1e-3,1e-5,1e-7,1e-9,1e-11)T. Table 3 and Figure 1 can show some convergence of the inversion algorithm. In Table 3, j denotes the number of iteration times, α is regularization parameter, and the solutions error is an average absolute error expressed by(4.1)Errsolu=‖rtrue-rinv‖2=(∑i=1717|ritrue-riinv|2)1/2,
where rinv denotes an inversion reconstruction solution corresponding to the true solution rtrue.Table 3
Solution errors with iteration numbers and regularization parameters.
jα=1e-3α=5e-4α=1e-4α=5e-5100.2603770.3073080.1722560.356620200.1467210.1211664.31409e-20.286733408.28802e-24.87719e-22.77992e-20.1978731003.71804e-21.04950e-21.50932e-26.41182e-22002.43731e-25.33821e-35.50282e-38.79467e-34001.91484e-23.53564e-37.33347e-41.60040e-46001.56244e-22.35879e-39.78162e-52.90766e-68001.27555e-21.57379e-31.30485e-55.39138e-810001.04143e-21.05009e-31.73988e-63.26469e-912008.50343e-37.00688e-42.31128e-73.11220e-9Figure 1
Solution errors with iteration numbers for different regularization parameters.By Table3 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 j=1000, the solution error is Errsolu=1.04143e-2 for α=1e-3, but it becomes Errsolu=1.73988e-6 for α=1e-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 aseps=1e-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 by(4.2)Errsolu(α)=15.4158α+0.0056,α∈[2e-5,3e-3].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α∈[2e-5,3e-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(4.3)‖rinv-rtrue‖2≤Eα,
where rinv represents the inversion solution, rtrue represents the true solution, and E 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 Section3.2, the inversion algorithm can be terminated if there is ∥ϵjα∥2≤eps, 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 α=2e-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.
eps1e-41e-61e-81e-101e-121e-141e-16Errsolu7.383e-27.206e-37.427e-47.294e-57.164e-67.382e-77.145e-8Figure 3
Solution errors with convergent precision.Finally, by choosing regularization parameter asα=2e-5, convergent precision as eps=1e-16, and other parameters as before, the true parameter vector can be reconstructed by the inversion algorithm given as (4.4)rinv=(2.00000,0.100000,0.500000,0.999998e-1,-0.500000,0.100000e-1,2.00000),
which is very close to the true solution and the solutions error is Errsolu=7.145e-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=5e-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:(5.1)rinv=(0.8649,1.753e-2,-7.103,19.77,4.863e-2,6.736e-4,1.475)T.
That is, the five coefficients in the dimensionless model (3.1)–(3.3) are given as(5.2)n=r1inv=0.8649,KF=r2inv=0.01753,(5.3)μ1=r3inv=-7.103,m=r7inv=1.475,(5.4)μ¯2(T)=19.77+0.04863T+0.0006736T2,
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 Z=1, which is plotted in Figure 4 compared with the actually measured breakthrough data.Figure 4
Reconstruction data and the measured breakthrough data.By Figure4, 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 as(5.5)Errinv=‖C(1,T;rinv)-Ĉ‖2=(∑k=1K1K|C(1,Tk;rinv)-Ĉk|2)1/2=0.008696,Errrel⁡=Errinv‖Ĉ‖2=0.001507,
respectively, where C(1,Tk;rinv) (k=1,2,…,K) denote reconstruction data by solving the forward problem with the inversion solution rinv, Ĉk (k=1,2,…,K) also denote the measured breakthrough data, and Ĉ=(Ĉ1,Ĉ2,…,ĈK)T.


## 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 oft=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 se=0.01753c0.8649. As for the source/sink coefficients, by the inversion results (5.4) and (5.3) and noting that μ¯2(T)=μ2(t) and T=vt/l, we can get the source coefficient in real-time dimension given as
(6.1)μ2(t)=19.77+0.01460t+0.00006095t2,
and the sink parameter is
(6.2)μ1=-7.102.
Since μ2(t)=γs,e-μs,e≫0 for t>0.5 (h) and μ1=γl-μl<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(t). By expression (6.1) in dimension form, we find that the coefficient of the third term of μ2(t) 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(t) takes higher-order polynomials, its expansion coefficients of higher-order terms also go to zero showing a numerical uniqueness of it.



---

*Source: 679531-2011-12-10.xml*