ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
……………………

Phan Thành Bắc

MÔ PHỎNG QUÁ TRÌNH LAN TRUYỀN VẬT CHẤT Ô
NHIỄM DƯỚI TÁC ĐỘNG CỦA CÁC YẾU TỐ ĐỘNG LỰC
TẠI VỊNH CAM RANH BẰNG MÔ HÌNH SỐ

LUẬN VĂN THẠC SĨ KHOA HỌC

Hà Nội - 2012


ĐẠI HỌC QUỐC GIA HÀ NỘI
TRƯỜNG ĐẠI HỌC KHOA HỌC TỰ NHIÊN
....................................

Phan Thành Bắc

MÔ PHỎNG QUÁ TRÌNH LAN TRUYỀN VẬT CHẤT Ô
NHIỄM DƯỚI TÁC ĐỘNG CỦA CÁC YẾU TỐ ĐỘNG LỰC
TẠI VỊNH CAM RANH BẰNG MÔ HÌNH SỐ
Chuyên ngành: Hải dương học
Mã số: 60.44.97

LUẬN VĂN THẠC SĨ KHOA HỌC

Người hướng dẫn khoa học: PGS. TS Nguyễn Minh Huấn


MỤC LỤC

MỞ ĐẦU..............................................................................................................................................1
CHƯƠNG 1. MÔ HÌNH SỐ TRỊ......................................................................................................4
1.1 TỔNG QUAN TÌNH HÌNH NGHIÊN CỨU

4

1.1.1 Tổng quan tình hình nghiên cứu trên thế giới..............................................................4
1.1.2 Tổng quan tình hình nghiên cứu trong nước................................................................6
1.2 MIKE 21 HD

8

1.2.1 Cơ sở toán học..............................................................................................................8
1.2.2 Phương pháp số..........................................................................................................12
1.3 MÔĐUN ECOLAB

16

1.3.1 Cơ sở lý thuyết............................................................................................................16
1.3.2 Ôxy hòa tan (DO) và nhu cầu ôxy sinh hóa (BOD)...................................................18
1.3.3 Các hợp phần của Nitơ...............................................................................................21
1.3.4 Hợp phần của Photpho...............................................................................................24
CHƯƠNG 2. TỔNG QUAN VÙNG NGHIÊN CỨU...................................................................25
2.1 TỔNG QUAN VỀ ĐIỀU KIỆN TỰ NHIÊN

25

2.1.1 Vị trí địa lí...................................................................................................................25

3.3 MỘT SỐ KẾT QUẢ TÍNH TOÁN

44

3.3.1 Kết quả tính toán cho mùa khô...................................................................................44
3.3.2 Kết quả tính toán cho mùa mưa..................................................................................74
3.4 ĐÁNH GIÁ TÁC ĐỘNG CỦA CÁC YẾU TỐ Ô NHIỄM

106

KẾT LUẬN

122

KIẾN NGHỊ

124

TÀI LIỆU THAM KHẢO

125


MỞ ĐẦU
Trong những năm gần đây, khu vực đầm thuỷ triều đang đứng trước nguy
cơ ô nhiễm nguồn nước. Đầm Thủy Triều nằm trong vịnh Cam Ranh, thuộc địa
bàn huyện Cam Lâm và thành phố Cam Ranh, tỉnh Khánh Hòa. Nơi đây phong phú
và đa dạng về số lượng cũng như trữ lượng thủy sản. Trong tương lai, đầm Thủy
Triều còn là mắt xích quan trọng trong việc phát triển du lịch của tỉnh Khánh Hòa
khi vịnh Cam Ranh đã được tỉnh này quy hoạch thành trung tâm du lịch biển tầm cỡ

