Meshless method using fundamental solution applied to computational simulation of groundwater flow of real aquifers: study case (Guariroba’s APA and Juazeiro do Norte)

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.


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, 2003). 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., 1984;McDonald & Harbaugh, 2003;Trefry & Muffels, 2007).
In late 1970s, Mathon & Johnston (1977) proposed the computational implementation of the Method of Fundamental Solutions and, in 1990s, Golberg & Chen (1998) outlined alternatives to apply the MFS to inhomogeneous and time-dependent problems, which was, until that moment, its major barrier . 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., 2019;Oh et al., 2019;Liu & Sarler, 2013;Liu, 2012;Barrero-Gil, 2012;Chen & Wang, 2010;Young et al., 2008;Alves & Chen, 2005;Golberg et al., 1999;Golberg & Chen, 1998;Bogomolny, 1985). Although the MFS was developed in 1960s , 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, 2016).
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, 1995;Wang & Zheng, 2016). 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, 2005). 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.

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, 1995). 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, 2019). 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, 1992).
where ( ) R x is the residual function applied in a collocation point x ∈ Γ, Γ is the boundary of a domain

3/15
Since approximated solution ( ) 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), Barrero-Gil (2012) and Liu & Sarler (2013).
We obtain Equation 6 writing Equations 3 and 4 in a matrixial notation, where f and g are the blocks of boundary conditions 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, 1998, Liu, 2012. Some strategies to circumvent ill-conditioning were considered in Chen et al. (2006), Young et al. (2008), Chen et al. (2016) and Hematiyan et al. (2016).
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, 1998;Wang & Zheng, 2016). 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, 2005;Golberg et al., 1999;Young et al., 2008). 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, 2019). 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, 1998;Liu, 2012), 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., 2019). 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), Companhia de Gestão dos Recursos Hídricos (2014) 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., 2012). 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, 2016). 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, 2010;Pino, 2014;Perme & Manevski, 2019).
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, 2014). Rodrigues Neto et al.

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 2 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, 1995Grande, , 2008Grande, , 2013Cavazzana et al., 2019).
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., 1993;Ribeiro, 2008).
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, 2008b; Campo Grande, 2008). 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, 2008;Oliveira et al., 2017;Capoane, 2019;Cavazzana et al., 2019).
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 2 . Piezometric values were calculated by difference between measured water level depth and altitude of the point. 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, 1992(Assine, , 2007Maia & Bezerra, 2014;Sales & Peulvast, 2007;Lima & Ribeiro, 2012;Scherer et al., 2014;Camacho & Sousa, 2017;Melo et al.., 2018). 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., 2014).
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, 2012;Fundação Cearense de Meteorologia e Recursos Hídricos, 2012).
Polygonal Domain has 36.88 km 2 and it is delimited by 4 observation wells. Static hydraulic level values were obtained from Companhia de Gestão dos Recursos Hídricos (2014) 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, 1995). 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, 2019). 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, 2019). 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, (Guariroba's APA and Juazeiro do Norte)  Hydraulic conductivity (kx, ky) were adopted equal to 6x10 -6 m/s (fine-grained sandstone) (Domenico & Schwartz, 1997). Specific storage (Ss), assuming incompressibility of water, can be adopted equal to 21% (Johnson, 1967), 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) to fine-grained sandstone. 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) and Liu (2012), 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, 1985;Golberg & Chen, 1998;Liu, 2012). Such behavior is observed in theoretical examples where analytical values are available for comparison purposes (Katsurada, 1994). 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., 2019). 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.

RESULTS AND DISCUSSION
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, 2005). 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, 1992). 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) 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, 2013;Araújo, 2016). Soares Junior et al. (2012) 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) 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  (1) Sd = Standard deviation, (2) n = Number of points, (3) U = Sum of ranks, (4) Z = z-Value.
watershed (Lacerda Filho et al., 2013). Souza & Castro (2013) 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.