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.
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 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]) 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]) where is Fredunlich’s adsorption coefficient and is called Freundlich’s exponent. By substituting expression (2.2) into (2.1), we can get where is called nonlinear retardation factor, and 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 , (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, , , and , , 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 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, (h). The values of the above-mentioned parameters are given in Table 1.
Before doing the experiment, the initial concentration of in the inflow was measured and denoted by [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 at the outflow are listed in Table 2.
By Table 2, we find that the solute concentration in the first outflow at (h) is (mg/L), which is double that of (mg/L) in the inflow, and the breakthrough data go down rapidly at the initial stage from (h) to (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 (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 given by (2.6) is a time-dependent function, that is, .
Now, let us give initial boundary value conditions for (2.3). Since the breakthrough data are available from the time of (h), we will take (h) as the initial moment. However, at the moment of (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 (mg/L), and (mg/L). According to the above information and by interpolation, we can get where 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
Obviously, if the retardation factor , Freundlich’s exponent , the adsorption coefficient , the source/sink coefficients and 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 and 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 , , and ; substituting them into (2.3), we have where , , and , , . Meantime, the initial condition (2.7) is transformed to and the boundary conditions are transformed to
The additional condition we will utilize is the measured breakthrough data listed in Table 2. Also by dimensionlessness to the real data, we have here 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 and 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 , and , and , we need a suitable approximate space to simulate 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 can be approximated by a quadratic polynomial given as follows: Then the model parameters we want to determine can be denoted by a vector where is Freundlich’s exponent, is the nonlinear adsorption coefficient, is the source/sink coefficient given by (2.5), and , , and are referred to as in (3.5), and is the initial index coming out in (3.2). Obviously, the vector belongs to Euclid space , and in concrete computations we will set a bounded ball as an admissible set for the unknown parameters, where is a positive constant.
For any prescribed (), set Thus, in order to get from the given , we only need to compute an optimal perturbation . 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 So, taking and () 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 : where the norm denotes ordinary Euclid norm, is regularization parameter, and is the additional data vector given by (3.4). Note that where () are basis functions of and () is numerical differential step. Furthermore, for given , denote and for , , denote Then the error function (3.9) can be transformed to the following form: 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 for , where the regularization parameter 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 , and convergent precision , and the additional measured data vector .
Step 2. Solve the forward problem (3.1)–(3.3) with Matlab to get and , and then obtain the vector and the matrix by formula (3.12),
Step 3. Choose suitable regularization parameter , and get an optimal perturbation vector by using formula (3.14).
Step 4. If there is , then the inversion algorithm can be terminated, and is taken as the coefficient solution, what we just want to determine; otherwise, go to Step 2 by replacing with .
4. Numerical Simulations
In order to verify numerical convergence of the inversion algorithm, several simulations will be presented by setting 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 , 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 and the numerical differential step vector as . 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 by where denotes an inversion reconstruction solution corresponding to the true solution .
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 , the solution error is for , but it becomes for . 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 . 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
By the computations, we find that numerical inversions should be performed for regularization parameters taking values in the interval of . 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 where represents the inversion solution, 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 , where 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 and initial iteration and differential steps as before, the solution errors are listed in Table 4 and plotted in Figure 3, respectively, where also denotes convergent precision, denotes solutions error defined by (4.1), abscissa represents logarithmic of , and longitudinal coordinates denote solutions errors in Figure 3.
Finally, by choosing regularization parameter as , convergent precision as , and other parameters as before, the true parameter vector can be reconstructed by the inversion algorithm given as which is very close to the true solution and the solutions error is .
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 and convergent precision as , 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: That is, the five coefficients in the dimensionless model (3.1)–(3.3) are given as 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 , which is plotted in Figure 4 compared with the actually 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 as respectively, where () denote reconstruction data by solving the forward problem with the inversion solution , () also denote the measured breakthrough data, and .
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 (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 . As for the source/sink coefficients, by the inversion results (5.4) and (5.3) and noting that and , we can get the source coefficient in real-time dimension given as and the sink parameter is Since for (h) and , 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 . By expression (6.1) in dimension form, we find that the coefficient of the third term of is 0.00006095, which is very small as compared with coefficients of the first two terms. Actually, by numerical computations we find that if takes higher-order polynomials, its expansion coefficients of higher-order terms also go to zero showing a numerical uniqueness of it.
Acknowledgment
This project is supported by the National Natural Science Foundation of China (nos. 11071148, 10926194, and 10471080).