chiến lược để quy hoạch, khai thác một cách hiệu quả tài nguyên khu vực đầm cũng
như việc kiểm soát và điều tiết các nguồn thải hợp lý hơn.
Nhận thức được mức độ cấp thiết của vấn đề môi trường vịnh Cam Ranh,
học viên lựa chọn hướng nghiên cứu với đề tài: “Mô phỏng quá trình lan truyền
vật chất ô nhiễm dưới tác động của các yếu tố động lực tại vịnh Cam Ranh bằng
mô hình số” để có thể mô phỏng một số vật chất có khả năng ảnh hưởng đến chất
lượng môi trường. Có nhiều kỹ thuật đánh giá mức độ ô nhiễm nước dựa vào giá trị
của các thông số chọn lọc. Các kỹ thuật này sử dụng các chỉ số để thực hiện mức độ
ô nhiễm. Trong đó có thể nêu một số chỉ số đang được công nhận như: Chỉ số ô
nhiễm dinh dưỡng (NPI) dựa vào các thông số NH 4+, NO3-, NO2-, tổng P, pH,
chlorophyll, độ dẫn điện và độ đục. Chỉ số ô nhiễm hữu cơ (OPI) dựa vào các thông
số BOD, COD, nhiệt độ và DO. Với nguồn số liệu có được từ một số đề tài được
thực hiện tại Viện Hải dương học như đề tài cấp Cơ sở phòng Vật lý biển, phòng
Thủy địa hóa, đề tài cấp Viện Khoa học và Công nghệ, Các Dự án hợp tác quốc tế,
tác giả sử dụng gói phần mềm MIKE 21 HD, ECO Lab để mô phỏng quá trình lan
truyền một số vật chất có thể gây ô nhiễm từ các nguồn thải của khu công nghiệp,
nuôi trồng thủy sản và khu dân cư trong 2 mùa: mùa mưa và mùa khô. Trong khuôn
khổ của luận văn, mục tiêu của học viên là có thể tính toán, mô phỏng, đưa ra được
bức tranh về quá trình động lực và quá trình truyền tải các vật chất gồm BOD, DO,
NO3-, PO4+, NH3+. Một kịch bản mô phỏng sự lan truyền các vật chất ô nhiễm với
giả thiết có sự gia tăng cực đại nồng độ các chất gây ô nhiễm từ số liệu thực đo tại
cống xả thải và công suất tính tại thời điểm khảo sát từ các nguồn thải của khu công
nghiệp, nuôi trồng thủy sản và khu dân cư để có thể đánh giá mức độ lan truyền và
ảnh hưởng của các vật chất này tới chất lượng nước các bãi tắm khu vực Cam Ranh.
2


Các kết quả nghiên cứu trong luận văn góp phần bổ sung thêm các thông
tin khoa học về những nghiên cứu, đánh giá vai trò và sự tác động của các từ các
nguồn thải của khu công nghiệp, nuôi trồng thủy sản và khu dân cư tác động ngược

là công cụ hữu hiệu cho các nhà môi trường học, sinh học, những nhà mô hình hóa
chất lượng nước và bất kỳ ai cần quan tâm tới việc đánh giá rủi ro và suy giảm các
hệ sinh thái thủy sinh.

4


