DOI QR코드

DOI QR Code

Groundwater Analysis for Jun Stream Basin in Donghae City using MODFLOW

MODFLOW를 이용한 동해시 전천유역의 지하수분석

  • Kim, Sam Eun (Department of Civil Engineering, Inha University) ;
  • Kim, Hung Soo (Department of Civil Engineering, Inha University)
  • 김삼은 (인하대학교 사회인프라공학과) ;
  • 김형수 (인하대학교 사회인프라공학과)
  • Received : 2021.10.06
  • Accepted : 2021.11.04
  • Published : 2021.11.30

Abstract

The objective of this study is to evaluate the change of groundwater and its effect to hydraulic head or parameter sensitivity in order to predict the effect of a new groundwater well using MODFLOW. With the downstream area of the Jun stream basin in the Donghae city, two condition were simulated; (i) First simulation condition is to simulate natural steady-state and the groundwater heads for observation wells are calibrated. After the calibration, existing wells are simulated to find out the pure effect of the next simulation condition, (ii) Secondly, the design wells are planed and simulated for the solution of water shortage problem at Donghae city. As a result of the simulation, the maximum drawdown which is approximately 0.35m occurs in layer 1. This represents that there is no intrusion of salty water and the design plan is acceptable. Finally, the sensitivity analysis of each parameter in MODFLOW is accomplished. We found out that the sensitive parameters are initial head, anisotropic factor, hydraulic conductivity, and general head boundary

이 연구의 목적은 MODFLOW를 이용하여 신규 지하수 관정 추가로 인한 영향을 예측하기 위해 지하수의 변화와 지하수두 또는 매개변수 민감도에 미치는 영향을 평가하는 것이다. 이 연구에서는 동해시의 전천유역 하류지역을 대상으로 두 가지 조건을 모의하였다. 첫 관측된 지하수위값을 이용하여 자연상태 정류 흐름을 모의 할 수 있는 매개변수를 결정한 후, 이미 설치되어 있는 관정까지 포함하여 모형을 구축하였다. 구축된 모형을 이용하여 지하수 관정의 추가 설치에 따른 영향을 모의하였다. 두 번째는 동해시의 물 부족 문제를 해결하기 위한 방법으로 관정을 계획하고 모의하였다. 모의 결과, 최대 지하수위 저하는 첫 번째 대수층에서 약 0.35m로 나타났으며, 이는 해수의 침입은 발생하지 않는다는 것과 관정의 계획이 물 부족의 해결책으로 무리가 없다는 것을 나타낸 것이다. 마지막으로 MODFLOW에서 사용되는 매개변수의 민감도 분석을 실시하였다. 민감도 분석 결과에서 초기 수두, 비등방성 계수, 수리전도도 및 일반수두경계가 민감한 매개변수로 판명되었다.

Keywords

1. Introduction

Groundwater refers to water in a saturated zone of geologic stratum. Groundwater constitutes the largest available source of fresh water far greater than all the lakes, reservoirs, and streams. In order to utilize a groundwater as water resource, it is necessary to understand the properties of groundwater flow. The characteristic of groundwater flow is very complex and it is usually modeled as a simplified mathematical form for the practical applications. MODFLOW is one of the simplified mathematical groundwater flow model.

MODFLOW has been widely used for groundwater simulations and practical applications. After the birth of MODFLOW, many researchers have been applied the model to real groundwater basins.

Kazemi and Rahnama (2006) simulated groundwater flow of Rafsanjan aquifer using MODFLOW and they investigated the decrease on groundwater table effects on subside of land in this plain. They study result indicated that decrease of water level with subside have a liner relationship. Rezaei and Sargezi (2009) studied the effects of artificial recharge on the aquifer of Goharkooh Plain and specified the best location for the implementation of artificial recharge using MODFLOW model. The results showed that the aquifer response to artificial recharge was positive and the artificial recharge didn’t have a destructive effect on the aquifer. Panagopoulos (2012) confirmed that MODFLOW tool has more capability to simulate the groundwater resources under the variable hydrodynamic characters of an aquifers system. Needhidasan and Nallanathel (2013) described that MODFLOW has used the input to construct and solve the equations of groundwater flow in aquifer system.

The purpose of this study is for the applications of a three-dimensional ground water flow model to the aquifer system in Jun stream watershed in Donghae city. Though based on the existing data, this research aims to simulate a real groundwater basin using MODFLOW in order to predict the effect of new hydraulic stress on the basin, and to examine how much the deviations of parameters used in MODFLOW contribute to the change of hydraulic head at the basin

