# Data Reconstruction for a Disturbed Soil-Column Experiment Using an Optimal Perturbation Regularization Algorithm

Authors: Gongsheng Li; De Yao; Yongzai Wang
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

This paper deals with data reconstruction problem for a real disturbed soil-column experiment using an optimal perturbation regularization algorithm. A purpose of doing the experiment is to simulate and study transport behaviors of Ca2+, Na+, Mg2+, K+, SO42−, NO3−, HCO3−, and Cl− when they vertically penetrating through sandy soils. By data analysis to breakthrough data of the eight kinds of solute ions, two kinds of models describing their transport behaviors in the column are given. One is the advection-dispersion equation with time-dependent reaction terms suitable for three ions of HCO3−, NO3−, and K+, the other is the ordinary advection-dispersion equation suitable for the rest ions. Furthermore, all the unknowns in each model are determined by utilizing the optimal perturbation regularization algorithm, respectively, and then the breakthrough data for each considered ion are reconstructed successfully. The inversion results show that the advection-dispersion model with suitable time-dependent reaction terms can be utilized to describe the experimental process and reconstruct the experimental data.


---

## 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 for a solute transport process in the soils and groundwater (see [1–4], e.g.). To understand transport behaviors of the soil in the presence of infiltrating contaminants, soil column experiments are often performed in the laboratory. There are disturbed and undisturbed soil-column infiltrating experiments, and the soil-column in a disturbed experiment is often loaded uniformly and orderly with a lucite tube. The primary purpose of doing a disturbed soil-column experiment is to simulate and study transport behavior of the considered ion when it vertically penetrating through the column with artificial chemical liquid.On the research of soil-column experiments in the known literatures we have, 1D case is often concerned with. In early studies, it is a main task to construct mathematical models for soil-column experiments. The typical works could belong to Nielsen and van Genuchten and their researching groups. Nielsen et al. (see [5], e.g.) gave a general advection-dispersion equation describing solute transport in unsaturated soils, and van Genuchten and Wagenet (see [6], e.g.) discussed solute transport models of two sites/two regions in the soils. With development of computational tools and techniques, numerical methods and software packages based on convection-dispersion equations are widely utilized in the research of soil-column infiltrating experiments (see [7–13], e.g.). Recently, the authors have ever considered an undisturbed soil-column experiment using a linear model (see [14]) and a nonlinear model (see [15]), respectively, and gave successful numerical inversions by utilizing the optimal perturbation algorithm based on the measured breakthrough data.As we know, for a disturbed soil-column infiltrating experiment, there always has no more complicated physical/chemical reactions in the solute transport process, and then a mathematical model seems to be easy to establish due to linearity of the transportation. However, it is still a trouble on how to uniquely determine transport parameters in the model which can not be measured directly by the experiment. Another thing is how to work out the unknown parameters with high precision as much as possible. Thus, inverse problem and inversion analysis methods are often encountered with for solutes transport problems in soil-column infiltrating experiments.In this paper, we will investigate a disturbed soil-column experiment carried out in Zibo, Shandong Province, China, and study transport behaviors of eight kinds of ions ofCa2+, Na+, Mg2+, K+, SO42-, NO3-, HCO3-, and Cl- by the inversion method. On the basis of advection-dispersion mechanism and hydrochemical analysis to the experimental data, two kinds of models describing their transport behaviors in the column are discussed respectively. By applying the optimal perturbation regularization algorithm, all unknown transport parameters for each ion are determined, and then the measured breakthrough data are reconstructed successfully. The inversion results show that the reconstruction data basically coincide with the measured breakthrough data, and the advection dispersion equation with suitable time-dependent reaction terms can be utilized to describe the experimental process and reconstruct the experimental data. The paper is arranged as follows.Section2 gives the soil-column experiment, breakthrough curves of the eight ions and preliminary analysis of their transport behaviors. In Section 3, two kinds of transport models are given respectively, one is suitable for HCO3-, NO3-, and K+, and the other is for the rest. With the help of the models, the parameter identification problem for each ion is set forth based on the breakthrough data, and an optimal perturbation regularization algorithm is introduced. In Section 4, all the unknowns for each considered ion are determined by the inversion algorithm, and the breakthrough data are reconstructed. Finally, conclusions and discussions on the real problem here are presented.


## 2. Experimental Analysis

The disturbed soil-column experiment was carried out in a simple apparatus in a lab in Shandong University of Technology. The device is installed with three systems which are infiltrating system, soil-column system, and sample collector system, respectively. The soil column is composed by fine sands with diameter of 2 [mm], the diameter of the lucite tube loading the column is 18.6 [cm], and the height of the column is about 62 [cm]. After infiltrating the column with distilled water for 24 hours, the experiment was performed by infiltrating into the column with artificial coalmine water at the temperature of 20 centigrade.By collecting water samples at the bottom of the column, 23 samples were collected, and all of them were immediately sent to tested and analyzed in Shandong General Inspecting Station of Geology and Environment. Thus, we get the so-called breakthrough data of positive ions ofCa2+, Na+, K+, and Mg2+ and negative ions of Cl-, SO42-, HCO3-, and NO3-, which are plotted in Figure 1 and Figure 2, respectively, where the transverse axis denotes the experimental time [h], and the vertical axis represents concentration of the solute ion at the bottom of the column [mg/L].Breakthrough data of positive solute ions in the experiment.
(a)(b)(c)(d)Breakthrough data of negative solute ions in the experiment.
(a)(b)(c)(d)In Figure1, we find that Ca2+, Na+, and Mg2+ in the liquid phase almost have same behaviors which go up quickly at the initial stage and then tend to steady situations gradually after 2-3 hours; however, there is a different performance for K+. Although K+ goes up rapidly during the initial stage of 3 hours, it is still in an increasing trend at a slow speed after t=3 [h].On the other hand, let us investigate behaviors of negative ions plotted in Figure2. Cl- and SO42- have the same tendency with Ca2+, Na+, and Mg2+, but the situations are complicated for HCO3- and NO3-. It is noticeable that the measured initial concentration of HCO3- in the inflow is very small, and the first sample (at t=2/3 [h]) gives a high concentration which is 152.7 [mg/L]. The reason may be ascribed to the actions of free CO2 in the inflow. From t=2/3 [h] to t=2 [h], the concentration goes down rapidly, it goes up quickly from t=2 [h] to t=6 [h], and approaches to an equilibrium after t=6 [h]. As for NO3-, its first breakthrough concentration at t=2/3 [h] is 1.4 [mg/L], which is lower than that of 1.7 [mg/L] in the inflow, and after t=2/3 [h], it goes up rapidly during t∈(2/3,3/2) [h], and then goes down rapidly until t=3 [h]. It is surprising that it still goes down after t=3 [h] but with small oscillations.In the following, on the basis of the above analysis, two kinds of mathematical models describing the solutes transport behaviors in the column will be given respectively, and an optimal perturbation algorithm will be introduced to determine the unknown parameters in the models.


