Integrated hydrological modeling of surface-groundwater interactions : the case of Denpasar-Tabanan Basin in the Southern Bali Island

INTEGRATED HYDROLOGICAL MODELING OF SURFACE - GROUNDWATER INTERACTIONS

The case of Denpasar-Tabanan Basin in the Southern Bali Island

Abenet Tadesse Teketel March, 2017

SUPERVISORS:

Dr. Maciek W. Lubczynski

Dr. Zoltan Vekerdy

INTEGRATED HYDROLOGICAL MODELING OF SURFACE - GROUNDWATER INTERACTIONS

The case of Denpasar-Tabanan Basin in the Southern Bali Island

Abenet Tadesse Teketel

Enschede, The Netherlands, March, 2017

Thesis submitted to the Faculty of Geo-Information Science and Earth Observation of the University of Twente in partial fulfilment of the requirements for the degree of Master of Science in Geo-Information Science and Earth Observation.

Specialization: Water Resource and Environment Management SUPERVISORS:

Dr. Maciek W. Lubczynski Dr. Zoltan Vekerdy

THESIS ASSESSMENT BOARD:

Dr. Ir. Christiaan van der Tol

Dr. W. van Verseveld (External Examiner, Deltares, The Netherlands)

DISCLAIMER

This document describes work undertaken as part of a programme of study at the Faculty of Geo-Information Science and Earth Observation of the University of Twente. All views and opinions expressed therein remain the sole responsibility of the author, and do not necessarily represent those of the Faculty.

ABSTRACT

The Denpasar-Tabanan (D-T) Basin (2270 km

) located in the southern part of Bali Island, Indonesia, is a densely populated area, also typical destination of international tourism. The basin is composed of unconsolidated volcanic materials storing groundwater in an unconfined aquifer. That aquifer provides the main water supply of the Bali Island although it often undergoes crises of water shortage. Therefore, it is crucial to understand its complex dynamic of surface-groundwater (SW-GW) interactions and to evaluate its groundwater resources to improve water management.

The dynamics of SW-GW- interaction was numerically simulated using three-dimensional steady-state and transient models. For that purpose, MODFLOW-NWT with stream flow routing (SFR2) and unsaturated zone flow (UZF1) packages were used. All data, including time series of rainfall, stream discharge, and potential evapotranspiration for the four-year period from 1

January 2009 to 31

December 2012, were simulated on a daily basis.

The steady-state model groundwater inflow consisted of: R

(95%) and q

(5%). The groundwater outflow consisted of: q

(47.8%), ET

(23.4%), Exf

(22.6%), and q

(6.1%) of total groundwater outflow. In the transient model simulation, groundwater inflow consisted of R

(75.4%), q

and ∆S

were 3.3% and 21.3%

of total groundwater inflow respectively. Regarding transient model groundwater outflow, q

(30.4%), E

f

(29.4%), ET

(14.1%) followed by ∆S

(18.9%) and q

(7.2%) of total groundwater outflow respectively. It was observed that in the years analysed the net recharge was largely positive and that streams largely gain groundwater, which reflects good groundwater potential of the D-T basin.

The calibrated transient model showed large spatiotemporal variability of groundwater fluxes. The R

ranged from 7.64 (January) to 2.30 mmday

(August) with the mean 3.36 mmday

, R

from 5.84 (January) to 0.26 mmday

(August) with the mean 1.34 mmday

, ET

from 1.05 (February) to 0.65 mmday

(July), E

f

from 1.48 (March) to 1.37 mmday

(December) and q

from 1.48 (January) to 1.37 mmday

(August). The temporal variability of fluxes was mainly due to seasonal variability of driving forces changing from dry to wet season. The large spatial variability of groundwater fluxes was primarily due to large variety of land covers and large spatial variability of rainfall.

Key Words: Surface-groundwater interactions, Bali, Volcanic aquifer, Water balance, MODFLOW-NWT

ACKNOWLEDGEMENTS

The preparation and completion of this thesis work would not have been possible without the involvement and contribution of kind persons around me. First, I thank the almighty God for the gift of life, good health and for His guidance to this far. Many thanks to the Netherlands Fellowship Program (NFP) for the generous financial support throughout my study period.

I would like to express my deep gratitude to my first supervisor, Dr. Maciek W. Lubczynski, for his supervision, patience, words of encouragement, constructive ideas and invaluable remarks throughout this research work. My deep gratitude also goes to my second supervisor Dr. Zoltan Vekerdy, for his supervision, useful comments, remarks, and words of encouragement, throughout the thesis study. It would be so difficult without their help. Sincere thanks to Dr. Tom H.M. Rientjes, Dr. ir. Suhyb Salama, and ir. Arno M.

van Lieshout for their valuable comments during the process of this research starting from thesis proposal.

You all patiently directed my wild imaginations to a worthwhile research item and ensured that I stayed within the scope of the study by frequently correcting my work. Sincere thanks to Miss Novi Rahmawati (Ph.D. candidate), who selflessly ensured I had all the necessary data, and for guiding me through the data that are written in Balinese language. My thanks are extended to government of Indonesia especially to Indonesian Agency for Meteorology, Climatology, and Geophysics Locally called DPPU; Biro of Geospatial Information, “BIG”; and Ministry of Energy and Mineral Resources of Indonesia, “KLPU” for sharing the meteorological and hydrological data for all station that are available in Bali Island.

Many thanks to the administrative and teaching staff of the ITC-WRM Program. It is here where I learned

most of my theoretical foundations in water resources management. To my fellow classmates, thank you

for the quality time spent together as one big ‘family’; in a way, you helped me in the completion of this

thesis. Last but not the least, my heartfelt gratitude to my family for their unconditional love, patience and

support to the seemingly endless journey. Special thanks to my dad and mom for this encouragement and

trust that forever keeps me strong. I am because you were. To my elder brother, my greatest source of

inspiration.

TABLE OF CONTENTS

1. INTRODUCTION ... 1

1.1. Background ...1

1.2. Problem statement ...2

1.3. Research setting ...3

1.3.1. Research objectives ... 3

1.3.2. Research question ... 3

1.3.3. Novelty of the study ... 3

1.3.4. Research hypothesis... 3

1.3.5. Assumptions ... 3

2. RESEARCH METHOD AND MATERIALS ... 4

2.1. Study area ...4

2.1.1. Location ... 4

2.1.2. Monitoring network ... 5

2.1.3. Climate ... 5

2.1.4. Topography and land cover ... 6

2.1.5. Hydrology ... 7

2.1.6. Hydrogeology ... 8

2.1.7. Previous studies in the area ... 8

2.2. Data processing...9

2.2.1. Watershed boundary ... 10

2.2.2. Precipitation ... 10

2.2.3. Potential evapotranspiration ... 12

