A High-Order Numerical Manifold Method for Darcy Flow in Heterogeneous Porous Media

One major challenge in modeling Darcy flow in heterogeneous porous media is simulating the material interfaces accurately. To overcome this defect, the refraction law is fully introduced into the numerical manifold method (NMM) as an a posteriori condition. To achieve a better accuracy of the Darcy velocity and continuous nodal velocity, a high-order weight function with a continuous nodal gradient is adopted. NMM is an advanced method with two independent cover systems, which can easily solve both continuous and discontinuous problems in a unified form. Moreover, a regular mathematical mesh, independent of the physical domain, is used in the NMM model. Compared to the conforming mesh of other numerical methods, it is more efficient and flexible. A number of numerical examples were simulated by the new NMM model, comparing the results with the original NMM model and the analytical solutions. Thereby, it is proven that the proposed method is accurate, efficient, and robust for modeling Darcy flow in heterogeneous porous media, while the refraction law is satisfied rigorously.


Introduction
Heterogeneity is the natural property of the geologic structure in groundwater systems.A real groundwater system can be regarded as a combination of several homogeneous subdomains.However, the contrast between the hydraulic conductivities of two adjacent subdomains presents a great challenge for numerical modeling.The two major limitations of modeling Darcy flow in heterogeneous porous media by conventional numerical methods are as follows.(1) The precision of the Darcy velocity (also known as specific discharge) is deficient, as the existing method usually solves the hydraulic head first, then computes the Darcy velocity by taking the derivation of the hydraulic head field.(2) The refraction law, which is expressed as the continuity of normal flux and discontinuity of tangential flux on the interface [1], can be hardly achieved.For the first limitation, the precision of the Darcy velocity decreases to a low order after derivation, and that leads to the discontinuity of the Darcy velocity around a node.For the second limitation, Cainelli et al. [2] summarized the accuracy of classic numerical schemes, such as the finite element method (FEM) [3,4] and finite volume method (FVM) [5][6][7], for modeling flow in heterogeneous formations.Results show that the above methods cannot satisfy the refraction law rigorously when the contrast of conductivities is large, hence the local computation error can be non-negligible.
In order to improve the first limitation (i.e., the lack of precision of the velocity solution), there are two different approaches.The first way is using some post-processing techniques to improve the velocity solution after the hydraulic head is solved [8][9][10].These methods are efficient, as they can be easily modified from the original FEM approach.For example, in Yeh's model [9], the head solution is substituted into the variational formula of Darcy's law, which is then directly solved by the same FEM approach.Although normal flux across the interface is continuous, the tangential component still cannot satisfy refraction law when facing contrasting hydraulic conductivities.The other way is to solve the head and velocity simultaneously, which is known as the mixed finite element method (MFEM) [11,12].However, iterations are still required for MFEM, making the calculation time-consuming [13].
Imposing the refraction law to the governing equation as an additional constraint condition is an efficient approach to dealing with the second limitation.For example, Zhou et al. [14] proposed a modified Yeh's method using a post-processing scheme to enforce the refraction law across material interfaces.Thus, iterations are required for the scheme to achieve convergence.Additionally, Xie et al. [15] developed a domain decomposed finite element method for solving the Darcy velocity in heterogeneous porous media without iterations.The basic idea is to separate the problem domain into several subdomains by material interfaces with coercive refraction law.Nevertheless, the above FEM-based methods are still mesh-dependent, which is inconvenient for pre-processing, especially when material interfaces are complex.
The numerical manifold method (NMM) first presented by Shi [16,17] can be regarded as a generalization to the conventional FEM.Due to the feature of solving both continuous and discontinuous problems in a unity frame, NMM is widely used in mechanical analysis, such as crack propagation [18][19][20][21], rock slopes failure [22][23][24], and dynamic problems [25][26][27][28].NMM also played an important role in fluid flow modeling.Ohnishi et al. [29] presented a node shift method based on NMM for 2D saturated-unsaturated unsteady seepage simulation.Zhang et al. [30] solved Navier-Stokes equations directly by NMM, modeling unsteady incompressible viscous flow.Furthermore, NMM models for unconfined seepage flow in homogeneous media were reported [31,32].However, it cannot exploit the properties of NMM to deal with homogeneous problems, compared with other continuous methods like FEM.Aiming to simulate Darcy flow in heterogeneous media, Hu et al. [33][34][35] developed a continuous approach with jump functions and a discontinuous approach with Lagrange multipliers, by enforcing the continuity of normal flux across material interfaces, while the refraction law for tangential flux is ignored.Zheng et al. [36] also proposed an NMM model for unconfined seepage flow by adopting the moving least squares (MLS) interpolation.Hence, NMM is an efficient approach for heterogeneous seepage problems.However, improving the accuracy of modeling, that is still a great challenge.
In this work, we developed a high-order NMM model for Darcy flow in heterogeneous porous media with a nodal-continuous weight function, as well as a local cover function derived by first order approximation of Taylor expansion.Then we introduced refraction law into the NMM model as an a posteriori constraint condition entirely to ensure that material interfaces are simulated correctly.The proposed NMM model has the following advantages.(1) The accuracy of the Darcy velocity solution is high and keeps continuous around a node, both of which are significant for modeling underground flow as well as solute transport.(2) The treatment of material interfaces in the NMM model is accurate and efficient, especially when facing various materials with interfaces intersecting each other, by taking full advantage of NMM's dual cover systems.More details of the proposed NMM model will be discussed in the following sections.
The framework of this paper is organized as follows.Section 2 briefly introduces some basic concepts of the original NMM theory.Then, the new NMM formation for Darcy flow in heterogeneous porous media is derived in Section 3.With the aim of testing the reliability of the proposed NMM model, numerical examples are simulated in Section 4. Finally, main conclusions and discussions are given in Section 5.

Basic Theory of NMM
Basic concepts of NMM are introduced in this section, including the dual cover systems and the manifold elements.The local approximation on physical covers and the global approximation on manifold elements are proposed next.This section shows that the NMM is suitable for modeling discontinuous problems because of the efficient meshing algorithm and the treatment of discontinuous boundaries.

