ABSTRACT

We investigated the influence of fictitious boundary distance, a parameter of MFS, to determine piezometric levels of two unconfined sedimentary aquifers assuming Dupuit-Forchheimer and steady-state flow hypothesis. Two study areas were modelled: Guariroba’s Environmental Protection Area, in Mato Grosso do Sul State, Brazil, and Juazeiro do Norte City, in Ceará State, Brazil. It was observed that in order to use the MFS as a numerical method in modeling groundwater flow, it is necessary to determine the best distance value of the fictitious boundary. This value can be chosen from the use of field data within the analyzed domain, where the relative error is a parameter to be minimized. Applying this methodology and comparing with the results of the MODFLOW application for the same set of initial data, we concluded that the MSF allows to estimate the piezometric level values within the analyzed domains and that the results of the statistical comparison between them point to the need to investigate the representativeness of both methods to determine which one is most appropriate for modelling the groundwater flow in each region.

Keywords:
Method of fundamental solutions; MFS; Ground-water flow; Numerical simulation; Meshless

RESUMO

Neste trabalho, investigou-se a influência da distância da fronteira fictícia, um dos parâmetros do MSF, para determinar o nível piezométrico de um aquífero sedimentar não confinado sob a adoção das hipóteses de Dupuit-Forchheimer e fluxo em regime permanente. Foram modeladas duas áreas de estudo: a APA do Córrego do Guariroba, no Mato Grosso do Sul, e a cidade de Juazeiro do Norte, no Ceará. Observou-se que para usar o MSF como método numérico na modelagem do fluxo de água subterrâneo é necessário determinar o melhor valor de distância da fronteira fictícia. Esse valor pode ser escolhido a partir da utilização de dados de campo dentro do domínio de análise, onde o erro relativo é um parâmetro a ser minimizado. Aplicando essa metodologia e comparando com os resultados da aplicação do MODFLOW para o mesmo conjunto de dados iniciais, concluímos que o MSF permite estimar os valores de nível piezométrico dentro dos domínios analisados e que os resultados da comparação estatística entre os dois modelos apontam para a necessidade de investigar a representatividade de ambos os métodos para determinar qual deles é mais apropriado para modelar o fluxo de água subterrâneo de cada região.

Palavras-chave:
Método das soluções fundamentais; MSF; Fluxo de água subterrâneo; Simulação numérica; Método livre de malha

INTRODUCTION

