## Aplikasi penggunaan Metode Finite Difference dalam pemodelan air tanah pada kasus pertambangan : literatur review

## Application of the Finite Difference Method in groundwater modeling in mining cases: literature review
Hana Trijayanti*, Tedy Agung Cahyadi, Rika Ernawati Magister Teknik Pertambangan Universitas Pembangunan Nasional “Veteran” Yogyakarta Jalan Padjajaran (SWK) No. 104 Sleman 55283 Indonesia e-mail: * trihananas@gmail.com
ABSTRAK Pada  area  kegiatan  penambangan  terdapat  berbagai  macam  kasus  yang  berkaitan  dengan  hidrogeologi. Permasalahan tersebut dapat diselesaikan dengan metode perhitungan numerik, yaitu finite difference . Finite difference digunakan kedalam Visual Modflow dengan membangun suatu model hidrogeologi. Dalam pemodelan Visual Modflow, tahapan yang dilakukan yaitu  persiapan  model  konseptual,  menentapkan  batas  area  pemodelan,  kalibrasi  dan  validasi. Pemodelan dapat diaplikasikan untuk menyelesaikan kasus air tanah pada pertambangan yaitu pemodelan air tanah setelah dilakukan penambangan, pemasangan drain hole pada dinding pit , pemodelan air tanah untuk kegiatan dewatering, intrusi air laut akibat adanya penambangan dan pemodelan sebaran kontaminan akibat kolam tailing penambangan. Pemodelan tersebut  dapat  menjadi  acuan  dalam  monitoring,  pengendalian,  dan pemanfaatan  terkait  dengan  pengelolaan  dan perlindungan air tanah.
Kata-kata kunci : air tanah, hidrogeologi, model, modflow, penambangan
ABSTRACT In the mining activity area there are various cases related to hydrogeology. These problems can be solved by numerical calculation method, namely finite difference. Finite difference is used into Visual  Modflow by building a hydrogeological model. In Visual Modflow modeling, the steps taken are the preparation of the conceptual model, setting the boundaries of the modeling area, calibration and validation. Modeling can be applied to solve groundwater cases in mining, namely groundwater modeling after mining, installation of drain holes on pit walls, groundwater modeling for dewatering activities, seawater intrusion due to mining and contaminant distribution modeling due to mining tailings ponds. The modeling can be used as a reference in monitoring, controlling, and utilization related to groundwater management and protection.
Keywords : groundwater, hydrogeology, model, modflow, mining

## PENDAHULUAN
Permasalahan pada air tanah seringkali dihadapkan dengan hal yang kompleks. Letaknya yang berada di bawah permukaan banyak dipengaruhi oleh berbagai faktor secara fisika maupun kimia  yang dapat mengubah  karakteristik dari  air tanah.  Dalam pengembangannya  sifat dan karakteristik  dari  air  tanah dapat  diketahui  dengan  menggunakan  berbagai  macam model  perhitungan,  salah  satunya yaitu dengan model numerik. Model numerik dapat dihitung berdasarkan Metode Perbedaan  Hingga (Finite Difference), dengan menetapkan nilai parameter hidrogeologi menjadi titik (node) dalam suatu domain. Finite Difference merupakan penjabaran dari persamaan Laplace dengan didasari konsep blok tiga dimensi (Anderson dkk, 2015).
Finite Difference membutuhkan penyederhaan data lapangan yang tinggi, maka diperlukan suatu simulasi model yang dibuat kedalam bentuk yang efisien dengan perangkat lunak khusus, yaitu Visual Modflow. Pemodelan dapat menjadi gambaran mengenai kondisi yang sebenarnya terjadi di lapangan  dan  dapat  memberikan
pemahaman  tentang  karakteristik  air  tanah.  Visual Modflow  untuk  pertama kali diperkenalkan oleh Waterloo Hydrogeologic pada bulan Agustus 1994, dengan mengembangkan konsep model air tanah yang ditemukan oleh Darcy ( Hariharan, 2017 ). Komponen yang diamati pada  konsep  air  tanah  yaitu  adanya  perbedaan head hidrolik yang  digunakan  untuk  mencari  arah  dan kecepatan aliran berdasarkan Hukum Darcy (Masud, 2002). Berdasarkan Hukum Darcy, konduktivitas hidrolik batuan merupakan sifat penting dalam pergerakan air tanah.
Pengembangan pemodelan air tanah telat banyak digunakan  pada  kasus  yang  terjadi  pada  pertambangan. Aplikasi pemodelan sangat berguna untuk mengevaluasi dampak kegiatan eksploitasi  hingga  penutupan  tambang terhadap air permukaan dari sudut pandang kuantitatif dan kualitatif (Rapantova, 2007). Area tambang yang nampak tergenang dapat dilakukan investigasi mendalam  melalui simulasi  pemodelan  berbentuk  perubahan  pola  aliran  air tanah untuk mengetahui asal mula dan pergerakan aliran air, hingga pembuatan simulasi pergerakan kontaminan yang berasal dari air asam tambang.

