10/6/2013
MÔ HÌNH NƯỚC NGẦM
Phần 1: Cơ sở lý thuyết
g y
g 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
Xác định và mô tả vấn đề
1. Xác định vấn đề cần mô phỏng
Cơ sở lý thuyết mô hình
Số liệu
– 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 – Độ chính xác
2. Xác định cơ sở lý thuyết và phát triển mô hình
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
Hiệu chỉnh mô hình & Đánh giá thông số
3. Hiệu chỉnh thông số mô hình
Kiểm định mô hình & Phân tích độ nhạy
– 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. được độ chính xác đặt ra
4. Kiểm định
Tổng hợp và hoàn chỉnh mô hì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)
Ứng dụng mô hình
– So sánh, đánh giá sự sai khác giữa kết quả tính toán và thực
đo
Trình diễn các kết quả
1
10/6/2013
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 ậ
www.epa.state.oh.us
– 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,
www.isws.illinois.edu
đ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 – Chuyển các phương trình vi phân từng 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.
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 giải thích các điều kiện địa chất và thủy văn. • Quá trình nào là
2
quan trọng đối với mô hình? • Các biên là gì? • Các giá trị thông số • 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
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
W
T
h
x
S
(
)
)
−
±
−
=
h 0
T x
y
xQ ( δ w
w
∑
K ' b '
∂ y ∂
h ∂ y ∂
∂ x ∂
h ∂ x ∂
h ∂ t ∂
w
=1
⎛ ⎜ ⎝
⎞ +⎟ ⎠
⎛ ⎜⎜ ⎝
⎞ +⎟⎟ ⎠
Dòng chảy
Thấm
Lượng trữ
Nguồn chảy vào/chảy ra
Bề mặt đất
Cột nước tầng có áp Cột nước tầng có áp
Lớp có áp (không thấm)
• Các phương trình chủ đạo Các phương trình chủ đạo • Các điều kiện biên • Các điều kiện ban đầu
h
Tầng có áp
Qx
b
z
K
y x
Tầng đá gốc
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
3
10/6/2013
Chuỗi Taylor
• Khai triển chuỗi Taylor đối với h(x) tại điểm
(x+∆x) lân cận điểm x ( ∆ ) lâ
ậ điể
n
2
(
(
x+h (
h
x )
++xh
n xh )(
′′ )(
)(
=
∆
′∆ +xhx+x )(
+
L
L
x ) ∆ ! n
x ) ∆ !2 • 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à: số của chuỗi số là n+ 1
n+
1
1
error
h
x )(
)
<
=
( ϑ
n +∆ x
max
x ( ) ∆ n+ ( )!1
Đạ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ị
2
( (
( ( x+h h
′ ′ )( )(
)( )(
) ) x
h
trí lân cận điểm x là: +xhx+xh h ∆ ⋅∆
=
′′ ′′ )( )( L+xh h
) ) x x ∆ ∆ !2 Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau:
x+h (
x )(
′ xh )(
=
≈
h x
∂ ∂
-hx ) ∆ x ∆
x∆
x∆
x
x
x ∆− x
x
x ∆+
4
10/6/2013
Đạ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à:
( (
( x-h
)(
′ )(
x ) =∆
+xhx-xh ⋅∆
′′ )( L-xh
2∆ 2 ) ) x ∆ !2
Trong đó đạo hàm bậc 1 được tính xấp xỉ như sau:
)
x
∆
′ )( xh xh )(
=
≈
h x
∂ ∂
)( ( -hxh x- x ∆
x∆
x∆
x
x
x
x ∆−
x
x ∆+
Đạo hàm bậc 2
2
(
( x+h x+h (
) ) x x
′ )( )(
∆ ∆
=
+xhx+xh )( )( h +xhx+x ⋅∆ ∆
′′ )( +xh L+xh )(
2
(
h
′′ xh )(
x+h (
x )
x )(
′ +xhx )(
∆
=
⋅∆−
L−
x ) ∆ !2 x ) ∆ !2
2
x+h (
(
)
∆
x ∆
′′ )( xh
=
≈
h ∂ 2 x ∂
x-+hxh-x )(2) 2 x ∆
5
Trong đó h”(x) là đạo hàm bậc 2 được tính xấp xỉ là:
10/6/2013
Xấp xỉ sai phân hữu hạn
)( xh
=
( ( xh h
) )
h i x =∆+
1 +
h h i h i
1 −
x∆
x∆
( xh ′ )( xh
) x =∆− ′= h i
x
h
i
i
−1
i
1 −
+
h
i
i =′
i =′ hi h
i′′ hi =h
-hh i i x ∆ Đạo hàm lùi bậc 1
-h +1 x ∆ Đạo hàm tiến bậc 1
1 2 +hh-h i i 2x 2 ∆ Đạo hàm trung tâm bậc 2
Lưới và rời rạc (Grids and Discretization)
y, j
Lưới Lưới
Vùng tính toán
i,j+1
∆y
i-1,j
i,j
i+1,j
• Quá trình rời rạc • p p Lưới được thiết lập để phủ kín vùng tính toán • Mục tiêu – dự báo mực nước tại
Điểm nút
i,j-1
x, i
∆x
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, …
– Sử dụng tốt cho vùng địa hình
đơn giản
Ô lưới
6
• Phương pháp sai phân hữu hạn – Được sử dụng phổ biến do tính Được sử dụng phổ biến do tính đơ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
10/6/2013
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ữ
• Lưới được tổ chức bởi các hàng (i), cột (j), lớp (k), mỗi một
(i) ộ (j) lớ (k)
ỗi ộ
nhật bởi các lưới. L ới đ ổ hứ bởi á hà 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)
Bề mặt đất (Ground surface)
• Dòng chảy đồng nhất,
Lớp không thấm
hA
Tầng ngậm nước (Aquifer)
Nút (Node)
∆x
hB
T
S
b
∂ x ∂
h ∂ x ∂
h ∂ t ∂
⎛ ⎜ ⎝
⎞ =⎟ ⎠
z
y
x
0 i = 0 i
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
Ô lưới (Grid Cell)
đẳng hướng, 1 chiều, và p có áp. • Phương trình mô tả cho dòng chảy xác định ở trên:
1.6
m
=
• Điều kiện ban đầu (initial diti ) condition): ( )0, xh
• Điều kiện biên (boundary
∆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
1.6 5.1
m m
= =
7
conditions): ),0( h t tLh ( ),
10/6/2013
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:
lt,
2
=
, +li 1
h ∂ 2 x ∂
h ∂ t ∂
t∆
S T • Đạo hàm bậc 2 đối với (WRT with respect
i
l
,1−
i
l
,1+
li,
to) x
l
ix,
2
x∆
≈
2
h ∂ 2 x ∂
l l l 2 +hh-h 1 1 i i i + − x ∆
i
, li , −li 1 1
• Đạo hàm bậc 1 đối với (WRT) t
l
+1
l
−1
l h i
l i
≈
≈
-h t
h ∂ t ∂
∆
i
l l -hh i i l ∆
h ∂ t ∂
i lùi
Tiến
Phương pháp sai phân sơ đồ hiện (Explicit Method)
lt,
1
, +li
•
t∆
li,
i
l
,1−
i
l
,1+
•
ix,
x∆
1
, −li
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 ở 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
1 +
2
ỏ ể
+
l h i
l h i
l h i 1 +
=
=
2
S T
l l 1 2 h h − i- i x ∆
− t ∆
S T
h ∂ 2 x ∂
h ∂ t ∂
(∆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, 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.
8
Từ pt vi phân chuyển sang sai phân (xấp xỉ)
10/6/2013
Phương pháp sai phân sơ đồ hiện (Explicit Method)
1 +
+
l h i
l h i
l h 1 i +
=
2
S T T
l l h 1 2 h − i- i x x ∆ ∆
− t t ∆ ∆
1
+
(
2
)
−
+
=
−
l h i-
1
l h i
1
l h i
l h i
l h i +
T S
t ∆ 2 x ∆
r
=
T S
−
+
t ∆ 2x ∆ ( ( l l hr h 1 i-
1 l l h 1 + + h i
l l h h += i
l l 2 h h 2 i
Đặt:
)l )l h h i
1 +
Bước thời gian l Đã biết Bước thời gian l+1 Chưa biết
Phương pháp sai phân sơ đồ hiện (Explicit Method)
Bề mặt đất (Ground surface)
1
+
h
h
hr (
2
h
h
)
=
+
−
+
l i
l i
l i-
1
l i
1
l i +
Lớp không thấm
hA
Tầng ngậm nước
Nút tính toán
∆x
hB
r
=
b
T S
t ∆ 2 x ∆
2
t
xr ∆=∆
S T
i = 0
1
2
3
4
5
6
7
8
9
10
Ô lưới Ô lưới
2
xr
t
.0
0128
d
5.18
min
∆=∆
=
≈
S T
Xem xét: r = 0.48
.0
0139
d
20
min
t =∆
≈
r = 0.52
9
∆x= 1m, L = 10m, b = 1.5m hA = 6.1m, hB = 1.5m, K = 0.5m/d, S = 0.02
10/6/2013
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
10/6/2013
Phân tích quá trình không ổn định đó diễn ra như thế nào?
Bề mặt đất
• Tại thời gian t = 0 (cid:198) ko có dòng chảy
Lớp ko thấm
hA
∆x
Tầng ngậm nước
hB
• Khi t > 0 (cid:198) có dòng chảy •
)1(
V
h
∆⋅
∆⋅
∆
=
1
b
∆x
S • Nước chảy ra khỏi một ô lưới trong thời khoảng ∆t là:
i = 0
1
2
… i-1
i
i+1 … 8
9
10
V
T
t
∆
=
∆
2
h ∆ x∆ 2/ x 2/ ∆
Ô lưới i
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: x
x
T
St
⋅
h ∆⋅∆⋅≤∆⋅
• Theo nguyên lý bảo toàn khối lượng: ∆V2 ≤ ∆V1
r
=
1 2
h ∆ x 2/ ∆ T t ∆ 2 ≤ S x ∆
Vậy nếu r > 0.5 (cid:198) khoảng thời gian qua lớn (cid:198) ô lưới nước ngầm không thể chứa đủ nước (cid:198) không ổn định
Phương pháp sai phân sơ đồ ẩn (Implicit Method)
t,
l
i
1
+ l
,1 +
i
1
− l
,1 +
i
1
, +l
t∆
• Sử dụng thông tin (H, Q) ở một nút thời gian trước tính toán nút thời gian trước tính toán cho 3 nút thời gian sau (sơ đồ ẩn).
i,
l
i
l
,1−
i
l
,1+
• Giải đồng thời cho tất cả các
ix,
x∆
i
1
− l
,1 −
i
1
+ l
,1 −
i
1
, −l
ệ
g
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). ộ ) g • PP này có bản chất luôn luôn ổn
định.
1 +
1 +
2
−
+
1 l + h 1 i-
l h i
l h i
1 l + h 1 i +
=
=
2
S T
h ∂ 2 x ∂
h ∂ t ∂
S T
l 2 h i x ∆
− t ∆
11
Chuyển từ pt vi phân sang pt sai phân hữu hạn
10/6/2013
Phương pháp sơ đồ ẩn (Implicit Method)
1 +
1 +
−
+
l 1 + h 1 i- i 1
l h i i
l 1 + h 1 i 1 i + +
=
2
l hS i i T
l 2 h i i x ∆
− t ∆
1
1
1
+
+
2
−
+
=
−
)
( l h i-
+ 1
l h i
l h i
l h i
1 l + h 1 i +
T S
t ∆ 2 x ∆
r
=
T S
t ∆ 2x ∆
1
1
+
rh
hr
rh
h
−
+
)21( +
−
=
l i-
+ 1
l i
l i
1 l + 1 i +
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)
Node No.
jy,
Unknown heads
Known heads
)5,1( )5,1(
W W
)5,2( )52(
)5,1( )4,5( )5,3( )53(
)5,4( )54(
T
T
Q
S
=
x
y
w
∑
)4,0(
)4,1(
)4,2(
)4,3(
)4,4(
)4,5(
∂ ∂ x ∂
h h ∂ ∂ x ∂
∂ ∂ y ∂
h h ∂ ∂ y ∂
h h ∂ ∂ t ∂
w
=1
⎛ ⎛ ⎜ ⎝
⎞ ⎞ +⎟ ⎠
⎞ ⎞ ±⎟⎟ ⎠
⎛ ⎛ ⎜⎜ ⎝
)3,0(
)3,5(
)3,1(
)3,2(
)3,3(
)3,4(
• Phương trình cơ bản
)2,0(
)2,1(
)2,5(
)2,2(
)2,3(
)2,4(
2
2
0
+
=
)1,0(
)1,1(
)1,2(
)1,3(
)1,4(
)1,5(
h ∂ 2 x ∂
h ∂ 2 y ∂
)0,1(
)0,2(
)0,3(
)0,4(
y∆ • Dòng chảy đồng nhất, đẳng hướng, ko có giếng
2
2
−
+
−
+
h i-
,j
h i+
,j
h i,j-
h i,j+
1
1
1
1
0
+
=
x∆
h i,j 2 x
h i,j 2 y
∆
∆
ix,
h
h
h
h
+
+
+
i-
1
,j
i+
1
,j
i,j-
1
i,j+
1
h
=
i,j
4
12
• Bình quân theo không gian (bình quân trong ô lưới)
10/6/2013
Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng (2‐D Heterogeneous Anisotropic Flow)
y
node (i,j)
i+1/2
i-1/2
∆x
x
y
T
T
cell
∂ x ∂
h ∂ x ∂
∂ y ∂
h ∂ y ∂
⎛ ⎛ ⎜ ⎝
⎞ ⎞ +⎟ ⎠
⎛ ⎛ ⎜⎜ ⎝
⎞ ⎞ 0=⎟⎟ ⎠
j+1
Qy,j+1/2
j+1/2
Qx,i-1/2
Qx,i+1/2
∆y
j
y
y
x
x
T
T
−
T
T
−
h y
h y
∂ ∂
∂ ∂
h x
∂ ∂
h ∂ x ∂
⎛ ⎜ ⎝
⎞ ⎟ ⎠
⎞ ⎟ ⎠
⎛ ⎜⎜ ⎝
⎞ ⎟⎟ ⎠
⎞ ⎟⎟ ⎠
i
/ 21
.j
i
/ 21
.j
i.j
/ 21
i.j
/ 21
+
−
+
−
0
+
=
⎛ ⎜ ⎝ x
⎛ ⎜⎜ ⎝ y
∆
∆
Tx và Ty là hệ số chuyển nước theo phương x và y
Dòng chảy 2 chiều, ko đồng nhất, ko đẳng hướng (2‐D Heterogeneous Anisotropic Flow)
h h
h h
h h
h h
−
1
1
i,j
i-
,j
i+
i,j
T
T
1
2
x 1
2
i-
/
,j
x i+
/
,j
x
− x
,j ∆
∆
⎞ ⎞ ⎟⎟ ⎠
⎛ ⎛ ⎜⎜ ⎝
⎛ ⎛ ⎜⎜ ⎝
⎞ ⎞ −⎟⎟ ⎠ ∆
h
h
h
h
−
i,j
i,j-
1
i,j+
i,j
T
T
1
2
1
2
y i,j-
/
y i,j+
/
− y
y
1 ∆
∆
⎛ ⎜⎜ ⎝
⎛ ⎜⎜ ⎝
⎞ ⎟⎟ ⎠
0
+
=
y
x ⎞ −⎟⎟ ⎠ ∆
T
2=
x i+
21 /
,j
T
x x T T 1 i+ i,j ,j x x +T 1 i+ i,j
,j
0
+
+
+
+
=
hA i,j
i
1
,j
hB i,j
i
1
,j
hC i,j
i,j
1
hD i,j
i,j
1
hE i,j
i,j
+
−
+
−
13
• Hệ số chuyển nước trung bình, điều hòa
10/6/2013
Các vấn đề chuyển đổi (Transient Problems)
x
y
T
T
S
∂ y y ∂
∂ x ∂
h ∂ x ∂
h ∂ y y ∂
h ∂ t ∂
⎛ ⎜ ⎝ ⎝
⎞ +⎟ ⎠ ⎠
⎛ ⎜⎜ ⎝ ⎝
⎞ =⎟⎟ ⎠ ⎠
1
−
)
( l + h i,j
l h i,j
1
S
+
+
+
+
=
hA i,j
1 l + 1 ,j i+
l hB i,j i-
+ 1
1 ,j
hC i,j
1 l + 1 i,j+
hD i,j
1 l + 1 i,j-
l + hE i,j i,j
i,j
t
∆
1
+
+
+
=
hA i,j
1 l + 1 ,j i+
hB i,j
l i-
+ 1
1 ,j
hC i,j
1 l + 1 i,j+
hD i,j
l + i,j-
1 1
′+ hE i,j
l + i,j
hF i,j
l i,j
S
E
E
−
=′,
j
i
i,j
i,j t ∆
MODFLOW
• Là mô hình toán được phát triển bởi Cục địa chất Mỹ
y) (USGS – US Geological Survey) (
g • 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: http://water.usgs.gov/nrp/gwsoftware/modflow.html • Gia diện đồ họa MODFLOW:
– GWV ‐ Groundwater Vitas: http://www.groundwater‐
vistas.com/ : vistas com/
– GMS – Ground Modeling System: http://www.aquaveo.com/ – PMWIN – Modflow for Window:
http://www.ifu.ethz.ch/publications/software/pmwin/index_EN
14
10/6/2013
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. 2. 3.
4.
5.
6.
7.
ậ
ầ
8.
9. 10. 11.
12.
Tầng ngậm nước không áp và có áp Đứt gãy và các biên không thấm Các địa tầng có áp hạn mịn và các lớp đáy liên kết Đại tầng có áp ‐ dòng chảy ngầm và sự thay đổi lượng trữ Sông – trao đổi nước ở các tầng ngậm nước Lưu lượng nước từ các kênh tiêu và suối Dòng chảy phù du – trao đổi nước ở các tầng ngậm nước ớ á Hồ chứa – trao đổi nước ở các tầng ngậm nước Bổ cập nước ngầm từ mưa và tưới Bốc thoát hơi nước Các giếng bơm hút và bổ cập nước ngầm Xâm nhập mặn
15