The specific objectives of this study include the simulation and calibration of a natural steady-state which has no well stresses, comparison of superiority of solution package in simulating complex groundwater basin, simulation study for the design wells provided for satisfying water shortage in Donghae city and, sensitivity analysis of the model parameters in the packages.

2. A Modular Three-Dimensional Finite-Difference Groundwater flow Model

Infancy of the study, Trescott (1975), Trescott and Larson (1976), and Trescott, Pinder, and Larson (1976) constructed the two- and three-dimensional finite-difference model. U.S. Geological Survey and others have studied the computer simulation of groundwater flow based on their previous works. The major purpose in designing a new groundwater flow model were to produce a program that could be readily modified, was simple to use and maintain, could be executed on a variety of computers with minimal changes, and was relatively efficient with respect to computer memory and execution time. In this study, majority of the explanation of MODFLOW is referenced to McDonald and Harbaugh (1988).

MODFLOW uses a modular structure wherein similar program functions are grouped together. New options can be added without the necessity of changing existing subroutines because of this structure. The program was originally written using FORTRAN 66 (McDonald and Harbaugh, 1984). It has subsequently been modified to use FORTRAN 77. This report documents the FORTRAN 77 version. The program is highly portable; it will run, without modification, on most computers. On some computers, minor modification may be necessary or desirable.

The major packages that are presently available include procedures to simulate the effects of wells, recharge, rivers, drains, evapotranspiration, and “general-head boundaries”. The solution algorithms available included two iteration techniques, the Strongly Implicit Procedure (SIP) and the Slice-Successive Overrelaxation method (SSOR).

The three-dimensional movement of ground water of constant density through porous earth material may be described by the partial-differential equation :

\(\begin{aligned}\frac{\partial}{\partial x}\left(K_{x x} \frac{\partial h}{\partial x}\right)+\frac{\partial}{\partial y}\left(K_{y y} \frac{\partial h}{\partial y}\right)+\frac{\partial}{\partial z}\left(K_{z z} \frac{\partial h}{\partial z}\right)-W=S_{s} \frac{\partial h}{\partial z}\end{aligned}\)

where Kxx, Kyy and Kzz are values of hydraulic conductivity along the x, y, and z coordinate axes, which are assumed to be parallel to the major axes of hydraulic conductivity (Lt-1); h is the potentiometric head (L); W is a volumetric flux per unit volume and represents sources and/or sinks of water (t-1); Ss is the specific storage of the porous material (L-1); and T is time (t).

Basic Package, Block-Centered Flow Package, River Package, Recharge Package, Well Package, Evapotranspiration Package, General Head Boundary Package are used in this study.

3. Applications

1) Study Area

Jun stream located at Donghae city in Kangwon province is a small size stream. The surface water basin area is 123.0km2 , the length of the stream is 18.42km, and the shape of the basin is a branching type. The surface water basin of Jun stream is located at 128°5¢ 332 of east longitude and 37°25¢ 242 of north latitude. The stream is generated from Samwhua Dong, flows to northeast side, joins with Shinhung stream, and flows to the East Sea in Donghae city. Precipitation of the site is 1,237mm. The slope of the stream bed is approximately 1/250 from the estuary to the point of the confluence and 1/50 at upstream.

The flow direction at upstream is approximately parallel to the direction of foliation of neighboring strata or strike of cracks. Tributaries are developed to the east-west direction, the upstream is a steep slope, and the longitudinal cross section shows its irregularity. Hydrogeologic units in this site are mainly formed as a granite, limestone, unconsolidated sediments, halfconsolidated sediments, crystalline metamorphic rock, metasedimentary rock, and early sedimentary rock.

As a result of the investigation of hydrogeology in this site, it was proved that there are two fault lines which are lying on the north-south direction in parallel and there are apparent differences in stratal distribution. Two fault lines divide the area into three sites which have the different hydraulic characteristics. Figure 1 shows the study area of Jun stream groundwater basin.

HKSJBV_2021_v23n4_342_f0001.png 이미지

Fig. 1. The study area of Jun stream groundwater basin.

2) Discretization and Boundary Conditions

