Abstract

To assess the effect of three grid cell properties (size, mean slope of the surface and distance between centre of grid and observation well) on groundwater models' performances, a tropical karst catchment characterized by monsoonal season in Rote Island, Indonesia was selected. Here, MODFLOW was used to develop models with five different spatial discretization schemes: 10 × 10 m, 20 × 20 m, 30 × 30 m, 40 × 40 m and 50 × 50 m. Using parameter estimation method, hydraulic conductivity and specific yield values over a selection of pilot points were estimated. The trends of the performances were calculated at each observation well in order to recommend the most appropriate location for observation well placement in terms of topographical characteristic. It is confirmed that the deterioration of model performance is mainly controlled by the increase of distance between well and centre of the cell, and the mean slope of the surface. Results reveal that model performance increases substantially for areas of low slope (<3%) and medium slope (3–10%) for a smaller grid cell size. Therefore, to improve model performance, it is recommended that the observations wells are placed in areas of low and medium slopes.

INTRODUCTION

Appropriate representation of spatial and temporal variability of hydrogeological parameters which determines model complexity is critical in model conceptualization (e.g., Li et al. 2014, 2015; Hailegeorgis & Alfredsen 2015; Pechlivanidis et al. 2017; Sun et al. 2016). In groundwater modelling the validity of both calibration and validation processes depends on how well the hydrogeological properties of the investigated area are conceptually, spatially and temporally represented in the model (Dehotin & Braud 2008). Discretization is one way to spatially and temporally assign physical properties of the real world to the model. A good discretization represented by optimal choice of size and number of cells would essentially improve the accuracy and efficiency of the model while reducing the complexity of the model through lesser time and cost consumption.

Spatial discretization comprises spatially dividing the model area into a finite number of cells to which model properties are characterized. In classic numerical methods, such as finite difference method, the nodes which are located in the centre of the cell represent the input components, such as precipitation and hydraulic conductivity. The node is also the place at which the numerical equation is solved to yield the output of the model, such as hydraulic head and discharge (Bundschuh & Arriaga 2010). It is well acknowledged that, centroid wise, the finer the grid size is, the better the node represents the output (e.g., Refsgaard 1997; Sciuto & Diekkrüger 2010; Wildemeersch et al. 2014). Figure 1 shows that the grid refinement progressively from 1-cell grid towards 16-cell grid results in a smaller distance between the node and observation well (dn). This obviously reduces the head disparity (Δhn) between simulated head (hs) and observed head (ho), which statistically improves model performance. However, on the other hand, due to the additional computational tasks as a result of increased grid numbers due to grid refinement, time consumption would exponentially increase. In some cases, much finer grid may only provide an insignificant model improvement at the huge expense of computational time.

Figure 1

Horizontal conceptualization of three examples of model grid representing observed head (ho), simulated head (hs), head disparity (Δhn) and distance between node and observation well (dn). The shadowed cells are the observed grid cells. Here, the denser grid represented by bigger cell number yields smaller distance (d3 < d2 < d1) and smaller head disparity (Δh3Δh2Δh1). In modelling, Δhn is assigned as statistics criterion.

Figure 1

Horizontal conceptualization of three examples of model grid representing observed head (ho), simulated head (hs), head disparity (Δhn) and distance between node and observation well (dn). The shadowed cells are the observed grid cells. Here, the denser grid represented by bigger cell number yields smaller distance (d3 < d2 < d1) and smaller head disparity (Δh3Δh2Δh1). In modelling, Δhn is assigned as statistics criterion.

A number of studies have explored the effect of grid cell size in model discretization using comparison of the performance of different models representing different cell sizes. Refsgaard (1997) developed four models (500, 1,000, 2,000 and 4,000 m finite difference grid) using MIKE SHE software for Karup catchment (440 km2) in Denmark. He analysed the statistics of discharge and hydraulic head of the models and proposed a guideline of using a maximum limit of 1,000 m grid model to simulate distributed hydrological models, as the coarser grid generated deteriorated model results. Li et al. (2008) studied the effect of topography representation as a result of coarsening of digital elevation model (DEM) resolution into two models (100 and 500 m grid) on evaporation rates in Yanqui Basin, China. They applied three methods (traditional, Taylor expansion and higher resolution techniques) on two resolution DEM inputs (500 × 500 and 100 × 100 m) to estimate subsurface evaporation with a MODFLOW-based groundwater model. They concluded that using finer resolution of surface topography gives better results in simulating evaporation calculation. Wildemeersch et al. (2014) compared the statistical performance of four models (250, 500, 750 and 1,000 m finite element grid) using a large synthetic catchment in Condroz, Belgium. Applying HydroGeoSphere, the effect of horizontal spatial discretization was investigated. They showed that although coarsening of grid size exponentially reduces computational time, the model accuracy with respect to hydraulic head and discharge decreased more than 50%. They also found that coarser models tend to overestimate discharge during high flow cycle as opposed to that during low flow periods.

