Simultaneous use of partial differential equations in conjunction with data analysis has proven to be an efficient way to obtain the main parameters of various phenomena in different areas, such as medical, biological, and ecological. In the ecological field, the study of climate change (including global warming) over the past centuries requires estimating different gas concentrations in the atmosphere, mainly CO2. The mathematical model of gas trapping in deep polar ice (firns) has been derived in [12,15,16,17], consisting of a parabolic partial differential equation that is almost degenerate at one boundary extreme. In this paper, we consider all the coefficients to be constants, except the diffusion coefficient that is to be reconstructed. We present the theoretical aspects of existence, uniqueness and simulation for such direct problem and consequently formulate the inverse problem that attempts at recovering the diffusion coefficients using given generated data Antarctic and Greenland Polar snow and ice constitute a unique archive of past climates and atmospheres. Based on a good understanding of the mechanisms controlling gas trapping in deep polar ice, and therefore of the processes of densification and pore closure in Firns (typically over the first hundred meters of the polar cap), several models have been derived as a result of the collaborations between the ICE3 team of the IGE and GIPSA Lab (24 publications [12] including 3 in Nature [15,16,17]). Considering the mass conservation equations, the concentration  o  of a gas  in open pores satisfies an initial-value, time-dependent advection-diffusion partial differential equation on a one-space dimension segment [0, z F ] with Dirichlet boundary condition at 0 and a mixed one at z F , for z  (0, z F ), t > 0: * American University of Beirut (AUB), Beirut, Lebanon. ()  American University of Beirut (AUB), Beirut, Lebanon. ()  Grenoble-Alpes University, Grenoble, France. () with  atm  (0) = 0. Moreover, D  (z) is the effective diffusion coefficient of the gas  in the Firn (m 2 /yr) and is given by with z eddy , r  , and c f are known constants, and D eddy (z), D CO2, air (z) diffusion coefficients. The remaining terms are considered constants in this paper, and summarized in Then, the inverse problem is formulated by defining the objective function in section 3, where its gradient is computed using directional derivatives (section 3.1). Testings are performed on the efficiency of the computed gradient using MATLAB's FminUnc and FminCon functions and Nonlinear Conjugate Gradient method (section 3.2). Finally, concluding remarks are given in section 4. We start by stating the assumptions on the diffusion coefficient D  and the the concentration of gas in the atmosphere  atm  . 1. The Diffusion Coefficient D  (z) satisfies the following properties:  Lipschitz continuous on [0, z F ), specifically for all  > 0, there exists L  such that: 2. Moreover, the boundary condition at z F can be satisfied in one of the following two ways:  Full Degeneration D  (z F ) = 0. In this case, since z F is a singular point for 1 D  (z) , we assume uniform convergence of its integral on [0, z F ], specifically:  Quasi-degeneration D  (z F ) = o F << 1, in which case  o  satisfies the Robin condition In this case, equation ( 5) is automatically satisfied as z F is no longer a singular point to 1 D  (z) . In this paper, we will adopt full degeneration on D  (z), for theoretical considerations. While for numerical computations, we consider both options. 3. The gas concentration  atm  is assumed to be Lipschitz continuous on [0, T e ]. Given that 1/D 1/2   L 2 (0, z F ), i.e. let H 1 (0, z f ) be the usual Sobolev subspace. However, for the Firn problem, since D  degenerates at z F , we use the following subspaces of H 1 (0, z f ): v 2 Note that the injection of H 1  in H 1 is continuous with: Accordingly, we define Lemma 2.1. H 1 ,d is a closed subspace of H 1  and therefore itself a Hilbert space. Proof. Let {v n }  H 1 ,d be a converging sequence with v its limit point, and let {v  n }  L 2 (0, z F ) be a uniformly converging sequence. We need to show that v  H 1 ,d , i.e. v  L 2 (0, z F ), v   L 2 (0, z F ) and v(0) = 0. Since {v n }  H 1 ,d , then v n  L 2 (0, z F ), v  n  L 2 (0, z F ) and v n (0) = 0 for all n. Moreover, n . Thus v  L 2 (0, z F ) and v   L 2 (0, z F ). It remains to prove that v(0) = 0 where lim Integrate (10) with respect to z from 0 to z F , then using Cauchy-Schwarz inequality In the sequel, we prove a more general estimate in H 1  . Specifically, one obtains the following result. Lemma 2.2. Given that and therefore on the basis of arguments that lead to (10) and (11), one has: In what follows we denote  o  by . Let   H 1 ,d (0, z F ); multiplying the pde in (1) and using integration by parts with respect to z, in addition to the initial and boundary conditions, then equation ( 12) is reduced to ( 14) where then, (14) becomes To deal with the issue of existence and uniqueness of ( 17), we use Lions theorem ([4], ), which is stated below. Then, we apply it to our problem. Theorem 2.3. Let V and H be 2 Hilbert spaces satisfying: with the injection from V to H is dense and continuous.
Assuming a bilinear form a(, then for u 0  H and F (t)  L 2 (0, T ; V * ), the initial value problem admits a unique solution u, satisfying: Application of Lions Theorem to (17) To define the Hilbert spaces H and V , we first make a change of variable: Then (17) becomes: with: Then, to be in line with Theorem 2.3, we let: , which for every t is a linear form in . u 0 = 0 Then (17) can be stated as follows: Specifically, we let: Naturally, for the above definitions of H and V , we have with continuous injection from V into H. To prove existence and uniqueness to (25), we proceed with a sequence of preliminary results. Theorem 2.4. The bilinear form A satisfies the following: where C, C 0 and C 1 are positive constants independent of v and w. Let 1. We start by checking the bi-continuity of A(., . ). Given that: with This proves the bi-continuity of A(., . ). We turn now to the coercivity of A(., .) on This implies that for Using the geometric inequality: , for all o > 0, then and one obtains: Thus choosing o = 1  min{G f , f 1 } > 0 such that: validates the weak coercivity. Last point to prove is the existence of a function f * (t)  L 2 (0, T ; V ), such that: Using the bi-continuity of A(., .) and Cauchy-Schwarz inequality on the inner product < ., . >, where C = z Taking the limit as n goes to infinity implies lim n Theorem 2.6. There exists an f * (t)  V , such that: By the Riesz-Frechet representation and Lemma (2.5), there exists f * (t)  V such that t and   H 1  , Thus, by isometry, f * (t) We can now state our main existence and uniqueness result. Theorem 2.7. Proof. By applying Lions theorem to the Firn problem (25) with the subspaces H = L 2 (0, z F ) and V = H 1 ,d (0, z F ) that satisfy (26), and using theorems 2.4 and 2.6 we get the result of this theorem. We start first by rescaling (1) to the unit square, then by discretizing the obtained problem in time using Finite Difference Euler-Implicit scheme, followed by space discretization using Finite Element. The Firn extracted data is discrete on some large interval, typically 100 meters, with a distance of multiple meters between the different measurements. Thus, we rescale our PDE in space from the interval [0, z F ] to the interval [0, 1] to test our discretization numerically. Moreover, given that we would like to go back several hundred years to approximate the solution, we also rescale our time interval [0, T e ] to the interval [0, 1]. Let t = t/T e , z = z/z F ,  ( t, z) =  (t, z), and D(z) = D(z), where by the chain rule we have: Then system (1) becomes for  ( t, z), z  [0, 1], t  [0, 1]: Replacing  by , t by t  [0, 1], z by z  [0, 1], D by D, we get the following semi-variational form where f > 0, M  = M  g RT  [1.9496 * 10 5 , 6.4826 * 10 4 ], G =  +  > 0 and F = v + w air > 0 are constants. Let the bilinear form A is given by: where By integrating equation (17) over the temporal interval [t, t + t], with 0  t  1  t, one reaches the following L 2 Integral Formulation: For the full discretization of the Firn equation, the term t+t t A((z, s), (z))ds is first discretized using an implicit right rectangular rule: leading to the following fully implicit scheme in time. Let N = {z i | i = 1, 2, ..., n} be the set of nodes based on the partition of (0, z F ) with ., n  1} the resulting set of elements. The P 1 finite element subspace X n of H 1 (0, 1) is given by: For that purpose, we let B n = { i | i = 1, 2, ...n} be a finite element basis of functions with compact support in (0, 1), i.e.,: where We hence obtain the following fully implicit Computational Model. Given  n (t)  X n,d + { atm }, one seeks: where Thus, (43) is equivalent to Given  n,d (t)  X n,d , one seeks t > 0,   X n,d : Since the two discrete systems (43) and (44) are equivalent, studying the existence, uniqueness, and stability of one of them is equivalent to that of the other. To prove the property of existence on system (44) and the consequent properties of uniqueness and stability, one introduces the map: with u and w satisfying: Define on X n,d the bilinear form: To prove the existence of w from u, we investigate the positive definiteness of B(t, ., .) on X n,d xX n,d . For that purpose, using the weak coercivity of A, as given in (31), specifically: Define now Gf min = min{G f , f 1 }. Letting then o = Gf min  , yields: Letting, then we can state the following lemma Lemma 2.8. Under the assumptions of theorem 2.7, if 0 < t < o |K| < 1 for o << 1, one has: where: Proof. Let  = w in (46), then by (47) we get Theorem 2.9. Under the assumptions of Lemma 2.8 system (45) has a unique solution, i.e. w = L t (u) exists and is unique. Proof. For w = 0, B(t, w, w) > 0 is positive definite, thus system (45) has a unique solution.
We look now into the issue of stability for the discrete system (45), in which we let  = w, yielding: leading, using Lemma 2.8, to the inequality: Looking now at the term A( atm  (t + t) 1 , w), one has: Then, there exists a constant C A , independent of h, such that: Proof. Note that: As  1 (1) = 0, and writing (for simplicity)  atm :=  atm  (t + t) then (52) reduces to: Thus: where To handle the right hand sides, for the first term we use 3 and the Cauchy-Schwarz inequality, to get For the second term, we use h 0  1 = h 2 in addition to the second mean value theorem: Using the finite element inverse inequality (cite Ciarlet, , inequality (3.3.32)): one obtains: A similar argument leads to: Hence, one obtains: Factoring h 1/2 ||w|| 2 on the right hand side and letting: yields the result of the lemma. Hence using the result of this lemma and on the basis of inequality (51), one obtains, using Cauchy-Schwarz inequality: This leads us to the following stability theorem. Theorem 2.11. Assume  atm  (.) C[0, 1] is Lipschitz, i.e., there exists L atm , such that: then, one has: Proof. To reach the result of this theorem, one simplifies ||w|| 2 in (51), following it with a division by c o , then use of the previous lemma, in addition to the Lipschitz property of  atm  . As a consequence of this stability result, let { n,d } be the sequence generated by the fully discrete scheme (43), then one derives the following estimate. Corollary 2.12. Let C A,atm = C A + L atm /3. Then, under the conditions of the stability theorem 2.11, one has: Proof. Starting with (48), one has with  0,d = 0: 1. for K  0, c o = 1 and therefore: For o sufficiently small such that 1/c o  (1 + o), with  1, one has: From 1. and 2., follows (54). Given the fully discrete system (44), we start by extracting its matrix form along with the properties of the matrices and the obtained algorithm in section 2.3.1. Then, we test the algorithm's robustness, performance and accuracy in sections 2.3.2 and 2.3.3. Let  =  j for j = 2, .., n in (44) and define the vector where by noting that for j = 2, .., n: where M is the  ( atm  (t + t)   atm  (t))  1 ,  j 2 is equivalent to the vector of length n  1, where is equivalent to where v 2 (t+t) = F(z n , t+t)e n1 = B(t+t) is an (n1)x1 vector of zeros except the last entry, B is an is equivalent to the (n  1) x 1 vector T e v 3 (t) where Assuming a uniform mesh in space, i.e. Similarly the matrices S(D) and A(D) are of size (n  1) x (n  1), where D is a vector of length n with entries To compute these matrices we approximate their integrals in D using the Mean Value Theorem or Trapezoidal rule, i.e. Similarly, the constant c 1 in (61) is approximated by The explicit extraction of these five matrices assuming a uniform and nonuniform mesh is detailed in [18]. In this section we summarize the properties of these matrices. Lemma 2.13. The mass matrix M is a symmetric positive definite tridiagonal matrix. Proof. The symmetry of M follows from the symmetry of the L2 inner product (M i,j = M j,i ). The tridiagonal nature follows from the 1D finite element integrals. As for the positive definiteness, let v = 0 be a vector of length n  1, then Lemma 2.14. The tridiagonal matrix K and the diagonal matrix B are both positive semi-definite. Moreover, Lemma 2.15. Assuming that the vector D > 0 is strictly decreasing The linearity of the matrix A(D) follows from the linearity its entries, given approximately by (63). Proof. Let v = 0 be a vector of length n  1, then since all the terms are positive, and at least one entry v j = 0, implying that one of the terms is nonzero. The linearity of the matrix S(D) follows from the linearity its entries, given approximately by (62). Input: Mass matrix M ; Matrix C  as defined in 56; End time T e ; Time step dt; Mesh size h Discrete initial condition vector  0 ; c 1 as defined in (64); The function  atm (t). Output: : n x m matrix with the computed solution vectors  i for (2 : n, i + 1) =  i+1 ; 8: end for The existence and uniqueness of the solution to the discrete problem (55) is proven in the appendix. Thus, system (55) is solved iteratively given (0), as detailed in Algorithm 1 which is written in MATLAB syntax. At each time step a system of linear equations has to be solved using either direct methods, or iterative methods such as Krylov Subspace methods. In our implementation we solve the linear system using MATLAB's backslash operator, which is based on direct solvers. We consider different end times T e = 1, 50, 100, 150 and firn depths z F = 1, 50, 100, 150. As for the mesh sizes, we consider h = 1/16, 1/32, 1/64, 1/128, 1/256 and set the time step  to be either h or h 2 . For testing purposes, we set the constants: We start first by testing the robustness of algorithm 1 in section 2.3.2, and then report the runtime of the algorithm in section 2.3.3. It is important to have a numerically stable fast Direct Problem, as it will be called many times from the inverse problem.
We test the robustness of the algorithm by decreasing the mesh size h from h = 1/2 4 to h = 1/2 8 for all the possible combinations of depth z F and end times T e that we are considering, with the time step dt set to either h or h 2 . We run algorithm 1 for test cases 1, 2a, 2b, 2c, and 2d with inputs described above, with dt set to h 2 , for all the considered h values. Then, we compute the L infinity and L 2 absolute and relative errors between solutions at time T e and common space points corresponding to h = 1/2 4 . These errors are computed between solutions for different h values less than 1/2 8 and the solution for h = 1/2 8 . We run algorithm 1 with same settings described above for all the test cases, but with dt set to h. We obtain very similar results to Tables 2 and 3 where the corresponding errors for dt = h are slightly smaller with a difference in the sixth fractional digit in normalized format, i.e. a relative difference of order 10 6 . In Tables 2 and 3 we only show the results for T e = 150 since for T e = 1, 50, 100 the corresponding relative errors are almost identical up to the fifth fractional digit in normalized format. However, the solution evolves with time as shown in figures 1, 2, 3, 4, 5, 6 as the concentration increases with time. When comparing the obtained concentrations for different z F 's, it is clear from figures 4 and 6 that we obtain the same solution up to numerical errors. These errors are due to the fact that even though we are using the same rescaled h = 1/256, the actual spacing in the original scale is Thus, in figure 4, the actual h 50 = 50/256, whereas in figure 6 the actual h 150 = 150/256 = 3 * h 50 . So, if we set h = 1/3 * 256 for Case 1 with z F = 150 we get the same solution, as shown in figures 7 and 8. The advantage of using the rescaled problem is that the obtained discrete problem's runtime is independent of the choice of z F and T e , and solely depends on the choice of the mesh size h and the time step t = dt. Thus, in As expected, setting dt = h 2 will slow down the algorithm, as compared to dt = h. The reason is not due to a difference in the iterations' complexity, as it is clear that the runtimes per time step are of the same order. However, the smaller dt is, the more time steps are needed to reach the end Time, implying that more linear systems are solved in total, leading to a slower runtime. Thus, in what follows we will be setting dt = h as the algorithm approximates the solution in less runtime. Moreover, the accuracy of the solution is not affected much as shown in figures 9, 10, 11 and 12 , where we plot the Case 1 obtained   (z, T e ) for z F = T e = 150, and dt = h or h 2 , with h = 1/16, 1/64, 1/128 and h = 1/256 respectively. The relative error between the solution for dt = h 2 and the solution for dt = h is of order 10 7 for h = 1/16, of order 10 8 for h = 1/64, and of order 10 9 for h = 1/128, and 1/256. As for the mesh size h, it is obvious that a finer mesh leads to a smoother more accurate solution but at the expense of requiring more runtime, as dt = h will be smaller. On the other hand, if the mesh size h, is not small enough with respect to z F then we may get a numerically wrong solution as shown in However, noting that the solution varies most in the neighborhood of the rescaled z = 0 for z F ranging between 50 and 150, then another alternative is using an adaptive mesh size that is finer near z = 0, and coarser away from z = 0. We consider the following adaptive space mesh over the rescaled interval [0, 1] for a given h: Similarly to the case of uniform mesh, setting dt = h doesn't affect the solution as shown in figures 13,14,15, and 16, where we plot the Case 1 obtained   (z, T e ) for z F = T e = 150, with a nonuniform mesh for h = 1/4, 1/8, 1/16, 1/64 respectively and dt = h or h 2 . Moreover, the relative error between the solution for dt = h and that of dt = h 2 is of order 10 8 . Thus, using a nonuniform mesh with h = 1/16 or 1/32, Algorithm 1 provides a smooth and accurate solution comparable to that of a uniform mesh with h = 1/128 or 1/256. 13 (70) Thus, it is sufficient to find D CO2,air , which can then be used to obtain all other D  's by (70). Consequently, the inverse problem is to find D CO2,air by minimizing the objective function where g  (z) =  o ,g (z, 1) is the given rescaled data measurement, and   (z) =  o ,c (z, 1) is computed using the direct problem (55) for the corresponding approximation of D  (z). subject to P (  , d) = 0 and boundary conditions where the continuous differential operator of (33) is There are different type of methods for solving constrained and unconstrained minimization problems, that could be generally categorized as gradient-free methods and gradient methods. The gradient-free methods require only the objective function evaluation.
Whereas gradient methods for solving a minimization problem require at least the objective function and gradient evaluation, such as the steepest descent, and the nonlinear Conjugate gradient method. Moreover, there are other methods that require the gradient and the hessian, such as SQP method, interior point method, and Newton's method for solving V = 0. Note that a MATLAB implementation of SQP and interior point methods within the "fmincon" function in the optimization toolbox, allows solving the minimization problem without providing the gradient and the hessian, where the first is approximated using finite differences and the hessian using finite differences or a positive definite quasi-Newton approximation using the BFGS method. However, if the gradient expression or some approximation of it is available, then it could be passed to fmincon. The computation of the gradient and the hessian may be time consuming, if possible or available. Thus, we will be considering methods that will require the gradient evaluation, and the hessian is approximated, if needed. For example, the gradient-based methods, such as steepest descent or nonlinear CG methods start with an initial guess d 0 , compute d n+1 = d n +  where 's computation depends on that of the gradient of V. One option is to approximate the gradient using finite differences as in "fmincon". We also consider another option where we get an approximation starting from the directional derivative of V (d) along the direction , denoted by D  V . The directional derivative of V (  , d) with respect to d is given by where D  and DV denote the Frechet derivatives of   and V with respect to d, and denotes the directional derivative of   along (z). In what follows, we refer to v , by v  . To define v  , we start by differentiating the PDE P (  , d) = 0 along the direction (z) to obtain Thus, v  is the solution of the following for z  (0, 1) In section 3.1 we define the procedure that gives the gradient of the objective function using the directional derivatives and optimize this procedure in terms of runtime as much as possible. Then, we test the obtained Algorithm 3 using different constrained and unconstrained optimization methods in section 3.2. In section 3.1.1 we discretize the directional derivative (77) and then approximate V in section 3.1.2. To solve (77) we first find its equivalent semi-variational form, then use Euler-Implicit discretization in time and Finite Element discretization in space. Let T = {  H 1 (0, 1) | (0) = 0} then (77) is given in variational form by (78). Using integration by parts with respect to z, in addition to the initial and boundary conditions, and (0) = 0, then equation ( 78) is reduced to (79) Let the bilinear form and Note that B(., ) = A(., ) for T = T e and D  = r  d. Thus the time and space discretization of the first two terms is similar to that of (17) with the exception that the vectors with  atm are zero since v  (0, t) = 0. By integrating (82) over the temporal interval [t, t + t], with 0  t  1  t, one reaches the following L 2 Integral Formulation: For the full discretization of equation ( 83), the term t+t t B(v  (z, s), (z)) ds is first discretized using an implicit right rectangular rule: and the second term t+t t U (  (z, s), (z)) ds using the trapezoidal rule leading to the following fully implicit scheme in time. 84) Applying Finite Element in space, and using definition (42) and the fact that v  (z 1 , t) = 0 and that   (z 1 , t + t) +   (z 1 , t) =  atm (t + t) +  atm (t), then (84) simplifies as follows Let  =  j for j = 2, .., n in (85) and define the vectors Note that the left-hand side and the first term on the right-hand side of (86) are obtained exactly the same way as their corresponding terms in (55). As for the remaining 2 terms, they are obtained as follows. For j = 2, 3, is equivalent to 2 T e J  (t). is equivalent to the (n  1) x 1 vector 2 T e p  (t). Similarly to c 1 in (64), the constant c 2, in (87) is approximated by 3.1.2 Computing the gradient using the Directional Derivative v  Given the directional derivative of V (  , d) with respect to d is with the directional derivative of   along (z) where D  and DV denote the Frechet derivatives of   and V with respect to d, then it is possible compute the gradient entry by entry, considering  as the canonical basis {e 1 , e 2 ,    , e n }. At a first glance this might appear computationally intense as it requires solving problem (86) n times for computing the gradient once. However, since the matrix M + T e t C is independent from  and only the right-hand side vector T e t J  ((t+t)+(t))+T e t (p  (t)+p  (t+t)) is -dependent, then it is possible to solve the n v ,e j systems simultaneously by solving a linear system with multiple right-hand sides. Moreover, it is possible to reduce the time needed to compute the right-hand side vectors by noting that the matrices S(e j ) and A(e j ) have at most 7 nonzero entries.
Thus, J e j ((t + t) + (t)) can be computed without performing matrix-vector multiplication, nor generating the matrices S(e j ) and A(e j ) for j = 1 : n. For that purpose, note that the multiplication of A(e j ) with any vector v costs at most 6 flops However, given the repeating patterns in (89), the n matrix-vector multiplications cost 4(n  1) flops where a j = v j+1  v j for j = 0 : n  2; and b j = v j+2  v j for j = 0 : n  3; and v 0 = 0 . Similarly for S(e j ), its multiplication with any vector v costs at most 8 flops Given the repeating patterns in (90), the n matrix-vector multiplications cost 2(n  1) + n  2 flops where a j = v j+1  v j for j = 0 : n  2; and c j = 2v j+1  v j+2  v j = a j  a j+1 for j = 0 : n  3; and v 0 = 0. Thus, all the right-hand side vectors J e j v can be computed simultaneously using O(n) flops in each time iteration of (86), to obtain the (n  1) x n block matrix J , where the vector v = (t + t) + (t) is given Moreover, given that Thus, it is possible to find v ,e j (t) for j = 1 : n and t  (0, 1] simultaneously by solving where the (n1)xn block R consists of adding (93) and (91); i.e. adding c 2,e 1 ( atm (t)+ atm (t+t)) to the first 2 entries in the first row of matrix (91); and the At the end of the time iterations of (94), the (n  1) x n block V  (1) is obtained. Note that a zero first row could be added to V  (1) to get an n x n matrix due to the boundary condition v  (0, t) = 0. This procedure is summarized in Algorithm 2. The above procedure is repeated for each gas . Then, This gradient calculation is summarized in Algorithm 3. Note that since the matrix B  = M + T e t C  is fixed throughout the time iterations of both Algorithms 1 and 2 for the same D  , we perform a PLU decomposition of the matrix once, and pass the obtained triangular matrices to the Algorithms, instead of the matrix C  . Algorithm Compute the block matrix J as defined in (91), and let R = J ; 5: Call the Block v  Algorithm 2 for M, P  , L  , U  , to get the (n  1) x n matrix V  8: end for We first start by briefly discussing different constrained and unconstrained minimization methods to solve our inverse problem in section 3.2.1. Then, we test these methods using our gradient algorithm 3 or with finite difference approximations, and compare the obtained results in sections .. In all the testings, we generate data using the direct problem for 3 gazes with r  = [0.5, 1, 1.5] and D CO2 corresponding to test cases described in the Direct problem part. Moreover, we consider T e = 1, 50, 100, 150, z F = 1, 5, 10, dt = h g , and h g = 1/65 or 1/128; where a data set is generated for each combination. Then, for each data set,   (z, 1) values for h = 1/16, 1/32, 1/64 are generated using linear splines, as the inverse problem is solved using these meshings (h = 1/16, 1/32, 1/64). We start by briefly introducing the steepest descent, and nonlinear Conjugate Gradient. As mentioned earlier, both methods aim at minimizing V (d), by starting with an initial guess d 0 , then iterating with where  n is the descent's direction and  n is the step length along this direction. The computation of  n depends on that of the gradient of V. The steepest descent method considers which is the steepest descent from d n . Whereas, the nonlinear Conjugate Gradient method considers There are different choices for  n1 based on the already computed g n = V (d n ), g n1 ,  n1 , and y n1 = (g n  g n1 ), such as: , Fletcher and Reeves, 1964 , Polak and Ribiere, 1969 Once the descent's direction is set, it remains to compute the step length  n . In both methods,  n is obtained by solving the following line search where  should satisfy the Armijo condition and curvature condition, or what is known as the Wolfe condition. For a survey of different nonlinear CG methods, refer to [7]. In our testings we use nonlinear Conjugate Gradient method as implemented in the Poblano's MAT-LAB Toolbox [6] for solving the unconstrained version of the inverse problem, where the gradient is computed via Algorithm 3. A line search satisfying the strong Wolfe conditions is used to guarantee global convergence of the Poblano optimizers. Note that we modified the algorithm to include the "Hager and Zhang" nonlinear CG version. We have also tested nonlinear Conjugate Gradient method, where the gradient is approximated using Finite Differences as implemented in "" Toolbox [5]. Matlab's Optimization Toolbox has several functions that solve unconstrained and constrained optimization problems, such as fminUn and fminCon with different methods such as SQP (sequential quadratic programming) and Interior Point methods. As mentioned, Matlab's optimization functions require the gradient as input, and if not available a Finite Difference approximation is used. We test these methods on our problem where the gradient is computed via Algorithm 3 and compare the results with the finite difference approximation version. We solve the unconstrained inverse problem using nonlinear Conjugate Gradient (NCG) method as implemented in the Poblano's MATLAB Toolbox [6], where the stopping criteria are set to 10 8 . We also solve it using MATLAB fminUnc's Quasi-Newton method with tolerance 10 6 .
The data sets are generated using a mesh size h g = 1/65. Moreover, the results for MATLAB fminUnc's Quasi-Newton method are shown, where the gradient is computed using MATLAB's built-in Finite Differences or Algorithm 3. Both versions converge in number of iterations of the same order with similar relative errors. It is clear that the version with Algorithm 3 is faster, specifically 10 times faster for h = 1/32, 1/64. All the methods fail to converge for h = 1/64, with the exception of NCG, Case1, z F = 5, T e = 50. Moreover, the relative error is of the order of 10 1 when the methods converge. Figure 17 shows a sample of obtained D CO2 for both NCG and FminUnc. At z = 0 we get a negative value of D CO2 (0) = 50. One could set these negative values to zero, since D CO2  0. Another observation is that the obtained D CO2 is not smooth. Thus, it is possible to apply some polynomial regression to the obtained data, to get a smoother solution. We test Matlab fmincon's SQP and interior point (IP) methods for two constrained inverse problems, the first with the constraint that d  0 (72), and the second with the constraints that d  0 and d is decreasing (d i  d i+1 ). We solve these inverse problems using all the generated data sets with a mesh size h g = 1/65. Table 7 shows the runtime of the algorithms for solving the contrained problem (72) in seconds needed till convergence in (iter) iterations for stopping criteria 10 6 and the L2 relative error with the exact D CO2 . The obtained results for FminCon's SQP and IP methods for Case 2d are shown, where the gradient is computed using MATLAB's built-in Finite differences or Algorithm 3. Similarly to FminUnc, the corresponding versions with Algorithm 3 are 10 times faster with the same order of relative errors. Moreover, the IP method has a more stable behavior than SQP as it converged for all cases with relative errors of order 10 1 . The top four plots in Figure 18 show the corresponding solution for IP method where the gradient is computed using MATLAB's built-in Finite differences or Algorithm 3 for Cases 2d. The obtained solutions have a very similar profile with nonsmooth behavior that is reduced for a smaller mesh. By imposing the extra condition that d is decreasing, we get better results as shown in In this paper, we study the theoretical and computational aspects of the Firn direct problem. Then, the inverse problem is formulated by defining the objective function in section, where its gradient is computed using directional derivatives in a computationally efficient way. Testings on the efficiency of the computed gradient using MATLAB's FminUnc and FminCon functions and Nonlinear Conjugate Gradient method validate that it gives similar results to the case where the gradient is computed using Finite difference but in less runtime (up to 10 times faster). This work lays the ground for constructing a robust inverse problem algorithm that should extract past history diffusion coefficients of different gases of interest for understanding climate changes, and which will be the content of a forthcoming work . To prove that the Discrete Problem (55), with the approximated matrices A(D) and S(D), has a unique solution, we seek the invertibility of the matrix M + T e tC, looking for conditions that would make it positive definite. Let v = 0  R n1 . Starting with the identity: one obtains: Moreover, associated with the mass matrix M , the norm: which is equivalent to the weighted l 2 norm: as shown in the following Lemma. Lemma 4.1. The M -norm, . M , and the weighted l 2 norm, . h , satisfy the inequalities: Proof. From Lemma 2.13 one has: By Geometric inequality (ab ) and definition of weighted l 2 norm we get We prove now the following lemma. Proof. This result is straightforwardly obtained from the expressions of the matrices S(D  ) and A(D  ). We focus now on Term 1 . For that purpose, let ) and introduce v 0 = 0, then Term 1 can be rewritten as: We derive now two inequalities that can provide lower bounds for Term 1 . Lemma 4.3. Under the continuity and positivity assumption of D  (z), (3), there exists an h 0 < 1 (n 0 > 1), such that for h < h 0 (n > n 0 ) one has and a generalized Poincare's inequality where and c D a positive constant independent from h. Proof. Let v 0 = 0, then based on the identity , and using Schwarz inequality, one deduces that for i = 1, .., n  1 : Consequently, one has for , i = 1, .., n  1: where As the sequence S n converges to I(D  ), i.e lim h0(n) then there exists an h 0 < 1 (n 0 > 1), such that for h < h 0 (n > n 0 ) one has, where 0 < c D < 2. Then, for i = 1, .., n  1: and consequently: We now consider Term 2 . We prove the following result. Lemma 4.4. Under the continuity, positivity (3), and Lipschitz continuity (4) assumptions on D  (z), there exists an h 0 < 1 (n 0 > 1), and and consequently from estimate (102), Proof. Let o > 0, the there exists  := (o) such that given that lim zz f D(z) = 0 and D(.) is continuous. Thus, Term 2 can be expressed as Then, using the triangle inequality, (105), and Lipschitz continuity where Note that o could be chosen very small, such that O = {n  2}, i.e. D n2 < o < D n3 .
In this case we get, To prove that the Discrete Problem (55), with the approximated matrices A(D) and S(D), has a unique solution, we seek the invertibility of the matrix M + T e tC, looking for conditions that would make it positive definite. Let v = 0  R n1 . Starting with the identity: one obtains: Moreover, associated with the mass matrix M , the norm: which is equivalent to the weighted l 2 norm: as shown in the following Lemma. Lemma 4.1. The M -norm, . M , and the weighted l 2 norm, . h , satisfy the inequalities: Proof. From Lemma 2.13 one has: By Geometric inequality (ab ) and definition of weighted l 2 norm we get We prove now the following lemma. Proof. This result is straightforwardly obtained from the expressions of the matrices S(D  ) and A(D  ). We focus now on Term 1 . For that purpose, let ) and introduce v 0 = 0, then Term 1 can be rewritten as: We derive now two inequalities that can provide lower bounds for Term 1 . Lemma 4.3. Under the continuity and positivity assumption of D  (z), (3), there exists an h 0 < 1 (n 0 > 1), such that for h < h 0 (n > n 0 ) one has and a generalized Poincare's inequality where and c D a positive constant independent from h. Proof. Let v 0 = 0, then based on the identity , and using Schwarz inequality, one deduces that for i = 1, .., n  1 : Consequently, one has for , i = 1, .., n  1: where As the sequence S n converges to I(D  ), i.e lim h0(n) then there exists an h 0 < 1 (n 0 > 1), such that for h < h 0 (n > n 0 ) one has, where 0 < c D < 2. Then, for i = 1, .., n  1: and consequently: We now consider Term 2 . We prove the following result. Lemma 4.4. Under the continuity, positivity (3), and Lipschitz continuity (4) assumptions on D  (z), there exists an h 0 < 1 (n 0 > 1), and and consequently from estimate (102), Proof. Let o > 0, the there exists  := (o) such that given that lim zz f D(z) = 0 and D(.) is continuous. Thus, Term 2 can be expressed as Then, using the triangle inequality, (105), and Lipschitz continuity where Note that o could be chosen very small, such that O = {n  2}, i.e. D n2 < o < D n3 . In this case we get,