2.2.4. Interception and infiltration rate ... 13

2.2.5. Stream discharge ... 14

2.2.6. Hydraulic properties ... 15

2.2.7. Head observation ... 15

2.2.8. Groundwater abstraction ... 15

2.3. Modeling flow chart ... 16

2.4. Conceptual model... 16

2.4.1. Defining Hydrostratigraphic units ... 17

2.4.2. Defining the flow system ... 17

2.4.3. Defining preliminary water balance ... 18

2.4.4. Defining the boundaries of the model ... 18

2.5. Numerical model ... 18

2.5.1. Software selection ... 18

2.5.2. Aquifer geometry and grid design ... 20

2.5.3. Driving forces ... 20

2.5.4. State variables ... 21

2.5.5. Parametric data ... 21

2.5.6. Boundary conditions ... 22

2.6. Model calibration ... 24

2.6.1. Steady-state model calibration ... 24

2.6.2. Warming-up period for transient model calibration ... 24

2.6.3. Transient model calibration ... 25

2.7. Error assessment and sensitivity analysis ... 25

2.8. D-T basin water balance ... 26

3. RESULTS AND DISCUSSION ... 28

3.1. Data processing calculation results... 28

3.1.1. Filling missed data for precipitation ... 28

3.1.2. Consistency of the precipitation records ... 28

3.1.3. Spatial data interpolation of rainfall ... 29

3.1.4. Interception and infiltration rate ... 31

3.1.5. Potential Evapotranspiration [PET] ... 32

3.1.6. Consistency of stream discharge ... 34

3.2. Steady-state model calibration ... 36

3.2.1. Calibrated head and error assessment ... 36

3.2.2. Calibrated stream discharges ... 38

3.2.3. Hydraulic conductivities ... 39

3.2.4. Water budget of the steady-state simulation using GHB conditions at the sea coast ... 39

3.2.5. Spatial variability of groundwater fluxes ... 42

3.2.6. Effects of changing GHB conductance upon lateral groundwater outflow to the ocean ... 43

3.2.7. Water budget of the steady-state simulation using CHD boundaries at the sea coast ... 43

3.2.8. Sensitivity analysis ... 45

3.3. Transient model calibration ... 46

3.3.1. Calibration heads and error assessment ... 47

3.3.2. Calibrated stream discharges ... 50

3.3.3. Hydraulic conductivities and specific yield ... 55

3.3.4. Water budget of the transient-state simulation ... 56

3.3.5. Temporal variability of groundwater fluxes ... 57

3.3.6. Spatial variability of groundwater fluxes ... 58

3.3.7. Yearly steady-state and transient variability of water fluxes ... 59

3.3.8. Sensitivity analysis ... 61

4. CONCLUSIONS AND RECOMMENDATIONS ... 63

4.1. Conclusions ... 63

4.2. Recommendations ... 64

Appendix I ... 69

Appendix II ... 69

Appendix III ... 70

Appendix IV ... 71

Appendix V ... 71

Appendix VI ... 74

Appendix VII ... 75

LIST OF FIGURES

Figure 1: Geological map and cross section across the study area (Modified after Purnomo & Pichler, 2015).

... 2

Figure 2: Location and elevation map of the D-T basin. ... 4

Figure 3: Monitoring stations of D-T basin ... 5

Figure 4: Mean monthly rainfall; maximum, minimum, and mean temperature at station Kuta of D-T basin. For the location of the station see Appendix V. ... 6

Figure 5: Land use map and percent area coverage of D-T basin. ... 7

Figure 6: Watershed boundaries, stream segments and gauging location of D-T basin. For the name of each station see Appendix V (B). ... 7

Figure 7: Geology of the study area. ... 8

Figure 8: Potentiometric surface of Bali Island (Modified after Nielsen & Widjaya, 1989). ... 9

Figure 9: Double mass curve for a precipitation data (after Gómez, 2007). ... 11

Figure 10: Methodology of flow chart. ... 16

Figure 11: Geological cross section across the study area (After Ministry of Energy and Mineral Resources of Indonesia)... 17

Figure 12: Proposed boundary conditions and locations in the D-T basin. ... 23

Figure 13: Schematic diagram of MODFLOW-NWT setup of D-T basin model. ... 26

Figure 14: Daily rainfall after filling missed data at station Kuta for the years from 2009 to 2013. For the location of station see Appendix V (A)... 28

Figure 15: Double mass curves of the precipitation gauges [units in mm]. The double mass curve in A & C shows consistency in the data but B & D shows inconsistency in the data. For the location of stations see Appendix V (A). ... 29

Figure 16: Sample significance test results of rainfall record on January 10, 2009. ... 30

Figure 17: Standard model variogram; distance in a unit of [m] and semi-variance in a unit of [m

]. ... 30

Figure 18: Kriged prediction and Kriging variance of D-T basin for long-term average rainfall from 01/01/2009 to 01/01/2012 [unit – mday

]. ... 31

Figure 19: Spatially variable interception (A) and infiltration rate (B) of D-T basin. ... 32

Figure 20: Temperature coefficient of determination for Sanglah and Kuta stations. ... 32

Figure 21: Spatially variable crop coefficient (A) and extinction depth (B) for D-T basin. ... 33

Figure 22: Average Rainfall (P), Infiltration rate (P

), Interception rate (I) and Potential evapotranspiration (PET) for four hydrological years from 2009 to 2012. ... 33

Figure 23: Sample double mass curve and frequency distribution for the stream gauge discharge data [Q- stream discharge in m

sec

]. For the location of station and log transform see Appendix V. ... 35

Figure 24: Relation between rainfall and stream discharge. The oval shape indicates uncertainty on the relation between measured stream flow and rainfall pattern. [RF – rainfall and Q – stream discharge]. For the location of stream gauging stations see Appendix V. ... 35

Figure 25: Relationship between simulated and observed heads in the D-T basin during steady-state IHM for the year from 2009 to 2012. ... 36

Figure 26: Potentiometric surface with location of heads, stream segments of the D-T basin during steady- state IHM. The stream segments with black lines indicates those that were not included during model calibration. Heads in m a.s.l. ... 37

Figure 27: Calibrated horizontal hydraulic conductivity (K

) distribution map of D-T basin after steady-state IHM [unit - mday

]. ... 39

Figure 28: Schematic representation volumetric water budget in case of steady-state IHM for the entire model of D-T Basin [ All units - mmyear

]. ... 41

Figure 29: Spatially variable ET

of D-T Basin for calibrated steady-state IHM [Unit – mday

]. ... 42

Figure 30: Spatially variable R

