Luận án 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ỏ

hấ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 101 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 102 - 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: STT Bộ thông số mô hình Giá trị 1 Hệ số xói mòn tách rời liên rãnh (ail) 5,3 e-05 2 Hệ số khả năng xói mòn đất (ket) 1,12 3 Hệ số vận chuyển bùn cát liên rãnh (acl) 2,5 e-04 4 Hệ số xói mòn rãnh (Ci) 0.321 5 Hệ số vận chuyển bùn cát rãnh (ac) 1.05 103 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 104 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 105 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 ,c cm nd theo phƣơng đứng của hai đƣờng phân phối xác suất. Giá trị tính toán của ,c cm nd sẽ so sánh với giá trị tiêu chuẩn, khi 106 ,c cm n d tính toán ≥ ,c cm nd 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: , _ .c c c c m n tieu chuan c c m n d m n    (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: 107 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à 108 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. 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: Ei= a.i 2 .Io b (3-2) 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: 109 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 110 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: 2 1,60,0002i oE i I   (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. Phƣơng trình cơ bản cho xói mòn rãnh có dạng tổng quát nhƣ sau: Er = c.i 2 .Io d (3-4) 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: 111 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: 1,670,00082 r o E i I   (3-5) 112 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ố. 113 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ì 114 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. 115 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. 116 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 hoc̣ “Sƣ̉ duṇg đấ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 117 đấ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. 118 [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. 119 [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 120 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. 121 [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. 122 [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. 123 PHỤ LỤC 1. Các bảng phân chia lƣu vực Bảng 1: Phân chia các tiểu lƣu vực trong lƣu vực Nậm Sập Lƣu vực Diện tích (km 2 ) Chiều dài lƣu vực (km) Độ dốc Ix (%) Độ dốc Iy (%) Chiều dài sông chính (m) Độ cao trung bình (m) 1 14.08 12.02 141.18 271.8 1.171 806 2 37.76 29.85 85.09 437.6 1.265 439 8 13.44 106.67 8.16 33.8 126 661 3 2.56 2.70 185.32 532.1 947 601 4 15.36 12.58 101.30 474.2 1.221 640 6 10.88 47.30 36.28 187.7 230 829 5 12.80 11.69 98.60 503.8 1.095 779 7 7.68 60.95 17.30 94.1 126 691 9 35.20 40.14 54.35 449.5 877 639 12 7.04 12.87 92.04 379.5 547 280 16 9.60 56.14 566.65 452.3 171 308 14 4.48 6.71 107.35 327.6 668 760 13 8.32 9.31 102.80 375.4 894 493 11 12.16 11.45 82.33 349.2 1.062 785 15 13.44 11.39 156.23 359.2 1.180 302 17 7.68 12.21 87.57 263.1 629 363 20 73.60 93.88 35.57 264.0 784 408 25 9.60 10.16 114.73 394.8 945 714 22 7.04 9.51 105.44 349.4 740 347 27 3.20 5.32 130.43 268.9 602 633 21 10.24 14.34 99.81 282.9 714 815 18 16.64 20.10 59.03 211.3 828 410 24 7.04 8.35 134.92 273.2 843 258 30 1.28 4.83 91.93 335.4 265 836 19 7.04 8.08 79.94 204.3 871 435 10 48.00 77.42 26.72 261.4 620 868 26 24.32 21.28 96.02 305.1 1.143 353 23 11.52 12.37 79.21 283.2 931 557 29 28.16 40.81 34.44 344.8 690 1223 31 8.32 21.72 43.55 353.2 383 1170 33 17.92 75.93 13.93 241.1 236 789 37 7.04 12.55 62.44 360.1 561 849 35 5.76 6.58 149.06 406.8 875 826 124 Lƣu vực Diện tích (km 2 ) Chiều dài lƣu vực (km) Độ dốc Ix (%) Độ dốc Iy (%) Chiều dài sông chính (m) Độ cao trung bình (m) 39 11.52 26.61 430.45 301.0 433 836 38 1.92 4.24 144.74 443.8 453 723 28 14.72 20.70 50.02 265.5 711 902 34 22.40 51.03 29.48 295.4 439 1151 41 14.72 33.00 45.93 277.0 446 1038 32 16.00 79.21 15.05 286.8 202 1001 40 10.24 16.73 66.30 238.4 612 941 36 97.28 94.08 30.80 344.4 1.034 1081 48 60.80 182.58 582.17 481.5 333 1080 44 7.04 33.36 34.48 242.3 211 978 46 9.60 19.35 64.11 352.0 496 1128 45 5.76 8.29 93.62 297.3 695 775 43 10.88 27.41 44.02 305.7 397 1184 51 7.04 13.26 70.63 407.6 531 855 47 5.12 13.58 59.76 399.3 377 843 55 8.32 15.70 60.94 225.8 530 788 49 7.68 11.21 100.23 463.4 685 583 50 3.84 7.27 124.89 305.7 528 617 42 22.40 29.55 57.80 334.3 758 1.127 59 4.48 11.49 61.91 205.4 390 853 52 7.68 15.30 74.76 320.7 502 1081 56 4.48 7.74 106.37 346.3 579 816 58 5.12 11.77 91.29 326.3 435 685 53 20.48 57.37 36.79 335.0 357 1088 63 10.24 20.81 39.81 374.5 492 883 61 1.28 4.92 100.74 306.6 260 722 64 1.70 65.60 122.91 426.2 130 749 60 8.96 19.27 72.10 236.2 465 923 54 12.80 20.55 61.95 316.4 623 1135 57 26.24 64.47 25.30 196.0 407 967 65 16.64 28.30 43.31 282.2 588 942 62 7.04 14.49 64.68 280.1 486 917 66 7.68 20.43 56.78 183.3 376 732 69 8.96 14.18 72.41 284.9 632 922 68 6.40 11.57 71.42 286.3 553 847 67 24.96 43.48 41.60 280.5 574 1064 71 19.84 31.54 54.78 287.3 629 596 70 28.16 40.52 55.10 316.1 695 491 72 10.24 9.57 173.89 375.5 1.070 890 125 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 (km 2 ) Chiều dài lƣu vực (km) Độ dốc Ix(%) Độ dốc Iy (%) Chiều dài sông chính (m) Độ cao trung bình (m) 1 3.00 2.40 66 114 1.251 564 2 2.00 2.47 270 113 809 440 4 2.00 2.06 242 111 969 442 3 3.00 2.18 364 133 1.379 1150 5 6.50 6.60 175 97 985 723 8 4.00 3.71 244 108 1.079 976 6 5.00 5.21 189 118 960 824 7 3.50 2.78 259 114 1.261 1050 11 2.50 2.85 244 102 876 944 10 2.50 2.65 214 105 945 732 9 3.50 2.90 299 95 1.207 1139 18 2.00 2.50 231 96 801 1188 14 4.00 3.25 229 96 1.232 1.265 13 2.50 3.42 176 72 730 1.103 17 0.50 1.31 253 91 382 1.111 15 1.50 1.62 236 88 926 1.255 19 2.50 2.81 239 99 889 1.212 22 22.50 57.69 278 76 390 980 23 3.50 3.05 224 89 1.148 1.016 20 3.00 3.46 176 95 866 1.481 21 2.00 1.97 348 113 1.015 1.366 16 2.00 2.11 222 94 950 1.521 26 4.00 4.15 169 101 963 1.343 25 1.00 2.30 149 72 435 971 29 3.00 3.40 163 111 882 1.202 27 2.50 2.70 218 106 927 1.236 12 12.00 8.33 138 99 1.440 1.433 33 3.00 4.11 155 71 730 1.714 34 3.50 3.76 196 68 931 1.488 24 6.00 6.06 175 93 990 1.051 35 2.00 2.13 235 80 940 1.484 28 3.00 2.28 341 113 1.314 1.455 37 0.50 1.31 194 74 381 1.097 38 3.00 3.32 152 81 904 1.416 36 4.00 3.70 202 110 1.082 907 40 0.50 2.42 119 74 207 1.021 32 5.00 4.76 160 69 1.050 1.800 39 2.50 3.07 185 77 814 1.229 30 7.50 4.87 188 102 1.541 1.829 31 5.00 3.49 285 111 1.431 1.695 41 2.00 2.22 173 71 901 1.347 42 3.50 3.41 184 87 1.026 1.079 45 34.50 157.53 278 64 219 833 44 3.00 3.53 198 91 850 1.193 126 Lƣu vực Diện tích (km 2 ) Chiều dài lƣu vực (km) Độ dốc Ix(%) Độ dốc Iy (%) Chiều dài sông chính (m) Độ cao trung bình (m) 46 3.50 2.89 293 80 1.209 2.248 43 11.50 7.30 164 88 1.575 1.370 47 2.00 3.39 167 92 590 1.068 49 19.00 863.64 62 44 22 697 50 4.50 4.20 171 74 1.071 1.178 51 19.50 80.58 307 91 242 711 52 2.50 3.24 140 73 771 1.189 48 4.00 2.43 225 95 1.647 1.255 54 11.50 5.66 182 96 2.033 1.733 58 3.00 1.71 324 88 1.751 1.377 57 6.00 5.72 141 85 1.049 683 53 4.50 4.56 223 93 986 583 60 1.00 2.19 216 94 457 384 56 15.00 7.67 191 90 1.956 1.087 55 2.00 1.98 283 104 1.011 775 61 5.50 6.54 143 92 841 792 59 5.50 3.35 273 90 1.643 1.039 62 9.50 8.23 130 102 1.155 727 64 4.00 3.89 204 114 1.028 692 63 5.00 4.73 213 101 1.057 573 66 2.00 1.91 256 76 1.049 806 67 2.00 1.61 319 98 1.245 580 65 4.00 3.68 187 77 1.087 812 69 5.00 3.98 250 95 1.257 656 71 0.50 1.34 217 100 374 505 70 1.50 1.24 420 90 1.208 442 68 4.00 3.78 208 87 1.057 411 76 1.50 4.59 116 93 327 291 73 2.50 4.48 169 84 558 361 77 2.00 2.71 190 84 737 447 75 0.50 1.01 231 83 495 226 72 4.50 5.88 174 78 765 782 79 1.50 2.34 248 80 641 213 81 1.00 2.15 201 64 465 170 82 3.00 3.38 241 104 888 519 74 2.50 2.81 203 82 891 801 78 11.50 9.64 156 89 1.193 301 80 2.50 2.66 178 82 940 878 84 7.00 9.20 130 82 761 322 87 6.00 74.07 210 47 81 154 85 2.00 2.01 223 92 993 621 86 2.00 2.36 322 113 849 652 89 5.00 8.88 98 71 563 681 83 2.50 2.83 214 85 882 456 92 1.00 1.79 263 76 560 327 88 2.00 2.13 233 105 940 599 93 4.50 4.69 221 110 960 779 127 Lƣu vực Diện tích (km 2 ) Chiều dài lƣu vực (km) Độ dốc Ix(%) Độ dốc Iy (%) Chiều dài sông chính (m) Độ cao trung bình (m) 90 5.00 6.51 116 68 768 571 91 4.50 4.57 202 86 984 510 94 4.50 5.81 178 83 774 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" ); 128 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) ); 129 } 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; 130 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 { } } } 131 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; 132 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; 133 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; 134 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; 135 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 ); 136 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 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 ) { 137 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 ); 138 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 ); 139 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; 140 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 ) { 141 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; 142 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; 143 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; 144 } 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; 145 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 ); 146 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