On Solving Nonlinear Moving Boundary Problems with Heterogeneity Using the Collocation Meshless Method

Cheng-Yu Ku 1,2 , Jing-En Xiao 1,* and Chih-Yu Liu 1 1 Department of Harbor and River Engineering, National Taiwan Ocean University, Keelung 20224, Taiwan; chkst26@mail.ntou.edu.tw (C.-Y.K.); 20452003@email.ntou.edu.tw (C.-Y.L.) 2 Center of Excellence for Ocean Engineering, National Taiwan Ocean University, Keelung 20224, Taiwan * Correspondence: 20452002@email.ntou.edu.tw; Tel.: +886-2-2462-2192 (ext. 6159)


Introduction
The free surface flow is a moving boundary problem governed by the Laplace equation but has nonlinear boundary conditions.The study of free surface seepage problem plays a crucial role in the analysis of hydraulic engineering.In the design of embankment, earth dams and rock-fill dams, finding the position of the moving boundary is of importance [1,2].The solution of the Laplace governing equation may be carried out by solving the boundary value problem.Because the solution of the free surface seepage flow is nonlinear, iterative techniques are often required in the solution process for matching the over-specified boundary conditions.Various mesh-based numerical methods [3 -13] have been used for the analysis of free surface flow.For the mesh-based methods, automatic grid regeneration [6] is commonly used to solve the free surface problems in mesh-based approaches.However, convergence problems are often raised due to the changing of element shapes and types in the process of the mesh generation.While the complexity of the boundary conditions is considered, mesh-based methods may become unstable since the automatic grid regeneration is likely to generate distorted meshes.
Recently, meshless methods have attracted much attention to solve free surface seepage problems [14].Compared to mesh-based methods, the discretization of the domain for meshless methods is relatively simple because only arbitrary collocation points need to be placed in the physical domain without using any elements.If the basic functions satisfy the governing equation, the collocation points may be conducted only on the boundary.Accordingly, the meshless method has advantages with problems involving complex geometry [15].Among several meshless methods, the collocation Trefftz method (CTM) may be regarded as an attractive boundary-type meshless method [16].In the past, study of the Trefftz method was less widespread because the ill-posedness of the CTM limits the applications of the method.However, using the characteristic length [17][18][19], the CTM has been adopted to obtain accurate solutions for solving the Laplace governing equation in three-dimensions enclosed by simply and multiply connected domains [20].The Trefftz method is first proposed by the German mathematician Erich Trefftz.Later, the CTM [21][22][23][24][25][26][27][28][29][30] is commonly used for solving partial differential equations.Since the CTM is categorized into the boundary-type meshless method, it approximates the solutions of the governing equation using the Trefftz basis functions where the solutions are described as the assembly of the Trefftz functions [31].The CTM requires the evaluation of the coefficients in which they may be obtained by solving the linear simultaneous equations assembled by using the boundary conditions at a number of collocation points.Applications of the CTM such as Laplace and modified Helmholtz equations [32,33] and the problem of boundary detection [34] has been studied.Due to the complexity, applications of the CTM are most limited to the homogeneous problems.In addition, the study of nonlinear moving boundary problems with heterogeneity using the CTM has not been reported yet.
This paper presents the study on solving nonlinear moving boundary problems in heterogeneous geological media using the CTM.The free surface flow is a moving boundary problem governed by Laplace equation but has nonlinear boundary conditions.We adopt the CTM to approximate the solution using Trefftz base functions satisfying the Laplace equation.An iterative scheme in conjunction with the CTM for finding the phreatic line with over-specified nonlinear moving boundary conditions is developed.To deal with flow in the layered heterogeneous soil, the domain decomposition method is used so that the hydraulic conductivity in each subdomain can be different.The method proposed in this study is verified by several numerical examples.The formulation of the proposed method is described as follows.

