On several ill-posed and ill-conditioned mathematical problems of soil physics

Several well-known mathematical models of concentration fields in the soil (both at the single aggregate and the profile scales) are considered. It is shown that the respective boundary value problems for steady-state profiles belong to the class of ill-posed problems, since their solution does not exist. It occurs because a certain set of processes (for example, diffusion transport + first-order kinetic of the consumption) restricts possible boundary conditions, which, therefore, can no longer be arbitrary. Ill-posed inverse problems are also briefly described as well as one ill-conditioned inverse problem of parameters identification for mathematical model of the soil organic matter concentration profile. Exact solution for this model is the sum of two exponents. For a certain input data it was shown that this problem belongs to the class of ill-conditioned, since a small bias in the input data causes a significantly larger error in the solution (i.e. in calculated parameters).


Introduction
Methods of mathematical processing of the experimental results and mathematical modeling in ecology and soil science are currently used very widely and continue to develop [1][2][3][4]. In addition, the basic mathematical tasks in soil science are similar (or even identical) to the tasks of mathematical physics.
It is well known that the problems of mathematical physics can be divided into forward and inverse. Tasks which start with the causes and then calculate the effects represent the forward tasks. The inverse problems start from observed effects and then calculate causes not yet known [5]. For example, the soil properties and the intensity of organic matter input are the causes for the formation of soil organic matter (hereinafter -SOM) concentration profile. Here calculation of such a profile based on known soil properties and litter input represents a forward problem, whereas a calculation of soil properties or litter input, based on the profile of organic matter represents an inverse problem (hereinafter -IP). Besides the concept of forward and IP, classification of problems as well-posed and ill-posed has been introduced.
If sets of valid input data F and possible solutions V are given, the computational problem of determining v Î V according to f Î F from the equation Âv = f (where Â -a continuous operator acting from the solution space to the input data space) is considered as well-posed if: (i) its solution v exists; (ii) for each f there is a unique v; (iii) solution continuously depends on f i.e. solution is stable (small changes in input data f correspond to small variations of the resulting v). In this case, instability means that the solution is not a continuous function of data, i.e. a small perturbation of the data corresponds to a large perturbation of the solution [5].
A problem is ill-posed (hereinafter -IPoP), if one of the above conditions (i)-(iii) is not satisfied [5]. However, IPoP should not be considered as something wrong, as a manifestation of mathematical unprofessionalism. This is just a term for a special class of problems. Moreover, the ill-posed formulation of the problem makes it often possible to realize some important properties of the object under study. IPoP can arise among both forward and inverse problems. In addition, there is a large class of ill-conditioned problems (hereinafter -ICoP), which, although theoretically are well-posed, in fact, do not differ from IPoP (from the point of view of practical calculations).
Let us assume that the problem is correct (its solution exists, the solution is unique and stable to the uncertainty in the input data). Theoretically, a problem has stability if a sufficiently small error in input data induces a small error in the solution. However, in practice, the errors of the input data cannot be made arbitrarily small: their accuracy is limited. The sensitivity of solution to small errors in the input data is called the conditionality of a computational problem. A problem is called illconditioned if small input errors (δf) cause strong changes in the solution (δv). The value α in the inequality δv ≤ α·δf is called condition number (here δv and δf are relative errors). For ICoP α >> 1. However, an exact answer to the question about the value α, at which a problem should be considered as ICoP, essentially depends on the solution accuracy requirements and on the level of the accuracy provided for the source data [6]. In addition, for non-linear models, we should expect α ¹ const, therefore, the problem can be well-conditioned with some values of δf, but become ICoP with others (usually high) δf. Note that the problem that is ill-posed due to instability, is an extreme case of ICoP (with α = ¥). However, in practice, there is no difference between such ill-posed problem and an ICoP, which is formally well-posed, but having so high α ¹ ¥ that (for a given δf) δv turns out to be completely unsatisfactory. The minimum possible values of δf are limited both by the accuracy of measurements (available at a given level of technological development) and by the minimum error of the model (available at a given level of the scientific progress).
In this paper, we consider some examples of both ICoP and IPoP that arose in soil science practice.

