Parameters Estimation in a Time-Fractiona Parabolic System of Porous Media

: The simultaneous estimation of coefﬁcients and the initial conditions for model fractional parabolic systems of porous media is reduced to the minimization of a least-squares cost functional. This inverse problem uses information about the pressures at a ﬁnite number of space time points. The Frechet gradient of the cost functional is derived. The application of the conjugate gradient method for numerical parameter estimation is also discussed. Computational results with noise-free and noisy data illustrate the efﬁciency and accuracy of the proposed algorithm


Introduction
Many industrial, environmental, and biological processes, related to solute transport in porous media and groundwater, see, e.g., [1][2][3], are modeled using a dispersion-reaction system of equations. The solving of such parabolic systems of equations of heterogeneous porous media usually requires the application of numerical methods.
Recently, there has been a sharp increase in research activities on modeling the flow of fluids in porous media with complex pore networks.
Many time-fractional diffusion equations and diffusion wave equations appear when the integer time derivatives are replaced by the fractional derivatives, and this technology is also used in porous media modeling. However, apart from such formal action, similar equations describe large classes of physical and chemical processes, such as anomalous diffusion, turbulent flow, chaotic dynamics, etc. A remarkable consequence of turbulence is the emergence of anomalous diffusion, see, e.g., [9]. We wrote the system of equations for flow in the triple continuum approach (4) in [8]. The first continuum describes a flow in the matrix of the porous media, the second continuum belongs to the network of a small highly connected fractured network, and the third continuum is related to the flow in a low-dimensional fractured network. We have the following system of equations for p = (p 1 , p 2 , p f ): where Ω ⊂ R d , Γ ∈ R d−1 and for the continuum index i = 1, 2, f , p denotes the pressure, k i is the permeability, k i ≥ k i0 > 0, i = 1, 2, k f ≥ 0 , and η is is the mass transfer term, which are proportional to the continuum permeabilities. The left time Caputo derivative [10][11][12][13] is used: for the function p i (x, t), i = 1, 2, f . All further results can easily be extended to the model for multi-continuum media: where i = 1, . . . , I and I is the number of contaminants.
The numerical method, based on the implicit finite difference approximation for the time approximation and the generalized multiscale finite element method for the spatial discretization is developed in [8] for solving the systems (1)-(4), k f = 0.
The one-dimensional case of the same model is studied in [14,15]. A priori estimates for the solutions of the first and third boundary value problems of the system (1)-(3), using energy inequalities, are obtained in [14]. In [15] is the constructed and analyzed positivity preserving numerical method for solving the second and third boundary value problems of the system (1)- (3). For the temporal discretization, the L1 formula on the non-uniform mesh is utilized while, for the discretization in space, a second-order monotone discretization on half-integer grid nodes is used.
The inverse problems, where it is required to determine some of the coefficients in the differential equations, initial values, boundary values, and some of the fractional orders α 1 , . . . , α I , α f by the observation data of solution p = (p 1 , . . . , p I , p f ), are important.
The inverse problem is indispensable for the precise and accurate modeling for analyzing many phenomena, such as heat-mass transfer, atmospheric pollution, and porous media, see, e.g., [5,9,10,16,17]. After relevant investigations of inverse problems, we can identify coefficients, boundary conditions, etc., to determine the equations themselves and can step into solving initial boundary value problems and initial value problems. Therefore, the inverse problem is the premise for studying the corresponding forward problem.
In the last decades, many research works indicated that the fractional differential equations are more relevant than those of classical models to describe non-Darcian flow or anomalous diffusion in some especial environment, in particular porous media, see, e.g., [18][19][20].
Actually, the aim of the current work is not only to establish a theory for direct problems related to time-fractional differential equations, but, first of all, to apply the theory to inverse problems. The inverse problems are highly various and, in this paper, we consider one inverse problem to illustrate how our framework for fractional derivative models operates. However, for inverse problems for the time-fractional diffusion equations, and especially for the systems of such equations, for example, to recover the initial data, source function, or diffusion coefficient, and so on, with some additional information, there are not so many works, see, e.g., [9,[21][22][23][24][25] and references therein. Even fewer are the results for simultaneously identifying more parameters in fractional order systems.
The rest of this paper is structured as follows. In the next section, we introduce some preliminaries concerning the fractional order calculus. We discuss the well-posedness of the direct (forward) problem in the one-dimensional case in Section 3. We also present some a priori estimates for the solution. Inverse problems are formulated in Section 4, while the quasi-solution of one of these problems is discussed in Section 5. The gradient of the cost functional is derived in the next section. The conjugate gradient method (CGM) is described in Section 7. The finite difference method (FDM) for the numerical solution of the inverse problem is presented in Section 8. The results from the numerical experiments are proposed and discussed in Section 9. The paper is finalized with some conclusions.

Preliminaries
In this section, we introduce useful notation and some formulas from fractional calculus without the precision of the smoothness of the functions, see, e.g., [12,13,26] for more details.
We consider another classical fractional derivative called the left Riemann-Liouville derivative: is the left Riemann-Liouville integral of order α. Next, the right Riemann-Liouville integral and derivative are defined: We also use the integration-by-parts formula: The relationship between the Caputo and Rieamann-Liouville derivatives is given as follows: is the right Caputo derivative. Further, we will use one-and two-parameter Mittag-Leffler functions [12,13]:

The Direct Problem
In this section, we introduce a model parabolic system (direct problem) with first or third boundary conditions and present two theorems that are useful for the following investigations.
Without the loss of generality, and for simplicity, in this paper, we develop our approach for the particular case of system (1)-(4), the triple continuum model in [27], in 1D geometry, as well as the models in [18][19][20]. Namely, in the rectangle Q T = {(x, t) : 0 ≤ x ≤ 1, 0 ≤ t ≤ T}, we consider the first boundary value problem: The existence and uniqueness of a weak solution to the initial boundary value problem (8)-(12) could be investigated following the scheme developed in [4]. Then, the stronger regularity of the weak solution of the direct problem (8)-(11) for more smooth input data can be established as well. Throughout the following, we assume that there exists a solution: of the problem (3), (4), and (8)- (11), where C m,n (Q T ) is the class of the functions that, together with their partial derivatives of order m with respect to x and order n with respect to t, are continuous in Q T = Ω × [0, T]. Therefore, throughout this paper we are interested in the classical solutions of (8)- (12), i.e., functions p = (p 1 , p 2 , p 3 ) whose derivatives ∂ α i p i ∂t α i , i = 1, 2, f , and ∂ 2 p i ∂x 2 , exist at all points in Q and satisfy (8)- (12) pointwise. The following results, obtained in [14], will be used in the quasi-solution and the analysis of the sensitivity problem in Sections 5 and 6.
Then, the initial-value problem (8)-(12) has a unique classical solution that satisfies the a priori estimate: where ρ 1 is a constant and E β , E β,β are Mittag-Leffler functions (7).
For the proof, see the conference paper (Theorem 1 in [14]) and Appendix A.
The proofs of Theorems 1 and 2 are based on the ideas in paper [28].

The Inverse Problems
In the fractional models, there are parameters, e.g., diffusion and convection, reaction coefficients, or boundary conditions, as well as initial conditions, and source terms that are difficult for observation purposes and have to be inferred indirectly from measurements. This gives rise to many different fractional differential equation inverse problems.
More precisely, we study the inverse problem (IP) of the identification of the permeabilityfluid viscosity constant coefficients k ≡ {k 1 , . . . , k I , k f }, the mass transfer coefficients and It can be rewritten in the operator form A(a) = g , where A : A → G is an injective operator and a ∈ A, g ∈ G is an Euclidean space of data g = {G i 11 , . . . , G i LM }. The inverse problem A(a) = g is ill-posed, i.e., its solution may not exist and/or its solution is nonunique and/or unstable to errors in measurements (15), see, e.g., [9,12,16,26,29].
Since the unknowns may not be uniquely identified due to the limited amount of inversion output, the inverse problem is often reformulated as an optimization problem with the minimizers considered as the approximate solution in some general sense. Here, the inverse problem is reduced to the minimization problem: where the functional J(a) characterizes the quadratic deviation of the model data from the experimental data, and we take it as follows: Of course, there are many other problems, which may be formulated for the determination of the coefficients, the right-hand side of the differential equations, or the initial conditions, for example: IP1 the reconstruction of the coefficients k 1 (x, t), . . . , k I (x, t) and k f (x, t) from the pressure measurements at several time moments 0 < t 1 < · · · < t m ≤ T; IP2 the reconstruction of the initial pressure ϕ i (x) from the final time observation IP3 the simultaneous reconstruction of the mass transfer coefficients η is (x, t) and the initial pressure ϕ i (x) at the above measurements; IP4 the simultaneous reconstruction of the permeability k 1 (x, t), . . . , k I (x, t) and k f (x, t) with the transfer coefficients η is (x, t).

Quasi-Solution of the Inverse Problem IP
This section is devoted to the formulation and existence of a quasi-solution to the IP. Further, we concentrate on the IP and, for more clarity, we will consider the triple continuum model (8)- (12), where d = 1, I = 2, k 1 , k 2 are positive constants, and k f is a non-negative constant, as well as η 12 i.e., we have three mass transfer coefficients. Additionally, we will take C 1 = C 2 = C f = 1.
Let us precise the admissible set: We rewrite the functional J(a) as follows: where δ(·) is the Dirac-delta function, reformulating the inverse problem as a minimization problem for the functional J(a), namely: Ifâ ∈ A , such that J(â) = 0, then the solution of the minimization problem (15) is called the exact solution of the inverse coefficient problem IP. Otherwise, the solution of the minimization problem is called the quasi-solution, see, e.g., [24,[27][28][29]. The proof follows the arguments as in (Chapter 6 in [24]), so we omit them.

The Gradient of the Cost Functional J(a)
In this section, we analyze the sensitivity problem and then derive the Frechet derivative of the functional (16). Now, we can introduce the matrices with the vector η = (e, b, c) and the following matrices: Next, if δa is an increment, we denote the deviation of solution p(p 1 , p 2 , p f ) by δp(x, t; δa) = p(x, t; a + δa) − p(x, t; a). Then, if it satisfies the following initial boundary value problem with accuracy up to terms of order o(|δa| 2 ), the sensitivity problem is as follows: Lemma 1. The mapping a → p(a) is Lipschitz continuous from A to H 1,0 (Q T ), i.e., for any a + δa ∈ A , the solution p ∈ H 1,0 (Q T ) to the problem (8)-(12) satisfies: Proof. The increment δp is a solution of the sensitivity problem (16). Next, we follow similar arguments as in Theorem 1 to obtain the estimate (18).
We have considered the first boundary value problem (8)- (12). In a similar way, we can derive the sensitivity problem when third boundary conditions (14) are posed. Then, the corresponding boundary conditions take the form: On the basis of Theorem 2, one can obtain an estimate analogously to (19). Now, we are in a position to formulate our main result.
Theorem 4. The gradient of the cost functional J(a) is given by the following: where • is the Hadamard product and the vector function R(x, t) is the solution of the adjoint problem (V = 3): for i = 1, 2, f , p 3 = p f , and δ(t − t m ) is the Dirac-delta function.
Proof. The proof is performed in two parts: Part 1.For the increment of the cost functional (17), we have: Part 2. We write the variational formulation of the sensitivity problem (18) using R(x, t) as a test function.
We organize the scalar product ·, · in Q T of the system (18) with the vector R(x, t) ≡ (r 1 (x, t), r 2 (x, t), r f (x, t)) and by its components, meaning we have: Now, we perform integration by parts (5), using the boundary and initial conditions in (18) and (21). For example, with respect to the first equation of (18), we have: Using the linearity of the scalar product, the product of the rest terms of (18) is written as follows: Collecting these expressions and using the adjoint problem (18), we obtain three equalities: Using (22), we find: (δp(x l , t m ; a)) 2 .
We have studied the IP at first boundary conditions (11). Similar results are valid for the IP when the boundary conditions (11) are replaced by the third boundary conditions (14).

Conjugate Gradient Method (CGM)
In this section, the CGM [4,24,29,30] will be described for the numerical identification of permeability-fluid viscosity coefficients k ≡ {k 1 , . . . , k I , k f }, the mass transfer coefficients η ≡ {η is , η i f , i, s = 1, . . . , I} , and the initial conditions: to the inverse problem IP. The following iterative process is used for the estimation of the triplet of vector functions a = (η, k, ϕ).
Step 1. Set an initial approximation vector a 0 and stopping parameter ε > 0. Suppose that we have a n . Next, we explain how to obtain the next approximation a n+1 .
Step 2. Check the stop condition : if |J(a n ) − J(a n−1 )| < ε and J(a n ) < J(a n−1 ), then a n is an approximate solution of the inverse problem IP . Otherwise, go to step 3.
Step 4. Solve the adjoint problem (21) with FDM (or FEM) and obtain the solution R(x, t).
Step 5. Determine the gradient of the cost functional J (a n ) with formula (20).
Step 6. Calculate the descent parameter α n = 2J(a n ) J (a n ) for the minimum errors gradient method.
Step 7. Calculate the next approximation: a n+1 = a n − α n J (a n ) and go to Step 2.

FDM Realization
In this section, we utilize the FDM in order to solve numerically the IP in the case of I = 2. We define a uniform partition of Ω = ω h × ω τ : The value of the pressure p i at the grid nodes is denoted by y n i,j , namely y n i,j = p i (x j , t n ), j = 0, 1, . . . , N x , n = 0, 1, . . . , N t , i = 1, 2, f . Similar notations are used for other functions e n = η 12 For the approximation of the Caputo fractional derivative of the function v n+1 , we apply the L1 formula [31][32][33] with accuracy O(τ 2−α i ): where The FDM approximation of the direct problem (8)- (12), with constant diffusion coefficients and time dependent functions η, regarding the assumptions in Section 5, is as follows: where G i n y n i,j := n ∑ s=1 ρ i n,s − ρ i n,s−1 y s i,j + ρ i n,0 y 0 i,j , n = 0, 1, . . . , N t − 1, Next, we approximate the adjoint problem (21). First, we apply the time inversion t = T − t and r i (x, t) = r i (x, T − t) = r i (x, t). Similarly, e( t) = e(T − t), b( t) = b(T − t), c( t) = c(T − t), and the terminal condition become the initial condition. For the fractional derivative, applying also the change of the variable λ = T − λ in the integral, we have: Further, from (6), taking into account that r i (x, T) = r i (x, 0) = 0, we obtain: Therefore, applying the L1 formula (23) for the approximation of the fractional derivative in the time-inverted adjoint problem, we obtain: Note that, for the values of b n+1 iν , we use the known values of b n+1 iν at the current iteration, iν . Similarly, the contribution S n+1 i,j is applied to the same measurement points (x l , t m ), l = 1, 2, . . . , L, m = 1, 2, . . . , M as in (15), where only the order is inverted in time. Finally, the numerical solution of (21) is obtained from r n i = r N t −n i , n = 0, 2, . . . , N t . The integrals in (20) are approximated by the second-order trapezoidal quadrature [34]: t) ,

Numerical Experiments
In this section, we verify the efficiency of the proposed numerical method. The numerical tests are performed as follows.
First, we choose functions η, ϕ, the values of k (called exact parameters), and the right-hand side g. Thus, the direct problem (8)- (12) is fully determined and we solve this problem numerically by (24) to find the solution y, which we refer to as the exact solution.
Then, we solve numerically the inverse problem with the initial guess a 0 , which is obtained from the exact parameters, adding perturbation with amplitude ρ 0 . Considering the observations, we take the points of measurements at ML number of inner grid nodes, and the values of the measurements are determined from the exact solution (discrete data) or by perturbing the exact solution with the perturbation amplitude ρ a (perturbed data).
The error of the pressure is estimated, comparing the solution p 1 , p 2 , p f , obtained by solving the inverse problem (for discrete or perturbed data) with the corresponding numerical solutions of the direct problem (24), computed for the exact a. The parametersk 1 , k 2 , k f , η 12 , η 1 f , and η 2 f , identified by the inverse problem, are compared with their exact values. We give the absolute (E a ), relative (E r ), and root mean square error (RMSE) errors, defined by the following: where k i is the exact value.
On Figures 1-3, we plot the recovered and exact functions η, ϕ and the solution p at the last time level for α 1 = α 2 = α f = 0.5. We observe a good fitting of the identified to the exact one.
In Table 1, we give the RMSE error for different values of the deviation ρ 0 and fractional order. As can be expected, more precise results are obtained for α 1 = α 2 = α 3 = 0.5 and a smaller deviation, but the influence of the value of ρ 0 is not so significant.     Example 2 (Perturbed data). Now, we consider another set of model parameters and functions: Numerical tests are performed for 36 points of measurements M = L = 6, the same as in Example 1, where N x = N t = 80, T = 1, ε = 10 −3 , ρ 0 = 0.15 in (25), ρ a = (ρ a,1 , ρ a,2 , ρ a, f ), and the following measurements: In Figure 4, we plot the exact solution and generated measurements for ρ a,1 = 0.035, ρ a,2 = 0.03, ρ a, f = 0.04, α 1 = α 2 = α f = 0.5. In Tables 2-6, we give errors of the numerical solution obtained by the inverse problem for different values of the fractional order. In Figures 5-10, we depict the corresponding exact and restored functions. In general, the numerical approach is efficient for recovering the diffusion coefficients, reaction term, initial functions, and solution p. We observe that more precise results are obtained for the moderate and equal values of the fractional orders α i , i = 1, 2, f . We may conclude that, for noisy data, the nonequal orders of the fractional derivatives affect the accuracy of the solution.
In Figures 11 and 12, we illustrate the changes of the errors in the maximal norm of the restored parameters and the functional and RMSE error at each iteration for α 1 For the considered set of measurements and deviations, the iterative process performs five iterations and stops since J 6 > J 5 .             Figure 11. Errors in maximal norm at each iteration of the restored parameters (η, k, φ), for i = 1 (line with circles), i = 2 (line with triangles), pi = f (line with squares), α 1 = α 2 = α f = 0.5, Example 2.

Summary and Future Work
In this paper, an inverse problem for the estimation of the coefficients and initial condition of a model time-fractional parabolic system in porous media was studied. The inverse problem uses information about the pressures at a finite number of points in the space-time domain.
This paper is devoted to the theoretical analysis and numerical implementation of the 1D inverse problem, which is reduced to the minimization of a least-squares cost functional.
The Frechet gradient of the functional was derived. Based on the explicit gradient formula, we described a CGM algorithm.
The numerical realization of the identification algorithm was proposed. The computational tests showed a good efficiency of the developed approach for simultaneously recovering the unknown diffusion coefficients, reaction terms, initial conditions, and the solution (pressure), even for noisy data. The different order of the fractional derivatives has a significant influence on the accuracy of the restored functions and solution. For α i < 0.25 or α i > 0.9, the algorithm fails to recover some of the parameters.
This present approach can be extended for multidimensional problems, but this will be undertaken in future work. Moreover, we plan to further improve the algorithm, applying the iterative regularization procedure of the CGM, such that it enhances the accuracy.   The embedding inequality is as follows: ∂p ∂x 2 0 and since η 12 ≥ 0, η 1 f ≥ 0, where ε = k 10 π 2 /l 2 , we obtain: Acting with the operator D −α 1 0t to both sides of this inequality, we find: In the same fashion, we obtain analogous estimates for p 1 and p f . Then, if we sum up all of these relations and apply the inequality [14], we have: and we derive: where M = M(m 12 , . . . , m f 2 , C 1 , C 2 , C f , α 1 , α 2 , α f , T) and N = N(l, k 10 , k 20 , α 1 , α 2 , α f , T) are constants. Next, neglecting the terms D −α i 0t ∂p i /∂x 2 0 ≥ 0, i = 1, 2, in the left-hand side of the inequality (A2), and applying Grönwall's inequality [28], we obtain (13).
Outline of the proof of Theorem 2. In the same manner, as in the proof of Theorem 1, we multiply Equation (8) by 2p 1 (x, t) and integrate the resulting expression with respect to x from 0 to l. Taking into account the boundary conditions (14), we derive: Further, if we take ε = k 10 b 1 +b 2 +2 , then there exists a constant M 1 > 0, such that we have: Acting on both sides of this inequality with the operator D −α 1 0t , and then applying inequality (A1), we derive: Similarly, we obtain the analogical inequalities for p 2 (x, t) and p f (x, t). Then, we have: where M 2 is a constant. By neglecting the third term of the left-hand side of inequality (A3), and applying Grönwall's inequality, we have: and we obtain the following inequality: where M 3 = Γ(β)E β,β (M 2 T β ). Finally, applying the inequality A1 for α := 2β and β := β, it follows from (A3) and (A4) that the a priori estimate holds.