Groundwater Modeling Irwan Iskandar, PhD

Groundwater Modeling Irwan Iskandar, PhD KK Eksplorasi Sumberdaya Bumi Teknik Pertambangan Fakultas Teknik Pertambangan dan Perminyakan ITB

Apakah yang dikerjakan dalam pemodelan What hidrogeologi? Hydrogeologist do? • Kuantifikasi Distribusi (keterdapatan) • Kuantifikasi aliran (recharge-discharge) • Prediksi dan simulasi pengambilan airtanah • Interaksi air permukaan dan airtanah • Interaksi air dengan batuan (hidro-geo-kimia) • Perencanaan dewatering

• Perencanaan eksplorasi - produksi (geothermal, CBM, Migas) [email protected]

Hidrogeologi dalam bidang pertambangan (tambang terbuka, tambang bawah tanah, geothermal, perminyakan ): Why do we need Hydrogeologist in our business?

[email protected]

Mine environment: Could we ask the hydrogeologist, “How”: could we control our mine water? so people can say our business is ‘green’ and eco-friendly Bagaimana pengelolaan air (drainase), desain pompa, dan kontrol kualitas air di tambang? Parameter apa saja yang perlu dipertimbangkan?

[email protected]

Tujuan (Pemodelan Hidrogeologi Tambang Terbuka) • Pola aliran airtanah dan air permukaan pra, selama, dan pasca tambang (kuantitatif) • Kuantifikasi jumlah air yang harus diatur dalam penyaliran

tambang • Environmental impact (prediksi, verifikasi lingkungan) • Desain sump, ditch, pump, pond

• Depressurization • “Water Management”

[email protected]

Air ke dalam Pit Tambang

[email protected]

Pola Aliran Air Permukaan Pra Tambang

Kursus Hidrogeologi IAGI Bandung, 25 – 26 April 2013

[email protected]

Pola Aliran Air Permukaan Pasca Tambang

[email protected]

Pola Aliran Air Airtanah dan Penurunan Air Tanah Akibat Tambang Head Equipotential Line

Perkiraan perimeter Pit