Failed condition for the solution existence
Apparently, in literature more attention was paid to inverse IPoP -see, for example, [5][6][7][8][9], but in soil science forward IPoP also occurs. In this section, some examples of both types are considered.
(another conditionf(h) = Co with a large but finite h -is also mentioned in [10]; however, since the original solution is given only for equations (2), we will consider the boundary value problem with these BoC). The "solution" to this problem is suggested in the form where , So However, it is obvious that such function is not a solution, because , whereas according to equations (2) f(¥) = Co.

Non-existence of problem (1)-(2) solution.
One can assume that the mistake in the solution mentioned above could be easily fixed. However, the origin of the mistake is much deeper: the solution of problem (1)-(2) does not exist (so this problem belongs to the class IPoP). It could be shown as follows. One could write well established (see, for example, [11]) general solution of equation (1) in this form: , where , C1, C2 -arbitrary constants. This solution can be rewritten as: but if x®¥ f(x) cannot have arbitrary value Co, and inevitably must be zero: with a root-associated carbon source. Example described above is not unique. In soil science calculation of spatial concentration profiles is often an IPoP. Let us give another example, which will be considered in general form. Smagin [10] developed the model (1) assuming that root systems can serve as the source of SOM: (4) where the exponential term describes the intensity of the organic matter input from roots (distributed vertically in accordance with the exponential law). BoC (2) were chosen again. This problem also belongs to the IPoP class, since its solution does not exist. Indeed, according to [12], a general solution of inhomogeneous linear equation (4) can be presented as the sum of its particular solution and its general solution of the corresponding homogeneous equation. In this case equation (1) will be a homogeneous equation needed and its general solution is equation (3). It is easy to check that function is a particular solution for equation (4). Therefore, its general solution is: This again shows that when x®¥ the only possible boundary condition for function f(x) (i.e. solution) is zero and equations (2) could not been satisfied with arbitrary value Co. It also means that the solution of the problem (2), (4) for Co ¹ 0 does not exist.

Oxygen concentration profile in a soil aggregate.
At the end of this section, we briefly mention the model of biogenic O2 consumption inside a moist spherical soil aggregate considered in [13]: , where C is O2 concentration; r is distance from the center of the aggregate (0 ≤ r ≤ rm; rm corresponds to the surface of the aggregate); complete formulation of the problem in [13] includes equation (5) and the following BoC: , where Co > 0 is treated as the concentration of O2 at the surface of the aggregate (and depends on the depth of its location in the soil). The "solution" of the boundary value problem (5)-(6) is given by the author in the form But this function could not be a solution, since after applying L'Hopital's rule it is easy to show that while according to equations (6) C(0) should be 0. In general, it can be proved that the solution of the problem (5)-(6) does not exist: for equation (5) and BoC C(rm) = Co the value of C(0) is always positive. Therefore, it is impossible to satisfy BoC C(0) = 0 and the problem of solving equation (5) with conditions (6) is IPoP. But, as mentioned above, this fact should not be considered as negative. IPoP arises when the author relies on a set of assumptions and experimental facts: (a) inside the ( ) aggregate only diffusion transport occurs; (b) consumption of O2 satisfies the first-order kinetics; (c) in the center of aggregate the oxygen sensor reported zero concentration of O2... Occurrence of IPoP (considering forward problem) indicates the need to carry out comprehensive analysis of all assumptions and "facts" to find those that most likely contain uncertainties (they are not real). For example, from the above items (a) -(c) the most reliable is (a). But the fact that the kinetics could differ from the first order seems to be quite probable, and here it is necessary to investigate how the change in the kinetics type (within acceptable limits) will affect the problem. Further, the fact that the concentration of O2 in the center of the aggregate will be small (outside the sensitivity of the sensor), but, nevertheless, not zero, in general seems very likely. Fortunately, a similar model was built by Gontchar-Zaykin et al. [14]. They showed that a slight modification gives a well-posed problem: instead of zero boundary condition, it is enough to claim the limitations: 0 < C(0) < Co. Thus, the analysis of ill-posed problem leads to the possibility of a reasonable choice of further ways to study the object by experiment: using more accurate sensors could show whether the O2 concentration in the center of the aggregate will differ from zero.

