Luận văn Đánh giá mức độ xâm nhập mặn do sử dụng quá mức nước ngầm trong vùng ven biển nam định bằng mô hình visual modflow

DĐ theo độ tổng khoáng hoá thể hiện ở hình 10. Độ tổng khoáng hoá của nước dưới đất biến đổi trong một khoảng rộng từ nhạt đến mặn. Ranh giới mặn nhạt (M=1g/l) xác định được khá rõ ràng, nội suy từ kết quả thí nghiệm giữa 2 lỗ khoan và kết quả đo địa vật lý. Độ tổng khoáng hoá của nước ở 35 đây biến đổi từ 0,29 đến 1g/l. Thành phần hoá học của nước biến đổi theo chiều tăng dần của độ tổng khoáng hoá như sau: - Khi M < 0,5g/l nước có thành phần bicacbonat clorua-natri đôi khi natri-canxi. - Khi M > 0,5g/l nước có thành phần clorua bicacbonat -natri Tức là khi độ tổng khoáng tăng thì hàm lượng clo tăng. - Diện phân bố của các vùng có độ tổng khoáng hoá lớn hơn 1g/l nằm ở phía bắc, phía đông nghiên cứu, trong đó lớn hơn 3g/l chiếm diện tích rất nhỏ ở phía đông bắc. Vùng nước có độ tổng khoáng hoá < 1g/l nằm ở phía tây nam, trong đó vùng giáp biển thuộc các huyện Hải Hậu, Nghĩa Hưng vẫn là nước nhạt. Diện phân bố của vùng nước nhạt trong tầng qp ngược với các tầng chứa nước trên nó. Sự tồn tại của vùng nước nhạt khá rộng dãi và sát đến tận bờ biển được tạm giải thích là có nguồn cung cấp từ xa theo các đứt gãy kiến tạo, các đới Karst hoá của các đá cổ cung cấp lên đẩy được nước mặn ra xa để tạo thành lưỡi nhạt khá lớn. Vùng nước lợ phân bố ở phía đông bắc tạo thành dải hẹp theo phương tây bắc- đông nam tạo thành vùng đệm gồm nước nhạt và nước mặn. Thành phần hoá học chuyển từ clorua bicacbonat sang clorua-natri. Vùng nước mặn phân bố ở góc đông bắc vùng nghiên cứu và còn phổ biến ở rộng dãi ở tỉnh Thái Bình. Thành phần hoá học ở đây là clorua-natri. - Nồng độ một số chất khác có trong nước: + Nồng độ sắt từ rất nhỏ đến 22mg/l, có khoảng 50% số mẫu có hàm lượng sắt >0,5mg/l. + Độ ôxy hoá từ 2 đến 12mg/l, trong đó 50% số mẫu có độ ôxy hoá lớn hơn 4mg/l. + Hàm lượng các chất độc hại, kim loại nặng, chất phóng xạ, các hợp chất nitơ đều thấp hơn tiêu chuẩn cho phép, một số mẫu ở một số lỗ khoan hàm lượng vi sinh lớn hơn tiêu chuẩn cho phép. 5. Tầng chứa nước Pliocen (n2) 36 - Độ tổng khoáng hoá nước dưới đất ở các lỗ khoan nghiên cứu biến đổi từ 0,58 (LK54) đến 16g/l (LK58) song do số lượng điểm nghiên cứu còn quá ít nên chưa vẽ được chính xác ranh giới mặn nhạt mà chỉ sơ bộ nhận định được là vùng phân bố nước nhạt rộng tương đối phù hợp với tầng chứa nước qp. - Thành phần hoá học của nước ở vùng nước nhạt là bicacbonat clorua natri còn ở vùng nước lợ và mặn là clorua natri. Hàm lượng các nguyên tố vi lượng, các chất độc hại đều nằm dưới tiêu chuẩn cho phép đối với nước ăn uống sinh hoạt. Tuổi của nước được xác định là 18 ± 5 năm. 6. Tầng chứa nước các trầm tích cacbonat trias giữa (t2) - Độ tổng khoáng hoá 0,46g/l thành phần là bicacbonat clorua -natri canxi. 7. Tính phân đới thuỷ địa hoá và nguồn gốc của nước dưới đất Vùng nghiên cứu trên mặt cắt địa chất thuỷ văn tồn tại 3 tầng chứa nước lỗ hổng (qh2; qh1; qp); 1 tầng chứa nước lỗ hổng- khe nứt vỉa (n2) và 1 tầng chứa nước khe nứt (t2). Tầng chứa nước trên cùng chủ yếu là nhạt, tầng chứa nước tiếp theo (qh1) bị mặn hoàn toàn. Hai tầng chứa nước tiếp theo là qp và n2 phần lớn diện tích là nhạt. Với đặc điểm như vậy có thể xếp vào kiểu phân đới thuỷ hoá nghịch tương đối phổ biến ở vùng ven biển nước ta. Về nguồn gốc của nước dưới đất: - Hai tầng chứa nước trên cùng rõ ràng có nguồn gốc từ khí quyển: được nước mưa, nước tưới cung cấp cho ở trên toàn bộ diện phân bố của tầng. - Nguồn gốc của nước tầng chứa nước qp và n2 hiện nay có nhiều giả thuyết. Giả thuyết có từ sớm nhất là nguồn gốc chôn vùi không đủ sức thuyết phục do đã khá rõ về đặc điểm địa chất thuỷ văn. Đa số ý kiến của nhà khoa học đều cho rằng có nguồn cung cấp từ xa và cung cấp từ dưới lên. Tôi cũng thống nhất ý kiến này với bằng chứng là kết quả phân tích tuổi của nước tầng qp là 25 ± 4 năm của tầng n2 là 18 ± 5 năm. Phân tích cấu trúc địa chất khu vực tôi cho rằng miền cung cấp cho nước của tầng qp và n2 vùng nghiên cứu là dải đá vôi t2a Ninh Bình chạy theo hướng tây bắc-đông 37 nam và dải đá vôi Hoà Bình, Hà Nam cũng t2a chạy theo hướng bắc tây bắc-nam đông nam. Nước vận động theo các đứt gãy kiến tạo, hệ thống Karst hoá của đá vôi đến vùng nghiên cứu thì đá vôi này chìm xuống, nước trong đó có áp lực cao hơn cung cấp cho các tầng chứa nước nằm trên.  Tóm lại đặc điểm thuỷ hoá vùng ven biển Nam Định có thể nói một nét khái quát là tồn tại các thấu kính, các đới nước nhạt xen kẽ nước mặn. 2.2. Phương pháp nghiên cứu 2.2.1. Phương pháp kế thừa số liệu Khu vực nghiên cứu có 2 tầng chứa nước đã được nghiên cứu khá chi tiết trong các đề tài nghiên cứu trong đó có đề tài “Xác định tiêu chí áp dụng tổ hợp các phương pháp địa chất thuỷ văn, địa vật lý, mô hình số để điều tra, đánh giá nhiễm mặn và tìm kiếm các thấu kính hoặc tầng chứa nước nhạt dải ven biển Nam Định” và có ý nghĩa quan trọng trong khai thác nước phục vụ cho nhu cầu sinh hoạt. Đó là tầng chứa nước lỗ hổng Holocen và tầng chứa nước Pleistocen. Chính vì vậy tôi đã sử dụng phương pháp kế thừa số liệu: kế thừa những kết quả điều tra ngoài thực địa và một số hoạt động xử lý trong phòng của các đề tài có liên quan để lấy số liệu phục vụ cho luận văn. Các số liệu được kế thừa từ Các phương pháp trong đề tài đã nghiên cứu bao gồm: Phân tích ảnh viễn thám Công tác phân tích ảnh viễn thám tức là ảnh chụp từ máy bay hoặc ảnh vệ tinh được thực hiện để tìm kiếm các thấu kính, tầng chứa nước nhạt vùng ven biển, có thể áp dụng có hiệu quả ở các dạng địa hình cao như cồn cát, dải cát, các đồi núi thấp ven biển, các nội dung chính cần phải thực hiện là: -Khoanh định các dải cát, cồn cát ven biển để xác định diện tồn tại các thấu kính, tầng nước nhạt. -Xác định cấu trúc địa chất để làm rõ miền cung cấp, hướng vận động của nước dưới đất. 38 Phương pháp thực địa Để có được nguồn số liệu để phục vụ cho công tác xây dựng và chỉnh lý mô hình cần kết quả của nhiều phương pháp khảo sát thực địa: 1. Phương pháp khảo sát thực địa Khảo sát thực địa được thực hiện tương tự như trong công tác điều tra đánh giá nước dưới đất nói chung được thực hiện ở tất cả diện tích nghiên cứu. Đối với công tác tìm kiếm các thấu kính, tầng nước nhạt vùng ven biển cần phải chú ý các vấn đề sau đây: -Quan sát mô tả địa hình, địa mạo để xác định hướng phát triển, kích thước, độ cao các cồn cát, dải cát, đụn cát. -Thành phần đất đá cấu tạo nên các cồn cát, dải cát, đụn cát. -Thảm thực vật phủ lên các dải cát, cồn cát, đụn cát trong đó cần phân định các loại thực vật ưa nước lợ, nước mặn. -Tình hình khai thác sử dụng nước dưới đất của nhân dân địa phương. 2. Công tác địa vật lý Xác định các thấu kính, tầng chứa nước nhạt gắn liền với xác định các cấu trúc địa chất, trong đó có việc phân chia mặt cắt địa chất ra các tầng lớp đất đá có thành phần thạch học khác nhau căn cứ vào điện trở suất. Sau đó việc tìm kiếm xác định các thấu kính tầng chứa nước nhạt được thực hiện trong các tầng chứa nước riêng biệt. -Đối với các tầng chứa nước nông (cồn cát, dải cát ven biển) thích hợp nhất để nghiên cứu cấu trúc địa chất là đo sâu điện trở có thiết bị phù hợp với chiều sâu nghiên cứu, ví dụ ABmax áp dụng trong đề tài là 300m có thể nghiên cứu đến chiều sâu khoảng 50. Các thấu kính nước nhạt thường tồn tại trong độ sâu ≤ 20m nguồn cung cấp cho chúng là nước mưa. Nước mặn thường xâm nhập theo chiều ngang từ phía biển hoặc vùng cửa sông tạo thành các nêm nước mặn. Áp dụng phương pháp Georada có thể giải quyết tốt nhiệm vụ xác định các thấu kính nước nhạt. 39 Phương pháp Georada đã được các nước phát triển trên thế giới áp dụng và đạt được những kết quả tốt. Ở nước ta cũng có 1 số nơi thực hiện đem lại kết quả tốt. Để phát hiện các thấu kính nước nhạt vùng ven biển cần bố trí tuyến đo Georada vuông góc với nguồn gây mặn như bờ biển. Các tuyến và điểm đo hình thành mạng lưới đo có độ dày phụ thuộc vào mức độ chi tiết của yêu cầu nghiên cứu. -Đối với các tầng chứa nước nằm sâu, việc xác định cấu trúc cũng được thực hiện bằng phương pháp điện trở có thiết bị phù hợp với chiều sâu nghiên cứu, tuy nhiên phương pháp này chỉ có thể áp dụng có hiệu quả đến chiều sâu 160m. Việc áp dụng phương pháp đo sâu điện trở được áp dụng với tổ hợp các phương pháp đo sâu điện từ (trường chuyển) phương pháp đo sâu phân cực kích thích. Để xác định độ tổng khoáng hoá của nước cần phải xác định mối quan hệ giữa độ tổng khoáng hoá và điện trở suất ςfK ở vùng nghiên cứu trên cơ sở đo karota hoặc đo tham số ở các lỗ khoan đã có. Trên cơ sở nghiên cứu ở vùng đồng bằng Bắc Bộ đã xác định mối quan hệ này như hình vẽ cho thấy ranh giới mặn nhạt ở chỗ có điện trở suất khoảng 10Ωm. Để nghiên cứu cấu trúc địa chất nên bố trí các tuyên đo vuông góc với phương cấu trúc, còn xác định ranh giới mặn nhạt trong tầng chứa nước thì bố trí tuyến đo vuông góc với ranh giới mặn nhạt. Độ dày của tuyến đo và điểm đo phụ thuộc vào mức độ chi tiết nghiên cứu được thiết kế cụ thể trong các đề án. Một biện pháp hữu hiệu để xác định ranh giới mặn nhạt theo chiều thẳng đứng là phương pháp đo địa vật lý lỗ khoan (karota). Ở các vùng có nhiều tầng, nhiều lớp chứa nước, điều kiện thuỷ địa hoá phức tạp, các lỗ khoan sau khi khoan xong nhất thiết phải đo karota mới thiết kế chống ống, ống lọc vào đoạn chứa nước nhạt. 3. Phương pháp khoan khai đào Công tác khoan khai đào được thực hiện thành các tuyến theo hướng vuông góc với dải cát, cồn cát vuông góc với ranh giới mặn nhạt. Mục đích của khoan khai đào là 40 xác định chiều dày của tầng chứa nước, xác định thành phần vật chất của đất đá chứa nước và tạo điểm để làm công tác thí nghiệm thấm. Vị trí khoan được xác định sau khi đo địa vật lý và nhất định phải đặt ở 1 điểm đo địa vật lý. Công tác chống ống kết cấu lỗ khoan được xác định sau khi đo karota lỗ khoan. 4. Phương pháp thí nghiệm thấm Công tác thí nghiệm thấm gồm có đổ nước trong hố đào, bơm nước thí nghiệm trong giếng đào, trong lỗ khoan. Công tác đổ nước thí nghiệm được thực hiện trong đới thông khí để xác định tính thấm của đất đá từ đó xác định khả năng thấm từ nước mưa cung cấp cho tầng chứa nước. Công tác bơm thí nghiệm được thực hiện để xác định các thông số địa chất thuỷ văn và lấy mẫu nước xác định độ mặn. 5. Công tác lấy mẫu phân tích hoá học và độ mặn Độ mặn của nước cần được xác định ở tất cả các điểm nước khảo sát bằng các dụng cụ đo độ mặn (Multiline P4 hoặc SenSion 5). Các điểm nước điển hình lấy mẫu để phân tích Ion Clo trong đó 1 số sẽ lấy mẫu để phân tích độ tổng khoáng hoá bằng cách xấy khô. Một số điểm cần lấy mẫu phân tích đồng vị, phân tích xác định tuổi của nước, nguồn gốc của nước dưới đất. 6. Phân tích cấu trúc và lịch sử phát triển địa chất Được thực hiện nhằm nội suy xác định nguồn cung cấp cho nước dưới đất, nguồn gốc nước dưới đất, phương vận động và tốc độ vận động dòng chảy dưới đất. 7. Phân tích mẫu nước ngoài thực địa và trong phòng Công tác phân tích mẫu nước được tiến hành ngoài thực địa và trong phòng thí nghiệm 41 Phân tích ngoài hiện trường được thực hiện ở tất cả các điểm khảo sát, các lỗ khoan. Các chỉ tiêu xác định là: nhiệt độ, độ pH; độ muối (đo bằng dụng cụ SenSion 5). Phân tích trong phòng. Các chỉ tiêu phân tích là: Các đại nguyên tố: Na+; K+; Ca+2; Mg+2; NH4+; Cl-; SO4-; HCO3-; CO3-; NO3-; No2-. Các chất khác: Cần sấy khô ở 105oC; 2.2.2. Phương pháp thống kê số liệu và xử lý số liệu Từ kết quả của các công tác khảo sát thực địa và công tác trong phòng để thống kê các loại số liệu phục vụ đầu vào cho mô hình: Cao độ và bề dày các tầng chứa nước thông qua các lỗ khoan khảo sát trong vùng nghiên cứu: Hình 11: Mạng lưới lỗ khoan quan trắc trong vùng nghiên cứu và vùng lân cận phục vụ điều tra địa chất và địa chất thủy văn 42 Mạng lưới lỗ khoan khảo sát địa chất và địa chất thủy văn trong vùng nghiên cứu và vùng lân cận phục vụ cho xây dựng mô hình Kết quả thu thập về, tọa độ vị trí các lỗ khoan nghiên cứu, cao độ địa hình, cao độ đáy, bề dày các tầng chứa nước tại vị trí các lỗ khoan được bố trí trong hình 11 được thể hiện trong bảng 6 phần phụ lục. Số liệu về độ tổng khoáng hóa, độ muối tại vị trí các điểm, lỗ khoan được thể hiện trong bảng 7 và 8 trong phần phụ lục. 2.2.3. Phương pháp mô hình số a, Mô hình Visual MODFLOW Hệ thống phần mềm mô hình Visual Modflow do công ty Waterloo Hydrogeologic Inc. (WHI) phát triển từ những năm 80 của thế kỷ trước. Cùng với sự lớn mạnh không ngừng của công nghệ máy tính và công nghệ phần mềm, hệ thống này cũng liên tục phát triển. Cho đến nay, hệ thống chương trình, mô hình mô phỏng và tính toán các quá trình nước dưới đất đã đạt đến một trình độ tương đối hoàn hảo. Các phiên bản mới nhất cho phép người sử dụng dễ dàng thiết lập được các mô hình NDĐ và quản lý chúng một cách hiệu quả. Các giao diện đầu vào, xử lý, và xuất dữ liệu của mô hình vì thế mà cũng thân thiện hơn với hầu hết các nhà nghiên cứu địa chất thủy văn. MT3D là một modul được phát triển lần đầu do nhóm chuyên gia Chunmiao Zheng, P. Patrick Wang, nó được tạo ra gần như song song với sự ra đời của phần mềm Modflow. MT3D là một mô hình số, ba chiều dùng để mô phỏng sự lan truyền vật chất trong nước dưới đất. Nó có thể mô phỏng được các quá trình đối lưu, phân tán, và khuyếch tán. b, Một số công thức tính toán được sử dụng trong mô hình Toàn bộ sự biến thiên độ cao mực nước duới đất được mô tả bằng một phương trình đạo hàm riêng duy nhất sau: 43 t hSsW z hT zy hT yx hT x zzyyxx ∂ ∂ =−      ∂ ∂ ∂ ∂ +      ∂ ∂ ∂ ∂ +      ∂ ∂ ∂ ∂ (1) Trong đó: - Txx, Tyy, Tzz là các hệ số dẫn nước theo phương x,y và z. Chiều z là chiều thẳng đứng. Ta có: T = K.m. (m: bề dày tầng chứa nước) - h là cốt cao mực nước tại vị trí (x,y,z) ở thời điểm t. - W là các giá trị bổ cập hay thoát của nước nước dưới đất tại vị trí (x,y,z) ở thời điểm t. W = W(x,y,z,t) là hàm số phụ thuộc thời gian t và không gian (x,y,z) - Ss là hệ số nhả nước - Ss = Ss(x,y,z), Kxx = Kxx(x,y,z), Kyy = Kyy(x,y,z), Kzz = Kzz(x,y,z) các hàm phụ thuộc vào vị trí không gian x,y,z (K: hệ số thấm của tầng chứa nước theo các chiều x, y , z) Phương trình (1) mô tả động thái mực nước trong điều kiện môi trường không đồng nhất và dị hướng. Phương trình (1) cùng với các điều kiện biên, điều kiện ban đầu của tầng chứa nước tạo thành một mô hình toán học về dòng chảy NDĐ. * Phương pháp giải Để giải phương trình trên, người ta phải tìm hàm số h(x,y,z,t), thoả mãn (1) và thoả mãn các điều kiện biên. Sự biến động của giá trị h theo thời gian sẽ xác định bản chất của dòng chảy, từ đó có thể tính được trữ lượng của tầng chứa nước cũng như tính toán các hướng của dòng chảy. Việc tìm lời giải giải tích h(x,y,z,t) của phương trình (1) chỉ khi nào miền nghiên cứu đưa về được mô phỏng bằng sơ đồ toán học. Thực tế, miền thấm có điều kiện rất phức tạp do đó người ta buộc phải giải bằng phương pháp gần đúng. Một trong các phương pháp giải gần đúng được áp dụng rộng rãi là phương pháp sai phân hữu hạn. 44 Khi áp dụng phương pháp này không gian nghiên cứu được phân ra hay rời rạc hóa thành nhiều ô. Ở mỗi ô, các giá trị tham gia vào phương trình được coi là không đổi. Giá trị này xấp xỉ với giá trị thực tế. Kết quả h(x,y,z,t) sẽ là một lưới ô các giá trị h. Bằng cách này người ta đưa phương trình đạo hàm riêng (1) về một hệ phương trình tuyến tính. Số lượng phương trình tương đương với số các ô lưới chia. Rõ ràng nếu bước lưới càng nhỏ thì kết quả thu được từ lời giải sai phân càng gần với lời giải đúng của phương trình (1). Thế nhưng khối lượng tính toán sẽ nhiều lên gấp bội, nên người ta phải tìm cách chọn ra bước lưới thích hợp. Nếu trong mỗi ô các giá trị tham gia tính toán trong phương trình không thay đổi đáng kể thì phép chia ô là hợp lý. Để hình dung được phương pháp sai phân áp dụng như thế nào, ta sẽ bắt đầu từ quá trình rời rạc hoá Hình 12: Ô lưới và các loại ô trong mô hình (tầng chứa nước K) 45 Hình 12 mô tả quá trình rời rạc hoá không gian. Không gian nghiên cứu được phân theo chiều thẳng đứng z thành các lớp chứa nước. Mỗi lớp chứa nước lại được chia thành các ô nhỏ hơn. Vùng hoạt động của nước dưới đất trong mỗi tầng chứa nước sẽ được đánh dấu là “ô trong miền tính”. Những ô cách nước hoặc không có dòng chảy thấm qua thì được đánh dấu là “ô ngoài miền tính”. * Phương trình vi phân Hệ phương trình sai phân nhận được từ phương trình (1) được thành lập trên cơ sở phương trình vi phân cân bằng của Buxines: Tổng dòng chảy đến và chảy đi từ một ô phải bằng sự thay đổi thể tích nước có trong ô. Giả thiết rằng khối lượng riêng của nước dưới đất là không đổi thì qui tắc cân bằng dòng chảy cho một ô được thể hiện bằng phương trình sau: V t hSsQ i i ∆∆ ∆ =∑ (2) - Qi là lượng nước chảy vào ô (nếu chảy ra thì Q lấy giá trị âm) - Ss là giá trị của hệ số nhả nước, nó chính là giá trị Ss(x,y,z) - ∆V là thể tích ô - ∆h là giá trị biến thiên của h trong thời gian ∆t tại ô lưới đang xét. 46 Hình vẽ mô tả cho một ô lưới (i,j,k) và 6 ô bên cạnh nó, (i-1,j,k), (i+1,j,k), (i,j- 1,k), (i,j+1,k), (i,j,k-1), (i,j,k+1) dòng chảy từ ô (i,j,k) sang các ô bên cạnh (nếu chảy vào mang dấu dương, chảy ra mang dấu âm). i,j,k i,j,k-1 i,j,k+1 i+1,j,k i,j-1,k i,j+1,k i-1,j,k Hình 13: Ô lưới i,j,k và 5 ô bên cạnh Nếu đặt Cri,j-1/2,k là sức cản thấm trong hàng thứ i, lớp thứ k giữa các nút lưới (i,j-1,k) và (i,j,k) được tính theo công thức: Cri,j-1/2,k=Kri,j-1/2,k∆ci∆vk/∆rj-1/2 (3) Trong đó: - Kri,j-1/2,k là hệ số thấm giữa các nút lưới (i,j,k) và (i,j-1,k) - ∆ci∆vk là diện tích bề mặt vuông góc với phương dòng chảy - ∆rj-1/2 là khoảng cách giữa các nút lưới (i,j,k) và (i,j-1,k) Và đặt lưu lượng cung cấp cho ô lưới từ biên theo phương trình tổng quát sau: ai,j,k,n = pi,j,k,n hi,j,k + qi,j,k,n (4) Trong đó: - ai,j,k,n biểu diễn dòng chảy từ nguồn thứ n vào trong nút lưới (i,j,k) - hi,j,k mực nước của nút (i,j,k) - pi,j,k,n , qi,j,k,n là các hệ số có thứ nguyên (L2t-1) và (L3t-1) tương ứng của phương trình. 47 Một cách tổng quát, nếu có N nguồn cấp vào trong ô lưới, lưu lượng tổng hợp Qsi,j,k có thể được viết như sau: Qsi,j,k = Pi,j,k hi,j,k + Qi,j,k (5) Trong đó: - Pi,j,k =Σ pi,j,k,n - Qi,j,k =Σ qi,j,k,n Viết cân bằng cho ô lưới (i,j,k) từ bước thời gian tm-1 đến tm ta có: Cri,j-1/2,k(hmi,j-1,k-hmi,j,k)+Cri,j+1/2,k(hmi,j+1,k-hmi,j,k)+ +Cci-1/2,j,k(hmi-1,j,k-hmi,j,k)+Cci+1/2,j,k(hmi+1,j,k-hmi,j,k)+ +Cvi,j,k-1/2(hmi,j,k-1-hmi,j,k)+Cvi,j,k+1/2(hmi,j,k+1-hmi,j,k)+ + Pi,j,khmi,j,k-1+Qi,j,k=Ssi,j,k(∆rj∆cj∆vk)( hmi,j,k-hm-1i,j,k)/(tm -tm-1) (6) Trong đó: - hmi,j,k là cốt cao mực nước tại bước thời gian m của ô (i,j,k) - Cri,j-1/2,k: là sức cản thấm trong hàng thứ i, lớp thứ k giữa các nút lưới (i,j-1,k) và (i,j,k). - Kri,j-1/2,k là hệ số thấm giữa các nút lưới (i,j,k) và (i,j-1,k) - ∆ci∆vk là diện tích bề mặt vuông góc với phương dòng chảy - ∆rj-1/2 là khoảng cách giữa các nút lưới (i,j,k) và (i,j-1,k) Phương trình trên sẽ được viết cho các ô mà mực nước thay đổi theo thời gian. Như vậy, ta sẽ lập được một hệ phương trình có số phương trình tương ứng với số ô lưới. Giải hệ phương trình này với điều kiện biết được mực nước hm-1i,j,k (điều kiện ban đầu) ta sẽ xác định được mực nước hmi,j,k. Cứ lần lượt như vậy, ta có thể xác định được mực nước cho bất kỳ thời điểm nào Hệ phương trình trên được giải bằng phương pháp lặp, người ta tiến hành chia nhỏ khoảng thời gian (tm-1,tm) kết quả nhận được là lời giải gần đúng của hệ phương trình. 48 Khi thời gian tăng lên thì h sẽ thay đổi. Khi h đạt được sự ổn định (chênh lệch h tính được giữa 2 bước thời gian kế cận nhau là nhỏ hơn một giá trị cho phép) thì mực nước đạt được sự cân bằng động và tại đây kết thúc quá trình tính toán. Để phương pháp lặp hội tụ, người ta chọn bước thời gian tăng theo cấp số nhân, khi đó thừa số 1/(tm-1 - tm) sẽ tiến nhanh tới 0 do đó các tổng có liên quan đến thừa số này hội tụ. * Điều kiện biên trong mô hình Có 3 loại điều kiện biên chính như sau: - Điều kiện biên loại I: là điều kiện biên áp lực trên đó được xác định trước (còn gọi là điều kiện biên Dirichlet). -Điều kiện biên loại II: là điều kiện biên lưu lượng được xác định trước (còn gọi là điều kiện biên Neumann). Trường hợp không có dòng chảy thì lưu lượng được xác định bằng không. - Điều kiện biên loại III: là điều kiện lưu lượng trên biên phụ thuộc vào sự thay đổi của áp lực (còn gọi là điều kiện biên Cauchy hoặc biên hỗn hợp). * Biên sông (River) Biên loại này được mô phỏng cho dòng chảy giữa tầng chứa nước và nguồn chứa nước thường là sông hay hồ... Nó cho phép dòng chảy từ tầng chứa vào trong nguồn chứa. Nước cũng có thể chảy từ nguồn chứa vào trong tầng chứa nước nhưng nguồn thấm này không phụ thuộc vào lưu lượng của sông. Hệ số sức cản thấm của biên sông được thể hiện trong công thức: Criv = KrLW/M (7) Trong đó : - Criv: giá trị sức cản thấm - Kr: hệ số thấm theo phương thẳng đứng của lớp trầm tích đáy lòng - L: chiều dài lòng sông trong ô - W: chiều rộng lòng sông trong ô 49 - M: chiều dày của lớp trầm tích đáy lòng Lưu lượng dòng thấm giữa sông và tầng chứa được tính theo công thức: Qriv = Criv (Hriv - h) khi h>Rbot (8) Trong đó: - Hriv: mực nước trong sông - h: mực nước của tầng chứa ngay dưới đáy lòng sông - Rbot: cốt cao đáy sông Trong trường hợp mực nước của tầng chứa nằm dưới đáy sông thì lúc đó lưu lượng dòng thấm sẽ đạt ổn định và tính theo công thức: Qriv = Criv (Hriv - Rbot) khi h ≤ Rbot (9) * Biên kênh thoát (Drain) Cơ chế hoạt động của biên loại này cũng không khác mấy so với biên sông ngoại trừ không cho phép nguồn thấm từ kênh vào tầng chứa. Điều này cũng có nghĩa rằng dòng thoát ra kênh Qd sẽ bằng không khi mực nước trong ô nhỏ hơn hoặc bằng cốt cao đáy kênh: Qd = 0 khi h ≤ d (10a) Khi mực nước nằm cao hơn đáy kênh thì lưu lượng dòng thoát ra kênh Qd sẽ được tính theo công thức: Qd = Cd (h-d) khi h > d (10b) Đối với kênh thoát, giá trị sức cản thấm Cd được tính như đối với sức cản thấm của biên sông Cr. * Biên mạch lộ (Spring) Loại biên này có thể mô phỏng bằng biên kênh thoát và chỉ hoạt động khi mực nước trong tầng chứa nằm cao hơn mặt đất. Sức cản thấm được đánh giá qua lưu lượng và mực nước của mạch lộ, mực nước của tầng chứa. * Biên bốc hơi (Evapotranspiration - ET) 50 Biên loại này đòi hỏi phải gán giá trị mô đun bốc hơi lớn nhất RETM cho các ô xảy ra quá trình bốc hơi. Giá trị này đạt được khi mực nước trong ô bằng với bề mặt địa hình (hs). Quá trình bốc hơi sẽ không xảy ra khi mực nước trong ô nằm dưới mực nước bốc hơi cho phép (d). Giữa hai giá trị này lượng bốc hơi (QET) sẽ được nội suy tuyến tính theo công thức: QET = QETM khi h>hs (11) Trong đó: - QETM = RETM .∆x.∆y QET = 0 khi h < (hs-d) (12) QET = QETM {h - (hs - d)}/d khi (hs-d) ≤ h ≤ hs (13) Hình 14: Điều kiện biên bốc hơi trong mô hình * Điều kiện biên tổng hợp (General head boundary - GHB) Điều kiện biên loại này cũng tương tự như điều kiện biên sông hoặc biên kênh thoát (hình 4). Lưu lượng dòng thấm qua biên được tính theo công thức: Qb = Cb (hb - h) (14) Sức cản thấm Cb cũng tương tự như sức cản thấm đáy lòng biểu thị sức cản dòng chảy giữa biên và tầng chứa nước. (hs - d) d h hs 51 Hình 15: Điều kiện biên tổng hợp (GHB) trong mô hình • Lỗ khoan hút nước hoặc ép nước (Well) Để mô phỏng các lỗ khoan hút nước trên mô hình, lưu lượng của các lỗ khoan trong ô lưới được đặt là lưu lượng tổng cộng QWT. QWT chính là bằng tổng lưu lượng của các lỗ khoan đặt trong các tầng chứa nước khác nhau (Σ Qi,j,k) (McDonald và Harbaugh,1988). Lưu lượng đơn lẻ cho các tầng chứa nước khác nhau đó được tính như sau: Qi,j,k = Ti,j,k (QWT/ΣTi,j,k) (15) Trong đó: - Ti,j,k là hệ số dẫn nước của tầng chứa nước - ΣTi,j,k là hệ số dẫn nước tổng cộng cho tất cả các lớp mà lỗ khoan khoan qua Tính hoàn chỉnh hay không hoàn chỉnh của lỗ khoan được mô phỏng bằng việc xác định vị trí đoạn ống lọc nằm trong tầng chứa nước mà lỗ khoan có trong thực tế Bán kính của lỗ khoan được mô phỏng trên mô hình lúc này sẽ là bán kính hiệu dụng re. Độ lớn của nó phụ thuộc vào kích thước của ô lưới và xác định theo công thức sau: Tầng chứa nước hi,j,k Nguồn cấp có mực nước không đổi Q hb,i,j,k Sức cản thấm (Cb,i,j,k) giữa nguồn và ô lưới i,j,k Ô i,j,k 52 re = 0.208a khi bước lưới đều a = ∆x = ∆y (16) c, Môdul vận chuyển nước dưới đất (MT3D) Môdul MT3D có thể mô phỏng cả 3 quá trình cơ bản của dịch chuyển chất đó là quá trình đối lưu, phân tán và khuyếch tán. Do đặc thù của chất gây mặn chủ yếu là do hàm lượng Nacl trong nước cao. Mà