Mô hình QUAL2K (hay Q2K) (River and Stream Water Quality Model)
được nâng cấp từ mô hình trước đó là QUAL2E (hay Q2E (Brown và Barnwell
1987)). Đây là mô hình mô phỏng chất lượng nước suối và sông một chiều có sự
tham gia của quá trình xáo trộn rối và bên. Một đặc điểm linh hoạt của mô hình này
là có thể chạy được trong môi trường Visual basic hoặc trong môi trường Excel. Mô
hình có những đặc điểm sau: có thể tính toán trên từng phân đoạn của sông và các
nhánh sông. Mô hình tính toán chu trình Nitơ. Thông qua các chu trình chuyển hóa
nitơ để biểu diễn các hợp chất cacbon (loại ôxy hóa nhanh và chậm), các loại
cacbon hữu cơ không sống (các phân tử cacbon, nitơ, phôtpho trong các hợp chất
hóa học). Các quá trình thiếu hụt ôxy gần tới giá trị không do các quá trình ôxy hóa,
trong đó quá trình khử nitơ như là bước tương tác đầu tiên. Tính toán thông lượng
trao đổi ôxy hòa tan và các dinh dưỡng giữa trầm tích và nước.
DELFT 3D của Viện nghiên cứu thuỷ lực Hà Lan cho phép kết hợp giữa
mô hình thuỷ lực 3 chiều với mô hình chất lượng nước. Ưu điểm của mô hình này
là việc kết hợp giữa các module tính toán phức tạp để đưa ra những kết quả tính mô
phỏng cho nhiều chất và nhiều quá trình tham gia.
SMS của Trung tâm nghiên cứu và phát triển kỹ thuật của quân đội Mỹ
xây dựng cho phép kết hợp giữa mô hình thuỷ lực 1, 2 chiều với mô hình chất lượng
nước, trong đó module RMA4 là mô hình số trị vận chuyển các yếu tố chất lượng
nước phân bố đồng nhất theo độ sâu. Nó có thể tính toán sự tập trung của 6 thành
phần bảo toàn hoặc không bảo toàn được tính toán theo lưới 1 chiều hoặc 2 chiều.
ECOHAM (phiên bản 1 và 2) là mô hình số 3D kết hợp giữa module thủy
lực với module sinh thái được phát triển bởi nhóm nghiên cứu của Trường đại học

được thương mại hoá. Một đặc điểm mạnh của MIKE rất dễ sử dụng với các giao
diện Windows, kết hợp chặt chẽ với GIS (hệ thống thông tin địa lý). MIKE tích hợp
các module thuỷ lực (HD) và chất lượng nước (ECO Lab), bao gồm: thuỷ lực,
truyền tải - khuếch tán chất lượng nước. MIKE là một mô hình với nhiều tính năng
mạnh, khả năng ứng dụng rộng rãi cho nhiều dạng thuỷ vực khác nhau.
1.1.2 Tổng quan tình hình nghiên cứu trong nước
Ở nước ta, trong những năm gần đây, hướng nghiên cứu, xây dựng và sử
dụng mô hình trong nghiên cứu thủy động lực – môi trường đang rất được quan
tâm. Trong đó những nghiên cứu, điều tra, tính toán ô nhiễm môi trường các vũng
6


vịnh và khu vực ven biển - khu vực tập trung chủ yếu các hoạt động kinh tế của con
người đã, đang được tiến hành. Chương trình hợp tác với Cơ quan hợp tác Quốc tế
Nhật Bản - JICA (1995 – 1998) của Viện Tài nguyên và Môi trường biển – Viện
Khoa học và Công nghệ Việt Nam, đã bước đầu sử dụng phương pháp tính dòng vật
chất bổ sung (Flux) và quỹ nguồn (Budget) chạy trên phần mềm chuyên dụng
CABARET of LOICZ (Mỹ) để đánh giá mức độ tích tụ và khuếch tán vật chất tại
một số điểm thuộc vịnh Hạ Long. Sau đó, phương pháp nghiên cứu này còn được sử
dụng tính toán mức độ dinh dưỡng của hệ đầm phá Tam Giang – Cầu Hai (Thừa
Thiên Huế). Tuy nhiên, phương pháp này chưa tính toán đến quá trình khuếch tán
vật chất trong không gian và chỉ giới hạn tại một số điểm nhất định.
Hoàng Dương Tùng (2004), trong phạm vi luận án tiến sĩ, đã sử dụng phần
mềm DELFT 3D - WAQ đánh giá khả năng chịu tải ô nhiễm của Hồ Tây với mục
đích xây dựng căn cứ khoa học trong việc xây dựng kế hoạch bảo vệ và phát triển
Hồ Tây. Nội dung đã xem xét đến khả năng biến động các yếu tố DO, BOD, COD,
NH4+, NO3-, PO4- theo không gian 2 chiều và thời gian.
Trong khuôn khổ đề tài cấp Bộ Thủy sản, Trần Lưu Khanh và các cộng sự
cũng đã tiến hành nghiên cứu sức chịu tải và khả năng tự làm sạch tại khu vực nuôi
cá lồng bè ở Phất Cờ (Quảng Ninh) và Tùng Gấu (Hải Phòng) dựa trên quá trình