## 3. Mathematical Model and Inversion Algorithm

### 3.1. The Mathematical Models

It is wellknown that, under suitable hypotheses, solute transportation in homogeneous porous media can be characterized using 1D advection-dispersion equation. For the disturbed soil-column experiment discussed in Section2, physical/chemical reactions between the liquid and solid phases can be ignored since the column is composed by clean silver sands, and all of the ions transport in the column could be dominated by hydrologic dispersion and chemical reactions/ions exchanges among the solutes ions in the liquid phase. On the other hand, on the basis of the analysis in Section 2, we can deduce that transport behaviors of the ions of Ca2+, Na+, Mg2+, and Cl- with SO42- could be described by the ordinary advection-dispersion equation, and transport model for the ions of HCO3-, K+, and NO3- could be different, but can be described by the advection-dispersion equation with different time-dependent reaction terms representing possible chemical reactions/ions exchanges depending on the experimental time.Based on the above analysis, the ions ofCa2+, Na+, Mg2+, and Cl- with SO42- almost have similar behaviors resulting from convection and dispersion mechanisms. Without loss of generality, taking SO42- as example, and denoting c1=c1(x,t) as its concentration in the liquid phase, then there is the ordinary advection-dispersion equation given below:(3.1)∂c1∂t=aLv∂2c1∂x2-v∂c1∂x,0<x<l,t0<t<Ttol,
where aL is the longitudinal dispersivity, v is the average flow velocity, l is the length of the column, t0 is the time at which the first sample collected, and Ttol is the total infiltrating time.Next, due to specialties of concentrations ofHCO3-, K+, and NO3- varying with the experimental time, we will utilize time-dependent reaction terms to describe their transport behaviors, respectively. Denote ci=ci(x,t) (i=2, 3, 4) as their concentrations of HCO3-, K+, and NO3-, respectively, then we have(3.2)∂ci∂t=aLv∂2ci∂x2-v∂ci∂x+βi-1(t)ci,0<x<l,t0<t<Ttol,i=2,3,4,
where βi-1(t) (i=2, 3, 4) represent the physical/chemical reactions depending on the experimental time for solutes ions of HCO3-, K+, and NO3-, respectively, which are unknown in the experiment.Finally, following the method as used in paper [15], initial and boundary conditions of (3.1) and (3.2) are given as follows.(3.3)ci(x,t0)=ci0+(ci(l,t0)-ci0)xmilmi,0≤x≤l,i=1,2,3,4,
here mi>0 (i=1, 2, 3, 4) are unknown initial indexes which will be determined later. In addition, the boundary conditions are given as:(3.4)ci(0,t)=ci0,∂ci∂x(l,t)=0,t0≤t≤Ttol,i=1,2,3,4.Then, two kinds of mathematical models describing solutes transport behaviors in the soil-column are given by (3.1) and (3.2) together with the initial boundary value conditions (3.3) and (3.4), respectively. Since the initial indexes in condition (3.3) and the reaction coefficient functions in (3.2) are all unknown, a suitable inversion algorithm must be applied to determine them numerically. Furthermore, thanks to the independence of (3.1) with (3.2), we can see in the next section that each ion can be reconstructed independently of the others.In the following, we will takeHCO3-(c2) as example and introduce an optimal perturbation regularization algorithm (see [14–16]), e.g.) to determine the unknown parameters based on the measured breakthrough data. For convenience of computation, we will transform the model to a dimensionless form.SetC2=c2/c20, Z=x/l and T=vt/l, then we have(3.5)∂C2∂T=aLl∂2C2∂Z2-∂C2∂Z+lvq1(T)C2,0<Z<1,T0<T<T¯,
where q1(T)=β1((l/v)T), and T0=vt0/l, T¯=vTtol/l.The initial and boundary value conditions are transformed to(3.6)C2(Z,T0)=1+(c2(l,t0)c20-1)Zm2,0≤Z≤1,(3.7)C2(0,T)=1,∂C2∂Z(1,T)=0,T0≤T≤T¯,
respectively.Now the problem left is to determine the initial parameterm2 and the reaction coefficient q1(T) which can lead to an inverse problem of parameter identification. Therefore, some additional information from the experiment is needed. The additional condition we will utilize is the breakthrough data as plotted in Figure 2. Also by dimensionless to the real data, we have(3.8)C2(1,Tk)=C2k,k=1,2,…,K.As a result forHCO3-, an inverse problem of determining the time-dependent reaction function and the initial index is formulated by (3.5), the initial boundary condition (3.6)-(3.7), and the over-posed condition (3.8).In addition, the initial concentration ofSO42-(c1), HCO3- (c2), K+(c3), and NO3- (c4) are given as c10=1023.2 [mg/L], c20=34.1 [mg/L], c30=17.2 [mg/L], and c40=1.7 [mg/L], respectively. By the experiment, we know that the length of the column is l=0.62 [m], the average flow velocity is v=2.4×10-5 [m/s], and the total infiltrating time is Ttol=21 [h]. Finally, noting that the soil medium is composed by sand with 2 [mm] diameter, the longitudinal dispersivity aL=8.8×10-3 [m], which is estimated by known experimental results (see [1], e.g.).