Dual Cover Systems and Manifold Elements
The basic concepts of NMM include mathematical cover (MC), physical cover (PC), and manifold element (ME).The mathematical cover system is a set of patches which completely cover the whole physical domain.Thus, each patch is an MC.PCs are generated from MCs, cutting via boundaries or cracks and material interfaces inside the physical area.More remarkably, the shape of MC may be arbitrary, which is not necessary in coinciding with the boundaries as long as the whole domain is covered.Here, standard FEM mesh is used to generate MCs for convenience, and the density of mesh decides computational accuracy.
A simple example is given to illustrate the above concepts.As shown in Figure 1, a rectangular domain is divided into two different materials by a vertical interface.A regular triangular mesh is used to form MCs in this example.Each node in the mesh is defined as a star, and MCs are formed by triangles adjacent to a star, denoted as P i , where the subscribe i is the number of the star.Thus, the standard shape of MC in this example is a hexagon (i.e., P 5 and P 6 ) combined by six triangles.However, quadrangular or triangular MCs may be made using fewer triangles on the boundaries (i.e., P 9 ).
Processes 2018, 6, x FOR PEER REVIEW 3 of 20 discontinuous problems because of the efficient meshing algorithm and the treatment of discontinuous boundaries.

Dual Cover Systems and Manifold Elements
The basic concepts of NMM include mathematical cover (MC), physical cover (PC), and manifold element (ME).The mathematical cover system is a set of patches which completely cover the whole physical domain.Thus, each patch is an MC.PCs are generated from MCs, cutting via boundaries or cracks and material interfaces inside the physical area.More remarkably, the shape of MC may be arbitrary, which is not necessary in coinciding with the boundaries as long as the whole domain is covered.Here, standard FEM mesh is used to generate MCs for convenience, and the density of mesh decides computational accuracy.
A simple example is given to illustrate the above concepts.As shown in Figure 1, a rectangular domain is divided into two different materials by a vertical interface.A regular triangular mesh is used to form MCs in this example.Each node in the mesh is defined as a star, and MCs are formed by triangles adjacent to a star, denoted as , where the subscribe is the number of the star.Thus, the standard shape of MC in this example is a hexagon (i.e., and ) combined by six triangles.However, quadrangular or triangular MCs may be made using fewer triangles on the boundaries (i.e., ).MCs are divided into independent parts by outside boundaries or material interfaces.Each part of MCs is regarded as a PC.This is denoted as , where the superscript is the number of the PC belonging to the i-th MC.As illustrated in Figure 2, is divided into and by boundaries and material interface.Similarly, is divided into and , and is divided into and .

Figure 2.
The dual cover systems: mathematical covers and physical covers.
After generating the cover systems, manifold elements-the basic integral region of the weak form of the problems-are defined as the overlapping areas of PCs.As exhibited in Figure 3, MEs are regarded as the common area of three layers of PCs when adopting triangular finite element mesh.Therefore, the star of the corresponding PCs is considered as a generalized node of the ME in NMM.MCs are divided into independent parts by outside boundaries or material interfaces.Each part of MCs is regarded as a PC.This is denoted as P j i , where the superscript j is the number of the PC belonging to the i-th MC.As illustrated in Figure 2, P 5 is divided into P 1 5 and P 2 5 by boundaries and material interface.Similarly, P 6 is divided into P 1 6 and P 2 6 , and P 9 is divided into P 1 9 and P 2 9 .
Processes 2018, 6, x FOR PEER REVIEW 3 of 20 discontinuous problems because of the efficient meshing algorithm and the treatment of discontinuous boundaries.

Dual Cover Systems and Manifold Elements
The basic concepts of NMM include mathematical cover (MC), physical cover (PC), and manifold element (ME).The mathematical cover system is a set of patches which completely cover the whole physical domain.Thus, each patch is an MC.PCs are generated from MCs, cutting via boundaries or cracks and material interfaces inside the physical area.More remarkably, the shape of MC may be arbitrary, which is not necessary in coinciding with the boundaries as long as the whole domain is covered.Here, standard FEM mesh is used to generate MCs for convenience, and the density of mesh decides computational accuracy.
A simple example is given to illustrate the above concepts.As shown in Figure 1, a rectangular domain is divided into two different materials by a vertical interface.A regular triangular mesh is used to form MCs in this example.Each node in the mesh is defined as a star, and MCs are formed by triangles adjacent to a star, denoted as , where the subscribe is the number of the star.Thus, the standard shape of MC in this example is a hexagon (i.e., and ) combined by six triangles.However, quadrangular or triangular MCs may be made using fewer triangles on the boundaries (i.e., ).MCs are divided into independent parts by outside boundaries or material interfaces.Each part of MCs is regarded as a PC.This is denoted as , where the superscript is the number of the PC belonging to the i-th MC.As illustrated in Figure 2, is divided into and by boundaries and material interface.Similarly, is divided into and , and is divided into and .

Figure 2.
The dual cover systems: mathematical covers and physical covers.
After generating the cover systems, manifold elements-the basic integral region of the weak form of the problems-are defined as the overlapping areas of PCs.As exhibited in Figure 3, MEs are regarded as the common area of three layers of PCs when adopting triangular finite element mesh.Therefore, the star of the corresponding PCs is considered as a generalized node of the ME in NMM.After generating the cover systems, manifold elements-the basic integral region of the weak form of the problems-are defined as the overlapping areas of PCs.As exhibited in Figure 3, MEs are regarded as the common area of three layers of PCs when adopting triangular finite element mesh.Therefore, the star of the corresponding PCs is considered as a generalized node of the ME in NMM.This is different to the concept of nodes in FEM.For simplicity, the generalized node will be termed as "node" for the rest of this paper.In addition, the shape of MEs are generally triangular, except for the MEs adjacent to boundaries or interfaces, which may have an arbitrary shape.As long as the MCs are based on triangular mesh, the number of generalized nodes of an ME is three.As shown in Figure 3, ME is overlapped by P 1 5 , P 1 6 , and P 1 9 .Thus, ME is termed as P This is different to the concept of nodes in FEM.For simplicity, the generalized node will be termed as "node" for the rest of this paper.In addition, the shape of MEs are generally triangular, except for the MEs adjacent to boundaries or interfaces, which may have an arbitrary shape.As long as the MCs are based on triangular mesh, the number of generalized nodes of an ME is three.As shown in Figure 3, ME○ 1 is overlapped by , , and .Thus, ME○ 1 is termed as , and ME○ 2 is similarly termed as .In this way, the discontinuity of the triangular element cut by boundaries or material interface is presented by NMM.Theoretically, NMM is a kind of partition of unity method (PUM).Its important feature is that the global approximation is established by combining local approximations using a PU function.According to the theory of NMM, local approximation is defined on each PC as a cover function, while global approximation on each ME is built through combining cover functions using weight functions, which are the same as PU functions.Details of NMM approximation is introduced in the following two subsections.