thể nghiên cứu các hệ sinh thái và nguồn lợi của vịnh Cam Ranh không nên và
không thể tách rời đầm Thủy Triều…” (GS.TS Mai Trọng Nhuận- 2008). Theo bản
đồ qui hoạch của tỉnh Khánh Hòa định hướng đến năm 2020 thì đầm Thủy Triều
ngày càng đóng vai trò quan trọng đến chất lượng môi trường nước toàn vịnh Cam
Ranh. Vì thế, tính toán lan truyền vật chất ô nhiễm vịnh Cam Ranh dựa trên công cụ
phần mềm MIKE là một hướng nghiên cứu mới mà học viên lựa chọn.

1.2 MIKE 21 HD
1.2.1 Cơ sở toán học
Mô hình MIKE 21 HD là gói công cụ trong bộ phần mềm DHI được xây
dựng bởi Viện Thủy Lực Hà Lan, đây là mô hình tính toán dòng chảy hai chiều
trong một lớp chất lỏng đồng nhất theo phương thẳng đứng.
Các phương trình nước nông
Các phương trình động lượng và liên tục tích phân trên toàn bộ cột nước h
= η+d trong các phương trình nước nông được viết lại như sau:

8


(1.1)

(1.2)

(1.3)

trong đó t là thời gian; x, y là tọa độ Đề Các; η là mực nước bề mặt; d là độ sâu của
nước tĩnh; h = η + d là độ sâu nước tổng cộng; u, v là các thành phần vận tốc theo
phương x và y; f = 2Ωsinθ là tham số Coriolis (Ω là vận tốc góc của Trái đất, θ là vĩ

độ địa lý);

môi trường xung quanh.
Biến số có đường gạch ngang biểu thị giá trị trung bình theo độ sâu. Ví dụ,
và

là các thành phần vận tốc trung bình theo độ sâu được xác định bởi:

(1.4)
Thành phần ứng suất bên Tij (i,j = x,y) bao gồm cả ma sát nhớt, ma sát rối và
chênh lệch bình lưu. Chúng được xác định bằng sử dụng công thức nhớt rối dựa
trên những biến đổi vận tốc trung bình theo độ sâu
(1.5)
Phương trình truyền tải nhiệt độ và độ muối
Các phương trình truyền tải nhiệt - muối tích phân trên toàn bộ cột nước
được viết dưới dạng:
(1.6)
(1.7)

trong đó,

và

tương ứng là nhiệt độ và độ muối trung bình theo độ sâu, F T và Fs

tương ứng là các hệ số khuếch tán ngang nhiệt độ và độ muối,

là nhóm nguồn

liên qua tới quá trình trao đổi nhiệt với khí quyển.
Phương trình truyền tải cho đại lượng vô hướng (scalar quantity)
Các phương trình truyền tải đại lượng vô hướng tích phân theo độ sâu có


sát đáy có thể được xác định từ hệ số Chezy, C, hoặc hệ số Manning, M
(1.11)
(1.12)
Ứng suất mặt
Ứng suất bề mặt

được xác định thông qua gió bề mặt. Ứng

suất mặt được tính toán dựa trên công thức thực nghiệm:
(1.13)

11


với

là mật độ không khí, cd là hệ số cản gió,

là tốc độ gió ở độ

cao 10m trên bề mặt biển. Vận tốc ma sát liên hệ với ứng suất bề mặt được cho bởi
công thức:
(1.14)

Hệ số cản cũng có thể là những giá trị không đổi hoặc phụ thuộc vào tốc
độ gió. Công thức bán thực nghiệm được đề xuất bởi Wu (1980, 1984) để xác định
giá trị của hệ số cản:
(1.15)



,

Tích phân phương trình 1.16 trên toàn bộ phần tử thứ i và sử dụng định lý
Gauss để viết lại tích phân thông lượng như dưới đây
(1.19)

13