### 3.2. The Optimal Perturbation Regularization Algorithm

In this subsection, we will also takeHCO3-(c2) as example to introduce the optimal perturbation regularization algorithm. Suppose that the unknown source parameter has an approximate polynomial expansion given as:(3.9)q1(T)≈∑i=0NaiTi,
where ai(i=0,1,…,N) are expansion coefficients. Then, the model parameters we are to determine can be denoted by a vector given as follows:(3.10)r=(a0,a1,…,aN,m2),
here ai(i=0,1,…,N) are coefficients referred to (3.9), and m2 is the initial index given in (3.6).Now, for any prescribedrj(j=0,1,…), set(3.11)rj+1=rj+δrj,j=0,1,….
Thus in order to get rj+1 from given rj, we need to compute an optimal increasing perturbation δrj. In the following, for convenience of writing, rj and δrj are abbreviated as r and δr, respectively.DenoteC2(Z,T;r) as the unique solution of the forward problem (3.5)–(3.7) for any prescribed vector r. Taking Taylor’s expansion for C2(1,T;r+δr) at r and ignoring higher order terms, we can get(3.12)C2(1,T;r+δr)≈C2(1,T;r)+∇TC2(1,T;r)⋅δr.
Furthermore, define a perturbation functional for δr as follows:(3.13)F(δr)=‖C2(1,Tk;r)+∇TC2(1,Tk;r)⋅δr-C2k‖2+α‖δr‖2,
where α>0 is regularization parameter, and(3.14)∇TC2(1,Tk;r)⋅δr≈∑i=0N+1C2(1,Tk;r+τi)-C2(1,Tk;r)τiδri,k=1,2,…,K,τ=(τ0,τ1,…,τN+1) is numerical differential steps vector.Denote(3.15)ξ=(C2(1,T1;r),C2(1,T2;r),…,C2(1,TK;r));η=(C21,C22,…,C2K),(3.16)gki=C2(1,Tk;r+τi)-C2(1,Tk;r)τi,G=(gki)K×(N+2).
Then, the functional expression (3.13) can be reduced to(3.17)F(δr)=‖Gδr-(η-ξ)‖2+α‖δr‖2.
It is easy to testify that minimizing (3.17) can be reduced to solve the following normal equation (see [17], e.g.)(3.18)GTGδr+αδr=GT(η-ξ).By the above discussions, we know that the best perturbationδr can be solved by the normal equation (3.18) via(3.19)δr=(αI+GTG)-1GT(η-ξ),
then the optimal coefficient vector r can be approximated by iterations (3.11) as long as the perturbation satisfying a given precision. This is the principal idea and procedure of the optimal perturbation regularization algorithm. The key points of performing the algorithm lie in appropriate choices of initial iteration r0, regularization parameter α, and differential steps vector τ=(τ0,τ1,…,τN+1).


## 3.1. The Mathematical Models