map for calibrated steady-state IHM in D-T Basin [Unit – mday

]. ... 43 Figure 31: The relationship between GHB conductance and head dependent boundary flow rate. The GHB conductance that masked by orange circle indicate the final value that was selected. In MODFLOW-NWT the GHB conductance is calculated based on polyline objects as in Section 2.4.6. ... 43 Figure 32: Sensitivity of model for horizontal hydraulic conductivity (A) & Vertical unsaturated zone hydraulic conductivity (B). ... 45 Figure 33: Sensitivity of model for UZF1 package parameter and driving forces: (C) extinction depth, (D) infiltration rate, (E) extinction water content, (F) potential evapotranspiration. ... 46 Figure 34: Sensitivity of model for (G) Brooks-Corey-Epsilon & (H) saturated water content. ... 46 Figure 35: Relationship between simulated and observed heads for the transient IHM of 11 observation points. ... 47 Figure 36: The potentiometric surface and stream segments of the D-T basin during transient model calibration at the last stress period, December 31, 2012. ... 48 Figure 37: Time series for the comparison of yearly observed and daily simulated heads for D-T basin. P – rainfall, H

– observed heads, and H

– simulated heads. ... 49 Figure 38: Relationship between observed and simulated discharge in the D-T basin for the transient-state model calibration of 13 stream gauge (2009 - 2012). For the location of gauges see Appendix V (A & B).

The oval shape in Figure 37 (J, K, and L) show that uncertainty in the measured rainfall and stream discharges. Because it is expected that at high rainfall records, stream discharge is higher and the opposite is true. Q – stream discharge and RF – rainfall. ... 54 Figure 39: Calibrated horizontal hydraulic conductivity (K

) distribution map of D-T basin after Transient- state IHM [unit - mday

]... 56 Figure 40: Temporal variability of groundwater fluxes in transient model calibration for gross recharge (R

), net recharge (Rn), surface leakage (E

f

), and groundwater evapotranspiration (ET

). ... 57 Figure 41: Temporal variability of rainfall, actual infiltration and PET ... 58 Figure 42: Spatially variable ET

map for calibrated transient IHM during dry (A) and wet (B) period in D- T Basin [Unit – mday

]. ... 59 Figure 43: Spatially variable Rg map for calibrated transient IHM during dry (A) and wet (B) period in D-T Basin. ... 59 Figure 44: Sensitivity of model for horizontal hydraulic conductivity [K

], maximum unsaturated zone vertical hydraulic conductivity [Kv

], extinction depth [EXTDP], extinction water content [EXTWC], Saturated water content [WC

], and Brooks-Corey-Epsilon [BC]. ... 61 Figure 45: Effects of changing unsaturated zone vertical hydraulic conductivity [Kvun] upon infiltration (A) and Exf

(B), effects of changing extinction depth [EXTDP] and GHB conductance upon ET

(C) and q

(D) respectively. ... 62

LIST OF TABLES

Table 1: Available data from the year 2009 to 2012 (T

- maximum temperature, T

- minimum temperature, RH – relative humidity, WS – wind speed, SS - sunshine duration, n.a - indicates that the data are not available within the study periods, EXTDP – Extinction depth, K

– saturated hydraulic

conductivity, and S

- specific yield). ... 10

Table 2: Interception loss rate that used in D-T Basin... 14

Table 3: Aquifer characteristics value of the basin (T - transmissivity, S

– specific yield, K

- horizontal hydraulic conductivity, SWL - surface water level). The location of bores can be seen in Appendix I. ... 15

Table 4: D-T extinction depth based on land cover. ... 21

Table 5: Observed and simulated head with calculated error assessment for 11 piezometers, H

– Observed head, H

– Simulated head [units – m]. ... 36

Table 6: Observed and simulated stream discharge with calculated error assessment for 16 gauges: Q

– observed stream discharge; and Q

– simulated stream discharge [unit - m

day

]. Stations that are highlighted by red colour indicate those station that are not used for model calibration and show the model response for those stations... 38

Table 7: Total water balance of D-T basin at steady-state IHM [mmday

]. ... 40

Table 8: water balance of land surface and unsaturated zone [mmday

]. ... 40

Table 9: Water balance of groundwater in steady-state IHM [mmday

]. ... 41

Table 10: Water balance of groundwater in steady-state condition [mmday

]... 44

Table 11: Observed, H

and simulated head, H

with calculated error assessment for 11 piezometers. [Units – m]. ... 48

Table 12: Observed and simulated stream discharge with GHB conductance of 0.1 m

day

per unit length calculated error assessment for 16 gauges in m

day

. Stations that are highlighted by red colour indicate those station that are not used for model calibration and shows the model performance for those stations. ... 55

Table 13: Final calibration output for model parameters and model variables in the D-T basin: EXTDP – extinction water content; EXTWC – extinction water content; THTS – saturated volumetric water content; THTI – initial volumetric water content; STRTOP – streambed top; STRTHICK – streambed thickness; SLOPE – stream slope; STRHC1 – streambed hydraulic conductivity; WIDTH1 – stream width; Kv

– maximum unsaturated zone vertical hydraulic conductivity; K

– horizontal hydraulic conductivity; S

– specific yield; and C – conductance. ... 56

Table 14: Long term average groundwater budget for entire model in transient-state IHM [mmday

] for the 2009-2012. IN – inflow to the aquifer system, OUT – outflow from the aquifer system, GW – groundwater. ... 57

Table 15: The yearly variability of driving forces and different groundwater balance components over the

three hydrological periods 1

January 2009 till 31

December 2012 MODFLOW-NWT simulation period

[All units in mm year

]. ... 60

LIST OF ABBREBATIONS

∆S Change in storage

BGI Biro of Geospatial Information

BMKG Indonesian Agency for Meteorology, Climatology, and Geophysics

CHD Time-Variant Specified-Head

DEM Digital Elevation Model

DPPU Ministry of Public Works

D-T Denpassar - Tabanan

ET

Groundwater Evapotranspiration

ET

Reference Evapotranspiration

ET

Unsaturated zone evapotranspiration

E

f

Groundwater exfiltration

EXTDP Extinction depth

EXTWC Extinction water content

FAO Food and Agriculture Organization

GHB General Head Boundary

GUI Graphical User Interface

h

Head in the aquifer

h

Head in streams

I Interception

IHM Integrated hydrological modelling

Kc Crop coefficient

KESDM Ministry of Energy and Mineral Resources of Indonesia

K

Horizontal hydraulic conductivity

Kv

Maximum unsaturated zone vertical hydraulic conductivity [UHC]

m a.s.l Meters above sea level

MAE Mean Absolute Error

ME Mean Error

MODFLOW Modular three dimensional finite-difference flow model

NS Nash-Sutcliffe efficiency

NWT Newtonian

P Precipitation

P

Actual infiltration

PET Potential evapotranspiration

P

Infiltration rate

q Stream discharge at the outlet of the catchment

Q Discharge

q

Dunnian saturated excess runoff

q