All the studies highlighted the importance of grid refinement principle in enhancing model performances and concluded that, to some extent, poor representation of the catchment surface attributed to coarser grid results in deterioration of model performance. However, little is known about the influence of cell properties on model performance. In the aforementioned studies, the model performance was analysed as an average for the whole catchment without considering the position of each calibration point and its respective individual cell properties. Although grid size refinement effect on model performances was calculated, none of them has carefully explored the impact of other grid cell properties, i.e., mean slope and distance between cell node and point of observation. It is essential to accurately evaluate the calibration capacity in each cell of the respective observation well. The size of catchment is another important factor to consider in modelling the hydrogeological process. Dehotin & Braud (2008) stressed that dealing with a large catchment involves a high level of uncertainty in parameterization due to the demands of a coarser model resolution and creating a practically less-complicated model at the expense of model accuracy. Compared to the above-mentioned studies, the present study used a median catchment area (20.11 km2) which is considered small enough to better analyse the effect of grid refinement and adequately large enough to assess the impact of slope and centroidal distance of the observation wells on model performance.

Therefore, in this study we would like to confirm that by assessing the effect of grid cell properties on model calibration we are able to efficiently select the most appropriate cell size during the discretization step of model development; also, to identify the most suitable location for observation well placement in terms of topographical characteristic. Hence, a guideline to select observation wells based on the effect of grid cell properties (size and distance between centre of grid cell and observation well) and the representation of topography is provided in this study. The findings of this study can not only assist modellers to construct models in spatial discretization steps but also provide practical guidelines for decision-makers on the allocation of investment by determining where observation wells should be placed to help increase model performance. Therefore, the aim of the study is to assess the impact of grid discretization properties in the calibration and validation performances of a groundwater model. The results are evaluated using standard statistical criteria, i.e., mean error (MEh), root mean square error (RMSEh) and Nash–Sutcliffe coefficient (NSEh) using the observed and simulated heads. Instead of using a synthetic catchment, the study employs a natural real small-scale catchment of a karst spring for the advantage of reducing parameter uncertainty due to better representation of hydrogeological properties and thus increasing model accuracy.

STUDY AREA AND DATA COLLECTION

The study area is located on Rote island, which is in the south-eastern part of Indonesia, geographically between latitudes 10°46′42.17″S – 0°43′36.91″S and longitudes 123°3′14.84″E – 23°9′17.64″E. The recharge area of the Mamar spring, locally known as Oemau Spring, is located around 3 km from Ba'a, the capital of Rote Island (Figure 2) and has a catchment area of 20.11 km2. The spring is one of the major springs in the island in terms of cultural values, the quantity of water discharge and the diversity of purposes it feeds. It is also very important for the inhabitants of the capital as it is the only water source for the capital and neighbouring areas. On the other hand, the significant increase of population in the capital at the expense of diminishing grass, forest and plantation areas shows that the anthropogenic demographic factors may contribute to the decrease of groundwater supply to the spring.

Figure 2

Location of Oemau Spring and its recharge area.

Figure 2

Location of Oemau Spring and its recharge area.

Among the vegetation, savannah-type vegetation dominates the area while settlements scatter around the plantation and farm locations. Plantations are mainly found around the streams flowing from the upper reach to the spring location, whereas the farms are situated in the middle of the recharge area.

The geological bedrock of the recharge area of the Oemau Spring is categorized as karst as different types of carbonates rock dominate the area. Permeable Holocene coralline limestones (63.93%) mostly cover the upper terrains and areas downstream around the spring. Bobonaro formations (27.83%), which mainly consist of a mixture of carbonate rocks and scaly clay, are largely found in the middle of the recharge areas where rice fields are located. The remaining 8.25% is alluvium deposits mainly formed from physical rocks weathering and sediments scatter the areas around the streams.