The grid is composed of 85 rows by 104 columns and it has a 50m wide by 50m long. The model vertically divide the study area into five layers. Layer 1 represents a unconfined aquifer. It contains continental deposit, fluvial deposit and surface layers. Layer 2 represents marine deposite layer that contains fine-grid silty sand and clay under the continental deposit layer near the East Sea. It plays role of a confining layer to the Bukpyung pudding layer. Layer 3 represents the Bukpyung pudding layer which is confined by the upper layer. Layer 4 represents thick mudstone layer which is impermeable. It has confining effect on the layer 5. Layer 5 is confined Pungchon limestone layer. Drill logs shows discrete geologic structures in this site. Based on the elevation and bottom elevation data from field drill log in each layer, the bottom elevations in whole basin are interpolated by using the Kriging method. The closer to the east coast, the deeper bottom elevation of each layer is and the thicker each layer is. Upper layers have the irregular distribution of bottom elevations of layer while lower layers have simple one. In MODFLOW, a layer cannot simply vanish when a pinch out occurs by the real stratification and limitation of data. Thus, this study uses 0.1m as a minimal thickness of the layer.

Boundary conditions are for the simulation on the effect of external flow in interested area. But it is impossible to determine distinct delineation of a groundwater basin without abundant geologic survey. Therefore, the general head boundary is utilized as a boundary condition to complement those problems. The border of groundwater basin is delineated along 30m-contour line of ground surface, which has small effect of upland, and the resident and drill logs are included within the boundary. Groundwater is assumed to flow from inland to the East Sea, namely, from the west to the east. Therefore, the outer groundwater inflows from the west inland. The northern and southern boundaries are applied as no flow boundaries and the eastern and western boundaries are mainly general head boundaries. Even though the cells are located at the northern and southern boundaries, if they are assumed as orthogonal groundwater flow lines, they may be general head boundary conditions. Initial boundary heads in General Head Boundary Package (GHB) were substituted by an initial observed head for steady-state simulation. Constant head boundaries are utilized as a boundary condition of the northeast part of the site where the East Sea bounds.

3) Hydraulic Characteristics

Aquifer tests for determining hydraulic characteristics were performed from 10th of February to 18th of May in 1996 (KOWACO, 1996). MODFLOW needs all hydraulic conductivities at the center of each model grids. A Kriging method is adopted in order to solve the problem with the limited measurements of hydraulic conductivities. The hydraulic conductivity around northeast side is greater than that of southeast part of the basin.

In MODFLOW, the vertical flow between layers is controlled by the vertical conductance coefficients, Vcont. In the process of calibration, it is necessary to change hydraulic conductivities for matching the calculated and observed groundwater levels. It takes a time to convert Vcont according to changes of hydraulic conductivities. Because of the inconvenience, Fortran program was developed to convert Vcont automatically.

4) Others

Jun stream functions as a source of recharge to the aquifer and a recipient of discharge from the aquifer. The data used in River Package is referenced from “A Report of Channel Improvement Plan for Jun Stream” (1994, Donghae city). Since it is considered that the effect on the aquifer system would be small, the tributaries of Jun stream are excluded.

Recharge rate is estimated as 0.0005m/day from the report by KOWACO(1996). Evapotranspiration estimate is 572mm per year when the equivalent method of soil moisture is used with rainfall data at Samchuk station. This value is referenced from the report by KOWACO (1996).

5) Well

Stresses from well are not considered in the first simulation condition in which model is calibrated because there is no dewatering in natural state. Though there is a personal small well, it is ignored because its effect is very slight and there is no exact data. Second simulation condition is to add the existing wells to the first calibrated steady-state simulation to represent the natural state. This simulation is inevitable since it makes it possible to examine the pure effect of third simulation of the design wells on groundwater flows at the basin. Row and column numbers and the discharge of existing wells are presented in table 1. Third simulation condition is for the prediction of the availability of groundwater by using design wells. Its objective is to investigate the change of groundwater system by the design wells. The location of wells are assumed based on the following considerations for Donghae city. First, Avoid the distributions of the densely grouped wells. Second, pumping discharge rate should be fine. Third, where it is easy to secure the site of pumping facility. Fourth, Groundwater should not be contaminanted.

Table 1. Existing Wells

HKSJBV_2021_v23n4_342_t0001.png 이미지

4. Sensitivity Analysis

Sensitivity of parameters which are used in the study is checked in order to investigate the effect of tested parameters on the Jun stream groundwater basin. Initial criteria for the sensitivity analysis is determined from the simulation study that includes the design wells for the predicted demand. That is to say, the hydraulic heads of layers in the last simulated results are chosen as the initial criteria of the sensitivity analysis. RMSE between hydraulic heads from the simulation study and changed heads from the sensitivity analysis of parameters is employed as the index of deviation.