Lateral groundwater outflow

q

Groundwater loss to the stream

q

Hortonian runoff

q

Groundwater gain from the stream

R

Gross recharge

RH Relative humidity

RMSE Root Mean Square Error

R

Net recharge

R

Unsaturated zone recharge

R

Total runoff

RVE Relative Volumetric Error

SFR2 Surface flow routing package

SS Sunshine duration

SW-GW Surface water and groundwater

S

Specific yield

T

Maximum temperature

T

Minimum temperature

UPW Upstream Weighting Package

UZF1 Unsaturated zone flow package

WCsat Saturated water content

WS Wind speed

Y Overall model performance

1. INTRODUCTION

1.1. Background

Surface water and groundwater and (SW-GW) are very crucial for the well-being of humans and for the nature in general (Sophocleous, 2002). The significance of these resources is increasing through time, as the number of population is rapidly increasing (Fitts, 2002 & Anderson et al., 2015). Recently, research effort has been extended to understand the interaction between SW-GW and to quantify the flow between these resources, since, understanding of the link between them is needed for effective management of the resource (Lubczynski & Gurwin, 2005; Sophocleous, 2005; Hassan et al., 2014 & Ala-aho et al., 2015). The hydrological interactions between SW-GW occur through the unsaturated zone and by infiltration into or exfiltration from the saturated zone (Anibas et al., 2009).

The groundwater basin of D-T is located in the most developed area of Bali Island, Indonesia. The Island is located at 8

south of the equator and has an area of 5,380 km

, while the area for D-T basin is 2270 km

. The island of Bali has a variety of volcanic topographic units such as eroded early Quaternary volcanoes, active stratovolcanoes, thick tephra deposits, pyroclastic flow slopes and closed caldera lakes (Nielsen &

Widjaya, 1989; Kayane et al., 1993 & Purnomo & Pichler, 2015). These topography types can be classified into two categories. The Quaternary upper volcanic sequence and the Pliocene lower calcareous sequence (Figure 1). The latter is composed of a sequence of limestone, the "prepatagung formation", and calcareous sandstone. The thickness of this aquifer material is unknown, whereas the thickness of the Quaternary upper formation is considered to be up to 150 m. This Quaternary upper formation is composed of different materials of volcanic origin. It includes mainly unconsolidated sand & gravel, volcanic ash, lava flow, breccia, lahar, "pumic", clay and tuff (Rai et al., 2015 & Purnomo & Pichler, 2015).

Cole (2012) stated that groundwater is the most widely used resource in Bali Island as aquifers are characterized by high permeability. Mainly the Quaternary volcanic aquifer is the most productive and it can produce up to 7862 m

day

. Nielsen & Widjaya (1989), strongly recommended for the expansion of well systems for irrigation and municipal demands, due to their finding that 75% of the recharge reaches to the Quaternary upper formation called production aquifer. However, no research has been done regarding integrated hydrological model [IHM].

This research focuses on assessing the interaction between SW-GW resources and estimate the groundwater budget of the D-T basin. The interaction of SW-GW is best quantified through IHM (Lubczynski & Gurwin, 2005 & Kumar, 2015). MODFLOW-NWT is one of the models developed that link surface, unsaturated and saturated zone in a reliable manner. The model is working under ModelMuse Graphical User Interface (GUI) and merge with unsaturated zone flow package [UZF1] and stream flow routing package [SFR2]

(Niswonger et al., 2011 & Hassan et al., 2014). Therefore, in this study MODFLOW-NWT was conducted

to quantify the exchange flux between surface, unsaturated and saturated zones for the simulation period of

four years from 2009 to 2012. Groundwater study in D-T basin was selected because: (i) the basin is located

in the most developed area on the Bali Island and having high potential of groundwater resources but

sometimes suffers from water scarcity; (ii) it is representative of unconsolidated aquifer around the world,

i.e. it consists of volcanic aquifer with intergranular porosity and medium to high permeability; (iii) there is

a four years of hydrological data available; (iv) groundwater recharge has been studied by Nielsen & Widjaya (1989) using analytical methods and Artabudi (2012), that can be used for comparing the final findings.

Geological and cross section map of Bali Island

Figure 1: Geological map and cross section across the study area (Modified after Purnomo & Pichler, 2015).

1.2. Problem statement

Groundwater is the most popular water resources in Bali Island. However, Bali is affected by frequent crises

of water shortage (Purnomo & Pichler, 2015 & Straub, 2011). There is also a lack of knowledge and

mismanagement of SW-GW resources which causes decline of the water table, saltwater intrusion and

deterioration of water quality. Cole (2012) stated that the Island lacks official government statistics and

documentation regarding the available amounts of water resources. Furthermore, he mentioned the need

for reliable studies on SW-GW resources for successful water utilization and management. Besides, the

author of this study found that Bali groundwater has only been studied by Nielsen & Widjaya (1989) who

did recharge assessment but only by analytical methods and Artabudi (2012) who estimated groundwater recharge using remote sensing application. The interaction of SW-GW and groundwater recharge using IHM has not been studied in Bali Island.

1.3. Research setting

The overall objective of the study is to develop an IHM of D-T basin for management purpose.

Specific objectives:

1. To set up the numerical model in D-T based on the data from 2009 till 2012.

2. To calibrate steady-state and transient IHM of the D-T basin.

3. To estimate the water balance of the D-T basin.

4. To characterize the dynamics of SW-GW interactions in the D-T basin.

1.3.2. Research question

1. How to integrate various sources of data in the IHM?

2. What are the key components of spatiotemporal variability of the water balances in the D-T basin?

3. How does the water balance of the D-T basin vary on daily and/or yearly basis?

4. How do the SW and GW resources interact?

1.3.3. Novelty of the study

The findings of this study will increase the understanding of an unconsolidated aquifer of the D-T basin and can be used as an asset for management purpose by including the following novelties:

1. First-time use of IHM in the D-T basin. No research has been done in the area related to SW-GW interactions.

2. Use of daily streamflow measurement for IHM calibration.

1.3.4. Research hypothesis

It is hypothesized that the calibration of the integrated transient numerical model of the D-T basin gives a reliable estimate of the SW-GW exchange flux and groundwater budget.

1.3.5. Assumptions

- Eventual leakages across the bottom boundary of the volcanic aquifer have a negligible impact on the flow system of the modeled aquifer.

- Eventual lateral fluxes across watershed boundaries are negligible.

- Variable density groundwater flow, advection and dispersive salt transport have a negligible impact on the flow system of the modeled aquifer.

- Due to lack of groundwater abstraction data both the steady-state and transient-state IHM were calibrated

without groundwater abstraction data. Therefore, the impact of groundwater abstraction in the water budget

of D-T basin have a negligible impact.

2. RESEARCH METHOD AND MATERIALS