It is wellknown that, under suitable hypotheses, solute transportation in homogeneous porous media can be characterized using 1D advection-dispersion equation. For the disturbed soil-column experiment discussed in Section2, physical/chemical reactions between the liquid and solid phases can be ignored since the column is composed by clean silver sands, and all of the ions transport in the column could be dominated by hydrologic dispersion and chemical reactions/ions exchanges among the solutes ions in the liquid phase. On the other hand, on the basis of the analysis in Section 2, we can deduce that transport behaviors of the ions of Ca2+, Na+, Mg2+, and Cl- with SO42- could be described by the ordinary advection-dispersion equation, and transport model for the ions of HCO3-, K+, and NO3- could be different, but can be described by the advection-dispersion equation with different time-dependent reaction terms representing possible chemical reactions/ions exchanges depending on the experimental time.Based on the above analysis, the ions ofCa2+, Na+, Mg2+, and Cl- with SO42- almost have similar behaviors resulting from convection and dispersion mechanisms. Without loss of generality, taking SO42- as example, and denoting c1=c1(x,t) as its concentration in the liquid phase, then there is the ordinary advection-dispersion equation given below:(3.1)∂c1∂t=aLv∂2c1∂x2-v∂c1∂x,0<x<l,t0<t<Ttol,
where aL is the longitudinal dispersivity, v is the average flow velocity, l is the length of the column, t0 is the time at which the first sample collected, and Ttol is the total infiltrating time.Next, due to specialties of concentrations ofHCO3-, K+, and NO3- varying with the experimental time, we will utilize time-dependent reaction terms to describe their transport behaviors, respectively. Denote ci=ci(x,t) (i=2, 3, 4) as their concentrations of HCO3-, K+, and NO3-, respectively, then we have(3.2)∂ci∂t=aLv∂2ci∂x2-v∂ci∂x+βi-1(t)ci,0<x<l,t0<t<Ttol,i=2,3,4,
where βi-1(t) (i=2, 3, 4) represent the physical/chemical reactions depending on the experimental time for solutes ions of HCO3-, K+, and NO3-, respectively, which are unknown in the experiment.Finally, following the method as used in paper [15], initial and boundary conditions of (3.1) and (3.2) are given as follows.(3.3)ci(x,t0)=ci0+(ci(l,t0)-ci0)xmilmi,0≤x≤l,i=1,2,3,4,
here mi>0 (i=1, 2, 3, 4) are unknown initial indexes which will be determined later. In addition, the boundary conditions are given as:(3.4)ci(0,t)=ci0,∂ci∂x(l,t)=0,t0≤t≤Ttol,i=1,2,3,4.Then, two kinds of mathematical models describing solutes transport behaviors in the soil-column are given by (3.1) and (3.2) together with the initial boundary value conditions (3.3) and (3.4), respectively. Since the initial indexes in condition (3.3) and the reaction coefficient functions in (3.2) are all unknown, a suitable inversion algorithm must be applied to determine them numerically. Furthermore, thanks to the independence of (3.1) with (3.2), we can see in the next section that each ion can be reconstructed independently of the others.In the following, we will takeHCO3-(c2) as example and introduce an optimal perturbation regularization algorithm (see [14–16]), e.g.) to determine the unknown parameters based on the measured breakthrough data. For convenience of computation, we will transform the model to a dimensionless form.SetC2=c2/c20, Z=x/l and T=vt/l, then we have(3.5)∂C2∂T=aLl∂2C2∂Z2-∂C2∂Z+lvq1(T)C2,0<Z<1,T0<T<T¯,
where q1(T)=β1((l/v)T), and T0=vt0/l, T¯=vTtol/l.The initial and boundary value conditions are transformed to(3.6)C2(Z,T0)=1+(c2(l,t0)c20-1)Zm2,0≤Z≤1,(3.7)C2(0,T)=1,∂C2∂Z(1,T)=0,T0≤T≤T¯,
respectively.Now the problem left is to determine the initial parameterm2 and the reaction coefficient q1(T) which can lead to an inverse problem of parameter identification. Therefore, some additional information from the experiment is needed. The additional condition we will utilize is the breakthrough data as plotted in Figure 2. Also by dimensionless to the real data, we have(3.8)C2(1,Tk)=C2k,k=1,2,…,K.As a result forHCO3-, an inverse problem of determining the time-dependent reaction function and the initial index is formulated by (3.5), the initial boundary condition (3.6)-(3.7), and the over-posed condition (3.8).In addition, the initial concentration ofSO42-(c1), HCO3- (c2), K+(c3), and NO3- (c4) are given as c10=1023.2 [mg/L], c20=34.1 [mg/L], c30=17.2 [mg/L], and c40=1.7 [mg/L], respectively. By the experiment, we know that the length of the column is l=0.62 [m], the average flow velocity is v=2.4×10-5 [m/s], and the total infiltrating time is Ttol=21 [h]. Finally, noting that the soil medium is composed by sand with 2 [mm] diameter, the longitudinal dispersivity aL=8.8×10-3 [m], which is estimated by known experimental results (see [1], e.g.).


## 3.2. The Optimal Perturbation Regularization Algorithm

In this subsection, we will also takeHCO3-(c2) as example to introduce the optimal perturbation regularization algorithm. Suppose that the unknown source parameter has an approximate polynomial expansion given as:(3.9)q1(T)≈∑i=0NaiTi,
where ai(i=0,1,…,N) are expansion coefficients. Then, the model parameters we are to determine can be denoted by a vector given as follows:(3.10)r=(a0,a1,…,aN,m2),
here ai(i=0,1,…,N) are coefficients referred to (3.9), and m2 is the initial index given in (3.6).Now, for any prescribedrj(j=0,1,…), set(3.11)rj+1=rj+δrj,j=0,1,….
Thus in order to get rj+1 from given rj, we need to compute an optimal increasing perturbation δrj. In the following, for convenience of writing, rj and δrj are abbreviated as r and δr, respectively.DenoteC2(Z,T;r) as the unique solution of the forward problem (3.5)–(3.7) for any prescribed vector r. Taking Taylor’s expansion for C2(1,T;r+δr) at r and ignoring higher order terms, we can get(3.12)C2(1,T;r+δr)≈C2(1,T;r)+∇TC2(1,T;r)⋅δr.
Furthermore, define a perturbation functional for δr as follows:(3.13)F(δr)=‖C2(1,Tk;r)+∇TC2(1,Tk;r)⋅δr-C2k‖2+α‖δr‖2,
where α>0 is regularization parameter, and(3.14)∇TC2(1,Tk;r)⋅δr≈∑i=0N+1C2(1,Tk;r+τi)-C2(1,Tk;r)τiδri,k=1,2,…,K,τ=(τ0,τ1,…,τN+1) is numerical differential steps vector.Denote(3.15)ξ=(C2(1,T1;r),C2(1,T2;r),…,C2(1,TK;r));η=(C21,C22,…,C2K),(3.16)gki=C2(1,Tk;r+τi)-C2(1,Tk;r)τi,G=(gki)K×(N+2).
Then, the functional expression (3.13) can be reduced to(3.17)F(δr)=‖Gδr-(η-ξ)‖2+α‖δr‖2.
It is easy to testify that minimizing (3.17) can be reduced to solve the following normal equation (see [17], e.g.)(3.18)GTGδr+αδr=GT(η-ξ).By the above discussions, we know that the best perturbationδr can be solved by the normal equation (3.18) via(3.19)δr=(αI+GTG)-1GT(η-ξ),
then the optimal coefficient vector r can be approximated by iterations (3.11) as long as the perturbation satisfying a given precision. This is the principal idea and procedure of the optimal perturbation regularization algorithm. The key points of performing the algorithm lie in appropriate choices of initial iteration r0, regularization parameter α, and differential steps vector τ=(τ0,τ1,…,τN+1).


## 4. Inversion Results