The area experiences an average annual rainfall between 1,000 and 2,300 mm, having two distinct seasons: dry (April–November) and rainy (December–March). In the rainy season, monthly rainfall amount reaches to around 400 mm in February, while its intensity then usually decreases in subsequent months and may only reach 4 mm/month in the dry season (August). The dry season ends in October or the middle of November followed by a sharp rise of rainfall in late November or December. The humidity increases during the wet months (December–February) to around 92%, and subsequently drops in the dry season in November to as low as 75%.

Generally, Oemau Spring responds moderately to precipitation. Figure 3 shows the relationship between the hydraulic head fluctuation at an observation well located next to the spring's pool and rainfall recorded in Lekunik rain gauge station. It shows that the water table began to rise in early January 2012, about two months after the beginning of the rain event in early November 2011. Meanwhile, the temporal delay of the dry period exemplified a similar time lag as the head began to recede in mid-July 2011, approximately two months after the last rain event in mid-May 2011.

Figure 3

Time series of groundwater levels at Oemau Spring and cumulative rainfall at Lekunik Station.

Figure 3

Time series of groundwater levels at Oemau Spring and cumulative rainfall at Lekunik Station.

METHODOLOGY

Groundwater model

The study employed MODFLOW (McDonald & Harbaugh 1988) developed by USGS, operated in groundwater modelling system (GMS) user interface (BYU 2014) to simulate the hydrogeological processes in the catchment. MODFLOW is a finite difference-based numerical groundwater model which has been widely used in different groundwater problems and schematization including its application in karst areas (Quinn et al. 2006; Martínez-Santos & Andreu 2010; Panagopoulos 2012; Abusaada & Sauter 2013; Gallegos et al. 2013; Mayaud et al. 2014; Vallner & Porman 2016). The governing equation for groundwater flow is expressed as:  
formula
(1)
where Kx,y,z is the hydraulic conductivity along the x, y and z coordinates in ‘m/day’, W is the volumetric flux that represents sources or sinks of water (1/s), h is the hydraulic head (m), Sy is the specific yield, and x, y and z are coordinate directions, and t is time (d).

The models were simulated under steady-state and transient conditions. Under the steady-state and transient simulations, a no-flow boundary was used along the catchment boundary of the spring outlet assuming no flux flowing through the boundary. The catchment boundary was delineated using ArcGIS (ESRI 1999). Meanwhile, a specified head boundary was set for the small stream downstream of the spring using Time-Variant Specified-Head (CHD) package provided in GMS MODFLOW.

The conceptual model consists of one horizontal layer which assumes the continuum representation of unconfined carbonate aquifer in the study area. The model was constructed using Layer Property Flow package. This modest geologic stratum representation is offset by a detailed spatial rendering of hydraulic conductivity to enhance the fitting accuracy (Voss 2011) using geostatistical tool, parameter estimation (PEST). Prior information of hydraulic conductivity and specific yield values were used and determined using the Boulton method (Boulton 1963), which analysed pumping test data at three locations in the catchment area. The porosity was set as a constant value of 0.3 to represent the typical nature of the recharge area (Johnson 1967). A uniform assumption for porosity was used due to the fact that its small values (fraction between 0 and 1) insignificantly affect the distribution of heads in an aquifer. Moreover, its function of releasing water was already represented by hydraulic conductivity specifically parameterized in the model. The drains were simulated using Drain (DRN) package (Harbaugh et al. 2000) to model the surface drainage network of the recharge area (Moiwo & Tao 2014) – see Figure 4(b) for the drainage network. The surface drain conductance assigned at the streams was set between 2,700 and 8,125 m/d. The values are dependent on the dimension of the drain and K values which were derived from local pumping test analysis. The input into the model is the recharge mainly characterized by surface (i.e., land use type) and atmospheric stresses (i.e., precipitation and evapotranspiration). The land use types used in the model represent the four main land use types in the recharge area (i.e., savannah grass/bush, plantation, rain-fed farm and settlement) (Figure 4(c)). In order to quantify the recharge, the water balance method (Bras 1990) was used. Using SCS-CN method (Viessman & Lewis 2003), each land use was assigned the run-off value as the input for the recharge. The calibration and validation processes used the daily observed groundwater heads for a 16-month period collected from deven dug wells scattered in the recharge area.

Figure 4

Conceptual model and boundary conditions of the groundwater models showing: (a) model grid, (b) source and sinks, (c) recharge and (d) observation wells.

Figure 4

Conceptual model and boundary conditions of the groundwater models showing: (a) model grid, (b) source and sinks, (c) recharge and (d) observation wells.