Failed condition of the solution uniqueness
Completely different reasons can lead to multiple (non-unique) solutions of a certain problem. In a short manuscript we are not able to cover this topic comprehensively and just give a few typical examples.

Cluster analysis
3.1.1. The multiple options of distance selection in classification of objects. In soil science in general and in soil physics specifically cluster analysis is widely used for the classification of any objects [15,16,17,18]. The cluster analysis is designed to classify observations (objects) into more or less homogeneous groups. Its governing idea is to calculate a certain measure of similarity or dissimilarity (for example, some distance) between each pair of objects. Low value of this distance indicates that objects are similar or close to each other, while a large value indicates a lack of similarity. But similarity estimates could be different [19]. Moreover, other approaches are based on a non-distance concept. Specifically, an approach associated with fundamental concepts of the graph theory can be used [20]. For example, let us consider only MATLAB environment which is not designed namely for cluster analysis and operated only with the most common distances. Even in this case there are 10 different distance options (Chebychev distance, City Block metric, Hamming distance, Mahalanobis distance, Minkowski metric, Euclidean and Standardized Euclidean distances etc.). Moreover, for example, in the Minkowski metric parameter of the exponent (fixed number) must be set by user since there is an indefinite number of possible exponents. Therefore, there is an indefinite number of options for cluster analysis even if we use Minkowski metric only.

The multiple options of object grouping algorithms.
The next problem is to obtain a hierarchical grouping of objects, where the objects with the highest similarity coefficients are placed together. Then the groups of objects are bound to form new groups (they are most closely connected too), and this process continues until a complete classification of the objects is obtained. Unfortunately, at this stage there are many methods for analyzing groups as well [19]. Moreover, along with the approach of binding single objects into clusters, it is possible to divide the entire set of objects into clusters. The difference in procedures may be caused by other reasons, for example, by choosing initial points of clustering. This can be a random selection of one or several points as well as not random "typical" point [20]. Considering only MATLAB environment again, there are 7 available algorithms ("nearest neighbor", "furthest neighbor", "average linkage" etc.). It may seem that this situation is not unusual because it is typical for almost any problems. For example, in the classic numerical analysis we have many algorithms for integrating of differential equations (see, for example, [6,12,21,22]). However, this difference is fundamental. Various algorithms of the differential equations numerical integration give solutions with error less than a small value ε comparing with the exact solution of a differential equation. Therefore, all the algorithms give the same solution with an error that can be considered arbitrarily small (within reasonable limits) while various algorithms of grouping in cluster analysis can give completely different solutions.

Multi-optional choice of the classes' number.
Finally, one of the most important issues in solving a cluster problem is the choice of the required number of clusters [23]. In the literature hundreds (!) cluster analysis methods are described. Their diversity origins from the fact that the identification of clusters in many ways is "art". Many of researchers who use the techniques of this analysis create a new method themselves [20]. For example, in MATLAB there are two fundamentally different possibilities for this choice: "Finding the Natural Divisions in the Data Set" and "Specifying Arbitrary Clusters". Each of them generates a set of variants of the clusters number (in the latter case this number is just defined by user). Thus, at each stage of cluster analysis there is a variety of criteria and/or algorithms, that results in the non-uniqueness of the solution of the whole classification task.

