Đồ Án Tốt Nghiệp
Mở Đầu
Vấn đề quản lý lưu vực sông hiện nay đang là vấn đề vô cùng
cần thiết và cấp bách, trên thế giới hiện nay để có thể quản lý
được tốt người ta dùng rất nhiều các mô hình tính toán, dưới đây
là một trong số đó, QUAL2K, đây là phần mềm mô hình của
EPA’s cục bảo vệ môi trường Mỹ. chúng ta sẽ nghiên cứu mô
hình này để xem chúng hữu ích và độ chính xác như thế nào để
có thể đưa ra các cách giải quyết chính xác.
Phần 1 Cơ sở lý thuyết của mô hình QUAl2K.
Cơ sở của phương pháp là hệ phương trình thủy lực, ổn định
một chiều. Hệ phương trình thủy lực Saint – Venant một chiều từ
hệ phương trình tính được cân bằng thủy lực.
Hệ phương trình phát tán
( . ) ( . )
j
j j j j
C
div u C div D gradC R
t
σ
σ
→
+ = +
Từ đó ta tìm được phương trình cơ bản của mô hình:
( . ) ( . )
j j
x j L j
C C
u C D R
t x x x

giới hạn đó là:
• Đoạn sông : max 25
• Phần tử tính toán :
≤
20% số đoạn sông hoặc tổng bằng 250
• Phần tử nguồn : max 7
• Phần tử nối : max 6
• Phần tử vào và ra : max25.
Mục đích đầu tiên của mô hình chất lượng nước là tạo ra công
cụ có khả năng mô phỏng tính chất thủy lực và chất lượng nước.
Mô tả công thức
Công thức cơ bản được giải bởi QUAL2E là một chiều, phát
tán dọc trục, công thức vận chuyển khối lượng bao gồm bình lưu,
phát tán, pha loãng, thành phần phản ứng, và sự tác động qua lại
giữa chúng, nguồn sông và lắng đọng. Công thức có thể được
viết như sau:
( )
( )
( )
x L
x
x
C
A D
A cu
M dc
x
dx dx A dx S
t x x dt
σ

S: nguồn sông hoặc lắng đọng (MT
-1
).
Bởi vì M = V.C, chúng ta có thể viết:
( )
.
M VC C V
V C
t t t t
σ σ σ σ
σ σ σ σ
= = +
Trong đó V = A
x
.dx thể tích phát tán (L
3
). Nếu chúng ta coi
dòng chảy là ổn định,
0
Q
t
σ
σ
=
, thì
0
V
t
σ
σ

những đặc điểm sau:
• Một chiều. Lòng sông là những nguồn nước trộn lẫn theo
chiều dọc và chiều sâu.
• Nhánh sông. Hệ thống có thể bao gồm một sông chính với
các sông nhánh.
• Khối nhiệt ngày đêm. Khối nhiệt và nhiệt độ được mô phỏng
như một công thức khí tượng học trên một mức độ thời gian.
• Tính chất thủy lực là ổn định. Đồng nhất, dòng chảy ổn định
được mô phỏng
• Động học chất lượng nước cả ngày đêm. Chất lượng nước
thay đổi mô phỏng theo các mức độ thời gian.
• Nhiệt và khối lượng đầu vào. Điểm và không điểm chịu tải
và nước chảy ra đều được mô phỏng.
QUAL2K còn bao gồm các phần tử mới:
• Phần mềm môi trường và giao diện. Q2K là một công cụ
trong môi trường Microsoft Windows. Số lượng tính toán dùng
chương trình Fortran 90. Excel được sử dụng để hiển thị đồ thị
trên giao diện cho người sử dụng. Tất cả các giao diện này có tác
dụng là chương trình trong Microsoft Office dùng ngôn ngữ:
Công cụ Visual Basic( VBA).
• Mô hình chia nhỏ. Q2E chia hệ thống thành các đoạn
sông gồm các phần tử có khoảng cách bằng nhau.
Q2K phân chia hệ thống thành các đoạn sông và các phần tử.
Thêm vào đó, khối lượng và các dòng chảy ra có thể vào nhiều
phần tử.
• Sự hình thành cacbon BOD. Q2K sử dụng 2 dạng
Cacbon BOD tượng trưng carbon hữu cơ. Hai dạng đó là dạng
Bùi Ngọc Hiếu_QLMT_K48
3
Đồ Án Tốt Nghiệp