## METODOLOGI
Tahapan  penelitian  dengan  metode  literatur  review. Literatur  review  merupakan  proses  pengumpulan informasi yang relevan dengan topik atau masalah yang menjadi  objek penelitian.  Informasi  tersebut  diperoleh dari beberapa jurnal, buku, ataupun tulisan yang berkaitan dengan  penelitian  dari  studi  pustaka yang dilakukan pengkajian.

## HASIL DAN DISKUSI
Pemodelan Visual Modflow didasarkan pada metode perhitungan beda hingga (Finite Difference).  Dalam metode Finite Difference disimulasikan dengan indek i, j,
k yang mewakili kolom, baris, lapisan, dan simbol node (titik) didalam ruang 3D (Gambar 1.). Jarak node di sepanjang baris  ditentukan  oleh ∆ x,  spasi sepanjang kolom  ditentukan oleh  ∆y, dan  spasi antar  lapisan ditentukan oleh ∆ z. Pada  gambar (a), ditampilkan  dalam bentuk 2D, node  terletak di  pusat  dalam sel atau blok Finite Difference , yang merupakan simbol yang mewakili head di blok tersebut. i merupakan kolom  dan j merupakan baris. Garis hitam yang melingkar mewakili domain masalah. Sel atau blok yang  berada diluar garis batas domain merupakan sel yang tidak aktif. Pada gambar (b), ditampilkan dalam  bentuk 2D, membentuk kelompok lima node dengan pusat node berada ditengah (i,j). Pada gambar (c),  ditampilkan  dalam  bentuk  3D,  dimana  ∆z mewakili jarak vertikal antar node, dan k adalah indeks vertikal. Sehingga dapat terlihat pusat node dalam 3D yaitu berada di tengah (h i,j,k ), akan dipengaruhi oleh faktor dari 7 node di sekelilingnya.
Bentuk persamaan Finite Difference dituliskan dengan Persamaan 1. dari persamaan diferensial  yang mewakili aliran air tanah 3D untuk kondisi heterogen dan anisotrop : ( Anderson dkk, 2015 ).
𝜕𝜕 𝜕𝜕𝑥𝑥 �𝐾𝐾𝑥𝑥 𝜕𝜕ℎ 𝜕𝜕𝑥𝑥 � + 𝜕𝜕 𝜕𝜕𝑦𝑦 �𝐾𝐾𝑦𝑦 𝜕𝜕ℎ 𝜕𝜕𝑦𝑦 � + 𝜕𝜕 𝜕𝜕𝑧𝑧 �𝐾𝐾𝑧𝑧 𝜕𝜕ℎ 𝜕𝜕𝑧𝑧 � ± 𝑊𝑊 = 𝑆𝑆𝑆𝑆 𝜕𝜕ℎ 𝜕𝜕𝑡𝑡 Dimana: Kx, Ky, Kz
= Konduktivitas hidrolik pada arah x,y,z (LT -1 ) h = Tinggi muka air tanah (L) W = Flux dengan menghitung faktor penambah, contoh : pemompaan atau penambahan dari sumber lain (L 3 T -1 ) Ss = Specific storage (LT -1 ) x, y, z = Arah koordinat t = Waktu (T)
Tahapan pemodelan Visual Modflow memerlukan persiapan awal, yang terdiri dari pembuatan model konseptual, penentuan kondisi batas model, kalibrasi dan validasi model.
A. Model Konseptual Model  konseptual  hidrogeologi  adalah  representasi  dari sistem  air  tanah  yang  digabungkan  dengan  interpretasi kondisi geologi, hidrologi, serta keseimbangan air ( water balanc e).  Model  tersebut  dapat  memberikan  info karakteristik bawah permukaan dan membantu mengurangi ketidakpastian yang mengatur aliran air tanah dan transport kontaminan di daerah penelitian ( Torres,  2020 ).  Contoh dari hasil pembuatan konseptual model dapat dilihat pada Gambar 2.
Gambar-1. Grid dan Notasi Beda Hingga ( Finite Difference ) ( Anderson dkk, 2015 )
Dalam konseptual model dapat terlihat berupa bangun tiga dimensi  yang  membentuk suatu  komponen ‘wadah’ atau perlapisan litologi yang ada pada daerah tersebut, sehingga hal inidapat menjadi salah satu acuan dalam mengasumsi kondisi hidrogeologi di daerah tersebut.
B. Kondisi Batas Model Kondisi batas simulasi numerik digunakan untuk menggambarkan kondisi fisik batas rembesan, yaitu kondisi muka air tanah yang dianggap konstan dari daerah rembesan. Infiltrasi curah hujan adalah sumber yang utama di wilayah penelitian (Xie dkk, 2020). Kondisi batas didapat berdasarkan  kondisi alam yang mendukung pada wilayah penelitian, seperti daerah pegunungan, batas lapisan  stratigrafi,  kontur,  danau  dan  sungai,  sebab  lebih mudah untuk diidentifikasi dan dikonsepkan.
C. Kalibrasi dan Validasi Model Hasil  perhitungan  numerik  dalam  software  memiliki perbedaan  jika  dibandingkan  dengan  kondisi yang sebenarnya di lapangan, sehingga akan ditemukan beberapa kesalahn pada saat proses  pembuatan model. Maka,
dilakukan kalibrasi model agar sesuai dengan kondisi di lapangan. Selama  proses kalibrasi dalam pemodelan ini, nilai konduktivitas hidrolik atau groundwater recharge akan  disesuaikan dalam model sesuai dengan tinggi air tanah simulasi dan air tanah terukur. Kalibrasi ini  akan membandingkan hydraulics  head hasil  perhitungan simulasi Modflow dengan Metode Root Mean Square Error (RMSE) yang  telah  tersedia  di software. Model  yang dikalibrasi harus memberikan hasil yaitu nilai RMSE yang minimal dan memiliki kecocokan antara head yang terukur dilapangan dengan yang dimodelkan (Tahershamsi, 2018). Nilai RMSE dapat dilihat pada grafik  scatter, seperti Gambar 3. Jika telah terkalibrasi, maka model telah tervalidasi dan siap untuk tahapan  selanjutnya  yaitu prediksi pola aliran air tanah.
Setelah  model  konseptual  telah  terkalibrasi  dan tervalidasi,  maka  selanjutnya  model  dapat  diaplikasikan untuk berbagai tujuan pengembangan research , seperti identifikasi pola aliran air tanah,  pembuatan drain hole , hingga simulasi penyebaran kontaminan.  Berikut merupakan contoh hasil  penelitian  dengan  menggunakan Metode Finite Different pada area pertambangan.
Gambar-3. Diagram Scatter Setelah Terkalibrasi ( Devy, 2018 )
Gambar-4 . Pola Aliran Airtanah dan Zona Neraca Airtanah Penambangan Aktif ( Devy, 2018 )
a. Pemodelan  Air  Tanah  Setelah  Terjadinya Penambangan Batubara Openpit
Tata guna lahan akan berpengaruh terhadap ketersedian air tanah secara kuantitas maupun kualitas. Perubahan pola aliran air tanah pada daerah penelitian dipengaruhi oleh faktor stuktur geologi dan kegiatan penambangan open pit . Setelah dilakukan kegiatan penambangan batubara open pit , terjadi  perubahan  morfologi  dan  geologi  sehingga mempengaruhi  aliran  air  tanah.  Setelah  dilakukan pemodelan  menggunakan  Visual  Modflow  didapatkan bahwa terjadi variasi terhadap nilai head hidrolik , dimana terjadi perubahan pola aliran air tanah menuju pit penambangan, muka air tanah  menurun, dan kenaikan neraca air tanah yang disebabkan perubahan tata guna lahan ( Devy, 2018 ) (Gambar 4).
b. Pemasangan Drain Hole pada Dinding Pit
Dinding pit penambangan selalu dihadapkan dengan kasus kelongsoran, terlebih pada kondisi dinding yang jenuh air dan terjadi rembesan. Air yang berada pada dinding tersebut dapat diminimalisir dengan pemasangan drain hole. Drain hole merupakan salah satu upaya dalam menurunkan muka airtanah dengan memanfaatkan gaya gravitasi. Dalam pemasangan drain hole pada dinding pit terdapat beberapa perkiraan yang dilakukan dengan suatu perancangan, agar fungsi dari drain hole dapat maksimal dan pemasangan yang seefisien mungkin.  Perencanaan desain dran hole dapat menggunakan perhitungan finite difference dengan model Visual Modflow. Dilakukan simulasi model terhadap nilai konduktivitas hidrolik, recharge , discharge , diameter pipa, dan panjang pipa. Maka, didapatkan hasil penurunan muka air tanah yang dihasilkan karena pemasangan drain hole dan  dapat menentukan jumlah dan jarak dari setiap drain hole yang akan dipasang pada dinding pit tambang (Alvando, 2017). Jika jarak antar drain hole semakin rapat, maka akan terjadi peningkatan penurunan muka air tanah. Namun disisi lain, semakin rapat jarak antar drain hole akan menurunkan  kekuatan lereng sebab daya dukung lereng tersebut akan semakin hialang akibat pemasangan drain hole tersebut (Fadhila, 2020).
c. Pemodelan  Air  Tanah  untuk Kegiatan  Dewatering Tambang
Kegiatan penambangan memerlukan kondisi area kerja yang tidak memiliki banyak genangan air.  Untuk tetap berada pada kondisi area tambang yang kering, maka diperlukan suatu upaya untuk mengeluarkan air yang telah terlanjur masuk kedalam pit tambang, dengan cara mine dewatering. Mine dewatering tidak terlepas dari fungsi pemompaan, sehingga diperlukan suatu perhitungan yaitu debit pompa dan dilakukan perhitungan lamanya waktu yang dibutuhkan untuk memompa (Adnyano, 2020). Kegiatan  dewatering  memiliki  dampak  yang  signifikan terhadap  keseluruhan  biaya  produksi  tambang. Maka, rencana kegiatan dewatering yang efektif dapat mengurangi biaya produksi dengan  jumlah yang cukup besar. Pemodelan aliran air tanah numerik dapat digunakan untuk membantu  dalam  membuat  rencana  desain  dewatering dengan  mensimulasikan  muka  air  tanah  dan  aliran rembesan. Hal ini bertujuan untuk mengidentifikasi lokasi dan jadwal pemompaan untuk dewatering sumur dan untuk menilai kinerja metode dewatering sebelumnya. Rencana dewatering dapat dirancang menggunakan Visual Modflow  dengan  memperhatikan  kondisi drainase dan fluktuasi muka air tanah pada sumur. Dengan
memperhatikan fluktuasi muka air tanah, kemudian dapat dilakukan suatu prediksi untuk waktu dan jumlah pompa pada saat proses  dewatering. Muka air tanah yang tinggi dapat memberi asumsi bahwa jumlah air yang akan masuk kedalam pit akan meningkat, begitupun sebaliknya. Maka, dibutuhkan rencana dengan dasar perhitungan dari Visual Modflow untuk dapat memperkirakan kebutuhan pompa yang akan digunakan (Alloisio, 2014).