Governing Equation and Boundary Conditions
The two-dimensional Laplace equation used to represent flow through a homogenous rectangular dam is expressed as ∆ϕ = 0 in Ω and where ϕ is the head, ∆ is the Laplacian, Ω represents the domain boundary of the problem, g and f denote the Dirichlet and Neumann boundary conditions, respectively.n denotes the normal vector.Γ D and Γ N denote the Dirichlet and Neumann boundary conditions.
The boundary conditions of the rectangular dam with a moving boundary, as depicted in Figure 1a, can be presented by Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 .The Dirichlet boundary conditions are imposed on the Γ 2 and Γ 5 , respectively.
where H 2 denotes the downstream elevation, H 1 denotes the upstream elevation.Neglecting the velocity head, the head is expressed as Water 2019, 11, 835 where Y(x) denotes the height above the sea level, p denotes the pore water pressure and γ denotes the water unit weight.On Γ 4 , the over-specified moving boundary conditions are given as On Γ 3 , the seepage face boundary condition is depicted as ϕ = Y(x) is unknown which can be solved by using the iterative scheme.On Γ 1 , the no-flow condition is given as Water 2019, 11, x FOR PEER REVIEW 3 of 21 where ) (x Y denotes the height above the sea level, p denotes the pore water pressure and  denotes the water unit weight.On 4  , the over-specified moving boundary conditions are given as On 3  , the seepage face boundary condition is depicted as is unknown which can be solved by using the iterative scheme.On 1  , the no-flow condition is given as

The Collocation Trefftz Method
The Laplace equation in polar coordinate system is expressed as where the radial coordinate is denoted by r and the angular coordinate is denoted by  .The solution of the Laplace governing equation is approximated by using the Trefftz basis functions satisfying the governing equation, as shown in Equation (10).The Trefftz basis functions are obtained by finding the general solutions using the separation of variables method [35].The Trefftz basis functions can be found to solve problems in a simply connected domain, as shown in Figure 2.

Seepage domain
Free surface

Separation point
Seepage face

The Collocation Trefftz Method
The Laplace equation in polar coordinate system is expressed as where the radial coordinate is denoted by r and the angular coordinate is denoted by θ.The solution of the Laplace governing equation is approximated by using the Trefftz basis functions satisfying the governing equation, as shown in Equation (10).The Trefftz basis functions are obtained by finding the general solutions using the separation of variables method [35].The Trefftz basis functions can be found to solve problems in a simply connected domain, as shown in Figure 2.

Formulation of T-Complete Basis Functions
We may apply the separation of variables [36].The solution may be in the following form: For simplicity, we let Then, Equation (10) can be rewritten as follows.
We divide on both sides in the above equation and the equation can be rewritten as two differential equations as follows.
Using the constant v to ensure positive or negative constants, we have    .Substituting Equations ( 16) and (17) into Equation (11), we have where 0 a and 0 b denote the coefficients.Considering the second scenario,

Formulation of T-Complete Basis Functions
We may apply the separation of variables [36].The solution may be in the following form: For simplicity, we let Then, Equation (10) can be rewritten as follows.
We divide U(r)V(θ) on both sides in the above equation and the equation can be rewritten as two differential equations as follows.
Using the constant v to ensure positive or negative constants, we have λ = 0, λ = v 2 and λ = −v 2 .Considering the first scenario λ = 0, we obtain the solutions as follows.
where D 1 , D 2 , D 3 and D 4 are constants.Using the boundary conditions of V(r, 0) = V(r, 2π), we may find that D 1 = 0. Substituting Equations ( 16) and (17) into Equation (11), we have where a 0 and b 0 denote the coefficients.Considering the second scenario, λ = v 2 , we obtain the following solutions.
Water 2019, 11, 835 where D 5 , D 6 , D 7 and D 8 are the coefficients.Inserting the above equations into Equation (11), we obtain ϕ = ar ν cos(νθ) + br ν sin(νθ) + cr −ν cos(νθ) + dr −ν sin(νθ), (21) where a, b, c and d denote the coefficients.Then, we may consider the last scenario λ = −v 2 .Since there is not able to find any non-zero periodic solutions of differential system for U(r), we may only find V(θ) = 0. Collecting all the solutions from the above results, the linearly independent solutions to Laplace equation can be obtained as follows.
In the numerical analysis, we approach the general solution in the form of infinite series of the Laplace equation in a simply connected domain by using a finite number of m.As a result, Equation ( 23) can be rewritten as where m represents the terms of the Trefftz order.The above Equation ( 24) can be used to match the Dirichlet boundary condition.We may also need to consider the Neumann boundary conditions as follows.
Equation ( 25) can be rewritten as where ∇ is the gradient and → n = (n x , n y ) denotes the normal vector.Equation (25) Using Equation ( 24), we may find the derivatives of ∂ϕ/∂r and ∂ϕ/∂θ as follows.

The Characteristic Length
The characteristic length plays a crucial role in controlling the proposed numerical approach in a stable way.Because the matrix assembled with Trefftz trial functions is a full matrix, the resultant system of linear equations may be ill-posed [17,18].The accuracy of the results from the Trefftz method depends sensitively on the order of the T-complete basis functions.Besides, the numerical solution may be unstable.Related to the CTM for solving two-dimensional Laplacian problems, Liu [17] proposed a characteristic length to mitigate the problems of the ill-posedness for the system of linear equations.Applying Dirichlet boundary condition, we obtain Using the CTM, we obtain the approximation solution of the Laplace equation as follows. where x is the coordinate of the collocation points and x ∈ Ω. Applying the Neumann boundary condition, we may obtain the following equations for simply connected domain using the characteristic length.
To mitigate the ill-posedness, the characteristic length [19], R, is adopted and is expressed as where maximum(r) denotes the maximum radial distance in the problem domain.After adopting the characteristic length in our numerical model, the ill-posed phenomenon is greatly reduced, and the accurate numerical solutions can be obtained.Collocating the numerical expansion from Equations ( 32) and ( 34) at boundary collocation points to match the given boundary conditions, we may obtain the following equation.
where T represents a l × M matrix, M = 2m + 1, b represents a M × 1 vector of unknown coefficients, B denotes a vector (size of l × 1) of given functions at boundary points, l represents the number of the boundary points and M represents the terms of the Trefftz order.i ≤ l and j ≤ l in which i and j are the number of boundary points for Dirichlet and Neumann boundary conditions, respectively.g 1 , g 2 , . . ., g i and f 1 , f 2 , . . ., f j denote the boundary values for Dirichlet and Neumann boundary conditions, respectively.In this article, we adopt the domain decomposition method (DDM) [37,38] to solve the nonlinear moving boundary problems in heterogeneous geological media.The DDM is commonly used to solve the problem with different physical characteristics in each subdomain.We first split the domain into two subdomains which are intersected only at the interface.Hence, each subdomain can be regarded as an independent soil layer with its own hydraulic conductivity.At the interface, the flux and the head must satisfy the continuity condition.For instance, we consider a rectangular domain, Ω, which can be split into two intersected subdomains, Ω 1 and Ω 2 .Figure 3 shows that the rectangular domain is divided into Γ 1 , Γ 2 , . . ., Γ 8 sub boundaries where Γ 1 , Γ 2 , . . ., Γ 4 and Γ 5 , Γ 6 , . . ., Γ 8 are sub boundaries of subdomains Ω 1 .and Ω 2 , respectively.At the interface, the sub boundaries, Γ 2 and Γ 6 , are overlapped at the same location.Therefore, additional boundary conditions are imposed on the boundary points to ensure the flux and the head at the interface must be the same. Water

The Iterative Scheme for Solving Free Surface
The nonlinearity of the moving surface flow is caused by the nonlinear boundary conditions.For solving free surface problems with nonlinear boundary conditions, the iterative scheme must be used in the numerical modeling.Due to the difficulty of the computation of the Jacobian matrix for Newton's method, the Picard scheme is adopted in this study.Applying boundary conditions, we obtain where k denotes the index of the collocation points on the free surface to be updated.The head at th J iteration is given as Matching all given boundary conditions, we may obtain a system of linear equations as where T Ω 1 with the size of l 1 × M 1 and T Ω 2 with the size of l 2 × M 2 are the T matrix shown in Equation (37) for Ω 1 and Ω 2 , respectively.The total head can be determined by collocating the inner points within subdomains, Ω 1 and Ω 2 .Consequently, the value of the total head, ϕ, can then be approximated by using Equation (33).

The Iterative Scheme for Solving Free Surface
The nonlinearity of the moving surface flow is caused by the nonlinear boundary conditions.For solving free surface problems with nonlinear boundary conditions, the iterative scheme must be used in the numerical modeling.Due to the difficulty of the computation of the Jacobian matrix for Newton's method, the Picard scheme is adopted in this study.Applying boundary conditions, we obtain where k denotes the index of the collocation points on the free surface to be updated.The head at J th iteration is given as where J denotes the number of iteration steps.We may obtain the following iterative equation [39].
where β is the under-relaxation factor and ϕ J (x k ) is the head to be updated.The value of β is ranging from 0 to 1.The iterative process begins from the first guess values for the moving surface and ceases while the stopping condition expressed in the following equation is achieved.
where ni is the collocation point number on the free surface.

Laminar Flow around a Cylinder
The first example under consideration is a laminar flow around a cylinder.The dimensions of the problem are depicted in Figure 4a.The radius of the cylinder at the center is 1 m.Because the geometry of the problem is clearly symmetric, we considered the half symmetry model.For a two-dimensional domain, Ω, enclosed by a boundary, the Laplace equation is expressed as  The order of the Trefftz basis functions, m, is 150.Totally, 900 collocation points including boundary points and sources are uniformly placed on the domain boundary, as shown in Figure 4a.In order to examine accuracy of the proposed method, 7786 inner points are collocated within the domain boundary.The computed results are validated with SVFLUX [40] which is a finite element seepage analysis program.Figure 4b shows the finite element mesh of SVFLUX.The results of the computed head are compared with the exact solution, as shown in Figure 5.It is found that the numerical solutions agree very well with those of the SVFLUX.

Nonlinear Moving Surface through a Rectangular Dam
The second example is a nonlinear moving surface through a rectangular dam, as shown in Figure 1a.This problem can be viewed as a benchmark problem in geotechnical engineering for finding the position of the moving boundary [3,11,14].The upstream and downstream heads are 24 m and 4 m, respectively.The dimensions of the problem are depicted in Figure 6.This rectangular dam is assembled with five boundary lines, including 1 Γ , 2 Γ , 3 Γ , 4 Γ and 5 Γ , as depicted in Figure 1b.The order of the Trefftz basis functions, m , is 75.We collocate 120 collocation points on moving boundary, as depicted in Figure 6.Since the process for finding the position of the nonlinear free surface is regarded as an inverse problem, the location of the separation point may also need to be examined.In this study, the initial guess of the separation point is at 14.2 m.
The numerical solutions of the free surface are shown in Figure 6 and the results are compared

Nonlinear Moving Surface through a Rectangular Dam
The second example is a nonlinear moving surface through a rectangular dam, as shown in Figure 1a.This problem can be viewed as a benchmark problem in geotechnical engineering for finding the position of the moving boundary [3,11,14].The upstream and downstream heads are 24 m and 4 m, respectively.The dimensions of the problem are depicted in Figure 6.This rectangular dam is assembled with five boundary lines, including Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 , as depicted in Figure 1b.criterion using the proposed iteration scheme.To further explore the accuracy of the computed results, we compare the computed location of the separation point with those from other numerical methods [3,11,14].As depicted in Table 1, the position of the free surface is almost identical with other methods.The results of the separation point calculated by three different methods [3,11,14] are 12.79, 12.68 and 12.84 m, respectively.It is found that the location of the separation point by using the CTM is 12.85 m which is close to that from previous studies.12.79 Chen, Hsiao, Chiu and Lee [11] 12.68 Xiao, Ku, Liu, Fan and Yeih [14] 12.84 The order of the Trefftz basis functions, m, is 75.We collocate 120 collocation points on moving boundary, as depicted in Figure 6.Since the process for finding the position of the nonlinear free surface is regarded as an inverse problem, the location of the separation point may also need to be examined.In this study, the initial guess of the separation point is at 14.2 m.
The numerical solutions of the free surface are shown in Figure 6 and the results are compared with those from other methods, such as Aitchison [3], Chen et al. [11] and Xiao et al. [14].Figure 7 shows that for solving the free surface the number of iteration step is 112 to reach the stopping criterion using the proposed iteration scheme.To further explore the accuracy of the computed results, we compare the computed location of the separation point with those from other numerical methods [3,11,14].As depicted in Table 1, the position of the free surface is almost identical with other methods.The results of the separation point calculated by three different methods [3,11,14] are 12.79, 12.68 and 12.84 m, respectively.It is found that the location of the separation point by using the CTM is 12.85 m which is close to that from previous studies.

Nonlinear Moving Surface through an Earth Dam
The third example is a nonlinear moving surface through an earth dam, as shown in Figure 8.The upstream and downstream hear are 18 m and 8 m, respectively.The dimensions of the problem are depicted in Figure 9a.The boundaries of the earth dam include 1  and 5  are the downstream and upstream boundaries, as depicted in Figure 8.The boundary values are given as 1  is the bottom of the dam where the no-flow condition is given as On 3  , the seepage face boundary condition is as follows.
 , the over-specified moving boundary conditions are as follows.
For this example, the Trefftz basis function order, m , is 150. 750points are collocated on the boundary.The first guess of the free surface is depicted in Figure 9a.The numerical solutions of the free surface are compared with those from other methods [8,14], as shown in Figure 9b.The results illustrate that the computed results are almost identical with other methods.Table 1.Comparison of computed result of the separation point with those from references.

Nonlinear Moving Surface through an Earth Dam
The third example is a nonlinear moving surface through an earth dam, as shown in Figure 8.The upstream and downstream hear are 18 m and 8 m, respectively.The dimensions of the problem are depicted in Figure 9a.The boundaries of the earth dam include Γ 1 , Γ 2 , Γ 3 , Γ 4 and Γ 5 in which Γ 2 and Γ 5 are the downstream and upstream boundaries, as depicted in Figure 8

Flow through Two Layered Soils
The modeling of two-dimensional heterogeneous isotropic subsurface flow in two layered soils is depicted in Figure 3.The coefficients of the permeability with extreme contrasts for two different soils, 1 k and 2 k , are adopted in which the

Flow through Two Layered Soils
The modeling of two-dimensional heterogeneous isotropic subsurface flow in two layered soils is depicted in Figure 3.The coefficients of the permeability with extreme contrasts for two different soils, 1 k and 2 k , are adopted in which the ) ( Γ 1 is the bottom of the dam where the no-flow condition is given as On Γ 3 , the seepage face boundary condition is as follows.
On Γ 4 , the over-specified moving boundary conditions are as follows.
Water 2019, 11, 835 15 of 20 For this example, the Trefftz basis function order, m, is 150.750 points are collocated on the boundary.The first guess of the free surface is depicted in Figure 9a.The numerical solutions of the free surface are compared with those from other methods [8,14], as shown in Figure 9b.The results illustrate that the computed results are almost identical with other methods.

Flow through Two Layered Soils
The modeling of two-dimensional heterogeneous isotropic subsurface flow in two layered soils is depicted in Figure 3.The coefficients of the permeability with extreme contrasts for two different soils, k 1 and k 2 , are adopted in which the k 1 = 10 −1 and k 2 = 10 −15 cm/s.The analytical solution of this example as shown in follows.
where L 1 is the width of the layer 1, L 2 is the width of the layer 2 and ϕ 2 is the head at the interface.The domain is split into two sub-domains which are simply connected.For each sub-domain, 112 boundary points are uniformly collocated.Figure 3 depicts that the rectangular domain boundary is split into eight sub-boundaries: Γ 1 , Γ 2 , . . ., Γ 8 .At the interface, we have the following additional boundary conditions.
Totally, 1520 interior points are collocated within the domain to approximate the head for examining the accuracy.The results of the computed head are compared with the exact solution, as shown in Figure 10.In addition, the accuracy of the results can reach to the order of 10 −11 , as shown in Figure 11.
Water 2019, 11, x FOR PEER REVIEW 16 of 21 where 1 L is the width of the layer 1, 2 L is the width of the layer 2 and 2  is the head at the interface.The domain is split into two sub-domains which are simply connected.For each subdomain, 112 boundary points are uniformly collocated.Figure 3 depicts that the rectangular domain boundary is split into eight sub-boundaries: 1  , 2  ,…, 8  .At the interface, we have the following additional boundary conditions.
Totally, 1520 interior points are collocated within the domain to approximate the head for examining the accuracy.The results of the computed head are compared with the exact solution, as shown in Figure 10.In addition, the accuracy of the results can reach to the order of  x (m)

Nonlinear Moving Surface through a Zoned-Earth Dam
The last case is a nonlinear moving surface through a zoned-earth dam, as shown in Figure 12.For the zoned-earth dam, the upstream and downstream heads are 18 m and 2 m, respectively.The dimensions of the problem are depicted in Figure 12.The values of the permeability for the upstream shell, the core and the downstream shell are For the zoned-earth dam, we divide the domain into three smaller sub-regions, as shown in Figure 13.For each sub-region, we collocate 400, 216 and 191 points on the sub-boundaries for the first, second and third sub-regions, respectively.Besides, we place 50, 66 and 66 collocation points on the moving boundaries, respectively.To validate the results, we compare the computed free surface with that from the finite difference seepage analysis program SEEP2D [41], as shown in Figure 14.It is found that the numerical results agree well with those obtained from the SEEP2D.

Nonlinear Moving Surface through a Zoned-Earth Dam
The last case is a nonlinear moving surface through a zoned-earth dam, as shown in Figure 12.For the zoned-earth dam, the upstream and downstream heads are 18 m and 2 m, respectively.The dimensions of the problem are depicted in Figure 12.The values of the permeability for the upstream shell, the core and the downstream shell are 1.43 × 10 −4 , 1.43 × 10 −5 and 1.43 × 10 −4 cm/s, respectively.

Nonlinear Moving Surface through a Zoned-Earth Dam
The last case is a nonlinear moving surface through a zoned-earth dam, as shown in Figure 12.For the zoned-earth dam, the upstream and downstream heads are 18 m and 2 m, respectively.The dimensions of the problem are depicted in Figure 12.The values of the permeability for the upstream shell, the core and the downstream shell are For the zoned-earth dam, we divide the domain into three smaller sub-regions, as shown in Figure 13.For each sub-region, we collocate 400, 216 and 191 points on the sub-boundaries for the first, second and third sub-regions, respectively.Besides, we place 50, 66 and 66 collocation points on the moving boundaries, respectively.To validate the results, we compare the computed free surface with that from the finite difference seepage analysis program SEEP2D [41], as shown in Figure 14.It is found that the numerical results agree well with those obtained from the SEEP2D.For the zoned-earth dam, we divide the domain into three smaller sub-regions, as shown in Figure 13.For each sub-region, we collocate 400, 216 and 191 points on the sub-boundaries for the first, second and third sub-regions, respectively.Besides, we place 50, 66 and 66 collocation points on the moving boundaries, respectively.To validate the results, we compare the computed free surface with that from the finite difference seepage analysis program SEEP2D [41], as shown in Figure 14.It is found that the numerical results agree well with those obtained from the SEEP2D.

Discussion
In this article, the CTM is adopted to solve the nonlinear moving boundary problems in layered heterogeneous media.Because of the characteristics of the non-linearity, solving nonlinear moving boundary problems with a moving surface remains a challenge.For modeling moving surface problems with nonlinear boundary conditions, the iterative scheme must be used.The sophisticated automatic mesh generation may be required using conventional mesh-based approaches.In addition, the complicated remeshing process in the iterative scheme may lead to problems of the convergence.In this study, we just need to place the boundary points on the domain boundary.Furthermore, only boundary nodes are required to adjust during the iteration process for find the moving boundary.Comparing with the traditional numerical methods, the proposed method is highly accurate.Therefore, the proposed method is advantageous for the nonlinear moving boundary analysis with a phreatic line.
To solve the flow through the layered soils, we adopt the CTM in conjunction with the DDM so that the hydraulic conductivity in each subdomain can be different.To verify the proposed method, numerical examples with nonlinear moving boundary are conducted.Besides, free surface flow through a zoned-earth dam is also carried out.The comparison of the results using the DDM with the exact solution depicts that the CTM with the use of the DDM can reach the accuracy in the order of 11 10 − .Although the CTM may be regarded as an attractive boundary-type meshless method, limitations for the applications may still remain due to the ill-posedness of the method.

Discussion
In this article, the CTM is adopted to solve the nonlinear moving boundary problems in layered heterogeneous media.Because of the characteristics of the non-linearity, solving nonlinear moving boundary problems with a moving surface remains a challenge.For modeling moving surface problems with nonlinear boundary conditions, the iterative scheme must be used.The sophisticated automatic mesh generation may be required using conventional mesh-based approaches.In addition, the complicated remeshing process in the iterative scheme may lead to problems of the convergence.In this study, we just need to place the boundary points on the domain boundary.Furthermore, only boundary nodes are required to adjust during the iteration process for find the moving boundary.Comparing with the traditional numerical methods, the proposed method is highly accurate.Therefore, the proposed method is advantageous for the nonlinear moving boundary analysis with a phreatic line.
To solve the flow through the layered soils, we adopt the CTM in conjunction with the DDM so that the hydraulic conductivity in each subdomain can be different.To verify the proposed method, numerical examples with nonlinear moving boundary are conducted.Besides, free surface flow through a zoned-earth dam is also carried out.The comparison of the results using the DDM with the exact solution depicts that the CTM with the use of the DDM can reach the accuracy in the order of 11 10  .Although the CTM may be regarded as an attractive boundary-type meshless method, limitations for the applications may still remain due to the ill-posedness of the method.

Discussion
In this article, the CTM is adopted to solve the nonlinear moving boundary problems in layered heterogeneous media.Because of the characteristics of the non-linearity, solving nonlinear moving boundary problems with a moving surface remains a challenge.For modeling moving surface problems with nonlinear boundary conditions, the iterative scheme must be used.The sophisticated automatic mesh generation may be required using conventional mesh-based approaches.In addition, the complicated remeshing process in the iterative scheme may lead to problems of the convergence.In this study, we just need to place the boundary points on the domain boundary.Furthermore, only boundary nodes are required to adjust during the iteration process for find the moving boundary.Comparing with the traditional numerical methods, the proposed method is highly accurate.Therefore, the proposed method is advantageous for the nonlinear moving boundary analysis with a phreatic line.
To solve the flow through the layered soils, we adopt the CTM in conjunction with the DDM so that the hydraulic conductivity in each subdomain can be different.To verify the proposed method, numerical examples with nonlinear moving boundary are conducted.Besides, free surface flow through a zoned-earth dam is also carried out.The comparison of the results using the DDM with the exact solution depicts that the CTM with the use of the DDM can reach the accuracy in the order of 10 −11 .Although the CTM may be regarded as an attractive boundary-type meshless method, limitations for the applications may still remain due to the ill-posedness of the method.

Conclusions
This paper presents a study on solving nonlinear moving boundary problems in heterogeneous soils using the CTM.The method is verified by several numerical examples.Application examples are also carried out.The findings are concluded.
The appearance of heterogeneous soils is often found in free surface flow problems.In the past, the CTM is usually limited to linear, homogeneous problems.In this article, we propose a novel idea for solving nonlinear moving boundary problems in layered heterogeneous soils using the collocation meshless method.The results show that the proposed method can be used to deal with nonlinear moving boundary problems in heterogeneity with large permeability contrasts.The method proposed in this study is capable of solving nonlinear moving boundary problems in layered heterogeneous media.However, it is still limited to the layered or zoned porous media in which the porous medium is still homogenous in each zoned domain.In addition, the proposed method can only be applied in saturated soils.

Figure 1 .
Figure 1.Nonlinear moving surface through a rectangular dam.(a) The cross section and boundary conditions and (b) collocation points on the boundary.

Figure 1 .
Figure 1.Nonlinear moving surface through a rectangular dam.(a) The cross section and boundary conditions and (b) collocation points on the boundary.

Figure 2 .
Figure 2. Illustration of the collocation scheme in the CTM.
obtain the solutions as follows.

and 4 D
are constants.Using the boundary conditions of

Figure 2 .
Figure 2. Illustration of the collocation scheme in the CTM.

Figure 3 .
Figure 3.The cross section and the collocation scheme of the two layered soil for the analysis.

2 ΩFigure 3 .
Figure 3.The cross section and the collocation scheme of the two layered soil for the analysis.
l 1 and l 2 are the number of boundary points; M 1 and M 2 are the T-complete basis function order for Ω 1 and Ω 2 , respectively.T I | Γ 2 of the boundary Γ 2 with the size of l I × M 1 and T I | Γ 6 of the boundary Γ 6 with the size of l I × M 2 are matrices at the interface.l I represents the boundary point number at the interface, 0 Ω 1 and 0 Ω 2 are matrices which all values are zero with the size of l 2 × M 1 and l 1 × M 2 , respectively.b Ω 1 denotes a M 1 × 1 vector of unknown coefficients of Ω 1 , b Ω 2 denotes a M 2 × 1 vector of unknown coefficients of Ω 2 .B Ω 1 and B Ω 2 denote vectors of given functions at boundary points of Ω 1 and Ω 2 , respectively.B I = [0 g 0 f ] T , 0 g and 0 f are vectors which all values are zero with the size of l I × 1.

Figure 5 .
Figure 5.Comparison of computed results of the CTM with those from SVFLUX.(a) Computed results of the CTM and (b) SVFLUX results.

Figure 5 .
Figure 5.Comparison of computed results of the CTM with those from SVFLUX.(a) Computed results of the CTM and (b) SVFLUX results.

Figure 6 .
Figure 6.Comparison of moving boundary results with other numerical methods.

Figure 6 .
Figure 6.Comparison of moving boundary results with other numerical methods.

Figure 7 .
Figure 7.The convergence of the iteration step for finding the free surface.

Figure 7 .
Figure 7.The convergence of the iteration step for finding the free surface.

Figure 8 .Figure 9 .
Figure 8.The cross section of the homogeneous earth dam for the analysis.
. The analytical solution of this example as shown in follows.

5 ΓFigure 8 . 21 Figure 8 .Figure 9 .
Figure 8.The cross section of the homogeneous earth dam for the analysis.
. The analytical solution of this example as shown in follows.

5 ΓFigure 9 .
Figure 9. Free surface flow through a homogeneous dam.(a) Initial guess and the computed free surface and (b) comparison of the results with those from other methods.

Figure 10 .
Figure 10.The computed head and the analytical solution in two-layered soils.

Figure 10 .
Figure 10.The computed head and the analytical solution in two-layered soils.

Figure 11 .
Figure 11.Absolute error of the example in two-layered soils.

5  8  6  4  1  7  1  2  3 Figure 12 .
Figure 12.The cross section and soil layer configuration of the zoned-earth dam for the analysis.

Figure 11 .
Figure 11.Absolute error of the example in two-layered soils.

Figure 11 .
Figure 11.Absolute error of the example in two-layered soils.

5 Γ 8 Γ 6 Γ 4 Γ 1 Γ 7 Γ 1 Ω 2 Ω 3 ΩFigure 12 .
Figure 12.The cross section and soil layer configuration of the zoned-earth dam for the analysis.Figure 12.The cross section and soil layer configuration of the zoned-earth dam for the analysis.

Figure 12 .
Figure 12.The cross section and soil layer configuration of the zoned-earth dam for the analysis.Figure 12.The cross section and soil layer configuration of the zoned-earth dam for the analysis.

Figure 13 .Figure 14 .
Figure 13.The collocation scheme of the zoned-earth dam for the DDM.

Figure 14 .
Figure 14.Comparison of the computed free surface with SEEP2D.

Table 1 .
Comparison of computed result of the separation point with those from references.