Groundwater Table Drawdown (max 50 in pit high wall

Drain Hole Desain • • • • • •

Metode drain Spasi Kedalaman lubang Time frame Evaluasi drain Efek ke kesetimbangan lingkungan airtanah di sekitarnya [email protected]

Slope Remediation

Infiltration Control

Crest Drain

Well Control

Toe Drain [email protected]

Pola Aliran Air Airtanah dan Penurunan Air Tanah Akibat Tambang

[email protected]

Pola Aliran Air Airtanah dan Penurunan Air Tanah Akibat Tambang

Kursus Hidrogeologi IAGI Bandung, 25 – 26 April 2013

[email protected]

Kuantifikasi Airtanah yang keluar dari dinding tambang Rate Dewatering Time [day] in m3/day

Rate Dewatering in liter/sec

Volume Cumulative [m3]

365

-196.73

-2.28

-71808

730

-392.42

-4.54

-286464

1095

-421.41

-4.88

-461440

1460

-403.90

-4.67

-589696

1825

-368.43

-4.26

-672384

2190

-340.22

-3.94

-745088

2555

-328.64

-3.80

-839680

2920

-312.28

-3.61

-911872

3285

-306.42

-3.55

-1006592

3650

-290.47

-3.36

-1060224

[email protected]

Kuantifikasi Airtanah dan Air Permukaan yang keluar dari dinding tambang

PIT

PIT (Kedalaman pit 150 m) PIT (Kedalaman pit 70 m)

Debit Air Limpasan

Debit Air Tanah

Total Debit Maksimum Air Masuk Pit

(m3/detik) 1.81

(m3/detik) 0.004

(m3/detik) 1.814

1.75

0.003

1.753

[email protected]

Desain Saluran Drainase Freeboard = 10-30 cm

X=lebar muka saluran

h = kedalaman saluran basah

Z = kemiringan saluran

5/3

1/ 2

A S Q 2/3 n.P

B = lebar dasar saluran

Debit MAX (liter/detik)

Kapasitas Debit Maks Saluran Air 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

Tinggi Muka Air (m) Kursus Hidrogeologi IAGI Bandung, 25 – 26 April 2013

0,8

0,9

1

[email protected]

Desain Saluran Di Sekitar Dinding Tambang dan Inpit Sump

Kursus Hidrogeologi IAGI Bandung, 25 – 26 April 2013

[email protected]

Tabel Balance Hasil Simulasi

Arah Aliran Air Tanah dan Kontur Head

[email protected]

Simulasi Dewatering Kursus Hidrogeologi IAGI Bandung, 25 – 26 April 2013

[email protected]

Pongsesa Infiltration Pond

Harapan Infiltration Pond

Concentration Cr6+ ppm

Simulasi contaminant transport

[email protected]

3D-Permeability Modeling • 3D-spatial distribution of permeability in a block modeling • Combination K value (primary variable) and RQD value as (secondary) • Replace previous layer (stratigraphical modeling) to grid based model

Field work

Conceptual model

Groundwater modeling

SIMULASI NUMERIK dalam HIDROGEOLOGI

Governing Equation 2h 2h 2h   2 0 2 2 x y z

•Steady state

 h  h  h h ( K x )  ( K y )  ( Kz )  S s x x y y z z t

•Transient

INTI SIMULASI NUMERIK • Teknik untuk mencari solusi persamaan differensial

• Merupakan pendekatan • Persamaan differensial didekati dengan persamaan linier simultan (persamaan matrix)

DUA JENIS UTAMA SIMULASI NUMERIK

1. Finite Difference (beda hingga)

2. Finite element (elemen hingga) Kelas kita: Finite Difference

FINITE DIFFERENCE • Introduced firstly by Richardson (1910) • The basic idea: to replace derivative at a point by ratio of change over a small but finite interval

LANGKAH-LANGKAH DALAM METODA FINITE DIFFERENCE

1. Pendekatan Terhadap Differensial h(x) B P

A

Deret Taylor:

x-h

x

x+h

x

dh 1 2 d 2 h 1 3 d 3h hx  x   h( x)  x  x  x  ... 2 3 dx 2 dx 6 dx

(1)

dh 1 2 d 2 h 1 3 d 3h hx  x   h( x)  x  x  x  ... 2 3 dx 2 dx 6 dx

(2)

hx  x   h( x  x)  2h( x)x 2

2

d h 4  O (  x )  ... 2 dx

+ (3)

Abaikan orde O(x4):  d 2h  14 Abaikan  h( x  x)  2h( x)  h( x  x) 2):  dx 2  0(h   x  x x

(4)

Error orde h2

(1)-(2):

1  dh  h( x  x)  h( x  x)     dx  x  x 2x

(5)

Error orde h2 Pendekatan slope pada P dengan menggunakan garis AB: CENTRAL DIFFERENCE

Dari Pers. (1): Abaikan orde h2:

1  dh  h( x  x)  h( x)     dx  x  x x

Forward Difference (CD)

1  dh  h( x)  h( x  x)     dx  x  x x

Backward Difference (BD)

2. Pembagian dalam sistem grid 1-D

2-D y

i-1

i+1

i h

i,j+1

i-1,j

i,j

i+1,j

h i,j-1

x

3. Penulisan pendekatan pers. differensial pada setiap titik dalam sistem grid 2-D

  2h  1  2   2 hi 1, j  2hi , j  hi 1, j   x  x  x x

y

i,j+1

i-1,j

i,j

1  h  hi1, j  hi1, j  CD :     x  x  x 2x

i+1,j

h

1  h  hi1, j  hi, j  FD :     x  x  x x

y i,j-1

x

x

1  h  hi, j  hii, j  BD :     x  x  x x

Dalam arah y   2h  1  2   2 hi , j 1  2hi , j  hi , j 1   y  x  x y

 h  1 hi, j 1  hi, j 1  CD :     y  y  y 2y  h  1 hi, j 1  hi, j  FD :     y  y  y y  h  1 hi, j  hi, j i  BD :     y  y  y y

5. Penyelesaian Persamaan Linier Simultan (Persamaan Matrix) • • •

Eliminasi Gauss Algoritma Thomas Iterasi: 1. Jacobi 2. Gauss - Seidel 3. SOR

Exercise 1

x 8,04

15 menit 8,18

8,36

8,53

Tugas: Cari h2,2, h3,2, h2,3, dan h3,3

h1,4

h2,4

h3,4

h4,4

7,68

?

?

8,41

h1,3

h2,3

h3,3

h4,3   2h  1  2   2 hi 1, j  2hi , j  hi 1, j  x  x  x x

y = x

7,19

 2h  2h  2 0 2 x x





?

?

h1,2

h2,2

h3,2

h4,2 2 h  1  2   2 hi , j 1  2hi , j  hi , j 1   y  x  x y

6,82

7,56

7,99

8,29

h1,1

h2,1

h3,1

8,33

h4,1 Steady state: hi,j sama pada tiap waktu n

Transient Condition   h    h    h  h    K zz  K xx    K yy   W  Ss x  x  y  y  z  z  t

Metode Beda Hingga Salah satu solusi pemecahan masalah persamaan dalam kondisi transient

hi,j

= (1/4) (hi-1,j + hi+1,j + hi,j-1 + hi,j+1)

Steady state: hi,j sama pada tiap waktu n

hin1,1j  hin1,1j  hin, j 11  hin, j 11  4hin, j 1  1/ T   Sa 2  1/ t   hin, j 1  hin, j  Transient: pencarian hi,j pada tiap waktu n

Pemodelan Lapisan (Pemodelan Parameter Hidraulik) 1. Layer Based  LPF (Later Package File)

2. Grid Cell Base

Physical Model (Layer Based) • Pembagian Unit Hidrostratigrafi

• Setiap unit lapisan hidrostratigrafi diterjemahkan sebagai “layer” • Setiap layer relatif homogen (K, S, θ dan parameter lainnya)

• Lapisan bisa datar, miring ataupun membentuk antiklinorium • Korelasi bisa manual oleh hydrogeologist atau dengan software  pembentukkan kontur struktur “layer”

Physical Model (Layer Based) • Easy when hydrostratigraphical unit has been defined

• Limited number of aquifer parameter e.g. 1 value in each unit. • Suitable in layered / sedimentary rock • Number of cells relatively low

• Time of simulation relatively short

Topo: [ID, x, y, z] Unit : akuitard, impermeabel, K<10-9 m/detik Unit : akuitard, impermeabel, K<10-7 m/detik

Unit : akuifer, permeabel, K> 10-5 m/detik

Unit : akuitard, impermeabel, K<10-9 m/detik

Unit : akuitard, impermeabel, K<10-7 m/detik

Bottom Layer 1: [ID, x, y, z, K, S, θ]

Bottom Layer 2: [ID, x, y, z, K, S, θ] Bottom Layer 3: [ID, x, y, z, K, S, θ] Bottom Layer 3: [ID, x, y, z, K, S, θ] Bottom Layer 4: [ID, x, y, z, K, S, θ]

Bottom Layer 5: [ID, x, y, z, K, S, θ]

Bottom Layer 5: [ID, x, y, z, K, S, θ]

Physical Model (Database) Hole ID

ID

X

Y

Z

Bottom 1 Bottom 2 Bottom 3 Bottom 4 Bottom 5 Bottom n

Data (Physical Model) • Topografi X, Y, Z (ASCII file) atau (.DXF) 3D polyline

• Log Bor ID, X, Y, Z (ASCII) atau (.DXF) 3D polyline

• Hydrologic Features (River, Lake, Sea, Pond, Stream) in dxf

• Hydrogeological Parameter (K, S, θ,) dalam 3 D data (X, Y, Z, K, S, θ) atau (ID, Layer, K, S, θ)

Physical Model (Grid Cell Based) • stratigraphy unit is not a must

• Each cell unit is translated to valued grid • Each cell is has one value (K, S, θ and other parameter) • Block model of the grid would hydrostratigraphical pattern

depend on the structure • Adjustment and interpolation were made by Hydrogeologist • A lot of Number of cells

• Time of simulation is relatively long

Hydraulics Conductivity and RQD Higher RQD

Lower RQD

4,00E-06

Higher Hydraulic Conductivity

Hydraulic Conductivity (m/s)

3,50E-06

y = 8E-06x - 5E-06 3,00E-06 2,50E-06 2,00E-06 1,50E-06 1,00E-06 5,00E-07

Lower Hydraulic Conductivity

0,00E+00 0,5

0,6 Linear (RQD vs K)

0,7

0,8

0,9

RQD Factor = (1-RQD/100)

Poly. (RQD vs K)

1

Contoh Grid Cell Based Model Dimension Easting (column) Northing (row)

Min Max Min Max Min Elevation Max

403500 405500 306000 308000 400 1200

Distance

Grid Size

Grid Sum

2000

10

200

2000

10

200

800

10

80

10 x 10 x 10 meter

Total Block Grid

3,200,000

Contoh Grid Cell Based Data Bor ID GW-01

GW-02 GW-03

GW-04

GW-05

GW-06

CGW-07

CGW-08

X (m)

Y (m)

24300 24300 24300 24300 24300 24300

-1200 -1200 -1200 -1200 -1200 -1200

23000 23000 23000 23800 23800 23800 23800 23800 25000 25000 25000 25000 25000 25000 22500 22500 22500 22500 22500 22500 22500 23200 23200 23200 23200 22700

1000 1000 1000 0 0 0 0 0 -1100 -1100 -1100 -1100 -1100 -1100 -2300 -2300 -2300 -2300 -2300 -2300 -2300 -2000 -2000 -2000 -2000 -1300

Depth (m) 6.4 12.5 20.6 28.0 33.1 35.6 37.2 3.3 6.5 7.7 12.4 18.1 24.4 36.4 39.4 6 11.1 21.2 23.2 29.1 47.1 6.5 16.7 21.5 26 30.5 35.9 39.5 3.5 12.5 32.9 40.3 6.2

k (m/s)

log k

1.14× 10-5 5.08× 10-6 2.77× 10-4 1.81× 10-5 5.47× 10-4 1.82× 10-5 2.89× 10-4 1.50× 10-5 1.85× 10-5 8.22× 10-5 2.89× 10-5 3.39× 10-4 7.29× 10-5 1.25× 10-4 4.78× 10-4 2.08× 10-5 3.24× 10-5 7.41× 10-5 1.85× 10-5 1.50× 10-6 2.08× 10-5 1.10× 10-4 2.55× 10-5 3.88× 10-5 6.74× 10-5 2.35× 10-5 1.04× 10-4 5.65× 10-6 3.44× 10-5 7.13× 10-4 1.27× 10-4 2.69× 10-4 3.89× 10-5 -6

-4.943 -5.294 -3.557 -4.742 -3.262 -4.739 -3.539 -4.822 -4.732 -4.085 -4.538 -3.469 -4.137 -3.903 -3.320 -4.681 -4.489 -4.130 -4.732 -5.822 -4.681 -3.957 -4.593 -4.411 -4.171 -4.628 -3.981 -5.247 -4.464 -3.147 -3.895 -3.570 -4.409

Data (Physical Model) Grid Cell Based

Statistical summary of log transformed k value Number of data

95

Minimum value

-5.57

Mean

-4.40

First quartile

-4.74

Standard deviation

0.68

Median

-4.46

Coefficient of variation

-0.15

Third quartile

-3.40

Skewness

-0.20

Maximum value

-3.15

Data (Physical Model) Grid Cell Based

• Estimasi parameter Hidraulik di daerah yang tidak ada data • 3D/4D space-time distribution modeling

Geostatistic-ordinary kriging

m

Zˆ ( xk )   i Z ( xi ) i 1

Zˆ ( xk )

The spatial correlation dianalisis menggunakan semivariogram γ(h) n

  ˆ(h )    ˆ(h j 1

ij

ik

1 n 2 ˆ (h)    Z ( X )  Z ( X  h) 2n i 1

Sill

Nugget (may be zero)

i

g

= Data Points = variogram model

Range

Lag or Separation Distance/data

)

3D Spatial Distribution

North

Co Located Co Kriging (CK) Grid kecil adalah RQD di batuan point/grid (size = 10 x 10) n1

n2

i 1

j 1

zˆ1 (x 0 )   λ1i z1 (x 1i )   λ2i z 2 (x 2 j ) Grid besar(merah) adalah data parameter hidraulik Size = 125 x 125

n1

 i 1

1i

 1,

n2

 i 1

2i

 0,

1 1 N (h) γ12 (h)   z(  z 2(x i )) (z( 1 xi  h ) 1 x i )) (z 2(x i  h ) 2 N (h) i 1

 (|h|) 450 400 350 300 250 200

γˆ(h ) sph

150 100 50 0 0

1000

2000 |h|

3000

4000

 3h  h 3  C 0  C1   3  2a  2a   C 0  C1 for h  a

  for 0  h  a  

3D Spatial Distribution Konduktivitas hidraulik (K)

Contoh Grid Cell Based Data dengan data struktur (RQD)

Conductivity Distribution (section view) A C

D

B

A

B

C

D

Conductivity Based on Lithology Distribution (section view) B

A

A

B

EXAMPLE • Case Study Underground Mine

Hydraulics Conductivity and RQD

Conductivity Distribution (section view) A C

D

B

A

B

C

D

Initial Head Condition (Based On 2011)

Mine Design N

N

N

Mining Drain Scenario (Assumption) • Mining Development assumption finished in two period, always open during mining activity. • Mine Production assumption finished in seven period. Closed every period with filling material, where filling material assumption is impermeable.

Drain Scenario (Assumption) File

Volume (m3)

Surface Area (m2)

Accumulatif Volume Opening (m3)

Dev1

243,342.00

190,504.00

243,342.00

Dev2

215,307.00

177,459.00

458,649.00

Mine1

244,681.00

92,227.00

703,330.00

Mine2

243,566.00

88,920.00

702,215.00

Mine3

246,932.00

109,066.00

705,581.00

Mine4

251,695.00

91,688.00

710,344.00

Mine5

241,731.00

98,695.00

700,380.00

Mine6

231,848.00

91,410.00

690,497.00

Mine7

234,107.00

120,843.00

692,756.00

Mining Drain Scenario (Assumption) Mining Development

Mine opening

Mine opening

Mine opening

Closed mining

Closed mining

Mine opening

Closed mining

Mine opening

Mine opening

Closed mining

Mine opening

Closed mining

Closed mining

Model Scenarios Model based on Conductivity Distribution (K values is not depend the lithology condition) 1. K values average 10-7 m/sec, rock condition is in low RQD mostly (worst scenario) 2. K values average 10-8 m/sec, rock condition is in moderate RQD mostly (moderate scenario) 3. K values average 10-9 m/sec, rock condition is in fair RQD mostly

MODEL TYPE- 1 K VALUES AVERAGE 10-7 M/SEC, ROCK CONDITION IS IN LOW RQD MOSTLY (WORST SCENARIO)

Model Calibration (Head Calibration on Steady State)

Observation Head (2011) vs Calculation Head (Model) – Steady State Condition

Model Calibration (Contour Head Calibration on Steady State)

Head contour Based On Observation Data (Env Field Measurements, 2011)

Head contour Based On Model Calculation (Steady State Condition)

Change of Groundwater Head and Flow (Plan View) Year: 1, 2, 3, 5,7 and 9

Change of Groundwater Head and Flow (Plan View) Year: 10, 20, 30, and 40

Head and Groundwater Flow (Section View) Year: 1, 2, 3, 5,7 and 9 N

N

N

N

N

N

Observation Well Location

2

4

1 3

Head vs Time (Observation Well) 880

Groundwater Head (m)

860

840

820

800

780

760 1

2

3

4

5

6

7

8

9

10

12

14

16

18

20

Year OBSERVATION 1/AInterpolated

OBSERVATION02/AInterpolated

OBSERVATION3/AInterpolated

OBSERVATION4/AInterpolated

22

24

26

28

30

Head Drawdown • Head Drawdown Radius Maximum: ± 139 m • Head Drawdown Depth Maximum: ± 127 m • Water table drawdown may not affected to unsaturated zone above the aquifer • All Drawdown will recover in ± 30 - 45 years after mine closure