trong đó, Ai là diện tích của phần tử thứ i, Ω là tích phân biến xác định trên A i , Гi là
biên của phần tử thứ i và ds tích phân biến dọc theo biên. n là véctơ pháp tuyến đơn
vị hướng ra ngoài biên. Các tích phân được tính bằng phương pháp cầu phương đơn
điểm, điểm cầu phương là điểm trọng tâm của phần tử, và tích phân biên được tính
dựa trên phép cầu phương tâm điểm, khi đó phương trình 1.19 được viết lại,
(1.20)

Ở đây Ui và Si tương ứng là các giá trị trung bình của U và S trên toàn bộ
phần tử thứ i và được đặt tại tâm của phần tử, NS là số cạnh của phần tử, n j véctơ
pháp tuyến ngoài đơn vị tại cạnh thứ j và Гj là chiều dài của giao diện thứ j.
Trong trường hợp 2D phép xấp xỉ Riemann được sử dụng để tính toán các
thông lượng đối lưu tại mặt phân cách của các phần tử. Sử dụng phép giải Roe để
ước lượng cho các biến phụ thuộc phía bên trái và bên phải của của giao diện. Độ
chính xác bậc hai theo không gian đạt được bằng cách sử dụng kỹ thuật tái cấu trúc
gradient tuyến tính. Các giá trị gradient trung bình được ước lượng thông qua phép
giải của Jawahar và Kamath, 2000.
Phương trình truyền tải
Các phương trình truyền tải xuất hiện trong mô hình nhiệt – muối, mô hình
rối và mô hình truyền tải. Tất cả các phương trình này đều có dạng chung. Trong
trường hợp 2D, các phương trình truyền tải có dạng tổng quát như phương trình

Biên mở
Các điều kiện biên mở có thể được đưa vào theo các dạng là lưu lượng
hoặc dao động mực nước mặt cho các phương trình thủy động lực. Với các phương
trình truyền tải, điều kiện biên có thể là các giá trị xác định hoặc giá trị gradient.
Điều kiện khô và ướt
Các giải pháp xử lý các vấn đề về biên di động (front khô và ướt) dựa trên
các nghiên cứu của Zhao và cộng sự (1994) và Sleigh và cộng sự (1998). Khi các
trường độ sâu nhỏ, vấn đề xảy ra là các phần tử được loại bỏ từ việc tính toán. Công
thức tính toán được xây dựng lại bởi sự giảm thông lượng động lượng tới giá trị
không và chỉ tính toán tới thông lượng khối lượng.

15


Độ sâu của mỗi phần tử biến đổi và các phần tử được sắp xếp thành các
loại khô, bán khô, ướt. Khi đó bề mặt các phần tử được kiểm tra để xác định các
điều kiện biên ướt.
Bề mặt của một phần tử được xác định là ngập nếu thỏa mãn hai tiêu
chuẩn: thứ nhất, độ sâu nước tại một cạnh của bề mặt phải nhỏ hơn độ sâu tới hạn
khô hdry, và độ sâu nước ở cạnh khác của bề mặt lớn hơn độ sâu độ sâu tới hạn ngập
hflood. Thứ hai, độ sâu tổng cộng của nước tĩnh tại cạnh có độ sâu nhỏ hơn h dry và
mực nước bề mặt tại cạnh khác đều phải lớn hơn giá trị 0.
Một phần tử được gọi là khô nếu độ sâu nước nhỏ hơn độ sâu giới hạn khô
hdry, và không một cạnh nào bị ngập. Phần tử này bị loại ra khỏi miền tính toán.
Một phần tử xem như là ngập một phần nếu nếu độ sâu nước lớn hơn h dry
và nhỏ hơn độ sâu giới hạn ướt, hoặc khi độ sâu nhỏ hơn h dry và một trong số các
cạnh khác là biên ngập nước. Trong trường hợp này thông lượng động lượng bằng
không và chỉ có thông lượng khối lượng được tính.
Một phần tử được gọi là ướt nếu độ sâu nước lớn hơn h wet. Trong trường
hợp này cả hai thành phần thông lượng khối lượng và thông lượng động lượng được