2.1. Study area

The D-T basin is located in the south-west of Bali Island, Indonesia with geographical coordinates of 8º39'S latitude and 115º13'E longitudes. It is part of the Bali Island with estimated area coverage of 2270 km

; whereas the whole Bali Island has an estimated area of 5,380 km

. The D-T basin has elevation variation from 0 to 2850 m a.s.l with the highest elevation located in the north and north-east of the basin and with the lowest elevation located in the south (Figure 2). The study mainly focused on the southern part of the area; where most population part of Bali lives in D-T basin, the study excluded Nusa Dua and Nusa Penida peninsula.

Location of D-T Basin

Figure 2: Location and elevation map of the D-T basin.

2.1.2. Monitoring network

The D-T basin, the southern Bali Island hydro-climatologically data such as rainfall, temperature, stream level, groundwater level, and pumping test data were monitored and available. These hydro-climatological data are available in the daily basis from 1

January 2009 to 31

December 2012 except for groundwater level. In D-T basin there are 21 meteorological stations, of all this 18 of them have sparsely located rain gauge stations and the remaining 3 are temperature, relative humidity, radiation and wind speed monitoring stations. Stream level at the southern basin outlet was monitored on a daily basis using 16 discharge gauge stations. The groundwater level was monitored by using 11 boreholes and dug wells (Figure 3).

Monitoring Stations

Figure 3: Monitoring stations of D-T basin 2.1.3. Climate

Based on average monthly rainfall, Bali has a pattern of monsoon type climate. Monsoon pattern occurs due to the air circulation changing direction every six months across the Indonesian region. In the area, month from April till the end of September is dry season and from October till the end of March is the wet season.

The air mass that brings the rain from the northwest equatorial wind in the wet season and the southeast wind from Australia in the dry season (Figure 2) cause the seasonal pattern in the area (Kayane et al., 1993).

Bali has an average annual rainfall of 2150 mm. Generally, during the rainy season, part of rainfall will

evaporate, other part will be taken up by the plant, some will run as overland flow and the remaining infiltrate

to the subsurface. About 30 to 50% of total rainfall has been estimated to infiltrate into the Quaternary

terrains. In fact, the percentage which infiltrate depends upon geological conditions, vegetation cover, land

use, and slope (Nielsen & Widjaya 1989). On average, the temperature ranges from 27-30

C (Figure 4) and

humidity 85-90 %. Moreover, in Bali soil and air temperature decreases at a lapse rate of 0.62

C/100 m

with elevation (Kayane et al., 1993).

Figure 4: Mean monthly rainfall; maximum, minimum, and mean temperature at station Kuta of D-T basin. For the location of the station see Appendix V.

2.1.4. Topography and land cover

The Island of Bali is a mountain chain that extends from the West to the East with volcanoes in Mount Batur (1717m) and Gunung Agung (3142 m) still active (Figure 1). The mountain chain that runs along the island of Bali causes morphological regions of Bali to be divided into several topographic and physiographic units (Purnomo & Pichler, 2015). Since the Northern part is high elevated, major drainage network is located in the south as well as the central part the D-T basin. In other respect, the most common land covers of D- T basin are Forest, 6%; bare soil, 9%; Agriculture, 68.3%; buildings, 13%; grassland and others, 3.4% (Figure 5).

20 22 24 26 28 30 32 34

0 4 8 12 16 20

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Temperature [oC]

Rainfall [mm month-1]

A. Mean monthy RF & T in 2009

Kuta_RF 2009 Max. Temp Min. Temp Mean. Temp

20 22 24 26 28 30 32 34

0 4 8 12 16 20

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Temperature [oC]

Rainfall [mm month-1]

A. Mean monthy RF & T in 2009

Kuta_RF 2010 Max. Temp Min. Temp Mean. Temp

20 22 24 26 28 30 32 34

0 4 8 12 16 20

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Temperature [oC]

Rainfall [mm month-1]

A. Mean monthy RF & T in 2009

Kuta_RF 2011 Max. Temp Min. Temp Mean. Temp

20 22 24 26 28 30 32 34

0 4 8 12 16 20

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Temperature [oC]

Rainfall [mm month-1]

A. Mean monthy RF & T in 2009

Kuta_RF 2012 Max. Temp Min. Temp Mean. Temp

Land use map Percent area ratio

Figure 5: Land use map and percent area coverage of D-T basin.

2.1.5. Hydrology

The D-T watershed boundaries are defined by topographic divides and delineate areas where surface water runoff drains into a common surface water body. The defined watershed boundaries are determined by science-based hydrologic principles, not favoring any administrative boundaries (Figure 6). Streams flow from the north side to the south direction following the valleys. Most of the streams are perennial (Rai et al., 2015), while their discharges originated from groundwater and surface runoff formed by precipitation.

The data from 16 stream gauges were available in a daily basis from 01/01/2009 to 31/12/2012.

Watershed boundaries of D-T basin

Figure 6: Watershed boundaries, stream segments and gauging location of D-T basin. For the name of each station see Appendix V (B).

0.48%

3.43%

5.99%

8.96%

12.82%

68.32%

100.00%

Water Grass Forest Bare Soil Buildings Aggriculture Total

2.1.6. Hydrogeology

Geologically, the study area is part of the “Sunda-Banda" volcanic islands arc. The arc is caused by the Indo- Australian and Eurasia plates. Since the late Tertiary, this process drives volcanism and produces a vast distribution of volcanic rocks (Figure 7). The Quaternary upper volcanic sequence which is also called unconsolidated layer and the Pliocene lower calcareous sequence also called consolidated layer, are the two dominant geologic formations of the area (Figure 1 & Figure 11). These volcanic rocks are rich in mafic minerals and exhibits considerable relief at the north (Purnomo & Pichler, 2015). Hills that surrounds Northern area forms surface water divides that coincide with the groundwater divides. Consequently, there is no flow contribution from outside of the basin and the outflow is the only through the streams discharge and lateral groundwater outflow towards the Indian ocean.

Pumping tests were conducted by Ministry of energy and mineral resources of Indonesia, locally "KESDM"

to evaluate the characteristics of the aquifer systems and groundwater potential of the basin. All the pumping tests were carried out in Quaternary layer and the Tertiary layer is considered to be an impervious layer.

Because of that, it was assumed that there is no hydraulic contact between the upper Quaternary and lower Tertiary layer.

D-T basin Geological map

Figure 7: Geology of the study area.

2.1.7. Previous studies in the area

Nielsen & Widjaya (1989), estimated groundwater recharge in southern Bali through five different

techniques for the year from 1980 to 1983. Based on analysis of well hydrographs they found that the

recharge value of 468 mm per annum, annual infiltration (≈25% of rainfall) gave 437 mm year

, base flow