The top elevations of the model were obtained from 2-m resolution DEM data provided by the Regional Spatial Planning Bureau on Rote Island. The conceptual model was constructed on the assumption that the aquifer is a single continuous and heterogeneous entity. The assumption was used mainly to reduce the complexity of the model, consider the limitation of local geological investigations, and signify the geological characteristics of the recharge area which is dominated by the carbonate formations. The model was then spatially discretized by a grid system comprising a finite difference mesh of cells. The grid cell of the model domain was homogenously discretized to five different sizes: 10 × 10 m, 20 × 20 m, 30 × 30 m, 40 × 40 m and 50 × 50 m, using finite difference method, to assess the impact of spatial discretization properties (cell size, mean slope and distance between observation well and centre of grid cell) (Figure 4(a)).

Model calibration and validation

Each of the five groundwater models was calibrated using two steps: (1) steady-state condition and (2) transient state condition, using recorded groundwater heads at seven observation wells. Steady-state calibration was performed to assist a smooth convergence of the model for the transient simulation. The main aim of the steady-state calibration in this study was to determine the distribution of horizontal hydraulic conductivity (K) across the recharge area in the model during the highest recharge as input to the model (Woessner & Anderson 1992). This was initially done by using parameter estimation method (Doherty 2001). In PEST, the pilot point technique (de Marsily 1984) was used to estimate the K values at 76 points evenly distributed across the model area by GMS. The distance among the points was equally set to 500 m. The estimated K values were then interpolated and distributed over the model using the inverse distance weighed method. The technique was used because it is able to generate a physically realistic K distribution in a heterogeneous system (Christensen & Doherty 2008). The maximal rainfall rate, 86.5 mm/day on 19 April 2011, was initially used as the recharge in the steady-state simulation to reflect the response of the highest groundwater head to maximal daily recharge. Transient simulations were performed to enhance the reliability of the simulation (Panagopoulos 2012) and aimed to calibrate K and Sy values for the duration of eight months from January to August 2011. In the validation step, the models were applied for another eight month period from September 2011 to April 2012. The simulation period was divided into 17 stress periods with a 10-day time step.

Model performance criteria

The performances of the developed models for the calibration and validation steps were evaluated using statistical goodness-of–fit measures: mean error (MEh), root mean squared error of head (RMSEh) and Nash–Sutcliffe coefficient (NSEh) as objective functions, respectively. The objective functions are described as follows:  
formula
(2)
 
formula
(3)
 
formula
(4)
where n is the number of observed wells, is the observed head (m) and is the simulated head (m).

The MEh was intended to measure the bias between the observed and simulated heads in the calibrated and validated models. The RMSEh was applied to indicate the global residual between observed and simulated heads. To evaluate how model results match with the observed ones in terms of assessing the relative degree of the residual discrepancy measured from the mean observed head, NSEh (Nash & Sutcliffe 1970) was used. A guideline from Moriasi et al. (2007) was used as criteria to evaluate NSEh values: very good (0.75 < NSEh ≤ 1.00), good (0.65 < NSEh ≤ 0.75), satisfactory (0.50 < NSEh ≤ 0.65) and unsatisfactory (NSEh ≤ 0.50).

RESULTS AND DISCUSSION

Grouping of slope cluster

The recharge area is dominated by a highly undulated topography with elevation ranging between 97.55 m and 339.55 m above sea level. Areas with mild slope are generally found in the middle of the catchment where rice farms are situated; and steeper slopes dominate the upper and lower areas of the catchment. Using ArcGIS the surface of the area studied was clustered into three slope groups (Figure 5), i.e., low, medium and high. The clustering is intended to categorize the area into distinct slope distribution and to explain the influence of each mean area slope to the performance of the model, in relation to different grid sizes.

Figure 5

Distribution of areas according to assigned slope cluster.

Figure 5

Distribution of areas according to assigned slope cluster.

The classification of the slope gradient is described in Table 1 and it is shown that the catchment is dominated by areas with low and medium clusters; 47.51% and 41.44% respectively, while high slope areas share only 11.04% of the total area. Table 2 shows where the seven observation wells are situated in the defined slope cluster. Four wells (OW5, OW6, OW13 and OW15), which are mainly located near the rice farm areas, fall in low slope cluster, while another two wells (OW2 and OW9) fall into medium slope cluster. Well OW17, situated in the hilly area upstream of the catchment, is categorized as being in high slope cluster.

Table 1

Slope cluster classification and its topology properties