## d. Intrusi Air Laut Akibat Penambangan
Ketersediaan air tanah di daerah pesisir memiliki kuantitas dan kualitas yang cenderung minim.  Aktivitas produksi sumur air tanah hingga adanya kegiatan penambangan di daerah pesisir dapat berpengaruh terhadap gerakan air laut yang semakin mendesak menuju daratan, yang disebut dengan intrusi air laut. Aktivitas penambangan terkhusus pada proses penggalian pit yang semakin dalam dan semakin dekat permukaan laut dan didukung dengan pemanfaatan sumur warga yang tinggi pada daerah pesisir, dapat berpotensi membuka jalur air laut akan masuk ke dalam sumur produksi yang digunakan warga sekitar untuk keperluan sehari (Gambar 5).Untuk mengatasi kebutuhan air tanah,  dapat dilakukan suatu simulasi intrusi air laut dengan menggunakan Visual Modflow. Simulasi dibuat dengan beberapa pertimbangan yaitu keselarasan antara nilai recharge, push back penambangan, dan  jumlah pemanfaatan air sumur pada sumur warga didukung dengan uji kualitas air. Simulasi dapat dirancang untuk membuat suatu prediksi jangka panjang dari jarak pergerakan intrusi air laut yang berada di air tanah (Sajeena, 2020).