separation gave 272 mm per annum, flow net analysis gave 492 mm per annum. Finally, they built an

analytical model based on land use and soil type and found recharge value of 645 mm per annum in light

soil, 538 mm per annum in medium soil and 376 mm per annuam in heavy clay soil. Based on the above

findings, they recommended the expansion of well systems for tourism, irrigation, and municipal demands.

They also showed the piezometric/potentiometric map of the area, in which groundwater flows from northern to the southern side of the basin (Figure 8). Additionally, Artabudi (2012) estimated the net recharge of Denpasar in two ways: (1) using “the Global Satellite mapping for precipitation” for the data from 2005 to 2009, the groundwater recharge value of 218 – 220 mmyear

was estimated; and (2) using in- situ rainfall data the groundwater recharge value of 650 – 660 mmyear

was obtained.

In the study area, the Quaternary volcanic aquifer which is also called the unconsolidated layer is the most productive and it can produce up to 7862 m

day

. This unconfined aquifer is characterized by high permeability and its volcanic rocks are characterized by intermediate hydraulic conductivity (K

). Sand and gravel of volcanic origin are highly permeable with transmissivity value of over 700 m

day

(Rai et al., 2015

& Purnomo & Pichler, 2015).

Southern Bali piezometric map

Figure 8: Potentiometric surface of Bali Island (Modified after Nielsen & Widjaya, 1989).

2.2. Data processing

Data collection is a prerequisite in groundwater modeling. Meteorological and hydro-geological data need to be collected as well as processed for effective model simulation of a given area (Kumar, 2015). In this study, the meteorological and hydrogeological data are available as shown in Table 1. Data such as groundwater depth, pumping test, hydrogeological and geological maps were collected from Ministry of Energy and Mineral Resources of Indonesia (locally "KESDM"). Most of the meteorological data such as daily rainfall, air temperature, relative humidity, streams level, and rating curves were collected from Ministry of Public Works (Locally called "DPPU"). Additionally, the land use map and 2012 daily rainfall data were collected from Biro of Geospatial Information (BIG) and Indonesian Agency for Meteorology, Climatology, and Geophysics (Locally called "BMKG") respectively.

Initially, in the data processing step, the missing rainfall data was filled using the coefficient of correlation

method (Teegavarapu & Chandramouli, 2005). Then, a double mass curve was used to check the consistency

of rainfall and stream discharge data (Searcy & Hardison, 1960). Afterward, all the available data (Table 1)

was converted in a way that MODFLOW-NWT can accept it. For instance, the available point observation of rainfall data was interpolated into rainfall map using kriging method. Kriged prediction was selected because it is a best linear unbiased predictor (Sterk & Stein, 1997; Hengl et al., 2007; Webster & Oliver, 2007; & Zhang et al., 2012). Spatially variable crop coefficient (K

), interception rate and infiltration rate were estimated using the 2009 land use map of the D-T basin. Using FAO Penman-Monteith method the available T

, T

, RH, WS, & SS together with spatially variable K

data were used to calculate potential evapotranspiration (PET). The hydraulic head was calculated from elevation and groundwater depths data.

Finally, the IHM was built and simulated in daily time steps for a four hydrological years from 1

January 2009 to 31

December 2012.

Table 1: Available data from the year 2009 to 2012 (Tmax - maximum temperature, Tmin - minimum temperature, RH – relative humidity, WS – wind speed, SS - sunshine duration, n.a - indicates that the data are not available within the study periods, EXTDP – Extinction depth, Kh – saturated hydraulic conductivity, and Sy - specific yield).

Required data Available data No. of

stations Required

Units Frequency of available data

Watershed boundary DEM - - -

Infiltration rate Rainfall 18 mday

daily

Potential

evapotranspiration T

, T

, RH,

WS, SS 3 mday

daily

Stream discharge Stream level &

rating curve 16 m³day

daily

Hydraulic head Groundwater

depth 11 m yearly

Groundwater

abstraction n.a. n.a. m³day

n.a.

Tidal head Tidal head 1 m daily

Model top elevation DEM - m a.s.l. -

Crop coefficient Land use map - - -

Interception Land use map - - -

EXTDP Land use map - - -

K

, S

K

S

14, 2 mday

, - n.a.

2.2.1. Watershed boundary