Slope cluster Degree Area Mean slope Standard deviation 
Low slope ≤3 9.55 0.83 0.57 
Medium 3 < slope ≤10 8.33 4.64 2.13 
High slope >10 2.22 17.26 7.39 
Slope cluster Degree Area Mean slope Standard deviation 
Low slope ≤3 9.55 0.83 0.57 
Medium 3 < slope ≤10 8.33 4.64 2.13 
High slope >10 2.22 17.26 7.39 
Table 2

Slope cluster category of selected observation wells

  Observation well
 
 OW2 OW5 OW6 OW9 OW13 OW15 OW17 
Slope cluster Medium Low Low Medium Low Low High 
  Observation well
 
 OW2 OW5 OW6 OW9 OW13 OW15 OW17 
Slope cluster Medium Low Low Medium Low Low High 

Model calibration and validation

First steady-state simulation results of groundwater heads were compared with the recorded heads in the observation wells (Domenico & Schwartz 1998) and were analysed using RMSEh. Table 3 summarizes the results of the steady-state simulations for the five models. From the results having RMSEh values ranging from 0.33 to 0.56 m, it can be concluded that the simulated heads satisfactorily matched with the observed heads. Moreover, the groundwater flow direction of the model, which confirms the pre-set flow direction, indicates a well-calibrated model (Domenico & Schwartz 1998). The performances of the five models in transient calibration and validation steps are presented in Table 4.

Table 3

Performance of steady-state simulations for the five models

Model Grid cell size RMSEh (m) 
10 × 10 m 0.40 
20 × 20 m 0.37 
30 × 30 m 0.33 
40 × 40 m 0.41 
50 × 50 m 0.56 
Model Grid cell size RMSEh (m) 
10 × 10 m 0.40 
20 × 20 m 0.37 
30 × 30 m 0.33 
40 × 40 m 0.41 
50 × 50 m 0.56 
Table 4

Models' performances during transient calibration and validation periods

Model Grid cell size Calibration
 
Validation
 
RMSEh (m) NSEh RMSEh (m) NSEh 
10 × 10 m 0.24 0.95 0.30 0.64 
20 × 20 m 0.22 0.95 0.29 0.65 
30 × 30 m 0.28 0.93 0.36 0.51 
40 × 40 m 0.37 0.73 0.42 −0.21 
50 × 50 m 0.57 0.48 0.53 −0.60 
Model Grid cell size Calibration
 
Validation
 
RMSEh (m) NSEh RMSEh (m) NSEh 
10 × 10 m 0.24 0.95 0.30 0.64 
20 × 20 m 0.22 0.95 0.29 0.65 
30 × 30 m 0.28 0.93 0.36 0.51 
40 × 40 m 0.37 0.73 0.42 −0.21 
50 × 50 m 0.57 0.48 0.53 −0.60 

Effect of discretization of grid size

Overall, the results of the calibration of the five models show that the increase of grid cell size significantly reduces model performance as represented by RMSEh and NSEh (Figure 6(a) and 6(b)). The effect of coarsening in weakening model performance becomes more visible for models bigger than 30 × 30 m grid, characterized by an exponential increase of RMSEh. As expected, the simulation time logarithmically reduced as an effect of grid coarsening (Table 5 and Figure 6(c)). Calculated RMSEh results for varying cell sizes at each observation well are presented in Figure 7. For all the wells, a significant increase in model performances (lower RMSEh) with regard to grid cell refinements was observed. The most notable RMSEh increases, ranging from 568.3 to 1,769.1% (Table 6), occurred for the wells OW6, OW9, OW13 and OW15, which are situated in low/medium slope cluster. On the other hand, well OW17 in the high slope cluster, experiences a much lower increase of RMSEh (269.9%). Interestingly, two observation wells (OW2 and OW5), although located in medium and low slope areas, showed a moderately low rise in RMSEh values (36.2% and 39.3%). This comparatively low increase might be due to the reason that they are located adjacent to the spring, where a communal pool is situated downstream. The almost steady pool surface water elevation would maintain a much smaller difference of water fluctuation during the simulation.

Table 5

Properties of five models, calibration results' performances and processing time