Local Approximation on Physical Covers
Local approximation is defined on each PC, representing the local field variables.That is, normally expressed as where is the total number of PCs, x are coordinates of points within the i-th PC, is a vector of unknown coefficients, and is the vector of polynomials bases.The latter can be expressed as , , , , , ⋯ , , ⋯ , for the first, second, and -th order approximation, respectively, where is a positive integer representing the approximation order.
However, in contrast with polynomial bases, specific bases can be selected as on each PC for achieving batter approximation.As a result, NMM has the ability to combine analytical solution with numerical method, that is to say, corresponding analytical solution can be chosen as the local approximation in the key position of problem domain, while polynomials bases are adopted on the rest of domain.Thus, a more accurate and efficient solution can be achieved.For example, Williams' displacement series are usually used as the local approximation near a crack tip for crack problems, making the stress solution around tip more accurate and reliable [18,19].Theoretically, NMM is a kind of partition of unity method (PUM).Its important feature is that the global approximation is established by combining local approximations using a PU function.According to the theory of NMM, local approximation is defined on each PC as a cover function, while global approximation on each ME is built through combining cover functions using weight functions, which are the same as PU functions.Details of NMM approximation is introduced in the following two subsections.

Local Approximation on Physical Covers
Local approximation is defined on each PC, representing the local field variables.That is, normally expressed as where n PC is the total number of PCs, x are coordinates of points within the i-th PC, a is a vector of unknown coefficients, and p T (x) is the vector of polynomials bases.The latter can be expressed as for the first, second, and n-th order approximation, respectively, where n is a positive integer representing the approximation order.
However, in contrast with polynomial bases, specific bases can be selected as p T (x) on each PC for achieving batter approximation.As a result, NMM has the ability to combine analytical solution with numerical method, that is to say, corresponding analytical solution can be chosen as the local approximation in the key position of problem domain, while polynomials bases are adopted on the rest of domain.Thus, a more accurate and efficient solution can be achieved.For example, Williams' displacement series are usually used as the local approximation near a crack tip for crack problems, making the stress solution around tip more accurate and reliable [18,19].

Global Approximation on Manifold Elements
Based on the NMM theory, global approximation is defined on each ME.The field variables ϕ on MEs are obtained by taking a weighted average of local approximations on corresponding PCs.Therefore, global approximation can be expressed as where N PC is the number of PCs on the ME and w i (x, y) is weight function of the i-th PC, which satisfies where U i is the geometric range of the i-th PC.

Steady Darcy Flow in Heterogeneous Porous Media by NMM
In this section, NMM formulation of steady Darcy flow in heterogeneous porous media is derived, including treatments for boundaries and material interfaces.

Governing Equations
Considering 2D steady Darcy flow in heterogeneous porous media, the continuity equation is derived using the conservation of mass, expressed as where ∇ is the gradient operator, v is the Darcy velocity or specific discharge, and Q is the source term.The constitutive equation for Darcy flow is Darcy's law, which means there is a linear relationship between the Darcy velocity and the hydraulic gradient.This is given as where K is the permeability tensor and h is the hydraulic head.Substituting Equation (6) into Equation ( 5), the governing differential equation can be rewritten as As shown in Figure 4, the associated boundary conditions can be written as for the Dirichlet condition on Γ h , q n = q (9) for the Neumann condition on Γ q , and for the refraction law across the material interface Γ m , where h is the given hydraulic head on Γ h , q n is normal flux on Γ q , q is the given flux, and v n and v s are normal and tangential components of the Darcy velocity, respectively.Superscript + and -represent the two sides of the material interface Γ m .In addition, v n and v s can be expressed as where n and s are the normal unit vector and tangential unit vector, respectively.where and are the normal unit vector and tangential unit vector, respectively.The first two conditions of refraction law can be easy deduced.They represent the continuity of the hydraulic head and normal flux.Here we briefly prove the third condition for tangential flux.As shown in Figure 5, a rectangular domain is divided into two different materials Ω and Ω , both of which are homogeneous and isotropic media.The permeability coefficients are and .Suppose fluid flows upside to downside, and and are the angles between the direction of the Darcy velocity and the normal direction of the material interface in Ω and Ω , respectively.From Equation (10a) it follows that (12) across the material interface, where is hydraulic gradient.Thus, the tangential component of the hydraulic gradient is continuous across the material interface.According to Darcy's law, the tangential component of the Darcy velocity can be expressed as Figure 5. Refraction law across the material interface.
Substituting Equation ( 13) into Equation ( 12) and considering the geometric condition, the tangential velocity in two sides of the interface satisfies The first two conditions of refraction law can be easy deduced.They represent the continuity of the hydraulic head and normal flux.Here we briefly prove the third condition for tangential flux.As shown in Figure 5, a rectangular domain is divided into two different materials Ω 1 and Ω 2 , both of which are homogeneous and isotropic media.The permeability coefficients are k + and k − .Suppose fluid flows upside to downside, and θ 1 and θ 2 are the angles between the direction of the Darcy velocity and the normal direction of the material interface in Ω 1 and Ω 2 , respectively.From Equation (10) across the material interface, where j is hydraulic gradient.Thus, the tangential component of the hydraulic gradient is continuous across the material interface.According to Darcy's law, the tangential component of the Darcy velocity can be expressed as Substituting Equation ( 13) into Equation ( 12) and considering the geometric condition, the tangential velocity in two sides of the interface satisfies which is equivalent to the third condition in Equation (10).
Through this subsection, NMM formation for Darcy flow in heterogeneous media is derived.The main purpose of this study is to solve the differential equation of Darcy flow, Equation (7) with constraint conditions of Equations ( 8)- (10) for the hydraulic head and the Darcy velocity using NMM.
velocity and the normal direction of the material interface in Ω and Ω , respectively.From Equation (10a) it follows that (12) across the material interface, where is hydraulic gradient.Thus, the tangential component of the hydraulic gradient is continuous across the material interface.According to Darcy's law, the tangential component of the Darcy velocity can be expressed as Figure 5. Refraction law across the material interface.
Substituting Equation ( 13) into Equation ( 12) and considering the geometric condition, the tangential velocity in two sides of the interface satisfies tan tan (14) which is equivalent to the third condition in Equation (10).Through this subsection, NMM formation for Darcy flow in heterogeneous media is derived.The main purpose of this study is to solve the differential equation of Darcy flow, Equation ( 7) with constraint conditions of Equations ( 8)- (10) for the hydraulic head and the Darcy velocity using NMM.

NMM Interpolations
In this paper, Taylor series expansion was applied to establish the local approximation on each PC.Thus, the base function on a PC can be adopted by first order approximation of the Taylor