vào và đầu ra của mô hình sẽ được thực hiện bằng công cụ trong
Excel, tất cả các công thức trong Excel dùng ngôn ngữ: Visual
Basic for Applications (VBA). Tất cả các công thức tính toán bằng
công cụ Fortran 90 được thi hành mau lệ. Tiếp sau đây là các
bước có bao nhiêu mô hình có thể cài đặt lên máy tính của bạn và
sử dụng chúng để làm mô phỏng .
Bước 1: copy the file, Q2Kv2_07.zip đến đường dẫn (ví dụ, C:\)
khi file được giải nén nó sẽ cho các file sau : file Excel
(Q2KMasterv2_07.xls), và một file chạy (Q2KFortran2_07.exe). Đầu
Bùi Ngọc Hiếu_QLMT_K48
4
Đồ Án Tốt Nghiệp
tiên trên giao diện của Q2K cho phép bạn chạy Q2K và biểu lộ
kết quả của nó. Thứ hai là Fortran có thể thực hiện được công
việc thực tế tính toán mô hình. Sẽ có hai file trong đường dẫn
giống nhau để mô hình có thể chạy chính xác. Chú ý sau khi bạn
chạy mô hình, một số file sẽ tự động được tạo ra bởi Fortran có
thể trao đổi thông tin với Excel.
Chú ý không xóa file .Zip. Nếu một vài lý do, bạn sửa Q2k,
bạn có thể sử dụng file zip để cài đặt lại mô hình.
Bước 2: tạo ra file theo đường dẫn C:\Q2Kv2_07 gọi là file dữ
liệu Datafiles.
Bước 3: mở Excel và chắc chắn macro security ở mức trung
bình (tranh 1) có thể yêu cầu sử dụng : Tools → Macro → Security.
Chắc chắn mức medium sẽ được chọn
Figure 1 The Excel Macro Security Level dialogue box. In order to run
Q2K, the Medium level of security should be selected.
Mở Q2KMasterFortranv2_07.xls. Khi bạn làm việc với nó hộp
thoại Macro Security sẽ hiện ra như sau:
Bùi Ngọc Hiếu_QLMT_K48

Figure 5 This window is displayed showing the progress of the model
computations as executed in Fortran. It allows you to follow the progress of a
model run.
Chương trình sẽ mô phỏng sông chính với hai nhánh sông.
Nếu chương trình làm việc chính xác hộp thoại sau đây sẽ xuất
hiện nếu bạn chạy thành công.
Bùi Ngọc Hiếu_QLMT_K48
7
Đồ Án Tốt Nghiệp
Ấn Ok, tiếp theo hộp thoại sau sẽ xuất hiện.
Hộp thoại trên sẽ cho phép bạn chọn phần của hệ thống bạn
muốn vẽ đồ thị. Như đã thấy, nó mặc định là sông chính. Ấn Ok
và nhìn thấy thời gian chạy của sông chính. Chú ý tất cả các đồ
thị đều được cập nhật khi nhấn OK.
Ngắt một lúc bạn nhìn thấy đồ thị của một nhánh sông, bạn
nhấn nút dưới bên trái bị che khuất.
Bùi Ngọc Hiếu_QLMT_K48
8
Đồ Án Tốt Nghiệp
Nguyên nhân là do đồ thị hộp thoại được chọn xuất hiện. Kéo
xuống bạn có thể chọn một nhánh khác.
Bước 7: Trên QUAL2K Worksheet click nút Open Old File. Mở
đường dẫn C:\Q2Kv2_07\DataFiles. Bạn nhìn thấy một file mới
được tạo ra với tên chỉ rõ ở cột 9 (trong trường hợp trên bức tranh
thứ 3 là Bogus062807.q2k). click nút hủy bỏ cacel quay trở lại
Q2K.
Chú ý trong thời gian Q2K chạy. Một file dữ liệu sẽ được tạo
ra với tên file chỉ rõ trong cột 9 trên QUAL2K Worksheet (Figure 3).
Chương trình tự động thêm vào phần mở rộng .q2k cho tên file.
Từ đó nó sẽ đè lên phiên bản của file trước, chắc chắn tạo ra sự