Shead represents initial hydraulic head of aquifer and is one of sensitive parameter. Upper thinner lines for each layer represent RMSE caused by change of Shead, and lower thicker lines represent deviations of RMSE and the difference between RMSE and deviation of Shead. That is to say, lower thicker lines show real sensitivity because Shead changes hydraulic head so directly that increase or decrease caused by deviation of Shead should be excluded. By the way, Shead is very sensitive parameter since it has direct effect.

TRPY is an anisotropy factor. It affects considerably to hydraulic head. In spite of sensitive characteristic of TRPY, MODFLOW assigns one constant anisotropy factor to a layer. Therefore, distributed input method is required for more exact simulation. HY is hydraulic conductivity. HY is multiplied by a factor in sensitivity analysis for preventing reasonable HY. It has a small effect on hydraulic head when deviation of hydraulic conductivity value is compared with the RMSE. However, hydraulic conductivity field data have so wide range that it needs caution to get the field data. BOT is the elevation of the aquifer bottom. The elevation of the layer 3 is most sensitive in the basin.

EVTR represents maximum evapotranspiration rate. It becomes more sensitive when it exceeds the standard model maximum evapotranspiration rate than when it doesn’t. EXDP is evapotranspiration extinction depth and SURF is evapotranspiration surface elevation. Both parameters are not sensitive. Though it is difficult to determine both parameters exactly, it has a small effect to hydraulic head.

GHB Head is the head on the boundary and it is sensitive parameter. The smallest change occurs in layer 5 since layer 5 has no general head boundary cells. RMSE of this analysis is linearly proportional to the deviation of GHB Head. GHB Cond is the hydraulic conductance of the interface between the aquifer cell and the boundary. The parameter is initially needs the by hydraulic conductivity multiplied by the thickness of assigned layer. Both parameters in General Head Boundary Package are sensitive. Therefore, a lot of field survey should be accomplished to evaluate those parameters.

RECH is the recharge flux. It is most sensitive in layer 1 where RECH affect directly because layer 1 is an unconfined aquifer. As RECH is increasing, the RMSE of hydraulic head in aquifers increase.

RIV Cond and Stage are parameters in RIV Package. RIV Cond is riverbed hydraulic conductance and Stage is the head in the river. They are insensitivity parameters. RMSEs caused by changes of the parameters are within 0.2m.

5. Calibration

Even though we simulate the packages in MODFLOW with the data sets obtained by field survey, the uncertainty of groundwater flow always exits due to the limitations of geologic distribution data, the heterogeneous characteristics of aquifer and so on. Therefore, the calibration process is inevitable for the reasonable simulation of groundwater basin. After the calibration, the changes of groundwater flow caused by the pumping of design wells can be simulated. This study assumes that hydraulic conductivity and boundary conditions may have the greatest uncertainty in the simulation according to the results of sensitivity analysis in the previous section. The calibrations for these two characteristics are mainly accomplished. The calibration for hydraulic conductivities is performed in a physically proper range based on the hydraulic conductivities interpolated with Kriging method. The reasonable result of calibration is evaluated from the comparisons of the calculated and observed groundwater heads and the groundwater budget in the basin.

The layer 2 and 4 are impermeable and the calibration is performed for the layers of 1 and 3. The observations of hydraulic head in layer 5 could not be found in the previous study (KOWACO, 1996) and it is not considered in this study for the calibration purpose. Tables 2 and 3 show the target head used for the calibration and the calculated head by calibration. Table 2 is for layer 1 and table 3 is for layer 3. There are several abnormal large residuals which represent the difference between the target and calculated heads in table 2 and 3 (table 2: MW-28, MW-29, JCW-6 and table 3: JCW-2). There may be two answers for the large residuals. First, the calibration is not enough to represent the characteristics of the basin. Second, the hydraulic characteristics of this abnormal point are so different that the additional field survey of aquifer may be needed. RMSE of the target and calculated heads at layer 1 is 0.20m with removing abnormal points and that at layer 3 is 0.35m. Approximately 1ft (30.48cm) is acceptable error of groundwater elevation in other researches (Restrepo at al., 1992).

Table 2. The target and calculated heads for layer 1

HKSJBV_2021_v23n4_342_t0002.png 이미지

Table 3. The target and calculated heads for layer 3

HKSJBV_2021_v23n4_342_t0003.png 이미지

