21
3 METODE PENELITIAN
3.1 Waktu dan Lokasi Penelitian Penelitian ini dalam dua tahapan, yaitu pengambilan data lapangan dan simulasi pemodelan. Pengambilan data lapangan dilaksanakan pada bulan Januari 2007 – Agustus 2007 dengan tahapan pengambilan contoh yang dilakukan pada bulan Januari- Februari, April – Mei dan Juli- Agustus. Pengambilan sampel dilakukan setiap dua minggu dan dilanjutkan dengan analisis contoh. Analisis contoh dilakukan di Laboratorium Kimia FMIPA Unila. Simulasi model hidrodinamika dan model ekosistem di perairan Teluk Lampung dilakukan pada posisi 5o24'2" - 5o46'26" LS dan 105o8'7"-105o37'12' BT (Gambar 4.). Simulasi pemodelan dilakukan pada bulan Juli 2009 – Desember 2010 yang meliputi pengumpulan dan analisis data sekunder, penyusunan model numerik dan parameter pemodelan, simulasi dan analisis hasil simulasi. Waktu simulasi disamakan dengan waktu pengambilan data lapangan. 3.2 Penentuan Lokasi Pengambilan Contoh Penentuan lokasi pengambilan contoh ditentukan dengan pertimbangan bahwa lokasi (stasiun) tersebut merupakan daerah sumber nutrien bagi perairan. Lokasi pengambilan contoh ditabulasikan dalam Tabel 5 dan disajikan pada Gambar 4. Tabel 5 Lokasi pengambilan contoh Stasiun 1 2 3 4 5 6 7 8 9 10
Posisi geografis Bujur Timur Lintang Selatan 105o16’12” 05o29’14’ o 105 15’15” 05o27’34” o 105 18’20” 05o27’12” o 105 15’8” 05o28’14” o 105 19’26” 05o29’23” o 105 21’46” 5o32’58” o 105 14’58” 5o32’27” o 105 15’02” 5o34’25” o 105 12’13” 5o40’00” o 105 20’00” 5o40’00”
Keterangan Pelabuhan Pendaratan Ikan Muara sungai Way Kahuripan Muara Sungai Way Lunik Pelabuhan Peti Kemas Kawasan Permukiman dan Industri Lokasi Wisata Pasir Putih Kawasan Tambak Hanura Kawasan Tambak Sidodadi Lokasi keramba jaring apung P. Pahawang Tengah laut
21
22
Gambar 4 Peta lokasi daerah model dan lokasi pengambilan sampel. 3.3 Alat dan Bahan Alat dan bahan yang digunakan dalam penelitian ini ditabulasikan pada Tabel 6. Tabel 6 Alat dan bahan yang digunakan dalam penelitian Alat dan Bahan Perangkat survei lapangan GPS (Global Positioning System) Botol sampel -
CTD Botol Van Dorn Spektrofotometer
DO meter dan pH meter Pompa vakum Stempel Pippete 0,1 ml Larutan Formalin dan Lugol Oven, Desikator dan Timbangan digital Perangkat analisis data Perangkat lunak komputer (Surfer 9, ArcGis 9.3, MS Excel, ODV 4.3, ARMSLite 2.12, Intel Visual Fortran 11)
Kegunaan -
Penentuan posisi Tempat sampel air dan sampel plankton Mengukur parameter fisik perairan Pengambilan sampel air Mengukur nilai absorbansi sampel
-
Pengukuran DO dan pH Penyedot sampel air untuk analisis kimia Pengambil fraksi fitoplankton Pengawet sampel plankton Analisis partikel padatan tersuspensi
-
Simulasi, visualisasi dan analisis data
23
3.4 Pengambilan Contoh 3.4.1 Pengambilan Contoh Air Laut Pengambilan sampel air dilakukan pada kedalaman 0 – 2 m menggunakan botol Van Dorn PVC. Pengambilan sampel air dilakukan untuk analisa kandungan nutrien terlarut dan kandungan padatan tersuspensi. Selain pengambilan sampel air juga dilakukan pengukuran secara langsung meliputi parameter oksigen terlarut menggunakan DO meter, pH menggunakan pH meter yang dikalibrasi dengan larutan standar dengan pH 4.00 dan 9.00 sebelum digunakan untuk pengukuran. Temperatur dan salinitas diukur menggunakan STD meter model YSI-30. 3.4.2 Pengambilan Contoh Plankton Pengambilan contoh fitoplankton dilakukan menggunakan plankton net dengan ukuran mata jaring 40 µm, dan untuk zooplankton menggunakan plankton net dengan ukuran mata jaring 80 µm. Sampel ditempatkan dalam wadah cool box untuk dianalisa kandungan klorofil-a dalam fitoplankton dan kandungan karbon dalam zooplankton. 3.4.3 Pengambilan Data Pasang Surut Data fluktuasi muka laut diperoleh dari stasiun pengamatan di Pelabuhan Peti Kemas PT. PELINDO II Panjang, dengan interval satu jam mulai pukul 00.00 sampai pukul 23.00 selama 15 hari pengamatan. Alat pengukuran menggunakan Tide Staff. Data pengukuran pasut selanjutnya digunakan sebagai data validasi model. 3.4.4 Pengumpulan Data Sekunder 3.4.4.1 Data Batimetri Data batimetri diperlukan untuk memprediksi variasi pola arus dan kecepatannya. Data batimetri merupakan hasil digitasi kedalaman dari peta batimetri yang dikeluarkan oleh Dishidros TNI-AL. Peta batimetri perairan Teluk Lampung disajikan pada Gambar 5.
24
Gambar 5 Peta batimetri Teluk Lampung (sumber: Dishidros TNI-AL). 3.4.4.2 Data Pasang Surut Parameter pasang surut digunakan sebagai data masukan di syarat batas terbuka untuk mengetahui proses yang membangkitkan proses hidrodinamika. Data pasang surut yang digunakan sebagai data masukan model adalah data prediksi pasang surut untuk daerah Bakauheni tahun 2007 yang dikeluarkan oleh Dinas Hidro Oseanografi (Dishidros) TNI AL. Data parameter pasang surut bervariasi terhadap waktu dan konstant sepanjang daerah syarat batas terbuka. Grafik data pasang surut yang digunakan sebagai masukan model disajikan pada Gambar 6.
Gambar 6 Pola pasang surut masukan model.
25
3.4.4.3 Data Debit Sungai Aliran sungai yang masuk ke perairan teluk sangat besar pengaruhnya dalam menentukan salinitas permukaan, selain sebagai sumber masukan nutrien ke dalam perairan. Data debit aliran sungai yang masuk ke perairan Teluk Lampung merupakan data debit rata-rata bulanan yang diperoleh dari Badan Pengelolaan Daerah Aliran Sungai (BAPEDAS) Dinas Kehutanan Propinsi Lampung. Data debit sungai juga digunakan untuk menghitung masukan nutrien dari daratan dengan persamaan: (1) dimana L adalah jumlah nutrien yang masuk dari sungai (ton/bulan), Q adalah debit sungai (l/detik), C adalah konsentrasi nutrien (mg/l), 2592000 adalah faktor konversi dari detk ke bulan dan 1000000 adalah faktor konversi dari mg ke ton 3.4.4.4 Data Salinitas Permukaan Dalam perhitungan simulasi model diperlukan data sebaran salinitas permukaan yang memiliki grid yang sama dengan grid batimetri. Data sebaran salinitas sebagai data masukan awal nilainya konstan terhadap ruang dan waktu. Data salinitas permukaan diperoleh dari www.nodc.noaa.gov/argo/float. dengan resolusi 2,5o x 2,5o. 3.4.4.5 Data Meteorologis Data meteorologis yang digunakan meliputi data arah dan kecepatan angin, temperatur udara, radiasi sinar matahari, kelembaban relatif, tekanan atmosfer, penutupan awan, dan presipitasi (hujan) disajikan pada Gambar 7 dan 8. Data meteorologis diperoleh dari ECMWF (European Center for Medium Range Forcasting) yang diunduh dari situs www.ecmwf.int. Data meteorologis ini merupakan data analisis ulang dan interpolasi dari data meteorologis yang diperoleh dari berbagai pusat pengamatan dan parameter meteorologi dunia. dengan resolusi spasial 2,5o x 2,5o dan interval setiap 3, 6, dan 12 jam pada ketinggian 10 m diatas permukaan laut dengan format NetCDF.
26
(a)
(c)
(e)
(b)
(d)
(f)
Gambar 7 Mawar angin masukan model tahun 2007 pada Bulan Januari (a), Februari (b), April (c), Mei (d), Juli (e) dan Agustus (f).
27
Gambar 8 Data meteorologis masukan model. 3.5 Analisis Contoh Untuk Verifikasi Model Sumber nutrien yang masuk perairan menggunakan data hasil analisa contoh yang diukur di stasiun pengambilan sampel. Analisis contoh dilakukan untuk mendapatkan data awal model yaitu dengan digunakan nilai terendah dari hasil analisis contoh.
28
3.5.1 Nutrien Konsentrasi nutrien yang berupa nutrien anorganik terlarut diperoleh dengan menganalisis sampel air menggunakan spektrofotometer. Analisa dilakukan untuk memperoleh konsentrasi PO4-P, NO3-N dan NH4-N. Sebelum dianalisa sampel air sebanyak 250 ml disaring menggunakan kertas saring nucleopore (diameter 47 mm dan porositas 0.2 µm). penyaringan dilakukan setidaknya 6 jam setelah sampling. 3.5.2 Karbon Organik Partikulat Konsentrasi karbon organik partikulat dihitung dengan pendekatan pengukuran total padatan tersuspensi (TSS) dengan persamaan (Bruce et al., 2006): (2) dimana LI adalah proporsi C yang hilang setelah pembakaran pada suhu 550oC selama 1 jam, nilainya 0.6 – 0.925 pada saat terjadi ledakan populasi diatom (Hipsey et al., 2006). Nilai 0.5 merupakan asumsi bahwa berat unsur C adalah 0.5 dari berat material organik (Bruce et al., 2006). 3.5.3 Konversi Klorofil-a dan Zooplankton ke Karbon Penentuan konsentrasi klorofil-a dilakukan dengan menggunakan metode spektrofotometer. Konsentrasi karbon fitoplankton diperoleh dengan mengkonversi klorofil-a ke karbon dengan asumsi 1 mg/L klorofil-a sama dengan 30 mg/L karbon (Chapra, 1997). Untuk keperluan simulasi konsentrasi karbon zooplankton dihitung menggunakan rasio 1 mg berat kering individu zooplankton setara dengan 1 µg C (Yanagi et al., 1997). Perbandingan antara berat kering dengan berat basah zooplankton menggunakan perbandingan 1 : 10 (Chapra, 1997). Berat basah dihitung dengan persamaan : (3)
29
dimana, BBS: berat zooplankton basah (mg); P1: berat sampel dengan kertas saring (mg/m3); P2 : berat kertas saring tanpa sampel; dan VS : volume sampel air tersaring (m3). 3.6 Model Hidrodinamika Simulasi hidrodinamik menggunakan kode pemodelan komputasi ELCOM , proses utamanya meliputi bahang permukaan, transpor massa dan momentum, dinamika lapisan tercampur dan limpasan aliran masuk (inflow). Data meteorologi digunakan untuk menentukan penetrasi bahang dari radiasi gelombang pendek dan fluks bahang permukaan yang disebabkan oleh evaporasi, bahang sensible, radiasi gelombang panjang dan tekanan angin. Angin permukaan menggerakkan momentum dan energi kinetik turbulen ke lapisan permukaan sehingga berperan dalam percampuran vertikal. Dasar persamaan transpor yang digunakan dalam ELCOM adalah persamaan Navier Stokes yang dirata-ratakan terhadap waktu dengan peratarataan Reynold yang ditulis sebagai (Hodges and Dallimore, 2010):
(4)
dimana u, v, w adalah komponen kecepatan pada sumbu horisontal x, y dan vertikal z;
adalah densitas air; p adalah tekanan dan
adalah
komponen Coriolis. Persamaan (4) dapat disimpulkan adanya pergerakan partikel yang tak terkendali dan transfer energi, sehingga untuk dua hal tersebut dilakukan perata-rataan dengan aturan sebagai berikut: (5) (6) dimana
adalah rata-rata pengukuran,
adalah deviasi dari rata-rata,
adalah nilai pengukuran, T adalah periode perata-rataan, dan a adalah konstanta.
30
Jadi dari persamaan (4) komponen u, v, w, dan p dapat diganti dengan ,
,
, dan
yang menyatakan nilai
pengukuran adalah nilai rata-rata ditambah deviasi atau komponen turbulen. Hasil proses perata-rataan dari persamaan gerak menghasilkan persamaan Navier Stokes dengan perata-rataan Reynold (RANS = Reynold Averaging Navier Stokes) sebagai berikut:
(7)
Hasil perata-rataan menghasilkan peubah baru ,
,
( , , dan
, dan
,
,
,
,
,
yang akan dinyatakan dengan derivatif aliran rata-rata
) yang dianalogikan dengan Hukum Viskositas Newton
,
dimana peubah ini dikenal sebagai Gesekan Reynold atau Gesekan Turbulen:
(8)
dimana
,
, dan
adalah viskositas Eddy masing masing pada sumbu x, y,
dan z. Pada sumbu x komponen Coriolis
<<
, sehingga
cukup kecil untuk ditiadakan. Jika persamaan (8) disubtitusikan ke persamaan (7) dan komponen Coriolis
diganti dengan simbol f akan diperoleh :
(9)
Untuk keperluan penyerderhanaan maka semua komponen simbol U, V dan W maka persamaan (9) menjadi:
diganti dengan
31
(10)
dimana U(x,y,z,t), V(x,y,z,t) dan W(x,y,z,t) adalah komponen kecepatan pada arah horizontal x, y dan vertikal z; t adalah waktu; g adalah percepatan gravitasi; adalah koefisien viskositas Eddy. Pada sumbu z semua komponen terlalu kecil jika dibandingkan
dan g
sehingga cukup kecil untuk ditiadakan, sehingga persamaan transpor pada sumbu z dapat ditulis : (11) Jika
;
dan komponen koefisien viskositas Eddy
dapat digantikan dengan simbol , maka secara sederhana persamaan transpor (10) dalam model ELCOM dapat ditulis sebagai: .
Persamaan transpor:
(12)
Kontinuitas: (13)
Kondisi batas momentum pada Permukaan bebas: (14) Dasar dan sisi: (15)
Transpor skalar (16)
32
Kondisi batas skalar (17)
Evolusi permukaan bebas (Free Surface Evolution) (18)
Gesekan angin pada permukaan bebas (19)
Input momentum oleh angin (20)
dimana, U: kecepatan Reynold yang dirata-ratakan atas waktu; i, j, k,m: komponen ruang; α,β: komponen ruang horizontal; η: tinggi permukaan bebas; ρ’: densitas anomali; ρo: densitas acuan; g: konstanta gravitasi; f : konstanta Coriolis; skalar;
: komponen permutasi tensor; ν: viskositas molekular; C: konsentrasi : kecepatan vektor angin dalam arah β;
pada ketinggian 10 m;
: koefisien bulk stress angin
: kecepatan wind shear pada arah α.
3.6.1 Penyelesaian Numerik Persamaan pembangun didiskretisasi pada grid solusi Cartesian dimana komponen kecepatan tunggal didefinisikan pada tiap sel yang berhadapan dan skalar didefinisikan pada pusat sel. Dalam persamaan diskretisasi sel yang berhadapan diwakili komponen subskrip i+1/2 sementara pada pusat sel diwakili nilai integer (i, j, k). Notasi dari Casulli and Cheng (1992) digunakan untuk penjelasan selanjutnya. Bentuk persamaan diskret akan digunakan subskrip untuk mewakili posisi dalam ruang diskret (i, j, k). Misalkan
mewakili nilai
vektor kecepatan di kolom air pada waktu n+1 pada posisi (i, j) yang ada pada solusi ruang dan waktu n* untuk seluruh k yang mencakup: (21) dimana
adalah ketinggian dari dasar domain pada titik (i, j);
ketingian permukaan bebas dan
adalah
adalah jumlah maksimum grid sel dalam
arah vertikal. Definisi yang sama juga digunakan untuk kuantitas vektor yang lain.
33
Dasar dari evolusi semi implisit dalam medan kecepatan pada komponen 2 dimensi dapat didiskretisasi pada solusi ruang * yang sama dengan pendekatan model TRIM (Tidal, Residual, Intertidal Mudflat) (Cheng et al., 1993) sebagai : (22) (23) dimana: G = sebuah sumber vektor eksplisit = implisitas permukaan bebas, nilai yang digunakan dalam ELCOM adalah diskretisasi langkah mundur euler ( Casulli and Cattani (1994) telah menunjukan bahwa metode langkah mundur Euler untuk solusi momentum dari persamaan hidrostatik dapat diperluas menjadi skema dua tingkat (persamaan 22 dan 23) yang secara formal akurat sampai ordo kedua (ketika
. Pemodelan
dengan grid kasar bagaimanapun
juga penambahan dalam akurasi numerik tidak selalu menghasilkan penambahan kemampuan model. Secara umum banyak simualasi yang dilakukan di daerah estuari dan reservoir, modus barotropik diselesaikan dengan kondisi CFL yang mungkin berkisar antara 5 - 10 atau lebih. Kondisi yang seperti ini diskretisasi semi implisit mungkin akan stabil tetapi ketepatan yang mewakili aliran fisik adalah fungsi dari aliran dengan pertimbangan karakter dari pemotongan galat adalah kritis untuk memahami kemampuan metode ini. Metode ordo pertama menyebabkan kesalahan komponen pada ordo kedua dan menghasilkan peredaman gelombang di permukaan bebas. Metode ordo kedua yang menyebabkan kesalahan adalah dispersif dan menghasilkan gelombang numerik pada permukaan bebas yang dibangkitkan sepanjang domain; pada umumnya menyebabkan gelombang barotropik linier meningkat pada lubang yang curam, menyebabkan kecepatan yang tinggi dalam wilayah yang dilokalisir sebagai gelombang permukaan yang dipengaruhi oleh topografi. Jadi, metode ordo pertama menghasilkan wakil yang bagus dari bentuk gelombang permukaan dan kecepatan barotropik lokal, tetapi menunjukkan perluasan peredaman respon inersia dari permukaan bebas. Sebaliknya metode ordo kedua memperbaiki energi dari gelombang permukaan dengan dispasi numerikal minimum, tetapi kurang
34
mewakili bentuk gelombang. Solusi hidrostatik dispersi gelombang mengakibatkan gaya lokal palsu yang melalui kolom air dan detrimental terhadap kemampuan solusi. Sebaliknya perluasan peredaman gelombang permukaan menyebabkan berkurangnya pergerakan dalam skala besar yang dihubungkan dengan respon barotropik ketika angin berkurang. Secara umum sistem yang sangat dianjurkan yang lebih baik adalah skema langkah mundur Euler sebagai energi gelombang diluar dua atau tiga periode yang sering tidak relevan pada ordo pertama. Menggunakan diskretisasi implisit dua tingkat (Casulli and Cheng, 1992), atau berbagai teknik eksplisit, matrik A dapat di gambarkan sebagai :
(24)
dimana
dalam matrik A adalah kondisi batas, sementara a, b dan c adalah
berturut-turut: (25) (26) (27) Koefisien
ditentukan dengan memilih teknik diskretisasi numerik, untuk
, komponen viskositas vertikal didiskretisasi menggunakan teknik Euler beda mundur, dan untuk model lapisan tercampur yang digunakan di ELCOM =0 dan diskretisasinya secara eksplisit dengan A adalah nol untuk semua komponen dalam diagonal utama. Komponen G dalam persamaan (22) dan (23) dapat digambarkan sebagai : (28) (29)
35
Operator L () mewakili diskretisasi advektif, B () mewakili diskretisasi baroklinik, D () mewakili diskretisasi difusi horizontal turbulen. Dalam ELCOM vertikal difusi dihitung menggunakan model percampuran vertikal. model percampuran diwakili oleh operator seperti: (30) dimana operator percampuran M akan dibahas kemudian. Dalam pendekatan TRIM (Casulli and Cheng, 1992) diterapkan diskretisasi difusi horizontal dalam jalur asal, penambahan kompleksitas tidak menemukan adanya keuntungan signifikan dalam akurasi. Selanjutnya dalam diskretisasi difusi horizontal (Dx, Dy) dari persamaan (28) dan (29) didiskretisasi menggunakan stensil ordo kedua seperti : (31) Komponen baroklinik (B) dalam arah x didiskretisasi sebagai, (32) dimana : F adalah sel yang memiliki permukaan bebas. Persamaan yang sama untuk difusi dan komponen baroklinik juga diperoleh untuk arah sumbu y. Penghilangan komponen difusi vertikal dalam persamaan transport momentum dan transport scalar dimungkinkan dengan menghilangkan inversi matrik tridiagonal untuk masing-masing komponen kecepatan horizontal dan transport scalar yang diperlukan untuk tiap kolom air (i, j) dalam skema TRIM. Persamaan evolusi permukaan bebas dapat di diskretisasi sebagai:
(33) dimana
adalah vektor grid spasi vertikal, dan operator
dan
mengindikasikan perbedaan diskret, seperti : (34) Substitusi persamaan (22) dan (23) dalam persamaan (33) menjadi persamaan pentadiagonal untuk tinggi permukaan bebas pada waktu n+1 yang
36
siap diselesaikan dengan menggunakan metode konjugasi gradien yang mirip dengan TRIM (Casulli dan Cheng, 1992). Salah satu kesulitan dalam model numerik sampai sekarang adalah skala kondisi aliran geofisik dengan kisaran yang lebar. Dalam hal tertentu gelombang internal mungkin sesekali menghasilkan pergerakan vertikal yang kuat pada wilayah yang relatif kecil. Ketika resolusi dari gelombang internal dianggap penting, pilihan metode numerik diarahkan pada kebutuhan untuk akurasi dan stabilitas dalam porsi yang kecil untuk keseluruhan medan aliran, sementara banyak metode diskretisasi eksplisit spasial stabil untuk CFL
O(0.5) secara umum kurang ketika arah aliran tidak sejalan dengan grid. Kelemahan metode kuadratik adalah komputasi yang memerlukan mesin dengan kemampuan tinggi. Bagaimanapun untuk daerah dengan CFL rendah (CFL<0.1), solusi dari metode kuadratik semi langrangian didominasi oleh komponen dalam tujuh titik stensil yang berlawanan yang dihasilkan oleh diskretisasi kuadratik berlawanan. Kesamaannya dapat dieksploitasi untuk mengurangi persyaratan komputasi pada daerah dengan CFL rendah tanpa secara signifikan mengurangi akurasi dari keseluruhan metode solusi. Konsep menerapkan skema yang berbeda pada daerah yang berbeda dapat disamaratakan dalam konsep metode numerik hibrida. Metode hibrida umum adalah satu dimana skema solusi yang berbeda diterapkan dalam daerah aliran yang berbeda didasarkan atas kriteria terukur dalam medan aliran. Untuk tujuan tersebut kriterianya adalah CFL, dan penerapan satu teknik diskretisasi untuk CFL rendah dan teknik yang berbeda untuk CFL tinggi. Metode hibrida sekarang telah diuji untuk dua tingkat menggunakan diskretisasi kuadratik semi Langrangian untuk daerah dimana 02 masih dipertanyakan, sehingga batas atas merupakan persyaratan yang tidak masuk akal. Untuk daerah dengan CFL>2 metode menerapkan diskretisasi semi Langrangian
37
linier sebagaimana diterapkan dalam TRIM untuk meminimalkan usaha mereposisi stensil. Untuk aliran terstratifikasi metode semi Langrangian digunakan dengan akurasi dan kestabilan ketika 0.11 dalam arah horizontal merupakan langkag penting berikutnya dalam memilih metode numerik. Langkah waktu maksimum secara umum dibatasi juga oleh pembangkitan kecepatan gelombang baroklinik pada daerah dengan stratifikasi terkuat atau difusi numerik maksimum yang dapat diterima dalam transpor skalar. Bentuk semi Langrangian dari adveksi diperoleh dengan menemukan pendekatan dalam ruang kontinum ( titik "Langrange") yang akan diadveksikan menjadi titik diskret (i, j, k) oleh medan kecepatan (U, V, W) sepanjang langkah waktu t. Posisi partikel (i, j, k) secara numerik dibariskan ulang sepanjang garis lurus yang diwakili oleh medan kecepatan U, V, W. medan U, V, W diperoleh dari tingkat waktu tunggal atau berganda, tergantung pada akurasi dan kompleksitas komputasi yang diinginkan. Pada tingkat waktu tunggal metode semi Langrangian linier, titik Langrange ditentukan menggunakan: (35) (36) (37) Nilai variabel trilinier berlawanan:
pada titik Langrangian ditentukan menggunakan interpolasi
38
=
(38)
Wakil diatas merupakan sebuah stensil delapan titik untuk metode semi Langrangian dengan interpolasi linier. Sebagai pendekatan medan aliran untuk aliran seragam 1D, teknik semi Langrangian mengurangi menjadi sebuah metode linier berlawanan untuk CFL<1. Memang kemungkinan lebih mudah untuk berfikir bahawa metode semi Langrangian sebagai stensil 3D liniear berlawanan yang dengan sukses direposisi untuk kondisi CFL tinggi. Selanjutnya bahwa metode semi Langrangian linier memperlihatkan tingkat penolakan dari difusi numerik yang merupakan ciri dari metode linier berlawanan. Hal ini dapat diperbaiki menggunakan metode kuadratik untuk interpolasi dari nilai pada titik Langrange. Metode semi Langrangian kuadratik meluas dari 8 titik stensil berlawanan dengan interpolasi trilinier menjadi 27 titik stensil berlawanan enggunakan interpolasi polinomial Langrangian kuadratik. Notasi yang digunakan Casulli and Cheng (1992) operator menginterpolasi medan kecepatan (yang berhadapan dengan x) ke posisi (i+1/2-a, j-b, k-d), dimana a, b, d adalah bilangan ril yang mewakili penggantian bentuk dari posisi (i+1/2, j, k). sebuah pendugaan asal dari jalur pergerakan partikel sepanjang medan kecepatan pada waktu n yang berakhir pada posisi (i+1/2,j ,k) setelah waktu t. Notasi yang sama digunakan pada yang menghadap y pada j+1/2. Untuk menyatakan jalur titik asal dengan pangkat (p) sehingga operator adveksi L() pada persamaan (28) dan (29) menjadi: (39)
39
Gambar 9
Ilustrasi garis komputasi Euler-Langrangian 2D menggunakan interpolasi kuadratik Langrangian (Hodges and Dallimore, 2010).
Gambar 9 menjelaskan bahwa vektor kecepatan A dengan komponen UA dan VA digunakan untuk melacak jalur partikel dari posisi (i, j) kembali ke basis B menggunakan momentum sub langkah waktu dt. Vektor kecepatan B dihitung dari sembilan titik grid berlawanan dari vektor kecepatan pada posisi (i, j). Vektor B digunakan untuk melacak jalur partikel kembali ke basis vektor kecepatan C yang diinterpolasi lagi dari 9 (sembilan) titik grid disekitarnya. Hal ini diulang sebanyak n waktu dimana ndt = t. jika basis posisi vektor kecepatan tidak diisi dengan sembilan titik berlawanan, stensil yang berlawanan harus direposisi. dalam kode yang sekarang, interpolasi linier digunakan untuk kejadia yang jarang ketika reposisi diperlukan. Vektor akhir dihasilkan dari operator Euler-Langrangian dari persamaan (24) dan (25) yaitu
. Jumlah dari sub langkah waktu
(n) dapat diatur seperlunya, dengan nilai tinggi yang menyediakan akurasi lebih besar dan biaya komputasi yang lebih tinggi. Dicatat bahwa n = 1 dimanapun berhubungan dengan diskretisasi kuadratik berlawanan dan akurasi yang rendah dengan ciri setidaknya CFLa<<1. Aturan umum, n minimum diatur sebagai fungsi lokal dari grid dan medan aliran seperti
.
40
Gambar 10 Skema interpolasi kuadratik Langrangian 3D dengan interpolasi berurutan dalam k, j kemudian i. Untuk Jelasnya, ilustrasi ini menunjukan interpolasi pada grid yang seragam, bagaimanapun metode ini dapat diterapkan untuk grid tak seragam tanpa adaptasi lebih lanjut (Hodges and Dallimore, 2010). Proses perhitungan jalur asal dalam 2D untuk interpolasi bilinier telah dibahas oleh Casulli and Cheng (1992) dan digambarkan pada Gambar 9 untuk interpolasi Langrangian kuadratik. Interpolasi Langrangian pada grid yang tidak seragam mempertimbangkan koordinat ruang (x, y, z) yang menghubungkan dengan indek komputasi (i, j, k). interpolasi Langrangian kuadratik 3D untuk menemukan nilai
pada berbagai titik (x(p)(p), y(p)(p), z(p)(p) dihitung dalam
tiga langkah seperti diilustrasikan dalam Gambar 10 yang memerlukan 9 interpolasi vertikal dari bentuk: (40) dimana subskrip (i, j) menyatakan posisi garis vertikal yang diinterpolasi, sementara
dan
adalah ±(0, 1, 2) dengan tanda ditentukan oleh arah berlawanan
dari stensil (sebagai tanda untuk kenaikan k pada posisi subskrip U) koefisien polinomial Langrangian untuk tiap garis (i, j) di hitung dari persamaan koefisien Langrangian standar (Al-Khafaji and Tooley, 1986 dalam Hodges and Dallimore, 2010) sebagai:
41
(41) dimana
adalah koordianat vertikal dari titik interpolasi dan tanda ditentukan
untuk mendapatkan stensil berlawanan. Interpolasi vertikal diikuti oleh 3 interpolasi horizontal pada arah y dalam bentuk: (42) Akhirnya, interpolasi tunggal pada arah x diterapkan sebagai: (43) Koefisien Langrangian dalam persamaan (40) dan (41) dihitung menggunakan persamaan (42) dengan y atau x digantikan untuk z yang sesuai. Stensil kuadratik digunakan untuk interpolasi Euler-Langrangian karena menguntungkan sebagai pengurang peredaman gelombang internal yang terjadi dengan 8 titik stensil linier; jadi meningkatkan kemampuan metode untuk memecahkan pergerakan bebas dari basin terstratifikasi. Sementara peningkatan ini perlu untuk dalam memodelkan reservoir terstratifikasi. Hal yang kemungkinan kurang penting dalam model estuari dimana yang bergerak dominan adalah aliran fisik. Komputasi ekstra memerlukan stensil kuadratik sedikitnya diperbaiki dengan kemampuan menghitung sumber komponen untuk aliran dalam kisaran 1
42
(46) Sebagai momentum percampuran dan komponen sumber, operator percampuran dalam persamaan (44) mewakili percampuran vertikal oleh komponen tekanan Reynold.
mewakili sumber skalar (yaitu transfer bahang
sepanjang permukaan bebas ke dalam lapisan tercampur oleh angin). Persamaan (45) adalah adveksi dari medan skalar yang dipecahkan medan aliran dan persamaan (46) dalah difusi horizontal karena pergerakan turbulen. Untuk lebih jelasnya, adveksi (persamaan 45) didefinisikan setiap langkah waktu t. Bagaimanapun ketika MAX(CFLα)>1 sub skala langkah waktu δt digunakan, dimana mδt = t dan persamaan (45) diiterasi sebanyak m kali. Derivasi berikut ini substitusi dari δt untuk t dan n+m t untuk n+1 membuat persamaan mencerminkan proses iterasi sub langkah waktu. Dalam notasi diferensial, bentuk konservatif persamaan transpor adalah: (47) dimana Ω adalah kontrol volum,
,
,
adalah area permukaan yang
berhadapan dengan kontrol volum. Untuk kejelasan dalam bentuk diskret, akan lebih berguna untuk menghilangkan notasi subskrip (i, j, k) untuk semua variabel dipusatkan sehingga ditulis sebagai
ditulis dengan bentuk sederhana sebagai C, dan . Konsentrasi skalar yang diadveksikan C* ditulis sebagai
fluks adveksi konservatif sepanjang rangkaian n* sel sebagai: (48) dimana operator dalam bentuk
didefinisikan dalam persamaan (32). Q adalah
fluks skalar yang melalui sel yang berhadapan, didefinisikan pada serangkaian sel n* untuk yang berhadapan dengan i+1/2 sebagai: (49) Serupa dengan definisi yang diterapkan pada yang berhadapan j+1/2 dan k+1/2. Tidak ada fluks pada sel yang berhadapan paling atas yang berisi permukaan bebas, sehingga
, dimana k = F adalah sel yang berisi
43
permukaan bebas. Hal ini mengikuti bahwa
dan
untuk setiap sel (i, j, F) dalam rangkaian n*: (50) Jadi, untuk semua sel dalam rangkaian n* (termasuk sel permukaan bebas): (51) Sejak konsentrasi skalar diperbaharui di pusat sel, adalah perlu untuk mendefinisikan metode interpolasi untuk nilai sel yang berhadapan seperti . Penyaringan pembatasan fluks ULTIMATE diterapkan dengan interpolasi ordo ketiga QUICKEST dilakukan dengan terutama sekali dalam menjaga medan skalar monotonik sementara membatasi difusi numerik. metode ULTIMATE QUICKEST konservatif dibatasi oleh CFL<1 dalam arah koordinat. Skema semiimplisit sekarang ini masih menyisakan kestabilan pada CFL tinggi, sehingga algoritma ULTIMATE QUICKEST harus dihitung dengan sukses melampaui sub skala langkah waktu sehingga maksimum CFLα kurang dari satu pada tiap sub langkah waktu. dalam aplikasinya model dengan resolusi kasar yang terstartifikasi kuat akan memiliki pembatasan langkah waktu berdasarkan pada mode baroklinik dalam solusi momentum dan CFLα>1 tidak pernah terjadi. Komponen difusi horizontal didiskretisasi untuk mengasilkan waktu skalar n+1 melampaui ruang solusi n*: (52) dimana Dx dan Dy adalah operator beda hingga untuk turunan kedua, sama dengan solusi kecepatan, kapanpun lokasi n+1 yang tidak pada waktu solusi ruang n* diperbarui menggunakan konsentrasi sel disekitarnya. Momentum dimasukkan oleh gaya angin dihilangkan pada lapisan batas dan proses turbulen di interior. pada model 2D perata-rataan kedalaman formulasi gesekan dasara Chezy-Manning telah diterapkan untuk menghitung hilangnya pada kolom air. Hal ini secara khusus berguna untuk menyediakan metode kalibrasi model perata-rataan kedalaman pantai/estuari untuk mereproduksi srangkaian data pasut yang diberikan. Dalam model 3D dengan stratifikasi data lapang yang detail untuk kalibrasi koefisien friksi dasar umumnya tidak tersedia.
44
Lebih lanjut tanpa transfer energi gelombang internal skala dasar ke gelombang skala subgrid, menjadi dipertanyakan apakah kondisi batas model dalam literatur dapat ditangkap secara aktual dinamika batas dan memprediksi lokasi yang tepat dan waktu menghilangnya serta fluks vertikal. kondisi batas sisi dinding (yaitu batas solid vertikal) jarang dimodelkan sebagai free-slip yang berdampak lebih sederhana implementasinya dalam metode numerik. Menegasikan gaya geser sisi dinding tidak akan menjadi lebih sederhana dalam pemodelan pergerakan. Masih ada banyak pekerjaan yang harus dilakukan dalam mengembangkan kondisi batas dasar dan sisi dalam model grid kasar. ELCOM menyediakan tiga bentuk dasar kondisi batas no-slip, free-slip, dan specified stress. Angin digunakan sebagai input momentum tekanan pada kondisi batas pada permukaan bebas dengan persamaan (Casulli and Cheng, 1992): (53) dimana
adalah viskositas Eddy dan
adalah tekanan angin. Kondisi batas ini
memerlukan solusi dari viskositas vertikal/difusi, (54) Lapisan tercampur yang dipengaruhi oleh angin yang termasuk permukaan bebas, dengan kedalaman (h) dihitung dalam bentuk diskret sebagai : (55) dimana ka dan kb adalah grid sel atas dan bawah yang mengindikasikan dari diskret lapisan tercampur karena pengaruh angin dalam kolom air (i,j) yang memiliki grid sel permukaan bebas
.
Pada ordo pertama dapat dilakukan pendekatan pendahuluan dari momentum angin sebagai distribusi yang seragam yang melalui lapisan tercampur (Imberger dan Peterson, 1990): (56) dimana
adalah ketinggian permukaan bebas dalam kolom air (i,j).
Perubahan bahang permukaan air dibangun dengan model standar bulk transfer (Hodges and Dallimore, 2010). Transfer energi yang melalui permukaan
45
bebas dipisahkan dalam komponen non penetratif dari radiasi gelombang panjang, transfer bahang sensible, dan kehilangan bahang melalui evaporasi. Pengaruh non penetratif dimasukkan sebagai sumber temperatur di lapisan permukaan dan lapisan tercampur, sehingga efek penetratif dimasukkan sebagai satu sumber dalam satu atau lebih lapisan grid pada basis penurunan eksponensial dan koefisien ekstingsi. Perubahan atau pertukaran pada bagian permukaan mencakup perpindahan panas yang disebabkan penetrasi gelombang pendek pada perairan dan fluks pada permukaan akibat penguapan, bahang sensibel (yaitu konveksi panas dari permukaan perairan menuju atmosfer) dan radiasi gelombang panjang. Radiasi gelombang pendek (280nm – 2800nm) pada umumnya diukur secara langsung. Radiasi gelombang panjang (>2800nm) diemisikan dari awan dan uap air, dapat diukur secara langsung atau dapat dihitung dari penutupan awan, suhu udara, dan kelembaban. Koefisien refleksi, atau albedo, dari radiasi gelombang pendek bervariasi diantara setiap perairan dan tergantung dari sudut matahari yang dibentuk dan kondisi gelombang permukaan saat pengukuran. Radiasi gelombang pendek dapat dibagi menjadi 4 komponen. ELCOM mengasumsikan nilai persen dari beberapa komponen Photosynthetically Active Radiation (PAR) 45% Near Infrared (NIR) 41% Ultra Violet A (UVA) 3.5% Ultra Violet B (UVB) 0.5% Jarak dari penetrasi radiasi ke dalam kolom perairan tergantung dari komponen peluruhan untuk setiap lebar band. ELCOM mengijinkan pengguna untuk mengatur koefisien ekstingsi untuk setiap band dalam run_elcom.dat file tetapi nilai pada umumnya adalah PAR 0.25/m NIR 1/m UVA 1/m UVB2.5/m Bila kualitas air disimulasikan melalui CAEDYM maka koefisien extinction PAR dihitung melalui CAEDYM. Kedalaman penetrasi dari radiasi gelombang pendek tergantung pada radiasi gelombang pendek bersih yakni penetrasi yang
46
sampai ke permukaan air dan koefisien ekstingsi. Persamaan yang diberikan untuk penetrasi radiasi matahari bersih dapat ditulis menjadi: (57) dimana air,
adalah radiasi gelombang pendek yang mencapai permukaan adalah penetrasi radiasi gelombang pendek bersih pada permukaan air
dan (sw) dan
adalah gelombang pendek albedo dipermukaan air yang
diberikan dengan persamaan:.
(58)
= 0.08,
= 0.02, D < standar banyaknya hari dalam satu tahun (365)
dan d adalah hari ke- dalam tahun. Penetrasi gelombang pendek pada kedalaman menurut hukum Beer-Lambert adalah : (59) dimana z adalah kedalaman dibawah permukaan air dan
adalah koefisien
atenuasi. Sehingga energi gelombang pendek per unit area yang memasuki lapisan k adalah : (60) atau (61) Untuk tujuan bahang diasumsikan bahwa semua Q dikonversi menjadi bahang. Jika terdapat energi gelombang pendek yang berlebihan pada dasar kolom perairan, ELCOM merefleksikan sebagian dari energy kembali ke dalam domain (mencapai 90%). Energi ini dibolehkan untuk mengalami propagasi kembali melalui kolom perairan, menurut Hukum Beer-Lambert. Radiasi gelombang panjang dihitung melalui salah satu dari tiga metode, tergantung dari asupan data. Tiga pengukuran yang dimungkinkan adalah: kejadian radiasi gelombang panjang, radiasi gelombang panjang bersih, dan penutupan awan.
47
Menghitung kejadian radiasi gelombang panjang, berarti memperhitungkan albedo dan radiasi gelombang panjang yang diemisikan dari lapisan permukaan perairan. Penetrasi Radiasi gelombang panjang pada lapisan permukaan berlaku: (62) dimana
adalah albedo untuk radiasi gelombang panjang, yang dianggap
konstan = 0.03 (Hendeson-Sellers, 1986 dalam Hodges and Dallimore, 2010). Radiasi gelombang panjang diemisikan dari lapisan permukaan perairan yang dideskripsikan Tennessee Valley Authority (1972) dalam Hodges and Dallimore (2010) adalah: (63) dimana
adalah emisivitas dari lapisan permukaan perairan (=0.96),
adalah
) dan Tw adalah
konstanta Stefan-Boltzmann (
temperatur mutlak dari lapisan permukaan perairan (yaitu temperatur dari lapisan permukaan). Densitas Energi dari radiasi gelombang panjang bersih yang dikumpulkan pada lapisan permukaan pada periode Δt sehingga menjadi (64) Dengan menggunakan radiasi gelombang panjang bersih, berarti memperhitungkan albedo pada lapisan permukaan perairan. Densitas energi Radiasi gelombang panjang yang dikumpulkan pada lapisan permukaan pada periode waktu Δt menjadi (65) Radiasi Gelombang panjang dapat juga diestimasi dari kondisi atmosfer, menggunakan fraksi penutupan awan
. Densitas energi kejadian
radiasi gelombang panjang pada permukaan perairan dapat diestimasikan menjadi (66) dimana: (67) Subskrib a mengacu pada sifat dari udara (68)
48
dimana
. Sebelumnya gelombang panjang emisi adalah (69)
Radiasi gelombang panjang bersih menjadi (70) Kehilangan bahang sensibel dari permukaan perairan dalam waktu periode dapat dituliskan sebagai: (71) dimana Cs adalah koefisien transfer kalor sensible dengan kecepatan angin diukur 10 meter, ketinggian referensi yang diukur dari permukaan air (= 1.3x10-3), ρa adalah densitas udara dalam kg/m3, Cp kalor spesifik dari udara pada tekanan yang tetap (=1003Jkg-1K-1), Ua adalah kecepatan angin yang diukur dengan menggunakan referensi ketinggian 10 m dalam m/s, dengan temperatur dalam Celsius ataupun Kelvin. Persamaan fluks bahang evaporatif (Fischer et al. (1979) dalam Hodges and Dallimore (2010) adalah: (72) dimana P adalah tekanan atmosfer dalam pascal, CL adalah koefisien kalor laten transfer (= 1.3 x 10-3) dengan pengukuran kecepatan angin dengan tinggi 10 m, ρa adalah densitas udara dalam kg/m3 , LE adalah kalor laten dari penguapan air (= 2.453 x 106J kg-1), Ua adalah kecepatan angin dalam m/s dengan ketinggian referensi 10 meter, ea tekanan uap air di udara, dan es adalah tekanan uap air jenuh pada suhu permukaan air Ts; kedua tekanan uap air diukur dalam pascal. Kondisi yang berlaku adalah
, sehingga tidak ada efek kondensasi yang
dimasukkan dalam perhitungan. Tekanan uap air jenuh es dihitung dengan menggunakan rumus MagnusTetens (Tennessee Valley Authority (1972) dalam Hodges and Dallimore (2010): (73) dimana Ts dalam derajat Celsius dan es dalam pascal. Sehingga, total densitas energi non-penetrative yang terkumpulkan pada lapisan permukaan selama periode
dapat ditulis sebagai berikut
49
(74) Perubahan massa pada sel lapisan permukaan (layer nomor ke N) terhadap flux kalor laten yang dihitung menjadi (75) dimana dx dan dy adalah ukuran grid dari sel lapisan permukaan dan Lv adalah kalor laten dari penguapan air. Bila diasumsikan temperatur air hujan sama dengan sel lapisan permukaan. Salinitas dan kualitas air diatur sama dengan nol. Perubahan massa lapisan permukaan (76) dimana r adalah kecepatan air turun dalam m/s Fluks massa permukaan didasarkan pada keseimbangan antara evaporasi dan hujan, ynag akan merubah massa dari sel lapisan permukaan. Pendekatan lapisan tercampur tiga dimensi pada percampuran neraca energi, transpor Turbulence Kinetic Energy (TKE) digunakan untuk menyediakan efek pergerakan dinamis 3 dimensi yang hasil percampuran merupakan keseimbangan dari : 1) TKE yang tersedia selama percampuran: TKEA 2) TKE yang dibutuhkan selama percampuran: Ereq 3) TKE yang keluar:
, dan
4) Residu energi percampuran: EM Hal ini sangat berguna untuk mencirikan dua tipe kejadian percampuran pada fluida terstratifikasi yaitu percampuran konvektif dari gradien densitas yang tak stabil yang menurunkan energi potensial fluida dan melepaskan TKE dan percampuran dari gradien densitas yang stabil yang menghilangkan TKE dan meningkatkan energi potensial. Gambar 11 menjelaskan bahwa pendekatan yang sekarang dilakukan adalah dengan mempertimbangkan sesuatu yang lebih sederhana dari kebanyakan model 1-D yang dalam persamaan deferensial menggabungkan ketebalan lapisan, laju pengiringarusan, dan transfer bahang melalui permukaan tidak digunakan. Dasar praduga pendekatan model ini adalah bahwa resolusi grid vertikal untuk 3D adalah terlalu kasar untuk mencukupi penghitungan solusi persamaan diferensial
50
untuk evolusi lapisan tercampur dan transfer bahang dengan gradien kondisi batas. Memang keperluan memodelkan dengan resolusi grid vertikal yang kasar adalah untuk menunda daya dorong persamaan diferensial difusi Eddy vertikal yang umumnya digunakan dalam percampuran.
Gambar 11 Perkembangan dari lapisan tercampur karena pendinginan permukaan dan stratifikasi tak stabil (a) stratifikasi stabil pada langkah waktu dimulai, (b) pendinginan permukaan menciptakan profil densitas tak stabil, (c) grid sel tak stabil yang telah tercampur (Hodges and Dallimore, 2010). Termodinamika permukaan ELCOM masuk kedalam sistem sebagai perubahan diskret terhadap perubahan struktur temperatur pada grid sel teratas dibandingkan sebagai kondisi batas pada transpor atau percampuran. Pada permulaan langkah waktu fluks bahang non penetratif (radiasi gelombang panjang, transfer bahang konvektif) ditambahkan pada grid sel paling atas sementara radiasi matahari ditambahkan ke kolom air menggunakan peluruhan eksponensial terhadap kedalaman. Transfer bahang merubah stratifikasi densitas juga menyediakan TKEA (gradien densitas takstabil yang dihasilkan oleh pendinginan bersih) atau meningkatnya Ereq (stratifikasi stabil yang dihasilkan oleh pemanasan bersih). Sekali medan densitas baru dihitung, proses percampuran dimodelkan berbasis pada lapisan demi lapisan tiap kolom air (i, j) dengan
51
membandingkan ketersediaan energi percampuran (TKEA) dari pembalikan konvektif, produksi geser, pengadukan angin dan energi yang disimpan (EM) terhadap peningkatan energi potensial yang diperlukan (Ereq) untuk mencampur grid sel diatas menjadi lapisan tercampur diatasnya. Percampuran dalam model dihitung dengan tahapan :
Menghitung input energi angin (77) dimana
adalah kecepatan gesekan angin dan Cn adalah koefisien
percampuran = 1,33
Jika kondisi batas digunakan, dihitung input energi dasar akibat pergeseran Edrag (78) dimana Cb adalah koefisien percampuran = 2,2
Dalam tiap kolom siklus dari sel permukaan ke sel dasar pada lapisan permukaan tercampur sampai Menghitung pembangkitan TKE karena shear, Eshear; (79) dimana Cs adalah koefisien percampuran = 0,15; S adalah gesekan=1 merujuk pada sel yang secara langsung dibawah lapisan tercampur, Menghitung energi yang diperlukan untuk percampuran Ereq (80) dimana g adalah koefisien gravitasi tereduksi; ml merujuk pada nilai di lapisan tercampur
Menghitung total energi yang tersedia jika dua sel total tercampur TKEmixed (81) dimana Cc adalah koefisien percampuran = 0,2
Menghitung waktu perkiraan selama percampuran TTKE (82) dimana Crr adalah koefisien = 50,0
52
Jika perhitungan waktu perkiraan tidak stabil didasarkan pada pembalikan kovektif Tconv (83)
Menghitung fraksi percampuran
f
(84)
Jika ada cukup energi dimana : Jika TKEmixed
akan terjadi percampuran
Jika TKEmixed
tidak terjadi percampuran
Jika
=1 semua skalar dan kecepatan dalam lapisan tercampur adalah
sama (tercampur total) .
(85)
Untuk tercampur sebagian (86)
dimana Cml adalah konsentrasi skalar dari lapisan tercampur; Cl adalah sel yang dicampur; tanda mengindikasikan nilai setelah percampuran; Ck adalah konsentrasi dari lapisan k dalam kolom k1 merujuk pada lapisan yang dicampur; dan kml-top adalah lapisan diatas lapisan tercampur.
Akhir siklus dari sel permukaan ke sel dasar, kelebihan energi percampuran yang keluar. (87)
dimana
adalah koefisien pengeluaran = 1,15; beberapa TKE tersisa setelah
keluaran transpor sebelum tersedia selama tahap percampuran berikutnya. Untuk mencampur sel dibawah lapisan tercampur (dari ketebalan kedalam lapisan tercampur (dari ketebalan
)
) fraksi percampuran
didefinisikan sebagai fraksi dari sel dibawahnya yang akan dicampur kedalam lapisan tercampur. kombinasi Potential Energy dari lapisan tercampur
53
(disimbolkan sebagai ml) dan sel akan dicampur setelah percampuran dengan persamaan: (88) dimana superskrip * mengindikasikan nlai setelah percampuran, jika dicatat bahwa: (89) dan (90) (91) kemudian,
(92)
(93)
(94)
(95)
(96)
(97)
(98)
54
3.6.2 Desain Simulasi Model Hidrodinamika Dengan menerapkan Radius Deformasi Rossby, maka dimensi lateral Teluk Lampung (50.000 m) jauh lebih kecil dari Radius Deformasi Rossby ( =54532768,4), sehingga efek Coriolis (f) dapat diabaikan. Untuk memenuhi kriteria stabilitas Courant-Friedichs-Lewy (CFL) dalam persamaaan momentum dengan berdasarkan pada kedalaman maksimum dan lebar sel maka didapatkan langkah waktu dengan persamaan sebagai berikut : ; t ≤15.65377 detik
(99)
dimana: : adalah lebar sel = x= y=500 meter : adalah kedalaman maksimum perairan daerah model (52 m) g
: adalah percepatan gravitasi bumi = 9.81 m/s2 Maka langkah waktu yang memenuhi syarat kestabilan CFL adalah 15, 65
detik, tetapi langkah waktu ( t) yang digunakan dalam simulasi adalah 15 detik. Daerah model dibagi dalam 85 x 109 sel dalam bentuk matrik dengan lebar sel (grid) x = y = 500 meter. Perubahan kedalaman diatur pada nilai konstan z = 2 meter, sehinggga jumlah grid vertikalnya akan bervariasi tergantung kedalaman perairan. Syarat batas tertutup merupakan daerah yang tidak memungkinkan massa air melewatinya, atau kecepatan dengan arah tegak lurus pantai adalah sama dengan nol. Syarat batas tertutup dapat dikatakan juga sebagai daerah yang mememiliki ketinggian lebih dari nol atau merupakan daerah daratan, sehingga berlaku persamaan : (100) Daerah model yang berbatasan dengan laut terbuka merupakan syarat batas terbuka, dimana pada simulasi ini syarat batas terbuka ditarik sebagai garis lurus antara daerah Tanjung Tikus disebelah barat hingga daerah Canti di sebelah timur.
55
3.7 Model Ekosistem 3.7.1 Persamaan Pembangun Model Ekosistem Model ekosistem CAEDYM telah digambarkan sebelumnya untuk model estuari (Robson and Hamilton, 2004 dalam Hipsey et al., 2009) dan untuk perairan pantai (Spillman et al., 2006 dalam Hipsey et al., 2009). Dalam model ini disimulasikan 2 kelompok fitoplankton dan satu kelompok zooplankton. Variabel yang disimulasikan meliputi karbon internal (IC), nitogen internal (IN), dan fosfor internal (IP) juga termasuk nutrien anorganik terlarut (PO4, NO3, NH4), materi organik terlarut (DOC, PON, POP), dan oksigen terlarut (DO). Proses yang menggambarkan interaksi dan perubahan antar variabel yang disimulasikan dalam CAEDYM disajikan pada Gambar 12.
Gambar 12 Interaksi yang mewakili model dalam CAEDYM (diadaptasi dan digambar ulang dari Hipsey et al., 2009).
56
3.7.2 Variabel yang Disimulasikan Keseluruhan variabel yang ditransportasikan oleh penggerak hidrodinamik di sajikan pada Tabel 7. dan deskripsi dari setiap parameter yang digunakan dalam pemodelan disajikan pada Lampiran 1, 2 dan 3. Parameter dasar sebagai kunci dari variabel dalam persamaan yang akan disimulasikan didefinisikan dalam persamaan 101 - 111, dimana j menyatakan pengidentikasi dari variabel dasarnya (Hipsey et al., 2009). Tabel 7
Variabel A AIN AIP a NA Z ZIN ZIP z NZ DOC POC DIC CO2 pCO2 pH TA TN PON DON NH4 NO3 TP POP DOP FRP DO I T ∆t ∆z ∆zbot ∆zsurf
Daftar dan deskripsi variabel yang disimulasikan dalam model CAEDYM Deskripsi
Unit satuan
Konsentrasi biomasa fitoplankton Konsentrasi internal nitrogen pada fitoplankton Konsentrasi internal fosfor pada fitoplankton Indek kelompok fitoplankton Jumlah kelompok fitoplankton yang disimulasikan Konsentrasi biomassa zooplankton Konsentrasi internal nitrogen pada zooplankton Konsentrasi internal fosfor pada zooplankton Indek kelompok zooplankton Jumlah kelompok zooplankton yang disimulasikan Konsentrasi karbon organik terlarut Konsentrasi karbon organik partikulat Konsentrasi karbon anorganik terlarut Konsentrasi karbon dioksida Tekanan parsial karbon dioksida Derajat keasaman, didasarkan pada konsentrasi ion H + Total alkalinitas Konsentrasi total nitrogen Konsentrasi nitrogen organik partikulat Konsentrasi nitrogen organik terlarut Konsentrasi amonium Konsentrasi nitrat Total konsentrasi fosfor Konsentrasi forfor organik partikulat Konsentrasi fosfor organik terlarut Fosfor reaktif terfilter Konsentrasi oksigen terlarut Intensitas cahaya Temperatur Langkah waktu Ketebalan lapisan sel Ketebalan lapisan sel diatas sedimen Ketebalan lapisan sel pada antar muka kolom air dan atmosfer
mg Chla L-1 mg N L-1 mg P L-1 mg C L-1 mg N L-1 mg P L-1 mg C L-1 mg C L-1 mg C L-1 mg C L-1 atm mg CaCO3 L-1 mg N L-1 mg N L-1 mg N L-1 mg N L-1 mg N L-1 mg P L-1 mg P L-1 mg P L-1 mg P L-1 mg O L-1 uE m-2 o C hari m m m
57
a) Fungsi yang mempunyai ketergantungan pada oksigen (101) (102) b) Fungsi yang mempunyai kergantungan pada temperatur (103) (104) c) Fluks sedimen terlarut (105) d) Pertukaran gas atmosfer (106) e) Penenggelaman (settling) (107) f) Mortalitas dan ekskresi biologis (108) g) Uptake biologis (109) h) Respirasi biologis (110) i) Grazing (111) 3.7.2.1 Cahaya Untuk produktivitas primer intensitas gelombang pendek pada permukaan di konversi menjadi komponen PAR (400 – 700 nm). PAR diasumsikan masuk kedalam kolom air mengikuti hukum Lambert-Beer persamaan (Hipsey et al., 2009): (112)
58
dimana
adalah PAR pada kedalaman z di bawah permukaan,
intensitas kejadian gelombang pendek dan
adalah
adalah koefisien ekstingsi PAR,
yang bervariasi secara spasial dan temporal yang mana dalam model akan mengikuti konsentrasi fitoplankton, partikel detritus dan anorganik, dan DOC (Hipsey et al., 2009) : (113) dimana
adalah atenuasi pada air murni (konstan), dan
adalah laju
peningkatan koefisien ekstingsi dengan pertambahan konsentrasi dari variabel yang berhubungan. 3.7.2.2 Oksigen Terlarut Dinamika DO dalam CAEDYM termasuk pertukaran dari atmosfer sedimen dan kolom air; konsumsi dari dekomposisi dan nitrifikasi bahan organik; produksi oksigen dari fotosintesis; dan konsumsi oksigen untuk respirasi zooplankton. Pertukaran atmosferik didasarkan pada model Wanninkhof (1992) dan persamaan fluks dari Rilley and Skirrow (1974) dalam Hipsey et al (2009).
(114) Keterangan 1: perubahan oksigen terlarut karena pertukaran atmosfer 2: perubahan oksigen terlarut dari fitoplankton fotosintesis dan respirasi dari fitoplankton dan zooplankton 3: perubahan oksigen terlarut selama proses nitrifikasi 4: perubahan oksigen karena kebutuhan reduksi dan oksidasi di sedimen 3.7.2.2 Siklus Karbon, Nitrogen dan Fosfor Kedua bentuk organik dan organik dari karbon, nitrogen dan fosfor baik yang terlarut maupun partikulat dimodelkan secara eksplisit mengikuti jalur perubahan POM-DOM-DIM. Siklus nitrogen termasuk proses tambahan dari denitrifikasi, nitrifikasi dan fiksasi N2 yang tidak ada dalam siklus karbon.
59
Keseimbangan DIC juga termasuk fluks CO2 dari atmosfer yang didasarkan pada perbedaan antara nilai pCO2 antara atmosfer dan kolom air. Perpindahannya dihitung mengikuti Wanninkhof (1992) dalam Hipsey et al (2009) dengan produk keterlarutan CO2 dari (Weiss, 1974 dalam Hipsey et al., 2009). Nilai fase gas dan cair CO2 dihubungkan dari hukum Henry untuk sirkulasi dari pCO2. Untuk menduga fraksi CO2 dari kolom DIC dan sistem penyangga karbonat, pH dan alkalinitas dimodelkan menurut Buttler (1982) dalam Hipsey et al (2009). Persamaan-persamaan pembangun dalam siklus karbon, nitrogen dan fosfor yang dimodelkan dalam CAEDYM adalah sebagai berikut (Hipsey et al., 2009): 1) Persamaan karbon -
Karbon dalam bentuk organik Perubahan oleh zooplankton
(115) Perubahan oleh fitoplankton
(116) Perubahan bentuk partikulat
(117) Perubahan bentuk terlarut
(118)
60 -
Karbon dalam bentuk anorganik Perubahan DIC (119) Total Alkalinitas (120) Konsentrasi CO2 (121) Tekanan Parsial CO2 (122) Derajat Keasaman (123)
2) Persamaan Nitrogen Nitrogen internal pada zooplankton (124) Nitrogen internal pada fitoplankton
(125) Perubahan nitrogen organik partikulat
(126) Perubahan nitrogen organik terlarut
(127)
61
Perubahan amonium
(128)
Perubahan nitrat
(129) Total nitrogen (130) 3) Persamaan fosfor Fosfor internal pada zooplankton (131) Fosfor internal pada fitoplankton (132)
Perubahan Fosfor organik partikulat
(133) Perubahan fosfor organik terlarut
(134) Perubahan fraksi fosfor reaktif (135) Total fosfor (136)
62
3.7.2.3 Dinamika Fitoplankton Nilai utama untuk biomassa fitoplankton, A, dikonfigurasi dalam satuan mgChlaL-1. Untuk tiap kelompok fitoplankton laju pertumbuhan potensial maksimum terjadi pada temperatur 25oC yang digandakan oleh nilai minimum dari pembatasan oleh cahaya, fosfor, dan nitrogen, sementara ada beberapa kemungkinan interaksi antara faktor pembatas (Rhee and Gotham, 1981 dalam Hipsey et al., 2009). Pembatasan cahaya pada pertumbuhan fitoplankton dikonfigurasi melalui pendefinisian karakter fitoplankton sebagai spesies fotoinhibitor atau nonfotoinhibitor. Pada jenis fotoinhibitor model yang dilakukan oleh Webb et al. (1974) dalam Hipsey et al (2009) digunakan untuk mengkuantifikasi fraksi pembatas dari laju potensial maksimum fiksasi karbon. Penghambatan pertumbuhan fitoplankton pada temperatur yang lebih tinggi fungsi temperatur fT1 digunakan dimana produktivitas maksimum terjadi pada temperatur TOPT dan jika temperatur meningkat maka produktivitas menurun pada temperatur maksimum yang diperbolehkan (TMAX). Temperatur di bawah temperatur standar (TSTD) maka pertumbuhan akan mengikuti hubungan eksponensial normal atau fungsi temperatur yang digunakan adalah fT2. Fosfor internal dan nitrogen dinamik dalam kelompok fitoplankton dimodelkan dengan model penyimpanan intra seluler dinamis yang dapat mengatur pertumbuhan melalui f(N) dan f(P). Fitoplankton akan memiliki konsentrasi variabel nutrien internal dengan pengambilan nutrien yang dinamis yang dibatasi oleh nilai maksimum dan minimum. Pengertian kehilangan melalui respirasi, mortalitas alami dan ekskresi dikelompokkan dalam istilah koefisien laju respirasi. Konstanta fres adalah fraksi respirasi murni yang hilang dan tidak termasuk mortalitas dan ekskresi. Kehilangan metabolisme nutrien dari mortalitas dan ekskresi adalah proporsional terhadap rasio penggandaan nutrien internal oleh koefisien laju kehilangannya. Model migrasi dan penenggelaman fitoplankton yang disimulasikan termasuk konstanta penenggelaman, hukum penenggelaman Stoke dan migrasi tanpa fotoinhibitor. Migrasi vertikal tanpa fotoinhibitor mengacu pada model
63
Kromkamp and Walsby (1990) dalam Hipsey et al (2009) yang diwakili oleh migrasi ke atas karena penyinaran dan migrasi ke bawah atau sedimentasi untuk mengisi penyimpanan nitrogen internal. Respon perpindahan nutrien diatur oleh penyimpanan nitrogen internal dalam sel dan asupan nitrogen eksternal ke dalam sel. Model yang mengacu pada hukum penenggelaman Stoke diasumsikan memilki respon yang beragam terhadap sintesis karbohidrat (fotosintesis) dan penggunaannya (respirasi), dan didasarkan pada perhitungan kepadatan secara dinamis. Hukum Stoke digunakan untuk menghitung kecepatan vertikal. Keseluruhan sub model dalam CAEDYM yang melibatkan fitoplankton di hitung berdasarkan persamaan-persamaan berikut ini (Hipsey et al., 2009). a) Fungsi uptake Uptake CO2 (137) Uptake FRP (138) Uptake NH4 (139) Uptake NO3 (140) Uptake N2 (141) dimana: (142) (143) (144) (145)
64
b) Respirasi (146) c) Ekskresi (147) (148) (149) (150) (151) (152) d) Sedimentasi dan migrasi vertikal
(153)
3.7.2.4 Zooplankton Model CAEDYM mengasumsikan bahwa zooplankton memiliki rasio C:N:P yang tetap dan tergantung pada rasio C:N:P dari berbagai sumber nutrien. Pertumbuhan zooplankton adalah fungsi dari temperatur dan laju memakan yang diuraikan dalam persamaan Michelis –Menten. Keseluruhan perhitungan sub model zooplankton yang dimodelkan dalam CAEDYM di uraikan dalam persamaan-persamaan berikut ini (Hipsey et al., 2009). a) Fungsi grazing -
Grazing karbon (154)
-
Grazing fosfor
(155)
65 -
Grazing nitrogen
(156) dimana: (157) (158) b) Respirasi (159) c) Mortalitas, engesti dan ekskresi (160) (161) (162)
(163) dimana: (164) (165)
(166) (167) 3.7.3 Desain Simulasi Model Ekosistem Komunitas fitoplankton yang disimulasikan menggunakan dua kelompok yaitu kelompok dinoflagelata dan kelompok diatom untuk menghitung secara kolektif fitoplankton yang lain. Total biomassa zooplankton yang disimulasikan
66
menggunakan satu kelompok yang diwakili oleh copepoda. Pemilihan kelompok yang mewakili fitoplankton maupun zooplankton berdasarkan hasil penelitian sebelumnya yang dilakukan oleh Damar (2002) yang menyatakan bahwa kelompok fitoplankton yang mendominasi di perairan Teluk Lampung adalah diatom dan dinoflagela, sedangkan kelompok zooplankton didominasi oleh copepoda dan protozoa. Fitoplankton dan zooplankton ditiadakan pada data inflow dan diatur pada nilai nol. Parameter yang akan digunakan untuk simulasi CAEDYM diperoleh dari berbagai literatur seperti yang disajikan pada Lampiran 1, 2, dan 3. Nilai konsentrasi nutrien pada data asupan digunakan nilai terendah dari data hasil observasi, nilai nutrien pada daerah batas terbuka mengacu pada penelitian LIPI di perairan Selat Sunda. Nilai nutrien internal sebagai asupan model diatur pada nilai konstan. Bakteri tidak dimodelkan secara eksplisit karena perannya dimasukkan dalam jalur mineralisasi bahan organik partikulat (POC, POP dan PON). Jadi ketersediaan bahan organik partikulat untuk pemangsaan zooplankton dianggap juga sebagai peran bakteri. Silika sebagai faktor pembatas tidak dimasukkan dalam model karena ketersediaan data. Nilai salinitas dan temperatur diperoleh dari model ELCOM yang juga merupakan penggerak dari model CAEDYM ini. Model CAEDYM dibangun dengan menggunakan langkah waktu 30 menit dan luaran model diatur dalam nilai harian. 3.8 Simulasi Model Data yang digunakan untuk input model sebagai data kondisi awal adalah data meteorologi dan data masukan aliran sungai. Data meteorologi termasuk di dalamnya adalah temperatur udara, arah dan kecepatan angin, radiasi matahari, tekanan atmosfer, penutupan awan, dan presipitasi. Data masukan aliran sungai meliputi meliputi volume harian, temperatur, salinitas, dan konsentrasi nutrien. Tahapan penggunaan model simulasi ELCOM-CAEDYM disajikan dalam Gambar 13.
67
Gambar 13
Skema tahapan penggunaan model simulasi ELCOM-CAEDYM (diadaptasi dan digambar ulang dari Hipsey et al., 2009).
3.9 Analisa Kecocokan dan Sensitivitas Model Untuk dapat menyatakan hasil simulasi berhasil baik atau dapat diterima dilakukan dengan membandingkan data lapangan dan data hasil simulasi dengan melakukan uji kecocokan dengan pendekatan normalisasi kesalahan mutlak dirata-ratakan (NMAE) (Alewell and Manderscheid, 1988; Bruce et al., 2006) terhadap nilai rata rata bulanan pada lapisan permukaan dengan persamaan : (168)
68
dimana, st = nilai simulasi pada waktu ke-t; ot = nilai pengamatan pada waktu ket; ō = nilai rata rata pengamatan selama periode simulasi; dan n = jumlah nilai hasil pengamatan. Metode NMAE mengukur simpangan absolut nilai simulasi dari nilai pengamatan, dinormalisasi terhadap rata-rata. Oleh karena kecocokan model tergantung dari sebaran data pengamatan, maka dihitung dari setiap variabel data lapangan. Kecocokan model dipertimbangkan dapat diterima ketika NMAE dari simulasi bernilai sama atau mendekati salah satu standar deviasi dari rata rata nilai data pengamatan. Analisis sensitivitas dilakukan pada parameter dalam model CAEDYM yang diasumsikan memiliki peran penting dalam perubahan variabel. Terdapat banyak metode untuk melakukan analisis sensitivitas yaitu mengkomputasi ulang dengan perubahan parameter dari nilai minimum kemudian dari nilai maksimum yang diperoleh dari referensi. Sumber referensi tidak selamanya memiliki nilai minimum atau maksimum, oleh karena iru bisa dilakukan dengan menambah atau mengurangi dengan setengah nilai parameter awal (Fasham et al., 1990). Analisis sensitivitas dalam penelitian ini dilakukan dengan menambahkan nilai masingmasing 10% dari nilai standar. Nilai 10% didasarkan pada perbandingan rata-rata selisih antara hasil observasi dan hasil model selurh variable, sehingga diasumsikan dengan penambahan nilai 10% dari parameter standar diharapkan akan memberikan hasil model mendekati hasil observasi. Koefisien sensitivitas (Sij) untuk menduga sensitivitas relatif dari variabel ke-i terhadap parameter ke-j dihitung dengan persamaan (Fasham et al., 1990; Chen et al., 1999; Bruce et al., 2006): (169) dimana
adalah koefisien sensitivitas,
terhadap nilai standar,
adalah perubahan variabel ke-i
adalah nilai standar variabel ke-i,
parameter ke-j terhadap parameter standar, dan ke-j.
adalah perubahan
adalah nilai standar parameter