Following the methods as used in paper [14, 15], we will implement the above optimal perturbation regularization algorithm to determine the unknown vector r given by (3.10) numerically, and then reconstruct the measured breakthrough data for the ions of HCO3-, SO42-, K+, and NO3-, respectively. The reaction coefficient can be determined in different approximate space, and one needs to choose a suitable approximate dimension N so as to reconstruct the measured data with high accuracy. For the four kinds of ions, we will realize the inversion with different dimension. All computations are performed in a PC of Dell Dimension 9200.


### 4.1. Data Reconstruction forHCO3-

SetN=7 and choose regularization parameter as α=0.08, initial iteration as r0=(0,0,0,0,0,0,1), and numerical differential steps vector as τ=(1e-1,1e-2,…,1e-9), then by 100 iterations which cost CPU time 295.4 [s], we obtain the inversion parameters given as follows:(4.1)rinv=(-1.7124,11.1543,-18.0604,5.9008,9.3842,-9.3487,3.1368,-0.3681,1.0727),
that is,
(4.2)q1(T)=-1.7124+11.1543T-18.0604T2+5.9008T3+9.3842T4-9.3487T5+3.1368T6-0.3681T7,(4.3)m2=1.0727,
and the relative inversion error is worked out as(4.4)Err=‖η-C2(1,T;rinv)‖2/‖η‖2=0.0113,
where η is the measured breakthrough data vector of HCO3- given in (3.15), and C2(1,T;rinv) is the reconstruction data. Furthermore, the inversion reconstruction data and the real measured breakthrough data are plotted in Figure 3.Figure 3
The breakthrough data and reconstruction data forHCO3-.


### 4.2. Data Reconstruction forSO42-

Let us consider the solute ion ofSO42-. As for (3.1) with corresponding initial boundary conditions, there is only one initial parameter m1 to be determined. Therefore, the optimal perturbation algorithm can be simplified as done in paper [14]. By utilizing the optimal perturbation algorithm without regularization term (choosing regularization parameter α=0), and taking initial iteration value as 1, numerical differential step as τ=1e-1, the initial index can be determined by 7 iterations costing CPU time 1.7 [s] which is m1=3.7588, and the relative inversion error is Err=0.0277. Figure 4 plots the measured breakthrough data and the reconstructed data of SO42-, respectively.Figure 4
The breakthrough data and reconstruction data forSO42-.


### 4.3. Data Reconstruction forK+andNO3-

As forK+, similarly as done for HCO3-, take N=3 here and choose regularization parameter as α=0.05, initial iteration as r0=(0,0,0,0,1), and numerical differential steps vector as τ=(1e-2,1e-3,1e-4,1e-5,1e-6), the unknown reaction coefficient, and the initial index for K+ are worked out by 31 iterations which costing CPU time 24.0 [s] given as follows:(4.5)q2(T)=-0.1413-0.1859T+0.1464T2-0.0305T3,m3=0.6558,
and the relative inversion error is Err=0.0141. Figure 5 plots the measured breakthrough data and the reconstructed data of K+, respectively.Figure 5
The breakthrough data and reconstruction data forK+.Now let us deal with data reconstruction forNO3- in which case the situation seems to be a little complicated. If still following the above method, there are no inversion results we need to appear. So, we perform the algorithm in a little different way as compared with the above computations by choosing a fixed initial index prescribed. Moreover, if setting an initial index in advance, then we only need to determine the time-dependent reaction coefficient.On the concrete computations, takingm4=3, N=9, regularization parameter as α=0.05, initial iteration as zero, and numerical differential steps vector as τ=(1e-1,1e-2,…,1e-10), an optimal reaction coefficient for NO3- is worked out by 500 iterations which costing CPU time 1636.7 [s] given as follows:(4.6)q3(T)=3.003-20.104T+36.064T2-10.092T3-21.755T4+6.761T5+17.050T6-14.783T7+4.523T8-0.491T9,
and the relative inversion error is Err=0.0285. The measured breakthrough data and the reconstruction data for NO3- are plotted in Figure 6, respectively.Figure 6
The breakthrough data and reconstruction data forNO3-.


## 4.1. Data Reconstruction forHCO3-

SetN=7 and choose regularization parameter as α=0.08, initial iteration as r0=(0,0,0,0,0,0,1), and numerical differential steps vector as τ=(1e-1,1e-2,…,1e-9), then by 100 iterations which cost CPU time 295.4 [s], we obtain the inversion parameters given as follows:(4.1)rinv=(-1.7124,11.1543,-18.0604,5.9008,9.3842,-9.3487,3.1368,-0.3681,1.0727),
that is,
(4.2)q1(T)=-1.7124+11.1543T-18.0604T2+5.9008T3+9.3842T4-9.3487T5+3.1368T6-0.3681T7,(4.3)m2=1.0727,
and the relative inversion error is worked out as(4.4)Err=‖η-C2(1,T;rinv)‖2/‖η‖2=0.0113,
where η is the measured breakthrough data vector of HCO3- given in (3.15), and C2(1,T;rinv) is the reconstruction data. Furthermore, the inversion reconstruction data and the real measured breakthrough data are plotted in Figure 3.Figure 3
The breakthrough data and reconstruction data forHCO3-.


## 4.2. Data Reconstruction forSO42-

Let us consider the solute ion ofSO42-. As for (3.1) with corresponding initial boundary conditions, there is only one initial parameter m1 to be determined. Therefore, the optimal perturbation algorithm can be simplified as done in paper [14]. By utilizing the optimal perturbation algorithm without regularization term (choosing regularization parameter α=0), and taking initial iteration value as 1, numerical differential step as τ=1e-1, the initial index can be determined by 7 iterations costing CPU time 1.7 [s] which is m1=3.7588, and the relative inversion error is Err=0.0277. Figure 4 plots the measured breakthrough data and the reconstructed data of SO42-, respectively.Figure 4
The breakthrough data and reconstruction data forSO42-.