Point withdrawal
Point withdrawal
Headwater boundary
Downstream boundary
Point source
Figure 6 QUAL2K segmentation scheme for a river with no tributaries.
Hệ thống gồm các sông nhánh (hình7). Số lượng các đoạn
sông được đánh số bắt đầu từ đoạn 1 và tăng dần ở thượng nguồn
của con sông chính. Khi đến chỗ nối với một nhánh sông là một
đoạn sông số thứ tự tiếp tục được đánh từ thượng nguồn từ nhánh
sông này. Quan sát cả thượng nguồn và các nhánh sông các số là
liên tiếp nhau theo một dãy sắp xếp tương tự đến các đoạn sông.
Chú ý các nhánh sông lớn của hệ thống đều được quy về như một
đoạn sông. Đặc biệt thực tế này rất quan trọng bởi vì phần mềm
cung cấp đồ thị của đầu ra mô hình trên một đoạn sông cơ bản.
Phần mềm tạo ra các đồ thị riêng biệt trên hệ thống sông chính
cũng như các sông nhánh.
Bùi Ngọc Hiếu_QLMT_K48
10
Đồ Án Tốt Nghiệp
19
18
17
16
19
18
17
16
1
5

6
8
7
6
9
11
10
9
11
10
24
23
22
25
HW#1
HW#2
HW#3
HW#4
(a) A river with tributaries (b) Q2K reach representation
Main stem
T
r
i
b
1
T
r
i
b
2

