Differential Quadrature Solution for One-dimensional Aquifer Flow

Differential Quadrature Method (DQM) has been applied to the solution of aquifer flow problems. Three examples from of each of the three one-dimensional aquifer flow equation problems, a confined aquifer flow with time dependent boundary conditions, a composite confined aquifer and an unconfined aquifer with seepage, were examined. The results of DQM solution were then compared with the results obtained from analytical solution, the Explicit Finite Differences Method and Implicit Finite Differences Method. Based on the comparison results, it was concluded that the DQM provides similar results but with relatively faster calculation speed, less nodes and memory usage.


INTRODUCTION
The one-dimensional unsteady flow equation for groundwater flow is a parabolic differential equation, and typically, Finite Differences Methods (FDM) and Finite Elements Method (FEM) are used for the numerical solution.In the solution of differential equations, DQM can be used as an alternative to these conventional methods.In the homogeneous or heterogeneous soil, confined or unconfined aquifer problems are common in Hydrogeology, Civil Engineering, Irrigation and Drainage Engineering.During recent years, the studies and research on the aquifer problems related to the initial and boundary conditions have been successfully carried out; meanwhile, many computational method and techniques have been developed for numerical solution of the governing equations.
Lockington [17] has studied the response of unconfined aquifer to sudden change in boundary head.In his paper, the analytical approximations to the solution of the one-dimensional Boussinesq equation were obtained using a weighted residual method.Wang and Anderson [32] have studied several groundwater flow problems, in a finite aquifer with recharge boundary, by using an explicit finite difference method for numerical solution.Onder [22] has found an analytical solution for one of the problems that Wang and Anderson had, by examining the flow resulting from a sudden rise or decline in the water stage of a flood channel in a composite aquifer.As "the case (1)", this problem was examined in the present study.
Homogeneous semi-infinite aquifer problems have been examined several times in the past, for instance, the solution for drawdown resulting from a step change in the drain is available in the literature [18,20,22,24].Non-uniform aquifer problem is "the case (2)" in the present study.
Bear [2] has addressed steady state flow in heterogeneous aquifers; Boonstra and Boehmer [5] have generated an approximate solution to the problem of groundwater flow in a composite dike aquifer system for the case that a well in the dike is pumped at a constant rate.Singh [31] has developed a semi-analytical model for computing drawdown in and around a partially penetrating large diameter well.Butter and Liu [6] have presented a semi-analytical solution for the analysis of drawdown data obtained from pump test performed in non-uniform aquifers, which represent a linear strip case.Guo [10] has examined a transient ground water flow between reservoirs and water table aquifers, and has given an analytical solution method.Ostfeld et al. [23] have proposed a general analytical solution scheme for determining groundwater levels for channel/group-water systems with recharge, and, use of Laplace transform method to solve a linearized form of the Boussinesq equation.A solution for leaky onedimensional flow problem with storage and skin effect in finite-width sink was offered [21].Li et al. [14] has presented the confined-unconfined flow in a horizontal confined aquifer.Surface water / groundwater interactions have been studied by several researchers using different techniques and approaches [15,19,23].As the last numerical example, surface and groundwater interact in an unconfined aquifer is "the case (3)" in this present study.
Also, Finite Element Method (FEM) has found a wide range of applications in groundwater investigations.Khebchareon and Saenton [13] used Crank-Nicolson and Galerkin Finite Element Method for 1D groundwater.Liou and Yeh [16] investigated a one-dimensional groundwater transport equation with two uncertain parameters, groundwater velocity and longitudinal dispersivity.A comparison, between these FEM solutions and DQM solutions that were obtained from the examples examined in the present study, was performed.