Model Grid cell size Number of active cell Average RMSEh (m) Average NSEh Processing time (minutes) 
10 × 10 m 201,139 0.24 0.95 4,137.13 
20 × 20 m 50,316 0.22 0.95 1,716.63 
30 × 30 m 22,412 0.28 0.93 729.08 
40 × 40 m 12,663 0.37 0.73 427.78 
50 × 50 m 8,061 0.57 0.48 135.25 
Model Grid cell size Number of active cell Average RMSEh (m) Average NSEh Processing time (minutes) 
10 × 10 m 201,139 0.24 0.95 4,137.13 
20 × 20 m 50,316 0.22 0.95 1,716.63 
30 × 30 m 22,412 0.28 0.93 729.08 
40 × 40 m 12,663 0.37 0.73 427.78 
50 × 50 m 8,061 0.57 0.48 135.25 

The simulations were performed on a desktop computer Intel Core i7-4770 CPU @ 3.40 GHz using 16 GB RAM.

Table 6

Simulation results of the five models at the seven observation wells with respect to three grid cell properties

Grid properties Increase of RMSEh (%) with grid cell coarseninga
 
OW 2 OW 5 OW 6 OW 9 OW 13 OW 15 OW 17 
Cell size 36.2% 39.3% 1,769.1% 1,670.3% 568.3% 1,760.0% 269.9% 
Slope 25.6% 40.0% 1,721.4% 1,473.2% 2,041.4% 1,252.6% 255.4% 
Distance 34.9% 29.3% 685.4% 636.6% 2,228.9% 164.3% 198.1% 
Grid properties Increase of RMSEh (%) with grid cell coarseninga
 
OW 2 OW 5 OW 6 OW 9 OW 13 OW 15 OW 17 
Cell size 36.2% 39.3% 1,769.1% 1,670.3% 568.3% 1,760.0% 269.9% 
Slope 25.6% 40.0% 1,721.4% 1,473.2% 2,041.4% 1,252.6% 255.4% 
Distance 34.9% 29.3% 685.4% 636.6% 2,228.9% 164.3% 198.1% 

aThe increase of RMSEh denotes a decrease in model performance.

Figure 6

Effects of grid cell coarsening on model performance and computational time: (a) RMSEh vs. cell size, (b) NSEh vs. cell size and (c) computational time vs. cell size.

Figure 6

Effects of grid cell coarsening on model performance and computational time: (a) RMSEh vs. cell size, (b) NSEh vs. cell size and (c) computational time vs. cell size.

Figure 7

RMSEh results for different cell size discretization.

Figure 7

RMSEh results for different cell size discretization.

Effect of surface slope of the grid

Figure 8 presents the model performance at each observation well with respect to mean slope of the surface area. In general, the increase of slope in each well as a result of modification of grid cell size results in deterioration of model performance represented by increased RMSEh values. The performance decrease due to increased slope differs variably, with the most notable decrease appearing at wells OW6, OW9, OW13 and OW15 (between around 1,250 and 2,040%), which are located in low and medium slope areas (Table 6). In contrast, OW17 representing high slope area responds moderately by around a 250% decrease to the increased slope. Meanwhile, a performance change is barely seen in OW2 and OW5, with only around 25% and 40% respectively, conceivably owing to their positions which are adjacent to almost steady surface water elevation of a pool downstream of the spring.

Figure 8

Mean slope of the area in cell vs. RMSEh.

Figure 8

Mean slope of the area in cell vs. RMSEh.

Therefore, it can be concluded that, generally, in an area of smooth topography the possibility of significantly enhancing model performance using grid refinement method is greater than in high slope areas. The rationale behind this is that grid size coarsening resulting in higher mean slope gradients consequently generates smaller modelled recharge volumes. This in the end would decrease simulated head. The relationship between slope augmentation and simulated head is shown in Figure 9, which confirms that the higher slopes contribute to lowered simulated head (hs) and increased MEh values. This is in agreement with the study of Moiwo et al. (2010) and Chaplot (2014), who suggested that increased slope contributes to higher runoff in contrast to smaller discharge.

Figure 9

The effect of slope augmentation on MEh.

Figure 9

The effect of slope augmentation on MEh.

The effect of grid cell slope in determining the model performance accordingly shows that the selection of observation wells' locations, with regard to their topographical gradient, to be included into a model is critical because it significantly contributes to the effectiveness of the model to achieve a better statistical performance while maintaining less computation time to do so.

Effect of distance between centre of grid cell and observation well