3.1.Cân bằng dòng chảy.
Như đã được miêu tả ở phần trước, đơn vị cơ bản của mô
hình Q2K là phần tử. Một dòng chảy ổn định cân bằng là phương
tiện cho mỗi mô hình phần tử.
ioutiinii
QQQQ
,,1
−+=
−
[1]
Trong đó Q
i
là lượng chảy ra từ phần tử i vào phần tử xuôi
dòng i + 1 [m
3
/d], Q
i–1
là lượng chảy vào từ phần tử ngược dòng i
– 1 [m
3
/d], Q
in,i
là tổng lượng chảy vào trong phần tử từ điểm
nguồn và không phải điểm nguồn [m
3
/d], và Q
out,i
là tổng lượng
chảy ra từ phần tử đó đến điểm chảy ra và không phải điểm chảy
ra [m

1
,,
1
,,,
[2]
Trong đó Q
ps,i,j
là lượng chảy vào từ điểm nguồn thứ j đến
phần tử i, psi tổng số điểm nguồn đến phần tử i, Q
nps,i,j
là lượng
chảy vào từ điểm không phải điểm nguồn chảy tới phần tử i, và
npsi là tổng số điểm không phải điểm nguồn chảy vào phần tử i.
Tổng lượng chảy ra từ các nguồn chảy ra được tính toán như
sau:
∑∑
==
+=
npai
j
jinpa
pai
j
jipaio
QQQ
1
,,
1
,,ut,
[3]

nước sẽ được chọn làm phương tiện tính toán.
• Nếu chiều rộng và chiều cao của đập bằng 0 và hệ số
đường cong ( a và
α
) được nhập vào. Phương tiện rating curves
được chọn làm phương tiện tính toán.
• Nếu không có quy định trước là mét. Q2K sử dụng công
thức Manning.
3.2.1 Đập nước
Bức tranh 11 cho thấy có bao nhiêu đập nước được miêu tả
trong Q2K. Chú ý một cái đập nước chỉ có thể xảy ra ở điểm cuối
của một phần tử đơn của một đoạn sông, bức tranh 11 cho thấy
các thông số sau H
i
là chiều sâu của phần tử ngược dòng của đập
nước [m], H
i+1
là chiều sâu của phần tử xuôi dòng của đập [m],
elev2
i
độ cao so với mực nước biển điểm cuối của phần tử ngược
dòng [m], elev1
i+1
độ cao so với mực nước biển điểm đầu của
phần tử xuôi dòng. H
w
độ cao của đập trên elev2
i
, H
d

i+1
Figure 11 A sharp-crested weir occurring at the boundary between two
reaches.
H
h
là độ cao ở đỉnh bên trên đập [m], B
w
là chiều rộng của đập
[m]. Chú ý là chiều rộng của đập khác với chiều rộng của phần
tử, B
i
Đây là một dạng đập trong đó H
h
/H
w
< 0.4, dòng chảy có liên
quan đến đầu nguồn (Finnemore and Franzini 2002)
2/3
83.1
hwi
HBQ =
[4]
Trong đó Q
i
là lượng chảy ra từ phần tử ngược dòng của đập,
m
3
/s, B
w
, H

iiiid
HelevHelevH
[7]
Chú ý độ hạ thấp có thể sử dụng để tính toán lượng Oxy và
CO
2
di chuyển qua đập ( xem trang 55 và 60).
Tại các khu vực mặt cắt ngang, chiều sâu, diện tích bề mặt và
thể tích phần tử i có thể được tính toán như sau.
iiic
HBA =
,
[8]
ic
i
i
A
Q
U
,
=
[9]
iiis
xBA ∆=
,
Bùi Ngọc Hiếu_QLMT_K48
15
Đồ Án Tốt Nghiệp
iiii
xHBV ∆=

c
=
[12]
H
A
B
c
=
[13]
Diện tích bề mặt và thể tích có thể tính như sau
xBA
s
∆=
xBHV ∆=
Số mũ b và
β
được đặc trưng trong bảng 1 chú ý tổng của b
và
β
phải kém hơn hoặc bằng 1. Nếu đây không phải là trường
hợp mà chiều rộng sẽ giảm với sự gia tăng dòng chảy. Nếu tổng
của chúng bằng 1 kênh sông là hình chữ nhật.
Table 1 Typical values for the exponents of rating curves used to
determine velocity and depth from flow (Barnwell et al. 1989).
Equation Exponent
Typical
value
Range
b
aQU =

c
=
[14]
Trong đó Q là lưu lượng dòng chảy [m
3
/s], S
0
độ dốc đáy sông
[m/m] , n là hệ số gồ ghề, A
c
diện tích mặt cắt ngang [m
2
] và P là
chu vi thấm ướt [m].
Q, U
B
0
1
1
s
s1
s
s2
H
S
0
B
1
Figure 12 Trapezoidal channel.
Diện tích mặt cắt ngang của một lòng sông hình thang được

10/3
5/2
2
21
2
110
5/3
)(5.0
11)(
−
−−
++






++++
=
kss
sksk
k
HssBS
sHsHBQn
H
[17]
Trong đó k = 1, 2, …n. n là số lần lặp. Ban đầu ước chừng
H
0

A
B
c
=
[20]
Chiều rộng bên trên B
1
[m] thể được tính toán như sau.
HssBB
ss
)(
2101
++=
Diện tích bề mặt và thể tích của phần tử có thể được tính toán
như sau:
xBA
s
∆=
1
xBHV ∆=
Đề xuất giá trị hệ số manning cho trong bảng 2, n đặc trưng
cho giá trị dòng chảy và chiều sâu (Gordon et al. 1992). Chiều sâu
giảm trong chiều dòng chảy thấp, liên quan đến sự dao động
thường xuyên được tăng lên. Giá trị của hệ số manning đã được
công bố từ 0.015 của lòng sông nhẵn nhịu đến 0.15 các lòng sông
gồ ghề nó miêu tả tình trạng dòng chảy có khả năng tạo thành bãi
ngầm (Rosgen, 1996). Điều kiện tới hạn của độ sâu ước lượng chất
lượng nước đại thể là kém hơn bãi ngầm sâu và nó liên quan đến
tính chất gồ ghề của độ cao.
Table 2 The Manning roughness coefficient for various open channel

i
elev1
i+1
Figure 13 A waterfall occurring at the boundary between two reaches.
Qual2k sẽ tính toán dòng chảy trong trường hợp độ cao so với
mực nước biển rất dốc trong ranh giới giữa hai đoạn sông , công
thức 7 dùng để tính toán sự hạ thấp dòng chảy. Chú ý sự hạ thấp
này chỉ tính toán khi độ cao so với mực nước biển xuôi dòng kết
thúc ở đoạn sông là lớn hơn điểm bắt đầu của đoạn sông xuôi
dòng tiếp theo nghĩa là elev2
i
> elev1
i+1
.
3.3Travel Time (Thời gian di chuyển)
Thời gian lưu của mỗi phần tử được tính toán như sau:
k
k
k
Q
V
=
τ
(1)
Trong đó
τ
k
là thời gian lưu của phần tử thứ k [d]. V
k
là thể

t,j
là thời gian di động.
3.4 Phát tán dọc trục
Bùi Ngọc Hiếu_QLMT_K48
19
Đồ Án Tốt Nghiệp
Hai lựa chọn được sử dụng để xác định sự phát tán dọc trục
cho ranh giới giữa hai phần tử. Đầu tiên, người sử dụng có thể
đánh giá, giá trị nhập vào trên Reach Worksheet. Nếu người sử
dụng không nhập giá trị, một công thức bên trong sẽ được dùng
tính toán phát tán cơ bản tính chất thủy lực lòng sông (Fischer et al.
1979),
*
22
,
011.0
ii
ii
ip
UH
BU
E =
[23].
Trong đó E
p,i
là hệ số phát tán dọc trục giữa phần tử i và phần
tử i+1 [m
2
/s], U
i

Độ phát tán của mô hình là E
i
( giá trị này sẽ được sử dụng
tính toán trong mô hình).
Nếu E
n,i
≤ E
p,i
, độ phát tán của mô hình, E
i
= E
p,i
− E
n,i
.
Nếu E
n,i
> E
p,i
, độ phát tán mô hình bằng Ei

= 0.
Trong một số trường hợp dưới đây, kết quả độ phát tán của
mô hình là lớn hơn độ phát tán vật lý. Như vậy, thì sự trộn lẫn
phát tán sẽ cao hơn thực tế, Chú ý dòng sông chảy ổn định, tập
trung độ dốc là không đáng kể. Nếu sự khác nhau là quan trọng,
chỉ lựa chọn các phần tử có chiều dài nhỏ hơn số phát tán mô
hình trở thành nhỏ hơn số phát tán vật lý.
3. Nhiệt độ của mô hình.
Như hình 14, cân bằng nhiệt cần tính toán sự trao đổi nhiệt từ

−+−+−−=
+−
−
−
−
cm 100
m
cm 100
m

cm 10
m

,,
36
3
,
1
'
1
'
1
,
1
1
ipww
is
ipww
ia
ipww

Q
T
V
Q
T
V
Q
dt
dT
ρρρ
[26].
Trong đó T
i
nhiệt độ phần tử i, [
0
C], t thời gian [d], E’
I
hệ số
phát tán lớn giữa phần tử i và phần tử i+1 [m
3
/d], W
h,i
mạng nhiệt
từ các điểm nguồn và không phải điểm nguồn trong phần tử
i[cal/d]
ρ
w
tỷ trọng của nước [g/m
3
], C

,
'
+
∆+∆
=
ii
ici
i
xx
AE
E
[27].
Chú ý hai loại điều kiện biên được sử dụng đến điểm cuối
cùng của dòng chảy xuôi dòng của sông, (1) điều kiện phát tán
bằn 0 và (2) điều kiện biên bắt buộc ở điểm cuối dòng chảy, cơ
hội lựa chọn tạo ra trên Downstream Worksheet.
Mạng nhiệt chiếu xuống từ nguồn được tính toán như sau
(recall Eq. 2)








+=
∑∑
==
npsi

JJJJIJ −−−+= )0(
[29].
I(0) là bức xạ sóng ngắn của mặt trời tại bề mặt nước, J
an
là
bức xạ sóng dài trong không khí, J
br
phản xạ sóng dài từ nước, J
c

là độ dẫn điện, và J
e
là sự bốc hơi.
Tất cả các dòng chảy đều biểu diễn bằng cal/cm
2
/d.
air-water
interface
solar
shortwave
radiation
atmospheric
longwave
radiation
water
longwave
radiation
conduction
and
convection

bức xạ ngoài khí quyển ( tầng trên của khí quyển trái đất),
[cal/cm
2
/d], a
t
không khí loãng, a
c
mây mỏng, R
s
suất phản chiếu
(phản xạ nhỏ), S
f
hệ số hiệu quả bóng tối bởi sinh vật và địa
hình).
Sự phát xạ ngoài khí quyển được đánh giá như sau:
α
sin
2
0
0
r
W
I =
[31].
W
0
là hằng số mặt trời [1367 W/m
2
or 2823 cal/cm
2






−=
imetrueSolarT
[33].
Trong đó:
timezoneLeqtimelocalTimeimetrueSolarT
lm
×−×−+= 604
[34].
Trong đó trueSolarTime là thời gian mặt trời xác định từ vị trí
thực tế của mặt trời trong bầu trời [minutes], localTime là thời gian
địa phương [thời gian chuẩn của địa phương], L
lm
kinh độ của địa
phương, timezone đới thời gian của địa phương liên quan đến giờ
chuẩn căn cứ theo kinh tuyến (GMT). Ví dụ như -8h ở đới Thái
Bình Dương là giờ chuẩn, thời gian địa phương ở các đới được
chọn trên QUAL2K Worksheet. Giá trị eqtime tượng trưng cho sự
khác nhau giữa thời gian mặt trời chính xác và thời gian mặt trời
trung bình.
QUAL2K tính toán độ nghiêng của mặt trời, múi giờ, độ cao
mặt trời và bán kính tiêu chuẩn ( khoảng cách giữa trái đất và mặt
trời), thời gian lúc mặt trời mọc và lúc mặt trời lặn sử dụng bởi
thuật toán Meeus (1999) như là một công cụ bởi nhánh nghiên
cứu bức xạ bề mặt NOAA’s.
NOAA sẽ xác định vị trí mặt trời dựa vào QUAL2K bao gồm

fac
ea
1
−
=
[36].
Trong đó n
fac
là hệ số mật độ không khí biến đổi từ xấp xỉ 2 ở
chỗ bầu trời sạch đến 4 hoặc 5 khu vực thành thị nhiều sương
mù. Hệ số phân tán phân tử a
1
tính toán như sau:
ma
101
log054.0128.0 −=
[37].
Trong đó m là khối lượng khí nhìn thấy được, tính toán như
sau:
253.1
)885.3(15.0sin
1
−
++
=
d
m
αα
[38].
α

biểu xấp xỉ 0.8), và elev là độ cao mặt đất tính bằng mét.
Phép đo nhằm xác định bức xạ mặt trời có thể áp dụng ở một
vài nơi. Ví dụ NOAA’s nghiên cứu suất phản chiếu bề mặt (ISIS)
có dữ liệu thay đổi từ Mỹ. ( Chọn
cả hai mô hình bức xạ mặt trời Bras or Ryan-Stolzenbach với hệ số
mật độ không khí thích hợp hoặc hệ số truyền không khí một
công cụ đặc biệt để lý tưởng hóa so sánh với bức xạ mặt trời
trước đó với giá trị được cân nhắc ở từng địa phương.
Sự suy giảm mây: sự giảm bức xạ mặt trời do bao phủ của
mây được tính toán với
2
65.01
Lc
Ca −=
[40].
C
L
là hệ số bao phủ bầu trời bởi mây.
Reflectivity. Reflectivity được tính toán như sau:
B
ds
AR
α
=
[41].
Trong đó A và B là hai hệ số liên quan đến bao phủ mây
( bảng 3).
Table 3 Coefficients used to calculate reflectivity based on cloud cover.
Cloudiness Clear Scattered Broken Overcast
C

d K
4
),
T
air
nhiệt độ khí quyển [
0
C], ε
sky
hệ số phát xạ không khí [không
thứ nguyên], R
L
hệ số phản xạ sóng dài [không thứ nguyên]. Độ
phát xạ là tỷ số giữa bức xạ sóng dài từ một vật đối với bức xạ
phát ra từ một vật hoàn toàn trong khoảng nhiệt độ như nhau.
Suất phản chiếu nói chung là nhỏ và được cho là bằng 0.03.
Mô hình bức xạ sóng dài không khí sẽ được chọn trên Light
and Heat Worksheet Qual2k. Ba phương pháp lựa chọn có thể sử
dụng trong qual2k tượng trưng cho hệ số phát xạ (ε
sky
):
1. Brunt ( mặc định ).
Công thức Brunt’s là một mô hình kinh nghiệm thường được
sử dụng trong mô hình chất lượng nước (Thomann and Mueller
1987),
airbaclear
eAA +=
ε
Trong đó A
a







=
a
air
clear
T
e
ε
Trong đó e
air
là áp suất khí quyển [mmHg], T
a
nhiệt độ khí
quyển
0
K, hệ số 1.333224 chuyển áp suất hơi từ mmHg sang
Bùi Ngọc Hiếu_QLMT_K48
25