Bài giảng Mô hình nước ngầm - Phần 1: Cơ sở lý thuyết - Nguyễn Mai Đăng

10/6/2013
1
MÔ HÌNH NƯỚC NGẦM
Phần 1: Cơ sở lý thuyết
Nguyễn Mai Đăng
Bộ môn Thủy văn & Tài nguyên nước
dang@wru.vn
0989.551.699
Quá trình ứng dụng mô hình
1. Xác định vấn đề cần mô phỏng
– Các yếu tố cần mô hình hóa
– Mối quan hệ và tương tác giữa chúng
Xác định và mô tả vấn đề
Cơ sở lý thuyết mô hìnhSố liệu
– Độ chính xác
2. Xác định cơ sở lý thuyết và phát triển mô hình
– Mô tả vềmặt toán học
– Loại mô hình
– Phương pháp số – lập trình trên máy tính
– Lưới tính toán, điều kiện biên, đ/k ban đầu
3. Hiệu chỉnh thông sốmô hình
– Tạm thời xác định các thông số của mô hình
– So sánh kết quả tính toán từmô hình với số liệu thực đo
– Hiệu chỉnh các thông số cho đến khi kết quảmô hình đạt
được độ chính xác đặt ra
Kiểm định mô hình &
Phân tích độ nhạy
Hiệu chỉnh mô hình &
Đánh giá thông số
Phát triển mô hình
.
4. Kiểm định
– Thiết lập bộ số liệu độc lập (so với giai đoạn hiệu chỉnh) 
– Chạy mô hình (giữ nguyên bộ thông số xác định được từ giai
đoạn hiệu chỉnh)
– So sánh, đánh giá sự sai khác giữa kết quả tính toán và thực
đo
Tổng hợp và hoàn
chỉnh mô hình
Ứng dụng mô hình
Trình diễn các kết quả
10/6/2013
2
Các công cụ để giải quyết các vấn đề của 
nước ngầm
• Các phương pháp tương tự và vật lý
– Some of the first methods used. 
• Các phương pháp phân tích
ì ú ả ậ ế– Những g  ch ng ta đã th o lu n cho đ n 
thời điểm này?
– Gặp khó khăn khi: các biên không đều, 
điều kiện biên khác nhau, tính chất không
đồng nhất, không đẳng hướng, nhiều giai
đoạn, phi tuyến
• Các phương pháp số
– Chuyển các phương trình vi phân từng
www.epa.state.oh.us
phần (đạo hàm riêng) của nước ngầm
sang hệ phương trình vi phân thường
(đạo hàm toàn phần) hoặc chuyển sang 
các phương trình đại số để giải ra nghiệm
bài toán.
www.isws.illinois.edu
Mô hình khái niệm
• Mô tả đại diện của 
hệ thống nước 
ngầm kết hợp với 
giải thích các điều       
kiện địa chất và 
thủy văn.
• Quá trình nào là 
quan trọng đối với 
mô hình?
• Các biên là gì?
• Các giá trị thông số         
có tồn tại không?
• Các giá trị thông số 
phải được thu thập 
là gì?
10/6/2013
3
Chúng ta thực sự muốn tính toán cái gì?
• Dòng chảy hướng ngang trong tầng ngậm 
nước có áp và thấm yếu
• Các phương trình chủ đạo
t
hSxxQhh
b
K
y
hT
yx
hT
x
W
w
wwyx ∂
∂=−±−+⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂
∂
∂+⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂ ∑
=1
0 )()('
' δ
Bề mặt đất
Cột nước tầng có áp
Dòng chảy Thấm Nguồn chảy
vào/chảy ra
Lượng trữ
• Các điều kiện biên
• Các điều kiện ban đầu
Tầng đá gốc
Tầng có ápQx
K
x
yz
h
Lớp có áp (không thấm)
b
Phương pháp sai phân hữu hạn
(Finite Difference Method)
• Phương pháp sai phân hữu hạn
ế– Thay th  các đạo hàm trong phương trình nước 
ngầm bằng các ‘xấp xỉ chuỗi Taylor’
– Thiết lập hệ phương trình đại số để giải tìm 
nghiệm
10/6/2013
4
Chuỗi Taylor
• Khai triển chuỗi Taylor đối với h(x) tại điểm 
( ∆ ) lâ ậ điểx+ x n c n  m x
• Nếu cắt chuỗi số từ các trị số n+1 trở đi, thì sai 
số của chuỗi số là
LL +∆′′∆′∆=∆ )(
!
)()(
!2
)()()()(
2
xh
n
x++xhx+xhx+x h xx+h n
n
        :
)()(
)!1(
)( 1
max
1
1
+∆=∆< nn+
n+
xxh
n+
xerror ϑ
Đạo hàm bậc 1 “tiến” (+∆x)
• Phát triển chuỗi Taylor “tiến” cho hàm mực nước h(x) tại vị 
trí lân cận điểm x là:
hxhhh )()()()()(
2
′′∆′
Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau:  
L+x+xx+x xx+
!2
⋅∆=∆
x
x-hxx+h 
x
hxh ∆
∆≈∂
∂=′ )()()(
x
xx ∆+
x∆
xxx ∆−
x∆
10/6/2013
5
Đạo hàm bậc 1 “lùi” (‐∆x)
• Phát triển chuỗi Taylor “lùi” cho hàm mực nước h(x) tại vị trí 
lân cận điểm x là:
)( 2∆
Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau:  
L-xhx+xhx-xhxx-h )(
!2
)()()( ′′′⋅∆=∆
xx--hxhhxh ∆≈∂=′ )()()(
xx ∆∂
x
xx ∆+
x∆
xxx ∆−
x∆
Đạo hàm bậc 2
+xhx+xhx+xhxx+h )()()()()(
2
′′∆′∆∆ L 
!2
⋅=
Trong đó h”(x) là đạo hàm bậc 2 được tính xấp xỉ là: 
L−′′∆′⋅∆−=∆ )(
!2
)()()()(
2
xhx+xhxx h xx+h
22
2 )()(2)()(
x
xx-+hxh-xx+h
x
hxh ∆
∆∆≈∂
∂=′′
10/6/2013
6
Xấp xỉ sai phân hữu hạn
i
hh
hxh =
)(
)(
-hhh iii =′ +1 2 11 2 +hh-h=h iiii′′ −+-hhh iii =′ −1
i
i
i
h xh 
hxxh
xx
′=′
=∆−
=∆+
−
+
)(
)( 1
1
x
x∆ x∆
x∆ x ∆x ∆
Đạo hàm 
lùi bậc 1
Đạo hàm 
tiến bậc 1
Đạo hàm 
trung tâm 
bậc 2
Lưới và rời rạc
(Grids and Discretization)
• Quá trình rời rạc
• Lưới được thiết lập để phủ kín 
y, j
Lưới
vùng tính toán
• Mục tiêu – dự báo mực nước tại 
các điểm nút của lưới
– Xác định ảnh hưởng của bơm
– Dòng chảy từ sông, 
• Phương pháp sai phân hữu hạn
– Được sử dụng phổ biến do tính
i,j
i,j+1
i+1,ji-1,j
i,j-1
Vùng tính toán
Điểm nút
∆y
đơn giản của phương pháp, dễ 
hiểu, dễ sử dụng và lại đạt kết 
quả  tốt
– Sử dụng tốt cho vùng địa hình 
đơn giản
x, i∆x
Ô lưới
10/6/2013
7
Lưới không gian 3 chiều
(Three‐Dimensional Grids)
• Một hệ thống tầng ngậm nước được chia thành các khối chữ 
nhật bởi các lưới. 
L ới đ ổ hứ bởi á hà (i) ộ (j) lớ (k) ỗi ộ• ư   ược t  c c   c c  ng  , c t  ,  p  , m  m t 
khối (block) được gọi là một “ô lưới“ (cell)
• Các loại lớp địa tầng
– Có áp
– Không áp
– Chuyển đổi
Có thể mô phỏng 
cho các lớp có vật 
liệu khác nhau
Dòng chảy tầng ngậm nước có áp 1 chiều
(1‐D Confined Aquifer Flow)
• Dòng chảy đồng nhất, 
đẳng hướng,  1 chiều, và 
có áp.
Bề mặt đất (Ground surface)
• Phương trình mô tả cho 
dòng chảy xác định ở trên:
• Điều kiện ban đầu (initial 
diti )
t
hS
x
hT
x ∂
∂=⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂
Tầng ngậm nước (Aquifer)
x
yz
hB
Lớp không thấm
b
hA
∆x
i = 0 1 2 3 4 5 6 7 8 9 10
Nút (Node)
con on :
• Điều kiện biên (boundary 
conditions):
Ô lưới (Grid Cell)
mtLh
mth
 5.1),(
 1.6),0(
=
=
∆x= 1 m, L = 10 m, b = 1.5 m
hA = 6.1 m, hB = 1.5 m, 
K = 0.5 m/d, S = 0.02
mxh 1.6)0,( =
10/6/2013
8
Tính toán xấp xỉ cho các đạo hàm
(Derivative Approximations)
• Phương trình nước ngầm 1 chiều, đồng 
nhất, đẳng hướng, có áp:
• Đạo hàm bậc 2 đối với (WRT with respect 
to) x
t
h
T
S
x
h
∂
∂=∂
∂
2
2
2
11
2
2 2
x
+hh-h 
x
h li
l
i
l
i
l
i ∆
≈∂
∂ −+
li ,1+
ix,
lt,
1, +li
li ,1−
1, −li
x∆
t∆
li,
• Đạo hàm bậc 1 đối với (WRT) t
t
-hh
t
h li
l
i
l
i ∆
≈∂
∂ +1
l
-hh
t
h li
l
i
l
i ∆
≈∂
∂ −1
Tiến lùi
Phương pháp sai phân sơ đồ hiện
(Explicit Method)
• Sử dụng thông tin (H, Q) tại các nút ở 
bước thời gian trước để tính cho nút ở  lt,
bước thời gian hiện tại (sơ đồ hiện).
• Tính toán cho từng nút một (tất cả các 
điểm không gian và thời gian) cho đến 
hết miền tính toán.
• Phải cẩn thận khi chọn bước thời gian 
(∆t), nếu lớn quá sẽ không ổn định trong 
tính toán, nếu nhỏ quá thì tính toán lâu, 
ủ ỏ ể ổ
li ,1+
ix,
1, +li
li ,1−
1, −li
x∆
t∆
li,
do vậy chọn vừa đ  nh  đ  bài toán  n 
định và cũng đáp ứng được sai số cho 
phép.
t
h
T
S
x
h
∂
∂=∂
∂
2
2
t
hh
T
S
x
hhh li
l
i
l
i
l
i
l
i-
∆
−=∆
+− ++ 1
2
11 2
Từ pt vi phân
chuyển sang sai phân (xấp xỉ)
10/6/2013
9
Phương pháp sai phân sơ đồ hiện 
(Explicit Method)
t
hh
T
S
x
hhh li
l
i
l
i
l
i
l
i-
∆
−=∆
+− ++ 1
2
11 2
( )lllll hhhhh 1 2+
2x
t
S
Tr ∆
∆=
l
i
l
i
l
i
l
i
l
i- hhhhhx
t
S
T −=+−∆
∆ +
+
1
112 )2(
Đặt:
iii-ii r 11 ++−+=
Bước thời gian l+1
Chưa biết
Bước thời gian l
Đã biết
Phương pháp sai phân sơ đồ hiện 
(Explicit Method)
Bề mặt đất (Ground surface)
Tầng ngậm nước
hB
Lớp không thấm 
b
hA
∆x
i = 0 1 2 3 4 5 6 7 8 9 10
Nút tính toán
Ô lưới
)2( 11
1 l
i
l
i
l
i-
l
i
l
i hhhrhh +
+ +−+=
T
Sxrt
x
t
S
Tr
2
2
∆=∆
∆
∆=
∆x= 1m, L = 10m, b = 1.5m
hA = 6.1m, hB = 1.5m, 
K = 0.5m/d, S = 0.02
Xem xét: 
r = 0.48
r = 0.52
min5.180128.02 ≈=∆=∆ d
T
Sxrt
min200139.0 ≈=∆ dt
10/6/2013
10
Kết quả tính theo sơ đồ hiện
(∆t = 18.5 min; r = 0.48 < 0.5)
Kết quả đường mực nước 
ổn định (ko bị dao động)
Kết quả tính theo sơ đồ hiện
(∆t = 20 min; r = 0.52 > 0.5)
Kết quả đường mực nước 
ko ổn định, do ∆t lớn
10/6/2013
11
Phân tích quá trình không ổn định đó 
diễn ra như thế nào?
• Tại thời gian t = 0 Æ ko có dòng 
chảy
• Khi t > 0 Æ có dòng chảy
Bề mặt đất
• Lượng nước có thể thoát ra từ 
lượng trữ trong một ô lưới trong 
khoảng thời gian ∆t:
• Nước chảy ra khỏi một ô lưới 
trong thời khoảng ∆t là:
Tầng ngậm nước
hB
Lớp ko thấm
b
hA
∆x
i = 0 1 2  i-1 i i+1  8 9 10
∆x
hxSV ∆⋅∆⋅=∆ )1(1
t
x
hTV ∆∆
∆=∆
2/2
• Theo nguyên lý bảo toàn khối 
lượng:     ∆V2 ≤ ∆V1
Ô lưới i
hxSt
x
hT ∆⋅∆⋅≤∆⋅∆
∆⋅ 
2/
Vậy nếu r > 0.5 Æ khoảng thời 
gian qua lớn Æ ô lưới nước ngầm 
không thể chứa đủ nước Æ không 
ổn định
2
1
2 ≤∆
∆=
x
t
S
Tr
Phương pháp sai phân sơ đồ ẩn
(Implicit Method)
• Sử dụng thông tin (H, Q) ở một 
nút thời gian trước tính toán lt,           
cho 3 nút thời gian sau (sơ đồ 
ẩn).
• Giải đồng thời cho tất cả các 
nút trong miền tính toán cùng 
một thời gian (khác với sơ đồ 
hiện là giải từng bước một).
li ,1+
ix,
1, +li
li ,1−
1, −li
x∆
t∆
li,
1,1 ++ li1,1 +− li
1,1 −− li 1,1 −+ li
• PP này có bản chất luôn luôn ổn 
định. 
t
h
T
S
x
h
∂
∂=∂
∂
2
2
t
hh
T
S
x
hhh li
l
i
l
i
l
i
l
i-
∆
−=∆
+− +++++ 1
2
1
1
11
1 2
Chuyển từ pt vi phân sang 
pt sai phân hữu hạn
10/6/2013
12
Phương pháp sơ đồ ẩn
(Implicit Method)
hhShhh li
l
i
l
i
l
i
l
i- −+− +++++ 111111 2
tTx ∆=∆ 2
2x
t
S
Tr ∆
∆=
( ) lilililili- hhhhhxtST −=+−∆∆ +++++ 1111112 2
l
i
l
i
l
i
l
i- hrhhrrh =−++− ++++ 11111 )21(
Bước thời gian l+1
chưa biết
Bước thời gian l
Đã biết
Dòng chảy ổn định 2 chiều
(2‐D Steady‐State Flow)
• Phương trình cơ bản
hhh W ∂⎞⎛ ∂∂⎞⎛ ∂∂
jy,
)5,1( )52( )53( )54(
)4,5(
)5,1(Node No. Unknown heads
Known heads
• Dòng chảy đồng nhất, đẳng hướng, ko có 
giếng
02
2
2
2
=∂
∂+∂
∂
y
h
x
h
t
SQ
y
T
yx
T
x w
wyx ∂=±⎟⎟⎠⎜
⎜
⎝ ∂∂+⎟⎠⎜⎝ ∂∂ ∑=1
y∆
)4,1( )4,2( )4,3( )4,4(
)3,1( )3,2( )3,3( )3,4(
)2,1( )2,2( )2,3( )2,4(
)1,1( )1,2( )1,3( )1,4(
)0,1( )0,2( )0,3( )0,4(
, , ,
)4,0(
)3,0(
)2,0(
)1,0(
)4,5(
)3,5(
)2,5(
)1,5(
• Bình quân theo không gian (bình quân trong 
ô lưới)
0
22
2
11
2
11 =∆
+−+∆
+−
y
hhh
x
hhh i,j+i,ji,j-,ji+i,j,ji- ix,x∆
4
1111 i,j+i,j-,ji+,ji-
i,j
hhhh
h
+++=
10/6/2013
13
Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng
(2‐D Heterogeneous Anisotropic Flow)
⎞⎛⎞⎛
i+1/2
y
node (i,j) i-1/2
0=⎟⎟⎠⎜
⎜
⎝ ∂
∂
∂
∂+⎟⎠⎜⎝ ∂
∂
∂
∂
y
hT
yx
hT
x
yx
j+1 
j
j+1/2 
Qx,i+1/2 Qx,i-1/2 
Qy,j+1/2 
∆x
∆y 
cell
Tx và Ty là hệ số chuyển nước theo 
phương x và y
021212121 =∆
⎟⎟⎠
⎞⎜⎜⎝
⎛
∂
∂−⎟⎟⎠
⎞⎜⎜⎝
⎛
∂
∂
+∆
⎟⎠
⎞⎜⎝
⎛
∂
∂−⎟⎠
⎞⎜⎝
⎛
∂
∂
−+−+
y
y
hT
y
hT
x
x
hT
x
hT
/i.j
y
/i.j
y
.j/i
x
.j/i
x
Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng
(2‐D Heterogeneous Anisotropic Flow)
• Hệ số chuyển nước trung bình, điều hòa
⎞⎛⎞⎛ hhhh
0
1
21
1
21
1
21
1
21
=∆
⎟⎟⎠
⎞
⎜⎜⎝
⎛
∆
−−⎟⎟⎠
⎞
⎜⎜⎝
⎛
∆
−
+
∆
⎟⎟⎠⎜
⎜
⎝ ∆
−−⎟⎟⎠⎜
⎜
⎝ ∆
−
y
y
hh
T
y
hh
T
x
x
T
x
T
i,j-i,jy
/i,j-
i,ji,j+y
/i,j+
,ji-i,jx
,j/i-
i,j,ji+x
,j/i+
x
i,j
x
,ji+
x
i,j
x
,ji+x
,j/i+ +TT
TT
T
1
1
21 2=
01111 =++++ −+−+ i,ji,ji,ji,ji,ji,j,jii,j,jii,j hEhDhChBhA
10/6/2013
14
Các vấn đề chuyển đổi
(Transient Problems)
t
hS
y
hT
yx
hT
x
yx
∂
∂=⎟⎟⎠
⎞
⎜⎜⎝
⎛
∂
∂
∂
∂+⎟⎠
⎞⎜⎝
⎛
∂
∂
∂
∂
( )
t
hh
ShEhDhChBhA
l
i,j
l
i,j
i,j
l
i,ji,j
l
i,j-i,j
l
i,j+i,j
l
,ji-i,j
l
,ji+i,j ∆
−=++++
+
+++++
1
11
1
1
1
1
1
1
1
l
i,ji,j
l
i,ji,j
l
i,j-i,j
l
i,j+i,j
l
,ji-i,j
l
,ji+i,j hFhEhDhChBhA =′++++ +++++ 111111111
t
S
EE i,ji,jji ∆−=′,
MODFLOW
• Là mô hình toán được phát triển bởi Cục địa chất Mỹ 
(USGS – US Geological Survey)
• Sử dụng phương pháp sai phân hữu hạn
• Một số Version đã phát triển tham khảo: 
• Gia diện đồ họa MODFLOW:
– GWV ‐ Groundwater Vitas: ‐
vistas com/.  :
– GMS – Ground Modeling System: 
– PMWIN – Modflow for Window: 
10/6/2013
15
MODFLOW có thể mô phỏng được cái gì?
Mô hình MODFLOW có thể mô phỏng 
được những yếu tố sau:
1. Tầng ngậm nước không áp và có áp
2. Đứt gãy và các biên không thấm
3. Các địa tầng có áp hạn mịn và các lớp 
đáy liên kết
4. Đại tầng có áp ‐ dòng chảy ngầm và 
sự thay đổi lượng trữ
5. Sông – trao đổi nước ở các tầng 
ngậm nước
6. Lưu lượng nước từ các kênh tiêu và 
suối
7. Dòng chảy phù du – trao đổi nước ở 
á ầ ậ ớc c t ng ng m nư c
8. Hồ chứa – trao đổi nước ở các tầng 
ngậm nước
9. Bổ cập nước ngầm từ mưa và tưới
10. Bốc thoát hơi nước
11. Các giếng bơm hút và bổ cập nước 
ngầm
12. Xâm nhập mặn