In general, the grid coarsening results in widening of the length between observation well and the centre of the cell, where the well is positioned in the model (Figure 10). Similar to the extent of RMSEh reduction as a result of grid cell expansion in wells OW6, OW9, OW13 and OW15, the drops in model performances due to increased distance implicated by cell size coarsening at those wells are up to approximately 23 times poorer (Figure 11 and Table 6). As such, due to the expansion of grid cell size, the increase of distance between well and the centre of the cell has a parallel effect on the model performance. This confirms the concept explained in the beginning (Figure 1), that the more distant an observation point is positioned due to grid cell enlargement, the poorer the model performance.

Figure 10

Relationships between cell size and the distance between centre of grid cell and observation well.

Figure 10

Relationships between cell size and the distance between centre of grid cell and observation well.

Figure 11

RMSEh results vs. distance between centre of grid cell and observation well.

Figure 11

RMSEh results vs. distance between centre of grid cell and observation well.

CONCLUSIONS

Overall, it is concluded from this study that the effect of three grid properties (size, mean slope of the surface and distance between centre of grid and observation well) in model discretization is obvious, and thus the accuracy of models depends on how well the spatial resolution of grid cells is predefined. As expected, the model simulation results using smaller cell size demonstrate a significant improvement to those of bigger cell sizes; however, the objective of this study was to quantify the magnitudes of improvements. The study confirms that the increase of distance between well and centre of the cell as a result of the coarsening of grid cells results in deterioration of model performance. The magnitude of model performance reduction as a result of increased distance increases up to almost 23 times. Also, it is found that, generally, model performance increases most significantly for areas of low slope (≤3%) and medium slope (3–10%). In these two areas, model performance greatly increases from approximately seven to 19 times for smaller grid sizes. However, in high slope areas (>10%), model performance increases only a maximum of four times. Another important finding is that the model performance is controlled by the change in grid cell surface geometry represented by the mean slope of the respective cell. Steeper slope would generate smaller flux into the cell thus it decreases simulated head, which results in decreased model performance represented by higher values of RMSEh. This is confirmed in the study; in low and medium slope areas the model performs much better (between around 14 and 21 times) compared to that in high slope area (four times) for smaller grid sizes. Thus, the better representation of model surface area would improve simulation of hydrological processes (i.e., such as run-off and recharge) related with topographic slope. This henceforth provides an improved insight and rationale of not only efficient grid cell selection but also the preference for selecting or placing observation well(s). The most optimal model performance was obtained in the low and medium slope areas given the cell size and slope are reduced. Hence, it is recommended to place observation wells within low and medium slope areas to help improve model performance.

ACKNOWLEDGEMENT

The authors are grateful to the Australian government who provided the first author with a scholarship to undertake a PhD program at Swinburne University of Technology, Australia through Australia Awards Scholarship (AAS) program.

REFERENCES