## e. Pemodelan  Sebaran  Kontaminan  Akibat  Kolam Tailing
Seiiring berjalannya waktu, kualitas air tanah akan semakin menurun oleh adanya suatu zat akibat  dari  kegiatan penambangan.  Setiap  komoditas  tambang  akan menghasilkan  sisa  ( tailing )  dengan karakteristik  yang bervariasi dan tergolong limbah. Tailing yang mengandung logam berat dikumpulkan dalam suatu kolam. Unsur logam berat pada kolam tailing dapat bereaksi secara fisika dan kimia, seperti difusi, dispersi, desorpsi dan disosiasi membentuk suatu kontaminan ( Zhang dkk,  2018 ). Kontaminan ini dapat bermigrasi keluar dari kolam penampung tailing menuju sungai maupun meresap ke air tanah dan dapat terbawa hingga puluhan kilometer ( Bortnikova dkk, 2013 ).  Hal ini  tidak ramah lingkungan apabila kontaminan dapat tembus dan terakumulasi ke dalam air tanah yang digunakan oleh masyarakat. Identifikasi mengenai karakteristik dan mekanisme pergerakan dari kontaminan yang ada pada air tanah akibat dari kolam tailing tersebut dapat diterapkan dengan metode numerik Beda  Hingga  membentuk suatu model yang dibantu dengan perangkat lunak Visual Modflow ( Xie dkk, 2020 ).  Model dapat dirancang untuk simulasi sebaran kontaminan pada jangka waktu kedepan sesaui dengan rencana  penelitian.  Simulasi  dibuat  dengan  beberapa pertimbangan dan asumsi terhadap model konseptual yang terdiri dari pengaruh karakteristik hidrostratigrafi, recharge , dan waktu produksi  tambang didukung dengan data uji kualitas pada air kolam tailing, air sungai, dan sedimen sungai  terdekat. Hasil dari simulasi dapat menunjukkan pola pergerakaan sebaran kontaminan yang dapat diketahui melalui waktu dan jarak yang dibutuhkan pada sebaran kontaminan tersebut ( Xie dkk, 2020 ) (Gambar 6 dan Gambar 7).
Gambar-5. Peta kontur muka air yang diprediksi setelah meningkatkan laju pemompaan ( Sajeena, 2020 )
Gambar-6. Prediksi Sebaran Mangan pada Air Tanah Dalam 20 Tahun (Tampak Horizontal) ( Xie dkk, 2020 )

