BỘ GIÁO DỤC VÀ ĐÀO TẠO
BỘ NÔNG NGHIỆP VÀ PTNT
TRƢỜNG ĐẠI HỌC THỦY LỢI
ĐÀO TẤN QUY
NGHIÊN CỨU XÂY DỰNG MÔ HÌNH TOÁN
MÔ PHỎNG DÕNG CHẢY VÀ VẬN CHUYỂN BÙN CÁT TRÊN LƢU VỰC VỪA VÀ NHỎ
LUẬN ÁN TIẾN SĨ KỸ THUẬT
HÀ NỘI, NĂM 2017
BỘ GIÁO DỤC VÀ ĐÀO TẠO
BỘ NÔNG NGHIỆP VÀ PTNT
TRƢỜNG ĐẠI HỌC THỦY LỢI
ĐÀO TẤN QUY
NGHIÊN CỨU XÂY DỰNG MÔ HÌNH TOÁN
MÔ PHỎNG DÕNG CHẢY VÀ VẬN CHUYỂN BÙN CÁT TRÊN LƢU VỰC VỪA VÀ NHỎ
\
Chuyên ngành: Thủy văn học
Mã số: 62 44 90 01
NGƢỜI HƢỚNG DẪN KHOA HỌC: 1. PGS.TS. Phạm Thị Hƣơng Lan
2. PGS.TS. Ngô Lê Long
HÀ NỘI, NĂM 2017
LỜI CẢM ƠN
Để hoàn thành đƣợc luận án, tác giả bày tỏ lòng biết ơn sâu sắc tới tập thể thầy cô hƣớng dẫn: PGS.TS Phạm Thị Hƣơng Lan và PGS. TS Ngô Lê Long (Trƣờng Đại
học Thủy Lợi) về sự hƣớng dẫn tận tình trong suốt quá trình nghiên cứu và viết luận án.
Nhân dịp này, tác giả trân trọng cảm ơn Bộ môn Mô hình toán và Dự báo khí tƣợng thủy văn, Khoa Thủy văn và Tài nguyên nƣớc, Phòng Đào tạo Đại học và sau Đại học, Ban Giám hiệu Trƣờng Đại học Thủy Lợi đã tạo điều kiện thuận lợi nhất để
luận án đƣợc hoàn thành.
Tác giả cũng xin trân trọng cảm ơn Bộ môn Toán học – Khoa Công nghệ thông
tin - Trƣờng Đại học Thủy lợi đã tạo điều kiện tốt nhất, quan tâm giúp đỡ về mọi mặt trong suốt quá trình nghiên cứu và hoàn thành luận án.
Cuối cùng, tác giả xin chân thành cám ơn bạn bè, đồng nghiệp và gia đình đã
giúp đỡ, khích lệ tinh thần trong suốt quá trình thực hiện luận án.
Hà Nội, tháng 02 năm 2017
Tác giả
Đào Tấn Quy
i
LỜI CAM ĐOAN
Tác giả xin cam đoan đây là công trình nghiên cứu của bản thân tôi. Các kết quả
nghiên cứu và các kết luận trong luận án là trung thực, không sao chép từ bất kỳ một
nguồn nào và dƣới bất kỳ hình thức nào. Việc tham khảo các nguồn tài liệu đã đƣợc
thực hiện trích dẫn và ghi nguồn tài liệu tham khảo đúng quy định.
Hà Nội, tháng 02 năm 2017
Tác giả
Đào Tấn Quy
ii
MỤC LỤC
LỜI CẢM ƠN ................................................................................................................ i LỜI CAM ĐOAN .......................................................................................................... ii DANH MỤC BẢNG BIỂU .......................................................................................... vi DANH MỤC HÌNH VẼ .............................................................................................. vii DANH MỤC CÁC TỪ VIẾT TẮT............................................................................. ix MỞ ĐẦU…… ................................................................................................................ 1 1. TÍNH CẤP THIẾT CỦA LUẬN ÁN ..................................................................... 1 2. MỤC TIÊU NGHIÊN CỨU .................................................................................. 2 3. ĐỐI TƢỢNG VÀ PHẠM VI NGHIÊN CỨU ....................................................... 2 4. NỘI DUNG NGHIÊN CỨU ................................................................................... 2 5. PHƢƠNG PHÁP NGHIÊN CỨU .......................................................................... 3 6. Ý NGHĨA KHOA HỌC VÀ THỰC TIỄN CỦA LUẬN ÁN ................................ 3 7. NHỮNG ĐÓNG GÓP MỚI CỦA LUẬN ÁN ....................................................... 4 8. CẤU TRÖC CỦA LUẬN ÁN ................................................................................ 4 CHƢƠNG 1. TỔNG QUAN VỀ MÔ HÌNH TOÁN MÔ PHỎNG VẬN CHUYỂN BÙN CÁT TRÊN LƢU VỰC VỪA VÀ NHỎ .............................................................. 5 1.1. Tổng quan về xói mòn và vận chuyển bùn cát trên lƣu vực ............................ 5 1.1.1. Một số khái niệm, thuật ngữ ..................................................................... 5 1.1.1.1. Xói mòn lƣu vực ............................................................................... 5 1.1.1.2. Bùn cát và bồi lắng ............................................................................ 6 1.1.2. Nguyên nhân chính gây xói mòn, ảnh hưởng đến xói mòn ....................... 6 1.1.2.1. Các yếu tố gây xói mòn ..................................................................... 6 1.1.2.2. Các yếu tố ảnh hƣởng đến xói mòn ................................................... 8 1.1.3. Phân loại xói mòn lưu vực ...................................................................... 11 1.1.4. Vận chuyển bùn cát trên lưu vực ............................................................ 13 1.2. Các nghiên cứu trên thế giới .......................................................................... 13 1.2.1. Nghiên cứu đánh giá chung về xói mòn .................................................. 13 1.2.2. Nghiên cứu về các mô hình mô phỏng quá trình xói mòn và vận chuyển
bùn cát ..................................................................................................... 15 1.2.2.1. Mô hình kinh nghiệm ...................................................................... 15 1.2.2.2. Mô hình nhận thức .......................................................................... 25 1.2.3. Các thuật toán giải trong các mô hình mô phỏng vận chuyển bùn cát trên lƣu vực .................................................................................................... 30 1.3. Các nghiên cứu trong nƣớc ............................................................................ 31 1.3.1. Nghiên cứu đánh giá chung về xói mòn .................................................. 31
iii
1.3.2. Các nghiên cứu ứng dụng mô hình toán mô phỏng xói mòn và vận chuyển bùn cát ........................................................................................ 32
1.4. Những khoảng trống trong nghiên cứu xói mòn và vận chuyển bùn cát –
2.1.1.1. 2.1.1.2. 2.1.1.3.
Hƣớng nghiên cứu chính của luận án ............................................................ 35 1.4.1. Những khoảng trống trong nghiên cứu xói mòn và vận chuyển bùn cát 35 1.4.2. Định hướng nghiên cứu của luận án....................................................... 36 Kết luận chƣơng 1: ..................................................................................................... 36 CHƢƠNG 2. NGHIÊN CỨU CƠ SỞ KHOA HỌC XÂY DỰNG MÔ HÌNH TOÁN MÔ PHỎNG VẬN CHUYỂN BÙN CÁT TRÊN LƢU VỰC VỪA VÀ NHỎ .......... 37 2.1. Cơ sở lý thuyết – Nghiên cứu đề xuất thuật toán giải .................................... 37 2.1.1. Cơ sở lý thuyết ........................................................................................ 37 Phƣơng trình thấm ........................................................................... 37 Phƣơng trình mô phỏng dòng chảy ................................................. 38 Phƣơng trình diễn toán xói mòn và vận chuyển bùn cát ................. 40 2.1.2. Nghiên cứu đề xuất thuật toán giải ......................................................... 44 2.1.2.1. ... Thuật toán giải hệ phƣơng trình mô phỏng dòng chảy trên bề mă ̣t lƣu vực .................................................................................................................... 44 2.1.2.2. Thuật toán giải hệ phƣơng trình mô phỏng dòng chảy trong kênh/sông: ......................................................................................................... 49 2.1.2.3. Thuật toán giải phƣơng trình mô phỏng xói mòn và vận chuyển bùn cát trên lƣu vực ................................................................................................. 50 2.1.2.4. Thuật toán giải phƣơng trình mô phỏng xói mòn và vận chuyển bù n cát trong kênh/sông ........................................................................................... 55 2.2. Xây dựng các thành phần của mô hình .......................................................... 56 2.2.1. Quá trình liên rãnh ................................................................................. 57 2.2.1.1. Khả năng xói mòn liên rãnh ............................................................ 58 2.2.1.2. Vận chuyển bùn cát liên rãnh .......................................................... 59 2.2.2. Quá trình xói mòn rãnh ........................................................................... 61 2.2.2.1. Xói mòn rãnh ................................................................................... 61 2.2.2.2. Vận chuyển bùn cát trong rãnh ....................................................... 63 2.2.3. Quá trình lòng kênh/sông ........................................................................ 64 2.2.3.1. Sức tải bùn cát trong kênh/sông ...................................................... 64 2.2.3.2. Vận chuyển bùn cát trong kênh/sông .............................................. 65 2.3. Phân tích lựa chọn ngôn ngữ xây dựng mô hình ........................................... 65 2.4. Cấu trúc và chức năng của một số chƣơng trình con ..................................... 66 2.4.1. Cấu trúc và một số mô đun chính ........................................................... 66 Cấu trúc chƣơng trình: Có hai thành phần cơ bản ........................... 66 2.4.1.1 2.4.1.2 Một số mô đun chính ....................................................................... 66 2.4.2. Chức năng của một số chương trình con ................................................ 72
iv
Phân tích tƣơng quan giữa xói mòn rãnh trên lƣu vực với độ dốc và
2.5. Giao diện sử dụng chƣơng trình ..................................................................... 72 Kết luận chƣơng 2 ....................................................................................................... 74 CHƢƠNG 3. THỬ NGHIỆM MÔ HÌNH ĐỂ MÔ PHỎNG DÕNG CHẢY VÀVẬN CHUYỂN BÙN CÁT CHO MỘT SỐ LƢU VỰC VỪA VÀ NHỎ ............................ 75 3.1. Số liệu đầu vào cho mô hình .......................................................................... 75 3.1.1. Tạo cơ sở dữ liệu ..................................................................................... 75 3.1.2. Chạy mô hình .......................................................................................... 75 3.2. Thử nghiệm mô hình để mô phỏng xói mòn và vận chuyển bùn cát ............ 75 3.2.1. Lưu vực Nậm Sập .................................................................................... 75 3.2.1.1. Vị trí địa lý ...................................................................................... 75 3.2.1.2. Đặc điểm địa hình ........................................................................... 76 3.2.1.3. Đặc điểm thổ nhưỡng ...................................................................... 77 3.2.1.4. Đặc điểm thảm phủ thực vật và hiện trạng sử dụng đất ................. 79 3.2.1.5. Đặc điểm khí hậu............................................................................. 80 3.2.1.6. Đặc điểm khí tượng thủy văn .......................................................... 81 3.2.1.7. Yêu cầu số liệu đầu vào mô hình .................................................... 85 3.2.1.8. Hiệu chỉnh, kiểm định mô hình ....................................................... 86 3.2.2. Lưu vực Phiêng Hiềng ............................................................................ 90 3.2.2.1. Vị trí địa lý ...................................................................................... 90 3.2.2.2. Đặc điểm địa hình ........................................................................... 91 3.2.2.3. Đặc điểm thổ nhƣỡng ...................................................................... 92 3.2.2.4. Đặc điểm thảm phủ thực vật và hiện trạng sử dụng đất .................. 94 3.2.2.5. Đặc điểm khí hậu............................................................................. 95 3.2.2.6. Đặc điểm khí tƣợng thủy văn .......................................................... 96 3.2.2.7. Yêu cầu số liệu đầu vào mô hình .................................................. 100 3.2.2.8. Hiệu chỉnh, kiểm định mô hình ..................................................... 102 3.2.3. Xây dựng phƣơng trình tƣơng quan ...................................................... 105 3.2.3.1. Phân tích độ nhạy của các thông số mô hình ................................ 105 Phân tích tƣơng quan giữa xói mòn liên rãnh trên lƣu vực với độ dốc 3.2.3.2. và cƣờng độ mƣa ............................................................................................ 108 3.2.3.3. cƣờng độ mƣa ................................................................................................. 110 Kết luận chƣơng 3 ..................................................................................................... 112 KẾT LUẬN VÀ KIẾN NGHỊ .................................................................................. 113 DANH MỤC CÁC CÔNG TRÌNH CỦA TÁC GIẢ ĐÃ CÔNG BỐ ................... 115 TÀI LIỆU THAM KHẢO ........................................................................................ 116 PHỤ LỤC .................................................................................................................. 123
v
DANH MỤC BẢNG BIỂU
Bảng 1-1. Bảng 3-1. Bảng 3-2. Bảng tra để tính C theo hội Khoa học đất Quốc tế ................................ 23 Loại hình sử dụng đất tại lƣu vực Nậm Sập ........................................... 79 Loại hình sử dụng đất tại lƣu vực Phiêng Hiềng.................................... 94
vi
DANH MỤC HÌNH VẼ
Các nhân tố chính ảnh hƣởng đến xói mòn đất ........................................ 8 Hình 1-1. Tỷ lệ xói mòn tại các châu lục trên thế giới ........................................... 14 Hình 1-2. Sử dụng lƣới hình chữ nhật giải phƣơng trình trong mô hình động lực 31 Hình 1-3. Sơ đồ sai phân Lax- Friedrichs .............................................................. 44 Hình 2-1. Sơ đồ hình thành dòng chảy và vận chuyển bùn cát trên lƣu vực .......... 56 Hình 2-2. Sơ đồ tính toán xói mòn và vận chuyển bùn cát .................................... 57 Hình 2-3. Sơ đồ khối tính toán xói mòn và vận chuyển bùn cát liên rãnh ............. 66 Hình 2-4. Sơ đồ khối tính toán xói mòn và vận chuyển bùn cát rãnh .................... 69 Hình 2-5. Hình 2-6. Sơ đồ khối tính lƣợng bùn cát vận chuyển trong sông đến cửa ra ......... 69 Hình 2-7. Màn hình khởi động chƣơng trình ......................................................... 72 Giao diện mở tệp dữ liệu ........................................................................ 73 Hình 2-8. Giao diện chạy chƣơng trình .................................................................. 73 Hình 2-9. Bản đồ vị trí lƣu vực Nậm Sập ............................................................... 76 Hình 3-1. Bản đồ địa hình lƣu vực Nậm Sập ......................................................... 77 Hình 3-2. Bản đồ thổ nhƣỡng lƣu vực Nậm Sập .................................................... 78 Hình 3-3. Bản đồ hiện trạng sử dụng đất lƣu vực Nậm Sập................................... 80 Hình 3-4. Bản đồ mạng lƣới sông ngòi lƣu vực Nậm Sập ..................................... 82 Hình 3-5. Bản đồ lƣới trạm quan trắc trên lƣu vực Nậm Sập ................................ 83 Hình 3-6. Bản đồ đẳng trị mô đun dòng chảy lƣu vực Nậm Sập ........................... 84 Hình 3-7. Các tiểu lƣu vực của Nậm Sập ............................................................... 86 Hình 3-8. Hình 3-9. Đƣờng quá trình lƣu lƣợng tính toán và thực đo .................................... 87 Hình 3-10. Đƣờng quá trình hàm lƣợng bùn cát tính toán và thực đo ..................... 88 Hình 3-11. Đƣờng quá trình lƣu lƣợng tính toán và thực đo .................................... 88 Hình 3-12. Đƣờng quá trình hàm lƣợng bùn cát tính toán và thực đo ..................... 89 Hình 3-13. Hàm lƣợng bùn cát trên lƣu vực Thác Mộc ........................................... 89 Hình 3-14. Bản đồ vị trí lƣu vực Phiêng Hiềng ....................................................... 90 Hình 3-15. Bản đồ địa hình lƣu vực Phiêng Hiềng .................................................. 91 Hình 3-16. Bản đồ thổ nhƣỡng lƣu vực Phiêng Hiềng ............................................. 93 Hình 3-17. Bản đồ hiện trạng sử dụng đất lƣu vực Phiêng Hiềng năm 2015........... 95 Hình 3-18. Bản đồ mạng lƣới sông ngòi lƣu vực Phiêng Hiềng .............................. 97 Hình 3-19. Bản đồ lƣới trạm quan trắc trên lƣu vực Phiêng Hiềng ......................... 98 Hình 3-20. Bản đồ đẳng trị mô đun dòng chảy lƣu vực Phiêng Hiềng .................... 99 Hình 3-21. Các tiểu lƣu vực của Phiêng Hiềng ...................................................... 101
vii
Hình 3-22. Đƣờng quá trình lƣu lƣợng tính toán và thực đo .................................. 103 Hình 3-23. Đƣờng quá trình hàm lƣợng bùn cát tính toán và thực đo ................... 103 Hình 3-24. Đƣờng quá trình lƣu lƣợng tính toán và thực đo .................................. 104 Hình 3-25. Đƣờng quá trình hàm lƣợng bùn cát tính toán và thực đo ................... 104 Hình 3-26. Đƣờng biểu diễn kết quả khi thay đổi thông số mô hình ..................... 107 Hình 3-27. Tƣơng quan giữa xói mòn liên rãnh với độ dốc và cƣờng độ mƣa ...... 109 Hình 3-28. Tƣơng quan giữa xói mòn rãnh với độ dốc và cƣờng độ mƣa ............. 111
viii
DANH MỤC CÁC TỪ VIẾT TẮT
FAO: Food and Ariculture Organization
GIS: Geographic information system
LFW: Lax Friedrich Weight
MUSLE: Modified Universal Soil Loss Equation
USLE: Universal Soil Loss Equation
UNEP: United Nation Enviromemt Program
ix
MỞ ĐẦU
1. TÍNH CẤP THIẾT CỦA LUẬN ÁN
Trong những năm gần đây, dƣới tác động của các nhân tố tự nhiên và hoạt
động khai thác tài nguyên của con ngƣời, đất đai đang ngày càng bị thoái hóa
nghiêm trọng và biến đổi một cách nhanh chóng. Sự tác động tiêu cực của tự nhiên,
sự biến đổi của khí hậu toàn cầu gây nên tình trạng xói mòn, vận chuyển bùn cát,
thoái hóa đất trên lƣu vực, đặc biệt là các vùng đất dốc. Theo báo cáo của Bộ Khoa
học và Công nghệ nƣớc ta có hơn 13 triệu ha đất bị suy thoái thành đất trống, đồi
núi trọc, trong đó có 1,2 triệu ha đã bị xói mòn trơ sỏi đá.
Việt Nam có điều kiện tự nhiên phong phú, thuận lợi cho việc phát triển sản
xuất nông nghiệp đa dạng, cùng với đó các hoạt động canh tác trên nền đất dốc dẫn
đến nguy cơ xói mòn và bạc màu đất rất cao, đặc biệt vào mùa mƣa. Hầu hết các
sông suối đều ngắn, độ dốc lớn nên khả năng giữ nƣớc kém, tốc độ dòng chảy lớn,
nên khi có mƣa lớn thƣờng gây ra các hiện tƣợng lũ quét, xói mòn và rửa trôi mạnh
làm tăng nguy cơ mất đất toàn khu vực. Rõ ràng xói mòn và vận chuyển bùn cát
đang là một vấn đề toàn cầu hiện nay và đang có xu hƣớng gia tăng. Trong khi đó
quỹ đất canh tác hết sức hữu hạn và dân số không ngừng phát triển. Theo các
chuyên gia của FAO – UNEP hàng năm trên thế giới có khoảng từ 5 đến 7 triệu ha
đất bị mất khả năng sản xuất do bị xói mòn đất. Ở Việt Nam, với 3/4 diện tích là đồi
núi và nằm trong vùng nhiệt đới ẩm gió mùa, nên xói mòn đƣợc xem là một hiểm
họa đối với đất dốc ở Việt Nam. Nếu không có biện pháp phòng chống xói mòn thì
hàng trăm tấn đất và dinh dƣỡng sẽ bị mất sau mỗi năm và đất trở nên thoái hóa
không còn khả năng canh tác.
Một trong những phƣơng pháp tốt nhất hiện nay để tính toán lƣợng đất mất đi
do quá trình xói mòn và vận chuyển bùn cát trên lƣu vực là mô hình toán. Ƣu điểm
của phƣơng pháp này: Ít chi phí; có thể mô phỏng đƣợc trên các lƣu vực với quy mô
1
khác nhau và có thể dự tính lƣợng bùn cát bị xói mòn theo các kịch bản khác nhau.
Các nhà khoa học đã mô phỏng xói mòn và mô phỏng quá trình vận chuyển bùn cát
trên lƣu vực dƣới dạng các phƣơng trình toán học và giải chúng trên máy tính điện
tử. Đây là một trong những hƣớng nghiên cứu tiếp cận với những công cụ, phần
mềm máy tính hiện đại. Đã có rất nhiều nghiên cứu trong nƣớc cũng nhƣ trên thế
giới về xói mòn với cách tiếp cận mô phỏng xói mòn theo các mô hình toán, tuy
nhiên chủ yếu theo hƣớng sử dụng phƣơng trình mất đất phổ dụng (USLE). Tuy
nhiên theo hƣớng nghiên cứu này thì phạm vi ứng dụng mang tính địa phƣơng, có
độ chính xác hạn chế khi áp dụng ở những lƣu vực khác nhau, chƣa đề cập đến quá
trình bồi lắng, chuyển tải hạt cát và không có khả năng tính toán cho từng trận mƣa
hay các bƣớc thời gian ngắn hơn. Chính vì vậy việc “Nghiên cứu xây dựng mô hình
toán mô phỏng dòng chảy và vận chuyển bùn cát trên lƣu vực vừa và nhỏ” là cần
thiết và cấp bách, áp dụng tính toán cho các lƣu vực sông ở Việt Nam.
2. MỤC TIÊU NGHIÊN CỨU
- Nghiên cứu cơ sở khoa học xây dựng mô hình toán mô phỏng xói mòn và
vận chuyển bùn cát trên lƣu vực vừa và nhỏ.
- Thử nghiệm mô hình đã xây dựng cho một số lƣu vực vừa và nhỏ.
3. ĐỐI TƢỢNG VÀ PHẠM VI NGHIÊN CỨU
- Đối tƣợng nghiên cứu: Mô hình toán mô phỏng xói mòn và vận chuyển bùn
cát.
- Phạm vi nghiên cứu: Các lƣu vực vừa và nhỏ.
4. NỘI DUNG NGHIÊN CỨU
- Nghiên cứu tổng quan về các mô hình toán mô phỏng xói mòn và vận
chuyển bùn cát trên lƣu vực vừa và nhỏ trên thế giới và ở Việt Nam. Đánh giá tồn
tại về kỹ thuật và chỉ ra vấn đề mà luận án tập trung giải quyết.
- Ứng dụng cơ sở lý thuyết về cơ chế xói mòn và vận chuyển bùn cát để phát
2
triển mô hình mô phỏng quá trình xói mòn và vận chuyển bùn cát trên lƣu vực vừa
và nhỏ.
5. PHƢƠNG PHÁP NGHIÊN CỨU
- Phƣơng pháp phân tích, tổng hợp tài liệu: Tổng hợp các tài liệu liên quan đến
nghiên cứu mô hình toán thủy văn, các nghiên cứu về xói mòn và vận chuyển bùn
cát, thu thập, bổ sung, cập nhật các số liệu khí tƣợng thủy văn, địa hình, điều kiện tự
nhiên, các hoạt động quản lý và khai thác trên lƣu vực sông.
- Phƣơng pháp kế thừa: Sử dụng có chọn lọc các kết quả nghiên cứu, điều tra
cơ bản trƣớc đây có liên quan đến nội dung của luận án.
- Phƣơng pháp viễn thám và GIS: Sử dụng các phần mềm GIS, ArcGIS,
Mapinfo để biên tập, trình bày các bản đồ và tính toán các thông số [1].
- Phƣơng pháp mô hình toán mô phỏng xói mòn và vận chuyển bùn cát: Là
phƣơng pháp đƣợc sử dụng để lƣợng hóa quá trình dòng chảy và vận chuyển bùn
cát trên lƣu vực.
- Phƣơng pháp chuyên gia và tham vấn ý kiến cộng đồng: Tiếp cận, học hỏi
các nhà khoa học có nhiều kinh nghiệm trong công tác nghiên cứu, các chuyên gia
và nhà quản lý có kinh nghiệm thực tế liên quan đến lĩnh vực nghiên cứu của luận
án.
6. Ý NGHĨA KHOA HỌC VÀ THỰC TIỄN CỦA LUẬN ÁN
Ý nghĩa khoa học
- Đƣa ra đƣợc luận cứ, cơ sở khoa học về việc xây dựng và thử nghiệm mô
hình để mô phỏng quá trình xói mòn và vận chuyển bùn cát trên lƣu vực vừa và nhỏ
áp dụng trong điều kiện của Việt Nam.
- Kết quả nghiên cứu của luận án sẽ góp phần khẳng định việc nghiên cứu xây
dựng mô hình toán mô phỏng xói mòn và vận chuyển bùn cát với ứng dụng công
nghệ viễn thám, GIS [1] là rất hiệu quả và cần thiết trong giai đoạn hiện nay.
3
Ý nghĩa thực tiễn
- Kết quả nghiên cứu của luận án là nguồn tƣ liệu tham khảo quan trọng góp
phần cung cấp thông tin về tình hình dòng chảy, vận chuyển bùn cát trên lƣu vực
sông suối miền núi phía Bắc trên địa bàn của tỉnh Sơn La.
- Kết quả nghiên cứu của luận án sẽ cung cấp cho địa phƣơng nguồn dữ liệu và
công cụ có thể giám sát, đánh giá, tra cứu thông tin, theo dõi tác động của xói mòn và
vận chuyển bùn cát đến các hoạt động sản xuất khai thác sử dụng đất và nƣớc để từ
đó đề xuất các giải pháp quy hoạch tài nguyên nƣớc, quy hoạch sử dụng đất phù hợp.
7. NHỮNG ĐÓNG GÓP MỚI CỦA LUẬN ÁN
- Xây dựng đƣợc mô hình toán mô phỏng xói mòn và vận chuyển bùn cát trên lƣu
vực vừa và nhỏ, với thuật toán sơ đồ sai phân Lax – Friedrich có thêm trọng số thời
gian, không gian để giải phƣơng trình dòng chảy và phƣơng trình vận chuyển bùn
cát trên lƣu vực.
- Xây dựng đƣợc phƣơng trình tƣơng quan giữa tính toán xói mòn liên rãnh,
xói mòn rãnh trên lƣu vực nghiên cứu: Lƣu vực Nậm Sập, lƣu vực Phiêng Hiềng
của tỉnh Sơn La, từ đó có thể dự báo lƣợng bùn cát bị xói mòn và vận chuyển trên
lƣu vực theo cƣờng độ mƣa.
8. CẤU TRÖC CỦA LUẬN ÁN
Ngoài phần mở đầu và kết luận, luận án đƣợc trình bày theo 3 chƣơng:
Chƣơng I: Tổng quan về mô hình toán mô phỏng vận chuyển bùn cát trên lƣu
vực vừa và nhỏ.
Chƣơng II: Nghiên cứu cơ sở khoa học xây dựng mô hình toán mô phỏng vận
chuyển bùn cát trên lƣu vực vừa và nhỏ.
Chƣơng III: Thử nghiệm mô hình để mô phỏng vận chuyển bùn cát cho một số
lƣu vực vừa và nhỏ.
4
CHƢƠNG 1. TỔNG QUAN VỀ MÔ HÌNH TOÁN MÔ PHỎNG VẬN CHUYỂN BÙN CÁT TRÊN LƢU VỰC VỪA VÀ NHỎ
1.1. Tổng quan về xói mòn và vận chuyển bùn cát trên lƣu vực
1.1.1. Một số khái niệm, thuật ngữ
1.1.1.1. Xói mòn lƣu vực
Theo từ điển bách khoa toàn thƣ về khoa học đất, xói mòn xuất phát từ tiếng
Latin là "erodere” chỉ sự ăn mòn dần, thuật ngữ xói mòn để chỉ các quá trình liên
quan đến các lớp đất đá, đá tơi ra và bị mang đi bởi các tác nhân nhƣ gió, nƣớc,
băng, tuyết tan hoặc hoạt động của sinh vật. Có rất nhiều khái niệm về xói mòn, cụ
thể nhƣ sau:
Xói mòn là hiện tƣợng di chuyển đất bởi nƣớc mƣa, bởi gió dƣới tác động của
trọng lực lên bề mặt của đất. Xói mòn đất đƣợc xem nhƣ là một hàm số với biến số
là loại đất, độ dốc địa hình, mật độ che phủ của thảm thực vật, lƣợng mƣa và cƣờng
độ mƣa [2].
Xói mòn là hiện tƣợng các phần tử mảnh, cục và có khi cả lớp bề mặt đất bị
bào mòn, cuốn trôi do sức gió và sức nƣớc [3].
Xói mòn là quá trình san bằng, trong đó có các hạt đất hay đá cứng bị nhào
lộn, rửa trôi và di chuyển dƣới tác dụng của trọng lực, gió và nƣớc là động lực
chính của quá trình này [4].
Quá trình xói mòn, trƣợt lở, bồi lấp thực chất là quá trình phân bố lại vật chất
dƣới ảnh hƣởng của trọng lực, xảy ra khắp nơi và bị chi phối chủ yếu bởi điều kiện
địa hình [5].
Xói mòn do nƣớc phụ thuộc vào năng lƣợng xói mòn của dòng nƣớc và sức
kháng của đất đá đối với dòng nƣớc chảy qua.
5
1.1.1.2. Bùn cát và bồi lắng
Quá trình xói mòn và lắng đọng bao gồm các quá trình phá vỡ, vận chuyển,
lắng đọng của các hạt bùn cát do tác động của các tác nhân xâm thực và tác nhân
vận chuyển của dòng chảy mặt [6].
Bùn cát lơ lửng là những hạt bùn cát có kích thƣớc nhỏ, nổi lơ lửng khắp trong
dòng nƣớc và chuyển động trôi theo dòng nƣớc. Tốc độ chuyển động của loại bùn
cát này bằng tốc độ chuyển động của dòng nƣớc.
Mật độ hay khối lƣợng riêng bùn cát là khối lƣợng của một đơn vị thể tích bùn
cát. Kí hiệu là ρs và đơn vị thƣờng dùng là kg/m3.
Bồi lắng là các hạt đất tách ra do xói mòn đƣợc lắng lại trong đất hoặc bên
trong các nguồn nƣớc nhƣ: hồ, suối và đất ngập nƣớc [7].
1.1.2. Nguyên nhân chính gây xói mòn, ảnh hưởng đến xói mòn
1.1.2.1. Các yếu tố gây xói mòn
Nguyên nhân chủ yếu của xói mòn đất là tác động của mƣa, việc nghiên cứu
quá trình xói mòn đƣợc chia thành hai vấn đề: (1) Xói mòn biểu hiện ra sao ở các
dạng mƣa khác nhau và (2) Xói mòn xảy ra nhƣ thế nào trên các loại đất khác nhau.
Vì vậy, quy mô hoạt động xói mòn phụ thuộc vào hai nhân tố: Cƣờng độ mƣa và
khả năng kháng xói của đất. Nói một cách toán học, xói mòn là một hàm của tác
động xói mòn do mƣa và tính xói mòn của đất.
Theo kết quả nghiên cứu của nhiều nhà khoa học đã chỉ ra các nhóm nhân tố
chính gây ra xói mòn và vận chuyển bùn cát trên lƣu vực nhƣ sau:
a) Nhóm nhân tố mưa:
Nhân tố mƣa gây xói mòn chủ yếu thông qua lƣợng mƣa và cƣờng độ mƣa.
Lƣợng đất mất do xói mòn tỷ lệ thuận với lƣợng mƣa, cƣờng độ mƣa. Điều này
đồng nghĩa với một khu vực nghiên cứu nhất định thì những trận mƣa nhỏ có lƣợng
6
xói mòn ít hơn những trận mƣa lớn; cùng một lƣợng mƣa, trận mƣa nào có cƣờng
độ mƣa lớn hơn (thời gian mƣa ngắn hơn) thì lƣợng đất xói mòn sẽ nhiều hơn.
Những trận mƣa có cƣờng độ mƣa lớn sẽ làm cho lƣợng nƣớc mƣa không kịp ngấm
xuống đất, lƣợng nƣớc mƣa này gần nhƣ chuyển toàn bộ thành dòng chảy mặt làm
cho vận tốc dòng chảy mặt tăng đáng kể. Dòng chảy mặt càng lớn thì sẽ tạo nên
xung lực lớn tách và cuốn trôi những hạt đất mặt dẫn đến lƣợng đất xói mòn lớn [8].
Lƣợng mƣa, cƣờng độ mƣa, và sự phân bố sẽ quyết định đến lực phân tán các
hạt của đất, đến lƣợng nƣớc và vận tốc của nƣớc chảy tràn. Thời gian mƣa ngắn
cũng hạn chế xói mòn do không đủ lƣợng nƣớc hình thành dòng chảy. Khi cƣờng
độ mƣa lớn, thời gian mƣa kéo dài thì xói mòn rất nghiêm trọng.
b) Nhóm nhân tố thành phần cơ lý của đất
Tính xói mòn của đất chủ yếu thể hiện qua thành phần cơ giới, độ xốp của đất
và tình trạng bề mặt đất. Đối với đất có thành phần cơ giới nặng (đất sét) thì kích
thƣớc các hạt nhỏ, mịn, liên kết chặt, khó bị phá vỡ nên nguy cơ xảy ra xói mòn là
không cao. Đối với đất có thành phần cơ giới trung bình (đất thịt) thì kích thƣớc hạt
nhỏ vừa phải, liên kết vừa phải, tơi xốp, dễ bị cuốn trôi khi xuất hiện dòng chảy mặt
nên nguy cơ bị xói mòn cao. Đối với đất có thành phần cơ giới nhẹ (đất cát), mặc dù
có kết cấu kém bền vững nhƣng có kích thƣớc hạt lớn khó vận chuyển nên nguy cơ
xảy ra xói mòn là không cao, loại đất này có khả năng thấm nƣớc tốt nhƣng giữ
nƣớc kém [9].
Đất bị xói mòn thành các hạt cơ bản (cát, bùn và sét) và các kết hạt (tổ hợp
của các hạt cơ bản). Kích thƣớc của các kết hạt trong khoảng từ 2μm tới 500μm với
tỉ trọng tƣơng đối khoảng 1,8. Bùn cát đƣợc tạo ra bởi sự xói mòn ở vùng cao
thƣờng là hỗn hợp các loại hạt này. Nhƣng các tỉ lệ là các hàm của các đặc trƣng
đất, của tỉ lệ giữa bùn cát đến từ rãnh với bùn cát đến từ vùng liên rãnh và sự phân
loại từ sự lắng đọng phía trên sƣờn [10, 11].
7
1.1.2.2. Các yếu tố ảnh hƣởng đến xói mòn
Có 5 nhân tố chính ảnh hƣởng tới xói mòn đất là địa hình, loại đất, thảm thực
vật, khí hậu và con ngƣời.
Hình 1-1. Các nhân tố chính ảnh hưởng đến xói mòn đất
a) Địa hình
Địa hình ảnh hƣởng đến xói mòn đất chủ yếu thông qua độ dốc và chiều dài
sƣờn dốc. Độ dốc là nhân tố quan trọng ảnh hƣởng đến xói mòn và dòng chảy mặt.
Độ dốc càng lớn thì xói mòn mặt càng lớn và ngƣợc lại. Cùng một cấp độ dốc, nếu
chiều dài sƣờn dốc càng lớn thì nguy cơ gây xói mòn đất càng cao. Một số kết quả
nghiên cứu đã chỉ ra rằng: Nếu chiều dài sƣờn dốc tăng lên hai lần thì lƣợng đất xói
mòn cũng tăng xấp xỉ hai lần (đối với đất sản xuất lâm nghiệp) và tăng lên gần ba
lần trên đất trồng cà phê. Trong điều kiện nhiệt đới thì ảnh hƣởng của chiều dài
sƣờn dốc rõ nét hơn so với (điều kiện) các nƣớc ôn đới [9, 10].
b) Lớp phủ
Lớp phủ bao gồm: Tán lá thực vật, thảm mục, xác thực vật. Tán cây ngăn các
hạt mƣa và nếu nó che kín mặt đất, giọt nƣớc rơi xuống từ các lá cây có ít năng
lƣợng hơn rất nhiều so với các giọt mƣa không bị chắn [12, 13, 14]. Tuy nhiên,
nhiều tán cây có những chỗ hở để cho giọt mƣa rơi thẳng xuống bề mặt đất và phá
8
vỡ các hạt đất. Vật chất tiếp xúc bề mặt đất có hiệu lực làm giảm sự xói mòn so với
một tán che. Nơi mà bề mặt đất đƣợc che phủ không xảy ra sự phá vỡ bởi giọt mƣa
vì không có khoảng cách rơi cho các hạt nƣớc để tạo ra năng lƣợng.
Lớp phủ thực vật trên bề mặt lƣu vực cũng làm giảm đáng kể lƣợng xói mòn
đất trên bề mặt lƣu vực và do đó làm giảm đáng kể lƣợng bùn cát gia nhập sông,
ảnh hƣởng này thông qua hai tác dụng:
- Cản trở quá trình hình thành dòng chảy mặt từ mƣa, tăng thời gian tiếp xúc
giữa nƣớc và đất, tăng lƣợng nƣớc thấm xuống đất, tăng dòng chảy ngầm và giảm
lƣợng dòng chảy mặt [15].
- Bảo vệ bề mặt đất khỏi sự tiếp xúc trực tiếp với gió và nƣớc, tăng sức kháng
đối với xói mòn của đất đá.
Lớp (tầng) tán rừng: là tầng hoạt động thứ nhất trong hiệu ứng thuỷ văn rừng,
tán rừng có tác dụng giữ lại một phần nƣớc mƣa: theo Lee MacDonald [16] thì
lƣợng nƣớc này thƣờng biến động trong khoảng 5 - 30% tùy thuộc vào lƣợng mƣa
lớn hay nhỏ. Bởi vậy, trong cùng một trận mƣa thì lƣợng mƣa thực tế dƣới tán rừng
luôn thấp hơn lƣợng mƣa ngoài đất trống. Mặt khác tán rừng cũng cản phần lớn
lƣợng nƣớc mƣa không cho chúng tác động trực tiếp vào bề mặt đất rừng. Tuy
nhiên, động năng giọt nƣớc mƣa dƣới tán rừng (chiều cao tán rừng lớn hơn 10m)
lớn hơn động năng giọt nƣớc mƣa nơi đất không [16, 17]. Trong rừng, mƣa bị tán
cây giữ lại trƣớc khi rơi xuống đất, với những cây có diện tích phiến lá lớn sẽ tích tụ
đƣợc những hạt nƣớc lớn hơn hạt mƣa tự nhiên. Tầng tán rừng thƣờng ở độ cao trên
10m đủ để những hạt nƣớc đạt tới vận tốc cuối. Điều đó có nghĩa là năng lƣợng của
toàn bộ những hạt nƣớc rơi xuống từ tán lá cây lớn hơn năng lƣợng của những hạt
mƣa tự nhiên rơi xuống đất không có che phủ. Do vậy nếu không có lớp che phủ
dƣới tán thì lƣợng xói mòn trong rừng sẽ rất lớn.
Hoạt động nông nghiệp đã làm tăng lƣợng đất xói mòn lên nhiều lần so với đất
có thảm thực vật tự nhiên che phủ. Hàng năm toàn thế giới ƣớc tính tổng lƣợng phù
9
sa từ các con sông đổ ra biển đã tăng từ 9 tỷ tấn lên đến 24 tỷ tấn [18].
Sông Nƣớc Lƣợng phù sa hàng năm
Hoàng Hà Trung Quốc 1600 (tỷ tấn)
Ganges Ấn Độ 1455 (tỷ tấn)
Amazon Một số nƣớc 363 (tỷ tấn)
Missisipi Mỹ 300 (tỷ tấn)
Kosi Ấn Độ 172 (tỷ tấn)
Mekong Một số nƣớc 170 (tỷ tấn)
Nile Một số nƣớc 111 (tỷ tấn)
Nhƣ vậy yếu tố ảnh hƣởng đến xói mòn bao gồm: Địa hình (L.S), S là độ dốc,
L là chiều dài sƣờn dốc; tác động của thực vật và các tác động của con ngƣời. Xói
mòn đất là hàm số tác động xói mòn của mƣa và tính xói mòn của đất (Xói mòn đất
= (tính gây xói mòn của mƣa) x (tính xói mòn của đất)). Các nhân tố ảnh hƣởng
đƣợc thể hiện trong hình vẽ sau:
Các hoạt động của con ngƣời nhƣ trồng rừng, phá rừng, canh tác đất, xây dựng
các công trình thủy lợi trên sông đều làm thay đổi các điều kiện mặt đệm lƣu vực và
10
điều kiện dòng chảy trong sông nên ảnh hƣởng đến xói mòn lƣu vực, lòng sông và
làm thay đổi hàm lƣợng bùn cát trong sông.
1.1.3. Phân loại xói mòn lưu vực
Căn cứ vào tác nhân gây xói mòn, ngƣời ta phân xói mòn đất thành 5 loại: Xói
mòn do nƣớc, do gió, do trọng lực, do tuyết tan và dòng bùn đá. Trong phạm vi
nghiên cứu của luận án tập trung nghiên cứu xói mòn do nƣớc. Xói mòn do nƣớc
còn đƣợc phân thành [10, 19]:
- Xói mòn bắn tóe (splash erosion): (1) Các thành phần của đất bị tách rời ra,
do lực va đập của hạt mƣa (cấu trúc đất bị phá vỡ, các hạt bị phân tán thành các hạt
riêng rẽ); (2) các hạt sau khi bị phân tán đƣợc vận chuyển đến nơi thấp hơn do dòng
chảy tràn; và (3) các hạt đƣợc bồi lắng, tích tụ ở các vị trí mới. Các hạt mƣa rơi trực
tiếp lên mặt đất với một lực đủ mạnh để phá vỡ cấu trúc đất, tách rời các hạt. Trên
đất dốc, phần lớn các hạt này có khuynh hƣớng dịch chuyển dần về phía cuối dốc.
Trên phạm vi nhỏ, các vật liệu này thƣờng tích tụ nơi trũng thấp xung quanh.
- Xói mòn bề mặt (sheet erosion): Nƣớc mƣa chảy tràn kéo theo các hạt vật
liệu theo đƣờng di chuyển của dòng chảy.
- Xói mòn rãnh nhỏ (rill erosion): Khi xói mòn lớp xảy ra kéo dài, nƣớc có
khuynh hƣớng tập trung vào các đƣờng, rãnh nhỏ (rill) và hình thành kiểu xói mòn
rãnh. Kiểu xói mòn này dễ nhận thấy khi nƣớc chảy theo các rãnh nhỏ trên mặt đất.
Khi mặt đất đồng đều, thƣờng gây ra xói mòn theo từng lớp mỏng, nhƣng thực tế,
trên mặt đất, dù trên một diện tích nhỏ, mặt đất ít khi đồng đều, luôn luôn có nơi
cao, nơi thấp.
- Xói mòn rãnh (gully erosion): Hình thành khi nhiều rãnh nhỏ phát triển thành
các rãnh lớn (gully).
Theo tổng kết các dạng xói mòn đất, nguyên nhân, các nhân tố môi trƣờng nhƣ
sau:
11
Xói mòn Nhân tố môi trƣờng Các dạng biểu hiện Nguyên nhân: các nguồn năng lƣợng khác nhau
Năng lƣợng gió
Xói mòn do (Wind gió erosion) Vết lằn, ụ đám đất, khí bụi.
1. Vận tốc gió, sự nhiễu động không khí. 2. Hƣớng gió thịnh hành. 3. Khả năng chống chịu môi trƣờng phụ thuộc vào độ nhám của đất, lớp phủ cây trồng. 4. Khả năng chống chịu đất phụ thuộc vào cấu trúc đất, thành phần cơ giới, và chất hữu cơ.
Lở đất
Lực trọng trƣờng, áp lực từ việc canh tác trên đất. 1. Phụ thuộc vào cƣờng độ canh tác (tần suất và loại hình canh tác).
Xói mòn khô (Dry mechanical erosion)
2. Phụ thuộc vào độ dốc, mức độ kết dính đất.
Tác động của giọt mƣa
Mảng cát, lớp bồi lắng.
Xói mòn bề mặt (Sheet erosion)
1. Lớp phủ cây trồng. 2. Độ dốc. 3. Đất. 4. Công nghệ và cấu trúc kiểm soát xói mòn.
mòn
Dòng suối rãnh nhỏ, nƣớc.
Xói dạng tuyến (Linear erosion)
Năng lƣợng dòng thuộc chảy phụ vào thể tích dòng bình chảy và tốc phƣơng vận dòng chảy.
1. Vận tốc dòng chảy phụ thuộc vào độ dốc, độ nhám bề mặt. 2. Thể tích dòng chảy phụ thuộc vào kích thƣớc lƣu vực và khả năng thấm. 3. Phẫu diện và tầng đá mẹ của đất.
chuyển
Trƣợt lở đất, lũ bùn.
Lực trọng trƣờng, mất cân bằng độ dốc.
Di khối (Mass movement)
1. Khối lƣợng lớp thổ nhƣỡng + nƣớc + cây trồng. 2. Độ ẩm của khối đất trƣợt lở. 3. Địa hình: đá, mức độ không thấm, độ dốc.
12
1.1.4. Vận chuyển bùn cát trên lƣu vực
Quá trình xói mòn và vận chuyển bùn cát bao gồm các quá trình phá vỡ, vận
chuyển, lắng đọng của các hạt bùn cát do tác động của các tác nhân xâm thực và tác
nhân vận chuyển của mƣa rơi và dòng chảy mặt [20]. Vận chuyển là sự cuốn trôi và
di chuyển của các hạt bùn cát từ các vùng đất cao qua sông ngòi và cuối cùng ra đại
dƣơng [21]. Quá trình vận chuyển bùn cát trên lƣu vực là một quá trình phức tạp
phụ thuộc vào lƣợng dòng chảy, quá trình xói mòn, chuyển tải và bồi lắng trên
lƣu vực. Các quá trình đó đƣợc các tác giả mô hình hoá thông qua các mô hình
toán [22, 25].
1.2. Các nghiên cứu trên thế giới
1.2.1. Nghiên cứu đánh giá chung về xói mòn
Quá trình xói mòn hiện nay đƣợc gắn với các hoạt động nông nghiệp. Nhiều
ngƣời đã cho rằng đất đai bị khai thác cạn kiệt có thể là nguyên nhân khiến các nền
văn minh quá khứ mất đi. Vì vậy, cùng với thoái hóa đất, xói mòn tồn tại nhƣ một
vấn đề trong suốt quá trình phát triển của nhân loại.
Xói mòn xảy ra trên phạm vi toàn cầu nhƣng với các mức độ khác nhau. Vấn
đề xói mòn đất đặc biệt nghiêm trọng tại các quốc gia nhiệt đới và cận nhiệt đới vì
áp lực dân số lớn, sự khan hiếm đất nông nghiệp màu mỡ và nguồn lao động nông
nghiệp nghèo chiếm đa số. Hiểm họa xói mòn đất đã là mối quan tâm của con ngƣời
kể từ bình minh của nền nông nghiệp. Tuy nhiên, quy mô và mức độ trầm trọng chỉ
gia tăng trong thế kỉ 20 theo sự bùng nổ dân số và tình trạng thiếu quản lý đất canh
tác [26]. Tỉ lệ xói mòn thay đổi ở các khu vực khác nhau từ 30 - 40 tấn/ha/năm [27].
Ở phạm vi toàn cầu, ƣớc lƣợng khoảng 1960 triệu ha đất có xu hƣớng bị xói
13
mòn, chiếm khoảng 15% tổng diện tích đất thế giới, trong số đó, có 50% bị xói mòn
trầm trọng, và đa số bị bỏ hoang [28], tỉ lệ xói mòn thay đổi từ 0.5 đến 350 tấn/ ha/
năm. Hàng năm, toàn cầu có khoảng 75 - 109 triệu tấn đất bị mất đi, tƣơng đƣơng
với thiệt hại 400 tỉ USD/ năm.
Hình 1-2. Tỷ lệ xói mòn tại các châu lục trên thế giới
Các nghiên cứu đầu tiên về xói mòn đất đƣợc các nhà khoa học ngƣời Đức
thực hiện vào những năm 1877. Năm 1907 tại Mỹ các chƣơng trình nghiên cứu về
xói mòn đất đƣợc bắt đầu khi Bộ Nông nghiệp nƣớc này tuyên bố chính sách về bảo
vệ nguồn tài nguyên đất [29]. Các nghiên cứu chi tiết đầu tiên về mƣa tác động đến
xói mòn [30] đƣợc tiến hành nghiên cứu và phân tích các tác động cơ học của hạt
mƣa lên đất và đƣa ra tiến trình xói mòn [2]. Công thức toán học đƣợc Zingg đƣa ra
để đánh giá ảnh hƣởng của độ dốc và độ dài của sƣờn dốc đến sự xói mòn [31].
Hiện nay xói mòn và vận chuyển bùn cát trên lƣu vực đƣợc nghiên cứu dƣới
nhiều loại hình và tính chất khác nhau. Xu hƣớng phổ biến trong nghiên cứu xói
mòn trên thế giới là mô hình hóa, diễn tả động lực của quá trình xói mòn và nghiên
cứu xói mòn kết hợp với các khoa học khác, chủ yếu để tìm hiểu quá trình cũng nhƣ
tác động của xói mòn lên môi trƣờng nhằm có đƣợc các biện pháp chống xói mòn
khả thi. Nghiên cứu theo hƣớng mô hình hóa theo hai trƣờng phái: Trƣờng phái
14
theo Bennett khẳng định rằng việc hình thành các rãnh trên bề mặt là do sự vận
chuyển vật chất rắn, do năng lƣợng dòng chảy [26]. Vì vậy, kiểm soát xói mòn nên
tập trung vào các biện pháp cơ giới để làm giảm tốc độ dòng chảy và năng lƣợng
xói mòn của nó, mà không cần giảm lƣợng dòng chảy trên lƣu vực. Trƣờng phái
dựa trên nghiên cứu của Ellison khẳng định rằng dòng chảy gia tăng theo sự suy
thoái cấu trúc bề mặt từ tác động của giọt mƣa [2]; do đó, kiểm soát xói mòn cần
tiến hành trên bề mặt lƣu vực, tập trung vào lớp phủ cây trồng, công nghệ trồng trọt.
1.2.2. Nghiên cứu về các mô hình mô phỏng quá trình xói mòn và vận chuyển
bùn cát
Mô hình xói mòn và vận chuyển bùn cát đƣợc phân thành hai loại [32]: Mô
hình kinh nghiệm; Mô hình nhận thức. Cụ thể các nghiên cứu xây dựng các mô hình
biểu diễn xói mòn, vận chuyển bùn cát trên lƣu vực nhƣ sau:
1.2.2.1. Mô hình kinh nghiệm
Các nghiên cứu về xây dựng mô hình kinh nghiệm tính toán xói mòn lƣu vực
phải kể đến các nghiên cứu nhƣ [31, 33, 36].
Mô hình kinh nghiệm đầu tiên đƣợc Zingg xây dựng để tính lƣợng bùn cát bị
xói mòn ở sƣờn dốc [31], trong đó lƣợng bùn cát bị xói mòn phụ thuộc vào độ dốc
và chiều dài sƣờn dốc . Sƣ̉ du ̣ng các dƣ̃ liê ̣u tƣ̀ các nhà nghiên cƣ́ u khác và các dƣ̃ liê ̣u tƣ̀ kết quả thí nghiê ̣m củ a mình , Zingg đã đề xuất mối quan hê ̣ giữa xói mòn với chiều dài, độ dốc của sƣờn dốc và đề xuất công thức tính xói mòn trên lƣu vực
nhƣ sau:
A = CS1.37L0.60 (1- 1)
trong đó:
A : Sự xói mòn đất trung bình trên mỗi đơn vị diện tích
L, S : Chiều dài, độ dốc của sƣờn dốc
C : Hằng số.
15
Smith đã đánh giá các ảnh hƣởng của các giải pháp bảo toàn cơ học cho 4 sƣ̣ kết hơ ̣p [35] gồm 1) sƣ̣ luân canh mù a vu ̣; 2) sƣ̣ xƣ̉ lý , canh tác đất; 3) lớ p phủ thƣ̣c
vâ ̣t (cây trồng, mùa vụ) (C); 4) các hệ số áp dụng hỗ trợ (P) cho phƣơng trình Zingg
để tính xói mòn trên lƣu vực.
A = CS7/5L3/5P (1- 2)
Ellison thông qua thí nghiê ̣m đã chỉ ra ảnh hƣở ng củ a năng lƣơ ̣ng đô ̣ng ho ̣c của mƣa đối với sự bong tách các phân tử đất [33], trong đó lƣơ ̣ng đất bi ̣ bong tách đƣợc tính toán trong quá trình thờ i đoa ̣n kéo dài 30 phút, mô tả trong phƣơng trình sau:
E = KV4.33D1.07I0.65 (1- 3)
trong đó:
E : Lƣợng đất bị bong tách trong quá trình bị bắn tóe , trong mô ̣t thờ i đoạn kéo
dài 30 phút (grams)
K : Mô ̣t hằng số
V : Vận tốc rơi củ a các ha ̣t mƣa (fps)
D : Đƣờng kính của các giọt mƣa (mm)
I : Cƣờng độ mƣa (inch/giờ ).
Musgrave đã đƣa ra phƣơng trình tính xói mòn lƣu vực [34], trong đó có xét
đến yếu tố khí hậu là cƣờng độ mƣa thời đoạn 30 phút tính toán cho lƣợng đất xói
mòn trong một mẫu đất Anh (0.4 ha) và lƣơ ̣ng nƣớ c mƣa lớ n nhất trong 30 phút, tần
1.75
suất 2 năm, (inch).
E = IRS1.35L0.35P30 (1- 4)
trong đó:
E : Lƣợng đất xói mòn trong mô ̣t mẫu đất Anh (0.4 ha)
16
I : Sự xói mòn riêng của đất (in)
R : Hệ số lớp phủ
S : Độ dốc mái (%)
L : Chiều dài mái dốc (ft)
P30 : Lƣợng nƣớc mƣa lớn nhất trong 30 phút, tần suất 2 năm, (inch).
Phƣơng trình Musgrave đƣơ ̣c sƣ̉ du ̣ng rô ̣ng rãi để ƣớ c tính lƣơ ̣ng xói mòn đất
tổng cô ̣ng tƣ̀ các lƣu vƣ̣c đầu nguồn.
Wischmeier và Smith đã xây dựng phƣơng trình mất đất phổ dụng USLE
(Universal Soil Loss Equation) để tính lƣợng bùn cát bị xói mòn [36], trong đó có
xét thêm các nhân tố ảnh hƣởng đến xói mòn nhƣ hệ số ảnh hƣởng của cây trồng và
ảnh hƣởng của các biện pháp canh tác đến xói mòn đất.
Xói mòn đƣợc xem nhƣ tích số của hệ số xói mòn của mƣa (nhân tố mƣa R),
hệ số xói mòn của đất (nhân tố đất K), hệ số chiều dài sƣờn dốc và độ dốc (nhân tố
địa hình LS), hệ số cây trồng và kĩ thuật canh tác (nhân tố thảm thực vật C), và hệ
số kiểm soát xói mòn (nhân tố thực hành bảo tồn P). Vì là một tích số nên nếu một
nhân tố có xu hƣớng bằng 0, xói mòn sẽ có xu hƣớng bằng 0. Phƣơng trình cơ bản
(mô hình kinh nghiệm) tính toán xói mòn trên bề mặt lƣu vực có dạng nhƣ sau:
A = R.K.LS.C.P (1- 5)
trong đó:
A là lƣợng đất mất trung bình hàng năm (tấn/ ha). Trong phƣơng trình USLE,
xói mòn đƣợc định nghĩa là tổng lƣợng đất đƣợc chuyển tới chân sƣờn dốc nơi các
quá trình lắng đọng bắt đầu diễn ra, hoặc các dòng chảy bắt đầu tập trung lại.
Hệ số R bằng E (động năng của lƣợng mƣa) nhân với I30 (cƣờng độ mƣa tối đa
trong 30 phút, đơn vị tính: cm/ giờ). Hệ số này tƣơng ứng với nguy cơ xói mòn tiềm
năng trong một khu vực nhất định mà xói mòn bề mặt (sheet) xuất hiện trên một
mảnh đất trống có độ dốc 9%.
17
Hệ số K phụ thuộc vào chất hữu cơ, kết cấu của đất, tính thấm của nó và cấu
trúc phẫu diện. Nó thay đổi từ 7/10 cho đất mỏng nhất đến 1/100 cho đất ổn định
nhất. Nó đƣợc đo trên các mảnh đất trống dài 22,2m và độ dốc 9%, canh tác theo
hƣớng dốc và không nhận đƣợc chất hữu cơ trong ba năm.
Hệ số LS phụ thuộc vào chiều dài sƣờn dốc và độ dốc.
Hệ số C là tỉ lệ đơn giản giữa xói mòn trên đất trống và xói mòn quan sát trên
một hệ thống cây trồng. Nhân tố này là sự kết hợp của lớp phủ cây trồng, trình độ
sản xuất và kỹ thuật thu hoạch liên quan. Nó thay đổi từ 1 trên đất trống đến 1/1000
dƣới rừng, 1/100 trên đồng cỏ và cây trồng, và 1 đến 9/10 dƣới gốc và củ cây trồng.
Hệ số P phản ánh thực hành kiểm soát xói mòn cụ thể nhƣ cày xới theo đƣờng
đồng mức, hoặc bậc thang.
Hệ số R của các trận mƣa đƣợc tính theo công thức của Brown và Foster [37]:
(1- 6)
trong đó:
R : Hệ số xói mòn do mƣa
j : Khoảng thời gian mƣa trong mỗi trận mƣa (phút)
n : Số khoảng thời gian trong trận mƣa
j : Tổng lƣợng mƣa trong khoảng thời gian j (mm)
Ij : Cƣờng độ của mƣa trong khoảng thời gian j (mm/h)
I30 : Cƣờng độ cực đại trong khoảng thời gian 30 phút mƣa (mm/h)
Vj : Vận tốc dòng chảy trong khoảng thời gian j (m/s).
Hệ số xói mòn do mƣa đƣợc xác định theo động năng mƣa và cƣờng độ mƣa
cực đại trong khảng thời gian 30 phút của trận mƣa.
Hệ số R này áp dụng trong những điều kiện có đủ thiết bị do mƣa ở độ phân
giải cao (các bƣớc thời gian ngắn, ví dụ: 1 phút, 2 phút, 5 phút, 10 phút….). Trong
18
trƣờng hợp những vùng nghiên cứu không có điều kiện thu thập các dữ liệu mƣa thì
có thể sử dụng công thức tính R theo lƣợng mƣa của tháng hoặc năm.
Trong những trƣờng hợp số liệu mƣa theo trận không sẵn có mà chỉ có lƣợng
mƣa theo năm thì có thể tính R từ lƣợng mƣa năm theo công thức thực nghiệm của
Renard và Freimund [38] nhƣ sau:
với F < 55mm
(1- 7) với F ≥ 55mm
Công thức Larionov (Liên Xô cũ):
(1- 8)
trong đó:
Xi : Lƣợng mƣa của trận mƣa thứ i
Ii,30 : Cƣờng độ mƣa lớn nhất trong 30 phút của trận mƣa thứ i.
Một số công thức tính hệ số R khác đƣợc thống kê theo Nguyễn Kim Lợi,
2005 trong bảng sau:
Tác giả Công thức
Roose (1975)
Chỉ số xói mòn tính theo lƣợng mƣa hàng năm (P) R=0,5 x P x 1,73
Morgan (1974)
Chỉ số xói mòn tính theo lƣợng mƣa hàng năm (P) R=9,28 x P-8,838
Foster et al (1981)
Chỉ số xói mòn tính theo lƣợng mƣa hàng năm (P) và I30 R=0,276 x P x I30
El-Swaify and others 1985
19
Tác giả Công thức
Chỉ số xói mòn tính theo lƣợng mƣa hàng năm (P) R=38,5 + 0,35 (P)
Wanapiryarat et al (1986)
Chỉ số xói mòn tính theo lƣợng mƣa hàng ngày (x) R=-3,2353+1,789ln(x)
Công thức của Nguyễn Trọng Hà (Đại học Thủy Lợi – Hà Nội)
Chỉ số xói mòn tính theo lƣợng mƣa hàng năm (P) R=0,548257P – 59,9
Việc xác lập công thức để tính toán cho hệ số R phụ thuộc vào từng khu vực nhất định do mỗi vùng đều có sự khác nhau về lƣợng mƣa, sự phân bố, tính chất mƣa v.v... Cƣờng độ mƣa càng lớn và thời gian mƣa càng lâu, tiềm năng xói mòn càng cao. Giá trị R thay đổi từ năm này qua năm khác nên việc xác định hệ số R chung là rất khó, muốn tính đƣợc hệ số R một cách chính xác phải dựa vào chế độ mƣa và số liệu thống kê của vùng nghiên cứu cụ thể qua nhiều năm. Khi tính toán hệ số R cho các khu vực khác nhau thì ta có thể áp dụng các công tính R của các khu vực đã nghiên cứu, nhƣng ta phải chọn công thức tính hệ số R phù hợp nhất với khu vực nghiên cứu.
K là hệ số thể hiện khả năng xói mòn của đất (soil erodibility). Nói cách khác đây là một nhân tố biểu thị tính dễ bị tổn thƣơng của đất với xói mòn và là đại lƣợng nghịch đảo với tính kháng xói mòn của đất. Đất có giá trị K càng lớn thì khả năng xói mòn càng cao. K phụ thuộc vào đặc tính của đất chủ yếu là sự ổn định về cấu trúc đất, thành phần cơ giới đất. Đặc biệt là ở các tầng đất trên mặt, thành phần cơ giới, hàm lƣợng hữu cơ có trong đất. Một số công thức tính hệ số K nhƣ sau:
Tác giả Công thức
100K=2,1.10-4M1,14(12-a) + 3,25(b-2) + 2,5(c-3) Wischmeier và Smith (1978)
Rosewell (1993) 100K=2,27M1,14(10-7)(12-a) + 4,28(10-3)(b-2) + 3,29(10-3)(c- 3)
ISSS (1995) 100K=2,241[2,1.10-4(12-a)M1,14 + 3,25(b-2) + 2,5(c-3)]
20
K : Hệ số xói mòn đất
M : Đƣợc xác định: (%) M = (% limon + % cát min)(100% - % sét)
a : Hàm lƣợng chất hữu cơ trong đất đo bằng phần trăm
b : Hệ số phụ thuộc vào hình dạng, sắp xếp và loại kết cấu đất
c : Hệ số phụ thuộc tiêu thấm của đất.
Hệ số LS là đại lƣợng biểu thị cho sự ảnh hƣởng của nhân tố độ dốc (S) và độ
dài sƣờn dốc (L) tới hoạt động xói mòn đất. Về mặt lý thuyết, khi tăng tốc độ dòng
chảy lên gấp đôi thì mức độ vận chuyển đối với các hạt có thể lớn hơn 64 lần, nó
cho phép mang các vật liệu hòa tan trong nƣớc lớn hơn gấp 30 lần và kết quả làm
tăng sức mạnh xói mòn gấp 4 lần [39]. S là độ dốc của sƣờn dốc , lƣợng mất đất lớn
khi độ dốc cao. Nó là tỷ lệ của sự mất đất từ độ dốc thực tế đối với độ dốc chuẩn
(9%) dƣới những điều kiện khác đồng nhất, sự liên hệ của sự mất đất đối với độ dốc
bị ảnh hƣởng bởi mật độ lớp phủ thực vật và kích thƣớc hạt đất. L là khoảng cách từ
đƣờng phân thủy ở đỉnh dốc đến nơi vận tốc dòng chảy chậm lại và vật chất bị trầm
lắng. Nó là tỷ số lƣợng mất đất ở các loại đất giống nhau có độ dốc giống nhau
nhƣng có chiều dài sƣờn khác nhau, so với chiều dài sƣờn của ô đất chuẩn (72.6
feet). --> L và S là 2 yếu tố đƣợc xem xét chung khi tính toán xói mòn. Tùy thuộc
vào từng khu vực mà ta có cách tính toán LS cho phù hợp.
Wischmeier và Smith đã đƣa ra công thức tính LS [36] nhƣ sau:
LS = (x/22,13)n (0,065 +0,045* s + 0,0065*s2) (1- 9)
trong đó:
x : Chiều dài sƣờn thực tế tính bằng đơn vị m
s : Độ dốc (%)
n : Thông số thực nghiệm: n = 0.5 khi S > 5%; n = 0.4 khi 3.5% < S < 4.5%
n = 0.3 khi 1% < S < 3.5%; n = 0.2 khi S < 1%.
C : Hệ số tỉ số giữa lƣợng đất mất trên một đơn vị diện tích có lớp phủ thực
vật và có sự quản lý của con ngƣời đối với lƣợng đất mất trên một diện tích trống
21
tƣơng đƣơng. Hệ số C là hệ số đặc trƣng cho mức độ che phủ đất của các lớp thực
phủ bề mặt, biện pháp quản lý lớp phủ, biện pháp làm đất, sinh khối đất. Giá trị của
C chạy từ 0 đến 1. Đối với vùng đất trống không có lớp phủ thực vật thì hệ số C
đƣợc xem là bằng 1.
Hệ số này đặc trƣng cho nhân tố lớp phủ thực vật đối với quá trình xói mòn.
Nhìn chung nhân tố cây trồng bao gồm:
- Tác dụng của lớp che phủ vòm lá và chiều cao của chúng (CI)
- Tác dụng của lớp phủ giữ nƣớc đối với xói mòn do dòng chảy bề mặt của đất
cũng nhƣ tốc độ thấm của đất (CII)
- Tác dụng của bộ rễ đối với quá trình thấm cũng nhƣ khả năng xói của đất
(CIII)
Hệ số C bằng tích của ba hệ số CI, CII và CIII. Trong đó CI đƣợc xác định
theo công thức sau:
CI = 1 – Fc.e(-0,1.H)
CII = 5,348.Sp2 – 3,711.Sp + 1,107 đối với Sp ≤ 30%
CII = 5,348.Sp2 – 1,116.Sp + 0,750 đối với Sp > 30%
CIII = 6600.106.X2 + (-1,6596.103) + 0,45 đối với các loại đất đƣợc
phủ bằng cây cỏ dại.
CIII = X2/54000 + X/180 + 0,45 đối với các loại đất đƣợc phủ bằng
các loại cây cỏ khác.
trong đó:
Fc : Phần trăm diện tích đất đƣợc che phủ bởi vòm lá (%)
H : Khoảng cách từ vòm lá đến mặt đất
Sp : Phần trăm diện tích đất đƣợc che phủ lớp mùn giữ nƣớc (%)
X : Phần trăm diện tích đất đƣợc che phủ lớp rễ cây cỏ (%)
22
Ngoài ra còn có thể xác định C theo bảng tra của hội Khoa học đất Quốc tế
Bảng 1-1. Bảng tra để tính C theo hội Khoa học đất Quốc tế
Cây và bụi cây với chiều cao khác nhau (không phủ kín đất) Cây hàng năm Độ che phủ (%) Rừng nhiệt đới lớp phủ (>50mm)
Bãi chăn thả, cây lâu năm thấp và có lớp phủ 1,00 0 4m 1,00 2m 1,00 1m 1,00 0,5m 1,00 1,00
0,55 0,30 10 20 0,97 0,95 0,95 0,90 0,93 0,83 0,92 0,83 0,009 0,55 0,30
0,17 30 0,92 0,85 0,79 0,75 0,17
0,09 0,05 40 50 0,89 0,87 0,80 0,75 0,72 0,65 0,66 0,58 0,003 0,09 0,06
0,027 0,015 60 70 0,84 0,81 0,70 0,65 0,58 0,51 0,50 0,41 0,001 0,056 0,053
0,008 80 0,78 0,60 0,44 0,33 0,050
0,005 0,002 90 100 0,76 0,73 0,55 0,50 0,37 0,30 0,24 0,16 0,0001 0,047 0,043
Trong phƣơng trình USLE thì yếu tố P đánh giá hiệu quả của các phƣơng thức
canh tác, phản ánh các hoạt động làm đất của con ngƣời nhằm bảo vệ đất trong việc
hạn chế xói mòn trên vùng đất dốc. Trong công thức mất đất phổ dụng, giá trị của
hệ số P đƣợc thành lập từ 3 yếu tố phụ và đƣợc tính theo công thức:
P = Pc * Pst * Pter
trong đó: Pc: Yếu tố phụ làm đất theo đƣờng đồng mức. Pst: Yếu tố phụ đƣờng viền
thức vật theo đƣờng đồng mức. Pter: Yếu tố phụ đắp bờ ngăn xói mòn.
Mô hình kinh nghiệm nói chung là đơn giản vì chủ yếu dựa trên phân tích các
quan sát và phân tích quan hệ đặc điểm từ các dữ liệu đo đạc [26]. Các giá trị tham
số trong mô hình kinh nghiệm có thể xác định thông qua việc hiệu chỉnh mô hình,
nhƣng thƣờng đƣợc xác định từ việc kiểm định các số liệu quan trắc thực tế. Các dữ
23
liệu này đặc biệt hữu ích nhƣ là một bƣớc quan trọng trong việc xác định nguồn gốc
bùn cát trên lƣu vực. Việc tính toán xói mòn kênh hay sự xói mòn mƣơng; không
chính xác cho trận mƣa đơn. Phƣơng trình tính xói mòn không đề cập đến vấn đề
bồi lắng trên lƣu vực và không tính toán cho một trận mƣa cũng nhƣ không ƣớc
lƣợng đƣợc sự lắng đọng, không ƣớc lƣợng đƣợc sự xói mòn rãnh sâu hay sự xói
mòn kênh. Chỉ giới hạn trong phạm vi lƣu vực nhỏ, với lƣu vực lớn đòi hỏi phải có
nhiều số liệu, dẫn đến sai số lớn khi tính toán. Việc xác định các hệ số trong phƣơng
trình USLE đƣợc xây dựng dựa trên số liệu của Mỹ, do đó việc áp dụng để tính cho
các lƣu vực khác sẽ không tránh khỏi những sai số. Tuy nhiên USLE có hiệu quả để
ƣớc lƣợng lƣợng đất mất trung bình trên một thời gian dài và lƣợng đất mất trung
bình hàng năm. Ƣu điểm của mô hình đơn giản, cơ sở dữ liệu lớn, các giá trị tham
số có sẵn, đƣợc sử dụng rộng rãi bởi các cơ quan nhƣ USDA thích nghi với các
vùng không đồng nhất nơi không xảy ra sự lắng đọng.
Williams đã nghiên cứu cải biên mô hình USLE và đƣa ra phƣơng trình
MUSLE tính lƣợng bùn cát sinh ra từ một trận mƣa [26]:
msf = 0,95(V.Qs)0,5.K.LS.C.P (1- 10)
trong đó:
msf : Lƣợng bùn cát tại mặt cắt cửa ra (kg) V : Tổng lƣợng lũ (m3) Qs : Lƣu lƣợng đỉnh lũ (m3/s)
K, LS, C, P: Nhƣ phƣơng trình USLE.
Mặc dù MUSLE giả thiết rằng sự lắng đọng xảy ra trong lƣu vực, nó chỉ đƣa
ra một ƣớc lƣợng của tổng sản lƣợng bùn cát và không đƣa ra ƣớc lƣợng của sản
lƣợng từng loại hạt riêng biệt. Sự lắng đọng tách riêng các hạt. Sự lắng đọng các hạt
xảy ra dễ dàng hơn ngay sau khi chúng rời khỏi nguồn, trong khi các hạt càng nhỏ,
càng nhẹ, di chuyển càng xa ngang qua lƣu vực trƣớc khi lắng đọng. Sự giảm theo
quy luật hàm mũ có thể đƣợc sử dụng biểu diễn vận chuyển bùn cát qua lƣu vực để
ƣớc lƣợng sự phân phối này. Mô hình MUSLE đơn giản, dễ áp dụng, sử dụng tập
24
hợp số liệu tham số của USLE, tuy nhiên phƣơng trình MUSLE không cung cấp
thông tin về phân bố theo thời gian của lƣợng bùn cát trong thời gian xảy ra dòng
chảy.
Các hạn chế của mô hình kinh nghiệm nhƣ sau:
Mô hình này chỉ áp dụng cho xói mòn bề mặt (sheet), vì thế nó không áp dụng
cho xói mòn dạng tuyến (linear) hoặc xói mòn khối (mass).
Mô hình đã đƣợc kiểm nghiệm và xác nhận tại bán bình nguyên và đất nƣớc
nhiều đồi núi với độ dốc sƣờn 1-20%. Mô hình không áp dụng cho các sƣờn dốc có
độ dốc lớn hơn 40%, nơi dòng chảy lớn hơn là một nguồn năng lƣợng lớn hơn so
với mƣa, và nơi có chuyển động khối của trái đất mạnh mẽ.
Các mối quan hệ giữa động năng và cƣờng độ mƣa thƣờng đƣợc sử dụng trong
mô hình này không áp dụng cho các vùng miền núi.
Mô hình này chỉ áp dụng cho dữ liệu trung bình trên 20 năm và không phù
hợp khi dùng cho các trận mƣa đơn.
Cuối cùng, một hạn chế lớn của mô hình là nó bỏ qua sự tƣơng tác nhất định
giữa các nhân tố. Ví dụ nó không tính đến ảnh hƣởng xói mòn độ dốc kết hợp với
che phủ thực vật, cũng nhƣ tác động của loại đất với ảnh hƣởng của độ dốc.
1.2.2.2. Mô hình nhận thức
Khác với mô hình kinh nghiệm, mô hình nhận thức đƣợc phát triển dựa vào
hiểu biết về các qui luật vận động và cơ chế vật lý của quá trình xói mòn, nghĩa là
dựa vào các hiểu biết đã đƣợc lý thuyết hoá dƣới dạng các định luật hay phƣơng
trình vật lý. Các quá trình vật lý của xói mòn có thể đƣợc kể ra gồm: quá trình tách
hạt đất (do năng lƣợng của hạt mƣa rơi hoặc một dạng năng lƣợng khác); quá trình
chuyển tải (với các định luật về dòng chảy mà quá trình này tuân thủ) và quá trình
bồi lắng của các hạt đất.
Các mô hình nhận thức mô phỏng xói mòn và vận chuyển bùn cát đƣợc xây
dựng dựa vào các phƣơng trình toán học mô phỏng các hiện tƣợng vật lý của quá
25
trình xói mòn rửa trôi đất. Cơ sở toán học của các mô hình toán là phƣơng trình liên
tục của Bennett [26]. Phƣơng trình liên tục thƣờng đƣợc sử dụng trong mô hình
động lực học của sự xói mòn đất dốc là:
(1- 11)
trong đó:
qs : Tải trọng bùn cát (khối lƣợng/1 đơn vị chiều rộng x 1 đơn vị thời gian)
x : Khoảng cách dọc theo chiều dốc
ρs : Mật độ khối của các hạt bùn cát
c : Nồng độ bùn cát trong dòng chảy (khối lƣợng bùn cát/một đơn vị thể tích
dòng chảy)
y : Độ dày lớp dòng chảy
t : Thời gian
Dr : Tốc độ xói mòn rãnh (khối lƣợng bị xói mòn /một đơn vị diện tích. một
đơn vị thời gian)
Dl : Tốc độ phân phối của bùn cát từ các vùng rãnh (khối lƣợng/một đơn vị
diện tích x một đơn vị thời gian).
Một số mô hình toán mô phỏng xói mòn và vận chuyển bùn cát trên lƣu vực
có thể kể đến nhƣ sau:
Mô hình ANSWERS đƣợc Beasley và các cộng sự đƣa ra vào năm 1980 (mô phỏng xói mòn đất trên lƣu vực có xem xét sƣ̣ khuếch tán chất lƣợng nƣớc) [41] bao gồm: mô ̣t quá trình dòng ch ảy và mô ̣t quá trình xói mòn dƣ̣a trên những hiện tƣợng vâ ̣t lý xảy ra trên lƣu vực . Quá trình xói mòn với giả thiết bùn cát bị bong tách do quá trình mƣa rơi, dòng chảy và dòng chảy cuốn trôi. Mô hình ANSWERS chia lƣu vƣ̣c thành các tiểu lƣu vực (ô lƣớ i) nhỏ và độc lập . Trong mỗi phần tử ô lƣới, các quá trình xói mòn và dòng chảy đƣợc xử lý nhƣ hai hàm độc l ập của cá c thông số
xói mòn và dòng chảy của mỗi phần tử đó . Trong mô hình, điều kiê ̣n bề mă ̣t và đô ̣ sâu dòng chảy tràn củ a mỗi phần (ô lƣớ i ) đƣơ ̣c xem là bằng nhau . Các khe rãnh
26
không đƣơ ̣c xem xét. Ảnh hƣởng của các khe dẫn nƣớc đƣợ c mô tả thông qua hê ̣ số . Trong mô nhám của phƣơng trình Mainning đã đƣợc sử dụng trong mô hình này
hình ANSWERS không xem xét thành phần bùn cát trong dòng chảy ngầm và
không xét đến vấn đề xói mòn trong kênh. Dƣ̃ liê ̣u đầu vào của mô hình ANSWERS khá lớn.
Mô hình CREAMS (Mô hình về hóa học, dòng chảy và xói mòn từ hệ quản lý
nông nghiệp) [27]. Các thành phần vận chuyển trầm tích của CREAMS (các chất
hóa ho ̣c, dòng chảy và xói mòn tƣ̀ các hê ̣ thống quản lý (canh tác) nông nghiê ̣p) đã phân tích sƣ̣ nƣ́ t bề mă ̣t và khe rãnh mô ̣t cách riêng biê ̣t . Sƣ̣ bong tách trong các khu vƣ̣c nƣ́ t bề mă ̣t và khe rãnh đƣơ ̣c xác đi ̣nh bở i mô hình USLE đã đƣơ ̣c chỉnh sƣ̉ a , thay đổi trong suốt quá trình dòng cải tiến . Quy trình này cho phép các thông số chảy và dọc theo các đƣờng dẫn nƣớc để miêu tả chính xác sự biến thiên theo không
gian.
Mô hình AGNPS (Mô hình nguồn ô nhiễm phân tán từ nông nghiệp) đƣợc xây
dựng do các nhà khoa học nông nghiệp Mỹ [42] nhằm dự đoán xói mòn và sự di
chuyển các chất dinh dƣỡng, hóa chất từ những lƣu vực nông nghiệp. Mô hình sẽ
chạy cho từng trận mƣa. Lƣu vực đƣợc chia ra thành những tiểu lƣu vực và lƣới ô
vuông có kích thƣớc bằng nhau. Mô hình bao gồm 3 mô hình con, đó là: Mô hình
xói mòn dựa trên phƣơng trình tính toán xói mòn mất đất phổ dụng USLE. Mô hình
dòng chảy dựa trên đƣờng cong thấm (SCS) mà mỗi giá trị đại diện cho một chế độ
đất, cây trồng, hệ số dẫn nƣớc, hệ số dòng chảy khác nhau. Quá trình vận chuyển
các chất dinh dƣỡng đất và hóa chất đƣợc mô hình hóa dựa trên các đặc tính của
đất, lƣợng các chất bón vào đất và khả năng vận chuyển của dòng chảy. Mô hình
làm việc trên cơ sở hệ thống lƣới ô vuông. Các ô vuông là những ô có diện tích
đồng nhất bằng nhau. Sự phân chia này cho phép ta phân tích mọi diện tích trong
lƣu vực. Các thành phần chính trong mô hình là phần thủy lực, xói mòn đất, vận
chuyển bùn cát và vận chuyển các chất dinh dƣỡng cây trồng nhƣ đạm, lân và các
chất hóa học. Trong phần thủy lực mô hình tính toán chủ yếu là thể tích dòng chảy
27
mặt (dòng lũ, đỉnh lũ). Tổng lƣợng xói mòn trên đất dốc. Xói mòn rãnh tính toán
trên 5 nhóm cấp hạt là sét, thịt, hạt kết nhỏ, hạt kết lớn và cát. Vận chuyển bùn cát
cũng đƣợc tính cho cả 5 cấp hạt cho từng ô vuông. Phần vận chuyển chất ô nhiễm
đƣợc chia nhỏ thành một phần tính toán các chất ô nhiễm gắn liền với bùn cát. Mô
hình không những mô phỏng cho từng trận mƣa đơn mà còn mô phỏng theo năm có
tính đến các yếu tố cân bằng nƣớc bao gồm: mƣa, bốc thoát hơi nƣớc, cây hút, chảy
tràn và thấm sâu. Vì thế động thái biến đổi của độ ẩm đất là yếu tố quan trọng ảnh
hƣởng đến khả năng sinh dòng chảy cho mỗi trận mƣa đƣợc tính toán một cách
khoa học. Tuy nhiên mô hình chủ yếu tập trung đánh giá sự di chuyển các chất dinh
dƣỡng, hóa chất từ các hoạt động sản xuất nông nghiệp trên lƣu vực.
Mô hình KINEROS (Mô hình xói mòn đô ̣ng lƣ̣c) là tập hợp các phần tử của
mô ̣t ma ̣ng lƣớ i nhƣ thƣ̣c vâ ̣t ho ̣c , các con kênh tự nhiên hoặc các ống dẫn nƣớc , ao
hồ, các khu chứ a mà chú ng liên kết vớ i nhau [43, 44]. Mô hình KINEROS có sƣ̣
hơ ̣p nhất củ a các phần vâ ̣n chuyển bù n cát và xói mòn . Thành phần vận chuyển bùn
cát của mô hình là dựa trên phƣơng trình liên tục không ổn địn h mô ̣t chiều. Mƣ́ c đô ̣
xói mòn, bồi lắng là sự kết hợp giữa xói mòn do mƣa rơi và mƣ́ c đô ̣ xói mòn /bồi
lắng thủ y lƣ̣c. Mƣ́ c đô ̣ xói mòn do mƣa rơi đƣơ ̣c đƣa ra bở i mô ̣t phƣơng trình kinh
nghiê ̣m mà trong đó mƣ́ c đô ̣ là tƣơng ƣ́ ng vớ i năng lƣơ ̣ng th ứ 2 của mƣa. Mƣ́ c đô ̣
xói mòn thủy lực (dòng nƣớc) đƣơ ̣c ƣớ c tính là tƣơng ƣ́ ng vớ i đô ̣ thiếu hu ̣t khả năng
vâ ̣n chuyển mà nó là sƣ̣ khác nhau giƣ̃a hàm lƣơ ̣ng bù n cát hiê ̣n ta ̣i trong dòng chảy
và nồng độ lớn nhất ở trạng thái ổn định.
Mô hình WEPP (Water Erosion Prediction Project) là mô hình tính toán xói
mòn dựa trên quá trình vật lý mô phỏng sự xói mòn và vận chuyển bùn cát trên lƣu
vực [45]. Mô hình WEPP tính toán lƣợng đất bị xói mòn dọc theo sƣờn dốc và
lƣợng chuyển tải bùn cát tại điểm cuối của sƣờn dốc. Sự xói bề mặt và xói rãnh nhỏ
cũng đƣợc mô phỏng trong mô hình này. Xói mòn bề mặt đƣợc mô phỏng nhƣ là
quá trình xói mòn do tác động của hạt mƣa, cƣờng độ mƣa, chuyển tải bùn cát do
lớp dòng chảy mặt tới các rãnh nhỏ. Xói mòn rãnh là hàm số của lớp dòng chảy
28
mặt, dung tích chuyển tải bùn cát và lƣợng bùn cát hiện có trong dòng chảy. Quá
trình hình thành dòng chảy mặt đƣợc coi là quá trình hình thành dòng chảy trên bề
mặt và quá trình tập trung dòng chảy trong rãnh nhỏ trên bề mặt đó. Quá trình diễn
toán dòng chảy đƣợc thực hiện bằng cách giải phƣơng trình sóng động học. Trong
lƣu vực, xói mòn trong kênh dẫn đƣợc xác định khi ứng suất kéo lớn hơn giá trị ứng
suất tiếp đáy. Sự chuyển động của bùn cát dọc theo sƣờn dốc hoặc trên bề mặt lƣu
vực đƣợc miêu tả trong mô hình WEPP dựa trên cơ sở của phƣơng trình liên tục
bùn cát [46] do Foster và Meyer xây dựng mà nó áp dụng cho điều kiện dòng chảy
trong rãnh.
Mô hình WEPP đã sử dụng phƣơng trình để tính toán xói mòn và trầm tích:
Di = Ki*Ie 2 *Ge*Ce*Sf (1- 12)
trong đó:
Di : Lƣợng trầm tích chuyển từ xói mòn mảng sang khu vực xói mòn dòng
(kg/m2 /s).
K i : Tính xói mòn mảng của mảng (kg/m3/s)
: Tác động của cƣờng độ mƣa (m/s) Ie
Ge : Nhân tố điều chỉnh lớp phủ
Sf = 1,05 – 0,85 exp(-4Sin): Nhân tố điều chỉnh độ dốc
Ce : Hệ số (1/m) (kg/(m2.h)).
Mô hình xói mòn đất Châu Âu EUROSEM là một mô hình dự báo xói mòn
đất bở i nƣớ c tƣ̀ các lƣu vƣ̣c nhỏ [47, 48]. EUROSEM có thể mô phỏng vâ ̣n chuyển bùn cát, xói mòn và bồi lắng bở i các quá trình trong các khe rãnh và trong các khe . Mô hình mô phỏng tổng lƣơ ̣ng dòng chảy , tổng nƣ́ t trên bề mă ̣t ở các sƣờ n dốc lƣơ ̣ng đất mất đi , quan hệ mƣa – trầm tích . Trong mô hình EUROSEM , sƣ̣ bong tách hạt đấ t do tác động mƣa rơi là tổng của mƣa rơi trực tiếp xuống và phụ thuộc
vào động năng của mƣa . Xói do nƣớc mƣa bắn tóe diễn ra trƣớ c khi dòng chảy bắt 0. Có 3 trƣờ ng đầu. Vì thế , hàm lƣợng bùn cát ban đầu nên đƣợc lấy bằng giá trị
29
hơ ̣p đƣơ ̣c xem xét cho vấ n đề xói mòn ở sƣờ n đồi bao gồm : Sƣờ n dốc không có các khe rãnh; sƣờ n dốc có khe rãnh và cá c khe nƣ́ t bề mă ̣t ; Sƣờ n dốc có các khe rãnh
dày ch ằng chịt . Quá trình xói mòn trong kênh đƣơ ̣ c xƣ̉ lý tƣơng tự nhƣ xói mòn
trong khe rãnh vớ i ngoa ̣i lê ̣ là tác đô ̣ng mƣa rơi là không đáng quan tâm và dòng
chảy bên (dòng chảy nhập lƣu) đối vớ i kênh tƣ̀ sƣờ n dốc trở nên quan tro ̣ng.
Mô hình “Công cụ đánh giá đất và nƣớc” SWAT (Soil and Water Assement
Tools) là một mô hình vật lý đƣợc xây dựng từ những năm 90 do tiến sĩ Jeff Arnold
thuộc trung tâm nghiên cứu đất nông nghiệp USDA- Agricultural Research Service
(ARS) xây dựng. Mô hình này đƣợc xây dựng để mô phỏng ảnh hƣởng của việc
quản lý sử dụng đất đến nguồn nƣớc, bùn cát và hàm lƣợng chất hữu cơ trong đất
trên hệ thống lƣu vực sông trong một khoảng thời gian nào đó. Mô hình SWAT tính
toán bồi lắng vào xói mòn cho những đoạn kênh có kích thƣớc giống nhau cho việc
mô phỏng toàn bộ hệ thống kênh. SWAT đƣợc phát triển từ mô hình mô phỏng tài
nguyên nƣớc lƣu vực nông thôn SWRRP cùng với một số mô hình khác nhƣ: Hệ
thống quản lý nông nghiệp về hóa chất, rửa trôi và xói mòn CREAMS mô hình
những ảnh hƣởng của sự tích trữ nƣớc ngầm GLEAMS và mô hình tính toán ảnh
hƣởng của hoạt động sản xuất đến xói mòn EPIC. Trong mô hình SWAT, lƣợng
bùn cát đƣợc tính bằng phƣơng trình mất đất phổ dụng cải tiến (MUSLE) [40].
MUSLE là một phƣơng trình sửa đổi từ phƣơng trình mất đất phổ dụng [25, 36]
(USLE) của Wishchmeier và Smith.
1.2.3. Các thuật toán giải trong các mô hình mô phỏng vận chuyển bùn cát
trên lƣu vực
Các mô hình toán mô phỏng xói mòn và vận chuyển bùn cát trên lƣu vực nêu
trên đều dựa trên phƣơng trình liên tục của Bennet. Các mô hình giải bằng phƣơng
pháp số trị nhƣ mô hình CREAM, EUROSEM... còn các mô hình giải bằng phƣơng
pháp giải tích nhƣ mô hình SWAT, WEPP... Với trƣờng hợp giải bằng phƣơng pháp
giải tích cho sƣờn dốc có độ dốc đồng nhất và với lƣợng mƣa vƣợt thấm.
Với trƣờng hợp tổng quát bài toán thƣờng đƣợc giải bằng phƣơng pháp số trị,
30
sử dụng các sơ đồ sai phân hữu hạn và một lƣới thời gian và không gian đặc trƣng
nhƣ đƣợc chỉ ra trong hình sau:
Hình 1-3. Sử dụng lưới hình chữ nhật giải phương trình trong mô hình động lực
1.3. Các nghiên cứu trong nƣớc
Việt Nam có địa hình chủ yếu là đồi núi, xói mòn đất diễn ra thƣờng xuyên
nên hiện tƣợng xói mòn đã đƣợc nghiên cứu từ rất sớm và đã đƣa ra kết luận
phƣơng pháp canh tác ruộng bậc thang sẽ giảm hiện tƣợng xói mòn.
1.3.1. Nghiên cứu đánh giá chung về xói mòn
Việt Nam nằm trong vùng nhiệt đới ẩm và có lƣợng mƣa tƣơng đối lớn (từ
1800 – 2000mm) nhƣng phân bố không đồng đều và tập trung chủ yếu vào mùa
mƣa. Lƣợng mƣa lớn tập trung tạo ra dòng chảy có cƣờng độ rất lớn, đây là nguyên
nhân chính gây ra hiện tƣợng xói mòn đất ở Việt Nam. Phù sa đổ ra biển Đông khoảng 200 triệu tấn/năm, ngƣời ta ƣớc tính trung bình 1m3 nƣớc trong sông chứa từ 50g - 400g bùn cát, riêng đồng bằng sông Hồng 1000g/m3 và có lúc đạt 2000g/m3.
Với tổng diện tích đất tự nhiên 33.121 triệu ha, cùng với khoảng 25 triệu ha
đất dốc, chiếm hầu hết lãnh thổ miền núi và trung du. Với biến động của môi
trƣờng, Việt Nam đang đứng trƣớc nguy cơ thoái hóa đất do xói mòn rửa trôi là rất
lớn.
Các số liệu thống kê hiện trạng sử dụng đất năm 2008 cho thấy, Việt Nam có
khoảng 25 triệu ha đất dốc có nguy cơ xói mòn và rửa trôi rất lớn. Theo các quan
31
trắc có hệ thống từ năm 1960 đến nay thì có khoảng 10% - 20% lãnh thổ bị ảnh
hƣởng xói mòn từ trung bình đến mạnh, đặc biệt là khu vực miền núi và trung du.
Theo Nguyễn Quang Mỹ [49] thì lịch sử nghiên cứu xói mòn của nƣớc ta có
thể chia làm 3 giai đoạn: Trƣớc năm 1954: Giai đoạn này chỉ mới bắt đầu xuất hiện
các biện pháp canh tác chống xói mòn nhƣ làm ruộng bậc thang, xây kè cống.v.v...
chứ xói mòn đất chƣa đƣợc nghiên cứu đƣa lên thành lý luận. Từ 1954 - 1975: Giai
đoạn này bắt đầu xuất hiện một số công trình nghiên cứu và nhiều biện pháp canh
tác chống xói mòn đƣợc đƣa ra hàng loạt trên các nông trƣờng miền núi phía Bắc.
Một số nghiên cứu đáng chú ý giai đoạn này nhƣ: Chu Đình Hoàng nghiên cứu sự
ảnh hƣởng của giọt mƣa đến xói mòn đất và chống xói mòn bằng biện pháp canh tác
[39]. Nguyễn Tử Xiêm và Thái Phiên với công trình nghiên cứu về đất đồi núi Việt
Nam [6]. Về mặt lý luận các tác giả đã đánh giá đƣợc năng lực phòng hộ của một số
dạng cấu trúc thảm thực vật rừng trong việc chống xói mòn và tiến hành các nghiên
cứu với quy mô rộng hơn.
1.3.2. Các nghiên cứu ứng dụng mô hình toán mô phỏng xói mòn và vận
chuyển bùn cát
Từ 1975 đến nay các công trình nghiên cứu trong nƣớc đã áp dụng phƣơng
trình mất đất phổ dụng của Wischmeier và Smith nhƣ sau:
Phạm Ngọc Dũng đã tiến hành nghiên cứu về ứng dụng phƣơng trình mất đất
phổ dụng vào dự báo tiềm năng xói mòn đất và đƣa ra các biện pháp chống xói mòn
cho các tỉnh Tây nguyên [9] và xác định giá trị của các yếu tố gây xói mòn đất theo
mô hình Wischmeier W.H và Smith D.D. Các nghiên cứu trên cũng chỉ tập trung
nghiên cứu xói mòn trên lƣu vực trên cơ sở ứng dụng phƣơng trình mất đất phổ
dụng để tính toán xói mòn trên lƣu vực mà chƣa xem xét đến quá trình dòng chảy
và vận chuyển bùn cát trên lƣu vực.
Nguyễn Trọng Hà trong nghiên cứu các yếu tố gây xói mòn và khả năng dự
báo xói mòn trên đất dốc [8] đã đƣa ra chỉ số xói mòn của mƣa và hệ số xói mòn đất
K để tiến tới áp dụng phƣơng trình mất đất phổ dụng USLE của Wischmeier và
32
Smith, dự báo đƣợc xói mòn đất và sơ bộ thể hiện tiềm năng xói mòn do mƣa trên
bản đồ tỷ lệ nhỏ. Nghiên cứu trên chƣa xem xét đến các yếu tố khác gây xói mòn
nhƣ độ dốc, độ dài sƣờn dốc, hệ thống cây trồng đa dạng, hệ số bảo vệ đất. Kết quả
của nghiên cứu mới chỉ ở mức độ dự báo xói mòn khái quát.
Phạm Hùng [11] đã phát triển mô hình toán kết hợp với phƣơng trình mất đất
phổ dụng của Wischmeier và Smith dựa trên cơ sở dữ liệu đƣợc hình thành từ hệ
thống thông tin địa lý (GIS) để đánh giá nguy cơ xói mòn lƣu vực. Bƣớc đầu tính
toán cho thấy nhiều khả năng ứng dụng tốt để đánh giá nguy cơ xói mòn lƣu vực
trong điều kiện của nƣớc ta. Trên cơ sở phân tích mƣa ngày lớn nhất và chỉ số xói
mòn do mƣa, nghiên cứu đã tính toán thử nghiệm ngƣỡng mƣa gây lũ quét và xác
lập các cấp khả năng xuất hiện lũ quét trên các lƣu vực điển hình. Tuy nhiên với kết
quả nghiên cứu trên chƣa đánh giá đƣợc phân bố theo thời gian của lƣợng bùn cát bị
xói mòn do các trận mƣa và chƣa xét đến cấp độ dốc trên lƣu vực là một trong
những nguyên nhân gây xói mòn lƣu vực và ảnh hƣởng đến khả năng vận chuyển
bùn cát trên lƣu vực.
Lại Vĩnh Cẩm đã sử dụng phƣơng trình mất đất phổ dụng để đánh giá tiềm
năng và mức độ xói mòn hiện tại của 4 lƣu vực lớn ở Miền Bắc Việt Nam [50]. Kết
quả nghiên cứu về xói mòn tại các lƣu vực đƣợc tích hợp với phân tích lƣu vực
nhằm chỉ ra các khu vực xói mòn nguy hiểm làm cơ sở cho đề xuất các biện pháp
phòng tránh hữu hiệu nhằm mục tiêu phát triển bền vững.
Phạm Hữu Tỵ đã mô phỏng rủi ro xói mòn vùng cảnh quan đồi núi trên cơ sở
sử dụng số liệu viễn thám và mô hình mất đất phổ dụng hiệu chính (RUSLE) [13],
trên cơ sở sử dụng dữ liệu địa hình, khí hậu, thảm thực vật v.v... thu đƣợc từ vệ tinh.
Quá trình nghiên cứu đã kiểm chứng công thức tính toán hệ số xói mòn đất của mƣa
R, hệ số địa hình LS, hệ số che phủ thực vật và biện pháp canh tác C của RUSLE
thông qua việc sử dụng ảnh vệ tinh Landsat ETM+, mô hình quan trắc Radar STM
và mô hình phân tích pha trộn quang phổ tuyến tính LSM. Kết quả đã xây dựng
đƣợc bản đồ đánh giá rủi ro xói mòn vùng và mô hình hỗ trợ quyết định quy hoạch
33
sử dụng đất lâm nghiệp vùng cảnh quan đồi núi lƣu vực sông Hƣơng thuộc huyện
Hƣơng Trà tỉnh Thừa Thiên Huế.
Nguyễn Văn Dũng đã đánh giá định lƣợng xói mòn đất đồi núi vùng Thanh –
Nghệ - Tĩnh [10] bằng phƣơng trình mất đất phổ dụng và hệ thống thông tin địa lý.
Kết quả tính hệ số K cho đồi núi vùng Thanh - Nghệ - Tĩnh cho thấy: Có sự khác
biệt với kết quả công bố hệ số K của tác giả Nguyễn Trọng Hà. Việc xác định hệ số
C cũng gặp khó khăn tƣơng tự nhƣ việc xác định hệ số K. Do điều tra, kiểm kê hiện
trạng sử dụng đất, thảm thực vật rừng không có thông tin chi tiết của lớp thảm thực
vật nhƣ độ che phủ của các tầng tán, tầng thảm mục v.v... Điều này dẫn tới không
xác định đƣợc các thông số cho việc xác định hệ số C. Từ thực tế trên rất cần có
công trình điều tra chuyên sâu về các thảm thực vật đặc trƣng ở Việt Nam để làm cơ
sở cho xác định hệ số C phục vụ nghiên cứu xói mòn - thoái hoá đất.
Lê Mạnh Hùng đã ứng dụng mô hình SWAT trong tính toán xói mòn lƣu vực
Mekong [51], từ đó ƣớc tính tải lƣợng bùn cát về Kratie, đầu nguồn của vùng đồng
bằng châu thổ. Kết quả tính toán cho thấy SWAT dự báo dòng chảy trên sông
Mekong với độ chính xác khá tốt (chỉ số NSE nằm trong khoảng 0.67 ÷ 0.86, chỉ số
PBIAS khoảng 11.96 ÷ 22.55). Kết quả cũng cho thấy SWAT có khả năng ƣớc tính
tải lƣợng bùn cát trên lƣu vực với độ tin cậy chấp nhận đƣợc. Trên phần hạ lƣu vực
Mekong, vùng có suất bùn cát lớn nhất là vùng từ Luang Prabang cho đến
Mudkahan với suất bùn cát trung bình khoảng 289.000 tấn/ngày. Đây là vùng có địa
hình dốc, lƣợng mƣa trung bình hàng năm lớn. Tải lƣợng bùn cát trung bình năm
tính toán tại Kratie trong giai đoạn 2007 - 2011 là khoảng 162 triệu tấn/năm.
Trần Quốc Vinh đã nghiên cứu sử dụng viễn thám (RS) và hệ thống thông tin
địa lý GIS để đánh giá xói mòn đất huyện Tam Nông tỉnh Phú Thọ [20]. Kết quả
của nghiên cứu chỉ ra các hệ số xói mòn đất theo phƣơng trình mất đất phổ dụng
phù hợp với vùng đất gò đồi huyện Tam Nông tỉnh Phú Thọ.
Trƣơng Văn Cảnh đã nghiên cứu ảnh hƣởng xói mòn đất của lƣu vực sông Cu
34
Đê đến sản xuất nông nghiệp [21]. Dựa trên cơ sở lí luận và thực tiễn về xói mòn
đất và khả năng ứng dụng GIS, viễn thám, nghiên cứu tập trung vào việc phân tích
hiện trạng xói mòn đất trong phạm vi lƣu vực sông Cu Đê bằng các ứng dụng của
GIS, từ đó có những đánh giá tác động đến đất sản xuất nông nghiệp. Đồng thời đề
xuất các biện pháp bảo vệ đất, chống xói mòn đất. Việc sử dụng phƣơng trình mất
đất phổ dụng USLE kết hợp với GIS và viễn thám trong nghiên cứu xói mòn đất và
phân tích tác động đến sản xuất nông nghiệp đem lại hiệu quả cao, tiết kiệm thời
gian và chi phí. Kết quả ở dạng dữ liệu GIS nên có tính trực quan, dễ sử dụng và
hiệu chỉnh. Chúng ta có thể nhận biết những thông tin về xói mòn đất và các ảnh
hƣởng của nó đến đất sản xuất nông nghiệp một cách nhanh chóng và chính xác.
Tuy nhiên các hệ số K và P trong nghiên cứu trên đƣợc lấy ở các tài liệu tham khảo
nên phần nào ảnh hƣởng đến kết quả nghiên cứu.
1.4. Những khoảng trống trong nghiên cứu xói mòn và vận chuyển bùn cát –
Hƣớng nghiên cứu chính của luận án
1.4.1. Những khoảng trống trong nghiên cứu xói mòn và vận chuyển bùn cát
– Các kết quả nghiên cứu tính toán xói mòn và vận chuyển bùn cát trong
nƣớc chƣa đánh giá đƣợc sự phân bố lƣợng bùn cát bị xói mòn do các trận mƣa theo
thời gian và chƣa xem xét đến cấp độ dốc trên lƣu vực mà cấp độ dốc là một trong
những nguyên nhân chính gây xói mòn lƣu vực và ảnh hƣởng lớn đến khả năng vận
chuyển bùn cát trên lƣu vực.
– Các nghiên cứu ở Việt nam chủ yếu là ứng dụng các mô hình có sẵn để
tính toán xói mòn và vận chuyển bùn cát.
– Các mô hình toán nhƣ mô hình ANSWER [41], WEPP [45, 52],
EUROSEM [47, 48], SWAT... là những mô hình mang tính thƣơng mại, có giá
thành cao và mã nguồn đóng nên khi áp dụng vào lƣu vực sông ở Việt Nam sẽ
không tránh khỏi những sai số nhất định vì bộ thông số trong mô hình đƣợc hiệu
chỉnh theo số liệu ở nƣớc ngoài.
35
– Mô hình ANSWER chƣa xét đến xói mòn trong kênh; Thuật toán giải của
mô hình WEPP, SWAT… bằng phƣơng pháp giải tích nên phải đồng nhất độ dốc
của sƣờn dốc là không phù hợp với thực tế; Số liệu đầu vào lớn.
1.4.2. Định hƣớng nghiên cứu của luận án
Để khắc phục những hạn chế nêu trên, luận án đi sâu vào định hƣớng nghiên cứu:
- Nghiên cứu cơ sở lý thuyết mô hình hóa dòng chảy, quá trình xói mòn và vận
chuyển bùn cát trên lƣu vực vừa và nhỏ. Từ đó xây dựng sơ đồ hình thành dòng chảy,
vận chuyển bùn cát và sơ đồ mô phỏng đƣờng quá trình dòng chảy bùn cát trên lƣu vực.
- Nghiên cứu đề xuất thuật toán giải mới ổn định và hội tụ đối với phƣơng
trình liên tục dòng chảy, phƣơng trình động lƣợng và phƣơng trình xói mòn và vận
chuyển bùn cát.
- Trên cơ sở nghiên cứu, luận án sử dụng thuật toán lƣợc đồ sai phân hữu hạn Lax – Friedrichs có thêm trọng số thời gian và không gian; Chọn ngôn ngữ C ++
xây dựng mô hình toán mô phỏng xói mòn và vận chuyển bùn cát trên lƣu vực vừa
và nhỏ nhằm khắc phục những nhƣợc điểm nêu trên.
Kết luận chƣơng 1:
Tác giả tổng quan về những mô hình mô phỏng vận chuyển bùn cát trên lƣu
vực trên thế giới và ở Việt Nam, từ đó đã chỉ ra những tồn tại của một số mô hình
nhƣ ở Việt Nam chỉ ứng dụng mô hình, những mô hình trên thế giới mang tính chất
thƣơng mại với giá thành rất cao, mã nguồn đóng, bộ thông số đƣợc xây dựng dựa
trên đặc điểm vật lý lƣu vực của thế giới nên khi ứng dụng vào lƣu vực ở Việt nam
sẽ có sai số lớn. Vậy để khắc phục những tồn tại nêu trên tác giả nghiên cứu xây
dựng mô hình toán mô phỏng vận chuyển bùn cát trên lƣu vực vừa và nhỏ ở Việt
Nam.
36
CHƢƠNG 2. NGHIÊN CỨU CƠ SỞ KHOA HỌC XÂY DỰNG
MÔ HÌNH TOÁN MÔ PHỎNG VẬN CHUYỂN BÙN CÁT
TRÊN LƢU VỰC VỪA VÀ NHỎ
2.1. Cơ sở lý thuyết – Nghiên cứu đề xuất thuật toán giải
2.1.1. Cơ sở lý thuyết
2.1.1.1. Phƣơng trình thấm
Lƣơ ̣ng mƣa sinh dòng chảy hoă ̣c dòng chảy trƣ̣c tiếp trên bề mă ̣t là thàn h phần
, xói mòn và vận
quan tro ̣ng trong viê ̣c xây dƣ̣ng mô hình mô phỏng dòng chảy chuyển bùn cát trên lƣu vực . Lƣơ ̣ng mƣa sinh dò ng chảy đƣơ ̣c tính toán bằng cách trƣ̀ trƣ̣c tiếp các tổn thất trong các yếu tố thủ y văn hoă ̣ c lƣơ ̣ng mƣa đầu vào. Tổn thất bao gồm:
1. Lƣu trữ
2. Bốc thoát hơi
3. Điền trũng
4. Thấm.
Có ba phƣơng pháp cơ bản để mô hình hóa lƣợng mƣa sinh dòng chảy gồm : 1) Các tổ n thất kể trên đƣơ ̣c mô hình hóa mô ̣t cách riêng biê ̣t và đƣơ ̣c liên kết vớ i nhau để xác định lƣơ ̣ng mƣa sinh dò ng; 2) Các tổn thất đƣơ ̣c tổng hơ ̣p toàn bộ trong mô ̣t mô hình; 3) Tổn thất do lƣơ ̣ng mƣa đo ̣ng la ̣i , sƣ̣ thoát bốc hơi nƣớ c và tổn thất do tổn thất do thấm xuống nƣớ c chƣ́ a trong các vũng , ao trũng đƣơ ̣c bỏ qua và chỉ có đất đƣơ ̣c sƣ̉ du ̣ng để tính toán lƣơ ̣ng mƣa sinh dòng . Trong phạm vi nghiên luận án sử dụng phƣơng pháp số 3.
Quá trình thấm đƣợc tính toán từ phƣơng trình tính toán lƣợng tổn thất thấm
trong thời gian mƣa. Phƣơng trình thấm của Green – Ampt Mein – Larson (1973):
(2- 1)
37
trong đó:
fi : Cƣờng độ thấm thực (m/s)
Ke : Độ dẫn thuỷ lực bão hoà (m/s) ϕ : Ảnh hƣởng trạng thái xốp của đất (m3/m3) 0 : Độ bão hoà ban đầu (m3/m3)
c : Ảnh hƣởng độ căng của ống mao dẫn (m)
I : Độ sâu thấm luỹ tích (m)
Thành phần cơ bản trong mô hình hóa quá trình xói mòn và vận chuyển bùn
cát trên lƣu vực là thành phần dòng chảy . Các hệ phƣơng trình sóng động học đã đƣơ ̣c sƣ̉ du ̣ng nhƣ là sƣ̣ xấ p xỉ dòng chảy mô ̣t chiều đƣợc rú t go ̣n từ hê ̣ phƣơng trình chuyển động tuần hoà n nƣớc dƣớ i hầu hết các điều kiện (trạng thái chảy) của
dòng chảy tràn bề mặt và đối với một số trƣờng hợp liên quan trực tiếp đến nƣớc do
mƣa trong các kênh.
2.1.1.2. Phƣơng trình mô phỏng dòng chảy
1. Hệ phương trình mô phỏng dòng chảy trên bề mặt lưu vực
Phƣơng trình liên tục:
(2- 2)
Phƣơng trình đô ̣ng lƣơ ̣ng [53]:
(2- 3)
trong đó:
h : Độ sâu dòng chảy (m)
u : Vâ ̣n tốc dòng chảy trung bình (m/s) r : Lƣợng mƣa hiệu quả hay lƣợng dòng chảy nhập bên (m/s)
t : Thờ i gian tính toán (s)
38
x : Khoảng cách theo hƣớng dòng chảy (m) g : Gia tốc tro ̣ng trƣờ ng (m/s2)
S0 : Độ dốc sƣờn dốc hoặc độ dốc lòng sông (%)
Sf : Độ dốc ma sát hay độ dốc cản (%).
Với giả thiết trong phƣơng trình động lƣợng độ dốc sƣờn dốc hoặc độ dốc
lòng sông bằng với độ dốc ma sát hay độ dốc cản và đô ̣ dốc đƣờ ng mă ̣t nƣớ c đƣơ ̣c giả thiết là bằng với độ dốc của mặt phẳng [54]. Tƣ̀ đó hê ̣ phƣơng trình đô ̣ng lƣơ ̣ng đƣơ ̣c rú t go ̣n:
S0 = Sf (2- 4)
Kết hợp với phƣơng trình Manning:
u = αhm-1 (2- 5)
Trong đó α là hệ số phụ thuộc vào độ dốc sƣờn dốc, hệ số nhám của bề mặt
lƣu vực, hệ số nhám có trị số phụ thuộc vào sử dụng đất , m là hệ số liên quan đến đi ̣a hình và đô ̣ nhám bề mă ̣t.
Thay (2-5) vào phƣơng trình (2-2) ta có:
(2- 6)
Với điều kiê ̣n biên:
h(0,t) = 0 khi t ≥ 0 (2- 7)
Điều kiê ̣n ban đầu:
h(x,0) = 0 khi x ≥ 0 (2- 8)
Hê ̣ số nhám Mainning củ a dòng chảy là :
39
1/2
2/3Sf
u = (1/n)RH (2- 9)
trong đó:
RH : Bán kính thủy lực
n : Hệ số sức cản Mainning của sức cản dòng chảy.
2. Hệ phương trình mô phỏng dò ng chảy trong kênh/sông:
Phƣơng trình liên tục:
(2- 10)
Phƣơng trình động lƣợng:
(2- 11)
Trong đó A(x,t) là diện tích mặt cắt ngang của dòng chảy (m2), và q A là lƣu
lƣợng dòng chảy bên trên mỗi đơn vi ̣ chiều dài lòng kênh (m3/s).
Vớ i sƣ̣ xấp xỉ đô ̣ng ho ̣c, phƣơng trình (2-11) có thể đƣợc viết:
(2- 12)
Với Q(x,t) đƣợc biểu diễn nhƣ là mô ̣t hàm củ a A (x,t), khi đó phƣơng trình (2-10) đƣơ ̣c biểu diễn nhƣ sau:
(2- 13)
2.1.1.3. Phƣơng trình diễn toán xói mòn và vận chuyển bùn cát
1. Phương trình mô phỏng xói mòn và vận chuyển bùn cát trên lưu vực
Mô hình mô phỏng xói mòn và vận chuyển bùn cát đƣợc thiết lập dựa vào
phƣơng trình toán học mô phỏng các hiện tƣợng vật lý của quá trình xói mòn và vận
40
chuyển bùn cát. Cơ sở toán học của các hiện tƣợng vật lý là phƣơng trình liên tục
của Bennett [26]. Phƣơng trình liên tục mô phỏng xói mòn và vận chuyển bùn cát
đƣợc biểu diễn nhƣ sau:
(2- 14)
, với , k: hằng vớ i: ;trong đó : q =
số lƣu trữ khi đó ta có:
Hay q cho bở i phƣơng trình:
và vớ i (2- 15)
trong đó:
qs : Lƣu lƣợng bùn cát qua một đơn vị mặt cắt (kg/s.m) x : Khoảng cách dọc sông (m) ρs : Độ đục bùn cát (kg/m3) Cs : Nồng độ tập trung bùn cát (m3/m3) h : Độ sâu dòng chảy (m) θ : Lƣợng bùn cát bị xói mòn hoặc bồi lắng (kg/s.m2).
Có nhiều phƣơng pháp giải gần đúng cho phƣơng trình (2-15) nhƣ Euler ;
Euler cải tiến hoă ̣c Runger Kutta RK 4.v.v… Luận án chọn phƣơng pháp Runger Kutta RK4 để giải gần đúng q(tn) ≈ qn của phƣơng trình (2-15).
(2- 16)
và Vớ i
Ta tính đƣợc:
41
2. Phương trình mô phỏng xói mòn và vận chuyển bù n cá t trong kênh/sông: Theo Bennett [26], ta có:
(2- 17)
trong đó:
er : Tốc độ cuốn trôi bùn cát (sƣ̣ bong tách ) bở i dòng chảy trong kênh /sông (kg/m/s)
d : Tốc độ bồi lắng bùn cát (sƣ̣ chìm lắng) (kg/m/s) qs : Dòng chảy bùn cát bên từ các dòng chảy tràn bề mặt kế bên (kg/m/s).
Các biến khác đã đƣợc xác định ở các phần trƣớc.
(2- 18) Điều kiện biên: C(0,t) = C0(t) t ≥ 0
Điều kiện ban đầu: C(x,0) = 0 x ≥ 0 (2- 19)
Trong đó C 0(t) là nồng độ bùn cát đến từ các phần bên của kênh. Trong giai
đoa ̣n xói mòn và bồi lắng trong kênh, các mặt cắt dọc đáy kênh ở các vùng đất trũng đƣơ ̣c giả thiết là ổn đi ̣nh vớ i các mă ̣t cắt do ̣c củ a các rãnh trên các khu vƣ̣c đồi nú i khi mà chú n g có thể thay đổi nhanh chóng trong các khoảng thờ i gian ngắn . Trầm , phần lớ n tích đáy trong các dòng chảy ở vù ng đất trũng đƣơ ̣c giả thiết là thô hơn các kích cỡ hạt bùn (0.062mm).
Mô ̣t phƣơng trình nói chung , ban đầu đƣơ ̣c ph át triển cho khả năng vận
42
(er) bở i dòng
chuyển bù n cát đáy , đã đƣơ ̣c sƣ̉ du ̣ng để mô hình hóa sƣ̣ cuốn trôi chảy trong kênh [28, 55].
er = a(η – ηc)n khi η ≥ ηc (2- 20)
er = 0 khi η < ηc (2- 21)
trong đó:
và (2- 22)
trong đó:
: Hệ số cho sƣ̣ cuố n trôi vâ ̣t liê ̣u (Kg.m2/N1.5s)
η : Ứng suất trƣơ ̣t trung bình (N/m2) ηc : Ứng suất tớ i ha ̣n trung bình cho kích thƣớ c ha ̣t đa ̣i diê ̣n (N/m2) Sf : Độ dốc ma sát hay độ dốc cản (%) n : Hệ số Manning RH : Bán kính thủy lực γ : Trọng lƣợng riêng của nƣớc (N/m3) γs : Trọng lƣợng riêng của bùn cát (N/m3) ds : Đƣờng kính hạt (mm) δ : Hệ số phụ thuộc vào các đặc tính của chất lỏng và trầm tích (không
thƣ́ nguyên).
Mƣ́ c đô ̣ bồi lắng trầm tích (dòng hƣớng xuống) là tỷ lệ với nồng độ bùn cát và vâ ̣n tố c lắng chìm hiê ̣u quả củ a ha ̣t cát. Mƣ́ c đô ̣ bồi lắng đƣơ ̣c tính t heo phƣơng trình sau:
d = εTWVSC (2- 23)
trong đó :
C : Nồng độ bùn cát (m3/m3)
ε : Hê ̣ số phu ̣ thuô ̣c vào các đă ̣c điểm củ a bù n cát và dòng chảy
Vs : Vận tốc lắng chìm của hạt cát (m/s)
TW : Chiều rô ̣ng phía trên dòng chảy (m) và các biến khác đã đƣợc xác định ở trên.
43
2.1.2. Nghiên cứu đề xuất thuật toán giải
Hiện nay có nhiều phƣơng pháp giải hệ phƣơng trình diễn toán dòng chảy và
hệ phƣơng trình diễn toán xói mòn và vận chuyển bùn cát, nhiều tác giả đã sử dụng
lƣợc đồ sai phân hữu hạn 4 điểm dạng ẩn, dạng hiện hình chữ nhật, lƣợc đồ sai phân
6 điểm, lƣợc đồ sai phân 4 điểm hình thoi v.v… để giải. Trong nghiên cứu này luận
án dùng thuật toán lƣợc đồ sai phân hữu hạn Lax - Friedrichs thêm trọng số thời
gian và trọng số không gian gọi là lƣợc đồ sai phân hữu hạn Lax - Friedrichs có
trọng số (LFW) để giải hệ phƣơng trình diễn toán dòng chảy và hệ phƣơng trình
diễn toán xói mòn và vận chuyển bùn cát. Thuật toán đã đƣợc kiểm tra tính hội tụ
kết quả khá tốt đối với những phƣơng trình nêu trên. Thuật toán giải đƣợc trình bày
chi tiết dƣới đây.
2.1.2.1. Thuật toán giải hệ phƣơng trình mô phỏng dòng chảy trên bề mă ̣t lƣu vực
1. Thuật toán giải:
Sơ đồ sai phân hữu hạn Lax - Friedrichs (LF) đƣợc mô tả nhƣ hình (2-1). Với
xj = jΔx, tn = nΔt công thức sai phân Lax - Friedrichs tính gần đúng đạo hàm cấp
một của hàm u = u(x,t) tại điểm (xj, tn):
(2- 24)
Hình 2-1. Sơ đồ sai phân Lax - Friedrichs
44
Lƣợc đồ sai phân hữu hạn Lax- Friedrich có trọng số (LFW) đƣợc viết:
Áp dụng LFW cho phƣơng trình (2-6). Khi đó ta đƣợc phƣơng trình sai phân
của (2-6) nhƣ sau:
(2- 25)
trong đó:
θ : Trọng số không gian
ω : Trọng số thời gian.
Với ω = 0 ta có công thức hiện giải :
2. Tính hội tụ của thuật toán
Nghiên cứu tính ổn định của thuật toán sai phân LFW đối với phƣơng trình
(2-6), tức là tính ổn định của thuật toán sai phân LFW với phƣơng trình:
(2- 26)
trong đó với là hằng số.
45
Ta giả sử tại giá trị gần đúng:
Với là sai số tƣơng ứng. Khi đó điều kiện để thuật toán sai phân ổn định là:
(2- 27)
Áp dụng thuật toán sai phân LFW cho (2-26), ta đƣợc phƣơng trình sai phân
cho sai số :
Hay
(2- 28)
Theo cách làm của Rezzolla, L., (2011) [56], ta có:
(2- 29)
Với và là hàm phức phụ thuộc số nguyên dƣơng và còn
đƣợc gọi là nhân tử khuyếch đại. Khi đó điều kiện (2-27) tƣơng đƣơng với:
Thay (2-29) vào (2.28) với Δx = xj – xj-1 = xj+1 - xj ta đƣợc
(2- 30)
Sử dụng công thức Euler eiφ = cosθ + isinθ ta có đƣợc từ (2-30):
46
Từ đó ta có:
.
Hay
(2- 31)
Khi đó
(2- 32)
Nhƣ vậy điều kiện để thuật toán sai phân LFW cho (2-26) ổn định là:
Tƣơng đƣơng
(2- 33)
47
Với ω = 0, θ = 1 thì LFW trở thành sai phân Lax - Friedrichs, khi đó (2-33) trở
thành:
(2- 34)
Tức là ta tìm đƣợc điều kiện đủ để (2-33) xảy ra là , ta đƣợc điều kiện
nhƣ ở các phƣơng pháp sai phân khác [57].
Ở đây khi , điều kiện để phƣơng pháp sai phân hội tụ khi đó là:
Điều đó chỉ ra rằng, ta cần có mối quan hệ giữa bƣớc thời gian và không gian là
.
Xét với . Ta chứng tỏ khi thì (2-33)
luôn đúng, tức là phƣơng pháp LFW ổn định.
Thật vậy:
Ta luôn có:
. (*)
Hơn nữa:
. (**)
Và
48
Do đó , ta có:
Từ (*), (**) và (***), ta nhận đƣợc một điều kiện đủ để (2-33) xảy ra là:
.
Vậy chúng ta có hai kết luận về điều kiện đủ để phƣơng pháp sai phân LFW ổn định
là:
Khi thì cần điều kiện .
Khi không cần điều kiện của
2.1.2.2. Thuật toán giải hệ phƣơng trình mô phỏng dòng chảy trong kênh/sông:
Từ phƣơng trình:
(2- 35)
Điều kiện biên:
khi t ≥ 0 Q(0,t) = Q0(t) (2- 36)
Điều kiện ban đầu :
Q(x,0) = 0 khi x ≥ 0 (2- 37)
49
Với , ở đó là chu vi thấm, thì là hàm theo
cho bởi (2-12) trở thành.
(2- 38)
Giải gần đúng cho tại theo sơ đồ sai phân LFW, nhƣng theo (2-
38) ta chỉ cần giải cho bởi phƣơng trình (2-35). Khi đó ta đƣợc phƣơng trình sai
phân cho (2-35):
(2- 39)
trong đó :
(2- 40)
Xét trƣờng hợp , với , ta giải đƣợc tƣờng minh theo
và nhƣ sau:
(2- 41)
2.1.2.3. Thuật toán giải phƣơng trình mô phỏng xói mòn và vận chuyển bùn cát trên
lƣu vực
1. Thuật toán giải
Có nhiều phƣơng pháp giải phƣơng trình (2-14), nhƣng để hạn chế sai số trong
bài toán thực tế, nhiều tác giả dùng phƣơng pháp sai phân hữu hạn sơ đồ hiện hoặc
sơ đồ ẩn. Tác giả dùng sơ đồ sai phân Lax - Friedrichs, nhƣng thêm trọng số thời
50
gian và không gian để tính gần đúng đạo hàm cấp một của hàm u = u(x,t) tại điểm
(xj, tn) với xj = jΔx, tn = nΔt.
Sử dụng sơ đồ sai phân LFW vào (2-14) ta đƣợc phƣơng trình sai phân:
(2- 42)
theo và các đại lƣợng gần đúng của : Ta giải đƣơ ̣c
(2- 43)
Với nên ta có:
(2- 44)
b) Tính ổn định của thuật toán
Thuật toán LFW cũng ổn định đối với phƣơng trình (2-14), ta dễ dàng kiểm tra
tính ổn định của thuật toán bằng cách hoàn toàn tƣơng tự nhƣ các bƣớc trong quá trình
kiểm tra điều kiện ổn định của thuật toán sai phân LFW cho phƣơng trình (2-26).
Quá trình vận chuyển bùn cát có thể sinh ra bùn cát
(sƣ̣ bong tách ) và sự bồi lắng bù n cát (sƣ̣ lắng đo ̣ng ). Sƣ̣ hình thành bù n cát tr ên các sƣờn dốc do hai quá trình: Động năng của hạt mƣa làm bong tách các phần tử đất ; Ứng suất trƣơ ̣t . Xói mòn và bồi lắng có thể xảy ra đồng thời ở các mức độ khác nhau.
51
Hàm θ của phƣơng trình (2-14) đƣợc biểu diễn nhƣ sau:
θ = er – d + eI (2- 45)
trong đó:
er : Lƣợng bù n cát sinh ra do ƣ́ ng suất trƣơ ̣t (kg/m2/s) d : Lƣợng bùn cát bồi lắng (chìm xuống) (kg/m2/s) ei : Lƣợng bù n cát sinh ra do tác đô ̣ng mƣa rơi (kg/m2/s).
Với er = KR(η)1.5 (2- 46)
trong đó :
KR : Hệ số bong tách đất của ứng suất trƣợt (Kg.m/N1.5.s)
η : Giá trị ứng suất trƣợt trung bình hiệu quả.
Sự bồi lắng bù n cát đƣơ ̣c tính theo phƣơng trình sau [58]:
d = εVSCa (2- 47)
trong đó:
ε : Hệ số phụ thuộc vào các đă ̣c điểm củ a bù n cát và dòng chảy
VS : Vận tốc lắng chìm của hạt cát (phần tƣ̉ ) (m/s)
Ca : Nồng độ bùn cát thực tế gần đáy (kg/m3).
d = εVSRCC (2- 48)
trong đó RC = Ca(x,t)/C(x,t) là hệ số mặt cắt nồng độ.
Khi không phân biê ̣t giƣ̃a bù n cát đáy và bù n cát lơ lƣ̉ ng (RC = 1.0) thì:
d = εVSC (2- 49)
đƣơ ̣c tính toán theo phƣơng trình Rubey:
(2- 50)
52
Với:
(2- 51)
trong đó:
γ : Trọng lƣợng riêng củ a nƣớ c (N/m3) γs: Trọng lƣợng riêng của bùn cát (N/m3) ν : Hệ số nhớt động học của nƣớc (m2/s)
dS : Đƣờng kính hạt (mm).
Hàm sinh ra bùn cát do mƣa là hàm số đơn giản biểu diễn sự bong tách do tác
đô ̣ng của mƣa kết hơ ̣p vớ i cƣờ ng đô ̣ mƣa i (x,t) là thƣớ c đo củ a sƣ̣ xói mòn [46, 55].
Nếu cƣờ ng đô ̣ mƣa phân bố đều trên toàn khu vƣ̣c thì i(x,t) = i(t).
eI = ai2 (2- 52)
Trong đó : a là hê ̣ số đƣơ ̣c xác định thông qua thí nghiê ̣m . Phƣơng trình (2-51) phản ánh tốc đô ̣ vâ ̣n chuyển trầm tích do dòng chảy nông ở các mái dốc [59]. Lane và Shirley đƣa ra phƣơng trình đơn giản cho tốc đô ̣ cuốn trôi bù n cát trên các mái dố c là:
(2- 53)
trong đó:
KI : Hệ số đo sự bong tách đất bở i tác đô ̣ng mƣa rơi (kg-s/m3)
(x,t) có thể r : Hê ̣ số củ a tốc đô ̣ sinh dòng chảy đối vớ i cƣờ ng đô ̣ mƣa i
đƣơ ̣c thể hiê ̣n nhƣ là mô ̣t thƣớ c đo củ a cƣờ ng đô ̣ dòng chảy chuẩn hóa sƣ̣ vâ ̣n
chuyển trầm tích bở i dòng chảy nô ng. Chú ý rằng khi r (x,t) = 0 (giai đoa ̣n trƣớ c khi
tạo ra lƣợng mƣa sinh dòng ) hoă ̣c khi i(x,t) = 0 (giai đoa ̣n sau khi ta ̣o ra lƣơ ̣ng mƣa
sinh dòng ) sẽ không có sƣ̣ cuốn trôi bù n cát do mƣa và khi r (x,t) = i(x,t), tốc đô ̣
cuốn trôi bù n cát do mƣa là không bị giới hạn bởi lƣợng mƣa sinh dòng.
53
Khi đó (2-14) đƣơ ̣c viết la ̣i:
(2- 54)
với: của tại (tức là đã giải gần đú ng trƣớ c đó . Cần tìm gần đú ng
).
Điều kiê ̣n biên trên:
khi t ≥ tP (2- 55)
Điều kiê ̣n ban đầu:
khi x ≥ 0 (2- 56)
trong đó tP là thời gian tạo thành các hồ , vũng nƣớc nhỏ.
Giải phƣơng trì nh (2-14) là để tìm ra nồng độ bùn cát tại bƣ ớc thời gian và
không gian đối vớ i các giá tri ̣ đã biết.
Sƣ̉ du ̣ng công thƣ́ c sai phân LFW cho phƣơng trình (2-54). Khi đó ta đƣợc
phƣơng trình sai phân của (2-54) nhƣ sau:
(2- 57)
Lấy theo và các đại lƣợng gần đúng của ta giải đƣơ ̣c
nhƣ sau:
54
(2- 58)
(2- 59)
2.1.2.4. Thuật toán giải phƣơng trình mô phỏng xói mòn và vận chuyển bù n cát
trong kênh/sông
Áp du ̣ng lƣợc đồ sai phân LFW vào phƣơng trình mô phỏng xói mòn và vận
chuyển bù n cát trong kênh/sông [26]:
(2- 60)
(2- 61)
theo và các đại lƣợng gần đúng của Từ đó ta giải đƣơ ̣c
nhƣ sau:
55
(2- 62)
2.2. Xây dựng các thành phần của mô hình
Trên cơ sở phân tích lý thuyết sự hình thành dòng chảy và vận chuyển bùn cát
trên lƣu vực hình 2-2. Luận án xây dựng sơ đồ khối mô phỏng quá trình xói mòn và
vận chuyển bùn cát trên lƣu vực sông nhƣ hình 2-3.
Hình 2-2. Sơ đồ hình thành dòng chảy và vận chuyển bùn cát trên lưu vực
56
Hình 2-3. Sơ đồ tính toán xói mòn và vận chuyển bùn cát
2.2.1. Quá trình liên rãnh
Quá trình phá vỡ và vận chuyển bùn cát đều xảy ra trên diện tích giữa các
57
rãnh. Các hạt đất bị bong tách do tác động của hạt mƣa đƣợc vận chuyển tới các
rãnh bởi lớp dòng chảy mỏng giữa các rãnh. Hầu hết lƣợng bùn cát tới các rãnh là
do dòng chảy giữa các rãnh mang đến, chỉ có một phần nhỏ bùn cát văng trực tiếp
tới các rãnh.
2.2.1.1. Khả năng xói mòn liên rãnh
Độ nhạy sự tách rời đất do tác động của mƣa sẽ khác nhau ứng với từng loại
đất vì chúng có các đặc trƣng vật lý, hóa học khác nhau. Các đặc trƣng đất ảnh
hƣởng tới khả năng xói mòn.
Phƣơng trình xói mòn liên rãnh là:
Ei= 0.0138 Kii2 (2- 63)
trong đó:
Ei : Tỷ lệ tách rời (kg/m2 h) Ki : Hệ số khả năng xói mòn đất (kg.h/Nm2) i : Cƣờng độ mƣa (mm/h).
Độ dốc liên rãnh đóng vài trò quan trọng trong việc xói mòn và vận chuyển
bùn cát, do vậy đã đƣa ra phƣơng trình độ dốc [60]:
Si = 2.96 (sinθ)0.79 + 0.56 (2- 64)
trong đó:
Si : Hệ số độ dốc liên rãnh θ : Góc dốc.
Sau đó Harmon và Meyer đƣa thêm thành phần thảm phủ [60], kết hợp với
phƣơng trình (2-63) xây dựng một phƣơng trình xói mòn liên rãnh:
eff [2.96 (sinθ)0.79 + 0.56] Ci
Di = 0.0138Ki i2 (2- 65)
trong đó:
eff : Biểu diễn sự ảnh hƣởng của tán cây tới khả năng xói mòn
i2
Ci : Hệ số thảm phủ - quản lý đối với phá vỡ và tách rời liên rãnh.
58
Mối quan hệ giữa ieff và Ci dựa vào loại ảnh hƣởng I, II và III theo
Wischmeier [61]. Loại I là tác động trên mặt đất chủ yếu từ tán cây; loại II là ảnh
hƣởng lớp phủ bề mặt đất, cây xanh còn lại và loại III là ảnh hƣởng sát mặt do rễ
cỏ, đã cải thiện cấu trúc đất và các thành phần còn lại. Sự ảnh hƣởng của loại I đƣợc
phản ánh trong ieff; trong khi đó sự ảnh hƣởng loại II và loại III phản ánh trong Ci.
Các giá trị nhân tố thành phần đƣợc nhân với nhau để nhận đƣợc giá trị Ci.
Hệ số C thay đổi theo thực hành nông nghiệp. Giá trị này bao gồm các tác
động của che phủ, năng suất, độ dài vụ mùa, tập quán canh tác. Hệ số C có thể tham
khảo trong bảng tính C của hội khoa học đất quốc tế đã xây dựng ở bảng 1-1.
Trong tính toán xói mòn và vận chuyển bùn cát trên lƣu vực mà xét đến toàn
bộ các hệ số trong các phƣơng trình trên là khó khả thi vì không có đủ số liệu. Để
đơn giản hóa tác giả chọn phƣơng trình sau:
Eil = ail.iket (2- 66)
trong đó:
ail : Hệ số xói mòn liên rãnh đƣợc xác định bằng cách hiệu chỉnh giữa thực
đo và tính toán trong mô hình
i : Cƣờng độ mƣa (mm/h)
ket : Hệ số khả năng xói mòn đất do tác động của hạt mƣa đƣợc xác định
bằng cách hiệu chỉnh giữa thực đo và tính toán trong mô hình
Eil : Xói mòn liên rãnh (kg/m2.h).
Hệ số ket đƣợc tính toán trên cơ sở tích hợp hệ số Cslr, mà Cslr là tỷ số giữa
lƣợng đất mất trên một đơn vị diện tích có lớp phủ thực vật và lƣợng đất mất trên
một diện tích trống tƣơng đƣơng với sự quản lý của con ngƣời trong phƣơng trình
USLE.
2.2.1.2. Vận chuyển bùn cát liên rãnh
Dòng chảy liên rãnh là dòng chảy tràn bề mặt. Dòng chảy tràn bề mặt có thể
59
dòng chảy tầng hoặc dòng chảy rối phụ thuộc vào hệ Reynol. Mức độ xói mòn bị
giới hạn bởi khả năng vận chuyển bùn cát liên rãnh. Các phƣơng trình biểu diễn khả
năng vận chuyển bùn cát dựa vào ứng suất trƣợt. Mối tƣơng quan cơ bản của
phƣơng trình khả năng vận chuyển bùn cát nhƣ sau:
Tcl = At (η-ηcr)1.5 (2- 67)
trong đó:
Tcl : Sức tải bùn cát liên rãnh (kg/m.h)
η : Ứng suất trƣợt (N/m2)
ηer : Ứng suất trƣợt tới hạn (N/m2).
Các hệ số At và er là các hàm của đƣờng kính và mật độ hạt.
Ứng suất trƣợt η do dòng chảy và năng lƣợng hạt mƣa tạo ra. Khi đó phƣơng trình
(2-67) trở thành:
eff Cit]
Tcl = γybs[(yb/yp). ai2 (2- 68)
trong đó:
a : Hệ số ƣớc lƣợng đƣợc
yb : Độ dày dòng chảy mặt
yb/yp : Tỷ số giữa độ dày của dòng chảy trên lƣu vực và độ sâu dòng chảy
trong các vùng trũng.
ieff : Cƣờng độ ảnh hƣởng tác động của hạt mƣa và Cit đƣợc tính nhƣ sau:
Cit = exp {-0.21[(yp/yb)-1]}1.18 (2- 69)
Phƣơng trình (2-68) không xét tới hệ số nhám, mà hệ số nhám ảnh hƣởng tới
sự phân bố bùn cát liên rãnh, do vậy trong mô hình này tác giả sử dụng phƣơng
trình của Meyer và Wischmeier [62] để xác định khả năng vận chuyển bùn cát liên
rãnh nhƣ sau:
60
5/3.q5/3
Tcl = acl. Io (2- 70)
trong đó:
Tcl : Sức tải bùn cát liên rãnh (kg/m.h)
I0 : Độ dốc sƣờn dốc (%) q : Lƣu lƣợng dòng chảy trên một đơn vị chiều dài (m3/m.s)
acl : Hệ số vận chuyển bùn cát liên rãnh.
2.2.2. Quá trình xói mòn rãnh
2.2.2.1. Xói mòn rãnh
Lƣợng bùn cát bị xói mòn trong rãnh trên một diện tích đƣợc mô hình hóa
bằng việc mô tả sự xói mòn trong từng rãnh riêng biệt [46]. Lúc đầu việc này dƣờng
nhƣ không khả thi bởi vì tồn tại rất nhiều các rãnh đặc trƣng và sự biến đổi về hình
dạng hình học của các rãnh. Tuy nhiên trong trƣờng hợp đặc biệt của các rãnh giữa,
phân tích sự xói mòn trong một số rãnh đƣợc chọn là khả thi và thích hợp với thực
tế.
Xói mòn rãnh và dòng chảy đƣợc giả thiết phân bố đồng đều, tuy nhiên về mặt
vật lý dòng chảy và xói mòn tập trung trong các kênh nhỏ, từ đó phƣơng trình xói
mòn rãnh đƣợc xác định:
Ei = a(η - ηcr)b (2- 71)
trong đó:
Ei : Tốc độ phá vỡ tiềm năng của xói mòn rãnh (kg/m2.h)
η : Ứng suất trƣợt của dòng chảy giả thiết là rộng và nông (N/m2)
ηcr: Ứng suất trƣợt tới hạn (N/m2).
Wischmeier và cộng sự đƣa ra phƣơng trình xói mòn rãnh [36] nhƣ sau:
Ei = 83.7 Krη1.5Cr (2- 72)
61
trong đó:
Ei : Tốc độ phá vỡ tiềm năng của sự xói mòn rãnh (kg/m2.h) Kr : Hệ số xói mòn đất cho xói mòn rãnh (kg.h/N.m2) η : Ứng suất trƣợt trung bình của dòng chảy rộng và nông (N/m2)
Cr : Hệ số lớp phủ tính trong xói mòn rãnh.
Ngoài hai công thức trên Foster và cộng sự giả thiết dòng chảy động học là ổn
định và sử dụng phƣơng trình dòng chảy đều của Darcy - Weishbach, xây dựng
phƣơng trình xói mòn rãnh cho sƣờn dốc có chiều dài x nhƣ sau:
qsrc = aγ1.5(f/8g)0.5sζx2KrCr/2 (2- 73)
trong đó:
qsrc : Tải trọng bùn cát của sƣờn dốc với độ dài x
γ : Khối lƣợng riêng của nƣớc (kg/m3)
f : Hệ số ma sát Darcy- Weishbach g : Gia tốc trọng trƣờng (m/s2)
s : Sin của góc dốc (độ)
ζ : Tốc độ mƣa vƣợt thấm (m/s).
Các công thức trên khá phức tạp, phụ thuộc nhiều vào nhân tố nên tác giả chọn
một phƣơng trình khác có mối quan hệ giữa các nhân tố đơn giản hơn để tính toán
trong mô hình nhƣ sau:
Ei = Ciq s K Cslr (2- 74)
trong đó:
Ei : Tốc độ phá vỡ tiềm năng của sự xói mòn rãnh (kg/m2.h)
Ci : Hệ số xói mòn rãnh đƣợc xác định từ việc hiệu chỉnh mô hình q : Lƣu lƣợng dòng chảy (m2/phút)
s : Sin của góc dốc (độ)
K : Hệ số xói mòn đất lấy từ USLE
Cslr : Tỷ lệ mất đất từ USLE.
62
2.2.2.2. Vận chuyển bùn cát trong rãnh
Bùn cát lơ lửng đƣợc phân bố đều trên toàn bộ độ sâu dòng chảy và ít bị lắng
đọng, bùn cát lơ lửng di chuyển xa hơn so với bùn cát đáy trƣớc khi nó lắng đọng,
đặc biệt là các hạt bùn cát mịn.
Phƣơng trình sức tải bùn cát trong rãnh rất quan trọng trong mô hình xói mòn
và vận chuyển bùn cát. Sự lắng đọng bùn cát phụ thuộc trực tiếp vào lƣợng bùn cát
bị phá vỡ ở liên rãnh.
Thông thƣờng bùn cát gồm các hạt có kích thƣớc và mật độ khác nhau. Với
mục đích đơn giản hóa, nên xem kích thƣớc và mật độ bùn cát là đồng nhất. Tuy
nhiên, trong các bài toán tính toán sự phân chia của bùn cát trong suốt quá trình
lắng đọng có vai trò quan trọng, do vậy việc lựa chọn phƣơng trình sức tải bùn cát
bùn cát cho phù hợp với tính toán là rất cần thiết cho mô hình.
Sức tải bùn cát trong rãnh phụ thuộc vào dòng chảy mặt, độ dốc của rãnh, lực
kéo của các hạt bùn cát nên Beasley đã xây dựng phƣơng trình sức tải bùn cát [41]
nhƣ sau:
Tc= 146.I0.q1/2 đối với q < 0.046 (m2/s)
Tc= 14600.I0.q2 đối với q > 0.046 (m2/s)
trong đó:
Tc : Sức tải bùn cát (kg/m.s)
q : Lƣu lƣợng dòng chảy bề mặt trên một đơn vị chiều dài (m2/phút)
I0 : Độ dốc của rãnh (%).
Hagen (1993) đƣa ra phƣơng trình tính nhƣ sau:
63
trong đó:
DEP : Tỷ lệ bồi lắng bùn cát (kg/s.m2)
Vx : Vận tốc chìm của các hạt (m/s) Tc : Sức tải bùn cát (kg/m.s) Q : Lƣu lƣợng dòng chảy bề mặt trên một đơn vị chiều dài (m2/s)
qs : Lƣu lƣợng bùn cát (kg/s.m).
Có nhiều phƣơng trình đƣợc xây dựng để tính sức tải bùn cát nhƣng luận án
chỉ nêu ra một số phƣơng trình tiêu biểu, mỗi phƣơng trình điều có tính ứng dụng
phù hợp với mỗi mô hình. Tác giả chọn phƣơng trình (2-75) [47] để đƣa vào mô
hình tính toán:
Tc = 3600.ac.q.ρs (2- 75)
trong đó:
Tc : Sức tải bùn cát (kg/m.s) q : Lƣu lƣợng dòng chảy bề mặt trên một đơn vị chiều dài (m2/s) ac : Hệ số sức tải bùn cát rãnh. ρs : Mật độ khối của các hạt bùn cát.
Cho đến nay, khi làm mô hình chƣa chỉ ra đƣợc phƣơng trình nào là ƣu điểm
nhất. Tất cả các phƣơng trình tính sức tải bùn cát đều đòi hỏi phải đƣợc hiệu chỉnh
và tối ƣu hóa thông qua tài liệu thực đo của dòng chảy tràn trên bề mặt.
2.2.3. Quá trình lòng kênh/sông
2.2.3.1. Sức tải bùn cát trong kênh/sông
Sức tải bùn cát đƣợc tính theo phƣơng trình của Parsons và nnk (1992) [63]:
(2- 76)
trong đó:
qsl : Sức tải bùn cát (kg/m.h)
64
Eil : Xói mòn liên rãnh (kg/m2.h) A : Diện tích (m2)
: Bƣớc thời gian tính toán (h) VQ : Tổng lƣợng dòng chảy (m3) q : Lƣu lƣợng dòng chảy bề mặt trên một đơn vị chiều dài (m2/phút).
2.2.3.2. Vận chuyển bùn cát trong kênh/sông
Vận chuyển bùn cát trong sông đƣợc mô phỏng theo phƣơng trình của
Engelund và Hansen [64] nhƣ sau:
(2- 77)
trong đó:
TCch : Sức tải bùn cát trong sông (kg/m.h)
I0 : Độ dốc đáy sông (%)
d : Đƣờng kính hạt trung bình (mm)
h : Độ sâu dòng chảy (m)
ρs : Mật độ khối của các hạt bùn cát.
Từ sơ đồ khối hình 2-3, tác giả sử dụng ngôn ngữ C++ để xây dựng mô hình
mô phỏng dòng chảy và vận chuyển bùn cát trên lƣu vực vừa và nhỏ.
Trong chƣơng trình tính toán xói mòn và vận chuyển bùn cát trên lƣu vực gồm
mô đun tính toán khả năng xói mòn liên rãnh, vận chuyển bùn cát liên rãnh, xói
mòn rãnh, vận chuyển bùn cát trong rãnh, sức tải cát trong kênh/sông. Nội dung các
mô đun đƣợc trình bày trong phần phụ lục của luận án.
2.3. Phân tích lựa chọn ngôn ngữ xây dựng mô hình
Ngôn ngữ lập trình C++ là ngôn ngữ lập trình cấp cao, đƣợc sử dụng rất phổ biến để lập trình hệ thống cùng với Assembler và phát triển các ứng dụng. C++ là một ngôn ngữ lập trình hệ thống rất mạnh và rất “mềm dẻo”, có thƣ viện chứa rất nhiều hàm (function) đã đƣợc tạo sẵn. Ngƣời lập trình có thể tận dụng các hàm này
65
để giải quyết các bài toán mà không cần tạo mới. Hơn thế nữa, ngôn ngữ C++ hỗ trợ rất nhiều phép toán nên phù hợp cho việc giải quyết các bài toán kỹ thuật có nhiều công thức phức tạp. Ngoài ra, C++ cũng cho phép ngƣời lập trình tự định nghĩa thêm các kiểu dữ liệu trừu tƣợng khác. Chính vì vậy, tác giả chọn C++ để xây dựng mô hình mô phỏng dòng chảy và vận chuyển bùn cát trên lưu vực vừa và nhỏ.
2.4. Cấu trúc và chức năng của một số chƣơng trình con
2.4.1. Cấu trúc và một số mô đun chính
2.4.1.1 Cấu trúc chƣơng trình: Có hai thành phần cơ bản
1. Thành phần thủy văn: Xƣ̉ lý dƣ̃ liê ̣u mƣa đầu vào và tính toán lƣợng mƣa sinh dòng qua phƣơng trình Green - Ampt - Lason (1973).
2. Thành phần xói mòn và bồi lắng : Tính lƣợng bùn cát sinh ra do tác đô ̣ng của
mƣa và dòng chảy; Tính lƣợng bùn cát bồi lắng trên lƣu vực và trong kênh/sông thông qua các phƣơng trình liên tu ̣c bù n cát vớ i thuâ ̣t toán lƣơ ̣c đồ sai phân hƣ̃u ha ̣n LFW.
2.4.1.2 Một số mô đun chính
1. Mô đun tính toán xói mòn và vận chuyển bùn cát liên rãnh:
Hình 2-4. Sơ đồ khối tính toán xói mòn và vận chuyển bùn cát liên rãnh
66
viết đoạn câu
Trên cơ sở sơ đồ khối hình 2-4 tác giả dùng ngôn ngữ C ++
lệnh chính nhƣ sau:
void CErosion::InterrillErosion( int nr )
{
if( pErosion[nr].OFRQ > 0 ) {
float RhoF = RhoW;
float RhoR2 = ( RhoS - RhoF ) / RhoF;
float Visco = pErosion[nr].Visco;
float LI = pErosion[nr].LI;
float ET = pErosion[nr].ET;
float OFRQ = pErosion[nr].OFRQ;
float Neff = pErosion[nr].INTPH / (1000*3600);
float D50 = pErosion[nr].D50;
pErosion[nr].VS = (float) ( 11.0f * Visco * ( sqrt( 1.0f + 0.01*D*D*D
) - 1.0f pErosion[nr].EI = ET * OFRQ / ( 2 * pErosion[nr].VS * LI );
pErosion[nr].EI = ET / ( 1 + pErosion[nr].VS/Neff );
if( pErosion[nr].EI > 1E16 || pErosion[nr].EI < -1E16 )
int t=1;
}
}
void CErosion::Deposition( int nr )
{
if( pErosion[nr].OFRQ > 0 ) {
float QSI = pErosion[nr].QSI;
float VS = pErosion[nr].VS;
float OFRQ = pErosion[nr].OFRQ;
float TC = pErosion[nr].TC;
float alpha= 0.5f * VS / OFRQ;
float alpha_eros = 0.5f;
float alpha = alpha_eros * VS / OFRQ;
pErosion[nr].DEP = alpha * ( TC - QSI );
67
if( pErosion[nr].DEP > 0 )
pErosion[nr].DEP = 0;
}
}
2. Mô đun tính toán xói mòn và vận chuyển bùn cát rãnh
viết đoạn câu lệnh
Trên cơ sở sơ đồ khối hình 2-5, tác giả dùng ngôn ngữ C ++
chính nhƣ sau:
void CErosion::RillErosion( int nr )
{
if( pErosion[nr].OFRQ > 0 ) {
float CS = pErosion[nr].CS;
float OFRQ = pErosion[nr].OFRQ;
float KFAK = pErosion[nr].Kfak;
float CFAK = pErosion[nr].Cfak;
float RR = pErosion[nr].RR;
float SinALPHAW = pErosion[nr].SinALPHAW;
pErosion[nr].ER = RR * OFRQ * SinALPHAW * KFAK * CFAK;
}
else
pErosion[nr].ER = 0;
}
void CErosion::Deposition( int nr )
{
float VS = pErosion[nr].VS;
float QSI = pErosion[nr].QSI;
float OFRQ = pErosion[nr].OFRQ;
float TC = pErosion[nr].TC;
pErosion[nr].DEP = VS * ( TC - QSI );
}
68
Hình 2-5. Sơ đồ khối tính toán xói mòn và vận chuyển bùn cát rãnh
3. Mô đun tính lƣợng bùn cát vận chuyển trong sông đến cửa ra của lƣu vực
Hình 2-6. Sơ đồ khối tính lượng bùn cát vận chuyển trong sông đến cửa ra
viết đoạn câu lệnh
Trên cơ sở sơ đồ khối hình 2-6, tác giả dùng ngôn ngữ C ++
chính nhƣ sau:
69
void CErosion::Erosion( int nr, int time, VARSTRUCT *pVar )
{
pErosion[nr].ER = pErosion[nr].EI = pErosion[nr].ET = pErosion[nr].QSI = 0; pErosion[nr].TC = pErosion[nr].DEP = pErosion[nr].Phi = 0; pErosion[nr].QSGR = pErosion[nr].WEIGHT = pErosion[nr].sedimenttransport = 0;
pErosion[nr].OFRQ = pVar[nr].qo1 / pErosion[nr].B; pErosion[nr].INTPH = pVar[nr].rain / delta_t;
pErosion[nr].PD = pVar[nr].rain; pErosion[nr].Qab = pVar[nr].qab / pVar[nr].itseg; pErosion[nr].VAUQ = pVar[nr].qo1 * 3600 * delta_t; if( pErosion[nr].OFRQ > 0 )
int t=1;
if( pErosion[nr].OFRQ < 0 ) int t=1;
if( pErosion[nr].Qab > 0 ) int t=1;
if( pErosion[nr].Qab < 0 ) int t=1;
if( pErosion[nr].QSI > 0 ) Sedimenttransport( nr, timestep );
if( pErosion[nr].Phi > 0 ) int t=1;
if( pErosion[nr].Phi < 0 ) int t=1;
float Sed = pErosion[nr].Phi * pErosion[nr].AREA; if( pErosion[nr].Qab > 0 &&Trans> 0 ) {
float sedtran = pErosion[nr].sedtran[0].QSGI; float time = m_pVar[nr].tchanel;
} if(sedtran> 0 ) { float ffactor = 1; float St = 0;
70
}
concentration= 0;
concentration = (float) 1E+12;
Tak[0].trans.ChangeAt( time, pVar[nr].qo1, concentration );
int TakNr = pStat->nr; pTak->time.ChangeAt( time, pVar[nr].qo1, concentration );
int tflrec = m_pTfl[nr].hydrologystations; CStation *pTak= GetTakPtr( m_pTfl[tflrec].hydrologystations Nr ); if( pStat ) { }
int TakNr = pTak->nr;
int t=1;
if( pTak ) { if( trans> 0 ) float trans= m_pVar[nr].qab / m_pVar[nr].itseg; float trans = concentration / ( chanelvolume*3600 ); }
Trans( nr ); float concentration = trans / (pVar[nr].qo1*3600); if(concentration<= 0 ) if(concentration> 1E+15 ) if( Tak != NULL ) if( Tak != NULL ) { } if( Tak != NULL && m_pTfl[nr].hydrology stations Nr > 0 ) { CStation *pTak = GetTakPtr( m_pTfl[nr].hydrologystations Nr ); } if( pErosion[nr].EI || pErosion[nr].ET || pErosion[nr].ER || pErosion[nr].DEP
pErosion[nr].S_ER += pErosion[nr].ER * delta_t;
pErosion[nr].S_EI += pErosion[nr].EI * delta_t; pErosion[nr].S_ET += pErosion[nr].ET * delta_t; pErosion[nr].S_DEP+= pErosion[nr].DEP * delta_t; pErosion[nr].S_PHI+= pErosion[nr].Phi * delta_t;
}
pErosion[nr].Channel += channelvolume; || pErosion[nr].Phi ) { }
71
2.4.2. Chức năng của một số chương trình con
Đo ̣c các dƣ̃ liê ̣u miêu tả các đă ̣c trƣng củ a lƣu vực. Tính toán các hệ số thấm trên các mặt phẳng bằng cách sử dụng phƣơng
trình Green - Ampt - Lason và ta ̣o ra thành phần lƣơ ̣ng mƣa sinh dòng.
Thƣ̣c hiê ̣n viê ̣c ghi chép (ghi lại thông tin ) vào nơi lƣu trữ thông tin tạm
thờ i củ a máy tính.
Giải các phƣơng trình ẩn phi tuyến tổng thể bằng việc sử dụng sơ đồ lặp
Newton - Rapson.
Tính lƣợng bùn cát sinh ra dƣới tác động của mƣa và ứng suấ
t trƣơ ̣t tƣ̀ dòng chảy mặt; Tính lƣợng bùn cát bồi lắng trên các tiểu lƣu vực bằng cách sƣ̉ du ̣ng lƣơ ̣c đồ sai phân hƣ̃u ha ̣n LFW với phƣơng trình liên tục bùn cát.
Tính lƣợng bùn cát sinh ra bởi ứng suất trƣợt do dòng chảy tập trung gây
ra; Và tính lƣợng bùn cát bồi lắng trong các kênh /sông bằ ng cách sƣ̉ du ̣ng lƣơ ̣c đồ sai phân hƣ̃u ha ̣n LFW.
Ghi ra các tâ ̣p tin đầu ra
2.5. Giao diện sử dụng chƣơng trình
Hình 2-7. Màn hình khởi động chương trình
72
Hình 2-8. Giao diện mở tệp dữ liệu
Hình 2-9. Giao diện chạy chương trình
73
Kết luận chƣơng 2
Trong chƣơng này, tác giả dùng lƣợc đồ sai phân hữu hạn Lax- Frierichs có
thêm trọng số thời gian và không gian (LFW) để giải phƣơng trình dòng chảy và
phƣơng trình vận chuyển bùn cát trên lƣu vực. Tác giả đã chứng minh tính hội tụ
của thuật toán đối với phƣơng trình dòng chảy và phƣơng trình vận chuyển bùn cát.
Trên cơ sở đó, mô hình đƣợc xây dựng dựa trên thuật toán LEW để mô phỏng
dòng chảy và vận chuyển bùn cát trên lƣu vực vừa và nhỏ.
74
CHƢƠNG 3. THỬ NGHIỆM MÔ HÌNH ĐỂ MÔ PHỎNG DÒNG CHẢY VÀVẬN CHUYỂN BÙN CÁT CHO MỘT SỐ LƢU VỰC VỪA VÀ NHỎ
3.1. Số liệu đầu vào cho mô hình
3.1.1. Tạo cơ sở dữ liệu
- Hệ thống thông tin địa lý (GIS) sẽ số hóa bản đồ thổ nhƣỡng, thảm phủ, địa
hình và hệ thống sông. Mô hình hóa độ cao của lƣu vực dựa vào bản đồ địa hình.
- Nhờ vào các phép phân chia, chồng ghép bản đồ theo sự đồng nhất về địa
hình, hƣớng dòng chảy, thổ nhƣỡng và tính chất thảm phủ của lƣu vực để phân chia
lƣu vực thành các tiểu lƣu vực, từ đó xác định đƣợc hƣớng dòng chảy của các tiểu
lƣu vực.
- Tạo ra tệp số liệu mƣa, lƣu lƣợng dòng chảy thực đo, lƣu lƣợng bùn cát lơ
lửng thực đo và các thông số của mô hình.
3.1.2. Chạy mô hình
- Khi nhập dữ liệu vào mô hình, chạy chƣơng trình ta có kết quả mô phỏng
dòng chảy, kết quả mô phỏng vận chuyển bùn cát trên lƣu vực dƣới dạng đồ thị.
- Trên cơ sở đó hiệu chỉnh các thông số của mô hình theo số liệu thực đo và
đƣợc bộ thông số của mô hình.
3.2. Thử nghiệm mô hình để mô phỏng xói mòn và vận chuyển bùn cát
3.2.1. Lưu vực Nậm Sập
3.2.1.1. Vị trí địa lý
Lƣu vực Nậm Sập thuộc địa phận các huyện Mộc Châu, Yên Châu, Bắc Yên, Mai Sơn, Vân Hồ tỉnh Sơn La, có diện tích tự nhiên 1085,52km2, chiếm 7,66% tổng diện tích toàn tỉnh. Nằm trong phạm vi địa lý: 104o11'09’’– 104o42'54 kinh độ Đông, 20o42'8 – 21o10'15 vĩ độ Bắc, phía Tây giáp lƣu vực suối Nậm Pàn, phía Đông giáp
huyện Vân Hồ, phía Nam giáp Lào.
75
Hình 3-1. Bản đồ vị trí lưu vực Nậm Sập
3.2.1.2. Đặc điểm địa hình
Nằm cách thành phố Sơn La 40,5km trên trục Quốc lộ 6 Hà Nội - Sơn La -
Điện Biên, lƣu vực Nậm Sập là một lƣu vực nằm sâu trong nội địa, có đặc điểm địa
hình phức tạp, bị chia cắt mạnh và độ dốc lớn, xen kẽ giữa những dãy núi.
Lƣu vực Nậm Sập có 2 hệ thống núi chính chạy qua: Hệ thống núi tả ngạn
sông Đà và hệ thống núi xen giữa sông Đà và sông Mã, hầu hết các dãy núi trong
tỉnh đều thấp dần theo hƣớng Tây Bắc - Đông Nam.
- Hệ thống núi tả ngạn sông Đà: Ranh giới giữa Sơn La và Yên Bái, bắt
nguồn từ Nậm Khan (Quỳnh Nhai) có độ cao 1.130m, chạy qua Mƣờng La, Bắc
Yên đến Phù Yên với các đỉnh cao từ 1.000m đến 2.500m, hình thành lƣu vực tả
ngạn sông Đà.
- Hệ thống núi xen giữa lƣu vực sông Đà và sông Mã: Bắt nguồn từ đỉnh Tà
Con (Thuận Châu) có độ cao 1.717m, chạy qua Mai Sơn, Yên Châu, Mộc Châu
76
gồm các đỉnh núi cao từ 1.000m đến 1.500m.
Hình 3-2. Bản đồ địa hình lưu vực Nậm Sập
3.2.1.3. Đặc điểm thổ nhưỡng
Lƣu vực Nậm Sập có 4 nhóm đất chính: Đất nâu vàng, đất xám feralit, đất
xám mùn trên núi và núi đá.
1. Đất nâu vàng
- Diện tích: 16.289,42ha, chiếm 15,01% diện tích lƣu vực.
- Vị trí phân bố: Phân bố chủ yếu tại huyện Mộc Châu và Yên Châu. Đất nâu
vàng chủ yếu đƣợc hình thành ở địa hình sƣờn dốc, phân bố ở địa hình bị chia cắt
mạnh nên thảm thực vật bị phá hủy. Phần diện tích đất tƣơng đối bằng phẳng.
2. Đất xám feralit
- Diện tích: 50.531,87ha, chiếm 46,57% diện tích lƣu vực.
- Vị trí phân bố: Phân bố chủ yếu tại huyện Yên Châu, Mộc Châu và một phần
thuộc huyện Bắc Yên, đất có độ phì cao.
77
3. Đất xám mùn trên núi
- Diện tích: 33.800,25ha, chiếm 31,15% diện tích lƣu vực.
- Vị trí phân bố: Phân bố chủ yếu tại huyện Mộc Châu và một phần thuộc
huyện Yên Châu. Thƣờng phân bố ở độ cao trên 900m. Khí hậu lạnh và ẩm hơn vùng dƣới, nhiệt độ bình quân năm khoảng 15 -200C. Địa hình cao, dốc, hiểm trở
nên xói mòn mạnh.
4. Núi đá
- Diện tích: 78.70,34ha, chiếm 7,2% diện tích lƣu vực.
- Vị trí phân bố: Phân bố chủ yếu tại khu vực giáp ranh giữa huyện Yên Châu
và Mộc Châu, một phần tại khu vực giáp giữa huyện Mộc Châu và Vân Hồ.
Hình 3-3. Bản đồ thổ nhưỡng lưu vực Nậm Sập
78
3.2.1.4. Đặc điểm thảm phủ thực vật và hiện trạng sử dụng đất
Lƣu vực Nậm Sập là lƣu vực có nhiều loại đất phù hợp với nhiều loại cây.
Diện tích rừng thuộc tiểu vùng là 27.573,5ha, chiếm 25,48% diện tích toàn tiểu
vùng, trong đó diện tích rừng tự nhiên giàu và trung bình là 8.847,06ha, diện tích
rừng tự nhiên nghèo là 11.752,06ha. Đất trồng cây công nghiệp chiếm 0,93%, đất
lúa màu chiếm 29,86%.
Bảng 3-1. Loại hình sử dụng đất tại lưu vực Nậm Sập
STT Loại rừng Diện tích (ha)
1 Đất trồng cây công nghiệp 1.005,56 So với diện tích lƣu vực (%) 0,93
2 Đất trống cỏ cây bụi 19.972,67 18,40
3 Đất trống cỏ cây gỗ rải rác 3.181,93 2,93
17,90 4 Đất trống có cỏ 19.434,16
0,51 5 Các loại đất khác 552,79
29,86 6 Lúa màu 32.408,59
0,78 7 Nƣơng rẫy 848,34
2,95 8 Núi đá 3.204,38
1,28 9 Rừng giàu 1.388,22
0,70 10 Rừng hỗn giao 756,87
1,05 11 Rừng lồ ô 1.135,94
12 Rừng nghèo 11.752,82 10,83
1,10 13 Rừng non có trữ lƣợng 1.195,12
1,57 14 Rừng non chƣa có trữ lƣợng 1.707,64
1,12 15 Rừng trồng 1.213,22
0,83 16 Rừng tre luồng 905,17
6,87 17 Rừng trung b́ ình 7.458,84
0,05 18 Rừng vầu 59,65
79
Hình 3-4. Bản đồ hiện trạng sử dụng đất lưu vực Nậm Sập năm 2015.
3.2.1.5. Đặc điểm khí hậu
Lƣu vực Nậm Sập có khí hậu nhiệt đới, mang đặc điểm khí hậu chung của vùng
Tây Bắc: Mùa đông lạnh, mùa hè nóng ẩm. Tuy nhiên, chế độ nhiệt, chế độ mƣa, số
giờ nắng có khác so với một số tiểu vùng.
Lƣợng bốc hơi trên lƣu vực Nậm Sập tại trạm Cò Nòi dao động từ 898 - 1487
mm/năm, tại trạm Mộc Châu từ 540 - 1.383mm/năm. Lƣợng bốc hơi trung bình trong
tháng tại trạm Cò Nòi dao động trong khoảng từ 74 - 124mm, tại trạm Mộc Châu là
45 - 115mm. Lƣợng bốc hơi lớn nhất vào tháng 4, nhỏ nhất vào tháng 7 (trạm Cò Nòi).
80
Tổng số giờ nắng trung bình các tháng tại trạm Cò Nòi từ 143 - 195 giờ/tháng,
tại trạm Mộc Châu từ 145 - 170 giờ/tháng, tháng có số giờ nắng ít nhất là tháng 1 và
tháng 2, nhiều nhất vào tháng 4, tháng 5, riêng cao nguyên Mộc Châu tháng 2 có số
giờ nắng lớn nhất so với các tháng còn lại trong năm. Tổng số giờ nắng trong năm
dao động từ 1.700 - 2.394 giờ/năm (trạm Cò Nòi), từ 1.354 - 2.022 giờ/năm (trạm
Mộc Châu).
Nhiệt độ trung bình năm tại trạm Cò Nòi dao động trong khoảng từ 19,70C - 22,20C, trạm Mộc Châu từ 17,00C - 20,30C, tháng 6 có nhiệt độ cao nhất tại trạm Cò Nòi là 25,80C, tại trạm Mộc Châu 23,40C vào tháng 7, nhiệt độ thay đổi lớn giữa
mùa đông và mùa hè, giữa ngày và đêm. Tuy nhiên, nhiệt độ có xu thế tăng trong
những năm gần đây.
Độ ẩm trung bình năm tại trạm Cò Nòi đạt khoảng từ 73% - 87%, tại trạm
Mộc Châu đạt khoảng từ 82% - 89%, mùa mƣa từ tháng 5 đến 10, độ ẩm tƣơng đối
cao, độ ẩm thấp nhất vào tháng 3, 4 và xuống đến 73% - 74%.
Trong năm đƣợc phân thành hai mùa rõ rệt: Mùa mƣa và mùa khô, mùa mƣa
kéo dài từ tháng 4 đến tháng 9 chiếm khoảng 75% - 80% tổng lƣợng mƣa cả năm,
tháng có lƣợng mƣa lớn nhất là tháng 7,8 đạt 250 - 280mm/tháng (trạm Cò Nòi) và
280 - 320mm/tháng (trạm Mộc Châu). Mùa khô kéo dài từ tháng 10 đến tháng 4
năm sau, lƣợng mƣa chỉ chiếm từ 20 - 25% tổng lƣợng mƣa năm, hai tháng có
lƣợng mƣa nhỏ nhất là tháng 12 và tháng 1.
Theo số liệu thống kê năm 2013 lƣợng mƣa bình quân năm của Sơn La là
111,4mm. Lƣợng mƣa trung bình tháng cao nhất tập trung vào tháng 7 với lƣợng
mƣa 335,7mm, tháng 12 có lƣợng mƣa thấp nhất với lƣợng mƣa trung bình 3,4mm.
(Niên giám thống kê tỉnh Sơn La năm 2013).
3.2.1.6. Đặc điểm khí tượng thủy văn
1. Hệ thống sông suối
81
Tiểu vùng Nậm Sập nằm trong lƣu vực của Nậm Sập Việt, gồm các phụ lƣu
chính nhƣ sau: Suối A Má, suối Ƣng, suối Môn, suối Co Păm, suối So Lung và suối
Vạt, Huổi Thƣơng, ngoài ra còn có rất nhiều các con suối nhỏ khác nhau, mật độ sông suối trên toàn lƣu vực khoảng 0,16km/km2. Mạng sông thƣa ở vùng đá vôi
(huyện Mộc Châu), dày hơn ở vùng Bắc Yên, Yên Châu.
Hình 3-5. Bản đồ mạng lưới sông ngòi lưu vực Nậm Sập
2. Mạng lưới trạm đo khí tượng thủy văn
+ Mạng lưới trạm đo mưa
82
Lƣu vực Nậm Sập có 5 trạm đo mƣa đƣợc xây dựng và hoạt động từ những
năm 60 của thế kỷ 20, hiện nay có 2 trạm đã ngừng hoạt động, chỉ còn trạm Mộc
Châu và Cò Nòi, Tà Làng đang hoạt động với chất lƣợng đo đạc đáng tin cậy.
+ Mạng lưới trạm Thủy Văn
Trên lƣu vực có 2 trạm thủy văn, nhƣng hiện nay chỉ còn trạm Thác Mộc
đang hoạt động.
Hình 3-6. Bản đồ lưới trạm quan trắc trên lưu vực Nậm Sập
a. Tổng lượng dòng chảy
83
Tài nguyên nƣớc mặt của lƣu vực Nậm Sập hàng năm vào khoảng 1,4 tỷ m3.
Tổng lƣợng dòng chảy trong 5 tháng mùa lũ chiếm khoảng 80% tổng lƣợng dòng
chảy năm, dòng chảy lớn nhất thƣờng tập trung vào tháng 8 hàng năm, các tháng
kiệt nhất thƣờng xảy ra vào tháng 3.
Dƣới tác động của các yếu tố khí hậu và mặt đệm, đặc biệt là mƣa và địa
hình, dòng chảy năm của sông suối cũng phân bố không đều trong lãnh thổ. Hình 3-
8 là bản đồ đƣờng đẳng trị mô đun dòng chảy năm trung bình thời kỳ 1963 -2012
trên lƣu vực Nậm Sập.
Hình 3-7. Bản đồ đẳng trị mô đun dòng chảy lưu vực Nậm Sập
84
b. Dòng chảy lũ
Mùa lũ trên lƣu vực Nậm Sập có thời gian xuất hiện trùng với mùa lũ trên các
sông, suối trên địa bàn tỉnh Sơn La, lƣu lƣợng dòng chảy bình quân tháng lớn nhất thƣờng rơi vào tháng 9, 10 và dao động trong khoảng từ 0,86 - 1,71m3/s.km2, lƣu lƣợng dòng chảy trung bình trong các tháng mùa lũ khoảng 0,96m3/s.km2.
Lƣu lƣợng đỉnh lũ dao động trong khoảng 340 - 365 m3/s, lũ lên nhanh, thời
gian lũ 2 - 3 ngày với thời gian đạt đỉnh 10 - 20 giờ. Biên độ lũ dao động từ 290 - 300m3/s.
c. Dòng chảy kiệt
Mùa khô tại lƣu vực Nậm Sập thƣờng bắt đầu từ tháng 10 đến tháng 4 năm
sau, khi mùa mƣa kết thúc lƣu lƣợng dòng chảy mặt giảm dần, lƣu lƣợng dòng chảy trung bình tháng thấp nhất thƣờng rơi vào tháng 3 với lƣu lƣợng nhỏ nhất 0,6m3/s. Đặc điểm tài nguyên nƣớc mặt trên lƣu vực Nậm Sập phụ thuộc nhiều vào sự điều
tiết dòng chảy trên sông Đà. Lƣu lƣợng dòng chảy tối thiểu trên suối Sập Việt là 1,32m3/s, suối Vạt là 0,54m3/s.
d. Dòng chảy bùn cát
Lƣợng bùn cát trên hệ thống sông, suối của lƣu vực Nậm Sập đƣợc hình
thành từ các tiểu lƣu vực do quá trình xói mòn và rửa trôi. Vào các tháng mùa lũ
hàm lƣợng bùn cát cao gấp nhiều lần so với các tháng mùa kiệt, thƣờng dao động trong khoảng 1000 - 5000g/m3. Tháng có hàm lƣợng bùn cát cao nhất thƣờng rơi vào tháng 8, 9 hàng năm dao động trong khoảng 2000 - 5000g/m3, các tháng có hàm lƣợng bùn cát thấp nhất thƣờng là tháng 1, 2 dao động trong khoảng 150 - 200g/m3.
3.2.1.7. Yêu cầu số liệu đầu vào mô hình
- Số liệu về khí tƣợng thủy văn của lƣu vực các năm 1962, 1973, 1980, 2010,
2011, 2012 và 2013 bao gồm các số liệu mƣa (giờ), nhiệt độ, bốc hơi và số liệu lƣu
lƣợng dòng chảy (giờ), lƣu lƣợng bùn cát (giờ) tại trạm Thác Mộc.
- Số liệu về tình hình sử dụng đất, về loại đất trên lƣu vực Nậm Sập tính đến
85
trạm Thác Mộc.
- Các bản đồ địa hình, bản đồ mạng lƣới sông suối, lƣới trạm khí tƣợng thủy
văn, bản đồ thổ nhƣỡng, bản đồ sử dụng đất.
Sử dụng công cụ GIS để phân chia lƣu vực thành các tiểu lƣu vực khác nhau.
Lƣu vực Nậm Sập đƣợc chia thành 72 tiểu lƣu vực. Với mỗi tiểu lƣu vực xác định đƣợc: Diện tích, độ dốc v.v... Chi tiết về đặc tính của từng tiểu lƣu vực đƣợc thể
hiện ở hình 3-8 và bảng 1 phần phụ lục.
Hình 3-8. Các tiểu lưu vực của Nậm Sập
3.2.1.8. Hiệu chỉnh, kiểm định mô hình
1. Hiệu chỉnh mô hình
Sau khi thiết lập đƣợc mô hình, tác giả chọn trận lũ tháng 9/1962 để hiệu
chỉnh bộ thông số của mô hình. Trong quá trình hiệu chỉnh luôn so sánh kết quả tính
toán lƣu lƣợng dòng chảy và lƣu lƣợng bùn cát với số liệu thực đo để hiệu chỉnh bộ
thông số của mô hình. Tiến hành hiệu chỉnh và so sánh tới khi nào đƣờng biểu diễn
86
thực đo và đƣờng biểu diễn tính toán tiệm cận với nhau, khi đó bộ thông số tìm
đƣợc là đạt và có thể sử dụng để tính toán lƣu lƣợng dòng chảy và lƣu lƣợng bùn
cát cho những năm không có số liệu thực đo. Kết quả hiệu chỉnh mô hình đƣợc thể
hiện (hình 3-9) dƣới dạng các biểu đồ so sánh kết quả tính toán và thực đo tại vị trí
các trạm thủy văn kiểm tra và chỉ số kiểm định NASH tƣơng ứng tại các trạm đó.
STT Bộ thông số mô hình Giá trị
4,8 e-05 1 Hệ số xói mòn tách rời liên rãnh (ail)
2 Hệ số khả năng xói mòn đất (ket) 1,22
2,15 e-04 3 Hệ số vận chuyển bùn cát liên rãnh (acl)
0.215 4 Hệ số xói mòn rãnh (Ci)
1.09 5 Hệ số vận chuyển bùn cát rãnh (ac)
Hình 3-9. Đường quá trình lưu lượng tính toán và thực đo
87
Hình 3-10. Đường quá trình hàm lượng bùn cát tính toán và thực đo
2. Kiểm định mô hình
Dùng bộ thông số của mô hình đã đƣợc hiệu chỉnh với trận lũ 1962 để kiểm định mô hình với trận lũ 1973. Kết quả mô phỏng đƣờng quá trình lƣu lƣợng dòng chảy,
hàm lƣợng bùn cát tính toán và thực đo tiệm cận nhau, đƣợc thể hiện ở hình sau:
Hình 3-11. Đường quá trình lưu lượng tính toán và thực đo
88
Hình 3-12. Đường quá trình hàm lượng bùn cát tính toán và thực đo
Kết quả tính toán hiệu chỉnh thông số của mô hình cho hệ số Nash đạt 81%, chênh lệch đỉnh lũ tính toán và thực đo là 13m3/s. Hệ số tƣơng quan giữa kết quả
tính toán và số liệu thực đo đạt 0.93.
Hình 3-13. Hàm lượng bùn cát trên lưu vực Thác Mộc
89
Do trạm Thác Mộc chỉ đo độ đục bùn cát đến năm 1980, do đó để mô phỏng
lƣu lƣợng bùn cát trên lƣu vực trong những năm gần đây, sử dụng bộ thông số mô
hình với số liệu mƣa giờ mô phỏng lƣu lƣợng bùn cát trên lƣu vực cho năm 2010,
2011, 2012 và 2013. Kết quả mô phỏng nhƣ hình 3-13.
Qua mô phỏng cho thấy lƣu lƣợng bùn cát trên lƣu vực phụ thuộc vào cƣờng
độ mƣa và tỷ lệ che phủ trên lƣu vực, độ dốc của lƣu vực v.v... Kết quả đƣợc chỉ
dẫn trong nội dung phân tích độ nhạy của các thông số của mô hình.
3.2.2. Lưu vực Phiêng Hiềng
3.2.2.1. Vị trí địa lý
Hình 3-14. Bản đồ vị trí lưu vực Phiêng Hiềng
90
Lƣu vực Phiêng Hiềng thuộc địa phận các huyện Bắc Yên, Phù Yên tỉnh Sơn La, có diện tích tự nhiên 431km2, chiếm 3,04% tổng diện tích toàn tỉnh. Nằm trong phạm vi địa lý: 104o24'41’ - 104o38'8’’ kinh độ Đông, 21o24'49’’- 21o05'38’’ vĩ độ
Bắc, phía Tây giáp lƣu vực suối Gạo huyện Bắc Yên, phía Bắc giáp với tỉnh Yên
Bái, phía Đông giáp với lƣu vực Suối Tốc huyện Phù Yên, phía Nam giáp với xã
Chiềng Sại huyện Bắc Yên.
3.2.2.2. Đặc điểm địa hình
Lƣu vực Phiêng Hiềng là một lƣu vực nằm sâu trong nội địa, có đặc điểm địa
hình phức tạp, bị chia cắt mạnh và độ dốc lớn, 100% diện tích tự nhiên thuộc lƣu
vực sông Đà, xen kẽ giữa những dãy núi.
Hình 3-15. Bản đồ địa hình lưu vực Phiêng Hiềng
91
Lƣu vực Phiêng Hiềng có 2 hệ thống núi chính chạy qua: Hệ thống núi tả
ngạn sông Đà và hệ thống núi thấp xen kẽ, hầu hết các dãy núi trong tiểu vùng đều
thấp dần theo hƣớng Tây Bắc - Đông Nam.
- Hệ thống núi tả ngạn sông Đà: Ranh giới giữa Sơn La và Yên Bái, bắt
nguồn từ Nậm Khan (Quỳnh Nhai) có độ cao 1.130m, chạy qua Mƣờng La, Bắc
Yên đến Phù Yên với các đỉnh cao từ 1.000m đến 2.500m, hình thành lƣu vực tả
ngạn sông Đà.
- Hệ thống núi thấp xen kẽ: Tập trung chủ yếu tại các xã Hồng Ngài, thị trấn
Bắc Yên, xã Suối Bau, xã Đá Đỏ, các đỉnh núi cao từ 1.000m đến 1.500m.
3.2.2.3. Đặc điểm thổ nhƣỡng
Lƣu vực Phiêng Hiềng có 3 nhóm đất chính: Đất xám feralit, đất nâu đỏ, đất
xám mùn trên núi.
1. Đất xám feralit
Diện tích: 14.797,85ha, chiếm 34,38% diện tích lƣu vực.
Vị trí phân bố: Phân bố chủ yếu tại xã Hồng Ngài, thị trấn Bắc Yên, Quang
Huy, Suối Tọ. Đất có độ phì cao, phù hợp cho trồng cây ăn quả tại vùng đất bằng
phẳng và trồng rừng tại những khu vực dốc.
2. Đất nâu đỏ
Diện tích: 7.180,03ha, chiếm 16,68% diện tích lƣu vực.
Vị trí phân bố: Phân bố chủ yếu tại xã Suối Bau, Sập Xa (huyện Phù Yên). Đất
nâu đỏ là loại đất tốt về mặt tính chất đất cũng nhƣ yếu tố địa hình tƣơng đối thuận
lợi cho việc canh tác, sản xuất nông lâm nghiệp.
3. Đất xám mùn trên núi
Diện tích: 21.058,30ha, chiếm 48,93% diện tích lƣu vực.
92
Vị trí phân bố: Phân bố chủ yếu tại các xã Phiêng Ban, Tà Xùa, Suối Tọ
(huyện Bắc Yên). Thƣờng phân bố ở độ cao trên 900m. Khí hậu lạnh và ẩm hơn vùng dƣới, nhiệt độ bình quân năm khoảng 15 - 200C. Địa hình cao, dốc, hiểm trở
nên xói mòn mạnh.
Hình 3-16. Bản đồ thổ nhưỡng lưu vực Phiêng Hiềng
93
3.2.2.4. Đặc điểm thảm phủ thực vật và hiện trạng sử dụng đất
Lƣu vực Phiêng Hiềng có nhiều loại đất phù hợp với nhiều loại cây, diện tích
đất trống chiếm 0,26%, đất trống có cỏ chiếm 29,60%, diện tích rừng nghèo chiếm
9,07%, diện tích rừng non có trữ lƣợng chiếm 1,57% v.v...
Bảng 3-2. Loại hình sử dụng đất tại lưu vực Phiêng Hiềng
STT Loại rừng
Diện tích (ha) 112.895,40 So với diện tích lƣu vực (%) 0,26 1 §Êt trèng
2 6.597,02 15,33 §Êt trèng cã c©y bôi
3 3.546,695 8,24 §Êt trèng cã c©y gç r¶i r¸c
4 12.738,70 29,60 §Êt trèng cã cá
5 7.438,01 17,28 Lóa mµu
6 53.821,67 0,13 C¸c lo¹i ®Êt kh¸c
7 2.385,28 5,54 Rõng giµu
8 158.603,10 0,37 Rõng hçn giao
9 3.902,69 9,07 Rõng nghÌo
2.675,11 6,22 10 Rõng non cã tr÷ lƣợng
1.850,81 4,30 11 Rõng non chƣa cã tr÷ lƣợng
30.203,17 0,07 12 Rõng trång
615.286,10 1,43 13 Rõng tre luång
616.941,80 1,43 14 Rõng trung b×nh
15 Núi đá 42.494,83 0,10
125.901,90 0,29 16 Rõng vÇu
140.759,50 0,33 17 Rõng lå «
0,01 18 Đất trồng cây công nghiệp 57.581,76
94
Hình 3-17. Bản đồ hiện trạng sử dụng đất lưu vực Phiêng Hiềng năm 2015
3.2.2.5. Đặc điểm khí hậu
Lƣu vực Phiêng Hiềng có khí hậu nhiệt đới, mang đặc điểm khí hậu chung của
vùng Tây Bắc: Mùa đông lạnh, mùa hè nóng ẩm. Tuy nhiên, chế độ nhiệt, chế độ
mƣa, số giờ nắng có khác so với một số tiểu vùng. Do lƣu vực không có trạm đo
nên sử dụng số liệu của các trạm Phù Yên, Bắc Yên là 2 trạm thuộc vùng phụ cận
tiểu lƣu vực.
Lƣợng bốc hơi tại trạm Phù Yên dao động từ 789 - 1.123mm/năm, tại trạm Bắc
Yên từ 770 - 1.291mm/năm. Lƣợng bốc hơi trung bình trong tháng tại trạm Phù Yên
dao động trong khoảng 61 - 106mm, tại trạm Bắc Yên là 67 - 133mm. Lƣợng bốc hơi
95
lớn nhất vào tháng 5, nhỏ nhất vào tháng 8 (trạm Bắc Yên) và lớn nhất vào tháng 6,
nhỏ nhất vào tháng 1 (trạm Phù Yên).
Tổng số giờ nắng trung bình các tháng tại trạm Phù Yên 90 - 170 giờ/tháng,
tại trạm Bắc Yên 119 - 173 giờ/tháng, tháng 1, 2 có số giờ nắng ít nhất, tháng 5, 6
có số giờ nắng nhiều nhất (trạm Phù Yên), tháng 1 có số giờ nắng ít nhất, tháng 4 có
số giờ nắng nhiều nhất (trạm Bắc Yên). Tổng số giờ nắng trong năm dao động từ
1319 - 2011giờ/năm (trạm Phù Yên), từ 1368 - 1922 giờ/năm (trạm Bắc Yên).
Nhiệt độ trung bình năm tại trạm Bắc Yên dao động trong khoảng từ 19,70C - 21,70C; tại trạm Phù Yên từ 22,60C - 24,60C; tháng 6 có nhiệt độ cao nhất là 28,60C (trạm Phù Yên), 26,10C ( trạm Bắc Yên). Tuy nhiên, nhiệt độ có xu thế tăng trong
những năm gần đây.
Độ ẩm trung bình năm tại trạm Phù Yên đạt khoảng 77% - 83%, trạm Bắc Yên
đạt khoảng 78% - 84%, độ ẩm tƣơng đối cao, độ ẩm thấp nhất vào tháng 3, 4 xuống
đến 78% - 79%.
Trong năm đƣợc phân thành hai mùa rõ rệt: Mùa mƣa và mùa khô, mùa mƣa
kéo dài từ tháng 4 đến tháng 9 chiếm khoảng 75% - 80% tổng lƣợng mƣa cả năm,
tháng 7, 8 có lƣợng mƣa lớn nhất đạt 278 - 279mm/tháng (trạm Phù Yên) và 250 -
258mm/tháng (trạm Bắc Yên). Mùa khô kéo dài từ tháng 10 đến tháng 4 năm sau,
lƣợng mƣa chỉ chiếm từ 23% đến 28% tổng lƣợng mƣa năm, tháng 12 có lƣợng
mƣa nhỏ nhất.
Theo số liệu thống kê năm 2013 lƣợng mƣa bình quân năm của Sơn La là
111,4mm. Lƣợng mƣa trung bình tháng cao nhất tập trung vào tháng 7 với lƣợng
mƣa 335,7mm, tháng 12 có lƣợng mƣa thấp nhất với lƣợng mƣa trung bình 3,4mm.
(Niên giám thống kê tỉnh Sơn La năm 2013).
3.2.2.6. Đặc điểm khí tƣợng thủy văn
1. Hệ thống sông suối
Lƣu vực Phiêng Hiềng có các phụ lƣu chính nhƣ sau: Suối Háng Đồng, suối
96
Ban, suối Bau, suối Làng, suối Tọ ngoài ra còn có rất nhiều các con suối nhỏ khác
nhau, mật độ sông suối trên toàn lƣu vực khoảng 0,16km/km2. Mạng lƣới sông thƣa
ở vùng đá vôi (huyện Phù Yên), dày hơn ở vùng Bắc Yên.
Hình 3-18. Bản đồ mạng lưới sông ngòi lưu vực Phiêng Hiềng
2. Mạng lưới trạm đo khí tượng thủy văn
+ Mạng lƣới trạm khí tƣợng
Do lƣu vực Phiêng Hiềng nhỏ nên không bố trí trạm đo mƣa, vùng phụ cận
lƣu vực có trạm đo mƣa Bắc Yên và Phù Yên với chất lƣợng đo đạc đáng tin cậy.
97
+ Mạng lƣới trạm Thủy Văn
Trên lƣu vực có trạm thủy văn Phiêng Hiềng đƣợc xây dựng từ năm 1960,
hiện nay đã ngừng hoạt động.
Hình 3-19. Bản đồ lưới trạm quan trắc trên lưu vực Phiêng Hiềng
a. Tổng lượng dòng chảy
Tài nguyên nƣớc mặt của lƣu vực Phiêng Hiềng hàng năm vào khoảng 0,57 tỷ m3. Tổng lƣợng dòng chảy trong 5 tháng mùa lũ chiếm khoảng 78% tổng lƣợng
98
dòng chảy năm, dòng chảy lớn nhất thƣờng tập trung vào tháng 8 hàng năm, các
tháng kiệt nhất thƣờng xảy ra vào tháng 3.
Hình 3-20. Bản đồ đẳng trị mô đun dòng chảy lưu vực Phiêng Hiềng
Dƣới tác động của các yếu tố khí hậu và mặt đệm, đặc biệt là mƣa và địa hình,
dòng chảy năm của sông suối cũng phân bố không đều trong lãnh thổ. Hình 3-18 là
99
bản đồ đƣờng đẳng trị mô đun dòng chảy năm trung bình thời kỳ 1963 - 2012 trên
lƣu vực Phiêng Hiềng.
b. Dòng chảy lũ
Mùa lũ trên lƣu vực Phiêng Hiềng có thời gian xuất hiện trùng với mùa lũ trên
các sông, suối trên địa bàn tỉnh Sơn La, mô đun dòng chảy bình quân tháng lớn nhất thƣờng rơi vào tháng 9, 10 và dao động trong khoảng 0,76 - 1,36m3/s.km2, mô đun dòng chảy trung bình trong các tháng mùa lũ khoảng 0,86m3/s.km2.
Lƣu lƣợng đỉnh lũ dao động trong khoảng 320 - 345m3/s, lũ lên nhanh, thời
gian lũ từ 2 - 3 ngày với thời gian đạt đỉnh từ 10 - 20 giờ. Biên độ lũ dao động từ 280 - 300m3/s.
c. Dòng chảy kiệt
Mùa khô tại lƣu vực Phiêng Hiềng thƣờng bắt đầu từ tháng 10 đến tháng 4
năm sau, khi mùa mƣa kết thúc lƣu lƣợng dòng chảy mặt giảm dần, mô đun dòng
chảy trung bình tháng thấp nhất thƣờng rơi vào tháng 3. Lƣu lƣợng dòng chảy tối thiểu trên Phiêng Hiềng là 1,94m3/s, suối Hàng Đồng là 0,58m3/s, suối Ban là 0,25m3/s
d. Dòng chảy bùn cát
Lƣợng bùn cát trên hệ thống sông, suối của lƣu vực Phiêng Hiềng đƣợc hình
thành từ các tiểu lƣu vực do quá trình xói mòn và rửa trôi. Vào các tháng mùa lũ
hàm lƣợng bùn cát cao gấp nhiều lần so với các tháng mùa kiệt. Tháng 8, 9 có hàm lƣợng bùn cát cao nhất dao động trong khoảng 2700 - 4900g/m3, tháng 1, 2 có hàm lƣợng bùn cát thấp nhất dao động trong khoảng 160 - 210g/m3.
3.2.2.7. Yêu cầu số liệu đầu vào mô hình
- Số liệu về khí tƣợng thủy văn của lƣu vực các năm 1973, 1976 bao gồm các
số liệu mƣa (giờ), nhiệt độ, bốc hơi tại trạm Phù Yên và số liệu lƣu lƣợng dòng
100
chảy (giờ), lƣu lƣợng bùn cát (giờ) tại trạm Phiêng Hiềng.
- Số liệu về tình hình sử dụng đất, về loại đất trên lƣu vực Phiêng Hiềng.
- Các bản đồ địa hình, bản đồ mạng lƣới sông suối, lƣới trạm khí tƣợng thủy
văn, bản đồ thổ nhƣỡng, bản đồ sử dụng đất.
Hình 3-21. Các tiểu lưu vực của Phiêng Hiềng
101
- Sử dụng công cụ GIS để phân chia lƣu vực thành các tiểu lƣu vực khác nhau.
Lƣu vực Phiêng Hiềng đƣợc chia thành 94 tiểu lƣu vực. Với mỗi tiểu lƣu vực xác
định đƣợc: Diện tích, độ dốc v.v... Chi tiết về đặc tính của từng tiểu lƣu vực đƣợc
thể hiện ở hình 3-21 và bảng 2 phần phụ lục.
3.2.2.8. Hiệu chỉnh, kiểm định mô hình
1. Hiệu chỉnh mô hình
Sau khi thiết lập đƣợc mô hình, tác giả chọn trận lũ tháng 8/1973 để hiệu
chỉnh bộ thông số của mô hình. Trong quá trình hiệu chỉnh luôn so sánh kết quả tính
toán lƣu lƣợng dòng chảy và lƣu lƣợng bùn cát với số liệu thực đo để hiệu chỉnh bộ
thông số của mô hình. Tiến hành hiệu chỉnh và so sánh tới khi nào đƣờng biểu diễn
thực đo và đƣờng biểu diễn tính toán phù hợp cả về dạng đƣờng và đỉnh thì khi đó
bộ thông số tìm đƣợc là đạt và có thể sử dụng để tính toán lƣu lƣợng dòng chảy và
lƣu lƣợng bùn cát cho những năm không có số liệu thực đo. Kết quả hiệu chỉnh mô
hình đƣợc thể hiện ở hình sau và bộ thông số của mô hình:
5,3 e-05
STT Bộ thông số mô hình Giá trị
1,12
1 Hệ số xói mòn tách rời liên rãnh (ail)
2,5 e-04
2 Hệ số khả năng xói mòn đất (ket)
0.321
3 Hệ số vận chuyển bùn cát liên rãnh (acl)
1.05
4 Hệ số xói mòn rãnh (Ci)
5 Hệ số vận chuyển bùn cát rãnh (ac)
102
Hình 3-22. Đường quá trình lưu lượng tính toán và thực đo
Hình 3-23. Đường quá trình hàm lượng bùn cát tính toán và thực đo
103
2. Kiểm định mô hình
Dùng bộ thông số của mô hình đã đƣợc hiệu chỉnh với trận lũ 1973 để kiểm
định mô hình với trận lũ 1976. Kết quả mô phỏng đƣờng quá trình lƣu lƣợng dòng
chảy, hàm lƣợng bùn cát tính toán và thực đo tiệm cận nhau, đƣợc thể hiện ở hình:
Hình 3-24. Đường quá trình lưu lượng tính toán và thực đo
Hình 3-25. Đường quá trình hàm lượng bùn cát tính toán và thực đo
104
Kết quả tính toán hiệu chỉnh thông số của mô hình tính toán lƣu lƣợng dòng
chảy lũ cho lƣu vực Phiêng Hiềng cho hệ số Nash đạt 83%, chênh lệch đỉnh lũ tính
toán và thực đo là 16 m3/s.
3.2.3. Xây dựng phƣơng trình tƣơng quan
3.2.3.1. Phân tích độ nhạy của các thông số mô hình
1. Phương pháp phân tích độ nhạy các thông số
Phân tích độ nhạy là sự nghiên cứu mối quan hệ giữa thông tin vào và ra của
mô hình. Độ nhạy có thể tính toán bằng nhiều phƣơng pháp hay phân tích định tính
hoặc định lƣợng. Có thể kể đến một số công trình phân tích độ nhạy nhƣ [18, 65].
Để phân tích độ nhạy các thông số trong mô hình toán, ngƣời ta đƣa ra khái
niệm véc tơ ứng với mỗi thông số (ví dụ véc tơ x nhận các giá trị trong khoảng a
b, ta viết x[a,b]. Số lƣợng các trị nhân tố thuộc [a,b] nhiều hay ít tùy thuộc yêu cầu
độ chính xác của phƣơng pháp.
Giả sử có véc tơ x[a,b], dùng phƣơng pháp Monte - Carlo [66] để lấy ngẫu
nhiên các giá trị nhân tố x k theo một hàm mật độ xác suất nhất định. Ứng với mỗi
giá trị nhân tố x k, ngƣời ta sẽ nhận đƣợc một giá trị kết quả tính toán của mô
hình, các giá trị của các hàm chỉ tiêu (ví dụ hàm Nash) đƣợc tính toán để so sánh
giữa kết quả tính toán và giá trị quan trắc. Sắp xếp theo chiều giảm dần các giá trị
của hàm chỉ tiêu, các giá trị x k đƣợc thay đổi theo. Các giá trị phân phối xác suất
của x k đƣợc chia thành hai phần: Phần thứ nhất là tập hợp các giá trị “tốt” (A),
phần thứ hai là tập hợp giá trị “không tốt” (B). Độ nhạy của thông số đang xét
đƣợc đánh giá theo khoảng cách lớn nhất theo phƣơng đứng của hai đƣờng
phân phối xác suất. Giá trị tính toán của sẽ so sánh với giá trị tiêu chuẩn, khi
105
tính toán ≥ tiêu chuẩn thì thông số x đang xét sẽ ảnh hƣởng nhiều đến
kết quả tính toán của mô hình, ngƣợc lại sẽ không ảnh hƣởng nhiều. Giá trị tiêu
chuẩn đƣợc tính theo công thức sau:
(3-1)
trong đó là hệ số phụ thuộc vào mức độ chính xác để lấy giá trị chỉ số lƣợng tính
toán “tốt” hay “không tốt”.
Trên cơ sở so sánh giữa kết quả tính toán của mô hình và giá trị quan trắc thực
tế, các giá trị thống kê đƣợc gán cho mỗi thông số. Thông qua các giá trị thống kê
đƣợc tính thông qua các hàm tiêu chuẩn nhƣ hàm Nash. Với tập hợp các giá trị (giá
trị khả dĩ) hàm tiêu chuẩn của các thông số, chúng ta có thể phân tích độ nhạy của
mỗi thông số trong cùng một mô hình tính toán.
2. Kết quả phân tích độ nhạy các thông số
Có nhiều thông số tham gia vào quá trình hình thành vận chuyển bùn cát trên
lƣu vực vừa và nhỏ mà đã đƣợc mô hình hóa thông qua mô hình mô phỏng vận
chuyển bùn cát trên lƣu vực. Phƣơng pháp phân tích độ nhạy có thể chỉ ra đƣợc
thông số nào có vai trò quan trọng trong mô hình, thông số nào không quan trọng.
Luận án phân tích độ nhạy tổ hợp các thông số: Cƣờng độ mƣa, hệ số xói mòn
bong tách liên rãnh, độ dốc lƣu vực, mật độ rãnh, độ che phủ đất và khả năng xói
mòn của đất. Kết quả phân tích độ nhạy của các thông số thể hiện nhƣ trên hình vẽ
sau:
106
Hình 3-26. Đường biểu diễn kết quả khi thay đổi thông số mô hình
Qua phân tích kết quả mô phỏng độ nhạy các thông số mô hình cho thấy,
cƣờng độ mƣa ảnh hƣởng lớn nhất đến sự vận chuyển bùn cát trên lƣu vực. Cƣờng
độ mƣa càng lớn thì mức vận chuyển bùn cát trên lƣu vực càng cao. Tiếp theo là
107
thông số độ dốc lƣu vực. Nếu độ dốc lƣu vực tăng lên 30% thì lƣợng bùn cát trên
lƣu vực tăng lên khoảng 12%. Nhƣ vậy độ dốc càng lớn thì xói mòn mặt càng lớn
và ngƣợc lại. Độ dốc càng lớn thì yêu cầu đối với cấu trúc thảm thực vật rừng
phòng hộ càng cao. Nếu giảm độ tàn che từ 0,7 - 0,8 xuống mức 0,3 - 0,4 thì xói
mòn đất sẽ tăng lên 42,2% và dòng chảy mặt tăng 30,4%. Tác giả đã thay đổi tỷ lệ
che phủ của lƣu vực, từ đó chỉ ra độ che phủ càng cao thì bùn cát càng giảm, khả
năng giữ nƣớc và đất càng tốt.
3.2.3.2. Phân tích tƣơng quan giữa xói mòn liên rãnh trên lƣu vực với độ dốc và
cƣờng độ mƣa
Lƣợng đất bị xói mòn liên rãnh trên lƣu vực là hàm số của cƣờng độ mƣa và
độ dốc lƣu vực. Thử nghiệm tính toán với các cƣờng độ mƣa khác nhau trên lƣu vực
Phiềng Hiềng và lƣu vực Nậm Sập cho ta kết quả tính toán xói mòn liên rãnh trên
lƣu vực.
b Ei= a.i2.Io
(3-2)
Phƣơng trình cơ bản cho sự bong tách liên rãnh có dạng tổng quát nhƣ sau:
trong đó:
a, b : Hệ số trong phƣơng trình tƣơng quan
i : Cƣờng độ mƣa (mm/h)
Io : Độ dốc lƣu vực (%).
Kết quả phân tích tƣơng quan giữa xói mòn liên rãnh trên lƣu vực với độ dốc
và cƣờng độ mƣa đƣợc chỉ ra trong hình vẽ sau:
108
Hình 3-27. Tương quan giữa xói mòn liên rãnh với độ dốc và cường độ mưa
109
Từ kết quả tính tác giả toán xây dựng đƣợc phƣơng trình cơ bản cho sự bong
tách liên rãnh trên lƣu vực vừa và nhỏ thuộc tỉnh Sơn La dựa trên số liệu tính toán
nhƣ sau:
(3-3)
trong đó:
i : Cƣờng độ mƣa (mm/h)
I0 : Độ dốc lƣu vực (%).
Nhƣ vậy với kết quả trên, trên cơ sở bản đồ địa hình lƣu vực, xác định đƣợc độ
dốc lƣu vực tại các tiểu lƣu vực, khi có cƣờng độ mƣa, ta sẽ xác định đƣợc lƣợng
đất bị xói mòn liên rãnh trên lƣu vực, từ đó đƣa ra các giải pháp để chống xói mòn,
giảm thiểu vận chuyển bùn cát trên lƣu vực.
3.2.3.3. Phân tích tƣơng quan giữa xói mòn rãnh trên lƣu vực với độ dốc và cƣờng
độ mƣa
Lƣợng bùn cát xói mòn rãnh trên lƣu vực là hàm số của cƣờng độ mƣa và độ
dốc lƣu vực. Thử nghiệm tính toán với các cƣờng độ mƣa khác nhau trên lƣu vực
cho ta kết quả tính toán xói mòn rãnh trên lƣu vực.
d
Phƣơng trình cơ bản cho xói mòn rãnh có dạng tổng quát nhƣ sau:
(3-4) Er = c.i2.Io
trong đó:
c, d : Hệ số trong phƣơng trình tƣơng quan
i : Cƣờng độ mƣa (mm/h)
I0 : Độ dốc lƣu vực (%).
Kết quả phân tích tƣơng quan giữa xói mòn rãnh trên lƣu vực với độ dốc và
cƣờng độ mƣa đƣợc chỉ ra trong hình vẽ sau:
110
Hình 3-28. Tương quan giữa xói mòn rãnh với độ dốc và cường độ mưa
Từ kết quả tính toán tác giả xây dựng đƣợc phƣơng trình cơ bản xói mòn rãnh
trên lƣu vực vừa và nhỏ thuộc tỉnh Sơn La dựa trên số liệu tính toán nhƣ sau:
(3-5)
111
Từ bản đồ địa hình lƣu vực, xác định đƣợc độ dốc lƣu vực tại các tiểu lƣu vực,
khi có cƣờng độ mƣa, ta sẽ xác định đƣợc lƣợng đất bị xói mòn rãnh trên lƣu vực
dựa vào phƣơng trình 3-5, từ đó đƣa ra các giải pháp để chống xói mòn, giảm thiểu
vận chuyển bùn cát trên lƣu vực.
Kết luận chƣơng 3
Tác giả chọn hai lƣu vực để kiểm tra độ tin cậy của mô hình, kết quả cho thấy
các đƣờng biểu diễn dòng chảy và đƣờng biểu diễn vận chuyển bùn cát trên hai lƣu
vực giữa tính toán và thực đo tiệm cận nhau. Ngoài ra, tác giả xây dựng đƣợc
phƣơng trình tƣơng quan giữa xói mòn liên rãnh với độ dốc và cƣờng độ mƣa,
phƣơng trình tƣơng quan giữa xói mòn rãnh với độ dốc và cƣờng độ mƣa dựa trên
quá trình phân tích độ nhạy của các thông số.
112
KẾT LUẬN VÀ KIẾN NGHỊ
Các kết quả nghiên cứu của luận án trong việc nghiên cứu cơ sở khoa học xây
dựng mô hình toán mô phỏng quá trình vận chuyển bùn cát trên lƣu vực vừa và nhỏ,
áp dụng cho các lƣu vực miền núi phía Bắc thuộc tỉnh Sơn La trên quan điểm tiếp
cận quá trình vật lý của quá trình vận chuyển bùn cát với các phƣơng pháp và công
cụ tính toán tiên tiến, đặc biệt là việc giải phƣơng trình toán học, kết hợp với việc lập trình bằng ngôn ngữ C++ đã rút ra một số kết luận nhƣ sau:
- Trên cơ sở tổng quan các mô hình mô phỏng quá trình vận chuyển bùn cát
trên lƣu vực vừa và nhỏ đang đƣợc áp dụng trên thế giới và ở Việt Nam, trong điều kiện các lƣu vực có địa hình dốc, việc lựa chọn ngôn ngữ C++ làm công cụ để mô
phỏng quá trình vận chuyển bùn cát trên lƣu vực vừa và nhỏ là hợp lý, thích hợp
trong điều kiện hiện nay, xây dựng công cụ
Để mô phỏng và phát triển công nghệ tính toán dòng chảy bùn cát trên lƣu
vực vừa và nhỏ là có cơ sở khoa học.
- Sơn La là tỉnh miền núi phía Bắc có địa hình dốc nên nguy cơ xảy ra xói
mòn, thoái hóa đất do vận chuyển bùn cát bề mặt lƣu vực là lớn, các yếu tố ảnh
hƣởng đến việc vận chuyển bùn cát trên lƣu vực gồm mƣa, thảm phủ thực vật, độ
dốc, loại đất và các hoạt động khai thác của con ngƣời. Các yếu tố tạo nguồn vật
chất cho việc vận chuyển bùn cát bề mặt lƣu vực là do thổ nhƣỡng.
- Nghiên cứu lý luận và thực nghiệm số các sơ đồ giải bằng phƣơng pháp sai
phân Lax - Friedrichs về quy mô không gian và thời gian nhằm nâng cao độ chính
xác và độ ổn định của mô hình tính toán, có độ ổn định cao nhất.
- Việc vận dụng mô hình tính toán mô phỏng vận chuyển bùn cát trên lƣu vực
cho thấy mƣa là yếu tố chính gây nên vận chuyển bùn cát trên lƣu vực. Độ dốc là
nhân tố quan trọng ảnh hƣởng đến xói mòn và dòng chảy mặt. Độ dốc càng lớn thì
113
xói mòn mặt càng lớn và ngƣợc lại. Nếu độ dốc lƣu vực tăng lên 30% thì lƣợng bùn
cát bị cuốn trôi trên lƣu vực tăng lên khoảng 12%.
- Từ lƣợng mƣa và độ dốc trên lƣu vực, thông qua phƣơng trình tƣơng quan
giữa xói mòn rãnh, liên rãnh với lƣợng mƣa và độ dốc, có thể tính toán đƣợc lƣợng
xói mòn rãnh và liên rãnh, từ đó tính đƣợc lƣợng bùn cát vận chuyển đến cửa ra của
lƣu vực.
Những đóng góp mới của luận án:
- Xây dựng đƣợc mô hình toán mô phỏng xói mòn và vận chuyển bùn cát trên lƣu
vực vừa và nhỏ, với thuật toán sơ đồ sai phân Lax - Friedrichs có thêm trọng số thời
gian, không gian để giải phƣơng trình dòng chảy và phƣơng trình vận chuyển bùn
cát trên lƣu vực.
- Xây dựng đƣợc phƣơng trình tƣơng quan giữa tính toán xói mòn liên rãnh,
xói mòn rãnh trên lƣu vực nghiên cứu: Lƣu vực Nậm Sập, lƣu vực Phiêng Hiềng
của tỉnh Sơn La, từ đó có thể dự báo lƣợng bùn cát bị xói mòn và vận chuyển trên
lƣu vực theo cƣờng độ mƣa.
Kiến nghị:
Việc tính toán mô phỏng vận chuyển bùn cát trong sông mới tính dựa trên cân
bằng bùn cát, tính toán cho hạt bùn cát đồng nhất. Vì vậy trong thời gian tới cần tiếp
tục nghiên cứu tính toán cho xói mòn đáy và xói mòn bờ.
Việc thử nghiệm mô hình mới chỉ trên hai lƣu vực nghiên cứu thuộc tỉnh Sơn
La, cần tiếp tục hoàn thiện mô hình, nghiên cứu thử nghiệm cho các lƣu vực khác.
Cần kết nối công nghệ GIS để thiết lập mô hình trên cơ sở mô phỏng địa hình
chi tiết của lƣu vực, nhƣ vậy sẽ tăng mức độ chính xác của mô hình.
114
DANH MỤC CÁC CÔNG TRÌNH CỦA TÁC GIẢ ĐÃ CÔNG BỐ
1. Đào Tấn Quy và Phạm Thị Hƣơng Lan. “Xây dựng bản đồ hiểm họa trƣợt
lở đất tỉnh Sơn La”. Tạp chí khoa học kỹ thuật Thủy Lợi và Môi trường. Số 51
T12/2015).
2. Đào Tấn Quy. “Xây dựng mô hình toán mô phỏng quá trình vận chuyển bùn
cát trên lƣu vực vừa và nhỏ. (Áp dụng cho lƣu vực suối Sập tỉnh Sơn La)”. Tạp chí
Khí tượng thủy văn. Số 657 T9/2015.
3. Đào Tấn Quy. “Xây dựng mô đun tính toán dòng chảy trong mô hình mô
phỏng quá trình vận chuyển bùn cát trên lƣu vực vừa và nhỏ (Áp dụng cho lƣu vực
suối Sập tỉnh Sơn La)”. Tạp chí Khoa học kỹ thuật Thủy Lợi và Môi trường. Số 48
T3/2015.
4. Đào Tấn Quy, Trần Kim Châu, Phạm Thị Hƣơng Lan và Nguyễn Thế Toàn.
“Ứng dụng hệ thông tin địa lý và quá trình phân tích cấp bậc để tiến hành xây dựng
bản đồ hiểm họa trƣợt lở đất tỉnh Sơn La”. Hội nghị khoa học thường niên ĐHTL
năm 2014. Trang 474- 476.
5. Đào Tấn Quy và Phạm Thị Hƣơng Lan. “Nghiên cứu cơ sở lý thuyết xây
dựng mô hình toán mô phỏng vận chuyển bùn cát trên lƣu vực vừa và nhỏ”. Hội
nghị khoa học thường niên ĐHTL năm 2013. Trang 175- 177.
6. Đào Tấn Quy và Phạm Thị Hƣơng Lan. “Nghiên cứu phân tích lựu chọn
ngôn ngữ phần mềm mô hình toán mô phỏng dòng chảy và vận chuyển bùn cát trên
lƣu vực vừa và nhỏ”. Tạp chí Khoa học kỹ thuật Thủy Lợi và Môi trường. Số 40
T3/2013.
115
TÀI LIỆU THAM KHẢO
[1] Đinh Văn Hùng., Ứng dụng Viễn thám và GIS đánh giá xói mòn đất khu vực
Yên Châu tỉnh Sơn La. Luận văn Thạc sĩ Khoa học, ĐHQG Hà Nội, 2009.
[2] Ellison W. D., Studies of raindrop erosion. Agricultural Engineering 25, p
181-182, 1944.
[3] FAO., The digital soil map of the world and derived soil properties. CDROM
Version 3.5. Rome, 1994.
[4] Hudson., Soil Conservation. London, Batsford, 1979.
[5] Cao Đăng Dƣ., Đánh giá tốc độ xói mòn lưu vực sông Hồng theo lượng bùn
cát trong sông. Báo cáo hội thảo khoa học, 2000.
[6] Nguyễn Tử Xiêm và Thái Phiên., Cơ sở khoa học kỹ thuật về canh tá c đất dốc .
Tài liệu hội thảo khoa ho ̣c “Sƣ̉ du ̣ng đất và bảo vệ rừng”, 1996.
[7] Bengt Carlsson., An introduction to sedimentation theory in wastewater
treatment. System and Control Group Uppsala University, 1998.
[8] Nguyễn Trọng Hà., Xác định các yếu tố gây xói mòn và khả năng dự báo xói
mòn trên đất dốc. Luận án Phó Tiến sĩ KH-KT, Trƣờng Ðại học Thủy lợi. Hà
Nội, 1996.
[9] Phạm Ngọc Dũng., Nghiên cứu một số biện pháp chống xói mòn trên đất đỏ
bazan trồng chè vùng Tây nguyên và xác định giá trị của các yếu tố gây xói
mòn đất theo mô hình Wischmeier W. H and Smith D.D. Luận án Phó tiến sĩ
khoa học Nông nghiệp. Hà Nội, 1991.
[10] Nguyễn Văn Dũng, Nguyễn Văn Dũng, Nguyễn Đình Kỳ., Đánh giá định
lượng xói mòn đất đồi núi vùng Thanh - Nghệ - Tĩnh bằng phương trình mất
116
đất phổ dụng và hệ thống thông tin địa lý. Tạp chí các khoa học về trái đất,
NXB Viện KH&CN Việt Nam, tập 34, số 1, tr. 31 – 37, 2012.
[11] Phạm Hùng., Nghiên cứu ứng dụng kỹ thuật mô hình toán trong tính toán xói
mòn lưu vực ở Việt Nam. Luận án tiến sĩ kỹ thuật, Trƣờng Ðại học Thủy lợi.
Hà Nội, 2001.
[12] Nguyễn Văn Khiết., Nghiên cứu đặc điểm xói mòn mặt khởi đầu dưới một số
thảm thực vật tại Lương Sơn - Hòa Bình. Luận văn thạc sỹ khoa học lâm
nghiệp, Đại học Lâm nghiệp, 2009.
[13] Phạm Hữu Tỵ và Hồ Kiệt., Mô phỏng rủi ro xói mòn vùng cảnh quan đồi núi
trên cơ sở sử dụng số liệu viễn thám và mô hình mất đất phổ quát hiệu chính
(RUSLE). Tạp chí khoa học đại học Huế số 48, 2008.
[14] ASCE., Sediment Transportation Mechanics: Introduction and Properties of
Sediment. ASCE. Journal of the Hydraulics Division, 1975.
[15] Rattan Lal., Soil erosion in the tropics. Principles and management, 1990.
[16] MacDonald L., Tài liệu tập huấn về thủy văn và quản lý lưu vực. Trƣờng Đại
học Lâm Nghiệp, 2009.
[17] Hudson., Bảo vệ đất và chống xói mòn (Đào Trọng Năng và Nguyễn Kim
Dung dịch). NXB Khoa học kỹ thuật. Hà Nội, 1981.
[18] Iman R. L and Helton J. C., An investigation of uncertainty and sensitivity
analysis techniques for computer models. Risk Analysis, 8 (1), 71-90, 1988.
[19] Nguyễn Văn Dũng và Trần Đức Viên., Ảnh hưởng của mưa và một số
phương thức sử dụng đất đến xói mòn đất và thu nhập của người dân ở vùng
đất dốc Tân Minh - Đà Bắc - Hoà Bình. Tạp chí Nông nghiệp và PTNT, Kỳ 1
tháng 12, Tr 36-38, 2003.
[20] Trần Quốc Vinh., Nghiên cứu sử dụng viễn thám RS và hệ thống thông tin địa
lý GIS để đánh giá xói mòn đất huyện Tam Nông tỉnh Phú Thọ. Luận án tiến
sỹ nông nghiệp, 2012.
117
[21] Trƣơng Văn Cảnh., Nghiên cứu ảnh hưởng xói mòn đất của lưu vực sông Cu
Đê đến sản xuất nông nghiệp. Đề tài cấp bộ Mã số: Đ2013-03-46-BS, 2014.
[22] Lane and Nearing., Profil Model documentation. USDA Water Erosion
Prediction project (WEEP) NSERL report NO. 2. Nation Soil erosion
Laboratory, USDA- ARS. W. Lafayette. IN, 1989.
[23] Rendon and Herredo., Unit sediment graph, 1978.
[24] Schramm., Ein Erosionsmodell mit räumlich und zeitlich veränderlicher
Rillenmorphologie. Mitteilung des Instituts für Wasserbauund Kulturtechnik,
Universität Fridericana zu Karlsruhe, Heft 190, 1994.
[25] Wischmeier W. H and Smith D. D., Relation of soil properties to its
erodibility. Soil Science Society of America Proc. 33(1), 131-137, 1960.
[26] Bennett., Concepts of mathematical modeling of sediment yield. Water
Resources Research . Vol. 10, No. 3, June, p. 485-492, 1974.
[27] CREAMS., Modeling developed coastal watersheds with theagricultural non-
point source mode, 1980.
[28] Croley T. E. II., Unsteady overland sedimentation. Journal of Hydrology,
Elsevier, 56, 325-346, 1982.
[29] Baver L. D., Soil Sci. Soc. Am. Proc. 3, 330 – 333, 1939.
[30] Laws J. O., Measurements of the fall-velocity of water-drops and raindrops.
Transactions, American Geophysical Union 22, 1941.
[31] Zingg A. W., Degree and Length of Land Slope as it Affects Soil Loss in
Runoff. Agric. Eng., 21(2), 59-64, 1940.
[32] Nguyễn Văn Khiết., Nghiên cứu xác định vai trò của một số yếu tố liên quan
đến xói mòn đất ở nước ta. Tạp chí khoa học Lâm nghiệp, Viện KH Lâm
Nghiệp Việt Nam, 1, trang 3145 - 3153, 2014.
118
[33] Ellison W. D., Soil erosion studies. AGRICULTURAL ENGINEERING 28:
145-146, 197-201, 245-248, 297-300, 349-351, 402-405, 1947.
[34] Musgrave., The causes and effects of sedimentation in lake Decatur. State
Water Survey Division, 1947.
[35] Smith D. D., A parameter-efficient hydrologic infiltration model. Water
Resources Research 14 (3), 533 – 538, 1941.
[36] Wischmeier, W. H and Smith D. D., Predicting rainfall erosion losses. A
guide to conservation planning, agricultural Handbook No 537, US Dep. Of
Agri. Washington, D. C, 1978.
[37] Brown, L. C and Foster G. R., Storm erosivity using idealized intensity
distributions. Trans. ASAE, 30, 379-386, 1987.
[38] Renard K. G and Freimund J. R., Using monthly precipitation data to estimate
the R-factor in the revised USLE. Journal of Hydrology, 157, 287–306, 1994.
[39] Chu Đình Hoàng., Tác dụng xói mòn của đất giọt mưa. Kỹ thuật Thủy Lợi, số
5, 1963.
[40] Williams J. R., Sediment Yield Predicting with Universal Equation Using
Runoff Energy Factor. Present and Prospective Technology for Predicting
Sediment Yield and Sources, USDA, ARS-S-40, 1975.
[41] ANSWERS., A model for watershed planning. Transaction of the ASAE, Vol.
23, No. 4, p. 938-944, 1980.
[42] AGNPS., Modelling sediment yield from agricultural watersheds. J. Soil
Water Conserv., 37, p.113-117, 1980.
[43] KINEROS., A kinematic runoff and erosion model. U. S. Department of
Agriculture, Agricultural Research Service. ARS-77, 130 pp. In: Arnold, J.G.,
Srinivasan, R., Muttiah, R. S., Williams J. R., 1998. Large area hydrologic
119
modeling and assessment Part I: Modeldevelopment. Journal of American
Water Resources Association 34 (1), 73–89, 1981.
[44] Woolhiser et al., KINEROS, a kinematics runoff and erosion model.
Documentation and user manual – U.S Dept. Of Agriculture, Agri. Res.
Service ARS - 77, March, 1990.
[45] WEPP., USDA – Water Erosion Prediction Project. Hillslope profile and
watershed model documentation NSERL Rep. No. 10, USDA – ARS, West
Lafayette, Indiana, 1995.
[46] Foster and Meyer L. D., A closed- form soil erosion equation for upland
areas. In: H.W. Shen (editor and publisher) sedimentation (Einstein)
Colorado State University, Fort Colling, CO. Kapitel 12, 1972.
[47] Morgan R. P. C., Quinton, J. N., Smith, R. E., Govers, G., Poesen, J. W. A.,
Auerswald, K., Chisci, G., Torri, D., Styczen, M. E., Folly, A. J. V., The
European soil erosion modell (EUROSEM): documentation and user guide.
Silsoe College, Cranfield University, 1998.
[48] Morgan P. C., Soil erosion modelling with EUROSEM at Embori and
Mukogodo catchments. Kenya, 2006.
[49] Nguyễn Quang Mỹ., Xói mòn đất hiện đại và các biện pháp chống xói mòn.
NXB Đại Học Quốc Gia, 2005.
[50] Lại Vĩnh Cẩm., Soil erosion study in North West region of Viet Nam by
intergrating waterhed analysis and universal soil loss equantion (USLE). Tạp
chí khoa học ĐHQG HN, KHTN số XI, 2000.
[51] Lê Mạnh Hùng và nnk., Kết quả ứng dụng mô hình SWAT trong tính toán xói
mòn lưu vực Mekong. Tạp chí khoa học công nghệ và thủy lợi số 12, 2012.
[52] Han F. et al., The WEPP Model Application in a Small Watershed in the Loess
Plateau. 2016.
120
[53] Hunink J. E. et al., Targeting of intervention areas to reduce reservoir
sedimentation in theTana catchment (Kenya) using SWAT. 2013.
[54] Woolhiser D. A and J. A. Liggett., Unsteady, one dimensional flow over a
plane – the rising hydrograph. Water Resources Research 3, 753-761, 1967.
[55] Foster et al., Modeling the erosion process. Chapter 8. In: Hydrologic
modeling of small watersheds, edited by C. T. Haan et al,. ASAE Monograph
No. 5, 1981.
[56] Rezzolla, L., Numerical methods for the solution of partial differential
equations. Lecture notes for the COMPSTAR School of Computational
Astrophysics, 8-13/02/2010. Caen, France, 2011.
[57] Kibler, D. F and Woolhiser, D. A., The kinematic cascade as a hydrologic
model. Hydrology Paper No. 39, Colorado State University, Fort Collins,
Colorado, USA, 1970.
[58] Einstein H. A., Deposition of Suspended Particles in a Gravel Bed, Jour.
Hydraulics Division, American Society of Civil Engineers, p 6102, 1968.
[59] Lane L. J. and Shirley E. D., Erosion and sediment yield equations: solutions
for overland flow. Paper presented at the Workshop on USLE Replacement.
NSEL, West Lafayette, IN, p. 22, 1985.
[60] Harmon, W. C. and Meyer, L. D., Multiple-Intensity Rainfall Simulator for
Erosion Research on Row Sideslopes. Transactions of the ASABE, 22, 100-
103, 1979.
[61] Wischmeler W. H., Estimating thesoil loss equation's cover and management
factor for undisturbed areas. In: Present and proslogy for predicting sediment
yields and sources. ARS-5-40. USDA Agricultural Research Service, p. 118-
124, 1975.
121
[62] Meyer L. D and Wischmeier W. H., Mathematical simulation of the process of
soil erosion by water. Trans. ASAE, 12, 754-758, 1969.
[63] Parsons A.J., Abrahams, A.D., and Simanton, J.R., Microtopography and soil-
surface materials on semi-arid piedmont hillslopes, southern Arizona. Journal
of Arid Environments, v. 22, p. 107-115, 1992.
[64] Engelund F. and Hansen E., Comparison Between Similarity Theory and
Regime Formulae. Basic Research - Progress Report No. 13, Jan. 1967.
Hydraulic Laboratory, Techn. Un. of Denmark, 1967.
[65] Werner M. G. F., Hunter N. M. and Bates. P. D., Identifiability of distributed
floodplain roughness values in flood extent estimation. Journal of Hydrology,
314, 139-157, 2005.
[66] Harvey G. and Jan T., An Introduction to Computer Simulation Methods. Part
2, Applications to Physical Systems, 1988.
122
PHỤ LỤC
1. Các bảng phân chia lƣu vực
Lƣu vực
Diện tích (km2)
Độ dốc Ix (%)
Độ dốc Iy (%)
Chiều dài sông chính (m) 1.171 1.265 126 947 1.221 230 1.095 126 877 547 171 668 894 1.062 1.180 629 784 945 740 602 714 828 843 265 871 620 1.143 931 690 383 236 561 875
Độ cao trung bình (m) 806 439 661 601 640 829 779 691 639 280 308 760 493 785 302 363 408 714 347 633 815 410 258 836 435 868 353 557 1223 1170 789 849 826
271.8 437.6 33.8 532.1 474.2 187.7 503.8 94.1 449.5 379.5 452.3 327.6 375.4 349.2 359.2 263.1 264.0 394.8 349.4 268.9 282.9 211.3 273.2 335.4 204.3 261.4 305.1 283.2 344.8 353.2 241.1 360.1 406.8
Chiều dài lƣu vực (km) 12.02 29.85 106.67 2.70 12.58 47.30 11.69 60.95 40.14 12.87 56.14 6.71 9.31 11.45 11.39 12.21 93.88 10.16 9.51 5.32 14.34 20.10 8.35 4.83 8.08 77.42 21.28 12.37 40.81 21.72 75.93 12.55 6.58
141.18 85.09 8.16 185.32 101.30 36.28 98.60 17.30 54.35 92.04 566.65 107.35 102.80 82.33 156.23 87.57 35.57 114.73 105.44 130.43 99.81 59.03 134.92 91.93 79.94 26.72 96.02 79.21 34.44 43.55 13.93 62.44 149.06
1 2 8 3 4 6 5 7 9 12 16 14 13 11 15 17 20 25 22 27 21 18 24 30 19 10 26 23 29 31 33 37 35
14.08 37.76 13.44 2.56 15.36 10.88 12.80 7.68 35.20 7.04 9.60 4.48 8.32 12.16 13.44 7.68 73.60 9.60 7.04 3.20 10.24 16.64 7.04 1.28 7.04 48.00 24.32 11.52 28.16 8.32 17.92 7.04 5.76
Bảng 1: Phân chia các tiểu lƣu vực trong lƣu vực Nậm Sập
123
Lƣu vực
Diện tích (km2)
Độ dốc Ix (%)
Độ dốc Iy (%)
39 38 28 34 41 32 40 36 48 44 46 45 43 51 47 55 49 50 42 59 52 56 58 53 63 61 64 60 54 57 65 62 66 69 68 67 71 70 72
11.52 1.92 14.72 22.40 14.72 16.00 10.24 97.28 60.80 7.04 9.60 5.76 10.88 7.04 5.12 8.32 7.68 3.84 22.40 4.48 7.68 4.48 5.12 20.48 10.24 1.28 1.70 8.96 12.80 26.24 16.64 7.04 7.68 8.96 6.40 24.96 19.84 28.16 10.24
Chiều dài sông chính (m) 433 453 711 439 446 202 612 1.034 333 211 496 695 397 531 377 530 685 528 758 390 502 579 435 357 492 260 130 465 623 407 588 486 376 632 553 574 629 695 1.070
Độ cao trung bình (m) 836 723 902 1151 1038 1001 941 1081 1080 978 1128 775 1184 855 843 788 583 617 1.127 853 1081 816 685 1088 883 722 749 923 1135 967 942 917 732 922 847 1064 596 491 890
301.0 443.8 265.5 295.4 277.0 286.8 238.4 344.4 481.5 242.3 352.0 297.3 305.7 407.6 399.3 225.8 463.4 305.7 334.3 205.4 320.7 346.3 326.3 335.0 374.5 306.6 426.2 236.2 316.4 196.0 282.2 280.1 183.3 284.9 286.3 280.5 287.3 316.1 375.5
Chiều dài lƣu vực (km) 26.61 4.24 20.70 51.03 33.00 79.21 16.73 94.08 182.58 33.36 19.35 8.29 27.41 13.26 13.58 15.70 11.21 7.27 29.55 11.49 15.30 7.74 11.77 57.37 20.81 4.92 65.60 19.27 20.55 64.47 28.30 14.49 20.43 14.18 11.57 43.48 31.54 40.52 9.57
430.45 144.74 50.02 29.48 45.93 15.05 66.30 30.80 582.17 34.48 64.11 93.62 44.02 70.63 59.76 60.94 100.23 124.89 57.80 61.91 74.76 106.37 91.29 36.79 39.81 100.74 122.91 72.10 61.95 25.30 43.31 64.68 56.78 72.41 71.42 41.60 54.78 55.10 173.89
124
Bảng 2: Phân chia các tiểu lƣu vực trong lƣu vực Phiêng Hiềng
Lƣu vực
Diện tích (km2)
Độ cao trung bình (m)
Độ dốc Ix(%)
Độ dốc Iy (%)
Chiều dài lƣu vực (km) 2.40 2.47 2.06 2.18 6.60 3.71 5.21 2.78 2.85 2.65 2.90 2.50 3.25 3.42 1.31 1.62 2.81 57.69 3.05 3.46 1.97 2.11 4.15 2.30 3.40 2.70 8.33 4.11 3.76 6.06 2.13 2.28 1.31 3.32 3.70 2.42 4.76 3.07 4.87 3.49 2.22 3.41 157.53 3.53
1 2 4 3 5 8 6 7 11 10 9 18 14 13 17 15 19 22 23 20 21 16 26 25 29 27 12 33 34 24 35 28 37 38 36 40 32 39 30 31 41 42 45 44
3.00 2.00 2.00 3.00 6.50 4.00 5.00 3.50 2.50 2.50 3.50 2.00 4.00 2.50 0.50 1.50 2.50 22.50 3.50 3.00 2.00 2.00 4.00 1.00 3.00 2.50 12.00 3.00 3.50 6.00 2.00 3.00 0.50 3.00 4.00 0.50 5.00 2.50 7.50 5.00 2.00 3.50 34.50 3.00
66 270 242 364 175 244 189 259 244 214 299 231 229 176 253 236 239 278 224 176 348 222 169 149 163 218 138 155 196 175 235 341 194 152 202 119 160 185 188 285 173 184 278 198
Chiều dài sông chính (m) 1.251 809 969 1.379 985 1.079 960 1.261 876 945 1.207 801 1.232 730 382 926 889 390 1.148 866 1.015 950 963 435 882 927 1.440 730 931 990 940 1.314 381 904 1.082 207 1.050 814 1.541 1.431 901 1.026 219 850
114 113 111 133 97 108 118 114 102 105 95 96 96 72 91 88 99 76 89 95 113 94 101 72 111 106 99 71 68 93 80 113 74 81 110 74 69 77 102 111 71 87 64 91
564 440 442 1150 723 976 824 1050 944 732 1139 1188 1.265 1.103 1.111 1.255 1.212 980 1.016 1.481 1.366 1.521 1.343 971 1.202 1.236 1.433 1.714 1.488 1.051 1.484 1.455 1.097 1.416 907 1.021 1.800 1.229 1.829 1.695 1.347 1.079 833 1.193
125
Lƣu vực
Diện tích (km2)
Độ cao trung bình (m)
Độ dốc Ix(%)
Độ dốc Iy (%)
Chiều dài lƣu vực (km) 2.89 7.30 3.39 863.64 4.20 80.58 3.24 2.43 5.66 1.71 5.72 4.56 2.19 7.67 1.98 6.54 3.35 8.23 3.89 4.73 1.91 1.61 3.68 3.98 1.34 1.24 3.78 4.59 4.48 2.71 1.01 5.88 2.34 2.15 3.38 2.81 9.64 2.66 9.20 74.07 2.01 2.36 8.88 2.83 1.79 2.13 4.69
46 43 47 49 50 51 52 48 54 58 57 53 60 56 55 61 59 62 64 63 66 67 65 69 71 70 68 76 73 77 75 72 79 81 82 74 78 80 84 87 85 86 89 83 92 88 93
3.50 11.50 2.00 19.00 4.50 19.50 2.50 4.00 11.50 3.00 6.00 4.50 1.00 15.00 2.00 5.50 5.50 9.50 4.00 5.00 2.00 2.00 4.00 5.00 0.50 1.50 4.00 1.50 2.50 2.00 0.50 4.50 1.50 1.00 3.00 2.50 11.50 2.50 7.00 6.00 2.00 2.00 5.00 2.50 1.00 2.00 4.50
293 164 167 62 171 307 140 225 182 324 141 223 216 191 283 143 273 130 204 213 256 319 187 250 217 420 208 116 169 190 231 174 248 201 241 203 156 178 130 210 223 322 98 214 263 233 221
Chiều dài sông chính (m) 1.209 1.575 590 22 1.071 242 771 1.647 2.033 1.751 1.049 986 457 1.956 1.011 841 1.643 1.155 1.028 1.057 1.049 1.245 1.087 1.257 374 1.208 1.057 327 558 737 495 765 641 465 888 891 1.193 940 761 81 993 849 563 882 560 940 960
80 88 92 44 74 91 73 95 96 88 85 93 94 90 104 92 90 102 114 101 76 98 77 95 100 90 87 93 84 84 83 78 80 64 104 82 89 82 82 47 92 113 71 85 76 105 110
2.248 1.370 1.068 697 1.178 711 1.189 1.255 1.733 1.377 683 583 384 1.087 775 792 1.039 727 692 573 806 580 812 656 505 442 411 291 361 447 226 782 213 170 519 801 301 878 322 154 621 652 681 456 327 599 779
126
Lƣu vực
Diện tích (km2)
Độ cao trung bình (m)
Độ dốc Ix(%)
Độ dốc Iy (%)
90 91 94
5.00 4.50 4.50
Chiều dài lƣu vực (km) 6.51 4.57 5.81
116 202 178
68 86 83
Chiều dài sông chính (m) 768 984 774
571 510 205
2. Trích dẫn đoạn chƣơng trình mô phỏng vận chuyển bùn cát trên lƣu vực
Mô đun mô phỏng vận chuyển bùn cát trên lƣu vực
#ifdef _DEBUG
void CErosion::AssertValid() const
{
CDocument::AssertValid();
}
void CErosion::Dump(CDumpContext& dc) const
{
CDocument::Dump(dc);
}
#endif //_DEBUG
BOOL CErosion::Init( ERGSETSTRUCT *pErgSet, VARSTRUCT *pVar,
int typ )
{
Empty();
int nr;
m_Typ = typ; // 386 oder 387
if( m_Typ == 181 || m_Typ == 387 )
m_MaterTyp = ST_SEDIMENT;
if( m_Typ == 289 )
m_MetabTyp = ST_PHOSPHORUS;
if( typ != 386 && typ != 387 && typ != 289 ) {
AfxMessageBox(“Error" );
127
return FALSE;
}
m_pErgSet = pErgSet;
delta_t = pErgSet->delta_t; // [thời gian h]
m_pTfl = pErgSet->pTFL->tfl;
m_pVar = pVar;
CPCDData *pPCD = pErgSet->pPCD;
m_NumbTfl = pErgSet->NumbTFL;
m_pPcd = pErgSet->pPCD;
CString strPfad = pErgSet->pfad;
strPfad += "Erosion.ers";
BOOL bAddToMRU = FALSE;
SetPathName( strPfad, bAddToMRU );
pErosion = new EROSSTRUCT[m_NumbTfl+1];
memset( pErosion, 0, (m_NumbTfl+1)*sizeof(EROSSTRUCT) );
for( nr=0; nr<=m_NumbTfl; nr++ ) {
pErosion[nr].NumbHeadSeg = (int) ceil(m_pVar[nr].tx / pErgSet->delta_t);
pErosion[nr].NumbNextSeg = (int) ceil(m_pVar[nr].txy / pErgSet->delta_t);
if( pErosion[nr].NumbNextSeg > 0 ) {
CGG::Initialize();
TFLNRLISTE *pList = new TFLNRLIST[pErgSet->numbTFL+1];
GGCreateFile( pErgSet, pList, FALSE, EROSION );
delete pList;
m_StartMin = pErgSet->req _minutes;
m_EndMin = pErgSet->end_minutes;
if(Tak!= NULL ) {
for( int i=0; i<=numb_Tak; i++ )
Tak[i].Time.FillEmpty( m_StartMin, m_EndMin, (MINUTES) (delta_t*60)
);
128
}
if( m_Typ == 181 || m_Typ == 182 ) {
for( nr=1; nr<=m_NumbTfl; nr++ ) {
int pcdnr = m_pTfl[nr].Proc.ErosNr;
if( pcdnr > 0 ) {
pErosion[nr].SLFAC = m_pTfl[nr].iy/1000;
pErosion[nr].SinALPHAW= (float) ( m_pTfl[nr].iy / sqrt(
m_pTfl[nr].iy*m_pTfl[nr].iy + 1000*1000 ) );
pErosion[nr].CosALPHAW= (float) ( 1000 / sqrt( m_pTfl[nr].iy*m_pTfl[nr].iy +
1000*1000 ) );
pErosion[nr].LI = m_pTfl[nr].area*500000 / m_pTfl[nr].flength;
pErosion[nr].AREA = m_pTfl[nr].area*1000000;
pErosion[nr].I = m_pTfl[nr].ix / 1000;
pErosion[nr].B = pErosion[nr].AREA / m_pTfl[nr].flength;
pErosion[nr].Qf = m_pVar[nr].qo1 + m_pVar[nr].qi1 + m_pVar[nr].qu1;
pErosion[nr].MQC = pPCD->GetPcdVar( 386, pcdnr, MP_MQC );
pErosion[nr].IRF = pPCD->GetPcdVar( 386, pcdnr, MP_IRF );
pErosion[nr].D50 = pPCD->GetPcdVar( 386, pcdnr, MP_D50 ) / 1000; ,
pErosion[nr].Visco = pPCD->GetPcdVar( 386, pcdnr, MP_VISCO );
pErosion[nr].CI = pPCD->GetPcdVar( 386, pcdnr, MP_CI );
pErosion[nr].RR = pPCD->GetPcdVar( 386, pcdnr, MP_RR );
pErosion[nr].KECG = pPCD->GetPcdVar( 386, pcdnr, MP_KECG );
pErosion[nr].Fs = pPCD->GetPcdVar( 386, pcdnr, MP_FS );
pErosion[nr].KET = pPCD->GetPcdVar( 386, pcdnr, MP_KET );
pErosion[nr].Cb = pPCD->GetPcdVar( 386, pcdnr, MP_RE );
pErosion[nr].Omega = pPCD->GetPcdVar( 386, pcdnr, MP_OMEGA );
pErosion[nr].Omega*= 0.01f * RhoW/(RhoS-RhoW);
pErosion[nr].fdep = pPCD->GetPcdVar( 386, pcdnr, MP_FDEP );
pErosion[nr].BBG = 0.5f;
129
pErosion[nr].Kfac = 0.5f; // [đơn vị kgh/Nm2]
pErosion[nr].Cfac = 0.5f; //
pcdnr = m_pTfl[nr].Proc.UsleNr;
if( pcdnr > 0 && pPCD->GetProc( 387, pcdnr ) ) {
pErosion[nr].BBG = pPCD->GetPcdVar( 387, pcdnr, MP_BBG );
pErosion[nr].Kfak = pPCD->GetPcdVar( 387, pcdnr, MP_KFAC);
pErosion[nr].Cfak = pPCD->GetPcdVar( 387, pcdnr, MP_CFAC );
pErosion[nr].KI = pPCD->GetPcdVar( 387, pcdnr, MP_KI );
}
float delta = 0.0028;
pErosion[nr].TAU0 = delta * G * ( RhoS - RhoW ) * pow(
pErosion[nr].D50, 0.33 ); }
}
}
if( m_Typ == 289 ) {
for( nr=1; nr<=m_NumbTfl; nr++ ) {
if( m_pTfl[nr].Proc.Matyp == 289 ) {
pErosion[nr].AREA =
m_pTfl[nr].area*1000000;
int pcdnr = m_pTfl[nr].Proc.MatNr;
}
{
int pcdnr = m_pTfl[nr].Proc.MatNr;
pErosion[nr].t_ow = m_pPcd->GetPcdVar( 286, pcdnr, MP_T_O );
pErosion[nr].t_dr = m_pPcd->GetPcdVar( 286, pcdnr, MP_T_I );
pErosion[nr].t_gw = m_pPcd->GetPcdVar( 286, pcdnr, MP_T_U ); }
else { }
}
}
130
m_bShowQC = TRUE;
m_bShowF = FALSE;
return TRUE;
}
Mô đun xói mòn liên rãnh
void CErosion::PotentialRillErosion( int nr )
{
float OFRQ = pErosion[nr].OFRQ;
float RE = pErosion[nr].RE;
pErosion[nr].EC = 1000 * 3600 * RE * RhoW * OFRQ;
}
void CErosion::DepositionRate( int nr )
{
float VS = pErosion[nr].VS;
float QSI = pErosion[nr].QSI;
float OFRQ = pErosion[nr].OFRQ;
float TC = pErosion[nr].TC;
pErosion[nr].DEP = -0.5f * ( VS / OFRQ ) * ( TC - QSI );
}
void CErosion::RillErosion( int nr )
{
float TK = pErosion[nr].TC;
float QSI= pErosion[nr].QSI;
float EC = pErosion[nr].EC;
float term1 = 1/sqrt( pErosion[nr].SLFAC );
float term2 = pErosion[nr].OFRQ * term1 * RFM;
float OFRH = pow( term2, 0.6 );
float GammaW = RhoW * G;
float TAUn = GammaW * OFRH * pErosion[nr].SLFAC;
131
float SLFAC = pErosion[nr].SLFAC;
float OFRQ = pErosion[nr].OFRQ;
float Kr = pErosion[nr].Kfac;
float C = pErosion[nr].Cfac;
float Cb = pErosion[nr].Cb;
pErosion[nr].ER = 1000 * 3600 * Cb * OFRQ * SLFAC * Kr * C;
if( pErosion[nr].ER < 0 )
int t=1;
}
void CErosion::sedimenttransport( int nr )
{
float EI = pErosion[nr].EI;
float ER = pErosion[nr].ER;
float DEP = pErosion[nr].DEP;
pErosion[nr].Phi = EI + ER + DEP;
if( pErosion[nr].Phi < 0 )
pErosion[nr].Phi = 0;
if( pErosion[nr].Phi > 0 )
int t=1;
pErosion[nr].QSI = pErosion[nr].Phi * m_pTfl[nr].flength;
}
void CErosion::ChannelTransport181( int nr )
{
if( pErosion[nr].OFRQ > 0 ) {
float SLG = pErosion[nr].SLG;
float D50 = pErosion[nr].D50;
float Z = 0.3f;
float term1 = (float) ( 1 / sqrt( SLG ) );
float term2 = pErosion[nr].OFRQ * term1 * RFM;
132
float OFRG = (float) pow( term2, 0.6 );
float VG = (float) ( (1 / RFM) * pow( OFRG, 0.666667 ) * sqrt( SLG )
pErosion[nr].WIDTH = 2.05 * pow( Z, -0.625) * pow( (1 + Z*Z), 0.125 ) *
pow( pErosion[nr].Qab * m_pTfl[nr].kst * term1, 0.375 );
float qch = pErosion[nr].Qab / pErosion[nr].WIDTH;
float VS = (float) ( RhoR * G * D50 * D50 / ( 18.0f * m_VISCO ) );
pErosion[nr].VS = VS;
float delta = 0.000285;
pErosion[nr].DEPTH = pErosion[nr].fdep * pow( qch *
m_pTfl[nr].kst * (1/sqrt(SLG)), 0.6 );
pErosion[nr].TAUch = RhoW * G * pErosion[nr].DEPTH * SLG;
float Ef = pErosion[nr].TAUch / ( D50 * (RhoS-RhoW)*G );
float nue = 0.7f * pow( Ef, 0.98 );
pErosion[nr].TCG = new * pErosion[nr].Omega * pErosion[nr].TAUch *
3.600 * (VG*VG/VS) * (1/G); }
else
pErosion[nr].TC = 0;
}
float D50 = pErosion[nr].D50;
float GammaW= RhoW * G;
float SLG = pErosion[nr].SLG;
float term1 = (float) ( 1 / sqrt( SLG ) );
float term2 = pErosion[nr].OFRQ * term1 * RFM;
float OFRG = (float) pow( term2, 0.6 );
float TAUch = pErosion[nr].TAUch;
float US = (float) ( sqrt( TAUch / RhoS ) );
float RE = (float) ( US * D50 / m_VISCO );
float TAU0 = pErosion[nr].TAU0;
133
float kch = 300;
pErosion[nr].EGR = (float) ( 3600 * pErosion[nr].WIDTH * kch * pow( (
1.35 * TAUch - TAU0 ), 1.05 ) );
}
{
float t = 0.1f;
float SLFAC = pErosion[nr].SLFAC;
float Width= (float) ( 2.05 * pow( t, -0.625) * pow( (1 + t*t), 0.125) *
pow( (QP*RFM)/sqrt(SLFAC), 0.375 ) );
}
else {
if( pErosion[nr].QSI > pErosion[nr].TC ) {
Deposition( nr );
}
else {
RillErosion( nr );
}
}
}
void CErosion::Erosion( int nr, int Timestep, VARSTRUCT *pVar )
{
pErosion[nr].ER = pErosion[nr].EI = pErosion[nr].ET =
pErosion[nr].QSI = 0;
pErosion[nr].TC = pErosion[nr].DEP = pErosion[nr].Phi = 0;
pErosion[nr].QSGR = pErosion[nr].WEIGHT =
pErosion[nr].OFRQ = pVar[nr].qo1 / pErosion[nr].B;
pErosion[nr].INTPH = pVar[nr]. rain / delta_t;
134
pErosion[nr].PD = pVar[nr].rain;
pErosion[nr].Qab = pVar[nr].qab / pVar[nr].itseg;
pErosion[nr].VAUQ = pVar[nr].qo1 * 3600 * delta_t;
if( pErosion[nr].OFRQ > 0 )
int t=1;
if( pErosion[nr].OFRQ < 0 )
int t=1;
if( pErosion[nr].Qab > 0 )
int t=1;
if( pErosion[nr].Qab < 0 )
int t=1;
if( m_pTfl[nr].Proc.ErosTyp == 386 )
AreasOder 386 ( nr, timestep, pVar );
else if ( m_pTfl[nr].Proc.ErosTyp == 387 )
AreasOder 387( nr, timestep, pVar );
if( pErosion[nr].QSI > 0 )
if( pErosion[nr].Phi > 0 ) // [đơn vị g/m2*h]
int t=1;
if( pErosion[nr].Phi < 0 )
int t=1;
float Freight = pErosion[nr].Phi * pErosion[nr].AREA; // [đơn vị g/h]
if( pErosion[nr].Qab > 0 && Freight > 0 ) {
if ( m_pTfl[nr].Proc.ErosTyp == 386 )
if ( m_pTfl[nr].Proc.ErosTyp == 387 )
}
}
Đoạn chƣơng trình mô tả vận chuyển bùn cát trong sông
if(Tak != NULL && m_pTfl[nr]. LevelNr > 0 ) {
CStation *pTak = GetTakPtr( m_pTfl[nr].LevelNr );
135
if( pTak ) {
int TakNr = pTak->nr;
int t=1;
}
}
if( pErosion[nr].EI || pErosion[nr].ET || pErosion[nr].ER || pErosion[nr].DEP ||
pErosion[nr].Phi ) {
pErosion[nr].S_EI += pErosion[nr].EI * delta_t; // đơn vị g/m2
pErosion[nr].S_ET += pErosion[nr].ET * delta_t; // đơn vị g/m2
pErosion[nr].S_ER += pErosion[nr].ER * delta_t; // đơn vị g/m2
pErosion[nr].S_DEP+= pErosion[nr].DEP * delta_t; // đơn vị g/m2
pErosion[nr].S_PHI+= pErosion[nr].Phi * delta_t; // đơn vị g/m2
}
}
BOOL IsValid( float value )
{
if( value < 0 || value > 1E32 )
return FALSE;
return TRUE;
}
void CErosion( int tflnr, int time )
{
float Flow=0.0f, Fldr=0.0f, Flgw=0.0f; // [đơn vị g/s]
float Qlow=0.0f, Qldr=0.0f, Qlgw=0.0f; // [đơn vị m3/s]
float Inflowfoc = 0;
float Inflow = 0;
float InflowFreight = 0;
if( pErosion[tflnr].MetabEntryLevelNr > 0 ) {
136
CStation *pTak = m_pErgSet->pOrgDoc->CGGC.GetTakPtr(pErosion[tflnr].
MetabEntryLevel Nr );
if( pTak ) {
int TakNr = pTak->nr;
MINUTES minute = (MINUTES) ( m_pErgSet->req_minutes + Timestep*
delta_t * 60 );
Infoww *= pErosion[tflnr]. MetabEntryshare;
}
}
if( time == 0 ) {
pErosion[tflnr].Fow = pErosion[tflnr].CStartFloor * m_pVar[tflnr].qo1;
pErosion[tflnr].Fdr = pErosion[tflnr].CStartFloor * m_pVar[tflnr].qi1;
pErosion[tflnr].Fgw = pErosion[tflnr].CStartFloor * m_pVar[tflnr].qu1;
Qlow = pErosion[tflnr].Qlow = pErosion[tflnr].Qow = m_pVar[tflnr].qo1;
Qldr = pErosion[tflnr].Qldr = pErosion[tflnr].Qdr = m_pVar[tflnr].qi1;
Qlgw = pErosion[tflnr].Qlgw = pErosion[tflnr].Qgw = m_pVar[tflnr].qu1;
}
float Nep = m_pVar[tflnr].Neo + m_pVar[tflnr].No + m_pVar[tflnr].
if( Nep > 0 ) {
Qlow = m_pVar[tflnr].Neo * m_pTfl[tflnr].area * 1000 / 3600;
Qldr = m_pVar[tflnr].No * m_pTfl[tflnr].area * 1000 / 3600;
Qlgw = m_pVar[tflnr].New * m_pTfl[tflnr].area * 1000 / 3600;
m_pVar[tflnr].New / Nep; }
float mFlow = ( Fzow + pErosion[tflnr].Flow ) / 2;
float mQlow = ( Qzow + pErosion[tflnr].Qlow ) / 2;
float Qow = pErosion[tflnr].Qow; // [đơn vị m³/s]
float Fow = pErosion[tflnr].Fow;
IsValid( pErosion[tflnr].Qow );
IsValid( pErosion[tflnr].Fow );
137
float kow = 1 - (float) exp( -delta_t / pErosion[tflnr].kow );
pErosion[tflnr].Fow = Fow + ( mFlow - Fow ) * kow;
pErosion[tflnr].Flow = Flow;
pErosion[tflnr].Qow = m_pVar[tflnr].qo1;
pErosion[tflnr].Qlow = Qlow;
pErosion[tflnr].S_Qlow += Qlow;
pErosion[tflnr].S_Flow += Flow;
pErosion[tflnr].S_Qow += pErosion[tflnr].Qow;
pErosion[tflnr].S_Fow += pErosion[tflnr].Fow;
float mFzdr = ( Fzdr + pErosion[tflnr].Fzdr ) / 2;
float mQzdr = ( Qzdr + pErosion[tflnr].Qzdr ) / 2;
float Qdr = pErosion[tflnr].Qdr;
float Fdr = pErosion[tflnr].Fdr;
IsValid( pErosion[tflnr].Qdr );
IsValid( pErosion[tflnr].Fdr );
float kdr = 1 - (float) exp( -delta_t / pErosion[tflnr].kdr );
pErosion[tflnr].Fdr = Fdr + ( mFzdr - Fdr ) * kdr;
pErosion[tflnr].Fzdr = Fzdr;
pErosion[tflnr].Qdr = m_pVar[tflnr].qi1;
pErosion[tflnr].Qzdr = Qzdr;
pErosion[tflnr].S_Qzdr += Qzdr;
pErosion[tflnr].S_Fzdr += Fzdr;
pErosion[tflnr].S_Qdr += pErosion[tflnr].Qdr;
pErosion[tflnr].S_Fdr += pErosion[tflnr].Fdr;
float mFzgw = ( Fzgw + pErosion[tflnr].Fzgw ) / 2;
float mQzgw = ( Qzgw + pErosion[tflnr].Qzgw ) / 2;
float Qgw = pErosion[tflnr].Qgw;
float Fgw = pErosion[tflnr].Fgw;
IsValid( pErosion[tflnr].Qgw );
138
IsValid( pErosion[tflnr].Fgw );
float kgw = 1 - (float) exp( -delta_t / pErosion[tflnr].kgw );
pErosion[tflnr].Fgw = Fgw + ( mFzgw - Fgw ) * kgw;
pErosion[tflnr].Fzgw = Fzgw;
pErosion[tflnr].Qgw = m_pVar[tflnr].qu1;
pErosion[tflnr].Qzgw = Qzgw;
pErosion[tflnr].S_Qzgw += Qzgw;
pErosion[tflnr].S_Fzgw += Fzgw;
pErosion[tflnr].S_Qgw += pErosion[tflnr].Qgw;
pErosion[tflnr].S_Fgw += pErosion[tflnr].Fgw;
float Freight = pErosion[tflnr].Fow + pErosion[tflnr].Fdr +
pErosion[tflnr].Fgw;
if( Freight < 0 )
int t=1;
pErosion[tflnr].QSI = 3600 * Freight / m_pTfl[tflnr].flength;
if(Tak != NULL )
if(Tak!= NULL ) {
int tflrec = m_pTfl[tflnr].LevelEzgRec;
CStation *pTak = GetTakPtr( m_pTfl[tflrec].LevelNr );
if( pTak &&Concentration > 0) {
int TakNr = pTak->nr;
}
}
if(Tak!= NULL && m_pTfl[tflnr].LevelNr > 0 ) {
CStation *pTak = GetTakPtr( m_pTfl[tflnr].LevelNr );
if( pTak ) {
int TakNr = pTak->nr;
int t=1;
139
else
MINUTES minute = (MINUTES) ( m_pErgSet->inq_minutes +
timestep * delta_t * 60 );
}
}
}
float D50 = pErosion[nr].D50;
float GammaW= RhoW * G;
float SLG = pErosion[nr].SLG;
float term1 = (float) ( 1 / sqrt( SLG ) );
float term2 = pErosion[nr].OFRQ * term1 * RFM;
float OFRG = (float) pow( term2, 0.6 );
float TAUch = pErosion[nr].TAUch;
float US = (float) ( sqrt( TAUch / RhoS ) );
float RE = (float) ( US * D50 / m_VISCO );
float TAU0 = pErosion[nr].TAU0;
float kch = 300;
pErosion[nr].EGR = (float) ( 3600 * pErosion[nr].WIDTH * kch * pow( (
1.35 * TAUch - TAU0 ), 1.05 ) );
}
{
float t = 0.1f;
float SLFAC = pErosion[nr].SLFAC;
float Width= (float) ( 2.05 * pow( t, -0.625) * pow( (1 + t*t), 0.125) *
pow( (QP*RFM)/sqrt(SLFAC), 0.375 ) );
}
else {
if( pErosion[nr].QSI > pErosion[nr].TC ) {
140
Deposition( nr );
}
else {
RillErosion( nr );
}
}
}
void CErosion::ChannelDeposition( int nr )
{
float FS = pErosion[nr].Fs;
float D50 = pErosion[nr].D50;
float BETA= 1.0f;
float TCGR= pErosion[nr].TCGR;
float QSGRN= pErosion[nr].QSGRN;
float RhoF=RhoW;
float Visco = pErosion[nr].Visco;
float RhoR2 = ( RhoS - RhoF ) / RhoF;
float D = D50 * (float) pow( (RhoR2 * G / (Visco*Visco) ), 0.333 );
float OFRQ = pErosion[nr].OFRQ;
float VSG = (float) ( 11 * Visco * ( sqrt(1+0.01*D*D*D) -1 / D50) );
float ALPHAGR = (float) ( sqrt( 1.9f * sqrt(D50) ) * VSG * BETA / (
sqrt(FS) * sqrt(OFRQ) ) );
pErosion[nr].DEPGR = ALPHAGR * ( TCGR - QSGRN );
}
void CErosion::ChannelErosion( int nr )
{
float D50 = pErosion[nr].D50;
float OFRQ = pErosion[nr].OFRQ;
float KECG = pErosion[nr].KECG;
141
float FS = pErosion[nr].Fs;
float Qab = pErosion[nr].Qab;
float TO;
if( B/h>= 30 )
TO = RhoW * G * h * J;
else
TO = RhoW * G * R * I;
float VS = pErosion[nr].VS;
float TOCDIMEN = (float) ( 0.4 * (VS*VS) / (G * D50) );
float TOC = (float) TOCDIMEN * (RhoS - RhoW) * G * D50;
pErosion[nr].EGR = (float) ( F * sqrt( 1.9 * sqrt(D50*1000.0f) ) *
sqrt(OFRQ) * KECG * ( TO-TOC ) / sqrt(FS) ); }
void CErosion::( int nr )
{
float QSI = pErosion[nr].QSI;
float EGR = pErosion[nr].EGR;
pErosion[nr].QSGRN = QSI + EGR;
}
void CErosion::( int nr )
{
float QSI = pErosion[nr].QSI;
float DEPGR = pErosion[nr].DEPGR;
float LI = pErosion[nr].LI;
pErosion[nr].QSGR = QSI + DEPGR * LI;
}
void CErosion::( int nr )
{
float QSGRN = pErosion[nr].QSGRN;
142
float DEPGR = pErosion[nr].DEPGR;
float LI = pErosion[nr].LI;
pErosion[nr].QSGR = QSGRN + DEPGR * LI;
}
void CErosion::( int nr )
{
pErosion[nr].QSGR = pErosion[nr].QSGRN;
}
void CErosion::( int nr )
{
float Foc = pErosion[nr]. Freight;
float Freight = pErosion[nr].Qab * Foc;
if( pErosion[nr]. Freight - Freight < 0 )
Freight = pErosion[nr].Freight;
pErosion[nr]. Freight -= Freight;
}
CErosionView* CErosion::GetActiveView()
{
CView *pView;
POSITION pos = GetFirstViewPosition();
if( pos )
pView = GetNextView(pos);
while( pos && !pView->IsWindowVisible() ) {
pView = GetNextView(pos);
}
if( pView )
if( pView->IsWindowVisible() )
return (CErosionView*) pView;
return ZERO;
143
}
BOOL CErosion::AddView( BOOL bVisible, CString WindowText )
{
class COrgDoc;
class CMapView;
class CZFD;
class CZFDView;
class CCompare;
class CCompareDlg;
class CErosion;
class CShapeView;
class CErosionView;
class CDlgLandMap;
class CMetab;
class CEventDoc : public CDocument
{
protected:
DECLARE_DYNCREATE(CEventDoc)
public:
CString m_EventName;
COrgDoc *m_pOrgDoc;
ERGSETSTRUCT *pErgSet;
CGG CERG, ColdERG;
CGG GGArea;
CVARDateCVAR;
CFormula vFormula;
CMemory MEMORY;
CNDate CN;
CMapView *grid_flood_view;
144
CMapView *grid_w_view;
CMapView *grid_damage_view;
CMapView *grid_ damagesum_view;
CGRID FloodGRID;
CGRID FloodHeightGRID;
CGRID WGRID;
CGRID Damage GRID;
CGRID DamageSumGRID;
CMyPtrList *Damage List;
CZFD *pZFD;
CErosion *pErosion;
CErosion *pMeta;
CShapeView *pShapeLandMap;
CDlgLandMap*pDlgLandMap;
PEGELSTRUCT *Levelq;
float delta_t;
float rain;
float rainsum;
float etp,evaporation;
public:
CEventDoc();
BOOL GetTimeSeriesRange( CDate *MinTotal, CDate *MinOverlap, CDate
*MaxOverlap, CDate *MaxTotal, CMyPtrList *pStatList=NULL );
BOOL createWGRID( CDate *Date, BOOL fMaxFlood );
BOOL createFloodGRID( CDate *Date, BOOL fMaxFlood );
void DeleteZFDWnd() { pZFD = NULL; };
void OnDeleteDisplayGrid( CMapView *View2delete );
BOOL UpdateViews();
void DeleteChildWindow( void **ptr );
145
virtual void OnCloseDocument();
virtual BOOL OnSaveDocument(LPCTSTR lpszPathName);
virtual BOOL OnOpenDocument(LPCTSTR lpszPathName);
virtual BOOL OnNewDocument();
public:
virtual ~CEventDoc();
#ifdef _DEBUG
virtual void AssertValid() const;
virtual void Dump(CDumpContext& dc) const;
#endif
protected:
DECLARE_MESSAGE_MAP()
};
CEventDoc_H
146