REFERENCES
Abusaada
,
M.
&
Sauter
,
M.
2013
Studying the flow dynamics of a karst aquifer system with an equivalent porous medium model
.
Groundwater
51
(
4
),
641
650
.
Boulton
,
N. S.
1963
Analysis of data from non-equilibrium pumping tests allowing for delayed yield from storage
.
Proceedings of the Institution of Civil Engineers, London, Great Britain
26
(
3
),
469
482
.
Bras
,
R. L.
1990
Hydrology: An Introduction to Hydrologic Science
.
Addison-Wesley
,
Reading, MA
,
USA
.
Bundschuh
,
J.
&
Arriaga
,
M. C. S.
2010
Introduction to the Numerical Modeling of Groundwater and Geothermal Systems: Fundamentals of Mass, Energy, and Solute Transport in Poroelastic Rocks
.
CRC Press
,
Boca Raton, FL
,
USA
.
BYU
2014
Groundwater Modeling System (GMS) Version 9.2.9
.
Environmental Modeling Research Laboratory, Brigham Young University
,
Provo, UT
,
USA
.
de Marsily
,
G.
1984
Spatial variability of properties in porous media: a stochastic approach
. In:
Fundamentals of Transport Phenomena in Porous Media
(
Bear
,
J.
&
Corapcioglu
,
M. Y.
, eds).
NATO ASI Series E
,
vol. 82
,
Martinus Nijhof Publishers
,
Boston, MA
,
USA
, pp.
719
769
.
Doherty
,
J.
2001
PEST Groundwater Data Utilities
.
Watermark Numerical Computing
,
Brisbane
,
Australia
.
Domenico
,
P. A.
&
Schwartz
,
F. W.
1998
Physical and Chemical Hydrogeology
,
2nd edn
.
Wiley
,
New York
,
USA
.
ESRI
1999
ArcView GIS Version 3.2
.
Environmental Systems Research Institute Inc.
,
Redlands, CA
,
USA
.
Gallegos
,
J. J.
,
Hu
,
B. X.
&
Davis
,
H.
2013
Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP
.
Hydrogeology Journal
21
(
8
),
1749
1760
.
Harbaugh
,
A. W.
,
Banta
,
E. R.
,
Hill
,
M. C.
&
McDonald
,
M. G.
2000
MODFLOW-2000, the U.S. Geological Survey Modular Ground-water Model: User Guide to Modularization Concepts and the Ground-water Flow Process
.
USGS Open-File Rep 00-92, 1–127
.
Johnson
,
A. I.
1967
Specific yield – Compilation of Specific Yields for Various Materials
.
USGS Water Supply Paper 1662-D, 1–74
.
Li
,
H. T.
,
Kinzelbach
,
W.
,
Brunner
,
P.
,
Li
,
W. P.
&
Dong
,
X. G.
2008
Topography representation methods for improving evaporation simulation in groundwater modeling
.
Journal of Hydrology
356
(
1–2
),
199
208
.
Li
,
H.
,
Xu
,
C.-Y.
&
Beldring
,
S.
2015
How much can we gain with increasing model complexity with the same model concepts
?
Journal of Hydrology
527
,
858
871
.
DOI:10.1016/j.jhydrol.2015.05.044.
Martínez-Santos
,
P.
&
Andreu
,
J. M.
2010
Lumped and distributed approaches to model natural recharge in semiarid karst aquifers
.
Journal of Hydrology
388
(
3
),
389
398
.
McDonald
,
M. G.
&
Harbaugh
,
A. W.
1988
A Modular Three-dimensional Finite-difference Ground-water Flow Model
,
USGS, Techniques of Water-Resources Investigations 06-A1, 1–258
.
Moiwo
,
J. P.
,
Lu
,
W.
,
Zhao
,
Y.
&
Yang
,
Y.
2010
Impact of land use on distributed hydrological processes in the semi-arid wetland ecosystem of Western Jilin
.
Hydrological Processes
24
(
4
),
492
503
.
Moriasi
,
D.
,
Arnold
,
J.
,
van Liew
,
M.
,
Bingner
,
R.
,
Harmel
,
R.
&
Veith
,
T.
2007
Model evaluation guidelines for systematic quantification of accuracy in watershed simulations
.
Transaction of the ASABE
50
,
885
900
.
Nash
,
J. E.
&
Sutcliffe
,
J. V.
1970
River flow forecasting through conceptual models, part 1: a discussion of principles
.
Journal of Hydrology
10
(
3
),
282
290
.
Pechlivanidis
,
I. G.
,
McIntyre
,
N.
&
Wheater
,
H. S.
2017
The significance of spatial variability of rainfall on simulated runoff: an evaluation based on the Upper Lee catchment, UK
.
Hydrology Research
48
(
4
),
1118
1130
.
DOI:10.2166/nh.2016.038.
Quinn
,
J. J.
,
Tomasko
,
D.
&
Kuiper
,
J. A.
2006
Modeling complex flow in a karst aquifer
.
Sedimentary Geology
184
,
343
351
.
Sun
,
J.
,
Li
,
Y. P.
,
Huang
,
G. H.
&
Wang
,
C. X.
2016
Analysis of interactive effects of DEM resolution and basin subdivision level on runoff simulation in Kaidu River Basin, China
.
Hydrology Research
DOI:10.2166/nh.2016.332.
Vallner
,
L.
&
Porman
,
A.
2016
Groundwater flow and transport model of the Estonian Artesian Basin and its hydrological developments
.
Hydrology Research
7
(
4
),
814
834
.
DOI:10.2166/nh.2016.104.
Viessman
,
W.
&
Lewis
,
G. L.
2003
Introduction to Hydrology
.
Pearson Education
,
Upper Saddle River, NJ
,
USA
.
Wildemeersch
,
S.
,
Goderniaux
,
P.
,
Orban
,
P.
,
Brouyère
,
S.
&
Dassargues
,
A.
2014
Assessing the effects of spatial discretization on large-scale flow model performance and prediction uncertainty
.
Journal of Hydrology
510
,
10
25
.
Woessner
,
W. W.
&
Anderson
,
M. P.
1992
Selecting calibration values and formulating calibration targets for ground-water flow simulations
. In:
Proceedings of the NWWA Conference on Solving Groundwater Problems with Models
,
Columbus, OH
,
USA
, pp.
199
212
.