The D-T watershed boundaries were defined by topographic divides and delineating areas where surface water runoff drains into a common surface water body. The defined watershed boundaries were determined upon science-based hydrologic principles, not favoring any administrative boundaries. The SRTM (Shuttle Radar Topographic Mission) 90 m DEM (digital elevation model) was used to delineate these watershed boundaries. The SRTM 90 m DEM provided by NASA, has significant importance in determining the flow direction and in delineating catchment areas according to their stream flow. This data is available on the internet with free of charge by USGS, (http://srtm.csi.cgiar.org/). Initially, the available DEM data was filled into remove small imperfections in the data using hydrology spatial analyst tool in ArcGIS, then using the filled DEM map, a raster map of flow direction from each cell to its steepest downslope neighbor was generated. Then, a raster map of accumulated flow into each cell was generated. Afterward, the steam gauges were used as a snap pour points to the cell of highest flow accumulation within a 1000 m distance. Finally, the contributing area above a set of cells or watershed boundaries was obtained (Figure 6).

2.2.2. Precipitation

Precipitation is one of the first essential input data in IHM. Precipitation is used in UZF1 package to be

partitioned into runoff, infiltration, evapotranspiration, unsaturated-zone storage and recharge. Precipitation

together with potential evapotranspiration is the input driving forces of UZF1 packages of MODFLOW- NWT (Niswonger et al., 2006). The D-T basin rainfall was measured using tipping bucket at 18 sparsely distributed rain gauge stations (Figure 3 and Appendix V). From the total number of stations, only the station called Kuta has missed rainfall data.

Estimation of the Missing Precipitation Records

Rainfall records can be missed due to instrument malfunction, tree growth or other reason. Filling the missed data is one of the most important tasks to be carried out in many hydrological studies. The most commonly applied methods that are recently used for filling missed data are inverse distance, inverse exponential, the coefficient of correlation weighted method, and others (ASCE, 1996). In this study, the “coefficient of correlation weighted method” (CCWM) was used to fill missing data as proposed by Teegavarapu &

Chandramouli (2005) & Gómez (2007). Reliable estimation of missed data using CCWM is strongly dependent on spatial autocorrelation in which the data nearby are more similar than that of far apart. Once the correlation coefficient is determined, the missed data at a given station is calculated using equation 2.1.





 n 1 i rmi

rmi n

1 i Pi

Pm

(2.1)

where P

- precipitation at the base station m; P

- precipitation at station i, n – number of nearby stations, and r

- correlation coefficient of station i with nearest n stations. The correlation coefficient, r

, is obtained from SPSS by using the method known as “Pearson - Moment Correlation Coefficient” between closest stations (Appendix III).

Consistency of the precipitation records

Hydrological data consist of time series data that are collected at a particular location. The use of raw data without checking for consistency in IHM can bring lots of error and uncertainty. Therefore, before hydrological data are used in such study, they should be tested by the double-mass curve technique to ensure their reliability (Searcy & Hardison, 1960). The Double-Mass curve is plotted as cumulative of one station (y-axis) against the average cumulative of nearby stations (x-axis) as in Figure 9.

Figure 9: Double mass curve for a precipitation data (after Gómez, 2007).

Deviation from or break in the slope of the double mass curve means that there is a change in the consistency of proportionality between the variables or simply it shows the degree of change in the relation.

Such deviation might be due to gauge location, observation method or exposure. The precipitation records can usually be adjusted by coefficients determined from the double-mass curve (Equation 2.2).

Pb δb δa

Pa

(2.2)

where P

- adjusted precipitation, P

- observed precipitation, δ

- the slope of the graph at the time was observed and δ

- the slope of the graph to which records are adjusted.

2.2.3. Potential evapotranspiration

McMahon et al., (2013) define potential evapotranspiration [PET] as "the rate at which evapotranspiration would occur from a large area completely and uniformly covered with growing vegetation which has access to an unlimited supply of soil water, and without advection and heating effects." PET is one of the driving forces in the UZF1 package. In the UZF1 package, the PET is applied at the land surface and decreases linearly with depth down to the assigned extinction depth where evapotranspiration no longer occurs (Allen et al., 1998). In most case PET calculated using FAO Penman-Monteith method (Equation 2.3). The method is recommended by the scientific community as the best estimate of evapotranspiration with minimum error compared to other methods (Wang et al., 2012). FAO Penman-Monteith method is applicable if data such as daily air temperature, wind speed, relative humidity, atmospheric pressure and relative sunshine duration are available. Consequently, the method was adapted for reference evapotranspiration computation, since, the required data such as daily air temperature, wind speed, relative humidity, atmospheric pressure and relative sunshine duration was available from the two stations namely called station Sanglah and Kuta.

Before ET

computation using the equation 2.3, the correlation between T

, T

, and T

for the available stations were constructed. A good coefficient of correlation between the station means that temperature is uniform spatially. Then, the computed ET

would be spatially uniform but temporally variable. In the case of low coefficient of correlation between T

, T

, and T

of the different stations kriging interpolation method would be applied to generate the spatially variable ET

values.

Note that, the computation equations to estimate evapotranspiration as shown in Equation 2.3 and Appendix II are the one that were applied in this study based on the available data. There are many possible ways of determining one parameter based on the availability of data in a given area. For detail study interested readers are referred to Allen et al.,(1998) & Raes & Munoz (2009).

2) U

* 0.34 γ(1 Δ

a) s e 2(e 273U T

* 900 γ n G) (R

* Δ

* 0.408 ETo





 







(2.3)

where ET

- reference evapotranspiration [mmday

], ∆ - slope vapour pressure curve [kPa°C

], R

- net

radiation at the crop surface [MJm

day

], G - soil heat flux density [MJm

day

], 𝛾 - psychrometric constant

[kPa

C

], T - mean daily air temperature at 2 m height [°C], U

- wind speed at 2 m height [ms

], e

vapour

pressure [kPa], e

- actual vapour pressure [kPa], e

-e

- saturation vapour pressure deficit [kPa], R

again can

be calculated using the equations as in Appendix II.

The computed ET

needs to be converted into PET since MODFLOW-NWT requires PET rather than ET

. The single and dual crop coefficient (K

) are the two approaches to converting ET

into PET (Allen et al., 1998). The former approach calculates as in equation 2.4; in which the K

depends on crop type and growing stage, whereas, in the latter case the K

splits it to two other factors. The first factor for evaporation and second factor for transpiration difference between crop reference surface (Zehairy, 2014 &

Weldemichael, 2016). Dual crop coefficient requires sufficient data about crop/vegetation and soil and hence for this study spatially variable but temporally invariable single crop coefficient approach was applied.

K c o * ET

PET  (2.4) where PET - potential evapotranspiration [mmday

], and Kc- crop coefficient [-].

The spatially variable K

in this study area was estimated based on the land use and vegetation cover. The D-T basin mainly covered by agricultural land, forest, bare soil, and grass. The K

value for bare soil which is about 9 % of the study area was assigned as 0.61. About 6 % of the area covered with forest and the K

of trees was assigned as 1.0 (Allen et al., 1998). Agriculture land that covers 68% of the total area and the majority of the agricultural land was covered by rice field. The K

values of rice during initial, crop development, mid-season, and late season stages were 0.62, 0.75, 1.16 and 0.67 respectively (Choudhury et al., 2013). It was assumed that the K

value for rice field to be during crop development stage. The K

value of 0.75 was used for agriculture land. After defining the crop coefficient for each land cover and vegetation type, the spatially variable K

was obtained for 500 * 500 grid cells. Then, the spatiotemporal variability of PET was calculated using Equation 2.4.

2.2.4. Interception and infiltration rate

Interception rate is the amount of rainfall that is retained by vegetation above the surfaces. The rate depends mainly on the rainfall duration, density, and morphology of vegetation cover. The interception loss from a given land cover can be calculated by using Equation 2.5 (Zehairy 2014 & Weldemichael 2016).

other ) A other * f I