## 4.3. Data Reconstruction forK+andNO3-

As forK+, similarly as done for HCO3-, take N=3 here and choose regularization parameter as α=0.05, initial iteration as r0=(0,0,0,0,1), and numerical differential steps vector as τ=(1e-2,1e-3,1e-4,1e-5,1e-6), the unknown reaction coefficient, and the initial index for K+ are worked out by 31 iterations which costing CPU time 24.0 [s] given as follows:(4.5)q2(T)=-0.1413-0.1859T+0.1464T2-0.0305T3,m3=0.6558,
and the relative inversion error is Err=0.0141. Figure 5 plots the measured breakthrough data and the reconstructed data of K+, respectively.Figure 5
The breakthrough data and reconstruction data forK+.Now let us deal with data reconstruction forNO3- in which case the situation seems to be a little complicated. If still following the above method, there are no inversion results we need to appear. So, we perform the algorithm in a little different way as compared with the above computations by choosing a fixed initial index prescribed. Moreover, if setting an initial index in advance, then we only need to determine the time-dependent reaction coefficient.On the concrete computations, takingm4=3, N=9, regularization parameter as α=0.05, initial iteration as zero, and numerical differential steps vector as τ=(1e-1,1e-2,…,1e-10), an optimal reaction coefficient for NO3- is worked out by 500 iterations which costing CPU time 1636.7 [s] given as follows:(4.6)q3(T)=3.003-20.104T+36.064T2-10.092T3-21.755T4+6.761T5+17.050T6-14.783T7+4.523T8-0.491T9,
and the relative inversion error is Err=0.0285. The measured breakthrough data and the reconstruction data for NO3- are plotted in Figure 6, respectively.Figure 6
The breakthrough data and reconstruction data forNO3-.


## 5. Discussions and Conclusions

### 5.1. Discussions

#### 5.1.1. On Reaction Coefficient

Firstly, let us investigate actions of the reaction coefficients. By the inversion result (4.2) for HCO3-, and noting  (3.2) (i=2) and its dimensionless form (3.5), we can get(5.1)β1(t)=-1.7124+1.5544t-0.3507t2+0.0160t3+0.00354t4-4.91E-4t5+2.30E-5t6-3.76E-7t7,
which is the reaction coefficient with dimensional time t for HCO3-. By (5.1) we can see that the expansion coefficients of β1(t) become smaller and smaller as N goes to large.Similarly, the reaction coefficientsβ2(t) and β3(t) in (3.2) (i=2,3) can be worked out which are given as(5.2)β2(t)=-0.1413-0.0259t+0.00284t2-8.26E-5t3,β3(t)=3.003-2.8016t+0.7004t2-0.0273t3-0.00820t4+3.55E-4t5+1.25E-4t6-1.51E-5t7+6.43E-7t8-9.73E-9t9,
respectively. Furthermore, βi(t) (i=1, 2, 3) are plotted in Figure 7 so as to study their changes with the experimental time, respectively.Figure 7
Reaction coefficients with time.From Figure7, we can find that the reaction to K+ almost has no change throughout the experiment; and the reaction to HCO3- goes up at the initial stage, and then keeps a little increase during t∈(5,15)[h]. As for the reaction to NO3-, it decreases during the initial stage, and then keeps unchanging during t∈(2,15)[h], and after t>15 [h], it becomes decreasing quickly.It is noticeable that the above arguments basically coincide with those hydrochemical analysis given in Section2 except for NO3- in the case of t>15 [h]. In addition, we also find that β2(t) and β3(t) both have similar asymptotic properties as β1(t), that is, the expansion coefficients of the reaction term go to zero as N becomes large.


#### 5.1.2. On Space Distribution at the First Outflow

Now, look at the following figure which plots space distributions at the first outflow time ofT=T0 for the four kinds of solute ions.From Figure8, we can see that at the first outflow time of t=2/3 [h], NO3-, SO42-, and K+ are in similar situations, and HCO3- has a special behavior.Figure 8
Space distributions atT=T0.ForNO3-, its concentration varies from 1 (t=0) to 1.4/1.7≈0.8235 (t=2/3 [h]), and goes down almost along with a straight line very slowly.ForSO42-, its concentration varies from 1 (t=0) to 10.0/1023.2≈0.0098 (t=2/3 [h]), and goes down with convexity but faster than NO3-.ForK+, its concentration varies from 1 (t=0) to 1.1/17.2≈0.0640 (t=2/3 [h]), and decreases also faster than NO3-, but with concavity.As forHCO3-, its concentration changes from 1 (t=0) to 152.7/34.1≈4.4780 (t=2/3 [h]), and goes up quickly with almost linearity.In summary, during this interval oft∈(0,2/3)[h], hydrodynamical advection dispersion actions may be play dominating roles in the solutes transport through the column. However, there are some chemical reactions occurring in HCO3- except for mechanical advection and dispersion. For example, the free CO2 in the inflow can produce HCO3- when penetrating through the column, and the amount of HCO3- reaches a maximum at the bottom of the column at t=2/3 [h].


### 5.2. Conclusions

