Mô hình nước ngầm - Phần 1: Cơ sở lý thuyết

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