Table 4 shows water budget related to each package used in the calibration process. Inflow of ‘Constant Head’ is zero. It means that there is no inflow of salty water from the East Sea. The water budget is mainly dependent on ‘Head Dependent Bounds’. That is to say, most of groundwater flow occurs around the border of groundwater basin. The steady-state calibration is accomplished relatively well, since ‘Percent Discrepancy’ is -0.01, near to zero.

Table 4. Water budget for each package used in the calibration process

HKSJBV_2021_v23n4_342_t0004.png 이미지

This study uses the heterogeneous hydraulic conductivities obtained by Kriging method and thus the hydraulic conductivities show high variations in the horizontal and vertical directions. Consequently, it is difficult to solve numerically. When the hydraulic properties show high variability, we can use SOR Package in MODFLOW for the simulation purposes. Figure 2 is the graph that shows convergence process of solution according to solution packages(SIP and SOR in MODFLOW) in the calibration. While the maximum head change for each iteration suddenly decreases and converges approximately in the 7th iterations when SOR Package is used (Figure 4.5), it severely vibrates when SIP Package is used. It can be seen that SIP solution Package is not proper to solve complicated multi-layer aquifer as shown in figure 4.5. The SOR solution technique works fairly well for this problem which has most of the heterogeneity in the vertical (Anderson, 1993).

HKSJBV_2021_v23n4_342_f0002.png 이미지

Fig. 2. Comparison of maximum head change for each iteration.

6. Simulation of Design Well

After the calibration for the natural state, the simulation for existing wells should be done before the simulation of design wells for the safe demands of water for domestic and industrial purposes of 20,000m3/day in Donghae city. The reason why the simulation of existing wells should be done before the simulation of design wells is to examine the pure effect of design wells on the Jun stream basin. Drawdowns caused by existing wells are calculated based on the calibrated results. The maximum drawdowns at layer 1, 3, and 5 are respectively about 0.07, 0.6, and 1.0m. Approximately the maximum 1m-drawdown at layer 5 occurs from the contours of drawdown compared to calibrated natural state. When we simulate for existing wells, there is no dried cell and so groundwater is not drying up for the operation of existing wells. Drawdowns near the east coast area is so insignificant that salty water doesn’t intrude into inland.

In order to choose the locations of well field sets, the most important things are the yield capacity and quality of groundwater among the commented conditions as afore mentioned. The data for yield capacity and quality of groundwater is referenced from the report by KOWACO (1996).

The existing morphology of groundwater is varying according to the positions. Yielding form of groundwater is divided into three types. The table 5 shows the distribution of aquifers and underground hydraulic characteristics. The specific pumpage in the second row represents that of JCW-3 observation well as a result of pumping test. It is smaller than any other specific pumpage in the basin. Water quality samples from seven drilled observation wells in Jun stream basin shows that water from JCW-3 is not proper.

Table 5. Aquifer distribution and hydraulic characteristic in Jun stream basin

HKSJBV_2021_v23n4_342_t0005.png 이미지

Since fault belts longitudinally exist from north to south around Jun stream bridge and the site around JCW-3 is excluded for choosing the design well sites. The fields of an design wells are determined in four places in the directions of the east and west sides of the basin. Each design well is arbitrarily distributed with the similar distance in the selected zones. The wells are sited near the Jun stream in order to put the wells around main stream. The fields of design wells are planned to be located like figure 3.

HKSJBV_2021_v23n4_342_f0003.png 이미지

Fig. 3. Field locations of design wells.​​​​​​​

The distributed discharges in each layer are calculated based on the transmissivities in assigned cells and the results are shown in the last column of table 6.

Table 6. Characteristics of design wells​​​​​​​

HKSJBV_2021_v23n4_342_t0006.png 이미지

The result of simulation for plan to construct design well is presented in figures 4. Contour lines of head drawdown are presented for the assigned layers of each well.

HKSJBV_2021_v23n4_342_f0004.png 이미지

Fig. 4. The drawdown of layers caused by the designed wells.​​​​​​​

Synthetically judged, there are no severe changes of drawdown contours caused by adding new well (i.e., design wells). Drawdown of southern part of the basin is larger than that of northern part. Drawdown at layer 3 is much smaller than drawdown at layer 1 because layer 3 is confined aquifer. Also, the drawdown affects the hydraulic head at the southern part of the basin comparatively far away. Drawdown occurs only at IW-11 point which is unique in layer 5. So, it can’t affect hydraulic head to the east coast. Because there is no dried cell when the design wells are simulated, it is considered that groundwater does not run dry for the simulation. The guess that groundwater does not dry up is supported by the water budget to be generated by MODFLOW, in which the inflow of groundwater at the basin equals the outflow. Additionally, the groundwater budget explains that there is no influx from the east coast caused by the design wells and we can know that there is no inflows in the constant head cells from MODFLOW outputs.