We give several concluding remarks to end this paper.(1)
By the model and data reconstruction results, we can see that the main transport behaviors of solutes ions in the soil-column are dominated by hydrodynamical advection and dispersion mechanisms. There are no complicated physical/chemical reactions between solid and liquid phases, especially for the solutes ions ofSO42-, Cl-, Ca2+, Na+, and Mg2+, whose transport behaviors through the column can be described by purely advection-dispersion equation (3.1) with suitable initial conditions.(2)
As forHCO3-, K+, and NO3-, their transport behaviors in the column could be described by the advection-dispersion equation (3.2) with suitable time-dependent reaction terms and initial conditions. It is noticeable that K+ has some properties similar to those five ions; but for HCO3- and NO3-, they have some different properties as compared with the other ions maybe resulting from chemical reactions or ions exchanges in the liquid phase.(3)
As stated in Section5.1.2, at the initial interval of t∈(0,2/3) [h], it is still a trouble to describe the solute transport behaviors in mathematics due to lack of measured data. Maybe better experimental devices, and other mathematical tools are needed. In addition, noting that nonlinear processes could occur quickly during this short time, perhaps it is feasible and prospective that utilizing the model in a smaller time scale of second instead of hour during this interval.(4)
By performing the inversion algorithm, we find that a numerical uniqueness on the inverse problem can be obtained in the sense of ignoring higher-order terms. Let us take the inversion forβ2(t) as example. For the reason of seeking for regularity, we perform the algorithm utilizing the initial index of m3=0.6558 determined in Section 4.3, and keeping the same inversion error. Table 1 gives the inversion results for β2(t) taking on higher-order polynomials including the case of N=3.Table 1
Inversion results forβ2(t) in higher-order polynomials.
Nβ2(t)3(−0.1410, −0.0260, 2.9e-3, −8.3e-5)4(−0.1400, −0.0265, 2.9e-3, −8.6e-5, 4.2e-8)5(−0.1446, −0.0191, 1.3e-4, 3.1e-4, −2.3e-5, 4.9e-7)6(−0.1409, −0.0222, 8.3e-4, 2.7e-4, −2.6e-5, 8.3e-7, −8.6e-9)7(−0.1412, −0.0201, 7.7e-6, 2.6e-4, 8.3e-6, −3.9e-6, 2.4e-7, −4.5e-9)From the table, we can find that if taking the first term as an approximation, the reaction coefficient should beβ2(t)=-0.14; if taking the front two terms as an approximation, the optimal reaction coefficient should be given by β2(t)=-0.14-0.02t. Another thing we should mention is that for different truncated levels N=3, 4, 5, 6 and N=7, all of the inversion results preserve the same relative error expressed by Err=0.0141. From this point of view, the numerical inversion is of uniqueness too. As to theoretical uniqueness of the inverse problem here, we will deal with it in our sequent papers.


## 5.1. Discussions

### 5.1.1. On Reaction Coefficient

Firstly, let us investigate actions of the reaction coefficients. By the inversion result (4.2) for HCO3-, and noting  (3.2) (i=2) and its dimensionless form (3.5), we can get(5.1)β1(t)=-1.7124+1.5544t-0.3507t2+0.0160t3+0.00354t4-4.91E-4t5+2.30E-5t6-3.76E-7t7,
which is the reaction coefficient with dimensional time t for HCO3-. By (5.1) we can see that the expansion coefficients of β1(t) become smaller and smaller as N goes to large.Similarly, the reaction coefficientsβ2(t) and β3(t) in (3.2) (i=2,3) can be worked out which are given as(5.2)β2(t)=-0.1413-0.0259t+0.00284t2-8.26E-5t3,β3(t)=3.003-2.8016t+0.7004t2-0.0273t3-0.00820t4+3.55E-4t5+1.25E-4t6-1.51E-5t7+6.43E-7t8-9.73E-9t9,
respectively. Furthermore, βi(t) (i=1, 2, 3) are plotted in Figure 7 so as to study their changes with the experimental time, respectively.Figure 7
Reaction coefficients with time.From Figure7, we can find that the reaction to K+ almost has no change throughout the experiment; and the reaction to HCO3- goes up at the initial stage, and then keeps a little increase during t∈(5,15)[h]. As for the reaction to NO3-, it decreases during the initial stage, and then keeps unchanging during t∈(2,15)[h], and after t>15 [h], it becomes decreasing quickly.It is noticeable that the above arguments basically coincide with those hydrochemical analysis given in Section2 except for NO3- in the case of t>15 [h]. In addition, we also find that β2(t) and β3(t) both have similar asymptotic properties as β1(t), that is, the expansion coefficients of the reaction term go to zero as N becomes large.


### 5.1.2. On Space Distribution at the First Outflow

Now, look at the following figure which plots space distributions at the first outflow time ofT=T0 for the four kinds of solute ions.From Figure8, we can see that at the first outflow time of t=2/3 [h], NO3-, SO42-, and K+ are in similar situations, and HCO3- has a special behavior.Figure 8
Space distributions atT=T0.ForNO3-, its concentration varies from 1 (t=0) to 1.4/1.7≈0.8235 (t=2/3 [h]), and goes down almost along with a straight line very slowly.ForSO42-, its concentration varies from 1 (t=0) to 10.0/1023.2≈0.0098 (t=2/3 [h]), and goes down with convexity but faster than NO3-.ForK+, its concentration varies from 1 (t=0) to 1.1/17.2≈0.0640 (t=2/3 [h]), and decreases also faster than NO3-, but with concavity.As forHCO3-, its concentration changes from 1 (t=0) to 152.7/34.1≈4.4780 (t=2/3 [h]), and goes up quickly with almost linearity.In summary, during this interval oft∈(0,2/3)[h], hydrodynamical advection dispersion actions may be play dominating roles in the solutes transport through the column. However, there are some chemical reactions occurring in HCO3- except for mechanical advection and dispersion. For example, the free CO2 in the inflow can produce HCO3- when penetrating through the column, and the amount of HCO3- reaches a maximum at the bottom of the column at t=2/3 [h].


## 5.1.1. On Reaction Coefficient