A f * (I

* P

I   (2.5)

where I - canopy interception per grid cell [mday

], P - precipitation [mday

], I

, and I

- interception loss rate by forest, and agriculture respectively in [%] of precipitation, and A

other A

- area ratios coverage per grid for forest and agriculture respectively.

The D-T basin interception loss was calculated as in Equation 2.5. The spatial variability of interception rate was calculated per grid cells and subtracted from spatially variable precipitation rate to get spatially variable infiltration rate. Several kinds of the literature suggested interception ratio, i.e. I

and I

based on land cover. For instance, Ghimire et al., (2012) determined rainfall interception by natural and planted forests in the middle mountains of central Nepal. According to their finding, interception loss for the evergreen natural forest was 22.4% of precipitation. Since similar forest type is present in the D-T basin (Heim, 2015), this study adapted the value of interception ratio 22.4% of precipitation as a value for forest cover. Van Dijk

& Bruijnzeel (2001) also studied rainfall interception in upland West Java, Indonesia for mixed agricultural

cropping system involving cassava, maize, and rice. According to their finding interception loss for the

agricultural land was 14.4% of precipitation. The final finding of Van Dijk & Bruijnzeel (2001) for mixed

crops was conducted in this study as interception loss rate for agriculture farmland. Finally, Corbett et al.,

(1968) determined rainfall interception by annual grass and chaparral and found that grass has an

interception loss rate of 6.5% of the total rainfall. This finding for grass was used in this study as interception loss from grassland. The values are summarized in Table 2.

Table 2: Interception loss rate that used in D-T Basin

D-T Land cover Interception loss Adapted literature

Forest cover 22.4% Ghimire et al., (2012)

Agriculture [mainly rice] 14.4% Van Dijk & Bruijnzeel (2001)

Grassland 6.5% Corbett et al., (1968)

As stated earlier the infiltration rate was calculated as the difference between rainfall and interception rate (Equation 2.6). Infiltration is the amount of water that percolate to the unsaturated zone. The infiltration rate highly depends on the vertical hydraulic conductivity and degree of saturation of the unsaturated zone.

The higher the vertical hydraulic conductivity, the easy through which water can pass through the soil, then the higher the infiltration rate would be and vice versa. When the infiltration rate is higher than the vertical hydraulic conductivity, the water hardly passes through the soil pores, then the excess rainfall will be redirected to the streams.

I P

P

  (2.6) where P

- infiltration rate per grid cell [mday

], P - precipitation [mday

]. This infiltration rate is one of the driving forces that was used in the UZF1 package to estimate the unsaturated zone storage and unsaturated zone evapotranspiration, and groundwater recharge (Niswonger et al., 2006).

2.2.5. Stream discharge

According to Braca (2008), the stage-discharge relation or rating curve in an open channel flow is used to convert the series of stage records into discharge records. Similarly, it is used to convert forecasted flow hydrographs into stage hydrographs. The relationship is highly affected by section and channel controls.

The latter is due to hydraulic properties and roughness of the surface downstream. Such hydraulic properties can be channel size, shape, slope, and curvature. On the other hand, the former can be caused by natural or man-made processes. Rock ledge, sand bar, and debris are grouped under natural processes. Dam, Weir, Flume, and Spillway are grouped under man-made processes. Because of the underlining reasons, the knowledge of channel features and a time series of stage and discharge together with low and high flow measurements are required to construct efficient rating curves. The stage-discharge relation can be calculated using the simplified forms of Manning equation (ISO, 2010) as shown in Equation 2.7.

)

( h a C

Q   (2.7)

where Q - discharge [m

day

], h - stage [m], a - gauge height of zero flow [m], h  a - effective depth of

water on the control, C - calibration coefficient [m

day

], α - slope of rating curve [-].

The D-T basin has 16 sparsely located discharge stations (Figure 6 and Appendix V). The data that was collected from "DPPU" was daily stream level for the year 2009 to 2012 together with rating curve for each station but it lacks stream discharges data. The availability of rating curve that was constructed using HYMOS (Sankhua & Srivastava, 2011) saves the time to convert the stages into stream discharges.

Therefore, the four-year daily stream level data were converted into stream discharge for each station and

then missing discharge data and its consistency was adjusted in the same manner as precipitation – section 2.1.1. Afterward, the final result was used as a state variable in SFR2 package (Niswonger & Prudic 2005).

2.2.6. Hydraulic properties

The pumping tests were used to obtain the aquifer characteristics such as transmissivity and storativity. Such data together with other data was used to build a numerical modeling and evaluate the resources in the given area. Table 3 shows the aquifer characteristics of D-T basin, this test value was used as an initial value in the model and then adjusted during the model calibration.

Table 3: Aquifer characteristics value of the basin (T - transmissivity, Sy – specific yield, Kh - horizontal hydraulic conductivity, SWL - surface water level). The location of bores can be seen in Appendix I.

ID Latitude Longitude

T [m

day

]

S

K

[mday

] SWL

[m]

Thickness [m]

Bottom elevation [m]

PZ1 8°25'29.21" 115°12'50.48" 1821.6 67.5 100.1 27.5 126.5 PZ2 8°30'56.37" 115°27'02.07" 2340.1 36.2 50.4 64.6 11.6

PZ3 8°20'9.69" 115°21'13.83" 134.8 2.7 89.8 50.2 140.1

PZ4 8°12'42.08" 115°17'17.63" 57.6 1.2 148.2 47.6 195.5 DP2 8°37'14.17" 115°14'0.91" 390.2 0.25 2.9 18.8 131.2 149.5 DP3 8°23'19.05" 115°13'20.55" 690.5 5.4 22.3 128.9 150.1 DP4 8°30'21.10" 115°12'33.51" 550.2 4.3 21.6 128.4 150.0 DP5 8°39'19.45" 115°11'46.89" 760.1 0.24 5.9 22.1 127.9 150.5 DP6 8°38'36.55" 115°13'56.06" 385.8 2.8 14.1 135.9 148.5 DP7 8°33'55.05" 115°10'50.25" 280.5 2.1 19.1 130.9 147.7 DP13 8°37'26.58" 115°14'09.05" 298.6 2.1 10.8 139.2 151.2

sbo1 8°31'44.41" 115°12'5.08" 7.1 0.1 62.0 62.2 130.2

sbo2 8°27'48.60" 115°11'15.25" 8.9 0.2 62.0 70.5 129.9

sbo3 8°32'18.33" 115°09'56.42" 5.4 0.1 4.2 111.0 115.2

sbo4 8°19'30.60" 115°22'00.83" 8.5 0.7 31.6 118.3 151.1 sbo9 8°22'33.38" 115°17'06.45" 110.1 0.8 9.9 140.1 151.0 sbe10 8°32'20.20" 115°06'54.27" 43.1 0.4 38.9 111.1 151.4

In the study area, all groundwater level measurements are undertaken in Quaternary deposits and show strong spatial variation to topographic altitude differences. However, such data were available neither daily nor monthly rather monitoring records show one record per year besides, the monitoring stations have low spatial coverage over the study area (Figure 3).

2.2.8. Groundwater abstraction

The D-T basin groundwater abstraction data could not be found in public services as well as from Ministry

of Energy and Mineral Resources of Indonesia, locally "KESDM". It was found out that the data is

confidential and the ministry is not willing to share the data even for study purposes. Because of the

underlining reason, the author builds the model without groundwater abstraction. Therefore, the result of

this model should be used with caution in case future model implemented using groundwater abstraction

data.

2.3. Modeling flow chart

The overall activities going to be applied to answer the research questions and to come up with the targeted objectives is summarized in the flow chart (Figure 10).

Figure 10: Methodology of flow chart.

2.4. Conceptual model

According to Anderson & Woessner, (1992), there are three essential steps in modeling protocol. These are

(1) “to establish the purpose of the model”, (2) “formulation of the conceptual model of the system”, and

(3) “formulation of the numerical model of the system”. The purpose of the model is for better

representation of the flow between surface, unsaturated zone and saturated zone. It is useful in

understanding the interaction and estimating the water balance in the D-T basin. A conceptual model is a