Finally, the plan to install new wells (design wells) near to the downstream of Jun stream can satisfy for predicted water demand, 20,000m3/day, and the groundwater system will not meet serious problem.

7. Conclusion

Even though the water can be obtained from the groundwater basin by existing wells, it is insufficient to satisfy the water demands for the domestic and industrial purposes in Donghae city. It was investigated that more water of 20,000m3/day is needed for satisfying water demand in Donghae city. Therefore, this study used MODFLOW to provide sufficient water demands to Donghae city from the groundwater basin through the processes of simulation, calibration, and sensitivity analysis using the design wells. We can obtain the water of 24,000m3/day from 12-design wells and it is sufficient to satisfy the water shortage problem of Donghae city. Also the study provided that the groundwater basin has no effect even though the basin losses the water as much as 24,000m3/day.

First step was to simulate natural steady-state which has no well pumping stresses. Hydraulic characteristics of aquifers were obtained at the points of drill logs for pumping tests. Using these characteristics, the calibration was performed and the simulated heads were compared with observed heads. Second step was to simulate the effect of existing wells and then the effects of design wells were investigated for the hydraulic heads. Third step was to simulate the design added wells for the increased water demand in the future. The required amount of water was assumed as a 24,000m3/day. New wells were distributed in the vicinity of a main stream of Jun stream by considering the condition of a basin. Finally, the sensitivity analysis for the parameters were performed based on the results from the third step.

Based on the results of this study, the following conclusion could be obtained. SOR Solution Package is better than SIP Solution Package to solve a complex groundwater basin for Jun stream. The maximum drawdown is approximately 0.35m after simulating the design wells to predict the water demand of 24,000m3/day. Therefore, the new stress will not cause a serious problem in the aquifer system. South part of the basin has larger drawdown and drawdown affects comparatively far away. Sensitivity parameters affecting hydraulic head in the aquifer at Jun stream basin are the parameters related to initial head, anisotropic factor, hydraulic conductivity, and general head boundary.

References

  1. Kazemi F, Rahnama MB (2006), "Modeling of groundwater in Rafsanjan aquifer and investigation of groundwater discharge effects", Thesis. Shahid Bahonar Kerman Iran University.
  2. KOWACO (1996), "Preliminary Report for Feasibility Study for the Link Development of Surface and Groundwater in the Northern Part of East Coast Area"
  3. McDonald, M. G. and Harbaugh, A. W. (1984), "A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model", U.S. Geological Survey.
  4. McDonald, M. G. and Harbaugh, A. W. (1988), "A Modular Three-Dimensional Finite-Difference Ground-Water Flow Model", U.S. Geological Survey.
  5. Needhidasan, S., & Nallanathel, M. (2013). "Application of Visual MODFLOW and GIS in Groundwater Modeling", International Journal of Civil Engineering (IJCE), 2(3), 41-50.
  6. Panagopoulos, G. (2012). "Application of MODFLOW for simulating groundwater flow in the Trifilia karst aquifer", Environ Earth Science, 67, 1877-1889. https://doi.org/10.1007/s12665-012-1630-2.
  7. Restrepo, J. I., Bevier, C. and Butler, D. (1992), "A Modular Three-dimensional Finite-Difference Groundwater Flow Model of the Surficial Aquifer System", South Florida Water Management District.
  8. Rezaei, M. Sargazi, A. (2009) "The Effects of Running Artificial Recharge on the Aquifer of Goharkooh Plain", Journal of Earth Science, 19, 99-106.
  9. Trescott, P. C. (1975), "Documentation of Finite-Difference Model for Simulation of Three-Dimensional GroundWater Flow", U.S. Geological Survey.
  10. Trescott, P. C. and Larson, S. P. (1976), "Documentation of Finite-Difference Model of Three-Dimensional Ground-Water Flow", Supplement to Open-File Report 75-438, U.S. Geological Survey.
  11. Trescott, P. C., Pinder, G. F., and Larson, S. P. (1976), "Finite-Difference Model for Aquifer Simulation in Two Dimensions with results of Numerical Experiments: U.S. Geological Survey.