Solving Transient Groundwater Inverse Problems Using Space–Time Collocation Tre ﬀ tz Method

: This paper presents a space–time meshfree method for solving transient inverse problems in subsurface ﬂow. Based on the transient groundwater equation, we derived the Tre ﬀ tz basis functions utilizing the method of separation of variables. Due to the basis functions completely satisfying the equation to be solved, collocation points are placed on the space–time boundaries. Numerical solutions are approximated based on the superposition theorem. Accordingly, the initial and boundary conditions are both regarded as space–time boundary conditions. Forward and inverse examples are conducted to validate the proposed approach. Emphasis is placed on the two-dimensional boundary detection problem in which the nonlinearity is solved using the ﬁctitious time integration method. Results demonstrate that approximations with high accuracy are acquired in which the boundary data on the absent boundary may be e ﬃ ciently recovered even when inaccessible partial data are provided.


Introduction
To model a real groundwater system, the forward problem and its inverse must be solved. The forward problems of groundwater modeling are concerned with the determination of head distribution in a domain when the domain geometric shape, boundary conditions, initial conditions, and hydraulic parameters are specified [1][2][3]. On the other hand, the inverse problems involve the findings of the missing geometric shape, unknown physical parameters, and other conditions of the system by fitting observed system states. In the inverse problems, the boundary identification problem is of importance, in which the geometric shape of inaccessible parts of the domain remains unknown and is found from the measured groundwater head or fluxes on the accessible boundaries [4][5][6]. The hyporheic flow distinguished from groundwater flowing near rivers, which is important in understanding recharge or discharge, may belong to the boundary detection problem. The study of forward problems has developed rapidly. However, research on inverse problems is still limited to the consideration of simple models [7][8][9].
In recent years, significant developments in the meshfree methods have been reported [10,11]. The meshfree methods, which do not require mesh construction and are especially advantageous for problems including complicated and irregular geometry, have attracted much attention. The radial basis function collocation method, the Trefftz method, the method of fundamental solutions, and the reproducing kernel collocation method have been widely used [12][13][14]. The Trefftz method developed by Erich Trefftz is one of the relatively simple meshfree methods due to the collocation points provided only on the boundary since the basis functions are solutions of the Laplace equation [15][16][17][18]. The Trefftz method is usually utilized only for solving the Laplace-type equations [19,20]. Prior research shows that this approach may be restricted to stationary and linear problems. The engineering applications of the Trefftz method for solving inverse problems in subsurface flow may be a challenging task [21].

The Governing Equation
The transient groundwater problem in two dimensions is [6] ∂ 2 h(r, θ, t) ∂r 2 + ∂h(r, θ, t) r∂r + ∂ 2 h(r, θ, t) where h represents the total head, r is the radius, t is the time, θ represents the polar angle, D represents the hydraulic diffusivity, and Ω t represents the space-time domain. The initial data is where h 0 is the initial total head. The boundary data are where B D (r, θ, t) represents the Dirichlet boundary condition, and B N (r, θ, t) represents the Neumann boundary condition. In the case of the forward problem, one of the above boundary conditions is often specified. On the other hand, for the inverse problem, over-specified boundary conditions are applied at the given location. In this study, the accessible data contaminated by different random noise are considered where s denotes the percentage of noise, B D denotes the noisy Dirichlet data, B N denotes the noisy Neumann data, and rand denotes the random number.

Space-Time Trefftz Functions
The space-time meshfree approach adopting Trefftz functions originates from the Trefftz method. We assumed that Water 2020, 12, 3580 3 of 18 in which R(r), φ(θ), and T(t) are independent functions. The following notations are adopted.
Substituting Equation (7) into Equation (1) yields Dividing R(r)φ(θ)T(t) in Equation (9) obtains where λ and ξ denotes the separation constants. To guarantee λ and ξ to be negative or positive constants, the value p and q are considered and assumed to be λ = p 2 or λ = −p 2 or λ = 0 and ξ = q 2 or ξ = −q 2 or ξ = 0. There are nine possible scenarios with the combination of positive or negative or zero value for λ and ξ. When λ = −p 2 and ξ = −q 2 , the following equation is obtained where c 1 ,c 2 ,c 3 ,c 4 , and c 5 are arbitrary constants to be determined. Based on the superposition theorem, we may acquire wherer = r/R 0 , R 0 is the characteristic length,t denotes the dimensionless time, which is described ast = tD/R 2 0 , A 1 , A 2p , · · · C 6pq denote arbitrary constants that have to be evaluated, I 0 denotes the modified Bessel function of the first kind, K 0 denotes the modified Bessel function of the second kind, J 0 denotes the Bessel function of the first kind, Y 0 denotes the Bessel function of the second kind, I q denotes the modified Bessel function of the first kind of the q order, K q represents the modified Bessel function of the second kind of the q order, J q denotes the Bessel function of the first kind of the q order, Y q represents the Bessel function of the second kind of the q order, s and m denote the order of the basis functions. In this study, only positive basis functions are considered [6].