DIFERENTIAL QUADRATURE METHOD
The DQM, providing a solution to the differential equations of any systems, was developed first by Richard Bellman [3].Afterwards, that was developed by Shu [27].DQM is an alternative approach to the standard methods such as finite difference and finite elements, for the initial value and boundary conditions encountered in physics and mathematics [28].In the DQM, a partial derivative of a function with respect to a space variable at a discrete point is approximated as a weighted linear sum of the function values at all discrete points in the region of that variable.The method can be written as, where x j are the discrete points in the variable, ) (r x u is the r th order derivative of the function, u(x j ) are the function values at these points, and A (r)  ij are the weight coefficients for the r th order derivative of the function.Considering the weight coefficients, Shu and Xue [30], made important studies, and they proposed some solutions in their studies.Determining the weight coefficients is the most crucial step in use of DQM.For calculation of the weight coefficients, a function is chosen.According to the chosen function, the method can take different names.In the Polynomial Differential Quadrature (PDQ), the function is approximated by a (N-1) th degree polynomial; in the Fourier Expansion Base Differential Quadrature (FDQ), the function is approximated by a Fourier series expansion; and another way to determine the weighting coefficients is to employ harmonic functions named the Harmonic Differential Quadrature (HDQ) [27,8].Civalek [8] had compared the methods of differential quadrature (DQ) and harmonic differential quadrature (HDQ) and, in his other studies, DQ and HDQ methods are presented for buckling, bending, and free vibration analysis of thin isotropic plates and columns [33,34].
At some boundary value problems, DQM performance is highly dependent on the boundary conditions and sampling grid points.The overall sensitivity of the model especially depends on the location and number sampling grid points.However, the boundary conditions can be easily implemented to DQ system.The boundary conditions, which are usually Dirichlet, Neuman and/or mixed type function, do not create any difficulty in this implementation process [4,7].Civalek [7] implies that the determination of the effective choice for any problem reduces the analysis time.For instance, previous studies show that, in the problems that have linear equations and homogeneous boundary conditions, selecting equal intervals are adequate for solution.In vibration problems, the choice of grid points through the Chebyshev-Gauss-Lobatto method is more reasonable [7].Similarly, in the water wave propagation and diffusion problems, the method can be used [11,12].

DQM IN GROUNDWATER FLOW APPLICATIONS
The basic equations of groundwater flow are obtained by means of mass conservation principle combined with Darcy law [2].In the one-dimensional unsteady flow, definition of mass conservation and use of notation in Fig. 1, Eq. ( 2) can be written.
In this equation; T=K.B is transmissivity capacity; h, h 1 and h 2 define the piezometric levels at main aquifer, unconfined aquifer and bottom aquifer respectively; K, K 1 and K 2 are permeability coefficients in main aquifer and semipervious stratums; B, B 1 and B 2 are width in main aquifer and semipervious stratums.According to h, h 1 and h 2 piezometric levels the leaking can be to inside or outside.In the unconfined aquifer as Fig. 2, use of u=h 2 transformation Eq. ( 2) is linearized according to h 2 , Eq. ( 3) can be written.
on condition that h is mean piezometric level, q S is seepage from soil surface.

Fig. 2. Leaky and unconfined aquifer
Determining the boundary conditions, for example, in the condition of h(x,0)=h i,1 =f(x), h(0,t)=h 1,s =g(x) and h(L,t)=h N,s ,=f(t), the Eq.2 and Eq.3, can be rewritten for the solution of DQM as seen in Eq. 4 and Eq.5  .In the present study, Polynomial Differential Quadrature Method (PDQM) was used and the coefficients of weight matrices were calculated by use of Quan and Chang's Approach [11,29,25,26].In the distribution, grid points are more frequently near boundaries (non-equally).Chebyshev-Gauss-Lobatto grid points distribution is more appropriate in the time dependent problems [7].Therefore, in the DQM solutions, Chebyshev-Gauss-Lobatto grid-points distribution was used, because the problem is depended on the time.

NUMERICAL SOLUTIONS
For the one-dimensional aquifer flow equation problems three examples (a confined aquifer flow with time dependent boundary conditions, a composite confined aquifer and an unconfined aquifer with seepage) were solved.The results obtained from the DQM solutions were compared with the results of Explicit Finite Differences Method (EFDM) and Implicit Finite Differences Method (IFDM).

Case 1
A confined aquifer having independent or time dependent boundary conditions in Fig. 3 was selected as Case 1.This aquifer has a transmissivity capacity (T) of 0.02 m 2 /min; a storage coefficient (S) of 0.002; and length (L) of 100 m [32].The initial and boundary conditions for the solution of Equation 7which is written for 1D flow are h(x, 0)=16m, h(0,t)=16 m and h(L,t) =11 m.The same example was solved use of analytical method and numerical method by Onder [22], and, the results were compared.
In this study, the piezometric level changes were calculated using DQM; which were compared with the EFDM, IFDM and analytical solution use of given by Onder [22].In the numerical solutions, different number of grid points were used.The differences between analytical solution and numerical solutions are calculated, and in the different times, Root Mean Square Error values for DQM, IFDM and EFMD are given in Table 1.For different number of grid points, e i error values and RMS Error can be written as follows: where h is the difference between analytical solution and numerical solution.Results of DQM, IFDM and EFDM solutions are given in Fig. 4 for some number of grid points.In the EFDM, some unstable results can occur depending on the soil hydraulic properties and the size of the spatial and temporal mesh.For the stable solution, CFL (Courant-Friedrichs-Lewy) condition may be agreed [9].Thus, some results can't be obtained in EFDM.accordingly.As seen in the table, results of DQM convergence rapidly as increasing N t values.DQM solution use of few number of grid points, the results are obtain closer to analytical solution.Also, RMS Error values at DQM are smaller than the other numerical methods.Furthermore, the less grid points are used in DQM.

Case 2
A heterogeneous semi infinite confined aquifer in Fig. 5 was selected as the case 2. The aquifer has two regions at the different soil properties.Left boundary condition is a sudden rise from 28 m to 25 m.The aquifer set a limit to impervious stratums at over and under.Consequently, in Eq(4) K 1 and K 2 is zero.The composite aquifer has the following properties: Transmissivity of region (1) T 1 = 800 m 2 /day and of region (2) T 2 =200 m 2 /day, storage coefficient of region (1) S 1 =0.0004, and of region (2) S 2 =0.0004.

Fig. 5. Composite aquifer
This problem is solved use of DQM, EFDM and IFDM and results are compared to analytical solution give by Onder [22].In Fig. 6, error values are given for some grid points.RMS Error values are seen in the Table 2.When the results are compare to analytical solution, it's seen that DQM results use of less grid points are very close to analytical solution than the other numerical methods.In the case as case 1, its seen that results of DQM convergence rapidly.

Case 3
In the case 3, an unconfined aquifer between two drainage canals was investigated (Fig. 7).From surface soil to the aquifer a seepage exist (q S >0).h 0 initial boundary value is 1 m, in the drainage canals water levels is 1m, and water level is a sudden rise to 2 m in the left canal boundary condition.Also, the seepage value from surface soil to the aquifer is supported q S =0.02 m/day [23].
In this study, this example was solved for a sudden rise in the drainage canal conditions use of DQM, IFDM and EFDM.Results of DQM, IFDM and EFDM are given in the Fig. 8 for different number of gird points.As seen in Case 1 and 2 that DQM results are closer to analytical solution than IFDM and EFDM.Therefore, the RMS differences between all numerical methods and DQM (N x =101, N t =51) are given in Table 3.

CONCLUSION
A DQM approach has been used first time in the groundwater hydraulics in this study.DQM has found increasing use in recent years in Hydraulic Engineering, because it is an alternative approach to the conservative methods.From the previous applications of DQM, it is seen that the results of DQM are converged rapidly and closer to analytical solutions than other numerical solutions [12].
In this study, with different boundary conditions, one dimensional groundwater problems were solved.The DQM results were then compared with analytical solution and the results obtained from IFDM and EFDM.The study has shown that the DQM is quickly converged and closer to analytical solution than the other methods.In the case 1, the results of DQM at 11x11 grid points are similar to results at 101x1001 grid points in IFDM and at 11x101 grid points in EFDM.In the case 2, the results of DQM at 51x51 are more closely to analytical solution results than 101x1001 grid points in IFDM and EFDM.The use of relatively fewer nodes with an acceptable accuracy is an advantage of the method, in order to save compile time and memory usage.
In the application of the method, choosing grid points and test function are important step to facilitate the stability of numerical solution.Parametric preferences in this study are appropriate for groundwater flow problems.

Fig. 4 .
Fig. 4. Results of DQM, IFDM and EFDM solution for some number of grid points

Table 1 .
RMS Error values

Table 2 .
RMS Error values for heterogeneous semi infinite confined aquifer problem.