Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng và biến động trầm tích đáy cho vùng biển Vịnh Hạ Long

Tạp chí Khoa học đhqghn, KHTN & CN, T.xxII, Số 1PT., 2006 Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng và biến động trầm tích đáy cho vùng biển vịnh Hạ Long Đinh Văn Ưu Trung tâm Động lực và Môi tr−ờng Biển Tóm tắt. Đã phát triển và ứng dụng hệ thống mô hình thuỷ động lực và vận chuyển vật chất lơ lửng đối với vùng biển vịnh Hạ Long. Hệ thống này bao gồm các mô hình 3 chiều (3D) thuỷ-nhiệt động lực và vận chuyển vật chất lơ lửng, mô hình lớp biên đáy. Những kết quả ban đầu cho thấy hệ thống các mô hình có thể ứng dụng cho các khu vực biển có điều kiện địa hình và thuỷ động lực phức tạp nh− vịnh Hạ Long cũng nh− các vùng cửa sông ven biển nhằm mục đích xây dựng hệ thống mô hình monitoring và dự báo môi tr−ờng biển. Từ khoá: hệ thống mô hình 3D,vật chất lơ lửng, Vịnh Hạ Long 1. Đặt vấn đề Sự hiện diện của các chất lơ lửng trong n−ớc đ−ợc nghiên cứu thông qua phân tách toàn bộ lớp n−ớc thành hai phần: phần n−ớc nằm trên có nồng độ t−ơng đối thấp và phần sát đáy có nồng độ cao (lớp đáy lỏng). Đối với phần trên, chúng ta có thể sử dụng hệ các ph−ơng trình bình l−u - khuếch tán vật chất áp dụng cho nồng độ chất lơ lửng, còn đối với phần sát đáy thì sử dụng ph−ơng trình biến đổi độ dày của toàn lớp trầm tích. Mô hình này có thể áp dụng đồng thời cũng nh− tách biệt vì vậy dễ dàng phát triển cho các loại bài toán khác nhau từ nghiên cứu vận chuyển phù sa, xác định nguồn gốc trầm tích đáy, bồi xói đáy sông, biển đến các bài toán lan truyền ô nhiễm. Trong khuôn khổ của bài báo này, chúng tôi giới thiệu các kết quả b−ớc đầu giải bài toán tính biến động lớp trầm tích đáy, và chất lơ lửng là phù sa. Các kết quả thu đ−ợc là cơ sở để phát triển ứng dụng cho bài toán vận chuyển các chất lơ lửng trong n−ớc biển và khả năng ảnh h−ởng của chúng lên chất l−ợng môi tr−ờng n−ớc và trầm tích đáy. Việc ứng dụng mô hình 3D ph−ơng trình bình l−u - khuếch tán nghiên cứu lan truyền vật chất cho phép giải quyết bài toán một cách chính xác hơn so với lớp các mô hình 2D và tựa 3D tr−ớc đây. Việc áp dụng ph−ơng pháp thể tích hữu hạn trong mô hình số cũng góp phần nâng cao khả năng ứng dụng của mô hình. 1.1. Mô hình vận chuyển chất lơ lửng trong lớp n−ớc Sử dụng ph−ơng trình bình l−u - khuếch tán đầy đủ đối với nồng độ chất lơ lửng: 11 Đinh Văn Ưu 12 z c zy c yx c xzyx cw z cv y cu xt c zyx zyx ∂ ∂ ∂ ∂+∂ ∂ ∂ ∂+∂ ∂ ∂ ∂+=∂ ∂+∂ ∂+∂ ∂−= =∂ ∂+∂ ∂+∂ ∂+∂ ∂ λλλφϕϕϕφ ][ )()()( (1) với hàm nguồn bao gồm suất nhập (sản sinh), xuất (tiêu hủy) và lắng đọng: ).( mcIS r∇−+=φ . (2) Đối với chất lơ lửng là phù sa - một hợp phần tựa bền vững với các nguồn xuất - nhập hầu nh− chỉ xẩy ra trên biên, vì vậy chỉ cần chú ý duy nhất đến quá trình lắng đọng. Tr−ớc mắt các điều kiện biên đối với các biên biển hở có thể cho giá trị không đổi về nồng độ hoặc l−u l−ợng phù sa dựa theo đánh giá thực tế thông th−ờng có giá trị không đáng kể. Các điều kiện t−ơng tự cũng có thể áp dụng đối với các biên cửa sông với những giá trị nhất định. Đối với mặt phân cách giữa lớp n−ớc và lớp đáy, các thông l−ợng trao đổi đ−ợc tính thông qua quá trình bứt xói và lắng đọng trầm tích. Suất lắng đọng qua biên này đ−ợc tính dựa vào vận tốc lắng đọng, nồng độ chất lơ lửng tại chỗ và giá trị vận tốc động lực t−ơng đối so với giá trị tới hạn cho phép lắng đọng xuống đáy. ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −≈⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛−= d bs d bs cwu ucwD τ τ11 2 * * , (3) với điều kiện u* < u*d hoặc t−ơng ứng τ < τd. Giá trị của vận tốc lắng đọng của chất lơ lửng phụ thuộc một cách phức tạp vào đặc tr−ng của trầm tích và yếu tố động lực học. Van Rijn (1984) đã đ−a ra một công thức thực nghiệm sau đây tính theo kích th−ớc hạt d, tỷ lệ giữa mật độ trầm tích và mật độ n−ớc, s, và độ nhớt động học, η: ( ) η18 12 −= sgdws . (4) thông th−ờng giá trị s ≈ 2,65 và η ≈ 1,5.10-6 m2/s. Giá trị d đ−ợc xác định theo công thức: ( )( )[ 50251011,01 dTd s ]−−−= σ , (5) trong đó ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −= 16 50 50 84 2 1 d d d d sσ là độ phân tán của kích th−ớc trầm tích, 2<sσ khi trầm tích t−ơng đối đồng nhất, crb crbcbT , ,,' τ ττ −= là tham số ứng suất phi thứ nguyên với cb,'τ Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng... 13 là ứng suất đáy do dòng chảy và crb,'τ là giá trị tới hạn đ−ợc tính theo công thức phụ thuộc vào kích th−ớc hạt và giá trị số Shields tới hạn: ( ) crscrb gd θρρτ 50, −= (6) Việc tính toán tham số Shields có thể sử dụng công thức của Van Rijn: , (7) ⎪⎩ ⎪⎨ ⎧ ≤< ≤<= − − 10*414,0 4*124,0 64,0 * 1 * DD DD crθ với biểu thức đối với kích th−ớc hạt phi thứ nguyên đặc tr−ng: ( ) 3/1 250* 1 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −= η sgdD . Nh− vậy kích th−ớc trầm tích 50dd = khi giá trị ứng suất t−ơng đối T > 25. Nghiên cứu vận tốc lắng đọng trong điều kiện dòng triều, Portela (1997) đã dẫn ra công thức đơn giản hơn chỉ phụ thuộc vào nồng độ chất lơ lửng: (8) ⎪⎩ ⎪⎨⎧ ≥= − − 31 31 /1,0).(02,0 /1,0).(002,0 mkgcsmc mkgcsm ws p Nồng độ chất lơ lửng gần đáy cb và c đ−ợc lấy theo nồng độ tại biên d−ới cùng lớp n−ớc hoặc giá trị trung bình cho lớp biên đáy. Giá trị của nồng độ này có thể tính theo quy luật hàm số mũ (Mayer, 1995): 1 0 0 1 − − ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −= HA w v s b v s eH A wcc (9) với H0 là độ dày lớp n−ớc có ảnh h−ởng, thông th−ờng độ dày này đ−ợc chọn bằng 1m và c là nồng độ trung bình trong lớp n−ớc đó. Hệ số khuếch tán Av đối với chất lơ lửng mịn đ−ờng kính nhỏ hơn 20 àm đ−ợc xem là giảm tuyến tính trong lớp sát đáy từ 3.10-10 m2/s đến 3.10-11 m2/s (Pohlmann, 1994). Cũng theo Pohlmann (1994) thì giá trị vận tốc tới hạn đối với lắng đọng có thể lấy nh− sau: (10) ( ) ⎪⎩ ⎪⎨ ⎧ ≤++ ≤ = − −−− −− smwsm smwsmwc smwsm u s ss s d /10.5/028,0 /10.510.5).(3,4log.02,0008,0 /10.5).(008,0 4 451 51 * f p Trong tr−ờng hợp đáy biển có sinh vật đáy, thì quá trình lắng đọng sinh học có thể tính nh− sau: bbiobio cwD = (11) với Đinh Văn Ưu 14 (12) smwbio /10.15,1 6−≈ Thông l−ợng bứt xói từ đáy có thể tính theo nhiều cách khác nhau đối với từng loại trầm tích đáy. Theo Pohlmann (1994), đối với đáy bùn thì suất bứt xói sẽ là: ( ) )/( 22*2* smtonuuCE ee −= (13) với hệ số Ce = 10-4 ton.s/m4 lấy theo kết quả thực nghiệm của Puls (1984) và Rodger et al (1985). Giá trị của vận tốc động lực bứt xói tới hạn đ−ợc lấy bằng 0,028 m/s. Nh− vậy đối với mô hình vận chuyển chất lơ lửng trong toàn lớp n−ớc, thông l−ợng vật chất qua biên sẽ là tổng đại số của hai hợp phần lắng đọng D và bứt xói E: FLCb = D - E (14) 1.2. Mô hình biến đổi độ dày lớp trầm tích đáy lỏng Sử dụng ph−ơng trình bảo toàn khối l−ợng để nghiên cứu biến đổi của độ dày lớp đáy lỏng: EDq t −+∇=∂ ∂ r.ξ (15) với là thông l−ợng vật chất vận chuyển ngang trong lớp đáy lỏng có thể xem đồng nhất với dòng vận chuyển trầm tích di đáy. Các hạng thức D và E t−ơng ứng suất lắng đọng và bứt xói trao đổi với lớp n−ớc nằm trên đã đ−ợc đề cập trong phần tr−ớc. Sau đây chúng ta sẽ tập trung phân tích các công thức tính toán đối với suất di đáy. qr Thông th−ờng suất di đáy qr.∇ đ−ợc tính theo véc tơ của dòng vận chuyển trầm tích di đáy , đại l−ợng này đ−ợc tính theo nhiều công thức khác nhau nh− Piter- Mayer, Van Rijn, v.v... phụ thuộc vào ứng suất đáy do tác động của sóng và dòng chảy và đặc tr−ng của trầm tích.. qr Công thức Piter-Mayer (1948) dạng tổng quát có thể viết nh− sau: ( 2/32/350 2/1 047,08 −⎥⎦ ⎤⎢⎣ ⎡ ⎟⎟⎠ ⎞⎜⎜⎝ ⎛ −= àθρ ρρ dgq s ) , (16) với ( ) 50 , gds cb ρρ τθ −= và 2/3 ' ⎟⎠ ⎞⎜⎝ ⎛= C Cà (17) trong đó ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= csk hC , 12log18 và ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛= 90 12log18' d hC , h là độ sâu lớp n−ớc, ks,c là tham số nhám. Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng... 15 Nếu bỏ qua hạng thức 0,047, công thức (16) trở về dạng đơn giản sau: ( ) ( ) [ 12 2/32 * 2 * . 8 − − −= sm g uuq s cr ρρ ], (18) Cùng với các giả thiết về hai lớp nh− trên, Van Rijn (1984) đã đ−a ra công thức tính vận chuyển di đáy trong dạng sau: cr cr s v vDdgq θθθ θ ρ ρρ >⎥⎦ ⎤⎢⎣ ⎡ ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ −⎥⎦ ⎤⎢⎣ ⎡ ⎟⎟⎠ ⎞⎜⎜⎝ ⎛ −= ,/1053,0 3,0*2/350 2/1 , (19) Giá trị ứng suất đáy phục vụ tính tham số Shields, θ, đ−ợc lấy từ kết quả mô hình thuỷ nhiệt động lực 3D: 2 *, ucb ρτ = . Nh− vậy, việc đ−a mô hình lớp trầm tích đáy lỏng vào mô hình 3D thuỷ động lực vừa cho phép cung ứng các điều kiện biên trao đổi chất qua biên cho phép giải bài toàn vận chuyển chất lơ lửng, đồng thời kết quả tính toán ứng suất đáy lại cung cấp đầu vào cần thiết cho mô hình lớp biên đáy lỏng. Trong tr−ờng hợp phát triển mô hình lớp biên đáy bao gồm các tầng trầm tích khác nhau, chúng ta có thể hoàn thiện mô hình biến động đáy đáp ứng các yêu cầu về nghiên cứu bồi, xói cũng nh− nguồn gốc, chất l−ợng môi tr−ờng. 2. Một số kết quả ứng dụng mô hình đối với vùng vịnh Hạ Long Với những kết quả phát triển mô hình 3D thuỷ động lực học n−ớc nông ven bờ, chúng tôi đã thử nghiệm kết hợp giải bài toán lan truyền chất lơ lửng dạng vệt cho vùng biển có địa hình phức tập nh− vịnh Hạ long (hình 1). Q u ả n g N i n h L ạ n g S ơ n B ắ c G ia n g T P . H ả i P h ò n g Q u ố c lộ 4 B Q u ố c lộ 3 1 Q uố c lộ 1 8 Đ − ờ n g x e g o ò n g Q u ố c l ộ 4 B Q u ố c l ộ 4 B S ô n g B a C h e T R U N G Q U O C Hình 1. Bản đồ vùng biển Quảng Ninh và vịnh Hạ Long Đinh Văn Ưu 16 Kết quả mô phỏng vệt loang do nguồn điểm đối với chất lơ lửng, áp dụng cho hai tr−ờng gió đối lập nhau trong 2 mùa, cho thấy vai trò của hoàn l−u d− có ý nghĩa rất quan trọng đối với khả năng lan truyền chất vào dải ven bờ và thoát khỏi vịnh (Đinh Văn Ưu và ctv, 2005). Với một đặc tr−ng thuỷ động lực nhất định, khả năng duy trì các chất lơ lửng tại các tầng sâu và lắng đọng xuống đáy là đáng kể. Cùng với dòng dự, các kết quả mô hình hoá dòng chảy tổng hợp cho thấy sự biến động của tr−ờng hoàn l−u hết sức mạnh mẽ, có thể thấy điều này trên hình 3 dẫn ra dòng chảy tổng hợp mùa hè trong 2 pha triều ng−ợc nhau (hình 2). Hình 2. Kết quả tính hoàn tổng hợp mùa hè trong hai pha triều ng−ợc nhau Hình 3. Kết quả tính vận chuyển chất lơ lửng trong mùa hè sau 36 h (trái) và 42 h Với sự biến động mạnh của dòng chảy tổng hợp, bên cạnh sự lan truyền chất lơ lửng theo dòng d− nh− đã nhận đ−ợc tr−ớc đây (Đinh Văn Ưu, 2003, Đinh Văn Ưu và ctv, 2005), sự biến đổi phân bố chất lơ lửng trong từng pha triều cũng cho thấy có sự biến đổi đáng kể của vị trí và hình dáng vệt. Chúng ta có thể nhận ra sự biến động phân bố của nồng độ chất lơ lửng cách nhau trong khoảng 6 giờ ven bờ vịnh Hạ long đ−ợc thể hiện trên hình 3, trong mùa hè. Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng... 17 Hình 4. Kết quả tính vận chuyển chất lơ lửng trong mùa hè (trái) và mùa đông sau 72 h Theo kết quả thu đ−ợc cho một khoảng thời gian dài với nhiều chu kỳ triều (5-7 ngày) có thể thấy rằng chất lơ lửng trong lớp n−ớc có xu thế bám sát bờ Hạ Long trong mùa hè và đi về cả hai phía cửa Lục và Cát Bà trong mùa đông (hình 4) Trên hình 5 dẫn ra độ dày của lớp trầm tích lắng đọng xuống đáy và khả năng bứt xói từ đáy đối với trầm tích bùn mịn cho thấy, trong các điều kiện thuỷ động lực cụ thể tồn tại những khu vực có khả năng xẩy ra hiện t−ợng lắng đọng hoặc bứt xói khác nhau đối với từng mùa. Điều này đồng nghĩa với khả năng xác định và dự báo các khu vực tập trung chất vệt trên đáy đối với bài toán chất l−ợng môi tr−ờng và khu vực có thể xẩy ra hiện t−ợng bồi, xói nếu chât lơ lửng đ−ợc xét ở đây là l−ợng phù sa trong n−ớc. Hình 5. Kết quả tính phân bố vùng lắng đọng và bứt xói trên đáy (mm) đối với trầm tích bùn mịn d = 10 àm sau 4 ngày đêm trong mùa hè Đinh Văn Ưu 18 Việc triển khai mô hình 3D kết hợp mô hình lớp biên đáy đã mở ra khả năng ứng dụng của mô hình này trong giải quyết các bài toán thuỷ thạch động lực và môi tr−ờng vùng biển ven bờ. Tuy nhiên để có thể ứng dụng rộng rãi trong t−ơng lai, bên cạnh từng b−ớc hoàn thiện mô hình 3D thuỷ động lực học cần phát triển mô hình lớp biên đáy theo h−ớng chi tiết hoá càng tầng trầm tích và cơ chế vận chuyển của chúng. Lời cảm ơn Các kết quả thu đ−ợc là một phần của Đề tài NC Cơ bản 722004, tác giả chân thành cảm ơn sự hỗ trợ đó. Tài liệu tham khảo 1. Đinh Văn Ưu, “Các kết quả phát triển và ứng dụng mô hình ba chiều (3D) thuỷ nhiệt động lực biển ven và n−ớc nông ven bờ Quảng Ninh”, Tạp chí Khoa học ĐHQG Hà Nội, XIX, 1, 2003, tr.108-117. 2. Đinh Văn Ưu, Đoàn Văn Bộ, Hà Thanh H−ơng, Phạm Hoàng Lâm, ứng dụng mô hình dòng chảy ba chiều (3D) nghiên cứu quá trình lan truyền các chất lơ lửng tại vùng biển ven bờ Quảng Ninh, Tuyển tập Hội nghị Cơ học thủy khí toàn quốc 2005 (xuất bản 1-2006) 3. Mayer B., Ein Dreinimensionales, Numerishes Schwebstoff-Transprtmodell mit Anwendung auf die Deutshe Bucht, Dissertation, GKSS 95/E/59, 1995. 4. Mayer-Piter E. and Muller, R., Formulas for Bed-Load Transport, Sec. Int. IAHR Congress, Stockholm, Sweden, 1948. 5. Pohlmann T. and Puls W., Currents and Transport In Water. In J. Sundermann (Editor), Circulation and Contaminant Fluxes In the North Sea. Springer Verlag Science, Michigan, 1994, 555-605. 6. Portela L.I., Effect of settling velocity on the modelling of suspended sediment transport. In Acinas, J.R. and C.A. Brebbia (Editors), Computer Modelling of Seas and Coastal Regions III, Computational Mechanics Publications, Sothampton, 1997, 381-390. 7. Puls W, Erosion characteristics of estuarine muds. Hydraulics Research Rep. IT265, Wallingford, UK, 1984. 8. Rodger J.G. and Odd N.V.M., Sludge disposal In coastal waters, Hydraulics Research Rep. SR 70, Wallingford, UK., 1985. 9. Van Rjin L.C., Sediment Transport, Journal of Hydraulic Engineering, Vol 110, N. 10, 1984, 11,12. Phát triển và ứng dụng mô hình tính toán vận chuyển chất lơ lửng... 19 VNU. JOURNAL OF SCIENCE, Nat., Sci., & Tech., T.xXII, n01AP., 2006 Development and application of the SPM transport and bed layer changes modelling system for the Halong Bay area Dinh Van Uu Marine Dynamics and Environment Center Vietnam National University A hydrodynamic and uspended matter transport modelling system was developed and applied to the Halong Bay area. The system is concluded: three- dimensional (3D) thermo-hydrodynamic model and suspended matter transport model, moving bottom layer model. The preliminary results show that the modelling system could be applied to the geological and hydrodynamical complexe region as Halong Bay and another estuarine areas to create a marine environmental monitoring and prediction system. Keywords: 3D modelling system, suspended matter, Halong Bay