SIMULASI INTRUSI AIR LAUT DENGAN MENGGUNAKAN SOFTWARE MATLAB R2010b SEAWATER INTRUSION SIMULATION USING MATLAB R2010b Agustina Cindy1 dan Idris Maxdoni Kamil2 Program Studi Teknik Lingkungan Fakultas Teknik Sipil dan Lingkungan Institut Teknologi Bandung Jl. Ganesha 10 bandung 40132 1
[email protected] dan
[email protected] Abstrak : Akuifer air tanah merupakan bagian yang penting di kawasan pantai. Di beberapa daerah dimana tidak ada air permukaan baik dari sungai maupun reservoir, air tanah merupakan satu-satunya alternatif yang dapat digunakan sebagai air baku selain dari air hujan. Namun di kawasan pantai rawan terjadi intrusi air laut akibat eksploitasi air tanah berlebih maupun kurangnya jumlah resapan air yang masuk ke dalam tanah. Untuk mencegah terjadinya intrusi air laut, banyak dilakukan studi untuk memodelkan intrusi. Dalam penelitian ini studi terbatas untuk pulau kecil, khususnya Pulau Pari (terletak di bagian utara kota Jakarta), dengan asumsi bahwa air asin dan air tawar merupakan dua buah fluida dengan densitas berbeda yang terpisahkan oleh bidang batas tegas sehingga tidak mengalami pencampuran. Tujuan penelitian ini untuk mengoptimalkan debit pengambilan air tanah dan mencari batas kritis bidang kontak. Persamaan diferensial parsial yang digunakan dalam pemodelan diselesaikan dengan metode beda hingga. Hasil simulasi menunjukkan bahwa semakin banyak diskretisasi, semakin luas area yang dapat dimodelkan. Dari hasil penelitian dengan beberapa ukuran grid, didapatkan pula bahwa nilai debit optimal tidak boleh lebih dari 0,03m3/detik. Apabila debit melebihi batas yang seharusnya akan menimbulkan dampak pada posisi bidang batas (interface). Dari hasil simulasi untuk dimensi grid 6x6, posisi kritis bidang batas terbentuk ketika debit pemompaan sebesar 0,036m3/detik. Kata Kunci : intrusi air laut, bidang batas tegas, metode beda hingga
Abstract : Groundwater aquifers are important resource in coastal regions. In many areas where there is no fresh surface water from rivers or reservoirs, the development of groundwater resources is practically the only alternative to storage of rainwater. However, the coastal aquifers are very vulnerable to the seawater intrusion through the overdraft of groundwater exploitation or insufficient recharge from upstream. In order to control seawater intrusion, seawater intrusion management models have been studied. In this study, study was limited for small islands, in particular the Pari Island (located in the Northern of Jakarta city), with the assumption that saltwater and freshwater were fluids with different density that are seperated by a sharp interface which mean these two liquids are not soluable in each other. This study addresses to optimize the groundwater pumping and to obtain the critical position for the interface. Partial differential equation used in this study was conducted by using finite difference method. The result shows that the more discretization was made, the greater the area to be modeled. From the results by using several grid size, it was also found that the optimal value of discharge should not be more than 0,03m3/second ;if the discharge exceeds the threshold value, it will make an impact on the position of the interface.From the simulation with grid dimension 6-by-6,the critical point of interface was formed when the pumping discharge rate was 0,036m3/second. Keywords : seawater intrusion, sharp interface, finite difference method
PENDAHULUAN Air tanah digunakan untuk hampir setengah bagian dari keseluruhan konsumsi air manusia. Pada daerah pantai, kelebihan pemompaan air akan menyebabkan penyusupan air asin ke dalam akuifer air tawar. Masuknya air asin ke dalam akuifer air tawar disebut sebagai 1
intrusi air laut. Kondisi ini menyebabkan penurunan kualitas air tanah. Kontrol terhadap beberapa proses intrusi telah dilakukan akan tetapi hanya mampu memberikan gambaran mengenai kondisi yang terjadi saat itu saja. Oleh sebab itu, simulasi numerik menjadi sangat penting untuk membuat prediksi pengaruh apa saja yang akan ditimbulkan pada masa yang akan datang. Dalam dunia pemodelan intrusi air laut, ada dua macam pendekatan yang dapat digunakan (Hinkelmann et al.,1999). Kondisi pertama mengasumsikan bahwa air tawar dan air asin merupakan dua buah fluida yang saling terpisahkan oleh bidang batas tegas. Hal ini menunjukkan bahwa kedua fluida tidak mengalami pencampuran. Sementara kondisi kedua mengasumsikan bahwa air asin dan air tawar mengalami pencampuran. Pencampuran ini terjadi pada zona transisi dimana gradien salinitas berubah terhadap dispersi hidrodinamik. Kondisi pertama merupakan bentuk penyederhanaan yang cukup masuk akal ketika zona transisi sempit sehingga dapat diabaikan. Namun apabila lebar zona transisi besar dan tidak dapat diabaikan, harus digunakan pendekatan kondisi kedua. Dalam makalah ini hanya akan dijelaskan mengenai pendekatan kondisi yang pertama. Untuk menggambarkan bidang batas yang tegas dalam simulasi intrusi air laut digunakan solusi numerik. Metode beda hingga digunakan untuk menyelesaikan persamaan aliran dengan melakukan diskretisasi menggunakan grid untuk mencari nilai head yang belum diketahui. Untuk permasalahan transisi seperti kasus intrusi air laut, dibutuhkan kondisi awal dan kondisi batas. Ada tiga kondisi batas yang dapat digunakan sebagai syarat batas domain model, yaitu: Tipe Dirichlet, tipe Neumann, dan tipe Cauchy (Notodarmojo,2005). Dalam pemodelan air tanah akan digunakan tipe Neumann sebagai kondisi batas.
METODOLOGI Metode penelitian yang diterapkan untuk studi mengenai simulasi intrusi air laut ditunjukkan dalam skema pada Gambar 1. Penentuan Objektif Model
Pengumpulan Data
Konseptualisasi
Pemrograman
Pemilihan Model
Penentuan Batas dan Domain Model
Verifikasi
Analisa Sensitivitas
Kalibrasi dan Uji Kemampuan Prediktif Model
Gambar 1. Skema metodologi penelitian
2
Penentuan Objektif Model Ada beberapa objek yang akan menjadi fokus dalam simulasi intrusi air laut, yaitu bentuk interface antara air tawar dengan air laut, pengaruh pumping terhadap tinggi tekan, serta ukuran grid untuk diskritisasi domain model. Pengumpulan Data Data yang digunakan untuk simulasi merupakan data sekunder yang didapat dari penelitian terdahulu (Hadi et al.,2003) untuk studi kasus bagian tengah Pulau Pari yang memiliki luas ±1km2 dan merupakan bagian dari Kepulauan Seibu yang terletak di sebelah utara Kota Jakarta. Konseptualisasi Formulasi model didahului dengan konseptualisasi daerah yang akan dimodelkan. Berikut adalah model konseptual aliran air dalam penampang melintang dua dimensi (cross section) dari bagian tengah Pulau Pari yang dapat dilihat pada Gambar 2.
Gambar 2. Pulau Pari (kiri) dan penampang melintang bagian yang dimodelkan (kanan) Penentuan Batas Domain Model Sebagai penyederhanaan dari kondisi yang sebenernya, pemodelan dilakukan dengan asumsi bahwa bidang kontak atau interface antara air tawar dan air asin merupakan batas tegas (sharp interface), berlaku aturan aproksimasi Dupuit-Forchheimer, pergerakan air laut diabaikan (pendekatan ini diterapkan pada keadaan dimana ketebalan air tawar sangat kecil bila dibandingkan dengan ketebalan air laut), arah aliran air horizontal, gradien hidrolik sama dengan muka air tanah, aliran terdistribusi secara seragam, akifer homogen dan isotropic, terdapat keseimbangan hidrostatis antara air tawar dan air asin, besaran permeabilitas dan storativitas bersifat konstan sepanjang arah vertical, sistem aliran secara keseluruhan dalam keadaan tunak (steady state), tipe kondisi batas adalah tipe Neumann dimana besar flux ditentukan. Batas domain model dapat dilihat pada Gambar 3.
3
Gambar 3. Batas Domain Model Pemilihan Model Pemodelan dilakukan dengan menggunakan metode numeric beda hingga. Persamaan aliran air tanah untuk dua fluida merupakan pengembangan dari konsep kekekalan massa dan Hukum Darcy. Berdasarkan penelitian yang terdahulu, aliran fluida dengan densitas yang berbeda dapat menghasilkan suatu larutan baru yang tercampur secara spontan atau dapat juga tetap dalam bentuknya masing-masing dan tidak tercampur (Tran,2004). Dari penelitian yang terdahulu (Hadi et al.,2003), persamaan-persamaan aliran air tanah untuk dua buah fluida dengan asumsi bahwa kedua fluida tidak saling bercampur diturunkan menjadi persamaan kontinuitas untuk sistem aliran air tanah yang dinyatakan sebagai berikut: ∁
(1)
Dimana adalah storativitas spesifik air tawar [1/L], adalah ketebalan akuifer air tawar [L], adalah tinggi tekan air tawar [L], t adalah waktu [T], adalah porositas efektif, adalah porositas kinematik (pada zona fluktuasi bidang batas), ∁ adalah posisi kedalaman bidang batas, adalah permeabilitas intrinsic air tawar [L/T], adalah berat jenis air tawar [M L-2 T-2], adalah viskositas air tawar, adalah discharge per satuan luas [L/T], adalah recharge per satuan luas [L/T], dan x,y menunjukkan arah horizontal dan vertikal. Letak notasi untuk persamaan (1) dilustrasikan dalam Gambar 4a.
4
Gambar 4a. Letak notasi untuk persamaan (1) Persamaan yang digunakan untuk membuat matriks dengan metode beda hingga adalah sebagai berikut: ,
,
,
,
(2)
,
Dengan penamaan grid yang dapat dilihat pada Gambar 4b. Aturan penamaan grid untuk dimensi yang lebih besar sama dengan yang terlihat pada Gambar 4b. Yang perlu dilakukan hanya memperbesar nomor pada arah x dan arah y.
Gambar 4b. Nomenclatur grid yang digunakan dalam diskretisasi Verifikasi Solusi numeric merupakan suatu hasil pendekatan sehingga tidak lepas dari akurasi dan stabilitas. Sebenarnya sangat sulit untuk mengetahui atau menguji akurasi solusi numeric karena solusi numeric umumnya digunakan untuk memecahkan persamaan diferensial yang tidak memiliki solusi analitik sebagai solusi eksaknya. Oleh sebab itu diperlukan dua persyaratan yang apabila dipenuhi akan memberikan semacam keandalan dari solusi tersebut, yaitu (Notodarmojo, 2005): Konvergensi error yang terjadi tidak boleh bertambah besar seiring berjalannya waktu
Stabilitas
, ,
∆ ∆
∆ ∆
; ∆
0 dan
,
∆
∆
Persyaratan di atas berlaku untuk semua node dan hanya terbatas untuk solusi eksplisit. Untuk solusi implicit tidak ada syarat stabilitas yang membatasi penggunaan ∆t. Analisis Sensitivitas Dalam simulasi intrusi air laut, komponen yang paling utama adalah besar grid yang digunakan, debit pemompaan, besar recharge, dan nilai permeabilitas intrinsik. Namun dalam makalah ini hanya akan dibahas mengenai besar grid dan debit pemompaan.
5
HASIL DAN PEMBAHASAN Perbandingan Left Hand Side Matriks Untuk Grid 5x5 dan 6x6 Mengacu pada persamaan (2) yang digunakan untuk membentuk matriks sebelah kiri dengan memasukkan batas domain tipe kedua (tipe Neumann), terbentuklah matriks sebelah kiri seperti yang terlihat pada Gambar 5. Perhatikan kotak berwarna kuning, hijau, dan merah. Kumpulan angka-angka tersebut membentuk matriks simetris yang disebut dengan matriks diagonal. Terlihat bahwa pola yang terbentuk baik untuk grid 5x5 maupun 6x6 adalah sama. Perbedaannya hanya terletak pada jumlah perulangan. Perulangan untuk grid 6x6 satu kali lebih banyak dari perulangan untuk grid 5x5. Semakin besar dimensi yang dibuat, semakin banyak perulangan yang dibutuhkan, semakin besar pula ukuran matriks yang dihasilkan. Besarnya ukuran matriks yang terbentuk dapat merepresentasikan dua hal. Pertama adalah perluasan domain model dimana luas area yang dimodelkan menjadi lebih besar dari semula. Kedua adalah diskretisasi yang semakin kecil dimana luas area yang dimodelkan tetap, tetapi segmen untuk tiap grid diperkecil sehingga jumlah segmen yang terbentuk menjadi semakin banyak.
3
-1
0
-1
0
0
0
0
0
3
-1
0
0
-1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-1
3
-1
0
-1
0
0
0
0
-1
3
-1
0
0
-1
0
0
0
0
0
0
0
-1
2
0
0
-1
0
0
0
0
-1
3
-1
0
0
-1
0
0
0
0
0
0
0
0
0
-1
0
0
4
-1
0
-1
0
0
0
0
-1
2
0
0
0
-1
0
0
0
0
0
0
0
0
0
-1
0
-1
4
-1
0
-1
0
-1
0
0
0
4
-1
0
0
-1
0
0
0
0
0
0
0
0
0
-1
0
-1
3
0
0
-1
0
-1
0
0
-1
4
-1
0
0
-1
0
0
0
0
0
0
0
0
0
-1
0
0
4
-1
0
0
0
-1
0
0
-1
4
-1
0
0
-1
0
0
0
0
0
0
0
0
0
-1
0
-1
4
-1
0
0
0
-1
0
0
-1
3
0
0
0
-1
0
0
0
0
-1
0
0
0
0
0
0
0
0
-1
0
-1
3
0
0
0
0
-1
0
0
0
4
-1
0
0
0
0
0
0
0
-1
0
0
-1
4
-1
0
0
-1
0
0
0
-1
0
0
0
0
0
0
0
-1
0
0
-1
4
-1
0
0
0
0
0
0
0
0
-1
0
0
-1
3
0
0
0
-1
-1
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
4
0
0
0
0
0
0
0
0
0
-1
0
0
-1
4
-1
0
-1
4
-1
0
-1
3
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
0
0
0
0
0
0
0
0
0
0
-1
0
Gambar 5. Left hand side matrix untuk grid dengan dimensi 5x5 (kiri) dan 6x6 (kanan) Untuk kasus aliran dalam kondisi tunak, besarnya diskretisasi tidak berpengaruh terhadap hasil yang didapatkan karena ∆x=∆y=0, sehingga ukuran grid hanya akan berpengaruh terhadap luas area yang dimodelkan. Pengaruh Besar Pemompaan Terhadap Head Air Tawar Di Atas Muka Air Laut Untuk melihat pengaruh pemompaan terhadap tinggi head, digunakan input debit sebesar 0 m3/detik, 0,02m3/detik, 0,03m3/detik, dan 0,04m3/detik. Keempat nilai parameter debit yang digunakan disimulasikan untuk model dengan dimensi 6x6. Berikut adalah hasil head yang didapatkan untuk setiap titik pada grid (lihat Gambar 6):
6
A
B
C
D
Gambar 6. Nilai head yang terbentuk untuk Q= 0 m3/detik (A), Q=0,02m3/detik (B), 0,03m3/detik (C), dan 0,04m3/detik (D) Kondisi (A) menunjukkan tinggi head pada saat tidak terjadi pemompaan. Batas 0m merupakan tinggi muka air laut yang menjadi datum dalam perhitungan. Apabila tidak terjadi pemompaan, head tertinggi berada pada ketinggian 250m di atas muka air laut. Ketika mulai terjadi pemompaan, head air tawar mengalami penurunan hampir sebesar 170m. Head tertinggi pada kondisi (B) berada pada ketinggian 120m di atas muka air laut. Ketinggian head terus mengalami penurunan (lihat kondisi (C)) hingga akhirnya terdapat nilai head yang negative seperti terlihat pada kondisi (D). Hal ini menunjukkan bahwa ketersediaan air dalam akifer air tawar tidak mencukupi sehingga terjadi penyusupan air asin ke dalam akifer air tawar untuk memenuhi kekosongan akifer air tawar akibat pemompaan. Artinya agar dalam kasus ini tidak terjadi intrusi air laut, debit pemompaan maksimal adalah sebesar 0,03m3/detik. Nilai batas pemompaan ini hanya berlaku untuk simulasi model dengan grid 6x6. Nilai batas pemompaan akan menjadi berbeda untuk luas area yang berbeda pula. Bentuk Interface Pada Titik Kritis Sebelumnya telah dibahas mengenai debit maksimal yang diperbolehkan untuk pemompaan. Sekarang akan dicari nilai perkiraan debit pemompaan yang menghasilkan kondisi stabil pada bidang batas. Rentang debit yang digunakan antara 0,03m3/detik – 0,04m3/detik dengan kenaikan sebesar 0,002m3/detik untuk setiap simulasi. Berikut adalah hasil yang didapatkan (lihat Gambar 7). Perhatikan bulatan berwarna biru dan hijau yang terdapat pada grafik (A) hingga grafik (F). Bulatan biru merupakan kondisi head batas 16m yang terus mengalami penurunan sebesar 4m. Fokus pada bulatan hijau yang terletak di bagian paling kiri dari setiap grafik. Pada kondisi (A), head bulatan hijau terletak pada ketinggian sekitar 11m. Pada kondisi B mengalami penurunan hingga ketinggian sekitar 8m. Pada kondisi (C), head menjadi semakin menurun yaitu setinggi 5m. Pada kondisi (D) ketinggian head adalah 7
sebesar 2m yang menunjukkan bahwa ketersediaan air dalam akifer air tawar dalam keadaan kritis. Adanya sedikit saja kenaikan debit pemompaan akan menyebabkan intrusi air laut. Hal ini ditunjukkan dengan nilai ketinggian head akibat pemompaan dengan head batas yang ditentukan yang memiliki nilai yang hampir sama (Perhatikan bulatan hijau keempat dan kelima pada kondisi(D)). A
B
C
D
E
FF
Gambar 7. Kondisi interface yang terbentuk untuk Q= 0,03 m3/detik(A), Q=0,032m3/detik(B), 0,034m3/detik(C), 0,036m3/detik(D), 0,038m3/detik(E), dan 0,04m3/detik(F) Sementara untuk kondisi (E) dan (F), ketinggian head sudah berada dibawah muka air laut. Dari keseluruhan grafik (A) – (F), dapat ditarik kesimpulan bahwa batas kritis interface terbentuk pada saat debit pemompaan (Q) = 0,036m3/detik.
8
KESIMPULAN Matriks diagonal pada left hand side memiliki pola yang sama untuk setiap dimensi grid yang dibuat. Perbedaannya hanya terdapat pada jumlah perulangan. Semakin besar dimensi grid yang dibuat, semakin besar perulangan yang dibutuhkan, dan semakin luas area studi yang dapat dimodelkan. Pada simulasi ini diketahui bahwa semakin besar debit pemompaan, semakin kecil tinggi head air tawar yang dihasilkan yang menunjukkan bahwa ketersediaan air dalam akifer air tawar semakin sedikit. Kekosongan akibat pemompaan menyebabkan air asin masuk ke dalam akuifer dan mengisi kekosongan tersebut. Pada akhirnya yang terambil dari pemompaan bukan air tawar melainkan air asin. Debit pengambilan maksimal yang diperbolehkan untuk pemodelan dengan grid 6x6 adalah sebesar 0,03m3/detik. Dari hasil simulasi untuk grid dengan dimensi 6x6, didapatkan bahwa posisi kritis bidang batas terbentuk pada saat debit pemompaan sebesar 0,036m3/detik. Hal ini akan berbeda apabila luas area yang dimodelkan berbeda pula.
DAFTAR PUSTAKA Hadi, S.,Kusuma, S.K., Bakti, H., Subardja, A.,(2003).Hidrodinamika Air Tanah Pulau Kacil Terumbu Karang: Kajian Berdasarkan Model Bidang Batas Tegas dengan Pendekatan Metoda Beda Hingga.RISET Geologi dan Pertambangan,Jilid 13 No.1.51-64. Herlambang, A., Indriatmoko, R.H.,(2005).Pengelolaan Air Tanah Dan Intrusi Air Laut.JAI Vol.1,No.2.214-217. Hinkelmann, R., Sheta, H., Class, Helmig,R.,(1999).A Comparison of Different Model Concepts For Saltwater Intrusion Processes.IAHS Publ. no. 263, 2000. Notodarmojo, Suprihanto.,( 2005). Pencemaran Tanah dan Air Tanah. Tran, T.M.,(2004).Multi Objective Management of Salwater Intrusion in Groundwater.Thesis of Environmental Engineering. Zhang, Y.,(2011).Groundwater Flow and Solute Transport Modeling.Dept. of Geology and Geophysics.University of Wyoming.
9
