Luận án Triển khai mô hình Delft 3D mô phỏng chất lượng nước hệ đầm phá Tam Giang - Cầu Hai và tính toán sức tải môi trường - Cao Thị Thu Trang

/tấn) BOD COD SS TDS Kg/tấn sản phẩm Rayon 42 30 52 55 100 Acetate 75 45 78 40 100 Nylon 125 45 78 30 100 Acrylic 210 125 216 87 100 Polyeste 100 185 320 95 150 Nguồn: [132]. 50 c. Nguồn thải chăn nuôi Thải lượng do chăn nuôi được tính dựa trên tổng đàn gia súc, gia cầm hàng năm và suất phát thải đơn vị cho các loại gia súc, gia cầm (bảng 2.4). Bảng 2.4. Hệ số phát thải do chăn nuôi (kg/con/năm) STT Thông số Gia cầm Trâu, bò Lợn 1 COD 2,57 262,4 52,64 2 BOD 1,61 164 32,9 3 N-T 3,6 43,8 7,3 4 P-T 0,156 11,3 2,3 5 NO3+NO2* 0,036 0,438 0,073 6 NH4* 0,864 10,512 1,75 7 PO4* 0,0702 5,09 1,035 8 Thể tích nước thải (m3/năm) 21,5 8 14,6 Nguồn: WHO, 1993 [128]; (*) - Tính theo San Diego - McGlone, M.L., S.V. Smith and V. Nicolas, 2000 [130] d. Nguồn thải từ nuôi thuỷ sản Nguồn thải từ nuôi thuỷ sản được tính dựa trên hệ số phát thải và sản lượng nuôi các loại thuỷ sản hàng năm của vùng. Chất thải thuỷ sản chủ yếu là các chất hữu cơ và dinh dưỡng. Lượng thải phát sinh tuỳ thuộc vào hình thức và đối tượng nuôi, trong đó nuôi tôm công nghiệp và nuôi cá lồng có lượng phát thải đáng kể nhất. Bảng 2.5. Hệ số phát thải từ nuôi thuỷ sản Các chất ô nhiễm Hệ số phát thải (kg/tấn/năm) Nuôi TC tôm sú Nuôi cá lồng CODc 28,4 15,9 BODc 8,1 4,5 N-T 5,2a 2,9b P-T 4,7a 2,6b NO3+NO2 c 0,05 0,03 NH4 c 1,25 0,70 PO4 c 2,12 1,17 Nguồn: (a)- [133]; (b)- [134]; (c) - [130] e. Nguồn ô nhiễm do rửa trôi đất 51 Tải lượng ô nhiễm do rửa trôi đất được tính dựa trên số liệu về diện tích sử dụng đất các loại, số ngày mưa trung bình năm trong khu vực và hệ số phát thải ô nhiễm do rửa trôi từ các kiểu sử dụng đất. Bảng 2.6. Hệ số phát thải ô nhiễm do rửa trôi đất (kg/km2/ngày mưa) [135] Thông số Đất rừng và đồng cỏ Đất nông nghiệp Đất trống Đất khu dân cư COD 20 28 26 42 BOD 14 18 16 38 N - T 10 36 32 20 P - T 4 8 6 12 2.2.2.2. Ước tính tải lượng ô nhiễm đưa vào khu vực đầm phá Tam Giang - Cầu Hai Tải lượng ô nhiễm đưa vào khu vực đầm phá của nhóm nguồn phát sinh ở khu vực ven bờ được tính dựa trên tình hình thực tế quá trình giảm thiểu chất thải trong khu vực hoặc quá trình xử lý chất thải của từng loại nguồn ô nhiễm. Ước tính tổng tải lượng ô nhiễm đưa vào đầm phá Tam Giang - Cầu Hai từ các nguồn khác nhau có thể sử dụng công thức sau: ∑Qij = ∑Qij phát sinh x Rij x (1 - Hij) (4) Trong đó: ∑Qij - Tổng thải lượng của chất i vào khu vực từ các nguồn j (4 nguồn) ∑Qij phát sinh - Tổng tải lượng ô nhiễm i phát sinh từ các nguồn j Rij - Tỷ lệ đưa nước thải vào sông tương ứng với i và j Hij - Hiệu suất xử lý tương ứng với i và j 2.2.2.3. Dự báo nguồn ô nhiễm Phương pháp dự báo nguồn ô nhiễm dựa trên phương hướng phát triển kinh tế - xã hội của khu vực đến năm 2020 và định hướng đến năm 2030. Cơ sở để dự báo là tốc độ phát triển kinh tế - xã hội của khu vực, tốc độ phát triển từng ngành và lĩnh vực cụ thể như phát triển dân số, phát triển du lịch - dịch vụ, công nghiệp, nuôi trồng thuỷ sản. 52 2.2.3. Phương pháp mô hình hóa Đánh giá STMT là một hướng tiếp cận có tính phòng ngừa, tính toán và dự báo cho tương lai. Mô hình là một công cụ hữu hiệu cho dự báo do có cơ sở khoa học và độ chính xác nhất định. Vì vậy, mô hình hóa là phương pháp chủ đạo, xuyên suốt trong nghiên cứu STMT của luận án. 2.2.3.1. Cơ sở lý thuyết của mô hình Delft -3D a) Mô hình thủy động lực [113] Phương trình cơ bản trong mô hình thuỷ động lực của Delft3d là phương trình Navier Stokes. Với giả thiết áp dụng phương trình Navier Stokes với chất lỏng không nén trong nước nông và phương pháp xấp xỉ Boussinesq. Sự biến đổi của thành phần vận tốc thẳng đứng trong phương trình động lượng được bỏ qua. Với mô hình 3 chiều, thành phần vận tốc thẳng đứng được tính toán từ phương trình liên tục. Phương trình liên tục Phương trình liên tục cho độ sâu trung bình (viết trong hệ toạ độ cong trực giao):     Q GVd GG GUd GGt                  )(1)(1 (5) Trong đó: , : là các tọa độ ngang trong hệ toạ độ cong trực giao;  GG , : là các hệ số chuyển đổi từ hệ toạ cong trực giao sang hệ tọa độ Đề các; d: là độ sâu tại điểm tính (độ sâu của nước dưới đường chuẩn (0 hải đồ)); : là mực nước tại điểm tính (so với 0 Hải đồ); U, V: lần lượt là các thành phần vận tốc trung bình theo các hướng , ; qin và qout : lần lượt là nguồn nước đưa vào và ra trên 1 đơn vị thể tích; H: là độ dày cột nước tại điểm tính (H = d + ); P, E: lần lượt là lượng mưa và bốc hơi. Trong phương trình (5): t  - biểu diễn sự thay đổi mực nước theo thời gian   53         GUd GG )(1 - biểu diễn sự thay đổi của vận tốc dòng chảy trung bình U theo hướng          GVd GG )(1 - biểu diễn sự thay đổi của vận tốc dòng chảy trung bình V theo hướng  Với Q thể hiện lượng thêm vào hay mất của nguồn nước, sự bốc hơi và mưa trên 1 đơn vị diện tích: EPdqqHQ outin     0 1 )( (6) b) Mô hình chất lượng nước [112]  Phương trình cơ bản của mô hình chất lượng nước Cơ sở toán học của mô hình lan truyền ô nhiễm là phương trình lan truyền và khuyếch tán vật chất: ),( tCFCu z C D z Cu y C D y Cu x C D xt C zzyyxx                                     (7) Trong các phương trình trên: Dx, Dy, Dz là các hệ số khuyếch tán theo các phương x, y, z F(C, t) là nguồn thêm vào hoặc mất đi. C: hàm lượng vật chất  Các quá trình chính xảy ra được mô phỏng bởi Delft-3D-WAQ liên quan đến các hợp chất C, N, P - Sự thoáng khí: Sự trao đổi oxy giữa khí quyển và bề mặt nước là một hàm tuyến tính về sự khác nhau giữa nồng độ oxy thực tế và bão hòa theo phương trình: Rrear=klrear x (Coxs-Cox)/H (8) klrear= ( a xvb Hc ) +(d x W2) (9) Coxs=f(T, Ccl)=f(T, SAL) (10) Trong đó: 54 a, b, c, d Các hệ số với các giá trị khác nhau của 11 lựa chọn thông khí Ccl Nồng độ clorua (gCl/m3) Cox Nồng độ oxy hòa tan thực tế (gO2/m3) Coxs Nồng độ oxy hòa tan bão hòa (gO2/m3) H Chiều cao cột nước (m) klrear Hệ số chuyển đổi thoáng khí (hệ số cấp khí) (1/ngày) Rrear Tốc độ thoáng khí (gO2/m3/ngày) SAL Độ muối (kg/m3) T Nhiệt độ nước (oC) V Vận tốc dòng (m/s) W Tốc độ gió ở độ cao 10 m (m/s) - Sản xuất sơ cấp: Sản xuất sơ cấp của TVPD trong Delft-3D - DYNAMO tuân theo các hàm sau: Rgpi=fnuti x flti x kgpi x Calgi (11) kppi=achfi x kgpi10x ftmpi (12) ftmpi=ktpgiT-10 (13) Calg= max (Calgi, Calgmin) (14) Trong đó: Calg Nồng độ sinh khối tảo (gC/m3) achf Nồng độ chlorophyll trong tảo (gChlof/gC) Calgmin Nồng độ sinh khối tảo thấp nhất, một giá trị ngưỡng (gC/m3) flt Hệ số giới hạn ánh sáng (-) fnut Hệ số giới hạn dinh dưỡng theo Monod (-) Rgp Tốc độ sản xuất sơ cấp (gC/m3/ngày) i Chỉ số cho nhóm tảo (tảo lục hoặc giáp (-)) T Nhiệt độ (oC) ftmp Hệ số giới hạn nhiệt độ đối với sự sản xuất (-) kgpi Tốc độ sản suất sơ cấp tiềm tàng (1/ngày) kgpi10 Tốc độ sản suất sơ cấp tiềm tàng ở 10oC (gC/gChlor/ngày) ktgp Hệ số nhiệt độ đối với sản xuất sơ cấp (-) Hệ số giới hạn dinh dưỡng trong DYNAMO chỉ xem xét đến dinh dưỡng bị giới hạn nhất: fnuti= min ( Cn Ksn+Cn , Cph Ksph+Cph , CSi KsSi+Csi ) (15) Trong đó: 55 Cn Nồng độ amoni + nitrat (gN/m3) Cph Nồng độ phosphat (gP/m3) CSi Nồng độ silicate vô cơ hòa tan (gSi/m3) Ksn Hằng số bán bão hòa cho dinh dưỡng nitơ (gN/m3) Ksph Hằng số bán bão hòa cho phôt phat (gP/m3) KsSi Hằng số bán bão hòa cho silicat (gP/m3) fnuti= fami + (1 – fami) x fnii (16) fami= Cam x Cph x Csi (Ksami xCph x Csi+Cam x Ksphi x Csi+Cam x Cph x Kssii+Cam x Cph xCsi) (17) fnti= Cni x Cph x Csi (Ksni x Cph x Csi+Cni x Ksphi x Csi+Cni x Cph x KsSii+Cni x Cph xCsi) (18) fam Hệ số giới hạn dinh dưỡng đặc trưng cho amoni (-) fni Hệ số giới hạn dinh dưỡng đặc trưng cho nitrat (-) Cam Nồng độ amoni (gN/m3) Cni Nồng độ nitrat (gN/m3) Cph Nồng độ phosphat (gP/m3) Csi Nồng độ silicate vô cơ hòa tan (gSi/m3) Ksam Hằng số bán bão hòa cho amoni (gN/m3) Ksni Hằng số bán bão hòa cho nitrat (gN/m3) Ksph Hằng số bán bão hòa cho phôt phat (gP/m3) KsSi Hằng số bán bão hòa cho silicat (gP/m3) Tốc độ hô hấp và tỷ lệ chết trong DYNAMO được mô tả theo: Rrspi=krspi x Calgi+frspi x Rgpi (19) Rmrti=kmrti x Calgi (20) kmrti=f(Ccl, T) (21) Trong đó: Ccl Nồng độ muối clorua (g/m3) frsp Phần sản xuất hô hấp (-) kmrt Tốc độ quá trình chết tổng cộng (1/ngày) krsp Tốc độ hô hấp duy trì (1/ngày) Rmrt Tốc độ chết tổng số (gC/m3/ngày) Rrsp Tốc độ hô hấp tổng (gC/m3/ngày) - Sự nitrat hóa: Amoni là đối tượng của quá trình nitrat hóa tuân theo phương trình hóa học: 56 NH4+ + 2O2  NO3- + 2H+ +H2O Quá trình nitrat hóa được xây dựng dưới dạng tổng của một tốc độ nền không đổi (giới hạn bậc 0) và tốc độ phụ thuộc vào nồng độ. Tốc độ này phụ thuộc vào lượng oxy hòa tan và amoni. Có hai lựa chọn thay thế cho sự phụ thuộc này. R=k0+k x fam x fox (22) fam= ( Cam Ks+Cam ) fox-( Cox (Ksox+Cox) ) Lựa chọn 1 (23) fam = Cam ^ fox -(1 - foxmin) x ( Cox-Coxc Coxo-Coxc ) 10a + foxmin Lựa chọn 2 (24) Trong đó R Tốc độ nitrat hóa (gN/m3/ngày) k0 Tốc độ nền không đổi (gN/m3/ngày) K hằng số động học (gN/m3/ngày) fam hệ số động học cho amoni (-) hoặc (gN/m3) fox hệ số động học cho oxy hòa tan (-) Cam Nồng độ amoni hòa tan (gN/m3) Ks Hằng số bán bão hòa (gN/m3) Cox Nồng độ oxy hòa tan (gO2/m3) Coxc Nồng độ oxy hòa tan tới hạn (gO2/m3) Coxo Nồng độ oxy hòa tan tối ưu (gO2/m3) Ksox Hằng số bán bão hòa cho oxy (gO2/m3) foxmin giá trị nhỏ nhất của chức năng giới hạn oxy (-) a hệ số độ cong (-) - Sự đề nitrat hóa: Nitrat là đối tượng của quá trình khử nitơ. Quá trình đề nitrat hóa được xây dựng dưới dạng tổng của một tốc độ nền không đổi (giới hạn bậc 0) và tốc độ phụ thuộc vào nồng độ. Tốc độ này phụ thuộc vào lượng oxy hòa tan và nitrat. Có hai lựa chọn thay thế cho sự phụ thuộc này. R=k0+k x fni x fox (25) fni= ( Cni Ks+Cni ) ^fox - 1.0- ( Cox Ksox+Cox ) Lựa chọn 1 (26) fni=Cni ^ fox- ( Coxc-Cox Coxc-Coxo ) Lựa chọn 2 (27) Trong đó R Tốc độ đề nitrat hóa (gN/m3/ngày) k0 Tốc độ nền không đổi (gN/m3/ngày) 57 k hằng số động học (gN/m3/ngày) fni hệ số động học cho amoni (-) hoặc (gN/m3) fox hệ số động học cho oxy hòa tan (-) Cni Nồng độ amoni hòa tan (gN/m3) Ks Hằng số bán bão hòa (gN/m3) Cox Nồng độ oxy hòa tan (gO2/m3) Ksox Hằng số bán bão hòa cho oxy (gO2/m3) foxmin giá trị nhỏ nhất của chức năng giới hạn oxy (-) Coxc Nồng độ oxy hòa tan tới hạn (gO2/m3) Coxo Nồng độ oxy hòa tan tối ưu (gO2/m3) - Sự khoáng hóa trong cột nước Tốc độ khoáng hóa của các thành phần BOD/COD phụ thuộc vào tuổi vật chất hữu cơ, với điều kiện là BOD và COD được mô phỏng đồng thời. Bằng cách mặc định hàm tuổi (fage) là 1,0. Tốc độ khoáng hóa được mô tả theo các phương trình động học bậc 1: Rmin=Fage x kmin x Cx (28) fage=fagemin+(fagemax-fagemin)x ef(AI) (29) AI= COD BOD5 (30) Trong đó AI Chỉ số tuổi của vật chất hữu cơ (-) Cx Giá trị của BOD/COD (gO2/m3) BOD5 Giá trị của BOD5 (gO2/m3) COD Giá trị tổng của COD (gO2/m3) fage Hàm suy giảm theo tuổi (-) kmin Hằngsố động học bậc 1 (1/ngày) Rmin Tốc độ khoáng hóa của BOD/COD (gO2/m3/ngày) Nhu cầu oxy trong trầm tích (SOD) được tính dưới dạng tổng của các quá trình bậc 0 và bậc 1, cho phép xác định nhu cầu oxy trầm tích cơ bản. Hai lựa chọn sẵn có đối với tốc độ khoáng hóa, một trong số đó tính đến sự thoát khí metan vào khí quyển. Tốc độ khoáng hóa theo lựa chọn 2 được tính theo: Rsod=(1-fch4)x ( k0sod H + ksod x SOD V ) (31) Trong đó SOD Khối lượng nhu cầu oxy trầm tích (gO2) fch4 phân đoạn vật chất hữu cơ thoát ra dưới dạng metan (-) 58 k0sod Nhu cầu oxy trầm tích cơ bản ổn định (gO2/m2/ngày) H Chiều sâu của cột nước (m) Rsod Nhu cầu oxy trầm tích (gO2/m 2/ngày) V Thể tích của khối nước (m3) - Sự khoáng hóa của các chất hữu cơ vẩn đục và chất hữu cơ khác ((DetC /OOC) Sự khoáng hóa của mảnh vụn khi sử dụng DetC/OOC (DetN/OON vv,) lựa chọn cho vật chất hữu cơ được xây dựng dựa trên giới hạn bậc 0 cộng bậc 1. Thông thường, giới hạn bậc 0 được thiết lập bằng 0 gC/m3/ngày. Tốc độ khoáng hóa bậc 1 phụ thuộc vào nhiệt độ nước như được mô tả trong phương trình 34. Rmini = k0 min + kmin x Cx (32) Trong đó Cx Nồng độ C hữu cơ (gC/m3 b) kmin Tốc độ khoáng hóa bậc 1 (1/ngày) k0min Tốc độ khoáng hóa bậc 0 (1/ngày) Rmin Tốc độ khoáng hóa đối với C, N hoặc P (gC/N/P/Si/m3b/ngày) i Chỉ số cho thành phần chất hữu cơ Công thức cho nitơ hữu cơ (DetN hoặc OON), phospho hữu cơ (Det P hoặc OOP), hoặc silic “hữu cơ” (DetSi hoặc OOSi) là giống hệt nhau. Tốc độ khoáng hóa bậc 1 có thể có giá trị tối đa hoặc tối thiểu phụ thuộc vào thành phần (tỷ lệ của vật chất hữu cơ). - Sự hấp phụ của phosphat Sự hấp phụ của phosphats được mô tả bởi công thức: Kads= Cpha x OHa Cphd x Cads (33) Cadst =fcor x ∑ 3 i=1 (ffei x Cimi)x 1 56x103x φ (34) Rsorp=ksorp x (Cphae - Cpha') (35) Với a Hằng số phản ứng thành phần (tỷ lệ nguyên tố), có giá trị = 0 với lựa chọn 1 Cads Nồng độ chất hấp phụ tự do (molFe/l nước) Cadst Tổng nồng độ mol của các vị trí hấp phụ (molFe/l nước) Cimi Nồng độ của phân đoạn vật chất vô cơ i = 1,2,3 (g D.W (trọng lượng 59 khô/m3 khối lượng (bulk) Cpha Nồng độ phosphat được hấp thụ (molP/l nước) Cpha’ Nồng độ phosphat được hấp thụ sau một bước thời gian (molP/m3 khối lượng) Cphae Nồng độ phosphat được hấp thụ ở trạng thái cân bằng (molP/l nước) Cphd Nồng độ phosphat hòa tan (molP/l nước) fcor Hệ số hiệu chuẩn của phần sắt (III) bị oxi hóa, có giá trị bằng 1 với lựa chọn 1 (-) ffei Phần Fe (III) hoạt hóa trong phần vật chất vô cơ i = 1, 2, 3 (gFe/g trọng lượng khô) ksorp Tốc độ phản ứng hấp phụ (1/ngày) Kads Hằng số cân bằng hấp phụ (mol (a-1)/ l(a-1) OH Hằng số cân bằng hấp phụ (mol/l nước) Rsorp Tốc độ hấp phụ hoặc giải hấp (gP/m3/ngày) φ Trạng thái lỗ rỗng (-) Việc giải phóng chậm phosphat giải hấp (PAP) được mô tả bằng tổng của các hằng số động học bậc 0 và bậc 1 theo: Rdes=k0+kl x Cphp (36) Trong đó: Cphp Nồng độ của phosphat giải hấp chậm (gP/m3) k0 Tốc độ nền không đổi (gP/m3/ngày) kl Hằng số động học bậc 1 (1/ngày) Rdes Tốc độ giải hấp phosphat chậm ((gP/m3/ngày) Tốc độ bậc 1 kl được thiết lập bằng 0 khi nhiệt độ giảm dưới nhiệt độ thấp tới hạn. - Sự hòa tan của opal silicat Opal silicat (OPAL) hòa tan chậm theo phương trình động học bậc 2 (giả định) sau: Rsol=ksol x Csip x (Cside - Csid φ ) (37) Trong đó Csid Nồng độ silicat hòa tan (gSi/m3) Cside Nồng độ silicat hòa tan cân bằng (gSi/m3) Csip Nồng độ silicat opal (gSi/m3) ksol Hằng số động học bậc hai (giả định) (m3/gSi/ngày) Rsol Tốc độ hòa tan silicat opal (gSi/m3/ngày) 60 φ Trạng thái lỗ rỗng (-) 2.2.3.2. Triển khai mô hình Delft 3D mô phỏng chất lượng nước hệ đầm phá Tam Giang - Cầu Hai và tính toán sức tải môi trường a) Thiết lập các dữ liệu đầu vào của mô hình  Địa hình Số liệu độ sâu và đường bờ của khu vực nghiên cứu (hình 2.3) được số hoá từ các bản đồ địa hình UTM tỷ lệ 1: 50.000 do Cục Đo đạc Bản đồ xuất bản, đây là những bản đồ mới xuất bản năm 2005 với hệ tọa độ nhà nước VN-2000. Những bản đồ đó đã được quét vào máy tính, số hoá và xử lý bằng các phần mềm AcrGis, MapInfo. Ngoài ra, địa hình vùng biển Thừa Thiên Huế và các vùng lân cận khi thiết lập mô hình dùng để NESTHD còn được tham khảo và bổ sung từ cơ sở dữ liệu địa hình ETOPO5 (Earth Topography - 5 Minute) của Trung tâm Tư liệu Địa vật lí Quốc gia Mỹ NGDC (National Geophysical Data Center) và GEBCO -1 (General Bathymetric Chart of the Ocean (GEBCO) one minute) của Trung tâm tư liệu hải dương học vương quốc Anh (British Oceanographic Data Centre-BODC). Hình 2.3.Trường độ sâu của mô hình thủy động lực 61  Số liệu khí tượng Những số liệu về khí tượng sử dụng trong mô hình thuỷ động lực bao gồm bức xạ mặt trời, độ ẩm tương đối, lượng mây, nhiệt độ không khí có giá trị trung bình theo mùa (mùa khô-tháng 5 và mùa mưa-tháng 11). Số liệu gió của tháng 11 năm 2011 và tháng 5 năm 2012 được quan trắc liên tục 6h/ốp ghi nhận tại Trạm Huế làm dữ liệu đầu vào cho mô hình tính toán. Các số liệu khí tượng thu thập từ trạm khí tượng Huế từ năm 2005 đến 2010 + Nhiệt độ không khí trung bình theo tuần + Độ ẩm không khí trung bình theo tuần + Số giờ nắng trong tuần + Lượng mưa theo tuần (số giờ mưa, số ngày mưa)  Số liệu thủy văn Các đặc trưng về dòng chảy của các sông (bảng 2.7) được sử dụng để thiết lập mô hình thủy động lực. Bảng 2.7. Những đặc trưng dòng chảy của các sông đổ vào hệ đầm phá TG – CH [59] Sông Chiều dài (km) Diện tích lưu vực (km2) Dòng chảy (l/s.km2) Lưu lượng bình quân (m/s) Tổng lượng nước (tỷ m3/năm) Hương Tả Trạch 54 717 70,2 50,3 1,580 Hữu Trạch 51 729 66,0 37,6 1,180 Bồ 64 720 63,0 45,3 1,420 Ô Lâu 66 900 57,8 17,3 0,545 Truồi 24 149 66,7 7,3 0,229 Nong 20 99 54,2 3,4 0,112 Cầu Hai 67,0 2,7 0,105 Số liệu dòng chảy đo đạc tại 3 trạm liên tục tại 3 điểm là phá Tam Giang, Đầm Sam, Đầm Cầu Hai mùa mưa 2011 và mùa khô 2012 đã được sử dụng để so sánh, hiệu chỉnh mô hình. Số liệu về dao động mực nước được cung cấp để hiệu chỉnh các tham số tính toán và làm số liệu đầu vào tại các biên mở ở phía biển của mô hình thuỷ động lực. 62 Dao động mực nước tại các trạm thủy văn thuộc tỉnh Thừa Thiên - Huế được dùng để hiệu chỉnh mô hình. Số liệu thuỷ triều tại các biên mở phía biển là những hằng số điều hoà thuỷ triều. Các hằng số điều hoà thuỷ triều được tính toán từ chuỗi số liệu quan trắc mực nước trong khoảng thời gian dài trong một số đề tài, dự án của Trung tâm Khí tượng Thuỷ văn biển, Viện Tài nguyên và Môi trường biển, Trường Đại học Thủy lợi. Những điểm biên lỏng phía biển không có số liệu quan trắc thì số liệu dùng cho mô hình là kết quả tính toán của mô hình NAO-TIDE [136] - chương trình dự báo thuỷ triều của Đài Thiên văn Quốc gia Nhật Bản NAO.  Lựa chọn miền tính, lưới tính Kích thước miền tính có phạm vi 70 km theo hướng bắc nam và 1 - 10 km theo hướng đông tây. Hệ thống đầm phá TG-CH có đường bờ khúc khuỷu nên hệ thống lưới cong trực giao được chọn sử dụng để tính toán cho khu vực. Toàn bộ khu vực tính toán bao gồm 95 x 513 ô lưới (hình 2.4). Miền tính có các biên lỏng phía biển là: cửa Thuận An và cửa Tư Hiền; ngoài ra còn có các biên lỏng sông: sông Ô Lâu, sông Hương và sông Truồi. Hình 2.4. Hình thái và phạm vi lưới tính của mô hình  Thời gian tính toán 63 Để phù hợp cho việc lựa chọn thông số đầu vào và hiệu chỉnh mô hình tính toán, mô hình được chạy theo hai mùa: mùa khô và mùa mưa. Mùa mưa bắt đầu tính toán từ 0h00, ngày 01/10/2011 và kết thúc lúc 23h00, ngày 30/11/2011. Mùa khô bắt đầu tính toán từ 0h00, ngày 01/4/2012 và kết thúc lúc 23h00, ngày 31/5/2012. Bước thời gian để chạy mô hình là 30 s. Các quá trình vật lý được lựa chọn trong mô hình: - Lực Coriolis. Vĩ độ trung bình được lấy là 16030’ độ vĩ bắc. - Trao đổi nhiệt với khí quyển qua lớp nước bề mặt: diện tích bề mặt khu vực tính (nơi diễn ra sự tương tác biển - khí quyển), lượng mây trên bầu trời, độ ẩm tương đối, nhiệt độ không khí, bức xạ mặt trời, bốc hơi, lượng mưa. - Ma sát đáy, sử dụng theo công thức tính toán hệ số Strickler. - Nhiệt độ nước biển và độ muối, sử dụng giá trị đặc trưng trung bình theo mùa nhiều năm tại các điểm trên biên lỏng.  Điều kiện ban đầu và điều kiện biên Mô hình thủy động lực - Tại thời điểm ban đầu: + Giá trị mực nước lấy bằng 0m + Độ muối lấy trung bình thấp nhất + Nhiệt độ nước biển cũng được lấy trung bình thấp nhất - Điều kiện biên: + Tại các biên mở phía biển (Thuận An, Tư Hiền): dùng các hằng số điều hoà thuỷ triều đã được tính toán và nội suy phù hợp với điều kiện địa phương. + Tại biên sông: sử dụng số liệu lưu lượng trung bình nhiều năm. + Tại tất cả các biên lỏng đều sử dụng và tham khảo các kết quả quan trắc nhiệt muối trong các đợt khảo sát năm 2011 và 2012 tại khu vực hệ đầm phá TG - CH. Bảng 2.8. Điều kiện biên của mô hình thủy động lực Thông số Sông Ô Lâu Sông Hương Sông Truồi Cửa Thuận An Cửa Tư Hiền Mùa mưa Nhiệt độ, oC 20 20 20 20 20 Độ muối, ‰ 0,2 0,2 1,4 15 15 64 Lưu lượng (m3/s) 27,79 72,3 11,8 Mùa khô Nhiệt độ, oC 25 25 25 25 25 Độ muối, ‰ 3 3 3 23 25 Lưu lượng (m3/s) 6,92 18 2,92 Mô hình chất lượng nước - Điều kiện ban đầu: Tại thời điểm ban đầu, các giá trị được thiết lập bằng 0. Sau mỗi lần chạy, điều kiện ban đầu được lấy là giá trị cuối cùng của file chạy trước đó - Điều kiện biên: có 3 biên sông và 2 biên phía biển. Lưu lượng của các biên sông và biên biển được chạy trên nền thủy động lực. Do các quá trình xảy ra trong thủy vực khá phức tạp, ngoài các quá trình vật lý như trao đổi nước, lắng đọng, khuyếch tán còn có các quá trình sinh học. Do đó, để mô phỏng chất lượng nước cũng cần mô phỏng các yếu tố sinh học như chlorrohyll a và các loại tảo lục, tảo giáp. - Dữ liệu nguồn thải: dữ liệu nguồn thải của mô hình chất lượng nước được tính toán và phân bổ chi tiết cho toàn vùng (xem mục b) về phương pháp tính toán và phân bổ điểm thải) - Các kịch bản mô phỏng: kịch bản cơ sở (thời điểm năm 2011 -2012) và các kịch bản dự báo cho các năm 2020, 2030, các kịch bản đột xuất được xây dựng trên cơ sở tình hình phát triển thực tế và các quy hoạch của các ngành, lĩnh vực (xem mục d) về cơ sở xây dựng các kịch bản) b) Thiết lập các dữ liệu nguồn thải Dựa trên đặc điểm lưu vực và kinh tế- xã hội, đã thiết kế 74 điểm nguồn thải đổ vào khu vực đầm phá TG - CH, được chỉ ra trong hình 2.5. 65 Hình 2.5. Vị trí các điểm thải khu vực TG - CH phục vụ chạy mô hình Lượng chất thải đưa vào hệ đầm phá TG - CH được thiết kế và phân bổ trên cơ sở kết quả tính toán tải lượng thải thời điểm 2011 - 2012 và thời điểm 2020, 2030. Mùa mưa được tính từ tháng 9 đến tháng 12 với khoảng 125 ngày, mùa khô được tính từ tháng 1 đến tháng 8 với khoảng 240 ngày. c) Hiệu chỉnh và đánh giá tính tương hợp của mô hình; thẩm định mô hình  Hiệu chỉnh mô hình thủy đông lực Bên cạnh các thông tin đầu vào của mô hình như nhiệt độ, độ muối, lưu lượng sông, địa hình, v.v. của khu vực thì hiệu chỉnh mô hình thủy động lực là một bước quan trọng ảnh hưởng đến các kết quả mô phỏng sau này. Kết quả tính toán từ mô hình được so sánh với số liệu quan trắc để đánh giá sự sai khác giữa tính toán và thực tế, từ đó kiểm tra, hiệu chỉnh các số liệu, các tham số nhằm mục đích tối ưu hoá quá trình tính toán và giảm sai khác đến một mức nhỏ nhất có thể. Tại khu vực đầm phá TG - CH, mô hình thủy động lực được hiệu chỉnh dựa trên số liệu quan trắc về dòng chảy trong hai mùa mưa và mùa khô năm 2011 -2012 và số liệu của mô hình. Với phần mềm Delft -3D tại khu vực, có 3 tham số chính ảnh hưởng đến kết quả mô 66 phỏng thủy động lực là hệ số nhám Manning, hệ số nhớt theo phương ngang và hệ số khuyếch tán theo phương ngang. Các hệ số này đều được gọi là các tham số hiệu chuẩn; giá trị của chúng phải được xác định trong quá trình hiệu chuẩn [112]. Cơ sở để hiệu chỉnh dựa vào sai số giữa mô hình và quan trắc, giá trị này càng nhỏ thì kết quả mô phỏng càng gần với thực tế. Sai số bình phương trung bình (RMSE- Root Mean Square Error) là một trong những chỉ tiêu quan trọng để đánh giá độ chính xác của mô hình.   N OP RMSE N i ii    1 2 Trong đó: i = 1, n là số lần quan trắc được thực hiện; Pi là giá trị dự báo của mô hình tại thời điểm i; Oi là giá trị quan trắc tại thời điểm i. Một chỉ số quan trọng để đánh giá độ tương hợp của mô hình là chỉ số Nash - Sutcliffe (1970) (NSE).         2 0,0 2 ,,0 1 XX XX NSE i iSi X0,i: Giá trị thực đo XS,i: Giá trị tính toán hoặc mô phỏng Xo : Giá trị thực đo trung bình Đánh giá độ tương hợp của mô hình qua chỉ số NSE như sau [137]: 0,75 <NSE<1 Rất tốt 0,65 <NSE<0,75 Tốt 0,5< NSE <0,65 Thỏa mãn NSE<0,5 Không thỏa mãn  Hiệu chỉnh mô hình chất lượng nước Trong thủy vực, có nhiều quá trình liên quan ảnh hưởng đến chất lượng nước. Việc hiệu chỉnh mô hình chất lượng nước là một công việc khó, liên quan đến rất nhiều tham số trong mô hình. Trong phần mềm mô phỏng chất lượng nước của Delft 3 D, có tới 85 tham số ảnh hư