Non-uniqueness of the model parameter identification problem.
When one identifies parameters of mathematical models, some distance between model prediction and experimental data is calculated [5,19,22,24]. Since this distance is a function of the model parameters, it is possible to solve the minimization problem: to find parameters for which the distance reaches minimal value. With these values of parameters the model fits the experimental data in the best way. Therefore, these values are considered as a solution of the parameters identifying problem. Obviously, non-uniqueness arises here from one of the reasons that we have already mentioned discussing non-uniqueness of cluster analysis. Indeed, in practice various distances are used: sum of squared deviations between model predictions and experimental data [5] (it corresponds to the sum of the absolute error squares); sum of weighted squared deviations [6,22,25]; sum of relative error squares [26]; some other types of distances are also used [19,27,28]. However, in practice nonuniqueness of solutions does not always lead to catastrophic consequences, which are usually associated with IPoP. Although we get formally different solutions (different sets of numerical values of parameters), but the values of the same parameter, (i) will not greatly differ from each other and (ii) assuming the adequacy of the model, will be close to the true values for models supplied with comprehensive set of experimental information. Therefore, in our opinion, the ill-posedness associated with the solution instability is much more dangerous, as well as the poor conditionality of the formally well-posed problem. Both occur when the model of the chosen level of complexity is not provided with adequate experimental data (for example if the complexity of the model exceeds the experimental data resolution).

Ill-conditionality of one inverse problem of the parameter identification
In soil science as well as in chemistry, biology, physics [20,22,29,30] a problem of representing experimental data as a sum of exponential functions (where coefficients have an important physical meaning) [10,31,32] is common.
Many models of SOM decomposition use only two exponential functions (one for labile components and one for recalcitrant): (7) where f is substrate (in % of the initial amount of substrate) at time x; C1 and C2 are initial amounts of substrates (%, therefore C1 + C2 = 100%) which are decomposed, respectively, with rate constants k1 and k2 [33]. If one needs to determine C1, C2, k1 and k2 from the residual substrate temporal dynamics, the solution of the problem can be presented (using the designation introduced in section "Introduction") as a vector ( ) ( ) ( ) The formula (7) is also a solution of the boundary value problem (2), (4) at Co = 0. In this case, the coefficients in equation (7) depend on physical parameters D, k, L, R, b: (8) Shape of the SOM depth profile represents SOM accumulation, decomposition and mass transfer processes rates. Therefore, it seems that it is theoretically possible to assess the intensity of these processes from the data on the SOM profile distribution [10]. But can it be used in practice?
So let us consider the problem of the SOM profile distribution model parameters identifying, i.e. we need to find the parameters C1, k1, C2 and k2 based on known SOM concentrations at various depths fi(xi) -and then calculate the values D, k, b and R: (9) As mentioned above, the problem of parameters identification is well studied and can be solved, for example, using the least squares method. However, we face the most difficult problem: how to find out that the correct results are obtained [34]. Apparently, the simplest solution is to test the properties of IP and its solution algorithm by following method: (i) set some values D, k, b, R and L; (ii) then calculate values fi(xi) by equation (7) for the depths xi where data on measured SOM concentration are available; (iii) then slightly "bias" the calculated values by adding random noise Δfi simulating the uncertainty of experimental data; and finally (iv) solve IP with "biased" values (fi(xi) + Δfi). The solution is the set of parameters that we denote D*, k*, R* and b* (since their values are slightly different from the original). Because we know the output values D, k, b and R exactly, the described algorithm gives an estimate of the error which arises from a certain IP solution method at given level of noise (if we perform the described algorithm N times). Therefore, it is possible to estimate the influence of the uncertainty in the "experimental" data on the error of the result, i.e. assess the conditionality of the problem or, in other words, its sensitivity to the source data uncertainty. According to [35] , where Np -number of model parameters, Ne -number of experimental points.
The idea underlying this sensitivity test is very simple and has a general character: to solve a problem with several datasets and to check how sensitive the solution is to a variation in input data [34]. Let us give a particular example. For calculations, we use values given in [10] for typical chernozem: D = 0.000251 m 2 year -1 , R = 0.227523 kg year -1 m -3 , k = 0.00278 year -1 , b = 2.558 m -1 and L = 0.008 kg year -1 m -2 . The values fi (given in table 1) were calculated according to the formula (7) with coefficients calculated by equations (8). These values (but with Gaussian noise Δfi) were used as "pseudo-experimental data" (further -PED). Several sets of PED for coefficient of variation CV = 5% are presented in table 1. For each ( ) ( ) were obtained using equations (9) and CVs for these parameters were evaluated (see MATLAB-code in Appendix).
The results are presented on a figure 1. Obviously, the identification problem is well-conditioned for parameter b because its value of CV approximately corresponds to CV of source data f. But for all other parameters we have typical ICoP. For example, if CV input data is 10% than CV for D will be 27.763·10 1.36 ≈ 636%, for R and k -about 1000% (k and R have almost identical curves of "CV -noise level" dependence, so on a figure 1 they are overlapping). In other words, if the experimental measurements (f) were obtained with an uncertainty of ≈10%, then the coefficients R and k obtained by solution of IP for model (7)-(9) could differ from the real ones by a factor of 10.
The reason of ill-conditional character of parameter identification is briefly illustrated on figure 2. As we see, the formula (7) with completely different values of the parameters (cf. C1 = -92.276 and C1 = -144.148) gives almost identical curves. Application of the method described here (conditional analysis through perturbation sensitivity) for other tasks in biology and soil science see, for example, in [36][37][38].