The Space-Time Collocation Scheme
The space-time collocation scheme is utilized to solve transient groundwater inverse problems in two dimensions as shown in Figure 1a. With the consideration of the space-time collocation scheme, the original Euclidean space domain, therefore, becomes a three-dimensional space-time, as depicted in Figure 1b. It is found that both initial and boundary data can be assigned on the space-time boundaries. Additionally, the space-time collocation scheme transforms the original transient groundwater inverse problem into an inverse BVP. The initial value problem used for solving inverse problem can be treated as the inverse boundary problem.  (15).
The space-time collocation scheme is utilized to solve transient groundwater inverse problems in two dimensions as shown in Figure 1a. With the consideration of the space-time collocation scheme, the original Euclidean space domain, therefore, becomes a three-dimensional space-time, as depicted in Figure 1b. It is found that both initial and boundary data can be assigned on the spacetime boundaries. Additionally, the space-time collocation scheme transforms the original transient groundwater inverse problem into an inverse BVP. The initial value problem used for solving inverse problem can be treated as the inverse boundary problem.

The Iterative Scheme for Modeling Boundary Detection Problem
The forward problems of groundwater modeling are concerned with the determination of head distribution in a domain. While the domain geometric shape, boundary conditions, initial conditions, and hydraulic parameters are specified, a linear system of simultaneous equations is yielded. The heads in the domain can be obtained using Equation (16) where the coefficients are solved by the linear system of simultaneous equations. On the other hand, the boundary identification problem is categorized into the inverse problem in which the missing geometric shape of the inaccessible part of the domain is the unknown. In this study, the missing geometric shape is described by a set of collocation points. The location ) , , ( of the collocation point is unknown, which has to be solved.
The coefficients are solved by the linear system of simultaneous equations from the given groundwater head or fluxes on the accessible boundaries.
Since the location ) , , ( of the collocation point is unknown, we may obtain the following equations. Equation (18)

The Iterative Scheme for Modeling Boundary Detection Problem
The forward problems of groundwater modeling are concerned with the determination of head distribution in a domain. While the domain geometric shape, boundary conditions, initial conditions, and hydraulic parameters are specified, a linear system of simultaneous equations is yielded. The heads in the domain can be obtained using Equation (16) where the coefficients are solved by the linear system of simultaneous equations. On the other hand, the boundary identification problem is categorized into the inverse problem in which the missing geometric shape of the inaccessible part of the domain is the unknown. In this study, the missing geometric shape is described by a set of collocation points. The location (r i , θ i , t i ) of the collocation point is unknown, which has to be solved. The coefficients are solved by the linear system of simultaneous equations from the given groundwater head or fluxes on the accessible boundaries.
Since the location (r i , θ i , t i ) of the collocation point is unknown, we may obtain the following equations.
Equation (18) is the nonlinear algebraic equation. Given (θ i , t i ), we can then solve r i using Equation (18). Most of the methods for solving nonlinear algebraic equations are based on the iteration scheme. Liu and Atluri (2008) proposed a time integration method named the Fictitious Time Integration Method (FTIM) [22,23]. The FTIM was first used to solve a nonlinear system of algebraic equations by introducing a fictitious time. The stationary point of these evolution equations is the solution for the original algebraic equation. Accordingly, we adopted the FTIM to solve Equation (18). The iterative scheme of the FTIM for solving the boundary detection problem is expressed as follows: where δ is the time step size, v is the non-zero coefficient, m is the value in the range of zero to one, and r i k is the value of r at the k th discrete time τ k = kδ.
The iterative scheme for the FTIM is utilized to deal with the system of nonlinear algebraic equations that are resulted from spatial discretization of the space-time collocation Trefftz method. It is obvious that the unknown in the iterative equation is the radius, which is the location of the boundary shape. Since the above equation is a nonlinear equation, the emphasis is placed on the two-dimensional boundary detection problem which the nonlinearity is solved using the FTIM. The iteration procedure begins by assigning a primary guess for the unknown spatial position and ends while the stopping criteria are reached.
where ε is the stopping criteria.

Accuracy and Convergence Analysis
A transient groundwater problem in two dimensions is conducted. The shape of irregular geometry is expressed as The initial data is h = sin(r cos θ) cos(r sin θ), (r, θ) ∈ Ω t .
The over-specified boundary conditions are applied using the following exact solution.
In this example, the hydraulic diffusivity is 10 −2 and the final elapsed time is 3 s. To validate our approach, a convergence analysis is firstly conducted. The following error is then evaluated where h E (r i , θ i , t i ) is the exact solution, h N (r i , θ i , t i ) is the numerical approximation, RMSE denotes the root mean square error, MAE denotes the maximum absolute error, and N I denotes the scattered inner points. Figure 2a shows the error versus the order of the basis functions. It appears that approximations with high accuracy are acquired while the order is superior to 14. The relationship between the absolute error and the number of boundary points is displayed in Figure 2b. Results obtained illustrate that

An Inverse Transient Groundwater Problem-A Cassini-like Boundary
The transient groundwater inverse problem is described as Equation (1). A two-dimensional Cassini-like boundary, as displayed in Figure 3, is defined as follows.

An Inverse Transient Groundwater Problem-A Cassini-like Boundary
The transient groundwater inverse problem is described as Equation (1). A two-dimensional Cassini-like boundary, as displayed in Figure 3, is defined as follows.
Water 2020, 12, x FOR PEER REVIEW 9 of 19 ( ) From Figure 3, not only the initial data are absent, but part of the boundary information is also not given. The specified boundary data are  From Figure 3, not only the initial data are absent, but part of the boundary information is also not given. The specified boundary data are h = e −2Dt sin(r cos θ) cos(r sin θ), (r, θ, t) ∈ ∂Ω given , The final time is 3 s, and hydraulic diffusivity is 0.1. The number of boundary points and the order of the basis functions are 1173 and 15, respectively. In addition, the parameters including δ = 0.8, v = 1.9, m = 0.1, and ε = 10 −12 are considered in this numerical implementation.
To obtain the recovered boundary data, 240 boundary points are assigned to the recovered boundary positions. To show the stability of our approach, the measured data are polluted by noise, where the noise levels are s = 0.001, s = 0.01, and s = 0.1. Figure 4 depicts the recovered boundary data by adding different noise levels. It seems that the recovered boundary data utilizing our method are consistent with the exact boundary data. However, Figure 4d shows that the discrepancy between numerical and analytical solution (at t = 0.2 s) with s = 0.1 may be found due to the large noise. Figure 5 displays the absolute error of our method at t = 0.2 s and t = 2.8 s. From Figure 5, the MAE (at t = 2.8 s) with s = 0, s = 0.001, s = 0.01, and s = 0.1 can reach the order of 10 −10 , 10 −4 , 10 −3 , and 10 −2 , respectively.

4.2.A Boundary Detection Problem-An Amoeba-like Boundary
The second example is to investigate the effectiveness and stability of the proposed method to solve this boundary detection problem. The transient groundwater inverse problem is described as Equation (1). We consider the given boundary positions defined by Equations (30)-(33) for four different cases in Figure 6.

A Boundary Detection Problem-An Amoeba-like Boundary
The second example is to investigate the effectiveness and stability of the proposed method to solve this boundary detection problem. The transient groundwater inverse problem is described as Equation (1). We consider the given boundary positions defined by Equations (30)-(33) for four different cases in Figure 6.
∂Ω given = (r, θ, t) r = e sin(θ) sin 2 (2θ) + e cos(θ) cos 2 (2θ), ∂Ω given = (r, θ, t) r = e sin(θ) sin 2 (2θ) + e cos(θ) cos 2 (2θ), ∂Ω given = (r, θ, t) r = e sin(θ) sin 2 (2θ) + e cos(θ) cos 2 (2θ), ∂Ω given = (r, θ, t) r = e sin(θ) sin 2 (2θ) + e cos(θ) cos 2 (2θ), Figure 6, not only the initial conditions are inaccessible, but partial boundaries are unspecified. Four different combinations of inaccessible space-time boundaries are conducted, as displayed in Figure 6. In case A, we consider a boundary detection problem, where 70% of the overall space-time boundary positions are given, as depicted in Figure 6a. In cases B, C, and D, not only the initial data but also parts of the space-time boundary positions are inaccessible. The percentage of the given boundary positions in cases B, C, and D is 50%, 30%, and 10%, respectively. As shown in Figure 6b-d, the specified boundary data are given at a partial portion of the space-time boundaries. The specified data are  From Figure 6, not only the initial conditions are inaccessible, but partial boundaries are unspecified. Four different combinations of inaccessible space-time boundaries are conducted, as displayed in Figure 6. In case A, we consider a boundary detection problem, where 70% of the overall space-time boundary positions are given, as depicted in Figure 6a. In cases B, C, and D, not only the initial data but also parts of the space-time boundary positions are inaccessible. The percentage of the given boundary positions in cases B, C, and D is 50%, 30%, and 10%, respectively. As shown in Figure 6b-d, the specified boundary data are given at a partial portion of the space-time boundaries. The specified data are h = (r cos(θ)) 2 + (r sin(θ)) 2 + 4Dt, (r, θ, t) ∈ ∂Ω given , h n = ∂h ∂n , (r, θ, t) ∈ ∂Ω given .
The boundary data on the unknown boundary positions are applied as h = (r i cos(θ i )) 2 + (r i sin(θ i )) 2 + 4Dt i , (r i , θ i , t i ) ∈ ∂Ω unknown (36) where (r i , θ i , t i ) are the positions of the unknown boundary. In this two-dimensional boundary detection problem, the final elapsed time is 3 s, and the hydraulic diffusivity is 0.1. The number of collocation points and the order of the basis functions to be 960 and 15, respectively. In addition, the parameters including δ = 0.8, v = 1.9, m = 0.1, and ε = 10 −12 are considered in this numerical implementation.
To recover the boundary positions at different times, 60 boundary collocation points are collocated. It seems that unspecified boundary positions can be completely recovered, as shown in Figure 7. To recover the boundary positions at different times, 60 boundary collocation points are collocated. It seems that unspecified boundary positions can be completely recovered, as shown in Figure 7.   Table 2 indicates the error with different noise levels. It appears that inaccessible boundary positions may be recovered even when noises are polluting the specified data.   Table 2 indicates the error with different noise levels. It appears that inaccessible boundary positions may be recovered even when noises are polluting the specified data.

A Boundary Detection Problem-A Flower-like Boundary
Three cases for the different combinations of unspecified boundary positions are conducted in this problem. For cases 1, 2, and 3 in Figure 9, the specified boundary locations are defined by Equations (37), (38), and (39), respectively. (c) From Equations (37)-(39), it is found that the accessible boundary data are only assigned at 1/8 portion of the space-time boundary.
The specified boundary data are assigned as  From Equations (37)-(39), it is found that the accessible boundary data are only assigned at 1/8 portion of the space-time boundary.
The specified boundary data are assigned as h = e −2Dt sin(r cos(θ)) sin(r sin(θ)), (r, θ, t) ∈ ∂Ω given , The boundary data on the unknown boundary positions are applied.
In this two-dimensional boundary detection problem, the final elapsed time is 3 s, and the hydraulic diffusivity is 0.1. The number of boundary points and the order of the basis functions are 960 and 15, respectively. In addition, the parameters including δ = 0.8, v = 1.9, m = 0.1, and ε = 10 −12 are considered in this numerical implementation.
To recover the boundary positions at different times, 240 boundary collocation points are collocated. It seems that the inaccessible boundary locations are completely recovered, as displayed in Figure 10. The absolute error of our approach for solving this boundary detection problem is also evaluated. Results obtained demonstrate that the MAE for case 1 to case 3 is 10 −3 , 10 −4 , and 10 −5 , respectively. It seems that the unknown boundary positions may be completely recovered, even 7/8 portion of the overall space-time boundaries are inaccessible.

Conclusions
We present a space-time meshfree method for solving transient groundwater inverse problems in this article. Key findings are described herein.
(1) For modeling of the two-dimensional transient groundwater inverse problems, an innovative concept, that the collocation points are placed on the space-time boundaries, is presented. Both the initial and boundary data are regarded as boundary conditions on the space-time boundary.

Conclusions
We present a space-time meshfree method for solving transient groundwater inverse problems in this article. Key findings are described herein.
(1) For modeling of the two-dimensional transient groundwater inverse problems, an innovative concept, that the collocation points are placed on the space-time boundaries, is presented. Both the initial and boundary data are regarded as boundary conditions on the space-time boundary. Using the space-time collocation, the initial value problem (IVP) for solving the transient groundwater equation is considered to be a problem of the inverse boundary value where the time marching for the IVP can be avoided. (2) Previous studies depict that the Trefftz method is limited to linear and stationary problems.
This study demonstrates that we solved transient groundwater inverse problems using the space-time meshfree method utilizing Trefftz functions with very high accuracy. The iterative scheme of the FTIM is utilized to deal with the system of nonlinear algebraic equations, which is resulted from spatial discretization of the space-time collocation Trefftz method. The advantages of the proposed meshfree method for solving inverse problems are also presented. Accurate approximations are yielded in which inaccessible boundary may be recovered, even 7/8 portion of the overall space-time boundaries are inaccessible.

Conflicts of Interest:
The authors declare that there is no conflict of interest.