Firstly, let us investigate actions of the reaction coefficients. By the inversion result (4.2) for HCO3-, and noting  (3.2) (i=2) and its dimensionless form (3.5), we can get(5.1)β1(t)=-1.7124+1.5544t-0.3507t2+0.0160t3+0.00354t4-4.91E-4t5+2.30E-5t6-3.76E-7t7,
which is the reaction coefficient with dimensional time t for HCO3-. By (5.1) we can see that the expansion coefficients of β1(t) become smaller and smaller as N goes to large.Similarly, the reaction coefficientsβ2(t) and β3(t) in (3.2) (i=2,3) can be worked out which are given as(5.2)β2(t)=-0.1413-0.0259t+0.00284t2-8.26E-5t3,β3(t)=3.003-2.8016t+0.7004t2-0.0273t3-0.00820t4+3.55E-4t5+1.25E-4t6-1.51E-5t7+6.43E-7t8-9.73E-9t9,
respectively. Furthermore, βi(t) (i=1, 2, 3) are plotted in Figure 7 so as to study their changes with the experimental time, respectively.Figure 7
Reaction coefficients with time.From Figure7, we can find that the reaction to K+ almost has no change throughout the experiment; and the reaction to HCO3- goes up at the initial stage, and then keeps a little increase during t∈(5,15)[h]. As for the reaction to NO3-, it decreases during the initial stage, and then keeps unchanging during t∈(2,15)[h], and after t>15 [h], it becomes decreasing quickly.It is noticeable that the above arguments basically coincide with those hydrochemical analysis given in Section2 except for NO3- in the case of t>15 [h]. In addition, we also find that β2(t) and β3(t) both have similar asymptotic properties as β1(t), that is, the expansion coefficients of the reaction term go to zero as N becomes large.


## 5.1.2. On Space Distribution at the First Outflow

Now, look at the following figure which plots space distributions at the first outflow time ofT=T0 for the four kinds of solute ions.From Figure8, we can see that at the first outflow time of t=2/3 [h], NO3-, SO42-, and K+ are in similar situations, and HCO3- has a special behavior.Figure 8
Space distributions atT=T0.ForNO3-, its concentration varies from 1 (t=0) to 1.4/1.7≈0.8235 (t=2/3 [h]), and goes down almost along with a straight line very slowly.ForSO42-, its concentration varies from 1 (t=0) to 10.0/1023.2≈0.0098 (t=2/3 [h]), and goes down with convexity but faster than NO3-.ForK+, its concentration varies from 1 (t=0) to 1.1/17.2≈0.0640 (t=2/3 [h]), and decreases also faster than NO3-, but with concavity.As forHCO3-, its concentration changes from 1 (t=0) to 152.7/34.1≈4.4780 (t=2/3 [h]), and goes up quickly with almost linearity.In summary, during this interval oft∈(0,2/3)[h], hydrodynamical advection dispersion actions may be play dominating roles in the solutes transport through the column. However, there are some chemical reactions occurring in HCO3- except for mechanical advection and dispersion. For example, the free CO2 in the inflow can produce HCO3- when penetrating through the column, and the amount of HCO3- reaches a maximum at the bottom of the column at t=2/3 [h].


## 5.2. Conclusions

We give several concluding remarks to end this paper.(1)
By the model and data reconstruction results, we can see that the main transport behaviors of solutes ions in the soil-column are dominated by hydrodynamical advection and dispersion mechanisms. There are no complicated physical/chemical reactions between solid and liquid phases, especially for the solutes ions ofSO42-, Cl-, Ca2+, Na+, and Mg2+, whose transport behaviors through the column can be described by purely advection-dispersion equation (3.1) with suitable initial conditions.(2)
As forHCO3-, K+, and NO3-, their transport behaviors in the column could be described by the advection-dispersion equation (3.2) with suitable time-dependent reaction terms and initial conditions. It is noticeable that K+ has some properties similar to those five ions; but for HCO3- and NO3-, they have some different properties as compared with the other ions maybe resulting from chemical reactions or ions exchanges in the liquid phase.(3)
As stated in Section5.1.2, at the initial interval of t∈(0,2/3) [h], it is still a trouble to describe the solute transport behaviors in mathematics due to lack of measured data. Maybe better experimental devices, and other mathematical tools are needed. In addition, noting that nonlinear processes could occur quickly during this short time, perhaps it is feasible and prospective that utilizing the model in a smaller time scale of second instead of hour during this interval.(4)
By performing the inversion algorithm, we find that a numerical uniqueness on the inverse problem can be obtained in the sense of ignoring higher-order terms. Let us take the inversion forβ2(t) as example. For the reason of seeking for regularity, we perform the algorithm utilizing the initial index of m3=0.6558 determined in Section 4.3, and keeping the same inversion error. Table 1 gives the inversion results for β2(t) taking on higher-order polynomials including the case of N=3.Table 1
Inversion results forβ2(t) in higher-order polynomials.
Nβ2(t)3(−0.1410, −0.0260, 2.9e-3, −8.3e-5)4(−0.1400, −0.0265, 2.9e-3, −8.6e-5, 4.2e-8)5(−0.1446, −0.0191, 1.3e-4, 3.1e-4, −2.3e-5, 4.9e-7)6(−0.1409, −0.0222, 8.3e-4, 2.7e-4, −2.6e-5, 8.3e-7, −8.6e-9)7(−0.1412, −0.0201, 7.7e-6, 2.6e-4, 8.3e-6, −3.9e-6, 2.4e-7, −4.5e-9)From the table, we can find that if taking the first term as an approximation, the reaction coefficient should beβ2(t)=-0.14; if taking the front two terms as an approximation, the optimal reaction coefficient should be given by β2(t)=-0.14-0.02t. Another thing we should mention is that for different truncated levels N=3, 4, 5, 6 and N=7, all of the inversion results preserve the same relative error expressed by Err=0.0141. From this point of view, the numerical inversion is of uniqueness too. As to theoretical uniqueness of the inverse problem here, we will deal with it in our sequent papers.



---

*Source: 732791-2012-02-26.xml*