Basic principles of ill-posed problems solution
The question of the IP solution existence consists of two aspects. One is the physical existence of the set of parameters generating the observed data; the second is the existence of a mathematical solution of operator equation Âv = f (see "Introduction"). A formal mathematical solution may not exist. To understand it better, we note that the measured data fe always have some uncertainty (Δf): fe = f+Δf. The question is whether it is possible to find a set of parameters ve, that strictly generates a dataset fe: fe = Âve. The answer is that sometimes we cannot find such parameters, and it is easy to understand why. Indeed, fe -are dataset f spoiled by noise. But this noise is not related to model parameters because it could be produced by phenomena that were described by concentration field equations given above. This is a reason why noise could not be predicted using the same operator equation, as for idealized dataset. It means that we should not expect that we always can find physically realistic model that exactly corresponds to observed dataset [5]. If a model fitted data perfectly, the following equality would have place: ||Âve -fe|| = 0.
But does it make sense to find such a model? It does not give meaningful results for noisy data (data are always noisy). So we should think about certain approximate approach to inversion, based on the search for a model that would fit observations within the given accuracy. Therefore, we are IOP Conf. Series: Earth and Environmental Science 368 (2019) 012011 IOP Publishing doi:10.1088/1755-1315/368/1/012011 9 looking for the pseudosolution of inverse problem, i.e. about such a solution from the certain class of models that fits the observed data in the best way. Pseudosolution of IP exists if there is such ve that , (10) where ||•|| means a certain measurement of discrepancy between observed and predicted noisy data [5]. Different formulations of this measurement for variable types of ve and fe -vectors, functions, etc.see, for example, in [5,34,39,40]). Based on a using measure one can choose Δ = Δ(Δf) but strictly speaking Δ ~ Δf and often simple rule of thumb Δ = Δf is used. This simple idea (inequality (10)) -is a cornerstone of a regularization theory [5]. However, this inequality provides an existence of the solution but not its uniqueness. Indeed, a variety of different ve can satisfy this inequality: for some of them left-hand side would be much less than Δ, for others -just slightly less, and for a certain value of а ve right-and left-hand sides would be strictly equal. The correct solution of IPoP and ICoP can be obtained using a variety of regularization methods (deterministic, statistical or descriptive regularization) -see, for example, [5,7,39,41] and references therein. The most commonly used method of regularization of ill-posed problems is Tikhonov regularization ("ridge regression", "weight decay", "Tikhonov-Miller method", "Tikhonov-Phillips regularization").
Tikhonov regularization is based on minimization of smoothing functional , where ξ > 0 -regularization parameter, chosen based on equation (10); ||Âve -fe|| -functional of discrepancy between predicted and observed data; Ω -regularization term [40]. Both functionals should be constructed in a way to be nonnegative. Let us describe briefly one of possible options of Tikhonov regularization method (where Morozov principle for choosing the regularization parameter [5,42,43] was used, but there are also different principles of this choice -see, for example [42][43][44][45]). For an arbitrary ξо corresponding value of ve(ξо) is calculated, providing a minimal value of J. If for calculated ve(ξо) it was obtained that ||Âve(ξо) -fe|| < Δ, next iteration (of functional minimization) is carried out using the value ξ1 > ξо. If it was obtained that ||Âve(ξо) -fe|| > Δ, for the next iteration value ξ1 < ξо should be used. For the obtained new value ve(ξ1) next check of discrepancy between functional and uncertainty level should be performed. If ||Âve(ξ1) -fe|| ≠ Δ, using the procedure described above new value of ξ2 should be chosen. Then a new level of ve(ξ2), minimizing functional, should be found, and again discrepancy should be compared with uncertainty level until it would be found a value of ξ* providing ||Âve(ξ*) -fe|| = Δ. Corresponding value of ve(ξ*), providing a minimum of functional J, would be considered as a regularized solution.
It is obvious, that the scheme described above can be directly used in a case when solution of operator equation (v) is a function. But here we want to discuss a simpler example, when v is a vector (containing a certain number of components). A problem of parameter identification for equation (7) described in section 4 belongs to this class. In this case J is a function and not a functional. So instead of functional minimization relatively simple problem of minimization of multi-argument function should be considered. The term ||Âve -fe|| is a sum of squared deviations between model predictions and observed data. If there are any prior estimates of parameter values in equation (7), that could be considered as components of vector va, term Ω = ||ve -va|| can be calculated as a sum of squared deviations between components of a resulting vector and components of a priory [46]. If components of the vector vary considerably, normalization is usually carried out in Ω. For example, these squared deviations are divided into squares of corresponding va vector components [41]. A priori estimates are usually based on data for systems similar in certain characteristics to systems of interest. For example   10 if we need to find parameters for chernozem and parameters for luvisol are known, one can use them for calculations.