Thành phần bình lưu - khuếch tán được xấp xỉ bằng công thức
(1.28)
trong đó, nồng độ tức thời c* được cho bởi quá trình truyền tải biến trạng thái trong
ECO Lab khi vật chất được bảo toàn trong suốt chu kỳ

sử dụng môđun AD.

Một lợi thế chính của phương pháp này là liên kết được phương pháp giải
hiện và các vấn đề phi tuyến từ các nguồn ECO Lab phức tạp, vì vậy ECO Lab và
thành phần bình lưu - khuếch tán có thể được giải một cách riêng lẻ.
Phương pháp giải số được sử dụng trong mô hình ECO Lab là phương
pháp Euler, Runge Kutta 4, Runge Kutta 5.

17


1.3.2 Ôxy hòa tan (DO) và nhu cầu ôxy sinh hóa (BOD)
a. Ôxy hòa tan (DO)
DO là lượng ôxy hoà tan trong nước cần thiết cho sự hô hấp của các sinh
vật nước (cá, lưỡng thê, thuỷ sinh, côn trùng v.v...) thường được tạo ra do sự hoà
tan từ khí quyển hoặc do quang hợp của tảo,... Nồng độ ôxy tự do trong nước phụ
thuộc vào nhiệt độ, sự phân huỷ hoá chất, sự quang hợp của tảo và v.v... Khi nồng
độ DO thấp, các loài sinh vật nước giảm hoạt động hoặc bị chết. Do vậy, DO là một
chỉ số quan trọng để đánh giá sự ô nhiễm nước của các thuỷ vực.
Quá trình cân bằng ôxy được xem xét theo các mức độ phức tạp khác nhau
của cân bằng tùy thuộc vào mục đích của người sử dụng. Có 4 mức độ khác nhau
mô tả cân bằng khối DO, trong phạm vi nghiên cứu của luận văn, chỉ tập trung vào
mức cân bằng bậc 3. Mức cân bằng này giả thiết rằng sự biến đổi của nồng độ ôxy
là tổng hợp của các quá tương tác nước - khí quyển (mặt phân cách), quá trình đạm
hóa, nhu cầu ôxy sinh hóa, quá trình quang hợp, quá trình hô hấp, nhu cầu ôxy trầm


Quá trình hô hấp của sinh vật (g O 2/m2/ngày). Sự suy giảm nồng độ ôxy
bởi quá trình hô hấp của sinh vật tự dưỡng và dị dưỡng thông qua biểu thức phụ
thuộc nhiệt độ.
(1.29e)
Nhu cầu ôxy cho phân hủy vật chất hữu cơ tại đáy (chỉ phụ thuộc vào hàm
lượng ôxy và nhiệt độ (g/m3/ngày). Lưu ý rằng các vật chất hữu cơ trầm tích trong

19


quá trình này không tính đến thành phần trầm tích có nguồn gốc từ các nguồn ô
nhiễm. Giá trị này chỉ phụ thuộc vào nồng độ ôxy và nhiệt độ.
(1.29f)

a. Nhu cầu ôxy sinh hoá (BOD)
BOD (Biochemical oxygen Demand - nhu cầu oxy sinh hoá) là lượng oxy
cần thiết để vi sinh vật oxy hoá các chất hữu. Trong môi trường nước, khi quá trình
oxy hoá sinh học xảy ra thì các vi sinh vật sử dụng oxy hoà tan, vì vậy xác định
tổng lượng oxy hoà tan cần thiết cho quá trình phân huỷ sinh học là phép đo quan
trọng đánh giá ảnh hưởng của một dòng thải đối với nguồn nước. BOD có ý nghĩa
biểu thị lượng các chất thải hữu cơ trong nước có thể bị phân huỷ bằng các vi sinh
vật.
Dạng cân bằng của nhu cầu ôxy sinh hóa (BOD) được mô tả bằng phương
trình:
(1.30)

Giải thích các từ ngữ:
S