NMM Interpolations
In this paper, Taylor series expansion was applied to establish the local approximation on each PC.Thus, the base function on a PC can be adopted by first order approximation of the Taylor expansion of the corresponding node.According to Equation ( 1), the local approximation of the hydraulic head on each PC, h i (x), can be written as where a is a vector of unknown coefficients, and bases are expressed as where x i and y i are coordinates of node i. Substituting Equation ( 16) into Equation ( 15), the local approximation on PCs can be rewritten as Comparing with conventional polynomials bases, the exact physical meanings of the unknown coefficients, the corresponding hydraulic head, and the gradient on each node are achieved.
Then, global approximation on each ME can be derived as where N is the shape function, expressed as and w i (x) is the weight function on the i-th PC, originally used in mechanics problems [37,38].This is expressed as where L i (i = 1, 2, 3) are area coordinates of the triangular finite element.The above weight function has four important features: (1) the Kronecker delta property (w i x j = δ ij (i, j = 1, 2, 3)), (2) the non-negative property (0 the partition of the unity condition (∑ 3 i=1 w i (x) = 1), ( 4) and a nodal gradient of zero of the weight function.
Therefore, the gradient of the hydraulic head solution at a generalized node is continuous, because of the last property, which is of great importance in this study.A continuous Darcy velocity can be then derived by Darcy's law at nodes.However, the discontinuity of velocity still exists on element edges.

Discrete Equations
The weak form of governing Equation ( 7) is derived by the weighted residual method, and then the functional expression derived by the variational principle can be written as [39] In each ME, the hydraulic head, h, is defined as where N is the shape function, and a is the unknowns.The hydraulic gradient, j, is defined as where B = ∇N.
The Darcy velocity, v, is expressed as where K is the permeability tensor.By substituting Equations ( 22)- (24) into Equation ( 21), the functional expression of each ME can be obtained.
where Ω e is the domain of each ME.Therefore, the functional expression of the whole problem domain, Equation (21), can be rewritten in a matrix form where n is the number of unknowns in the problem, C is the conductivity matrix with n × n components, and F is the flux term with n components, expressed as respectively.Their components can be derived by taking the derivatives with the unknowns, expressed as According to the variational principle, the solution of the weak form of the differential equation is equivalent to a function a which makes Π stationary with respect to arbitrary changes of δa, that is where δa 1 , δa 2 , . . ., δa n are arbitrary changes.Therefore, Equation ( 7) reduces to a standard linear form as which can be rewritten as a matrix form Ca = F (33) and Therefore, the governing differential equation, Equation ( 7) is transformed a system of linear equations (Equation ( 34)), which can be easily solved.In addition, the hydraulic head and Darcy velocity can be obtained using Equations ( 22) and (24), respectively.

Treatments for Boundaries and Material Interfaces
In Section 3.3, the Dirichlet boundary condition and refraction law across material interfaces are not involved.The treatments for boundaries and material interfaces are introduced in this subsection.Generally, two basic methods to impose additional constraints are the Lagrange multiplier method and the penalty method.Although the Lagrange multiplier method can achieve an accurate solution, additional degrees of freedom must be introduced into the governing function, which usually leads to computational difficulties.The penalty method is more convenient as it does not require a change of the size of the coefficient matrix and is free of the above drawbacks of the Lagrange multiplier method.Note that the penalty number should be neither too small nor too large to enforce the constraints and avoid computational difficulties.In this study, the penalty method is adopted, and penalty number α is determined as 10 6 .The accuracy and robustness are verified by a series of numerical tests.
In addition, the treatments of material interfaces by existing methods are inadequate.In this paper, additional constraints of the Dirichlet boundary condition and refraction law across the material interface are introduced by forming a new functional in which Π D is the Dirichlet boundary condition term and Π m is material interface term.
For the Dirichlet boundary condition, Π D is expressed as where α is the penalty number.Thus, the components of the conductivity matrix, C ij , and the flux term, F i , for the equilibrium equation can be derived as For the material interface, Π m can be written into three parts Each part of Π m can be expressed as Thus, the components of the conductivity matrix can be respectively derived as Because the constraints of refraction law are fully introduced into the NMM model as a posteriori conditions, the material interface can be simulated accurately.The following section will present some numerical examples to verify the new NMM model for Darcy flow in heterogeneous porous media.

Numerical Examples
To verify the proposed NMM model, a set of numerical examples are carried out.First, two rectangular homogeneous models with different boundaries are studied.They are compared with analytical solutions in order to demonstrate the accuracy in modeling homogeneous Darcy flow.Then, a heterogeneous model, with hydraulic conductivity varying smoothly within the computational domain, is simulated to verify the capability and accuracy of the model for continuous-heterogeneous problems.Lastly, the new model is applied to a heterogeneous problem with material interfaces, showing the accuracy and efficiency for Darcy flow in heterogeneous porous media.

Example 1: Verification of Homogeneous Darcy Flow with an Analytical Solution
In order to verify the accuracy for a Darcy flow simulation before modeling heterogeneous problems, a homogeneous model is studied.As shown in Figure 6a, a 1 × 1 m rectangular domain with Dirichlet boundary conditions is as follows: h(0, y) = 0, h(x, 0) = 0, h(1, y) = y, and h(x, 1) = x.The model is homogeneous and isotropic with a permeability coefficient of k(x, y) = 1 m/s.The analytical solution of the hydraulic head is h(x, y) = xy.Thus, the Darcy velocity can be derived by Darcy's law or Equation ( 6) as follows:

Example 1: Verification of Homogeneous Darcy Flow with an Analytical Solution
In order to verify the accuracy for a Darcy flow simulation before modeling heterogeneous problems, a homogeneous model is studied.As shown in Figure 6a, a 1 1 m rectangular domain with Dirichlet boundary conditions is as follows: 0, 0, , 0 . The model is homogeneous and isotropic with a permeability coefficient of , 1 m/s.The analytical solution of the hydraulic head is , .Thus, the Darcy velocity can be derived by Darcy's law or Equation ( 6) as follows: , .
(a) (b) The numerical results solved by the proposed method with a 10-layers mesh are shown in Figure 7. Figure 7a shows the distribution of the hydraulic head in the problem domain, and Figure 7b,c are the distributions of the Darcy velocity component and , respectively.It is obvious that the hydraulic head and Darcy velocity are both continuous in the problem domain.To realize the advanced feature of the new model, a contrast simulation using the conventional NMM model is studied.Figure 8a is the result of the hydraulic head, which is similar to Figure 7a.Compared with the continuous Darcy velocity in Figure 7b,c, the velocity field obtained though the conventional model is discontinuous both around nodes and across element edges (Figure 8b,c).Therefore, it is proven that the proposed NMM model can achieve a continuous nodal velocity.
The numerical results solved by the proposed method with a 10-layers mesh are shown in Figure 7. Figure 7a shows the distribution of the hydraulic head in the problem domain, and Figure 7b,c are the distributions of the Darcy velocity component v x and v y , respectively.It is obvious that the hydraulic head and Darcy velocity are both continuous in the problem domain.To realize the advanced feature of the new model, a contrast simulation using the conventional NMM model is studied.Figure 8a is the result of the hydraulic head, which is similar to Figure 7a.Compared with the continuous Darcy velocity in Figure 7b,c, the velocity field obtained though the conventional model is discontinuous both around nodes and across element edges (Figure 8b,c).Therefore, it is proven that the proposed NMM model can achieve a continuous nodal velocity.Since the analytical solution is known in this problem, we can compare the numerical result with its exact result.Figure 9a shows the comparison of the hydraulic head solution along a profile located at x = 0.3 m, and Figure 9b,c   Since the analytical solution is known in this problem, we can compare the numerical result with its exact result.Figure 9a shows the comparison of the hydraulic head solution along a profile located at x = 0.3 m, and Figure 9b,c  Furthermore, in order to analyse the accuracy and convergence of the numerical solution, the relative error of the hydraulic head in the L 2 norm can be written as  Since the analytical solution is known in this problem, we can compare the numerical result with its exact result.Figure 9a shows the comparison of the hydraulic head solution along a profile located at x = 0.3 m, and Figure 9b,c  Furthermore, in order to analyse the accuracy and convergence of the numerical solution, the relative error of the hydraulic head in the L 2 norm can be written as  Furthermore, in order to analyse the accuracy and convergence of the numerical solution, the relative error of the hydraulic head in the L 2 norm can be written as where h is the numerical solution of the hydraulic head and ĥ is the analytical solution.Similarly, the relative error of the Darcy velocity in the L 2 norm is defined as where v is the numerical solution of the Darcy velocity and v is the analytical solution.
After convergence analysis for both the new NMM model and the conventional NMM model, the convergence curves of the hydraulic head and Darcy velocity are plotted in Figure 10a,b, respectively.As shown in Figure 10a, the two curves are almost coincident, indicating that they have the same accuracy and convergence rate for the solution of the hydraulic head.However, Figure 10b shows that the velocity solution of the new NMM model is more accurate than the convention NMM model, and their convergence rates are still the same.The values of relative error and convergence rate of both hydraulic head and Darcy velocity by the two NMM models are shown in Tables 1 and 2, respectively.
After convergence analysis for both the new NMM model and the conventional NMM model, the convergence curves of the hydraulic head and Darcy velocity are plotted in Figure 10a,b, respectively.As shown in Figure 10a, the two curves are almost coincident, indicating that they have the same accuracy and convergence rate for the solution of the hydraulic head.However, Figure 10b shows that the velocity solution of the new NMM model is more accurate than the convention NMM model, and their convergence rates are still the same.The values of relative error and convergence rate of both hydraulic head and Darcy velocity by the two NMM models are shown in Tables 1 and  2, respectively.To demonstrate the capability of handling complex boundaries by the proposed method, the same model with a quadratic Dirichlet boundary condition is simulated, as shown in Figure 6b.The boundary conditions are h(0, y) = 0, h(x, 0) = 0, h(1, y) = y 2 , and h(x, 1) = x, and the permeability coefficient remains unchanged.Note that the source term in this model is Q = 2x.The analytical solutions of the hydraulic head and Darcy velocity are h(x, y) = xy 2 , v x = −y 2 , and v y = −2xy, respectively.
First, the problem is simulated by the proposed NMM model with a 60-layers mesh.As shown in Figure 11, the distributions of the hydraulic head and Darcy velocity are presented.From Figure 11a, an obvious difference in hydraulic head distribution caused by the change of boundary conditions can be observed compared to Figure 8a.Furthermore, the velocity field is also continuous in the problem domain (Figure 11b,c).
boundary conditions are 0, 0 , , 0 , and the permeability coefficient remains unchanged.Note that the source term in this model is 2 .The analytical solutions of the hydraulic head and Darcy velocity are , , , and 2 , respectively.
First, the problem is simulated by the proposed NMM model with a 60-layers mesh.As shown in Figure 11, the distributions of the hydraulic head and Darcy velocity are presented.From Figure 11a, an obvious difference in hydraulic head distribution caused by the change of boundary conditions can be observed compared to Figure 8a.Furthermore, the velocity field is also continuous in the problem domain (Figure 11b,c).By using different mesh densities, the accuracy of the proposed method for complex boundary conditions is verified.As shown in Figure 12a, the hydraulic head solution along a profile located at x = 0.3 m, in spite of different mesh densities, almost achieves its analytical solution.The Darcy velocity solution is also reliable compared with its analytical result (Figure 12b,c).In addition, when the mesh density reaches 60 layers, the Darcy velocity solution almost achieves its analytical solution.Convergence analysis is performed both for hydraulic heads and velocities.Figure 13a shows the convergence curve of the hydraulic head, and Figure 13b shows the convergence curve of the  By using different mesh densities, the accuracy of the proposed method for complex boundary conditions is verified.As shown in Figure 12a, the hydraulic head solution along a profile located at x = 0.3 m, in spite of different mesh densities, almost achieves its analytical solution.The Darcy velocity solution is also reliable compared with its analytical result (Figure 12b,c).In addition, when the mesh density reaches 60 layers, the Darcy velocity solution almost achieves its analytical solution.
same model with a quadratic Dirichlet boundary condition is simulated, as shown in Figure 6b.The boundary conditions are 0, 0 , , 0 0 , 1, , and , 1 , and the permeability coefficient remains unchanged.Note that the source term in this model is 2 .The analytical solutions of the hydraulic head and Darcy velocity are , , , and 2 , respectively.
First, the problem is simulated by the proposed NMM model with a 60-layers mesh.As shown in Figure 11, the distributions of the hydraulic head and Darcy velocity are presented.From Figure 11a, an obvious difference in hydraulic head distribution caused by the change of boundary conditions can be observed compared to Figure 8a.Furthermore, the velocity field is also continuous in the problem domain (Figure 11b,c).By using different mesh densities, the accuracy of the proposed method for complex boundary conditions is verified.As shown in Figure 12a, the hydraulic head solution along a profile located at x = 0.3 m, in spite of different mesh densities, almost achieves its analytical solution.The Darcy velocity solution is also reliable compared with its analytical result (Figure 12b,c).In addition, when the mesh density reaches 60 layers, the Darcy velocity solution almost achieves its analytical solution.Convergence analysis is performed both for hydraulic heads and velocities.Figure 13a shows the convergence curve of the hydraulic head, and Figure 13b shows the convergence curve of the  Convergence analysis is performed both for hydraulic heads and velocities.Figure 13a shows the convergence curve of the hydraulic head, and Figure 13b shows the convergence curve of the Darcy velocity.However, there is no big difference with Figure 10, except that the relative errors of both the hydraulic head and Darcy velocity are larger than last examples, since the problem has more complex boundary conditions.
Darcy velocity.However, there is no big difference with Figure 10, except that the relative errors of both the hydraulic head and Darcy velocity are larger than last examples, since the problem has more complex boundary conditions.

Example 2: Accuracy Darcy Velocity Solution for Continuous-Heterogeneous Media
A continuous-heterogeneous model, as shown in Figure 14, is simulated in this example.The boundary conditions are as follows: , 5 , 5 25 and 5, 5, 25 .The hydraulic conductivity varies smoothly within the problem domain, expressed as k , 1e . In addition, a source term 4e , is included in this model.The analytical solution of hydraulic head is , , and the solution of the velocity field can be derived as 2e and 2e .Figure 15a shows the distribution of hydraulic head, whereas Figure 15b,c show the distribution of the Darcy velocity and , respectively.The contour plots show that the distribution of numerical results is symmetric and continuous by the proposed model.In order to study the robustness of the cover system, numerical results along a profile at y = 0 m with 20, 40, and 60 layers of mesh are compared with the analytical solution.Figure 16a shows the comparison of the hydraulic head solution between the proposed NMM model and the analytical solution along a profile located at y = 0 m, indicating high precision for the hydraulic head solution.Furthermore, the Darcy velocity solution also shows good agreement with the analytical solution in this continuous-heterogeneous media in Figure 16b,c.

Example 2: Accuracy Darcy Velocity Solution for Continuous-Heterogeneous Media
A continuous-heterogeneous model, as shown in Figure 14, is simulated in this example.The boundary conditions are as follows: h(x, −5) = h(x, 5) = x 2 − 25 and h(−5, y) = h(5, y) = 25 − y 2 .The hydraulic conductivity varies smoothly within the problem domain, expressed as k(x, y) = 1e −3 × x 2 + y 2 .In addition, a source term Q = 4e −3 × x 2 − y 2 , is included in this model.The analytical solution of hydraulic head is h(x, y) = x 2 − y 2 , and the solution of the velocity field can be derived as v x = −2e −3 × x 3 + xy 2 and v y = 2e −3 × x 2 y + y 3 .
Processes 2018, 6, x FOR PEER REVIEW 14 of 20 Darcy velocity.However, there is no big difference with Figure 10, except that the relative errors of both the hydraulic head and Darcy velocity are larger than last examples, since the problem has more complex boundary conditions.

Example 2: Accuracy Darcy Velocity Solution for Continuous-Heterogeneous Media
A continuous-heterogeneous model, as shown in Figure 14, is simulated in this example.The boundary conditions are as follows: , The hydraulic conductivity varies smoothly within the problem domain, expressed as k , 1e . In addition, a source term 4e , is included in this model.The analytical solution of hydraulic head is , , and the solution of the velocity field can be derived as 2e and 2e .Figure 15a shows the distribution of hydraulic head, whereas Figure 15b,c show the distribution of the Darcy velocity and , respectively.The contour plots show that the distribution of numerical results is symmetric and continuous by the proposed model.In order to study the robustness of the cover system, numerical results along a profile at y = 0 m with 20, 40, and 60 layers of mesh are compared with the analytical solution.Figure 16a shows the comparison of the hydraulic head solution between the proposed NMM model and the analytical solution along a profile located at y = 0 m, indicating high precision for the hydraulic head solution.Furthermore, the Darcy velocity solution also shows good agreement with the analytical solution in this continuous-heterogeneous media in Figure 16b,c.
(-5,- Figure 15a shows the distribution of hydraulic head, whereas Figure 15b,c show the distribution of the Darcy velocity v x and v y , respectively.The contour plots show that the distribution of numerical results is symmetric and continuous by the proposed model.In order to study the robustness of the cover system, numerical results along a profile at y = 0 m with 20, 40, and 60 layers of mesh are compared with the analytical solution.Figure 16a shows the comparison of the hydraulic head solution between the proposed NMM model and the analytical solution along a profile located at y = 0 m, indicating high precision for the hydraulic head solution.Furthermore, the Darcy velocity solution also shows good agreement with the analytical solution in this continuous-heterogeneous media in Figure 16b,c.( For this continuous-heterogeneous problem, convergence curves of both the hydraulic head and velocity are shown in Figure 17.According to Figure 17a, the hydraulic head solution found using the proposed NMM model is more accurate than that of the conventional NMM model, and the convergence rate is still the same.Moreover, the proposed high-order NMM model also shows better accuracy of the Darcy velocity solution than the conventional model, while their convergence rates are almost the same.

Example 3: Verification of Refraction Law for Discontinuous-Heterogeneous Media
The last example is to verify the refraction law across material interfaces [14], as shown in Figure 18a ( For this continuous-heterogeneous problem, convergence curves of both the hydraulic head and velocity are shown in Figure 17.According to Figure 17a, the hydraulic head solution found using the proposed NMM model is more accurate than that of the conventional NMM model, and the convergence rate is still the same.Moreover, the proposed high-order NMM model also shows better accuracy of the Darcy velocity solution than the conventional model, while their convergence rates are almost the same.

Example 3: Verification of Refraction Law for Discontinuous-Heterogeneous Media
The last example is to verify the refraction law across material interfaces [14], as shown in Figure 18a For this continuous-heterogeneous problem, convergence curves of both the hydraulic head and velocity are shown in Figure 17.According to Figure 17a, the hydraulic head solution found using the proposed NMM model is more accurate than that of the conventional NMM model, and the convergence rate is still the same.Moreover, the proposed high-order NMM model also shows better accuracy of the Darcy velocity solution than the conventional model, while their convergence rates are almost the same.For this continuous-heterogeneous problem, convergence curves of both the hydraulic head and velocity are shown in Figure 17.According to Figure 17a, the hydraulic head solution found using the proposed NMM model is more accurate than that of the conventional NMM model, and the convergence rate is still the same.Moreover, the proposed high-order NMM model also shows better accuracy of the Darcy velocity solution than the conventional model, while their convergence rates are almost the same.The last example is to verify the refraction law across material interfaces [14], as shown in Figure 18a.Considering a rectangular domain of 120 × 120 m with material interfaces AB and CD, the permeability coefficients are k 1 = 100 m/d and k 2 = 10 m/d, respectively.The hydraulic heads are fixed with 1 m on the left side and 0 m on the right side, while the top boundary has a specified inflow flux of q = 1 m/d.The bottom of the model is impermeable.Due to the dual cover systems, the NMM mesh (Figure 18b) does not need to conform with the material interfaces, unlike the conventional FEM mesh (Figure 18c).In order to study the influence of the material interface, numerical results with two different mesh densities along profiles at x = 100 m and y = 80 m are plotted in Figures 20 and 21.As shown in Figure 20, the hydraulic head across the material interface is continuous, while the Darcy velocity is discontinuous.Similar results are shown by Figure 21, however, the Darcy velocity in the y direction is continuous as the material interface happens to be perpendicular to the y-axis.
The comparison of the hydraulic head distribution simulated by the proposed NMM model and Zhou's FEM model with same mesh size is shown in Figure 19.The results show that the hydraulic head is continuous across the material interfaces AB and CD.fixed with 1m on the left side and 0 m on the right side, while the top boundary has a specified inflow flux of 1 m/d.The bottom of the model is impermeable.Due to the dual cover systems, the NMM mesh (Figure 18b) does not need to conform with the material interfaces, unlike the conventional FEM mesh (Figure 18c).In order to study the influence of the material interface, numerical results with two different mesh densities along profiles at x = 100 m and y = 80 m are plotted in Figures 20 and 21.As shown in Figure 20, the hydraulic head across the material interface is continuous, while the Darcy velocity is discontinuous.Similar results are shown by Figure 21, however, the Darcy velocity in the y direction is continuous as the material interface happens to be perpendicular to the y-axis.In order to study the influence of the material interface, numerical results with two different mesh densities along profiles at x = 100 m and y = 80 m are plotted in Figures 20 and 21.As shown in Figure 20, the hydraulic head across the material interface is continuous, while the Darcy velocity is discontinuous.Similar results are shown by Figure 21, however, the Darcy velocity in the y direction is continuous as the material interface happens to be perpendicular to the y-axis.
In order to study the influence of the material interface, numerical results with two different mesh densities along profiles at x = 100 m and y = 80 m are plotted in Figures 20 and 21.As shown in Figure 20, the hydraulic head across the material interface is continuous, while the Darcy velocity is discontinuous.Similar results are shown by Figure 21, however, the Darcy velocity in the y direction is continuous as the material interface happens to be perpendicular to the y-axis.For the purpose of verifying the accuracy of modeling refraction law across material interfaces, the velocity components of element nodes on the different sides of the material interface AB are plotted in Figures 22 and 23.The results show that the requirement of the refraction law cannot be strictly met using the conventional NMM model (Figure 22).Thus, the numerical results may be unreliable.As shown in Figure 23, the velocity components solved by the proposed model fully satisfy the refraction law, as well as Zhou's model [14].For the purpose of verifying the accuracy of modeling refraction law across material interfaces, the velocity components of element nodes on the different sides of the material interface AB are plotted in Figures 22 and 23.The results show that the requirement of the refraction law cannot be strictly met using the conventional NMM model (Figure 22).Thus, the numerical results may be unreliable.As shown in Figure 23, the velocity components solved by the proposed model fully satisfy the refraction law, as well as Zhou's model [14].
As a result, the proposed NMM model can achieve the same precision of the Darcy velocity as Zhou's method without iterations, while the dual cover systems of NMM model are more flexible and convenient without the need of conforming to material interfaces.Moreover, the new model is proven to be insensitive to mesh size, so that reliable results can be achieved through a sparse mesh.Therefore, the proposed NMM model is very suitable for modeling the Darcy flow in heterogeneous porous media.
For the purpose of verifying the accuracy of modeling refraction law across material interfaces, the velocity components of element nodes on the different sides of the material interface AB are plotted in Figures 22 and 23.The results show that the requirement of the refraction law cannot be strictly met using the conventional NMM model (Figure 22).Thus, the numerical results may be unreliable.As shown in Figure 23, the velocity components solved by the proposed model fully satisfy the refraction law, as well as Zhou's model [14].For the purpose of verifying the accuracy of modeling refraction law across material interfaces, the velocity components of element nodes on the different sides of the material interface AB are plotted in Figures 22 and 23.The results show that the requirement of the refraction law cannot be strictly met using the conventional NMM model (Figure 22).Thus, the numerical results may be unreliable.As shown in Figure 23, the velocity components solved by the proposed model fully satisfy the refraction law, as well as Zhou's model [14].As a result, the proposed NMM model can achieve the same precision of the Darcy velocity as Zhou's method without iterations, while the dual cover systems of NMM model are more flexible and convenient without the need of conforming to material interfaces.Moreover, the new model is proven to be insensitive to mesh size, so that reliable results can be achieved through a sparse mesh.Therefore, the proposed NMM model is very suitable for modeling the Darcy flow in heterogeneous

Discussion and Conclusions
In this paper, a new high-order NMM model for the Darcy flow in heterogeneous porous media is presented.Due to the lack precision of measuring the velocity solution by ordinary methods, a new weight function with a continuous nodal gradient was adopted into the original NMM model.In order to achieve best precision when modeling material interfaces, the refraction law is fully introduced into the NMM model as an a posteriori requirement.Compared to other numerical schemes, the proposed NMM model has three major advantages as follows: (1) the solution of Darcy velocity is accurate, and continuous at each node; (2) the material interfaces exactly satisfy the refraction law, and thus the simulation results can be accurate and reliable; and (3) the pre-processing work is efficient, due to the two independent cover systems of NMM, in contrast with other numerical schemes.Especially for complex heterogeneous problems even with intersecting material interfaces, the proposed model provides a significant advantage.
Additionally, a set of numerical examples are modeled to verify the accuracy, efficiency, and robustness of the proposed high-order numerical manifold method.First, Darcy flow in homogeneous porous media with different boundary conditions is simulated.The numerical results show good precision compared with analytical solutions, which is the foundation for modeling heterogeneous problems.Then, a continuous-heterogeneous example with hydraulic conductivity varying smoothly in the problem domain is demonstrated.Compared with the analytical solution, the developed method shows the ability to model heterogeneous problem without material interfaces.Moreover, convergence analysis shows that the new model has better accuracy for the Darcy velocity solution, especially in heterogeneous problems.Finally, the new high-order NMM model is applied into a heterogeneous problem with material interfaces and different boundary conditions.Results show that refraction law is exactly satisfied by the proposed model, while the computational cost, especially the mesh effort, is less than that of the existing models.
Furthermore, the proposed high-order NMM model can be extended to 3D cases and applied into fluid flow and contaminant transportation in heterogeneous porous media.

Figure 1 .
Figure 1.Finite element mesh and physical domain.

Figure 1 .
Figure 1.Finite element mesh and physical domain.

Figure 1 .
Figure 1.Finite element mesh and physical domain.

Figure 2 .
Figure 2. The dual cover systems: mathematical covers and physical covers.

Figure 3 .
Figure 3. Manifold elements generated by physical covers.

Figure 3 .
Figure 3. Manifold elements generated by physical covers.

Figure 4 .
Figure 4. Physical domain and boundaries of a Darcy flow problem.

Figure 4 .
Figure 4. Physical domain and boundaries of a Darcy flow problem.

Figure 5 .
Figure 5. Refraction law across the material interface.

Figure 7 .Figure 8 .
Figure 7. Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 10-layer mesh.

Figure 7 .
Figure 7. Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 10-layer mesh.
are the comparison of the Darcy velocity component v x and v y , respectively.It is indicated that the accuracy of the new NMM model is close to the analytical solution without a refined mesh.

Figure 7 .Figure 8 .
Figure 7. Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 10-layer mesh.

Figure 9 .
Figure 9.The comparison of solutions computed by the proposed numerical manifold method (NMM) model using a 10-layer mesh with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 8 .Figure 8 .
Figure 8. Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, (c) Darcy velocity component in the y direction using the conventional NMM model with a 10-layer mesh.

Figure 9 .
Figure 9.The comparison of solutions computed by the proposed numerical manifold method (NMM) model using a 10-layer mesh with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 9 .
Figure 9.The comparison of solutions computed by the proposed numerical manifold method (NMM) model using a 10-layer mesh with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 10 .
Figure 10.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.(NDF: number of degrees of freedom).

Figure 10 .
Figure 10.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.(NDF: number of degrees of freedom).

Figure 11 .
Figure 11.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the conventional NMM model with a 60-layers mesh.

Figure 12 .
Figure 12.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 11 .
Figure 11.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the conventional NMM model with a 60-layers mesh.

Figure 11 .
Figure 11.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the conventional NMM model with a 60-layers mesh.

Figure 12 .
Figure 12.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 12 .
Figure 12.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 0.3 m: the (a) hydraulic head, (b) Darcy velocity in the x direction, and (c) Darcy velocity in the y direction.

Figure 13 .
Figure 13.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Figure 14 .
Figure 14.A continuous-heterogeneous model for Darcy flow: model geometry and boundary conditions.

2 Figure 13 .
Figure 13.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Figure 13 .
Figure 13.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Figure 14 .
Figure 14.A continuous-heterogeneous model for Darcy flow: model geometry and boundary conditions.

Figure 15 .
Figure 15.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 60-layer mesh.

Figure 16 .
Figure 16.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at y = 0 m: the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction.

Figure 17 .
Figure 17.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: the (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Figure 15 .Figure 15 .
Figure 15.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 60-layer mesh.

Figure 16 .
Figure 16.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at y = 0 m: the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction.

Figure 17 .
Figure 17.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: the (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Figure 16 .
Figure 16.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at y = 0 m: the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction.

Figure 15 .Figure 16 .
Figure 15.Computed contour map of the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction using the proposed NMM model with a 60-layer mesh.

Figure 17 .
Figure 17.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: the (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

4. 3 .Figure 17 .
Figure 17.The comparison of convergence curves in the L 2 norm by the proposed NMM model and the conventional NMM model: the (a) the relative error of the hydraulic head and (b) the relative error of the Darcy velocity.

Processes 2018, 6 ,Figure 18 .Figure 19 .
Figure 18.A discontinuous-heterogeneous model with material interfaces: (a) model geometry and boundary conditions, (b) non-conforming mesh based on the proposed NMM model, and (c) conforming mesh based on Zhou's FEM model [14].

Figure 20 .Figure 18 .
Figure 20.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 100 m: the (a) hydraulic head, (b)

Figure 18 .Figure 19 .
Figure 18.A discontinuous-heterogeneous model with material interfaces: (a) model geometry and boundary conditions, (b) non-conforming mesh based on the proposed NMM model, and (c) conforming mesh based on Zhou's FEM model [14].

Figure 19 .
Figure 19.Numerical results of hydraulic head distribution simulated by (a) the proposed NMM model and (b) Zhou's FEM model [14].

Figure 20 .
Figure 20.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 100 m: the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction.

Figure 20 .Figure 21 .
Figure 20.The comparison of solutions computed by the proposed NMM model using different mesh densities with an analytical solution along a profile located at x = 100 m: the (a) hydraulic head, (b) Darcy velocity component in the x direction, and (c) Darcy velocity component in the y direction.Processes 2018, 6, x FOR PEER REVIEW 17 of 20

Figure 22 .Figure 21 .
Figure 22.Verification of the refraction law on the material interface using the conventional NMM model: the (a) normal component on AB, (b) tangential component on AB, (c) normal component on

Figure 22 .Figure 22 .
Figure 22.Verification of the refraction law on the material interface using the conventional NMM model: the (a) normal component on AB, (b) tangential component on AB, (c) normal component on CD, and (d) tangential component on CD.

Figure 22 .
Figure 22.Verification of the refraction law on the material interface using the conventional NMM model: the (a) normal component on AB, (b) tangential component on AB, (c) normal component on CD, and (d) tangential component on CD.

Figure 23 .
Figure 23.Verification of the refraction law on the material interface using the conventional NMM model: the (a) normal component on AB, (b) tangential component on AB, (c) normal component on CD, and (d) tangential component on CD.

Figure 23 .
Figure 23.Verification of the refraction law on the material interface using the conventional NMM model: the (a) normal component on AB, (b) tangential component on AB, (c) normal component on CD, and (d) tangential component on CD.
In this way, the discontinuity of the triangular element cut by boundaries or material interface is presented by NMM.

Table 1 .
Error in the L 2 norm of hydraulic head and Darcy velocity by the proposed NMM model.

2 Error of Hydraulic Head L 2 Error of Darcy Velocity Relative Error Convergence Rate Relative Error Convergence Rate
1kh means the half number of mesh layers in the vertical direction.

Table 2 .
Error in the L 2 norm of hydraulic head and Darcy velocity by the convention NMM model.
kh L

Table 1 .
Error in the L 2 norm of hydraulic head and Darcy velocity by the proposed NMM model.
1kh means the half number of mesh layers in the vertical direction.

Table 2 .
Error in the L 2 norm of hydraulic head and Darcy velocity by the convention NMM model.