Concluding remarks
Apparently Chudnovskii [47] was one of the first scientists who raised the question about ill-posed inverse problem in soil science. He considered IPoP for soil thermal conductivity identification and methods for solving of this problem. Some inverse IPoP of ecology and soil science were solved in [48][49][50][51] etc. However, the ideology of IPoP in soil science unfortunately has not been developed comprehensively, despite the fact that IPoP is common in this area. Some researchers treated IPoP as a usual (i.e. in this case -well-posed) problem of numerical analysis. Of course, the results obtained with this approach, are ambiguous, since their uncertainty is high.
In this brief article we were able to consider only a few issues related to ill-posed and illconditioned problems, as well as only a very limited number of examples of such problems in the soil science. Thus, we absolutely did not discuss the classical problems of experimental data processing. But the number of IPoP is very high among them. Specifically, the problem of function differentiation (if the function has an uncertainty), probability density estimation for empirical distributions, solving some types of integral equations, etc. lead to IPoP. On the other hand, in recent years, the processes of automatic collection and processing of experimental data in soil science and ecology have become increasingly widespread [52][53][54][55]. So solving different ill-posed and ill-conditioned problems could become much more important in the near future.
Finally, it should be noted that due to the limited volume of the paper we briefly considered the most important topic: how to solve IPoP and ICoP. Fortunately, these methods are well developed in the mathematical literature.