After development of digital computers, traditional finite difference method (FDM) and finite element method (FEM) became popular among researches and were applied to approximate functional values to discretized systems in some manifold physical problems (McDonald & Harbaugh, 2003McDonald, M. G., & Harbaugh, A. W. (2003). The history of MODFLOW. Ground Water, 41(2), 280-283. PMid:12656294. http://dx.doi.org/10.1111/j.1745-6584.2003.tb02591.x.
http://dx.doi.org/10.1111/j.1745-6584.20... ). Extensive amount of scientific work (theoretical and applied) generated in 1960s and 1970s provided background which allowed the construction of corporative softwares, such MODFLOW and FEFLOW. In that same time, other numerical techniques were investigated with less enthusiasm, such as Boundary Element Methods and other meshless methods (Brebbia et al., 1984Brebbia, C. A., Telles, J. C. F., & Wrobel, L. C. (1984). Boundary element techniques. Berlin: Springer-Verlag. http://dx.doi.org/10.1007/978-3-642-48860-3.
http://dx.doi.org/10.1007/978-3-642-4886... ; McDonald & Harbaugh, 2003McDonald, M. G., & Harbaugh, A. W. (2003). The history of MODFLOW. Ground Water, 41(2), 280-283. PMid:12656294. http://dx.doi.org/10.1111/j.1745-6584.2003.tb02591.x.
http://dx.doi.org/10.1111/j.1745-6584.20... ; Trefry & Muffels, 2007Trefry, M. G., & Muffels, C. (2007). FEFLOW: a finite-element ground water flow and transport modelling tool. Ground Water, 45(5), 525-528. http://dx.doi.org/10.1111/j.1745-6584.2007.00358.x.
http://dx.doi.org/10.1111/j.1745-6584.20... ).

In late 1970s, Mathon & Johnston (1977)Mathon, R., & Johnston, R. L. (1977). The aproximate solution of elliptic boundary-value problems by fundamental solutions. SIAM Journal on Numerical Analysis, 14(4), 638-650. http://dx.doi.org/10.1137/0714043.
http://dx.doi.org/10.1137/0714043... proposed the computational implementation of the Method of Fundamental Solutions and, in 1990s, Golberg & Chen (1998)Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press. outlined alternatives to apply the MFS to inhomogeneous and time-dependent problems, which was, until that moment, its major barrier (Chen et al., 2008Chen, C. S., Karageorghis, A., & Smylirs, Y. S. (2008). The method of fundamental solutions: a meshless method. Atlanta: Dynamic Publishers.). Since then, the Method of Fundamental Solutions (MFS) was analyzed under its particularities, such the needing of a fictitious boundary, the conditioning of matrixes, location of external points and alternatives to improve its limitations, providing an excellent background to its application in different physical phenomena (Gu et al., 2019Gu, Y., Fan, C. M., & Xu, R. P. (2019). Localized method of fundamental solutions for large-scale modeling of two-dimensional elasticity problems. Applied Mathematics Letters, 93, 8-14. http://dx.doi.org/10.1016/j.aml.2019.01.035.
http://dx.doi.org/10.1016/j.aml.2019.01.... ; Oh et al., 2019Oh, J., Zhu, H., & Fu, Z. (2019). An adaptive method of fundamental solutions for solving the Laplace equation. Computers & Mathematics with Applications, 77(7), 1828-1840. http://dx.doi.org/10.1016/j.camwa.2018.11.021.
http://dx.doi.org/10.1016/j.camwa.2018.1... ; Liu & Sarler, 2013Liu, Q. G., & Sarler, B. (2013). Non-singular method of fundamental solutions for two-dimensional isotropic elasticity problems. Computer Modeling in Engineering & Sciences, 91(4), 235-266.; Liu, 2012Liu, C. S. (2012). An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Engineering Analysis with Boundary Elements, 36(8), 1235-1245. http://dx.doi.org/10.1016/j.enganabound.2012.03.001.
http://dx.doi.org/10.1016/j.enganabound.... ; Barrero-Gil, 2012Barrero-Gil, A. (2012). The method of fundamental solutions without fictitious boundary for solving Stokes problems. Computers & Fluids, 62, 86-90. http://dx.doi.org/10.1016/j.compfluid.2012.03.024.
http://dx.doi.org/10.1016/j.compfluid.20... ; Chen & Wang, 2010Chen, W., & Wang, F. Z. (2010). A method of fundamental solutions without fictious boundary. Engineering Analysis with Boundary Elements, 34(5), 530-532. http://dx.doi.org/10.1016/j.enganabound.2009.12.002.
http://dx.doi.org/10.1016/j.enganabound.... ; Young et al., 2008Young, D. L., Tsai, C. C., Chen, C. W., & Fan, C. M. (2008). The method of fundamental solutions and condition number analysis for inverse problems of Laplace equation. Computers & Mathematics with Applications, 55(6), 1189-1200. http://dx.doi.org/10.1016/j.camwa.2007.05.015.
http://dx.doi.org/10.1016/j.camwa.2007.0... ; Alves & Chen, 2005Alves, C. J. S., & Chen, C. S. (2005). A new method of fundamental solutions applied to nonhomogeneous elliptic problems. Advances in Computational Mathematics, 23(1-2), 125-142. http://dx.doi.org/10.1007/s10444-004-1833-5.
http://dx.doi.org/10.1007/s10444-004-183... ; Golberg et al., 1999Golberg, M. A., Chen, C. S., & Muleshkov, A. S. (1999). The method of fundamental solutions for time-dependent problems. Transactions on Modelling and Simulation, 22, 377-386.; Golberg & Chen, 1998Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press.; Bogomolny, 1985Bogomolny, A. (1985). Fundamental solutions method for elliptic boundary value problems. SIAM Journal on Numerical Analysis, 22(4), 644-669. http://dx.doi.org/10.1137/0722040.
http://dx.doi.org/10.1137/0722040... ). Although the MFS was developed in 1960s (Chen et al., 2008Chen, C. S., Karageorghis, A., & Smylirs, Y. S. (2008). The method of fundamental solutions: a meshless method. Atlanta: Dynamic Publishers.), its implementation was mostly oriented to theoretical analysis and other areas than groundwater flow, mainly because the difficult to implement time-dependent problems and, consequently, transient states, as well as the implementation of collocation points inside the domain (Wang & Zheng, 2016Wang, F. Z., & Zheng, K. H. (2016). The method of fundamental solutions for steady-state groundwater flow problems. Zhongguo Gongcheng Xuekan, 39(2), 236-242. http://dx.doi.org/10.1080/02533839.2015.1082936.
http://dx.doi.org/10.1080/02533839.2015.... ).

Under the Dupuit-Forchheimer, water incompressibility, isotropy of the sedimentary aquifer and steady-state flow assumptions, groundwater flow can be described by Laplace differential operator (Haitjema, 1995Haitjema, H. M. (1995). Analytic element modeling of groundwater flow. London: Academic Press.; Wang & Zheng, 2016Wang, F. Z., & Zheng, K. H. (2016). The method of fundamental solutions for steady-state groundwater flow problems. Zhongguo Gongcheng Xuekan, 39(2), 236-242. http://dx.doi.org/10.1080/02533839.2015.1082936.
http://dx.doi.org/10.1080/02533839.2015.... ). This allows the application of the MFS to compute numerically the piezometric value of points inside a domain, what suggests its availability to be used as a numerical tool in groundwater flow modelling. The application of MFS to groundwater study incorporates the advantages of a meshless method and provides an alternative to situations where traditional models present difficulties to solve (Liu & Gu, 2005Liu, G. R., & Gu, Y. T. (2005). An introduction to meshless methods and their programming. Netherlands: Springer.). Thereby, the objective of this work is to apply the MFS in a real groundwater flow examples to evaluate its behavior and results.

Two study areas were simulated, the first is located in Brazil’s Northeast region (Juazeiro do Norte city) and second in Midwest region (Environmental Protected Area of Guariroba river). To reach this objective, we compared the numerical values obtained with field data and MODFLOW results to the same areas. MODFLOW result was used for comparison purposes in function of its availability as a free-to-use and open source software and its reliability acquired as a groundwater modelling tool among researches around the world. The main contribution of this investigation is to provide the initial steps of groundwater modelling with the MFS, a meshless method, and to verify the influence of its main characteristics, the need of a fictitious boundary to perform such task.

MATERIAL AND METHODS

Method of Fundamental Solutions (MFS) applied to ground-water flow

Groundwater flow can be described by Laplace equation in two dimensions under Dupuit-Forchheimer assumptions: water incompressibility, isotropy of the sedimentary aquifer and steady-state flow (Haitjema, 1995Haitjema, H. M. (1995). Analytic element modeling of groundwater flow. London: Academic Press.). To identify equipotentials curves for hydraulic head in a domain under those assumptions requires to solve the Laplace differential operator to a set of boundary conditions. These boundary conditions can be of Dirichlet, Neumann, mixed, Robin or Cauchy types, whereas they provide a well-posed problem (Jazayeri & Werner, 2019Jazayeri, A., & Werner, A. D. (2019). Boundary condition nomenclature confusion in groundwater flow modeling. Ground Water, 57(5), 664-668. PMid:31001804. http://dx.doi.org/10.1111/gwat.12893.
http://dx.doi.org/10.1111/gwat.12893... ). Analytical solution to such equations systems is not always easy to determine. In most cases, numerical methods are used to perform an approximation.

The Method of Fundamental Solutions (MFS) is a truly meshless method which uses the fundamental solution of a differential operator as a shape function in a weighted-residual integral, wherein the weighting function is the Dirac’s Delta. This integral formulation is also called Collocation Method and can be summarized as Equations 1 and 2 to MFS (Brebbia & Dominguez, 1992Brebbia, C. A., & Dominguez, J. (1992). Boundary elements: an introductory course. Boston: WIT Press Computational Mechanics Publications.).

where $R\left(x\right)$ is the residual function applied in a collocation point $x\in \text{}\text{}\Gamma \text{}\text{}$, $\text{}\text{}\Gamma \text{}\text{}$ is the boundary of a domain $\text{}\text{}\Omega \text{}\text{}\subset R^{\text{2}}$ and $\delta \left(x,x^{\prime }\right)$ is Dirac’s Delta function in a source point x' applied in a collocation point $x$ (Figure 1).

Since approximated solution $\tilde{u}\left(x\right)$ of a point $x\in \text{}\text{}\Omega \text{}\text{}$ is written as Equation 3, where $F\left(x,x_i\right)$ is the fundamental solution for Laplace differential operator in two-dimensions, ${\lambda }_i$ are constant coefficients and $u\left(x\right)$ is the analytical solution, the relationship in Equation 4 can be outlined. $\text{}\text{}\Gamma \text{}\text{}=\text{}\text{}\text{}\text{}\text{}\text{}\Gamma \text{}\text{}{}_1+\text{}\text{}\Gamma \text{}\text{}{}_2$, $\text{}\text{}\Gamma \text{}\text{}{}_1$ is the part of $\text{}\text{}\Gamma \text{}\text{}$ with Dirichlet boundary condition and $\text{}\text{}\Gamma \text{}\text{}{}_2$, with Neumann boundary condition.

Fundamental solution of Laplace differential operator (Equation 5), where $r$ is the distance vector between $x$ and x', provides an indetermination when coincident source and collocation points are chosen $(r=0)$ during matrixes construction. A fictitious boundary is created to avoid such indetermination (Figure 1). This fictitious boundary is an independent parameter and alternatives to avoid its need were investigated in Chen & Wang (2010)Chen, W., & Wang, F. Z. (2010). A method of fundamental solutions without fictious boundary. Engineering Analysis with Boundary Elements, 34(5), 530-532. http://dx.doi.org/10.1016/j.enganabound.2009.12.002.
http://dx.doi.org/10.1016/j.enganabound.... , Barrero-Gil (2012)Barrero-Gil, A. (2012). The method of fundamental solutions without fictitious boundary for solving Stokes problems. Computers & Fluids, 62, 86-90. http://dx.doi.org/10.1016/j.compfluid.2012.03.024.
http://dx.doi.org/10.1016/j.compfluid.20... and Liu & Sarler (2013)Liu, Q. G., & Sarler, B. (2013). Non-singular method of fundamental solutions for two-dimensional isotropic elasticity problems. Computer Modeling in Engineering & Sciences, 91(4), 235-266..

We obtain Equation 6 writing Equations 3 and 4 in a matrixial notation, where $f$ and $g$ are the blocks of boundary conditions (Equation 4). Algebraically isolating vector $\left[\lambda \right]$ (Equation 7), the numerical values of interested points $s_i\in \text{}\text{}\Omega \text{}\text{}$ can be computed by Equation 8.

The position of the fictitious boundary influences numerical results. For an off-set or a circle around domain, quality of approximation increases as the fictitious boundary moves away until the point that matrixes become ill-conditionate and computers abort operation accusing a non-inversible system. This same result is found as number of source points increases (Golberg & Chen, 1998Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press., Liu, 2012Liu, C. S. (2012). An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Engineering Analysis with Boundary Elements, 36(8), 1235-1245. http://dx.doi.org/10.1016/j.enganabound.2012.03.001.
http://dx.doi.org/10.1016/j.enganabound.... ). Some strategies to circumvent ill-conditioning were considered in Chen et al. (2006)Chen, C. S., Cho, H., & Golberg, M. A. (2006). Some comments on the ill-conditioning of the method of fundamental solutions. Engineering Analysis with Boundary Elements, 30(5), 405-410. http://dx.doi.org/10.1016/j.enganabound.2006.01.001.
http://dx.doi.org/10.1016/j.enganabound.... , Young et al. (2008)Young, D. L., Tsai, C. C., Chen, C. W., & Fan, C. M. (2008). The method of fundamental solutions and condition number analysis for inverse problems of Laplace equation. Computers & Mathematics with Applications, 55(6), 1189-1200. http://dx.doi.org/10.1016/j.camwa.2007.05.015.
http://dx.doi.org/10.1016/j.camwa.2007.0... , Chen et al. (2016)Chen, C. S., Karageorghis, A., & Li, Y. (2016). On choosing the location of the sources in the MFS. Numer Algo, 72(1), 107-130. http://dx.doi.org/10.1007/s11075-015-0036-0.
http://dx.doi.org/10.1007/s11075-015-003... and Hematiyan et al. (2016)Hematiyan, M., Arezou, M., & Haghighi, A. (2016). Reducing the condition number of the coefficient matrix in the method of fundamental solutions for 2D Laplace equation. In Proceedings of 17th Advances in Boundary Element & Meshless Techniques. Ankara, Turkey..

The Method of Fundamental Solutions cannot be applied directly to inhomogeneous problems related to Laplace operator. To solve an inhomogeneous case, a particular solution must be integrated to remove inhomogeneity from boundary, and later superposed. Basic formulation of MFS to Laplace differential operator allows only the implementation of time-independent problems. The implementation of time-dependent systems involves the application of specific mathematical strategies to each situation and to each type of boundary data (Fairweather & Karageorghis, 1998Fairweather, G., & Karageorghis, A. (1998). The method of fundamental solutions for elliptic boundary value problems. Advances in Computational Mathematics, 9(9), 69-95. http://dx.doi.org/10.1023/A:1018981221740.
http://dx.doi.org/10.1023/A:101898122174... ; Wang & Zheng, 2016Wang, F. Z., & Zheng, K. H. (2016). The method of fundamental solutions for steady-state groundwater flow problems. Zhongguo Gongcheng Xuekan, 39(2), 236-242. http://dx.doi.org/10.1080/02533839.2015.1082936.
http://dx.doi.org/10.1080/02533839.2015.... ). Some works performed previous approaches and contributed to application of the MFS to diverse physical phenomena, such heat transference, diffusion problems, time-dependent models and theoretical analysis (Alves & Chen, 2005Alves, C. J. S., & Chen, C. S. (2005). A new method of fundamental solutions applied to nonhomogeneous elliptic problems. Advances in Computational Mathematics, 23(1-2), 125-142. http://dx.doi.org/10.1007/s10444-004-1833-5.
http://dx.doi.org/10.1007/s10444-004-183... ; Golberg et al., 1999Golberg, M. A., Chen, C. S., & Muleshkov, A. S. (1999). The method of fundamental solutions for time-dependent problems. Transactions on Modelling and Simulation, 22, 377-386.; Young et al., 2008Young, D. L., Tsai, C. C., Chen, C. W., & Fan, C. M. (2008). The method of fundamental solutions and condition number analysis for inverse problems of Laplace equation. Computers & Mathematics with Applications, 55(6), 1189-1200. http://dx.doi.org/10.1016/j.camwa.2007.05.015.
http://dx.doi.org/10.1016/j.camwa.2007.0... ). However, application of the MFS in the groundwater flow determination is limited and few works are found in the literature. The difficult to implement time-dependent problems and, consequently, transient states, as well as the implementation of collocation points inside domain (which represent numerical constraints, such lakes, pumping wells, rivers and others) makes the MSF unattractive compared to traditional Finite Difference and Finite Element Methods, which deal with such inhomogeneities with simple adjustments. Nevertheless, MFS provides smoother equipotentials curves near the analyzed boundary, as will be shown later in this work, requests less input parameters than MODFLOW and it does not need the building of a previous mesh to perform a numerical approximation.

To implement the MFS requires a set of data-points, which allows the construction of a boundary around a domain, where differential operator rules the analytical values. Those data-points must have scalar values associated to Dirichlet or Neumann boundary conditions. Since boundary is well-defined and fundamental solution of differential operator is known, Equations 6-8 can be used to compute numerical values. In a two-dimensional analysis of a sedimentary aquifer ruled by Laplace equation, a set of observation wells can determine a boundary and its static piezometric levels are associated to Dirichlet boundary conditions (Jazayeri & Werner, 2019Jazayeri, A., & Werner, A. D. (2019). Boundary condition nomenclature confusion in groundwater flow modeling. Ground Water, 57(5), 664-668. PMid:31001804. http://dx.doi.org/10.1111/gwat.12893.
http://dx.doi.org/10.1111/gwat.12893... ). Assuming that there are no inhomogeneities inside the analyzed area, the MFS is theoretically able to estimate values to any point inside the domain.

Methodology

We compared the capacity of the MFS predicts field hydraulic heads at points inside domain varying off-set distance. Predicted hydraulic heads were compared with predicted values provided by MODFLOW model, which uses Finite Difference Method as numerical tool. Numerical values obtained by MFS with best off-set distance were compared to MODFLOW result through piezometric equipotentials, in order to explicit the difference between them. Flowchart (Figure 2) illustrates the methodology applied in this work.

The MFS was implemented in SCILAB and Java programming languages. Routines were created to apply the method to a variable off-set distance and to a set of boundary points, respectively. MODFLOW model was operated through UFC-FLOW interface.

Since off-set position (working as a fictitious boundary) influences computed values (Golberg & Chen, 1998Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press.; Liu, 2012Liu, C. S. (2012). An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Engineering Analysis with Boundary Elements, 36(8), 1235-1245. http://dx.doi.org/10.1016/j.enganabound.2012.03.001.
http://dx.doi.org/10.1016/j.enganabound.... ), the best off-set distance was determined after its value varied between 100m and 200,000m in steps of 100m. This range was chosen because we verified that those values evidence the representative behavior of MFS applied to this proposed analysis. Then we compared computed values with field data at observation wells location to each step. The best off-set value was chosen as the value which proportionated the minimal relative error between numerical result and field data. Once best off-set distance was determined, it was used to compute piezometric values to a grid of points inside domain. Both generated piezometric equipotentials (by MFS and MODFLOW) were created by an interpolation of values in a 200x200m mesh grid. This mesh resolution was determined in order to draw equipotential curves. Although the MFS does not require a mesh of points, it was performed for comparison purposes.

Two study areas were analyzed. First is located in Middle-west Brazil region in Campo Grande City, an area known as Guariroba’s Environmental Protected Area (Figure 3A) and second, in Northeast region in Juazeiro do Norte City (Figure 3B). They were chosen in function of the availability of observation wells data (static piezometric level) and geological formation reports. Both areas are composed by a set of observations wells (PO) which delimits a polygonal domain (Figure 3). Some wells are located in boundary and others, inside domain. Wells located in boundary were used as boundary points to MFS and MODFLOW models. Domain points were used to compare field measurement of hydraulic heads with computed numerical values to its same coordinates. Indexes 2PO refers to Guariroba’s APA case and 1PO, to Juazeiro do Norte.

Hydraulic heads of Guariroba’s APA wells were obtained in same expedition during an interval less than 24 hours and special care was taken to avoid pumping operations in last 2 days before measurement (Cavazzana et al., 2019Cavazzana, G. H., Lastoria, G., & Gabas, S. G. (2019). Surface-groundwater interaction in unconfined sedimentary aquifer system in the Brazil’s tropical wet region. Revista Brasileira de Recursos Hídricos, 24, e8. http://dx.doi.org/10.1590/2318-0331.241920180136.
http://dx.doi.org/10.1590/2318-0331.2419... ). In Juazeiro do Norte city, hydraulic heads were provided by an annual Hydrological Resources Management Company of Ceará State (Companhia de Gestão de Recursos Hídricos do Estado do Ceará – COGERH) report. To simulate the same situation as Araújo (2016)Araújo, R. S. A. (2016). Metodologias baseadas no mínimo resíduo para a escolha do parâmetro de forma utilizado em RBF’s aplicadas no método assimétrico livre de malhas (MESHLESS) (Tese de doutorado). Centro de Tecnologia, Universidade Federal do Ceará, Fortaleza., Companhia de Gestão dos Recursos Hídricos (2014)Companhia de Gestão dos Recursos Hídricos – COGERH. (2014). Boletim de monitoramento dos poços com DataLogger no Cariri-CE (ano de 2013). Fortaleza: COGERH. report was used, it is related to the measurement period of report’s previous year, 2013. Correspondent watersheds were delimited using TOPODATA digital elevation model (DEM) to allow a representative characterization of study areas. These criteria ensure the implementation of a non-confined aquifer under a time-equivalent measurement period, what prevents simulation of an inconsistent groundwater flow state. Error analysis associated to input data was not performed in this work.

The MFS requests boundary values, its coordinates and fundamental solution of the differential operator, while MODFLOW model requests, besides boundary values and its coordinates, horizontal and vertical permeabilities, storage, porosity and thickness of the permeability layer as input parameters (Akram et al., 2012Akram, F., Rasul, M. G., Khan, M. M. K., & Amir, M. S. I. I. (2012). A comparative view of groundwater flow simulation using two modelling software - MODFLOW and MIKE SHE. In Proceedings of the 18th Australasian Fluid Mechanics Conference. Red Hook: Curran Associates.). Conceptual model adopted to simulate the groundwater flow with MFS was a bidimensional domain, delimited by boundary points (wells) associated to Dirichlet boundary condition (piezometric level). In MODFLOW, one permeable layer was implemented. Parameters values (hydraulic conductivity, specific storage and porosity) were adopted from literature and thickness of modelled layer was extracted from digital elevation model (DEM).

Qualitative analysis of results was performed. We discussed about how main characteristic of MFS - the need of a fictitious boundary - influences numerical approximation in real scenarios (study areas). Results were discussed evocating literature to corroborate and endorse the general behavior of simulated areas. A graphic was plotted to show relative error between numerical values and field data for a variable off-set distance. Best off-set distance was chosen as the value which provided less summation of errors relative to field data. Best off-set distance was used to compute numerical values of a grid inside domain. This allowed the building of piezometric level contours and correspondent raster files.

Statistical analysis was performed over these data. First, we calculated central measure tendency parameters (media, median, standard deviation) for each raster file. Since the collection of surface data does not represent a population with a normal curve distribution, because piezometric values are not aleatory, but obtained from numerical methods, we chose Mann-Whitney U test (non-parametric) to compare results from MFS and MODFLOW models. H0 was stated as both group of values do not differ between themselves, and H1, both group of values do differ between themselves with a level of significance of 5% (α = 0.05, two-tailed) (Umeh & Eriobu, 2016Umeh, E. U., & Eriobu, N. O. (2016). Comparison of two sample tests using both relative efficiency and power of test. Open Journal of Statistics, 6(2), 331-345. http://dx.doi.org/10.4236/ojs.2016.62029.
http://dx.doi.org/10.4236/ojs.2016.62029... ). Mann-Whitney U test is similar to t-test, but it is applicable in situations where the shape of data distribution is not specified. This test can be interpreted as an investigation if first data distribution values is greater than second in an ordinal counting (McKnight & Najab, 2010McKnight, P. E., & Najab, J. (2010). Mann-Whitney U Test. In I. B. Weiner & W. E. Craighead (Eds.), The Corsini encyclopedia of psychology. Hoboken: John Wiley & Sons. http://dx.doi.org/10.1002/9780470479216.corpsy0524.
http://dx.doi.org/10.1002/9780470479216.... ; Pino, 2014Pino, F. A. (2014). A questão da não normalidade: uma revisão. Revista de Economia Agrícola, 61(2), 17-33.; Perme & Manevski, 2019Perme, M. P., & Manevski, D. (2019). Confidence intervals for the Mann-Whitney test. Statistical Methods in Medical Research, 28(12), 3755-3768. PMid:30514179. http://dx.doi.org/10.1177/0962280218814556.
http://dx.doi.org/10.1177/09622802188145... ).

We also generated a difference raster between results of each model for each study area. Root Mean Square Deviation was computed as well as the Kolmogorov-Smirnov test for normality (since difference represents error between models results). H0 was stated as data obeys a normal distribution and H1, data does not obey a normal distribution with a significance level of 5% (Pino, 2014Pino, F. A. (2014). A questão da não normalidade: uma revisão. Revista de Economia Agrícola, 61(2), 17-33.).

Study areas

Guariroba’s APA - MS - Brazil

The Environmental Protected Area (Área de Proteção Ambiental - APA) of the Guariroba River was created in 1995 by the Campo Grande Municipal Decree nº 7.183/1995, it is delimited by Guariroba’s river watershed and has about 360 km² of area. Guariroba River is a tributary of Botas River and belongs to Paraná River Basin. Guariroba’s APA is located in Brazilian Midwest Region, in Mato Grosso do Sul State, in Campo Grande City. Its Köppen climatic classification is Aw, i.e. rainy summer and dry winter. Figure 4 shows annual distribution of medium temperature and accumulated rainfall between 1996 and 2012 (Campo Grande, 1995Campo Grande. Prefeitura. (1995, 22 de setembro). Decreto nº 7.183, de 21 de setembro de 1995. Diário Oficial do Município de Campo, Campo Grande., 2008Campo Grande. Secretaria Municipal do Meio Ambiente e Gestão Urbana – SEMADUR. (2008). Plano de manejo da área de proteção ambiental dos mananciais do Córrego Guariroba – APA do Guariroba. Campo Grande., 2013Campo Grande. Prefeitura. (2013). Plano municipal de saneamento básico de Campo Grande. Campo Grande.; Cavazzana et al., 2019Cavazzana, G. H., Lastoria, G., & Gabas, S. G. (2019). Surface-groundwater interaction in unconfined sedimentary aquifer system in the Brazil’s tropical wet region. Revista Brasileira de Recursos Hídricos, 24, e8. http://dx.doi.org/10.1590/2318-0331.241920180136.
http://dx.doi.org/10.1590/2318-0331.2419... ).

Sedimentary Caiuá Group, created during the Early Cretaceous (133-120 Million Years Ago - Ma) after consolidation of the Serra Geral Group, occurs in almost entire Guariroba’s APA. Caiuá Group was formed by deposition of wind sediments (sand sea) and posteriori accommodation due sedimentary establishment of Bauru Group in Late Cretaceous (80-60 Ma) (Figure 5). This geological history created what is called Bauru Aquifer System (Sistema Aquífero Bauru – SAB, composed by Caiuá and Bauru Groups), an extensive outcropped free aquifer that occupies about 37% of the Mato Grosso do Sul State and where this study area is located (Figure 6). Sedimentologic analysis of geological units of Caiuá Group shows the presence of fine sandstone, well selected sediment and depths from dozens up to 280 meters, while units of Serra Geral Group are composed by fractured basaltic shedding. Deeper rivers beds reach the formations of Serra Geral Group in Guariroba’s APA, since they are superposed (Fernandes & Coimbra, 1994Fernandes, L. A., & Coimbra, A. M. (1994). O Grupo Caiuá (Ks): revisão estratigráfica e contexto deposicional. Revista Brasileira de Geociencias, 24(3), 164-176. http://dx.doi.org/10.25249/0375-7536.1994164176.
http://dx.doi.org/10.25249/0375-7536.199... ; Batezelli, 2010Batezelli, A. (2010). Arcabouço tectono-estratigráfico e evolução das Bacias Caiuá e Bauru no Sudeste brasileiro. Revista Brasileira de Geociências, 40(2), 265-285. http://dx.doi.org/10.25249/0375-7536.2010402265285.
http://dx.doi.org/10.25249/0375-7536.201... ; Mato Grosso do Sul, 2010Mato Grosso do Sul. Governo do Estado. Secretaria de Estado de Meio Ambiente, do Planejamento, da Ciência e Tecnologia. (2010). Plano estadual de recursos hídricos de Mato Grosso do Sul. Campo Grande: Editora UEMS.; Uechi et al., 2016Uechi, D. A., Gabas, S. G., & Lastoria, G. (2016). Analysis of the heavy metals in the Bauru Aquifer System in Mato Grosso do Sul. Engenharia Sanitaria e Ambiental, 22(1), 155-167. http://dx.doi.org/10.1590/s1413-41522016142430.
http://dx.doi.org/10.1590/s1413-41522016... ; Campo Grande, 2008Campo Grande. Secretaria Municipal do Meio Ambiente e Gestão Urbana – SEMADUR. (2008). Plano de manejo da área de proteção ambiental dos mananciais do Córrego Guariroba – APA do Guariroba. Campo Grande.; Cavazzana et al., 2019Cavazzana, G. H., Lastoria, G., & Gabas, S. G. (2019). Surface-groundwater interaction in unconfined sedimentary aquifer system in the Brazil’s tropical wet region. Revista Brasileira de Recursos Hídricos, 24, e8. http://dx.doi.org/10.1590/2318-0331.241920180136.
http://dx.doi.org/10.1590/2318-0331.2419... ).

Beneath Serra Geral Group, Corumbá Group is identified as a confined sedimentary set of geological unities composed by carbonate and silicate sediments deposited between 600 and 490 Ma, this geological group is usually known as Guarani Aquifer and reaches depths up to 1,800 meters in some regions and its boundaries extend for four South-American countries: Brazil, Argentina, Uruguay and Paraguay (Boggiani et al., 1993Boggiani, P. C., Fairchild, T. R., & Coimbra, A. M. (1993). O grupo Corumbá (Neoproterozóico-Cambriano) na região central da Serra da Bodoquena (Faixa Paraguai), Mato Grosso do Sul. Revista Brasileira de Geociencias, 23(3), 301-305. http://dx.doi.org/10.25249/0375-7536.1993233301305.
http://dx.doi.org/10.25249/0375-7536.199... ; Ribeiro, 2008Ribeiro, W. C. (2008). Aqüífero Gauarani: gestão compartilhada e soberania. Estudos Avançados, 22(64), 227-238. http://dx.doi.org/10.1590/S0103-40142008000300014.
http://dx.doi.org/10.1590/S0103-40142008... ).

Guariroba’s APA morphologic surface is composed by river plains separated by soft hills and the altitude varies between 665 and 455 m (Instituto Nacional de Pesquisas Espaciais, 2008bInstituto Nacional de Pesquisas Espaciais – INPE. (2008b). TOPODATA 20S555ZN. Retrieved in 2019, December 6, from http://www.dsr.inpe.br/topodata/index.php
http://www.dsr.inpe.br/topodata/index.ph... ; Campo Grande, 2008Campo Grande. Secretaria Municipal do Meio Ambiente e Gestão Urbana – SEMADUR. (2008). Plano de manejo da área de proteção ambiental dos mananciais do Córrego Guariroba – APA do Guariroba. Campo Grande.). Red, Red-Yellow Latosols and Neosol Quartzarenic present in higher areas, and Fluvial and Hydromorphic Neosols in fluvial valleys characterize the soil type of the watershed. Soil occupation activity is mostly cattle breeding and natural vegetation is predominantly of Cerrado biome type, shrubby vegetal specimens distributed in denser and sparser adjacent areas. Guariroba’s river flow is dammed downstream in a superficial reservoir where administration of Campo Grande City produces potable water to sustain a diversity of urban, industrial, rural and others anthropic activities (Campo Grande, 2008Campo Grande. Secretaria Municipal do Meio Ambiente e Gestão Urbana – SEMADUR. (2008). Plano de manejo da área de proteção ambiental dos mananciais do Córrego Guariroba – APA do Guariroba. Campo Grande.; Oliveira et al., 2017Oliveira, A. K. M., Fernandes, V., Pirajá, R. V., & Silva, M. H. S. (2017). Multitemporal evaluation of landscape of the Area of Environmental Protection (APA) of the water catchment Guariroba Stream, Campo Grande, Mato Grosso do Sul, by means of satellite images. Ra’eGa, 42, 8-20.; Capoane, 2019Capoane, V. (2019). Susceptibility to erosion in the Corrego Guariroba watershed, Campo Grande, Mato Grosso do Sul. Acta Brasiliensis, 3(2), 49-55. http://dx.doi.org/10.22571/2526-4338180.
http://dx.doi.org/10.22571/2526-4338180... ; Cavazzana et al., 2019Cavazzana, G. H., Lastoria, G., & Gabas, S. G. (2019). Surface-groundwater interaction in unconfined sedimentary aquifer system in the Brazil’s tropical wet region. Revista Brasileira de Recursos Hídricos, 24, e8. http://dx.doi.org/10.1590/2318-0331.241920180136.
http://dx.doi.org/10.1590/2318-0331.2419... ).

The Polygonal Domain is located in medium-lower part of the APA and it is delimited by 8 observation wells, an area of 181.88 km². Piezometric values were calculated by difference between measured water level depth and altitude of the point. Measurement expedition was accomplished in the same day and its values are an extract of data provided by Cavazzana et al. (2019)Cavazzana, G. H., Lastoria, G., & Gabas, S. G. (2019). Surface-groundwater interaction in unconfined sedimentary aquifer system in the Brazil’s tropical wet region. Revista Brasileira de Recursos Hídricos, 24, e8. http://dx.doi.org/10.1590/2318-0331.241920180136.
http://dx.doi.org/10.1590/2318-0331.2419... , who monitored water level of wells in Guariroba’s APA for two consecutive years (2016-2017). Statics levels used were collected in 02/12/2015 - rainy season.

Juazeiro do Norte - CE - Brazil

Juazeiro do Norte City is located in Brazilian Northeast Region, in Ceará State, has a population of 274,207 hab, 248,832 km² of area and belongs to metropolitan region of the Cariri. Köppen climatic classification is Aw - summer is the rainy season and winter is dry. Climatic properties of the region are widely affected by the annual displacement of Intertropical Convergence Zone (ITCZ) (Lima et al., 2017Lima, M. T. V., Meireles, A. C. M., Oliveira, C. W., & Nascimento, M. T. B. (2017). Koppen-Geiger and Thornthwaite climatic classification for the metropolitan region of the Cariri, Ceará. Revista Geama, 3(3), 136-143.; Lima & Ribeiro, 2012Lima, G. G., & Ribeiro, S. C. (2012). Geomorfologia e paisagem do município de Juazeiro do Norte/CE: relações entre a natureza semi-árida e os impactos antrópicos. Revista Geonorte, 2(4), 520-530.; Fundação Cearense de Meteorologia e Recursos Hídricos, 2012Fundação Cearense de Meteorologia e Recursos Hídricos – FUNCEME. (2012). Levantamento de reconhecimento de média intensidade dos solos: mesorregião do Sul Cearense (280 p.). Fortaleza: FUNCEME.; Instituto Brasileiro de Geografia e Estatística, 2018Instituto Brasileiro de Geografia e Estatística – IBGE. (2018). Área territorial brasileira. Rio de Janeiro: IBGE., 2019Instituto Brasileiro de Geografia e Estatística – IBGE. Diretoria de Pesquisas, Coordenação de População e Indicadores Sociais. (2019). Estimativas da população residente com data de referência 1o de julho de 2019. Rio de Janeiro: IBGE.).

Polygonal Domain is located in Cariri Valley, pediplan of the Araripe sedimentary basin (an interior basin), and its geological configuration is characterized by the following sequence of sedimentary geological units: Abaiara, Missão Velha and Brejo Santo formations. They were created by deposition of materials from erosion of higher crystalline exposed regions of Borborema Province during succession of rift processes in Late Jurassic and Early Cretaceous (Assine, 1992Assine, M. L. (1992). Análise estratigráfica da Bacia do Araripe, Nordeste do Brasil. Revista Brasileira de Geociencias, 22(3), 289-300. http://dx.doi.org/10.25249/0375-7536.1992289300.
http://dx.doi.org/10.25249/0375-7536.199... , 2007Assine, M. L. (2007). Bacia do Araripe. Boletim de Geociências da Petrobras, 15(2), 371-389.; Maia & Bezerra, 2014Maia, R. P., & Bezerra, F. H. R. (2014). Conditioning structural of relief in Northeast Brazilian. Mercator, 13(1), 127-141. http://dx.doi.org/10.4215/RM2014.1301.0010.
http://dx.doi.org/10.4215/RM2014.1301.00... ; Sales & Peulvast, 2007Sales, V. C., & Peulvast, J. P. (2007). Evolução morfoestrutural do relevo da margem continental do Estado do Ceará, Nordeste do Brasil. Caminhos de Geografia, 7(20), 1-21.; Lima & Ribeiro, 2012Lima, G. G., & Ribeiro, S. C. (2012). Geomorfologia e paisagem do município de Juazeiro do Norte/CE: relações entre a natureza semi-árida e os impactos antrópicos. Revista Geonorte, 2(4), 520-530.; Scherer et al., 2014Scherer, C. M. S., Sá, E. F., & Córdoba, V. C., Sousa, D. C., Aquino, M. M., & Canelas Cardoso, F. M. (2014). Tectono-stratigraphic evolution of the Upper Jurassic-Neocomian rift succession, Araripe Basin, Northeast Brazil. Journal of South American Earth Sciences, 49, 106-122. http://dx.doi.org/10.1016/j.jsames.2013.10.007.
http://dx.doi.org/10.1016/j.jsames.2013.... ; Camacho & Sousa, 2017Camacho, C. R., & Sousa, F. R. F. R. O. (2017). The structural framework of the Araripe sedimentary basin, Borborema Province, based on aeromagnetic data. Revista do Instituto de Geociências, 17(3), 140-161.; Melo et al.., 2018Melo, B. G. V., Carvalho, I. S., & Bertolino, L. C. (2018). Sedimentary enviroment of Brejo Santo Formation, Late Jurasisc of Araripe Basin, Brazil. Anuário do Instituto de Geociências, 41(1), 265-282. http://dx.doi.org/10.11137/2018_1_265_282.
http://dx.doi.org/10.11137/2018_1_265_28... ). Sediments of upper layer, Abaiara formation, are composed by fine up to medium granulometry sandstones with the thickness of 400m in some areas (Figure 7) (Scherer et al., 2014Scherer, C. M. S., Sá, E. F., & Córdoba, V. C., Sousa, D. C., Aquino, M. M., & Canelas Cardoso, F. M. (2014). Tectono-stratigraphic evolution of the Upper Jurassic-Neocomian rift succession, Araripe Basin, Northeast Brazil. Journal of South American Earth Sciences, 49, 106-122. http://dx.doi.org/10.1016/j.jsames.2013.10.007.
http://dx.doi.org/10.1016/j.jsames.2013.... ).

Watershed (Figure 8) is part of contribution area of Salgado River, an intermittent water body, tributary of Jaguaribe River. Natural vegetation is predominantly of the Caatinga Biome type, with notorious presence of deciduous and thorny vegetal specimens, adapted to low precipitation and high temperatures. Geomorphology is associated to soft hills and a shallow drainage system with the occurrence of inselbergs in east in contrast of the high declivities provided by Araripe geological formation (residual crystalline areas) in west portion. Soil occupation is characterized by intensive urban activity on the Juazeiro do Norte and adjacent cities, due to metropolis conurbation. Fluvial and Litolic Neosols associated to Latosols and Argisols characterize recent fluvial deposition and older deposition areas, respectively (Lima & Ribeiro, 2012Lima, G. G., & Ribeiro, S. C. (2012). Geomorfologia e paisagem do município de Juazeiro do Norte/CE: relações entre a natureza semi-árida e os impactos antrópicos. Revista Geonorte, 2(4), 520-530.; Fundação Cearense de Meteorologia e Recursos Hídricos, 2012Fundação Cearense de Meteorologia e Recursos Hídricos – FUNCEME. (2012). Levantamento de reconhecimento de média intensidade dos solos: mesorregião do Sul Cearense (280 p.). Fortaleza: FUNCEME.).

Polygonal Domain has 36.88 km² and it is delimited by 4 observation wells. Static hydraulic level values were obtained from Companhia de Gestão dos Recursos Hídricos (2014)Companhia de Gestão dos Recursos Hídricos – COGERH. (2014). Boletim de monitoramento dos poços com DataLogger no Cariri-CE (ano de 2013). Fortaleza: COGERH. and represent consisted values of annual monitored series. Piezometric values were calculated by difference between topographic altitude and static level depth.

Modelling parameters

Each observation well is associated to a Dirichlet boundary condition. Piezometric level is the dependent variable of the differential operator, which rules analytical values in domain, and has the same value of hydraulic head in a non-confined aquifer (Haitjema, 1995Haitjema, H. M. (1995). Analytic element modeling of groundwater flow. London: Academic Press.). Flow rate through boundary in a normal direction is associated to a Neumann boundary condition. Its value can be equal to zero in a lateral impermeable layer or another scalar value in a pumping well on the boundary. (Jazayeri & Werner, 2019Jazayeri, A., & Werner, A. D. (2019). Boundary condition nomenclature confusion in groundwater flow modeling. Ground Water, 57(5), 664-668. PMid:31001804. http://dx.doi.org/10.1111/gwat.12893.
http://dx.doi.org/10.1111/gwat.12893... ). Since analyzed polygons have no lateral confining layers, because domains are surrounded by sedimentary formations and wells data were collected in a stationary state, no Neuman boundary condition was implemented. However, a set of Dirichlet boundary conditions is sufficient to make the equation system well-determined, i.e. an equation system with only one solution or a well-posed boundary value problem (Jazayeri & Werner, 2019Jazayeri, A., & Werner, A. D. (2019). Boundary condition nomenclature confusion in groundwater flow modeling. Ground Water, 57(5), 664-668. PMid:31001804. http://dx.doi.org/10.1111/gwat.12893.
http://dx.doi.org/10.1111/gwat.12893... ).

Table 1 summarizes adopted values to parameters in MODFLOW simulation. Thickness of permeable layer was extracted from hypsometry and from geological information of each area. In Guariroba’s APA, Caiua Group (permeable) is located above Serra Grande Group (adopted as impermeable). Difference between highest and lowest topographic elevation values provides a gross approximation of modelled layer’s thickness. Since we are analyzing piezometric level in absolute altimetric value, the superposing of topographic surface only matters to an optional analysis, the interaction between groundwater and surface drainage system. Therefore, previous approximation is appropriated to compute piezometric levels inside domain. So, thickness can be estimated as 665-455 = 210 m. In Juazeiro do Norte City, the value adopted is equal to the thickness of the Abaiara formation, about 400m (Figure 6), the supposed geological formation where the water table occurs.

Hydraulic conductivity (kx, ky) were adopted equal to 6x10^-6 m/s (fine-grained sandstone) (Domenico & Schwartz, 1997Domenico, P. A., & Schwartz, W. (1997). Physical and chemical hydrogeology (2nd ed.). New York: John Wiley & Sons.). Specific storage (Ss), assuming incompressibility of water, can be adopted equal to 21% (Johnson, 1967Johnson, A. I. (1967). Specific yield-compilation of specific yields for various materials (Water Supply Paper, No. 1662-D). Washington: U.S. Government Printing Office.), since this value is lower than porosity (n) and corresponds the amount of that porosity which can be drained (or produced) in correspondent water volume. Porosity value (n) equal to 25% was adopted from Domenico & Schwartz (1997)Domenico, P. A., & Schwartz, W. (1997). Physical and chemical hydrogeology (2nd ed.). New York: John Wiley & Sons. to fine-grained sandstone.

RESULTS AND DISCUSSION

Figure 9 shows the comparison between error related to field values for MODFLOW and MFS results for a variable off-set distance (100m to 200,000m in steps of 100m). As affirmed by Golberg & Chen (1998)Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press. and Liu (2012)Liu, C. S. (2012). An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Engineering Analysis with Boundary Elements, 36(8), 1235-1245. http://dx.doi.org/10.1016/j.enganabound.2012.03.001.
http://dx.doi.org/10.1016/j.enganabound.... , the quality of numerical approximation is affected by the position of fictitious boundary. As off-set distance increases, error related to field data abruptly reduces to a minimum and then increases again until a stable value. After the 1500^th step (off-set distance = 150,000m) for both examples, computational routine reached ill-conditioned matrixes, however, it was able to finish the procedure until last step.

Results in Figure 9 contrasts with expected behavior of MFS applied in mathematical analysis and examples. Theoretically, numerical error describes an asymptotic curve tending to zero as the fictitious boundary moves away from domain limits. (Bogomolny, 1985Bogomolny, A. (1985). Fundamental solutions method for elliptic boundary value problems. SIAM Journal on Numerical Analysis, 22(4), 644-669. http://dx.doi.org/10.1137/0722040.
http://dx.doi.org/10.1137/0722040... ; Golberg & Chen, 1998Golberg, M. A., & Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems. In M. A. Golberg (Ed.), Boundary integral methods: numerical and mathematical aspects (pp. 103-176). Boston: Wit Press.; Liu, 2012Liu, C. S. (2012). An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation. Engineering Analysis with Boundary Elements, 36(8), 1235-1245. http://dx.doi.org/10.1016/j.enganabound.2012.03.001.
http://dx.doi.org/10.1016/j.enganabound.... ). Such behavior is observed in theoretical examples where analytical values are available for comparison purposes (Katsurada, 1994Katsurada, M. (1994). Charge simulation method using exterior mapping functions. Japan Journal of Industrial and Applied Mathematics, 11(1), 47-61. http://dx.doi.org/10.1007/BF03167213.
http://dx.doi.org/10.1007/BF03167213... ). Results plotted in Figure 9 suggest that field data, in both cases, do not represent an analytical solution of Laplace operator to that set of boundary conditions. This is acceptable considering that any modelling procedure is associated with uncertainties. Errors related to conceptual assumptions and measurement procedure can lead to such results (Middlemis et al., 2019Middlemis, H., Walker, G., Peeters, L., Richardson, S., Hayes, P., & Moore, C. (2019). Groundwater modelling uncertainty: implications for decision making: Summary of the Groundwater Modelling Uncertainty Workshop - Australasian Groundwater Conference. Australia: Flinders University.). Besides that, convergence of MFS was preserved, since all error tended to a stable value. These results show that the choice of an off-set distance in those study areas does not follow the general idea of how the MFS behaves in theoretical analysis. Wells 2PO-19 and 2PO-05 are located close to each other and this suggest that they are under the influence of the same set of errors, what reflects in similar shape of its respective curves in Figure 9.

It is important to notice that MODFLOW requests more input data and its numerical approach occurs propagating boundary and initial conditions cell-to-cell in a connected mesh through finite-difference form of continuity equation (Harbaugh, 2005Harbaugh, A. W. (2005). MODFLOW-2005, the U.S. Geological Survey modular ground-water model: the Ground-Water Flow Process: U.S. Geological Survey Techniques and Methods 6 – A16, variously p. Reston: U.S. Geological Survey.). The MFS uses a residual, integral formulation to transform integral into a summation at collocation points in a boundary of a well-posed problem (Brebbia & Dominguez, 1992Brebbia, C. A., & Dominguez, J. (1992). Boundary elements: an introductory course. Boston: WIT Press Computational Mechanics Publications.). While MODFLOW incorporates more details about hydrogeological characteristics of modeled areas, here adopted as literature values (hydraulic conductivity, porosity and specific storage), MFS only request information from boundary conditions and differential operator in the form of its fundamental solution, what is an advantage because it is less susceptible to input errors (associated to multiples parameters) and a disadvantage because it has less parameters to perform further calibration procedures. In addition, MFS results are not sensible to a refined mesh, which leads to less computational effort when we are interested in particular points inside the domain. MFS, however, in its basic formulation is limited to steady-state and homogeneous domain simulations.

In Juazeiro do Norte case, other meshless method, Kansa’s Method with Radial Basis Function, applied by Araújo (2016)Araújo, R. S. A. (2016). Metodologias baseadas no mínimo resíduo para a escolha do parâmetro de forma utilizado em RBF’s aplicadas no método assimétrico livre de malhas (MESHLESS) (Tese de doutorado). Centro de Tecnologia, Universidade Federal do Ceará, Fortaleza. computed best numerical value equal to 363m for the 1PO-67 well (5.60% of error related to field data), more inaccurate than MFS and MODFLOW models. Unlike MFS, Kansa’s Method requires a parameter “c” also know as shape factor to implement the numerical procedure. Despite the shape factor adds a source of error into analysis, this meshless method allows the positioning of collocation points inside the domain and does not require a fictitious boundary (Parand & Rad, 2013Parand, K., & Rad, J. A. (2013). Kansa method for the solution of a parabolic equation with an unknown spacewise-dependent coefficient subject to an extra measurement. Computer Physics Communications, 184(3), 582-595. http://dx.doi.org/10.1016/j.cpc.2012.10.012.
http://dx.doi.org/10.1016/j.cpc.2012.10.... ; Araújo, 2016Araújo, R. S. A. (2016). Metodologias baseadas no mínimo resíduo para a escolha do parâmetro de forma utilizado em RBF’s aplicadas no método assimétrico livre de malhas (MESHLESS) (Tese de doutorado). Centro de Tecnologia, Universidade Federal do Ceará, Fortaleza.). Soares Junior et al. (2012)Soares Junior, D., Godinho, L., Pereira, A., & Dors, C. (2012). Frequency domain analysis of acoustic wave propagation in heterogeneous media considering iterative coupling procedures between the method of fundamental solutions and Kansa’s method. International Journal for Numerical Methods in Engineering, 89(7), 914-938. http://dx.doi.org/10.1002/nme.3274.
http://dx.doi.org/10.1002/nme.3274... applied both methods, Kansa’s and MFS, to simulate acoustic wave propagation in heterogeneous media. The MFS was applied to model homogeneous part of the domain and Kansa’s method was applied to model heterogeneities. This same strategy can be further investigated to evaluate its application in groundwater modeling, mainly because it can propose alternatives to insert constraints such pumping wells and lakes in a meshless analysis, what the MFS cannot provide alone in its basic formulation.

Some extracted values from Figure 9 are displayed in Table 2. It shows that values computed by MFS and MODFLOW do not diverge more than 1% in relation to themselves. However, related to field values, MODFLOW reached best approximation. Once best off-set distance value was chosen (20,000m to Juazeiro do Norte case and 5,000m to Guariroba’s APA), Figure 10 was generated by interpolation of results for each study area for a grid (200x200m) of points inside domain. It shows that piezometric equipotentials have a similar shape for both methods and for both study areas.

The application of Mann-Whitney U test (Table 3) to compare numerical values in grid obtained by MODFLOW and MFS leads to rejection of H0 in Guariroba’s APA case, which means that ratings for both group of values suggest significant differences between them (U = 1.03+E07, n1 = n2 = 4480, p = 0.0237<0.05, two-tailed). In Juazeiro do Norte case, H0 was failed to reject, so both datasets are not significantly different (U = 419029, n1 = n2 = 906, p = 0.4394>0.05, two-tailed). Convergence of results provided by both methods indicates suitability between adopted hypothesis and numerical propagation of boundary values to points inside domain by finite-difference and fundamental solution methods to Juazeiro do Norte case. Although MFS result was proved significantly different from MODFLOW result in Guariroba’s APA case, the shape of piezometric levels allow the inference of similar groundwater flow behavior. This result does not point the inconsistency of one or other model but suggests the need for further investigation of representativeness of each model for this study area.

In Figure 10A, convergence of equipotential levels endorses the interpretation that two sub-basins share groundwater flow until they join themselves in same direction and contribute to surface flow in drainage system. Lastoria et al. (2017)Lastoria, G., Gabas, S. G., Cavazzana, G. H., Casadei, J., & Souza, T. A. (2017). Potencialidade dos recursos hídricos na bacia docórrego Guariroba, município de Campo Grande MS. In E. L. Cardoso (Ed.), Geologia ambiental: tecnologias para o desenvolvimento sustentável (Vol. 1). Ponta Grossa: Atena. applied the software SURFER 9.0 to draw piezometric levels to a set of observation wells and data collected in 01/12/2015 in Guariroba’s APA, rainy season. They also concluded that groundwater flow contributes to maintenance of superficial flow, what corroborates with the assumption of a non-confined aquifer for modelling purposes.

In Figure 10B, equipotentials show groundwater flowing obliquely to drainage system. Since field values are the result of an annual monitored series which was consisted, boundary values propagate this characteristic and indicates that general groundwater flow do not follow the superficial drainage during entire year. This is acceptable when considering intermittent behavior of drainage beds in semiarid regions and the existence of adjacent sedimentary geological formations which do not share necessarily same topographic boundaries with superficial watershed (Lacerda Filho et al., 2013Lacerda Filho, J. V., Valente, C. R., Lopes, R. C., Oliveira, I. W. B., Sachs, L. L. B., Silva, V. A., Lachaal, F., Mlayah, A., Bédir, M., Tarhouni, J., & Leduc, C. (2013). Groundwater flow modelling in semi-arid regions: the Zéramdine-Béni Hassen Miocene aquifer system (east-central Tunisia). In Colloque International Eau et Climat: Regards Croisés Nord/Sud. Fès: Faculté des Lettres et des Sciences Humaines Saïs.). Souza & Castro (2013)Souza, C. D., & Castro, M. A. H. (2013). Simulação do fluxo hídrico subterrâneo por estimativa de parâmetros usando cargas hidráulicas observadas: caso do Cariri cearense, Brasil. Recursos Hídricos, 34(1), 43-61. http://dx.doi.org/10.5894/rh34n1-4.
http://dx.doi.org/10.5894/rh34n1-4... investigated groundwater flow of an area that englobes the study area of Juazeiro do Norte. They simulated a more detailed model using daily data, MODFLOW model and MIGHA as a calibration process. They observed groundwater flow coherence, since flux pointed to lower areas and towards drainage system. In contrast, we used consisted values of annual monitored series, what leads to groundwater flow results representative for an entire year instead of a range of a particular day.

Although the MODFLOW was more accurate computing the piezometric levels of observation wells, it demonstrated inconsistencies working close to boundary. The MFS provided smoother values and, consequently, smother piezometric equipotentials. MFS also allows the determination of numerical values without a mesh, saving computational resources when there is a list of points of interest. Moreover, the approximation quality does not improve as domain is discretized for more points. Numerical value is the best value available by the method even if only one domain point is estimated to same set of boundary data.

Absolute error distribution (Figure 11) for both study cases reject H0 with Kolmogorov-Smirnov normality test (A: RMSE = 0.85, P-value = 1.21E-82 < 0.05, B: RMSE = 3.30, P-value = 6.61E-10 < 0.05), which means both error distributions are not described by a normal curve with significance of 5%. Error values associated to spatial points are not randomly observed. Further correlation analysis between datasets is encouraged. Figure 11 also shows that errors have lower values near boundary wells. This is expected because boundary data is a fixed parameter and both methods are forced to reach its values at those points.

CONCLUSIONS

We proposed to apply the MFS to model the groundwater flow of two study areas. It was required to choose the best off-set distance value as input parameter. So, minimizing the error related to field values of observation wells was the strategy adopted. Then we applied the same set of initial data into MODFLOW for comparison purposes.

Both models provided numerical solutions to estimate piezometric levels inside polygonal domains. However, limited field data induces the application of inaccurate methods for assessing the representativeness of results. Convergence between both results is a sign that both numerical methods provide a good approximation to the same set of initial values, but the lack of statistical convergence between them does not imply necessarily that one or other method is invalid. Only a specific study of representativeness with enough density of field data is able to assign this role to each model to each case. Besides that, and in general perspective, we expect that areas where Dupuit-Forchheimer assumption is representative and under a steady-state analysis it is possible to obtain significant numerical results with the MFS working as a modeling tool, since an off-set optimization method is implemented.

ACKNOWLEDGEMENTS

This research project was funded by a scholarship provided by Fundação Cearense de Apoio à Pesquisa (FUNCAP).

REFERENCES

Publication Dates

History