## KESIMPULAN
Finite Difference merupakan metode dasar yang digunakan pada perangkat lunak Visual Modflow. Visual Modflow dapat merancang sebuah model hidrologi, dimana terdapat beberapa tahapan dalam  membangun model tersebut. Model konseptual adalah hal yang paling mendasar dalam rancangan model yang dilakukan dengan Visual Modflow. Dalam konseptual model terdiri atas lapisan titologi  yang  dapat  dikelompokkan  menjadi komponen  lapisan  unit  hidrostratigrafi  yang  memiliki karakteristik  yang  berupa  konduktivitas  hidrolik  batuan, porositas, dan struktur geologi yang mempengaruhi posisi lapisan batuan. Dibutuhkan pula batasan dalam pemodelan yang dapat berupa dataran tinggi atau sungai yang disebut dengan constant head . Selanjutnya, dilakukan kalibrasi yang bertujuan untuk menyesuaikan kondisi model dengan keadaan yang sebenarnya di lapangan. Hasil kalibrasi dapat dilihat pada persentase RMSE pada grafik scatter yang memiliki nilai rendah. Jika nilai RMSE rendah, maka model tersebut telah terkalibrasi dan dinyatakan valid untuk dilakukan berbagai simulasi model selanjutnya.
Modflow dapat  diaplikasikan  untuk  merancang simulasi  dan  prediksi  pada  kasus  hidrogeologi  yang terdapat di tambang, sehingga langkah tersebut dapat dijadikan acuan sebagai upaya terbaik  dalam  monitor, pengendalian,  serta  pemanfaatan  terkait  pengelolaan  dan perlindungan airtanah.

## UCAPAN TERIMA KASIH
Penulis  mengucapkan  terima  kasih  kepada Lembaga Penelitian dan Pengabdian Masyarakat (LPPM) Universitas Pembangunan Nasional “Veteran” Yogyakarta yang telah mendanai penelitian ini.

## DAFTAR PUSTAKA
[1]   Anderson, M. P., Woessner, W. W., & Hunt, R. J.
(2015). Applied Groundwater  Modeling: Simulation Of Flow And Advective Transport . Academic Press . London.
[2]   Devy, S. D. (2018). Pemodelan airtanah dan Neraca Airtanah Dampak Penambangan Batubara Open Pit pada Lipatan Sinklin di Daerah Muara Lawa, Kabupaten Kutai Barat, Provinsi Kalimantan Timur. SPECTA Journal of Technology , 2 (2), 69-83.
[3]   Alvando, P. V., Darul, A., & Irawan, D. E. (2017). Pemodelan Dampak Pemasangan Inclined Drain Hole  Pada Wall  Pit  Pt  Xxx  Terhadap Penurunan Muka Airtanah.
[4]   Fadhilah, R., Widodo, L. E., & Iskandar, I. (2020). Pengaruh Jarak Antar Drain Hole Terhadap Penurunan Muka Air Tanah pada Lereng Tambang Terbuka Batubara. Journal of Science and Applicative Technology , 4 (2), 116-120.
[5]   Alloisio, S. A. R. A. H., Douglas, B., McKittrick, R., & Prigneau, P. (2004). Groundwater modelling for
large-scale mine dewatering in Chile: MODFLOW or FEFLOW.
[6]   Adnyano, A. I. A., & Bagaskoro, M. (2020). Kajian Teknis Dewatering System Tambang Pada Pertambangan Batubara. PROMINE.
[7]   Sajeena, S., & Swathy, P. S. (2020). Studies on saline water intrusion in the coastal stretch of Kadalundi river basin, Malappuram district, Kerala using visual MODFLOW–A case study.
[8]   Hariharan, V., & Shankar, M. U. (2017). A review of visual MODFLOW applications in groundwater modelling. In IOP Conference Series: Materials Science and Engineering (Vol. 263, No. 3, p. 032025). IOP Publishing. India.
[9]   Rapantova, N., Grmela, A., Vojtek, D., Halir, J., & Michalek, B. (2007). Ground water  flow modelling applications in mining hydrogeology. Mine water and the environment , 26 (4), 264-270. Italy.
[10] Alloisio, S. A. R. A. H., Douglas, B., McKittrick, R., & Prigneau, P. (2004). Groundwater modelling for large-scale mine dewatering in Chile: MODFLOW or FEFLOW.
Water Management Consultants .Chile.
[11] Zhang, X., Yang, H., & Cui, Z. (2018). Evaluation And Analysis Of Soil Migration And Distribution Characteristics  Of  Heavy  Metals  In  Iron Tailing s. Journal Cleaner Production Vol. 172 . USA.
[12] Bortnikova, S., Yurkevich, N., Bessonova, E., Karin, Y., & Saeva, O. (2013) : The Combination Of Geoelectrical Measurements And Hydro-
Geochemical  Studies  For  The  Evaluation  Of Groundwater  Pollution  In  Mining Tailing s  Areas. Threats To The Quality Of Groundwater Resources . Springer, Berlin .
[13] Xie, W., Ren, B., Hursthouse, A. S., Wang, Z., & Luo,  X. (2020).  Simulation  Of  Manganese Transport In Groundwater Using Visual Modflow: A Case Study From Xiangtan Manganese Ore Area In Central  China .  Polish Journal Of  Environmental Studies Vol. 30 . Polandia.
[14] Torres, H. (2020). Assessing Groundwater Contamination Risk And Detection Of Unknown Sources  Using A Multi-Component Reactive Transport Model. Journal Of Geoscience And Environment Protection Vol.8. Australia.
[15] Masud,  A.,  &  Hughes,  T.  J.  (2002).  A  stabilized mixed  finite  element  method  for  Darcy flow. Computer  methods  in  applied mechanics  and engineering , 191 (39-40), 4341-4370. USA.
[16] Tahershamsi, A., Feizi, A., & Molaei, S. (2018). Modeling groundwater surface by MODFLOW math code  and geostatistical  method. Civil  Engineering Journal , 4 (4), 812. Iran.